diff --git a/.clang-format b/.clang-format
new file mode 100644
index 0000000000000000000000000000000000000000..542c97e89adce40066f0b9a96cc8be9746640d65
--- /dev/null
+++ b/.clang-format
@@ -0,0 +1,5 @@
+---
+Language: Cpp
+BasedOnStyle: Google
+KeepEmptyLinesAtTheStartOfBlocks: true
+...
diff --git a/INSTALL.swift b/INSTALL.swift
index c18142c1a62c7c8a011ead2ee8d0d869c9e072ec..a4a2a9208c9443aacc796c87dfe1b96ee88fbb34 100644
--- a/INSTALL.swift
+++ b/INSTALL.swift
@@ -30,7 +30,7 @@ or:
./configure CC=icc
to use an Intel compiler. The main "programs" can be found in the "examples/"
-directory.
+directory. See README for run parameters.
SWIFT has been successfully built and tested with the following compilers:
@@ -105,3 +105,14 @@ among the different cores on each computing node.
DOXYGEN: the doxygen library is required to create the SWIFT API
documentation.
+
+
+
+ SWIFT Coding style
+ ==================
+
+The SWIFT source code is using a variation of the 'Google' style. The
+script 'format.sh' in the root directory applies the clang-format-3.8
+tool with our style choices to all the SWIFT C source file. Please
+apply the formatting script to the files before submitting a merge
+request.
diff --git a/Makefile.am b/Makefile.am
index 9d4f6371b7aa37a1750c5fb8dbed17f9ff48442e..fb4eb5f6d6b63a7d0e034e0a3202ac61066e6e25 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -1,7 +1,6 @@
-
# This file is part of SWIFT.
-# Copyright (c) 2012 pedro.gonnet@durham.ac.uk.
-# 2015 matthieu.schaller@durham.ac.uk.
+# Copyright (c) 2012 pedro.gonnet@durham.ac.uk
+# 2015 matthieu.schaller@durham.ac.uk
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@@ -23,4 +22,4 @@ ACLOCAL_AMFLAGS = -I m4
SUBDIRS = src examples doc tests
# Non-standard files that should be part of the distribution.
-EXTRA_DIST = INSTALL.swift
+EXTRA_DIST = INSTALL.swift .clang-format format.sh
diff --git a/README b/README
index 59c362415a9199d08fb6dfbb1ed044c66e647254..a010c4595a6e3a113879659390b373553c6e07d3 100644
--- a/README
+++ b/README
@@ -14,6 +14,7 @@ See INSTALL.swift for install instructions.
Usage: swift [OPTION] PARAMFILE
Valid options are:
+ -a Pin runners using processor affinity
-c Run with cosmological time integration
-d Dry run. Read the parameter file, allocate memory but does not read
the particles from ICs and exit before the start of time integration.
diff --git a/configure.ac b/configure.ac
index 59eea9c4e55d5aa22fd5a2f41e7f9e0a97ad9571..497107121753f212dd8b07f5a8e8eed7acdf82b5 100644
--- a/configure.ac
+++ b/configure.ac
@@ -329,7 +329,7 @@ fi
AM_CONDITIONAL([HAVEPARALLELHDF5],[test "$have_parallel_hdf5" = "yes"])
# Check for setaffinity.
-AC_CHECK_FUNC(pthread_setaffinity_np, AC_DEFINE([HAVE_SETAFFINITY],[true],
+AC_CHECK_FUNC(pthread_setaffinity_np, AC_DEFINE([HAVE_SETAFFINITY],[1],
[Defined if pthread_setaffinity_np exists.]) )
AM_CONDITIONAL(HAVESETAFFINITY,
[test "$ac_cv_func_pthread_setaffinity_np" = "yes"])
diff --git a/examples/CosmoVolume/cosmoVolume.yml b/examples/CosmoVolume/cosmoVolume.yml
index 6075ad99799ddb671e34a5b590713e8441be56f6..45d32ae081823374d523a6105b196a4e1e7aa3c6 100644
--- a/examples/CosmoVolume/cosmoVolume.yml
+++ b/examples/CosmoVolume/cosmoVolume.yml
@@ -1,4 +1,3 @@
-
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
@@ -17,9 +16,9 @@ Scheduler:
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
- time_end: 1. # The end time of the simulation (in internal units).
+ time_end: 1e-4 # The end time of the simulation (in internal units).
dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
@@ -32,6 +31,10 @@ Snapshots:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml
index 6d05cb5b43c541490432c8f283b662c9d7c8d184..cfc6f0ae6d15a50471900e0bd09113e942fc116d 100644
--- a/examples/ExternalPointMass/externalPointMass.yml
+++ b/examples/ExternalPointMass/externalPointMass.yml
@@ -1,4 +1,3 @@
-
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.9885e33 # Grams
@@ -32,6 +31,10 @@ Snapshots:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/SedovBlast/sedov.yml b/examples/SedovBlast/sedov.yml
index b5ad8d855535e70133464fa6a117db37a2c9fb8b..1f491c248ffbd590dda5b5dcdf86b3eca27c8ce7 100644
--- a/examples/SedovBlast/sedov.yml
+++ b/examples/SedovBlast/sedov.yml
@@ -1,4 +1,3 @@
-
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
@@ -32,6 +31,10 @@ Snapshots:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-3 # Time between statistics output
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/SodShock/sodShock.yml b/examples/SodShock/sodShock.yml
index 7bb1895c15c80c6a4e54248d1d2bf85dea278ce1..d61c593f10387ad1a6bfc9cbbe9e89f05f84c52c 100644
--- a/examples/SodShock/sodShock.yml
+++ b/examples/SodShock/sodShock.yml
@@ -1,4 +1,3 @@
-
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
@@ -32,6 +31,10 @@ Snapshots:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/UniformBox/uniformBox.yml b/examples/UniformBox/uniformBox.yml
index ce989d22ad8a2310ef86303e1b7b6788fdcc427f..15d956a5d752262db91208f689005a7dc85f30d5 100644
--- a/examples/UniformBox/uniformBox.yml
+++ b/examples/UniformBox/uniformBox.yml
@@ -1,4 +1,3 @@
-
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
@@ -32,6 +31,10 @@ Snapshots:
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/main.c b/examples/main.c
index 19d2d7a7871200ba6e4902d5918a59511bf32b47..9f371d64dde9c66f118477beb145f426946fc9d9 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -27,10 +27,10 @@
/* Some standard headers. */
#include <fenv.h>
-#include <unistd.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
+#include <unistd.h>
/* MPI headers. */
#ifdef WITH_MPI
@@ -52,6 +52,7 @@ void print_help_message() {
printf("\nUsage: swift [OPTION] PARAMFILE\n\n");
printf("Valid options are:\n");
+ printf(" %2s %8s %s\n", "-a", "", "Pin runners using processor affinity");
printf(" %2s %8s %s\n", "-c", "", "Run with cosmological time integration");
printf(
" %2s %8s %s\n", "-d", "",
@@ -88,6 +89,7 @@ void print_help_message() {
* @brief Main routine that loads a few particles and generates some output.
*
*/
+
int main(int argc, char *argv[]) {
struct clocks_time tic, toc;
@@ -121,24 +123,16 @@ int main(int argc, char *argv[]) {
fflush(stdout);
#endif
+/* Let's pin the main thread */
#if defined(HAVE_SETAFFINITY) && defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
- if ((ENGINE_POLICY) & engine_policy_setaffinity) {
- /* Ensure the NUMA node on which we initialise (first touch) everything
- * doesn't change before engine_init allocates NUMA-local workers.
- * Otherwise, we may be scheduled elsewhere between the two times.
- */
- cpu_set_t affinity;
- CPU_ZERO(&affinity);
- CPU_SET(sched_getcpu(), &affinity);
- if (sched_setaffinity(0, sizeof(cpu_set_t), &affinity) != 0) {
- error("failed to set entry thread's affinity");
- }
- }
+ if (((ENGINE_POLICY)&engine_policy_setaffinity) == engine_policy_setaffinity)
+ engine_pin();
#endif
/* Welcome to SWIFT, you made the right choice */
if (myrank == 0) greetings();
+ int with_aff = 0;
int dry_run = 0;
int dump_tasks = 0;
int with_cosmology = 0;
@@ -153,7 +147,10 @@ int main(int argc, char *argv[]) {
/* Parse the parameters */
int c;
- while ((c = getopt(argc, argv, "cdef:gGhst:v:y:")) != -1) switch (c) {
+ while ((c = getopt(argc, argv, "acdef:gGhst:v:y:")) != -1) switch (c) {
+ case 'a':
+ with_aff = 1;
+ break;
case 'c':
with_cosmology = 1;
break;
@@ -393,6 +390,7 @@ int main(int argc, char *argv[]) {
message("%zi parts in %i cells.", s.nr_parts, s.tot_cells);
message("%zi gparts in %i cells.", s.nr_gparts, s.tot_cells);
message("maximum depth is %d.", s.maxdepth);
+ fflush(stdout);
}
/* Verify that each particle is in it's proper cell. */
@@ -419,8 +417,9 @@ int main(int argc, char *argv[]) {
/* Initialize the engine with the space and policies. */
if (myrank == 0) clocks_gettime(&tic);
struct engine e;
- engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, engine_policies,
- talking, &prog_const, &hydro_properties, &potential);
+ engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
+ engine_policies, talking, &prog_const, &hydro_properties,
+ &potential);
if (myrank == 0) {
clocks_gettime(&toc);
message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
@@ -431,10 +430,6 @@ int main(int argc, char *argv[]) {
/* Write the state of the system before starting time integration. */
if (!dry_run) engine_dump_snapshot(&e);
- /* Now that everything is ready, no need for the parameters any more */
- free(params);
- params = NULL;
-
/* Init the runner history. */
#ifdef HIST
for (k = 0; k < runner_hist_N; k++) runner_hist_bins[k] = 0;
@@ -521,19 +516,20 @@ int main(int argc, char *argv[]) {
int count = 0;
for (int l = 0; l < e.sched.nr_tasks; l++)
if (!e.sched.tasks[l].skip && !e.sched.tasks[l].implicit) {
- fprintf(file_thread, " %03i %i %i %i %i %lli %lli %i %i %i %i %i\n",
- myrank, e.sched.tasks[l].last_rid, e.sched.tasks[l].type,
- e.sched.tasks[l].subtype, (e.sched.tasks[l].cj == NULL),
- e.sched.tasks[l].tic, e.sched.tasks[l].toc,
- (e.sched.tasks[l].ci != NULL) ? e.sched.tasks[l].ci->count
- : 0,
- (e.sched.tasks[l].cj != NULL) ? e.sched.tasks[l].cj->count
- : 0,
- (e.sched.tasks[l].ci != NULL) ? e.sched.tasks[l].ci->gcount
- : 0,
- (e.sched.tasks[l].cj != NULL) ? e.sched.tasks[l].cj->gcount
- : 0,
- e.sched.tasks[l].flags);
+ fprintf(
+ file_thread, " %03i %i %i %i %i %lli %lli %i %i %i %i %i\n",
+ myrank, e.sched.tasks[l].last_rid, e.sched.tasks[l].type,
+ e.sched.tasks[l].subtype, (e.sched.tasks[l].cj == NULL),
+ e.sched.tasks[l].tic, e.sched.tasks[l].toc,
+ (e.sched.tasks[l].ci != NULL) ? e.sched.tasks[l].ci->count
+ : 0,
+ (e.sched.tasks[l].cj != NULL) ? e.sched.tasks[l].cj->count
+ : 0,
+ (e.sched.tasks[l].ci != NULL) ? e.sched.tasks[l].ci->gcount
+ : 0,
+ (e.sched.tasks[l].cj != NULL) ? e.sched.tasks[l].cj->gcount
+ : 0,
+ e.sched.tasks[l].flags);
fflush(stdout);
count++;
}
@@ -552,8 +548,8 @@ int main(int argc, char *argv[]) {
FILE *file_thread;
file_thread = fopen(dumpfile, "w");
/* Add some information to help with the plots */
- fprintf(file_thread, " %i %i %i %i %lli %lli %i %i %i %lli\n", -2, -1, -1, 1,
- e.tic_step, e.toc_step, 0, 0, 0, cpufreq);
+ fprintf(file_thread, " %i %i %i %i %lli %lli %i %i %i %lli\n", -2, -1, -1,
+ 1, e.tic_step, e.toc_step, 0, 0, 0, cpufreq);
for (int l = 0; l < e.sched.nr_tasks; l++)
if (!e.sched.tasks[l].skip && !e.sched.tasks[l].implicit)
fprintf(
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 761e2f67c16766bc341b2bcb47f40ebbc4ea30ad..b3b9b52ad73fad2bde9b40b983f88acc3e985fbb 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -1,4 +1,3 @@
-
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1 # Grams
@@ -32,6 +31,10 @@ Snapshots:
UnitCurrent_in_cgs: 1 # Unit system for the outputs (Amperes)
UnitTemp_in_cgs: 1 # Unit system for the outputs (Kelvin)
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/plot_tasks.py b/examples/plot_tasks.py
index ffabdf06f9324fba35770e7fbba4cfadc8add770..e070bdb565d6d8305b7974137292910dd720afe2 100755
--- a/examples/plot_tasks.py
+++ b/examples/plot_tasks.py
@@ -56,8 +56,8 @@ pl.rcParams.update(PLOT_PARAMS)
# Tasks and subtypes. Indexed as in tasks.h.
TASKTYPES = ["none", "sort", "self", "pair", "sub", "init", "ghost", "drift",
- "kick", "send", "recv", "grav_pp", "grav_mm", "grav_up",
- "grav_down", "grav_external", "part_sort", "gpart_sort",
+ "kick", "kick_fixdt", "send", "recv", "grav_pp", "grav_mm",
+ "grav_up", "grav_down", "grav_external", "part_sort", "gpart_sort",
"split_cell", "rewait", "count"]
TASKCOLOURS = {"none": "black",
@@ -69,6 +69,7 @@ TASKCOLOURS = {"none": "black",
"ghost": "cyan",
"drift": "maroon",
"kick": "green",
+ "kick_fixdt": "green",
"send": "yellow",
"recv": "magenta",
"grav_pp": "mediumorchid",
diff --git a/examples/plot_tasks_MPI.py b/examples/plot_tasks_MPI.py
index cf591028b6dc3f724847b52c9efac4355d27f87e..44705a0164034f8a22f1f3aa860a64d3e7656d2e 100755
--- a/examples/plot_tasks_MPI.py
+++ b/examples/plot_tasks_MPI.py
@@ -62,8 +62,8 @@ pl.rcParams.update(PLOT_PARAMS)
# Tasks and subtypes. Indexed as in tasks.h.
TASKTYPES = ["none", "sort", "self", "pair", "sub", "init", "ghost", "drift",
- "kick", "send", "recv", "grav_pp", "grav_mm", "grav_up",
- "grav_down", "grav_external", "part_sort", "gpart_sort",
+ "kick", "kick_fixdt", "send", "recv", "grav_pp", "grav_mm",
+ "grav_up", "grav_down", "grav_external", "part_sort", "gpart_sort",
"split_cell", "rewait", "count"]
TASKCOLOURS = {"none": "black",
@@ -75,6 +75,7 @@ TASKCOLOURS = {"none": "black",
"ghost": "cyan",
"drift": "maroon",
"kick": "green",
+ "kick_fixdt": "green",
"send": "yellow",
"recv": "magenta",
"grav_pp": "mediumorchid",
diff --git a/format.sh b/format.sh
new file mode 100755
index 0000000000000000000000000000000000000000..5bf78c5bb663d5cf22178bff868912a059c5c4fe
--- /dev/null
+++ b/format.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+clang-format-3.8 -style=file -i src/*.[ch] src/*/*.[ch] src/*/*/*.[ch] examples/main.c tests/*.[ch]
diff --git a/src/Makefile.am b/src/Makefile.am
index 322a4823527f078da9142229c6e2b8d100379610..6d96fd01de5b4be002ad5dc03c65d5145ada7177 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -48,7 +48,8 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
# Include files for distribution, not installation.
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
- vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h \
+ vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h kick.h \
+ timestep.h drift.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
diff --git a/src/cell.c b/src/cell.c
index 29072173b060a936fe6f4f43c091581d522985b3..97a67583a7872999d7b7e7c5a0eb9787c83b4bd3 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -50,8 +50,8 @@
#include "atomic.h"
#include "error.h"
#include "gravity.h"
-#include "hydro_properties.h"
#include "hydro.h"
+#include "hydro_properties.h"
#include "space.h"
#include "timers.h"
diff --git a/src/cell.h b/src/cell.h
index a69a1a74648d76aeb13126d37a3c1265397dc2bf..65665c9dcfcd652cc688f652008044d618b10790 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -136,13 +136,13 @@ struct cell {
int hold, ghold;
/* Spin lock for various uses. */
- lock_type lock, glock;
+ swift_lock_type lock, glock;
/* ID of the previous owner, e.g. runner. */
int owner;
/* Momentum of particles in cell. */
- float mom[3], ang[3];
+ double mom[3], ang_mom[3];
/* Mass, potential, internal and kinetic energy of particles in this cell. */
double mass, e_pot, e_int, e_kin;
diff --git a/src/common_io.c b/src/common_io.c
index 143c1ca94160f2a43c40aa44384803ce721e4dbf..4399c83f7ebd2bb5355b7df9cc1fdecd301e36aa 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -147,8 +147,8 @@ void readAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data) {
*
* Calls #error() if an error occurs.
*/
-void writeAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data,
- int num) {
+void writeAttribute(hid_t grp, const char* name, enum DATA_TYPE type,
+ void* data, int num) {
hid_t h_space = 0, h_attr = 0, h_err = 0;
hsize_t dim[1] = {num};
@@ -186,7 +186,8 @@ void writeAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data,
*
* Calls #error() if an error occurs.
*/
-void writeStringAttribute(hid_t grp, char* name, const char* str, int length) {
+void writeStringAttribute(hid_t grp, const char* name, const char* str,
+ int length) {
hid_t h_space = 0, h_attr = 0, h_err = 0, h_type = 0;
h_space = H5Screate(H5S_SCALAR);
@@ -225,7 +226,7 @@ void writeStringAttribute(hid_t grp, char* name, const char* str, int length) {
* @param name The name of the attribute
* @param data The value to write
*/
-void writeAttribute_d(hid_t grp, char* name, double data) {
+void writeAttribute_d(hid_t grp, const char* name, double data) {
writeAttribute(grp, name, DOUBLE, &data, 1);
}
@@ -235,7 +236,7 @@ void writeAttribute_d(hid_t grp, char* name, double data) {
* @param name The name of the attribute
* @param data The value to write
*/
-void writeAttribute_f(hid_t grp, char* name, float data) {
+void writeAttribute_f(hid_t grp, const char* name, float data) {
writeAttribute(grp, name, FLOAT, &data, 1);
}
@@ -246,7 +247,7 @@ void writeAttribute_f(hid_t grp, char* name, float data) {
* @param data The value to write
*/
-void writeAttribute_i(hid_t grp, char* name, int data) {
+void writeAttribute_i(hid_t grp, const char* name, int data) {
writeAttribute(grp, name, INT, &data, 1);
}
@@ -256,7 +257,7 @@ void writeAttribute_i(hid_t grp, char* name, int data) {
* @param name The name of the attribute
* @param data The value to write
*/
-void writeAttribute_l(hid_t grp, char* name, long data) {
+void writeAttribute_l(hid_t grp, const char* name, long data) {
writeAttribute(grp, name, LONG, &data, 1);
}
@@ -266,7 +267,7 @@ void writeAttribute_l(hid_t grp, char* name, long data) {
* @param name The name of the attribute
* @param str The string to write
*/
-void writeAttribute_s(hid_t grp, char* name, const char* str) {
+void writeAttribute_s(hid_t grp, const char* name, const char* str) {
writeStringAttribute(grp, name, str, strlen(str));
}
diff --git a/src/common_io.h b/src/common_io.h
index 70ed25993cb110a1faf42a41c08d316c81b54271..ed1c96801c904d9952c7b3ca995b847afdc3fb43 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -52,10 +52,7 @@ enum DATA_TYPE {
*start a run or optional.
*
*/
-enum DATA_IMPORTANCE {
- COMPULSORY = 1,
- OPTIONAL = 0
-};
+enum DATA_IMPORTANCE { COMPULSORY = 1, OPTIONAL = 0 };
/**
* @brief The different particle types present in a GADGET IC file
@@ -88,14 +85,14 @@ void duplicate_hydro_gparts(struct part* const parts,
void readAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data);
-void writeAttribute(hid_t grp, char* name, enum DATA_TYPE type, void* data,
- int num);
+void writeAttribute(hid_t grp, const char* name, enum DATA_TYPE type,
+ void* data, int num);
-void writeAttribute_d(hid_t grp, char* name, double data);
-void writeAttribute_f(hid_t grp, char* name, float data);
-void writeAttribute_i(hid_t grp, char* name, int data);
-void writeAttribute_l(hid_t grp, char* name, long data);
-void writeAttribute_s(hid_t grp, char* name, const char* str);
+void writeAttribute_d(hid_t grp, const char* name, double data);
+void writeAttribute_f(hid_t grp, const char* name, float data);
+void writeAttribute_i(hid_t grp, const char* name, int data);
+void writeAttribute_l(hid_t grp, const char* name, long data);
+void writeAttribute_s(hid_t grp, const char* name, const char* str);
void createXMFfile(const char* baseName);
FILE* prepareXMFfile(const char* baseName);
@@ -108,7 +105,6 @@ void writeXMFline(FILE* xmfFile, char* fileName, char* partTypeGroupName,
char* name, size_t N, int dim, enum DATA_TYPE type);
void writeCodeDescription(hid_t h_file);
-void writeSPHflavour(hid_t h_file);
void writeUnitSystem(hid_t h_file, struct UnitSystem* us);
#endif
diff --git a/src/cycle.h b/src/cycle.h
index 1278c83e8b43324662bdeb0de75eec08faf4fd82..4925808f5a4c4ea7828bc1a6b7a9490d2d2ca255 100644
--- a/src/cycle.h
+++ b/src/cycle.h
@@ -334,7 +334,7 @@ typedef unsigned __int64 ticks;
extern "C"
#endif
ticks
- __getReg(int whichReg);
+ __getReg(int whichReg);
#pragma intrinsic(__getReg)
static __inline ticks getticks(void) {
@@ -481,9 +481,9 @@ INLINE_ELAPSED(inline)
/* MIPS ZBus */
#if HAVE_MIPS_ZBUS_TIMER
#if defined(__mips__) && !defined(HAVE_TICK_COUNTER)
+#include <fcntl.h>
#include <sys/mman.h>
#include <unistd.h>
-#include <fcntl.h>
typedef uint64_t ticks;
diff --git a/src/debug.c b/src/debug.c
index 53a03d66aee2c169a555ed00a2efa2d5b984066a..a3647915d96a9456cb6d8177d144dab56fd0b97e 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -24,8 +24,8 @@
#include "config.h"
#include "const.h"
-#include "part.h"
#include "debug.h"
+#include "part.h"
/* Import the right hydro definition */
#if defined(MINIMAL_SPH)
@@ -100,8 +100,9 @@ void printgParticle(struct gpart *gparts, long long int id, size_t N) {
void printParticle_single(struct part *p, struct xpart *xp) {
- printf("## Particle: id=%lld", p->id);
+ printf("## Particle: id=%lld ", p->id);
hydro_debug_particle(p, xp);
+ printf("\n");
}
#ifdef HAVE_METIS
diff --git a/src/drift.h b/src/drift.h
new file mode 100644
index 0000000000000000000000000000000000000000..05b09bb7910e8cddaf9fb24bb248f120f9db9eea
--- /dev/null
+++ b/src/drift.h
@@ -0,0 +1,101 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DRIFT_H
+#define SWIFT_DRIFT_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+#include "debug.h"
+#include "hydro.h"
+
+/**
+ * @brief Perform the 'drift' operation on a #gpart
+ *
+ * @param gp The #gpart to drift.
+ * @param dt The drift time-step
+ * @param timeBase The minimal allowed time-step size.
+ * @param ti_old Integer start of time-step
+ * @param ti_current Integer end of time-step
+ */
+__attribute__((always_inline)) INLINE static void drift_gpart(
+ struct gpart* gp, float dt, double timeBase, int ti_old, int ti_current) {
+ /* Drift... */
+ gp->x[0] += gp->v_full[0] * dt;
+ gp->x[1] += gp->v_full[1] * dt;
+ gp->x[2] += gp->v_full[2] * dt;
+
+ /* Compute offset since last cell construction */
+ gp->x_diff[0] -= gp->v_full[0] * dt;
+ gp->x_diff[1] -= gp->v_full[1] * dt;
+ gp->x_diff[2] -= gp->v_full[2] * dt;
+}
+
+/**
+ * @brief Perform the 'drift' operation on a #part
+ *
+ * @param p The #part to drift.
+ * @param xp The #xpart of the particle.
+ * @param dt The drift time-step
+ * @param timeBase The minimal allowed time-step size.
+ * @param ti_old Integer start of time-step
+ * @param ti_current Integer end of time-step
+ */
+__attribute__((always_inline)) INLINE static void drift_part(
+ struct part* p, struct xpart* xp, float dt, double timeBase, int ti_old,
+ int ti_current) {
+ /* Useful quantity */
+ const float h_inv = 1.0f / p->h;
+
+ /* Drift... */
+ p->x[0] += xp->v_full[0] * dt;
+ p->x[1] += xp->v_full[1] * dt;
+ p->x[2] += xp->v_full[2] * dt;
+
+ /* Predict velocities (for hydro terms) */
+ p->v[0] += p->a_hydro[0] * dt;
+ p->v[1] += p->a_hydro[1] * dt;
+ p->v[2] += p->a_hydro[2] * dt;
+
+ /* Predict smoothing length */
+ const float w1 = p->h_dt * h_inv * dt;
+ if (fabsf(w1) < 0.2f)
+ p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+ else
+ p->h *= expf(w1);
+
+ /* Predict density */
+ const float w2 = -3.0f * w1;
+ if (fabsf(w2) < 0.2f)
+ p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
+ else
+ p->rho *= expf(w2);
+
+ /* Predict the values of the extra fields */
+ hydro_predict_extra(p, xp, ti_old, ti_current, timeBase);
+
+ /* Compute offset since last cell construction */
+ xp->x_diff[0] -= xp->v_full[0] * dt;
+ xp->x_diff[1] -= xp->v_full[1] * dt;
+ xp->x_diff[2] -= xp->v_full[2] * dt;
+}
+
+#endif /* SWIFT_DRIFT_H */
diff --git a/src/engine.c b/src/engine.c
index 598f8e8ecb2a0540d5b82c7c3ce2be7a69e78463..59d83757830b62cd5c229b9937f9530462ee7773 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -29,11 +29,11 @@
#include <float.h>
#include <limits.h>
#include <sched.h>
+#include <stdbool.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
-#include <stdbool.h>
/* MPI headers. */
#ifdef WITH_MPI
@@ -56,22 +56,35 @@
#include "error.h"
#include "hydro.h"
#include "minmax.h"
+#include "parallel_io.h"
#include "part.h"
#include "partition.h"
-#include "parallel_io.h"
#include "serial_io.h"
#include "single_io.h"
#include "timers.h"
-const char *engine_policy_names[13] = {
- "none", "rand", "steal", "keep",
- "block", "fix_dt", "cpu_tight", "mpi",
- "numa_affinity", "hydro", "self_gravity", "external_gravity",
- "cosmology_integration"};
+const char *engine_policy_names[13] = {"none",
+ "rand",
+ "steal",
+ "keep",
+ "block",
+ "fix_dt",
+ "cpu_tight",
+ "mpi",
+ "numa_affinity",
+ "hydro",
+ "self_gravity",
+ "external_gravity",
+ "cosmology_integration"};
/** The rank of the engine as a global variable (for messages). */
int engine_rank;
+#ifdef HAVE_SETAFFINITY
+/** The initial affinity of the main thread (set by engin_pin()) */
+static cpu_set_t entry_affinity;
+#endif
+
/**
* @brief Link a density/force task to a cell.
*
@@ -112,6 +125,7 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c,
const int is_with_external_gravity =
(e->policy & engine_policy_external_gravity) ==
engine_policy_external_gravity;
+ const int is_fixdt = (e->policy & engine_policy_fixdt) == engine_policy_fixdt;
/* Am I the super-cell? */
if (super == NULL && (c->count > 0 || c->gcount > 0)) {
@@ -130,9 +144,14 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c,
c->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0, 0,
c, NULL, 0);
- /* Add the kick task. */
- c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0, 0, c,
- NULL, 0);
+ /* Add the kick task that matches the policy. */
+ if (is_fixdt) {
+ c->kick = scheduler_addtask(s, task_type_kick_fixdt, task_subtype_none,
+ 0, 0, c, NULL, 0);
+ } else {
+ c->kick = scheduler_addtask(s, task_type_kick, task_subtype_none, 0, 0,
+ c, NULL, 0);
+ }
if (c->count > 0) {
@@ -1799,8 +1818,9 @@ void engine_barrier(struct engine *e, int tid) {
/**
* @brief Mapping function to collect the data from the kick.
+ *
+ * @param c A super-cell.
*/
-
void engine_collect_kick(struct cell *c) {
/* Skip super-cells (Their values are already set) */
@@ -1808,15 +1828,13 @@ void engine_collect_kick(struct cell *c) {
/* Counters for the different quantities. */
int updated = 0, g_updated = 0;
- double e_kin = 0.0, e_int = 0.0, e_pot = 0.0;
- float mom[3] = {0.0f, 0.0f, 0.0f}, ang[3] = {0.0f, 0.0f, 0.0f};
- int ti_end_min = max_nr_timesteps, ti_end_max = 0;
+ int ti_end_min = max_nr_timesteps;
/* Only do something is the cell is non-empty */
if (c->count != 0 || c->gcount != 0) {
/* If this cell is not split, I'm in trouble. */
- if (!c->split) error("Cell has no super-cell.");
+ if (!c->split) error("Cell is not split.");
/* Collect the values from the progeny. */
for (int k = 0; k < 8; k++) {
@@ -1828,36 +1846,197 @@ void engine_collect_kick(struct cell *c) {
/* And update */
ti_end_min = min(ti_end_min, cp->ti_end_min);
- ti_end_max = max(ti_end_max, cp->ti_end_max);
updated += cp->updated;
g_updated += cp->g_updated;
+ }
+ }
+ }
+
+ /* Store the collected values in the cell. */
+ c->ti_end_min = ti_end_min;
+ c->updated = updated;
+ c->g_updated = g_updated;
+}
+
+/**
+ * @brief Collects the next time-step by making each super-cell recurse
+ * to collect the minimal of ti_end and the number of updated particles.
+ *
+ * @param e The #engine.
+ */
+void engine_collect_timestep(struct engine *e) {
+
+ int updates = 0, g_updates = 0;
+ int ti_end_min = max_nr_timesteps;
+ const struct space *s = e->s;
+
+ /* Collect the cell data. */
+ for (int k = 0; k < s->nr_cells; k++)
+ if (s->cells[k].nodeID == e->nodeID) {
+ struct cell *c = &s->cells[k];
+
+ /* Make the top-cells recurse */
+ engine_collect_kick(c);
+
+ /* And aggregate */
+ ti_end_min = min(ti_end_min, c->ti_end_min);
+ updates += c->updated;
+ g_updates += c->g_updated;
+ }
+
+/* Aggregate the data from the different nodes. */
+#ifdef WITH_MPI
+ {
+ int in_i[1], out_i[1];
+ in_i[0] = 0;
+ out_i[0] = ti_end_min;
+ if (MPI_Allreduce(out_i, in_i, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD) !=
+ MPI_SUCCESS)
+ error("Failed to aggregate t_end_min.");
+ ti_end_min = in_i[0];
+ }
+ {
+ unsigned long long in_ll[2], out_ll[2];
+ out_ll[0] = updates;
+ out_ll[1] = g_updates;
+ if (MPI_Allreduce(out_ll, in_ll, 2, MPI_LONG_LONG_INT, MPI_SUM,
+ MPI_COMM_WORLD) != MPI_SUCCESS)
+ error("Failed to aggregate energies.");
+ updates = in_ll[0];
+ g_updates = in_ll[1];
+ }
+#endif
+
+ e->ti_end_min = ti_end_min;
+ e->updates = updates;
+ e->g_updates = g_updates;
+}
+
+/**
+ * @brief Mapping function to collect the data from the drift.
+ *
+ * @param c A super-cell.
+ */
+void engine_collect_drift(struct cell *c) {
+
+ /* Skip super-cells (Their values are already set) */
+ if (c->drift != NULL) return;
+
+ /* Counters for the different quantities. */
+ double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
+ double mom[3] = {0.0, 0.0, 0.0}, ang_mom[3] = {0.0, 0.0, 0.0};
+
+ /* Only do something is the cell is non-empty */
+ if (c->count != 0 || c->gcount != 0) {
+
+ /* If this cell is not split, I'm in trouble. */
+ if (!c->split) error("Cell has no super-cell.");
+
+ /* Collect the values from the progeny. */
+ for (int k = 0; k < 8; k++) {
+ struct cell *cp = c->progeny[k];
+ if (cp != NULL) {
+
+ /* Recurse */
+ engine_collect_drift(cp);
+
+ /* And update */
+ mass += cp->mass;
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot += cp->e_pot;
mom[0] += cp->mom[0];
mom[1] += cp->mom[1];
mom[2] += cp->mom[2];
- ang[0] += cp->ang[0];
- ang[1] += cp->ang[1];
- ang[2] += cp->ang[2];
+ ang_mom[0] += cp->ang_mom[0];
+ ang_mom[1] += cp->ang_mom[1];
+ ang_mom[2] += cp->ang_mom[2];
}
}
}
/* Store the collected values in the cell. */
- c->ti_end_min = ti_end_min;
- c->ti_end_max = ti_end_max;
- c->updated = updated;
- c->g_updated = g_updated;
+ c->mass = mass;
c->e_kin = e_kin;
c->e_int = e_int;
c->e_pot = e_pot;
c->mom[0] = mom[0];
c->mom[1] = mom[1];
c->mom[2] = mom[2];
- c->ang[0] = ang[0];
- c->ang[1] = ang[1];
- c->ang[2] = ang[2];
+ c->ang_mom[0] = ang_mom[0];
+ c->ang_mom[1] = ang_mom[1];
+ c->ang_mom[2] = ang_mom[2];
+}
+
+void engine_print_stats(struct engine *e) {
+
+ const struct space *s = e->s;
+
+ double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
+ double mom[3] = {0.0, 0.0, 0.0}, ang_mom[3] = {0.0, 0.0, 0.0};
+
+ /* Collect the cell data. */
+ for (int k = 0; k < s->nr_cells; k++)
+ if (s->cells[k].nodeID == e->nodeID) {
+ struct cell *c = &s->cells[k];
+
+ /* Make the top-cells recurse */
+ engine_collect_drift(c);
+
+ /* And aggregate */
+ mass += c->mass;
+ e_kin += c->e_kin;
+ e_int += c->e_int;
+ e_pot += c->e_pot;
+ mom[0] += c->mom[0];
+ mom[1] += c->mom[1];
+ mom[2] += c->mom[2];
+ ang_mom[0] += c->ang_mom[0];
+ ang_mom[1] += c->ang_mom[1];
+ ang_mom[2] += c->ang_mom[2];
+ }
+
+/* Aggregate the data from the different nodes. */
+#ifdef WITH_MPI
+ {
+ double in[10] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
+ double out[10];
+ out[0] = e_kin;
+ out[1] = e_int;
+ out[2] = e_pot;
+ out[3] = mom[0];
+ out[4] = mom[1];
+ out[5] = mom[2];
+ out[6] = ang_mom[0];
+ out[7] = ang_mom[1];
+ out[8] = ang_mom[2];
+ out[9] = mass;
+ if (MPI_Allreduce(out, in, 10, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
+ MPI_SUCCESS)
+ error("Failed to aggregate stats.");
+ e_kin = out[0];
+ e_int = out[1];
+ e_pot = out[2];
+ mom[0] = out[3];
+ mom[1] = out[4];
+ mom[2] = out[5];
+ ang_mom[0] = out[6];
+ ang_mom[1] = out[7];
+ ang_mom[2] = out[8];
+ mass = out[9];
+ }
+#endif
+
+ const double e_tot = e_kin + e_int + e_pot;
+
+ /* Print info */
+ if (e->nodeID == 0) {
+ fprintf(e->file_stats,
+ " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e\n",
+ e->time, mass, e_tot, e_kin, e_int, e_pot, mom[0], mom[1], mom[2],
+ ang_mom[0], ang_mom[1], ang_mom[2]);
+ fflush(e->file_stats);
+ }
}
/**
@@ -1868,7 +2047,6 @@ void engine_collect_kick(struct cell *c) {
* @param mask The task mask to launch.
* @param submask The sub-task mask to launch.
*/
-
void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask) {
@@ -1913,18 +2091,18 @@ void engine_init_particles(struct engine *e) {
if (e->nodeID == 0) message("Initialising particles");
+ engine_prepare(e);
+
/* Make sure all particles are ready to go */
/* i.e. clean-up any stupid state in the ICs */
if (e->policy & engine_policy_hydro) {
- space_map_cells_pre(s, 1, cell_init_parts, NULL);
+ space_map_cells_pre(s, 0, cell_init_parts, NULL);
}
if ((e->policy & engine_policy_self_gravity) ||
(e->policy & engine_policy_external_gravity)) {
- space_map_cells_pre(s, 1, cell_init_gparts, NULL);
+ space_map_cells_pre(s, 0, cell_init_gparts, NULL);
}
- engine_prepare(e);
-
engine_marktasks(e);
/* Build the masks corresponding to the policy */
@@ -1987,13 +2165,7 @@ void engine_init_particles(struct engine *e) {
*/
void engine_step(struct engine *e) {
- int updates = 0, g_updates = 0;
- int ti_end_min = max_nr_timesteps, ti_end_max = 0;
- double e_pot = 0.0, e_int = 0.0, e_kin = 0.0;
- float mom[3] = {0.0, 0.0, 0.0};
- float ang[3] = {0.0, 0.0, 0.0};
double snapshot_drift_time = 0.;
- struct space *s = e->s;
TIMER_TIC2;
@@ -2002,66 +2174,11 @@ void engine_step(struct engine *e) {
e->tic_step = getticks();
- /* Collect the cell data. */
- for (int k = 0; k < s->nr_cells; k++)
- if (s->cells[k].nodeID == e->nodeID) {
- struct cell *c = &s->cells[k];
-
- /* Recurse */
- engine_collect_kick(c);
-
- /* And aggregate */
- ti_end_min = min(ti_end_min, c->ti_end_min);
- ti_end_max = max(ti_end_max, c->ti_end_max);
- e_kin += c->e_kin;
- e_int += c->e_int;
- e_pot += c->e_pot;
- updates += c->updated;
- g_updates += c->g_updated;
- mom[0] += c->mom[0];
- mom[1] += c->mom[1];
- mom[2] += c->mom[2];
- ang[0] += c->ang[0];
- ang[1] += c->ang[1];
- ang[2] += c->ang[2];
- }
-
-/* Aggregate the data from the different nodes. */
-#ifdef WITH_MPI
- {
- int in_i[1], out_i[1];
- in_i[0] = 0;
- out_i[0] = ti_end_min;
- if (MPI_Allreduce(out_i, in_i, 1, MPI_INT, MPI_MIN, MPI_COMM_WORLD) !=
- MPI_SUCCESS)
- error("Failed to aggregate t_end_min.");
- ti_end_min = in_i[0];
- out_i[0] = ti_end_max;
- if (MPI_Allreduce(out_i, in_i, 1, MPI_INT, MPI_MAX, MPI_COMM_WORLD) !=
- MPI_SUCCESS)
- error("Failed to aggregate t_end_max.");
- ti_end_max = in_i[0];
- }
- {
- double in_d[5], out_d[5];
- out_d[0] = updates;
- out_d[1] = g_updates;
- out_d[2] = e_kin;
- out_d[3] = e_int;
- out_d[4] = e_pot;
- if (MPI_Allreduce(out_d, in_d, 5, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
- MPI_SUCCESS)
- error("Failed to aggregate energies.");
- updates = in_d[0];
- g_updates = in_d[1];
- e_kin = in_d[2];
- e_int = in_d[3];
- e_pot = in_d[4];
- }
-#endif
+ /* Recover the (integer) end of the next time-step */
+ engine_collect_timestep(e);
/* Check for output */
- while (ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0) {
+ while (e->ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0) {
e->ti_old = e->ti_current;
e->ti_current = e->ti_nextSnapshot;
@@ -2082,7 +2199,7 @@ void engine_step(struct engine *e) {
/* Move forward in time */
e->ti_old = e->ti_current;
- e->ti_current = ti_end_min;
+ e->ti_current = e->ti_end_min;
e->step += 1;
e->time = e->ti_current * e->timeBase + e->timeBegin;
e->timeOld = e->ti_old * e->timeBase + e->timeBegin;
@@ -2094,15 +2211,15 @@ void engine_step(struct engine *e) {
if (e->nodeID == 0) {
/* Print some information to the screen */
- printf(" %6d %14e %14e %10d %10d %21.3f\n", e->step, e->time, e->timeStep,
- updates, g_updates, e->wallclock_time);
+ printf(" %6d %14e %14e %10zd %10zd %21.3f\n", e->step, e->time,
+ e->timeStep, e->updates, e->g_updates, e->wallclock_time);
fflush(stdout);
+ }
- /* Write some energy statistics */
- fprintf(e->file_stats, "%d %f %f %f %f %f %f %f %f %f %f %f\n", e->step,
- e->time, e_kin, e_int, e_pot, e_kin + e_int + e_pot, mom[0], mom[1],
- mom[2], ang[0], ang[1], ang[2]);
- fflush(e->file_stats);
+ /* Save some statistics */
+ if (e->time - e->timeLastStatistics >= e->deltaTimeStatistics) {
+ engine_print_stats(e);
+ e->timeLastStatistics += e->deltaTimeStatistics;
}
/* Re-distribute the particles amongst the nodes? */
@@ -2114,11 +2231,17 @@ void engine_step(struct engine *e) {
/* Build the masks corresponding to the policy */
unsigned int mask = 0, submask = 0;
- /* We always have sort tasks and kick tasks */
+ /* We always have sort tasks and init tasks */
mask |= 1 << task_type_sort;
- mask |= 1 << task_type_kick;
mask |= 1 << task_type_init;
+ /* Add the correct kick task */
+ if (e->policy & engine_policy_fixdt) {
+ mask |= 1 << task_type_kick_fixdt;
+ } else {
+ mask |= 1 << task_type_kick;
+ }
+
/* Add the tasks corresponding to hydro operations to the masks */
if (e->policy & engine_policy_hydro) {
@@ -2390,23 +2513,58 @@ void engine_dump_snapshot(struct engine *e) {
(float)clocks_diff(&time1, &time2), clocks_getunit());
}
-#if defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
-static bool hyperthreads_present(void) {
-#ifdef __linux__
- FILE *f =
- fopen("/sys/devices/system/cpu/cpu0/topology/thread_siblings_list", "r");
+#ifdef HAVE_SETAFFINITY
+/**
+ * @brief Returns the initial affinity the main thread is using.
+ */
+static cpu_set_t *engine_entry_affinity() {
+
+ static int use_entry_affinity = 0;
+
+ if (!use_entry_affinity) {
+ pthread_t engine = pthread_self();
+ pthread_getaffinity_np(engine, sizeof(entry_affinity), &entry_affinity);
+ use_entry_affinity = 1;
+ }
+
+ return &entry_affinity;
+}
+#endif
+
+/**
+ * @brief Ensure the NUMA node on which we initialise (first touch) everything
+ * doesn't change before engine_init allocates NUMA-local workers.
+ */
+void engine_pin() {
- int c;
- while ((c = fgetc(f)) != EOF && c != ',')
+#ifdef HAVE_SETAFFINITY
+ cpu_set_t *entry_affinity = engine_entry_affinity();
+ int pin;
+ for (pin = 0; pin < CPU_SETSIZE && !CPU_ISSET(pin, entry_affinity); ++pin)
;
- fclose(f);
- return c == ',';
+ cpu_set_t affinity;
+ CPU_ZERO(&affinity);
+ CPU_SET(pin, &affinity);
+ if (sched_setaffinity(0, sizeof(affinity), &affinity) != 0) {
+ error("failed to set engine's affinity");
+ }
#else
- return true; // just guess
+ error("SWIFT was not compiled with support for pinning.");
#endif
}
+
+/**
+ * @brief Unpins the main thread.
+ */
+void engine_unpin() {
+#ifdef HAVE_SETAFFINITY
+ pthread_t main_thread = pthread_self();
+ pthread_setaffinity_np(main_thread, sizeof(entry_affinity), &entry_affinity);
+#else
+ error("SWIFT was not compiled with support for pinning.");
#endif
+}
/**
* @brief init an engine with the given number of threads, queues, and
@@ -2418,6 +2576,7 @@ static bool hyperthreads_present(void) {
* @param nr_nodes The number of MPI ranks.
* @param nodeID The MPI rank of this node.
* @param nr_threads The number of threads per MPI rank.
+ * @param with_aff use processor affinity, if supported.
* @param policy The queuing policy to use.
* @param verbose Is this #engine talkative ?
* @param physical_constants The #phys_const used for this run.
@@ -2427,7 +2586,7 @@ static bool hyperthreads_present(void) {
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
- int nr_threads, int policy, int verbose,
+ int nr_threads, int with_aff, int policy, int verbose,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential) {
@@ -2456,6 +2615,7 @@ void engine_init(struct engine *e, struct space *s,
e->ti_current = 0;
e->timeStep = 0.;
e->timeBase = 0.;
+ e->timeBase_inv = 0.;
e->timeFirstSnapshot =
parser_get_param_double(params, "Snapshots:time_first");
e->deltaTimeSnapshot =
@@ -2467,12 +2627,16 @@ void engine_init(struct engine *e, struct space *s,
e->dt_min = parser_get_param_double(params, "TimeIntegration:dt_min");
e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
e->file_stats = NULL;
+ e->deltaTimeStatistics =
+ parser_get_param_double(params, "Statistics:delta_time");
+ e->timeLastStatistics = e->timeBegin - e->deltaTimeStatistics;
e->verbose = verbose;
e->count_step = 0;
e->wallclock_time = 0.f;
e->physical_constants = physical_constants;
e->hydro_properties = hydro;
e->external_potential = potential;
+ e->parameter_file = params;
engine_rank = nodeID;
/* Make the space link back to the engine. */
@@ -2483,74 +2647,117 @@ void engine_init(struct engine *e, struct space *s,
if (nr_queues <= 0) nr_queues = e->nr_threads;
s->nr_queues = nr_queues;
+/* Deal with affinity. For now, just figure out the number of cores. */
#if defined(HAVE_SETAFFINITY)
const int nr_cores = sysconf(_SC_NPROCESSORS_ONLN);
- int cpuid[nr_cores];
+ cpu_set_t *entry_affinity = engine_entry_affinity();
+ const int nr_affinity_cores = CPU_COUNT(entry_affinity);
+
+ if (nr_cores > CPU_SETSIZE) /* Unlikely, except on e.g. SGI UV. */
+ error("must allocate dynamic cpu_set_t (too many cores per node)");
+
+ char *buf = malloc((nr_cores + 1) * sizeof(char));
+ buf[nr_cores] = '\0';
+ for (int j = 0; j < nr_cores; ++j) {
+ /* Reversed bit order from convention, but same as e.g. Intel MPI's
+ * I_MPI_PIN_DOMAIN explicit mask: left-to-right, LSB-to-MSB. */
+ buf[j] = CPU_ISSET(j, entry_affinity) ? '1' : '0';
+ }
+
+ if (verbose && with_aff) message("Affinity at entry: %s", buf);
+
+ int *cpuid = malloc(nr_affinity_cores * sizeof(int));
cpu_set_t cpuset;
- if ((policy & engine_policy_cputight) == engine_policy_cputight) {
- for (int k = 0; k < nr_cores; k++) cpuid[k] = k;
- } else {
- /* Get next highest power of 2. */
- int maxint = 1;
- while (maxint < nr_cores) maxint *= 2;
- cpuid[0] = 0;
- int k = 1;
- for (int i = 1; i < maxint; i *= 2)
- for (int j = maxint / i / 2; j < maxint; j += maxint / i)
- if (j < nr_cores && j != 0) cpuid[k++] = j;
+ int skip = 0;
+ for (int k = 0; k < nr_affinity_cores; k++) {
+ int c;
+ for (c = skip; c < CPU_SETSIZE && !CPU_ISSET(c, entry_affinity); ++c)
+ ;
+ cpuid[k] = c;
+ skip = c + 1;
+ }
+
+ if (with_aff) {
#if defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
- /* Ascending NUMA distance. Bubblesort(!) for stable equidistant CPUs. */
- if (numa_available() >= 0) {
- if (nodeID == 0) message("prefer NUMA-local CPUs");
-
- const int home = numa_node_of_cpu(sched_getcpu());
- const int half = nr_cores / 2;
- const bool swap_hyperthreads = hyperthreads_present();
- bool done = false;
- if (swap_hyperthreads && nodeID == 0)
- message("prefer physical cores to hyperthreads");
-
- while (!done) {
- done = true;
- for (int i = 1; i < nr_cores; i++) {
- const int node_a = numa_node_of_cpu(cpuid[i - 1]);
- const int node_b = numa_node_of_cpu(cpuid[i]);
-
- /* Avoid using local hyperthreads over unused remote physical cores.
- * Assume two hyperthreads, and that cpuid >= half partitions them.
- */
- const int thread_a = swap_hyperthreads && cpuid[i - 1] >= half;
- const int thread_b = swap_hyperthreads && cpuid[i] >= half;
-
- bool swap = thread_a > thread_b;
- if (thread_a == thread_b)
- swap = numa_distance(home, node_a) > numa_distance(home, node_b);
-
- if (swap) {
- const int t = cpuid[i - 1];
- cpuid[i - 1] = cpuid[i];
- cpuid[i] = t;
- done = false;
+ if ((policy & engine_policy_cputight) != engine_policy_cputight) {
+
+ if (numa_available() >= 0) {
+ if (nodeID == 0) message("prefer NUMA-distant CPUs");
+
+ /* Get list of numa nodes of all available cores. */
+ int *nodes = malloc(nr_affinity_cores * sizeof(int));
+ int nnodes = 0;
+ for (int i = 0; i < nr_affinity_cores; i++) {
+ nodes[i] = numa_node_of_cpu(cpuid[i]);
+ if (nodes[i] > nnodes) nnodes = nodes[i];
+ }
+ nnodes += 1;
+
+ /* Count cores per node. */
+ int *core_counts = malloc(nnodes * sizeof(int));
+ for (int i = 0; i < nr_affinity_cores; i++) {
+ core_counts[nodes[i]] = 0;
+ }
+ for (int i = 0; i < nr_affinity_cores; i++) {
+ core_counts[nodes[i]] += 1;
+ }
+
+ /* Index cores within each node. */
+ int *core_indices = malloc(nr_affinity_cores * sizeof(int));
+ for (int i = nr_affinity_cores - 1; i >= 0; i--) {
+ core_indices[i] = core_counts[nodes[i]];
+ core_counts[nodes[i]] -= 1;
+ }
+
+ /* Now sort so that we pick adjacent cpuids from different nodes
+ * by sorting internal node core indices. */
+ int done = 0;
+ while (!done) {
+ done = 1;
+ for (int i = 1; i < nr_affinity_cores; i++) {
+ if (core_indices[i] < core_indices[i - 1]) {
+ int t = cpuid[i - 1];
+ cpuid[i - 1] = cpuid[i];
+ cpuid[i] = t;
+
+ t = core_indices[i - 1];
+ core_indices[i - 1] = core_indices[i];
+ core_indices[i] = t;
+ done = 0;
+ }
}
}
+
+ free(nodes);
+ free(core_counts);
+ free(core_indices);
}
}
#endif
+ } else {
+ if (nodeID == 0) message("no processor affinity used");
+
+ } /* with_aff */
- if (nodeID == 0) {
+ /* Avoid (unexpected) interference between engine and runner threads. We can
+ * do this once we've made at least one call to engine_entry_affinity and
+ * maybe numa_node_of_cpu(sched_getcpu()), even if the engine isn't already
+ * pinned. Also unpin this when asked to not pin at all (!with_aff). */
+ engine_unpin();
+#endif
+
+ if (with_aff) {
#ifdef WITH_MPI
- printf("[%04i] %s engine_init: cpu map is [ ", nodeID,
- clocks_get_timesincestart());
+ printf("[%04i] %s engine_init: cpu map is [ ", nodeID,
+ clocks_get_timesincestart());
#else
- printf("%s engine_init: cpu map is [ ", clocks_get_timesincestart());
+ printf("%s engine_init: cpu map is [ ", clocks_get_timesincestart());
#endif
- for (int i = 0; i < nr_cores; i++) printf("%i ", cpuid[i]);
- printf("].\n");
- }
+ for (int i = 0; i < nr_affinity_cores; i++) printf("%i ", cpuid[i]);
+ printf("].\n");
}
-#endif
/* Are we doing stuff in parallel? */
if (nr_nodes > 1) {
@@ -2570,8 +2777,10 @@ void engine_init(struct engine *e, struct space *s,
if (e->nodeID == 0) {
e->file_stats = fopen("energy.txt", "w");
fprintf(e->file_stats,
- "# Step Time E_kin E_int E_pot E_tot "
- "p_x p_y p_z ang_x ang_y ang_z\n");
+ "# %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n",
+ "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "p_x", "p_y",
+ "p_z", "ang_x", "ang_y", "ang_z");
+ fflush(e->file_stats);
}
/* Print policy */
@@ -2597,6 +2806,7 @@ void engine_init(struct engine *e, struct space *s,
/* Deal with timestep */
e->timeBase = (e->timeEnd - e->timeBegin) / max_nr_timesteps;
+ e->timeBase_inv = 1.0 / e->timeBase;
e->ti_current = 0;
/* Fixed time-step case */
@@ -2687,19 +2897,24 @@ void engine_init(struct engine *e, struct space *s,
if (pthread_create(&e->runners[k].thread, NULL, &runner_main,
&e->runners[k]) != 0)
error("Failed to create runner thread.");
- if (e->policy & engine_policy_setaffinity) {
+
+ /* Try to pin the runner to a given core */
+ if (with_aff &&
+ (e->policy & engine_policy_setaffinity) == engine_policy_setaffinity) {
#if defined(HAVE_SETAFFINITY)
/* Set a reasonable queue ID. */
- e->runners[k].cpuid = cpuid[k % nr_cores];
+ int coreid = k % nr_affinity_cores;
+ e->runners[k].cpuid = cpuid[coreid];
+
if (nr_queues < e->nr_threads)
- e->runners[k].qid = cpuid[k % nr_cores] * nr_queues / nr_cores;
+ e->runners[k].qid = cpuid[coreid] * nr_queues / nr_affinity_cores;
else
e->runners[k].qid = k;
/* Set the cpu mask to zero | e->id. */
CPU_ZERO(&cpuset);
- CPU_SET(cpuid[k % nr_cores], &cpuset);
+ CPU_SET(cpuid[coreid], &cpuset);
/* Apply this mask to the runner's pthread. */
if (pthread_setaffinity_np(e->runners[k].thread, sizeof(cpu_set_t),
@@ -2713,9 +2928,23 @@ void engine_init(struct engine *e, struct space *s,
e->runners[k].cpuid = k;
e->runners[k].qid = k * nr_queues / e->nr_threads;
}
- /* message( "runner %i on cpuid=%i with qid=%i." , e->runners[k].id , */
- /* e->runners[k].cpuid , e->runners[k].qid ); */
+ if (verbose) {
+ if (with_aff)
+ message("runner %i on cpuid=%i with qid=%i.", e->runners[k].id,
+ e->runners[k].cpuid, e->runners[k].qid);
+ else
+ message("runner %i using qid=%i no cpuid.", e->runners[k].id,
+ e->runners[k].qid);
+ }
+ }
+
+/* Free the affinity stuff */
+#if defined(HAVE_SETAFFINITY)
+ if (with_aff) {
+ free(cpuid);
+ free(buf);
}
+#endif
/* Wait for the runner threads to be in place. */
while (e->barrier_running || e->barrier_launch)
diff --git a/src/engine.h b/src/engine.h
index a830a37538d7356d0d5de7c67778e8ef7fcc1ed0..392cd9f25c975213340497cacc13451a4b5618b9 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -39,16 +39,15 @@
/* Includes. */
#include "hydro_properties.h"
#include "lock.h"
+#include "parser.h"
+#include "partition.h"
+#include "physical_constants.h"
+#include "potentials.h"
#include "proxy.h"
#include "runner.h"
#include "scheduler.h"
#include "space.h"
#include "task.h"
-#include "parser.h"
-#include "partition.h"
-#include "physical_constants.h"
-#include "potentials.h"
-#include "threadpool.h"
#include "units.h"
/* Some constants. */
@@ -111,6 +110,9 @@ struct engine {
/* The task scheduler. */
struct scheduler sched;
+ /* Common threadpool for all the engine's tasks. */
+ struct threadpool threadpool;
+
/* The minimum and maximum allowed dt */
double dt_min, dt_max;
@@ -133,6 +135,13 @@ struct engine {
/* Time base */
double timeBase;
+ double timeBase_inv;
+
+ /* Minimal ti_end for the next time-step */
+ int ti_end_min;
+
+ /* Number of particles updated */
+ size_t updates, g_updates;
/* Snapshot information */
double timeFirstSnapshot;
@@ -141,8 +150,10 @@ struct engine {
char snapshotBaseName[200];
struct UnitSystem *snapshotUnits;
- /* File for statistics */
+ /* Statistics information */
FILE *file_stats;
+ double timeLastStatistics;
+ double deltaTimeStatistics;
/* The current step number. */
int step;
@@ -194,8 +205,8 @@ struct engine {
/* Properties of external gravitational potential */
const struct external_potential *external_potential;
- /* Common threadpool for all the engine's tasks. */
- struct threadpool threadpool;
+ /* The (parsed) parameter file */
+ const struct swift_params *parameter_file;
};
/* Function prototypes. */
@@ -204,7 +215,7 @@ void engine_compute_next_snapshot_time(struct engine *e);
void engine_dump_snapshot(struct engine *e);
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
- int nr_threads, int policy, int verbose,
+ int nr_threads, int with_aff, int policy, int verbose,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential);
@@ -226,5 +237,7 @@ void engine_redistribute(struct engine *e);
struct link *engine_addlink(struct engine *e, struct link *l, struct task *t);
void engine_print_policy(struct engine *e);
int engine_is_done(struct engine *e);
+void engine_pin();
+void engine_unpin();
#endif /* SWIFT_ENGINE_H */
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 4dd8d5b7c7194ae072f90e7703907500c208b9b9..0f62511eced181fbf3b5b781a50314dd08f7c0ef 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -31,11 +31,10 @@
* @param phys_const The physical constants in internal units.
* @param gp Pointer to the g-particle data.
*/
-__attribute__((always_inline))
- INLINE static float gravity_compute_timestep_external(
- const struct external_potential* potential,
- const struct phys_const* const phys_const,
- const struct gpart* const gp) {
+__attribute__((always_inline)) INLINE static float
+gravity_compute_timestep_external(const struct external_potential* potential,
+ const struct phys_const* const phys_const,
+ const struct gpart* const gp) {
float dt = FLT_MAX;
@@ -55,10 +54,9 @@ __attribute__((always_inline))
* @param phys_const The physical constants in internal units.
* @param gp Pointer to the g-particle data.
*/
-__attribute__((always_inline))
- INLINE static float gravity_compute_timestep_self(
- const struct phys_const* const phys_const,
- const struct gpart* const gp) {
+__attribute__((always_inline)) INLINE static float
+gravity_compute_timestep_self(const struct phys_const* const phys_const,
+ const struct gpart* const gp) {
float dt = FLT_MAX;
@@ -73,8 +71,8 @@ __attribute__((always_inline))
*
* @param gp The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void gravity_first_init_gpart(struct gpart* gp) {}
+__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
+ struct gpart* gp) {}
/**
* @brief Prepares a g-particle for the gravity calculation
@@ -84,8 +82,8 @@ __attribute__((always_inline))
*
* @param gp The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void gravity_init_part(struct gpart* gp) {
+__attribute__((always_inline)) INLINE static void gravity_init_part(
+ struct gpart* gp) {
/* Zero the acceleration */
gp->a_grav[0] = 0.f;
@@ -100,8 +98,8 @@ __attribute__((always_inline))
*
* @param gp The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void gravity_end_force(struct gpart* gp) {}
+__attribute__((always_inline)) INLINE static void gravity_end_force(
+ struct gpart* gp) {}
/**
* @brief Computes the gravitational acceleration induced by external potentials
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
index 72ae5fd55f399448d78c6f2d1533c7f007f16dd1..21d775d703c8bd7000ad03386fe818124642a8f1 100644
--- a/src/gravity/Default/gravity_debug.h
+++ b/src/gravity/Default/gravity_debug.h
@@ -17,8 +17,8 @@
*
******************************************************************************/
-__attribute__((always_inline))
- INLINE static void gravity_debug_particle(struct gpart* p) {
+__attribute__((always_inline)) INLINE static void gravity_debug_particle(
+ struct gpart* p) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
diff --git a/src/hydro.h b/src/hydro.h
index aacbb6ac1d16b38133ee573ee2b7ad95918fc9e5..b2ae9d57c399ecea818e9f3dc7db238e01487a9a 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -19,20 +19,24 @@
#ifndef SWIFT_HYDRO_H
#define SWIFT_HYDRO_H
-#include "./const.h"
+/* Includes. */
+#include "const.h"
+#include "hydro_properties.h"
+#include "kernel_hydro.h"
+#include "part.h"
/* Import the right functions */
#if defined(MINIMAL_SPH)
-#include "./hydro/Minimal/hydro_iact.h"
#include "./hydro/Minimal/hydro.h"
+#include "./hydro/Minimal/hydro_iact.h"
#define SPH_IMPLEMENTATION "Minimal version of SPH (e.g. Price 2010)"
#elif defined(GADGET2_SPH)
-#include "./hydro/Gadget2/hydro_iact.h"
#include "./hydro/Gadget2/hydro.h"
+#include "./hydro/Gadget2/hydro_iact.h"
#define SPH_IMPLEMENTATION "Gadget-2 version of SPH (Springel 2005)"
#elif defined(DEFAULT_SPH)
-#include "./hydro/Default/hydro_iact.h"
#include "./hydro/Default/hydro.h"
+#include "./hydro/Default/hydro_iact.h"
#define SPH_IMPLEMENTATION "Default version of SPH"
#else
#error "Invalid choice of SPH variant"
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index b4ed2daf385691958d5849b0ac2d651ac6f98d34..4a6b1900374eb6b422e6fa7422901293b36fd5eb 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -27,7 +27,7 @@
*
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
- struct part* p, struct xpart* xp,
+ const struct part* p, const struct xpart* xp,
const struct hydro_props* hydro_properties) {
const float CFL_condition = hydro_properties->CFL_condition;
@@ -53,9 +53,8 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_first_init_part(struct part* p, struct xpart* xp) {
-}
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+ struct part* p, struct xpart* xp) {}
/**
* @brief Prepares a particle for the density calculation.
@@ -65,8 +64,8 @@ __attribute__((always_inline))
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_init_part(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+ struct part* p) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
@@ -86,8 +85,8 @@ __attribute__((always_inline))
* @param p The particle to act upon
* @param time The current time
*/
-__attribute__((always_inline))
- INLINE static void hydro_end_density(struct part* p, float time) {
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+ struct part* p, float time) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -171,8 +170,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_reset_acceleration(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+ struct part* p) {
/* Reset the acceleration. */
p->a_hydro[0] = 0.0f;
@@ -218,8 +217,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_end_force(struct part* p) {}
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+ struct part* p) {}
/**
* @brief Kick the additional variables
@@ -239,16 +238,17 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_convert_quantities(struct part* p) {}
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part* p) {}
/**
* @brief Returns the internal energy of a particle
*
* @param p The particle of interest
+ * @param dt Time since the last kick
*/
-__attribute__((always_inline))
- INLINE static float hydro_get_internal_energy(struct part* p) {
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part* p, float dt) {
return p->u;
}
diff --git a/src/hydro/Default/hydro_debug.h b/src/hydro/Default/hydro_debug.h
index 2e7c3d683aa18eee7a550ac89517e3bd01e42107..79ee392d46ca75a2c097bf045b2d82c9f3dc96c0 100644
--- a/src/hydro/Default/hydro_debug.h
+++ b/src/hydro/Default/hydro_debug.h
@@ -17,8 +17,8 @@
*
******************************************************************************/
-__attribute__((always_inline))
- INLINE static void hydro_debug_particle(struct part* p, struct xpart* xp) {
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+ struct part* p, struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index 4f85299b9d61b3a66389bac3527a63068ab96db9..4d651b61bd934388414d54fa813820111d26e682 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -268,11 +268,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
* @brief Density loop (non-symmetric vectorized version)
*/
-__attribute__((always_inline))
- INLINE static void runner_iact_nonsym_vec_density(float *R2, float *Dx,
- float *Hi, float *Hj,
- struct part **pi,
- struct part **pj) {
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
+ struct part **pi, struct part **pj) {
#ifdef VECTORIZE
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index 09a7f4c7caa6a2a9fcf77c8d83b0ffcab2ef8b54..de71963b05149fa0b2df6222424478a6ad9b1f44 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/src/hydro/Default/hydro_io.h
@@ -104,10 +104,6 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Kernel function description */
writeAttribute_s(h_grpsph, "Kernel", kernel_name);
- // writeAttribute_f(h_grpsph, "Kernel eta", const_eta_kernel);
- // writeAttribute_f(h_grpsph, "Weighted N_ngb", kernel_nwneigh);
- // writeAttribute_f(h_grpsph, "Delta N_ngb", const_delta_nwneigh);
- // writeAttribute_f(h_grpsph, "Hydro gamma", const_hydro_gamma);
/* Viscosity and thermal conduction */
writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
@@ -123,11 +119,6 @@ void writeSPHflavour(hid_t h_grpsph) {
writeAttribute_f(h_grpsph, "Viscosity decay length", const_viscosity_length);
/* Time integration properties */
- // writeAttribute_f(h_grpsph, "CFL parameter", const_cfl);
- // writeAttribute_f(h_grpsph, "Maximal ln(Delta h) change over dt",
- // const_ln_max_h_change);
- // writeAttribute_f(h_grpsph, "Maximal Delta h change over dt",
- // exp(const_ln_max_h_change));
writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
const_max_u_change);
}
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 9fa1aae9dbb13ac3aa7e2271f30ef57e91ae85ae..0973acb0fb46411c778f2551fbe91621825f0278 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -25,7 +25,7 @@
*
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
- struct part* p, struct xpart* xp,
+ const struct part* p, const struct xpart* xp,
const struct hydro_props* hydro_properties) {
const float CFL_condition = hydro_properties->CFL_condition;
@@ -46,9 +46,8 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_first_init_part(struct part* p, struct xpart* xp) {
-}
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+ struct part* p, struct xpart* xp) {}
/**
* @brief Prepares a particle for the density calculation.
@@ -58,8 +57,8 @@ __attribute__((always_inline))
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_init_part(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+ struct part* p) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
@@ -79,8 +78,8 @@ __attribute__((always_inline))
* @param p The particle to act upon
* @param ti_current The current time (on the integer timeline)
*/
-__attribute__((always_inline))
- INLINE static void hydro_end_density(struct part* p, int ti_current) {
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+ struct part* p, int ti_current) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -148,8 +147,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_reset_acceleration(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+ struct part* p) {
/* Reset the acceleration. */
p->a_hydro[0] = 0.0f;
@@ -193,8 +192,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_end_force(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+ struct part* p) {
p->entropy_dt *=
(const_hydro_gamma - 1.f) * powf(p->rho, -(const_hydro_gamma - 1.f));
@@ -230,8 +229,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_convert_quantities(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part* p) {
p->entropy = (const_hydro_gamma - 1.f) * p->entropy *
powf(p->rho, -(const_hydro_gamma - 1.f));
@@ -241,10 +240,13 @@ __attribute__((always_inline))
* @brief Returns the internal energy of a particle
*
* @param p The particle of interest
+ * @param dt Time since the last kick
*/
-__attribute__((always_inline))
- INLINE static float hydro_get_internal_energy(struct part* p) {
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part* p, float dt) {
+
+ const float entropy = p->entropy + p->entropy_dt * dt;
- return p->entropy * powf(p->rho, const_hydro_gamma - 1.f) *
+ return entropy * powf(p->rho, const_hydro_gamma - 1.f) *
(1.f / (const_hydro_gamma - 1.f));
}
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index 4f00d8e020f9eaac30fe7f272ca1248d7e4eec58..31e89d438bc96c0f0f2ba56249664d28036cb607 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -17,8 +17,8 @@
*
******************************************************************************/
-__attribute__((always_inline))
- INLINE static void hydro_debug_particle(struct part* p, struct xpart* xp) {
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+ struct part* p, struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 41f3348db10172c1d0d5813debb50365252a413b..b977f25386fab0787c925635297a48fa85a8df24 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -83,8 +83,8 @@ __attribute__((always_inline)) INLINE static void hydro_write_particles(
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "SmoothingLength",
FLOAT, N, 1, parts, N_total, mpi_rank, offset, h, us,
UNIT_CONV_LENGTH);
- writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Entropy",
- FLOAT, N, 1, parts, N_total, mpi_rank, offset, entropy, us,
+ writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "Entropy", FLOAT, N,
+ 1, parts, N_total, mpi_rank, offset, entropy, us,
UNIT_CONV_ENTROPY_PER_UNIT_MASS);
writeArray(h_grp, fileName, xmfFile, partTypeGroupName, "ParticleIDs",
ULONGLONG, N, 1, parts, N_total, mpi_rank, offset, id, us,
@@ -104,10 +104,6 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Kernel function description */
writeAttribute_s(h_grpsph, "Kernel", kernel_name);
- // writeAttribute_f(h_grpsph, "Kernel eta", const_eta_kernel);
- // writeAttribute_f(h_grpsph, "Weighted N_ngb", kernel_nwneigh);
- // writeAttribute_f(h_grpsph, "Delta N_ngb", const_delta_nwneigh);
- // writeAttribute_f(h_grpsph, "Hydro gamma", const_hydro_gamma);
/* Viscosity and thermal conduction */
writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
@@ -116,11 +112,4 @@ void writeSPHflavour(hid_t h_grpsph) {
"Legacy Gadget-2 as in Springel (2005)");
writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
-
- /* Time integration properties */
- // writeAttribute_f(h_grpsph, "CFL parameter", const_cfl);
- // writeAttribute_f(h_grpsph, "Maximal ln(Delta h) change over dt",
- // const_ln_max_h_change);
- // writeAttribute_f(h_grpsph, "Maximal Delta h change over dt",
- // exp(const_ln_max_h_change));
}
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index f3553a009c22a2f6353796e8a278f0db7d66d294..f69dc3f1798f014e895c4a63760805b1739cec94 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -39,17 +39,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_first_init_part(struct part* p, struct xpart* xp) {
-}
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+ struct part* p, struct xpart* xp) {}
/**
* @brief Prepares a particle for the volume calculation.
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_init_part(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+ struct part* p) {
#ifdef SPH_GRADIENTS
/* use the old volumes to estimate new primitive variables to be used for the
@@ -127,8 +126,8 @@ __attribute__((always_inline))
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_end_volume(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_end_volume(
+ struct part* p) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -387,8 +386,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_end_gradient(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_end_gradient(
+ struct part* p) {
#ifndef SPH_GRADIENTS
float h, ih, ih2, ih3;
@@ -531,8 +530,8 @@ __attribute__((always_inline))
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_prepare_fluxes(struct part* p, struct xpart* xp) {
+__attribute__((always_inline)) INLINE static void hydro_prepare_fluxes(
+ struct part* p, struct xpart* xp) {
/* initialize variables used for timestep calculation */
p->timestepvars.vmax = 0.0f;
@@ -546,8 +545,8 @@ __attribute__((always_inline))
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_reset_acceleration(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+ struct part* p) {
/* figure out what to put here */
}
@@ -559,8 +558,8 @@ __attribute__((always_inline))
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_end_fluxes(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_end_fluxes(
+ struct part* p) {
/* do nothing */
}
@@ -572,8 +571,8 @@ __attribute__((always_inline))
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_convert_quantities(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part* p) {
float volume;
GFLOAT m;
@@ -605,17 +604,17 @@ __attribute__((always_inline))
}
// MATTHIEU
-__attribute__((always_inline))
- INLINE static void hydro_end_density(struct part* p, float time) {}
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+ struct part* p, float time) {}
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* p, struct xpart* xp, int ti_current, double timeBase) {}
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {}
-__attribute__((always_inline))
- INLINE static void hydro_end_force(struct part* p) {}
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+ struct part* p) {}
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt, float half_dt) {}
-__attribute__((always_inline))
- INLINE static float hydro_get_internal_energy(struct part* p) {
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ struct part* p) {
return 0.f;
}
diff --git a/src/hydro/Gizmo/hydro_debug.h b/src/hydro/Gizmo/hydro_debug.h
index 2cc957ed883436ce57e9d53d00a073693c9495df..365d85a2f651cf98b0713e8d82f11ae70fa9beaa 100644
--- a/src/hydro/Gizmo/hydro_debug.h
+++ b/src/hydro/Gizmo/hydro_debug.h
@@ -17,8 +17,8 @@
*
******************************************************************************/
-__attribute__((always_inline))
- INLINE static void hydro_debug_particle(struct part* p, struct xpart* xp) {
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+ struct part* p, struct xpart* xp) {
printf(
"x=[%.16e,%.16e,%.16e], "
"v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], volume=%.3e\n",
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 4fe875d3d07315051ef8b3051665a9ea0ef261b8..30a8d6cbebc851b44a5ee2339950aec9e15057c0 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -194,11 +194,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_hydro_loop1(
}
/* this corresponds to task_subtype_hydro_loop1 */
-__attribute__((always_inline))
- INLINE static void runner_iact_nonsym_hydro_loop1(float r2, float *dx,
- float hi, float hj,
- struct part *pi,
- struct part *pj) {
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_hydro_loop1(float r2, float *dx, float hi, float hj,
+ struct part *pi, struct part *pj) {
float r;
float xi;
@@ -487,11 +485,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_hydro_loop2(
#endif
}
-__attribute__((always_inline))
- INLINE static void runner_iact_nonsym_hydro_loop2(float r2, float *dx,
- float hi, float hj,
- struct part *pi,
- struct part *pj) {
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_hydro_loop2(float r2, float *dx, float hi, float hj,
+ struct part *pi, struct part *pj) {
#ifndef SPH_GRADIENTS
@@ -1025,11 +1021,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_hydro_loop3(
}
/* this corresponds to task_subtype_fluxes */
-__attribute__((always_inline))
- INLINE static void runner_iact_nonsym_hydro_loop3(float r2, float *dx,
- float hi, float hj,
- struct part *pi,
- struct part *pj) {
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_hydro_loop3(float r2, float *dx, float hi, float hj,
+ struct part *pi, struct part *pj) {
runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0);
}
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 01957065cf5b1d89bb839449ac60b9ec965c2c47..4222daafe82e7dc977cd87f57a5b9a235d505f00 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -31,7 +31,7 @@
*
*/
__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
- struct part* p, struct xpart* xp,
+ const struct part* p, const struct xpart* xp,
const struct hydro_props* hydro_properties) {
const float CFL_condition = hydro_properties->CFL_condition;
@@ -53,8 +53,8 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
* @param p The particle to act upon
* @param xp The extended particle data to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_first_init_part(struct part* p, struct xpart* xp) {
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+ struct part* p, struct xpart* xp) {
xp->u_full = p->u;
}
@@ -68,8 +68,8 @@ __attribute__((always_inline))
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_init_part(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+ struct part* p) {
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
@@ -88,8 +88,8 @@ __attribute__((always_inline))
* @param p The particle to act upon
* @param time The current time
*/
-__attribute__((always_inline))
- INLINE static void hydro_end_density(struct part* p, float time) {
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+ struct part* p, float time) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -143,8 +143,8 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_reset_acceleration(struct part* p) {
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+ struct part* p) {
/* Reset the acceleration. */
p->a_hydro[0] = 0.0f;
@@ -172,14 +172,7 @@ __attribute__((always_inline))
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, int t0, int t1, double timeBase) {
- const float dt = t1 - t0;
-
- /* Predict internal energy */
- const float w = p->u_dt / p->u * dt;
- if (fabsf(w) < 0.2f)
- p->u *= approx_expf(w); /* 4th order expansion of exp(w) */
- else
- p->u *= expf(w);
+ p->u = xp->u_full;
/* Need to recompute the pressure as well */
p->force.pressure = p->rho * p->u * (const_hydro_gamma - 1.f);
@@ -194,8 +187,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_end_force(struct part* p) {}
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+ struct part* p) {}
/**
* @brief Kick the additional variables
@@ -228,8 +221,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*
* @param p The particle to act upon
*/
-__attribute__((always_inline))
- INLINE static void hydro_convert_quantities(struct part* p) {}
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part* p) {}
/**
* @brief Returns the internal energy of a particle
@@ -239,9 +232,10 @@ __attribute__((always_inline))
* energy from the thermodynamic variable.
*
* @param p The particle of interest
+ * @param dt Time since the last kick
*/
-__attribute__((always_inline))
- INLINE static float hydro_get_internal_energy(struct part* p) {
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part* p, float dt) {
return p->u;
}
diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
index 85fdbac0a7893d575e52629600a5407c1bf77fcc..127ba75e99418b6a5dc197a44ccdc77de3cdef15 100644
--- a/src/hydro/Minimal/hydro_debug.h
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -17,8 +17,8 @@
*
******************************************************************************/
-__attribute__((always_inline))
- INLINE static void hydro_debug_particle(struct part* p, struct xpart* xp) {
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+ struct part* p, struct xpart* xp) {
printf(
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], "
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 3427ec538613842f8fbcf0d8ba5f9ba5a0b8d540..f453d8fa1c495472af27ccf98356a3dab0894e98 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -131,18 +131,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute dv dot r. */
- float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2];
- dvdr *= r_inv;
-
- /* Compute the relative velocity. (This is 0 if the particles move away from
- * each other and negative otherwise) */
- const float omega_ij = fminf(dvdr, 0.f);
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
/* Compute sound speeds */
const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
- const float v_sig = ci + cj + 3.f * omega_ij;
+ const float v_sig = ci + cj;
/* SPH acceleration term */
const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
@@ -157,8 +153,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->a_hydro[2] += mi * sph_term * dx[2];
/* Get the time derivative for u. */
- pi->u_dt += P_over_rho_i * mj * dvdr * wi_dr;
- pj->u_dt += P_over_rho_j * mi * dvdr * wj_dr;
+ pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
+ pj->u_dt += P_over_rho_j * mi * dvdr * r_inv * wj_dr;
/* Get the time derivative for h. */
pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
@@ -208,18 +204,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute dv dot r. */
- float dvdr = (pi->v[0] - pj->v[0]) * dx[0] + (pi->v[1] - pj->v[1]) * dx[1] +
- (pi->v[2] - pj->v[2]) * dx[2];
- dvdr *= r_inv;
-
- /* Compute the relative velocity. (This is 0 if the particles move away from
- * each other and negative otherwise) */
- const float omega_ij = fminf(dvdr, 0.f);
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
/* Compute sound speeds */
const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
- const float v_sig = ci + cj + 3.f * omega_ij;
+ const float v_sig = ci + cj;
/* SPH acceleration term */
const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
@@ -230,7 +222,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->a_hydro[2] -= mj * sph_term * dx[2];
/* Get the time derivative for u. */
- pi->u_dt += P_over_rho_i * mj * dvdr * wi_dr;
+ pi->u_dt += P_over_rho_i * mj * dvdr * r_inv * wi_dr;
/* Get the time derivative for h. */
pi->h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 43926d526b47fa8f1b92bb862edabefb7417de3c..1bbfe1065358017e0b27e2131ff717848221fe9c 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -104,9 +104,6 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Kernel function description */
writeAttribute_s(h_grpsph, "Kernel", kernel_name);
- // writeAttribute_f(h_grpsph, "Kernel eta", const_eta_kernel);
- // writeAttribute_f(h_grpsph, "Weighted N_ngb", kernel_nwneigh);
- // writeAttribute_f(h_grpsph, "Delta N_ngb", const_delta_nwneigh);
writeAttribute_f(h_grpsph, "Hydro gamma", const_hydro_gamma);
/* Viscosity and thermal conduction */
@@ -115,11 +112,6 @@ void writeSPHflavour(hid_t h_grpsph) {
writeAttribute_s(h_grpsph, "Viscosity Model", "No model");
/* Time integration properties */
- // writeAttribute_f(h_grpsph, "CFL parameter", const_cfl);
- // writeAttribute_f(h_grpsph, "Maximal ln(Delta h) change over dt",
- // const_ln_max_h_change);
- // writeAttribute_f(h_grpsph, "Maximal Delta h change over dt",
- // exp(const_ln_max_h_change));
writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
const_max_u_change);
}
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 4bedad7e0ccbca7656694c070c6140d1163db9cc..b4e37d672408ad3c0adf1f37fcd4fcbd50095d92 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -26,8 +26,8 @@
/* Local headers. */
#include "error.h"
-#include "kernel_hydro.h"
#include "hydro.h"
+#include "kernel_hydro.h"
void hydro_props_init(struct hydro_props *p,
const struct swift_params *params) {
@@ -55,7 +55,7 @@ void hydro_props_print(const struct hydro_props *p) {
message("Hydrodynamic kernel: %s with %.2f +/- %.2f neighbours (eta=%f).",
kernel_name, p->target_neighbours, p->delta_neighbours,
p->eta_neighbours);
- message("Hydrodynamic integration: CFL parmeter: %.4f.", p->CFL_condition);
+ message("Hydrodynamic integration: CFL parameter: %.4f.", p->CFL_condition);
message(
"Hydrodynamic integration: Max change of volume: %.2f "
"(max|dlog(h)/dt|=%f).",
diff --git a/src/intrinsics.h b/src/intrinsics.h
index 21b8c8e68bc45d8799db496ff30ac0cfb289acea..27e0bcc729b58493aed8c7eae7dfcdfc8f0855aa 100644
--- a/src/intrinsics.h
+++ b/src/intrinsics.h
@@ -26,8 +26,8 @@
* This is a wrapper for the GCC intrinsics with an implementation (from
* Hacker's Delight) if the compiler intrinsics are not available.
*/
-__attribute__((always_inline))
- INLINE static int intrinsics_clz(unsigned int x) {
+__attribute__((always_inline)) INLINE static int intrinsics_clz(
+ unsigned int x) {
#ifdef __GNUC__
/* Use GCC intrinsics if possible */
@@ -66,8 +66,8 @@ __attribute__((always_inline))
* This is a wrapper for the GCC intrinsics with an implementation (from
* Hacker's Delight) if the compiler intrinsics are not available.
*/
-__attribute__((always_inline))
- INLINE static int intrinsics_popcount(unsigned int x) {
+__attribute__((always_inline)) INLINE static int intrinsics_popcount(
+ unsigned int x) {
#ifdef __GNUC__
/* Use GCC intrinsics if possible */
diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h
index cca449efa6dde73eac428d1e20fa5c28156360fa..fedc046eed3bc64ed09120537cf6863179b171a6 100644
--- a/src/kernel_gravity.h
+++ b/src/kernel_gravity.h
@@ -26,11 +26,11 @@
#include "vector.h"
#define const_iepsilon (1. / const_epsilon)
-#define const_iepsilon2 (const_iepsilon *const_iepsilon)
-#define const_iepsilon3 (const_iepsilon2 *const_iepsilon)
-#define const_iepsilon4 (const_iepsilon2 *const_iepsilon2)
-#define const_iepsilon5 (const_iepsilon3 *const_iepsilon2)
-#define const_iepsilon6 (const_iepsilon3 *const_iepsilon3)
+#define const_iepsilon2 (const_iepsilon * const_iepsilon)
+#define const_iepsilon3 (const_iepsilon2 * const_iepsilon)
+#define const_iepsilon4 (const_iepsilon2 * const_iepsilon2)
+#define const_iepsilon5 (const_iepsilon3 * const_iepsilon2)
+#define const_iepsilon6 (const_iepsilon3 * const_iepsilon3)
/* The gravity kernel is defined as a degree 6 polynomial in the distance
r. The resulting value should be post-multiplied with r^-3, resulting
@@ -42,18 +42,28 @@
#define kernel_grav_degree 6
#define kernel_grav_ivals 2
#define kernel_grav_scale (2 * const_iepsilon)
-static float kernel_grav_coeffs
- [(kernel_grav_degree + 1) * (kernel_grav_ivals + 1)] = {
- 32.0f * const_iepsilon6, -192.0f / 5.0f * const_iepsilon5,
- 0.0f, 32.0f / 3.0f * const_iepsilon3,
- 0.0f, 0.0f,
- 0.0f, -32.0f / 3.0f * const_iepsilon6,
- 192.0f / 5.0f * const_iepsilon5, -48.0f * const_iepsilon4,
- 64.0f / 3.0f * const_iepsilon3, 0.0f,
- 0.0f, -1.0f / 15.0f,
- 0.0f, 0.0f,
- 0.0f, 0.0f,
- 0.0f, 0.0f,
+static float
+ kernel_grav_coeffs[(kernel_grav_degree + 1) * (kernel_grav_ivals + 1)] = {
+ 32.0f * const_iepsilon6,
+ -192.0f / 5.0f * const_iepsilon5,
+ 0.0f,
+ 32.0f / 3.0f * const_iepsilon3,
+ 0.0f,
+ 0.0f,
+ 0.0f,
+ -32.0f / 3.0f * const_iepsilon6,
+ 192.0f / 5.0f * const_iepsilon5,
+ -48.0f * const_iepsilon4,
+ 64.0f / 3.0f * const_iepsilon3,
+ 0.0f,
+ 0.0f,
+ -1.0f / 15.0f,
+ 0.0f,
+ 0.0f,
+ 0.0f,
+ 0.0f,
+ 0.0f,
+ 0.0f,
1.0f};
/**
@@ -76,8 +86,8 @@ __attribute__((always_inline)) INLINE static void kernel_grav_eval(float x,
* version).
*/
-__attribute__((always_inline))
- INLINE static void kernel_grav_eval_vec(vector *x, vector *w) {
+__attribute__((always_inline)) INLINE static void kernel_grav_eval_vec(
+ vector *x, vector *w) {
vector ind, c[kernel_grav_degree + 1];
int j, k;
@@ -179,8 +189,8 @@ __attribute__((always_inline)) INLINE static void blender_eval_vec(vector *x,
* distance x (Vectorized version). Gives a sensible answer only if x<2.
*/
-__attribute__((always_inline))
- INLINE static void blender_deval_vec(vector *x, vector *w, vector *dw_dx) {
+__attribute__((always_inline)) INLINE static void blender_deval_vec(
+ vector *x, vector *w, vector *dw_dx) {
vector ind, c[blender_degree + 1];
int j, k;
diff --git a/src/kernel_hydro.h b/src/kernel_hydro.h
index 06f950b0e2f3aec68055ba7976b32dc00980589e..b1774d8f35b7eddb6c2fdb0c341fa6299de74582 100644
--- a/src/kernel_hydro.h
+++ b/src/kernel_hydro.h
@@ -137,22 +137,22 @@ static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
/* Ok, now comes the real deal. */
/* First some powers of gamma = H/h */
-#define kernel_gamma2 ((float)(kernel_gamma *kernel_gamma))
-#define kernel_gamma3 ((float)(kernel_gamma *kernel_gamma *kernel_gamma))
+#define kernel_gamma2 ((float)(kernel_gamma * kernel_gamma))
+#define kernel_gamma3 ((float)(kernel_gamma * kernel_gamma * kernel_gamma))
#define kernel_gamma4 \
- ((float)(kernel_gamma *kernel_gamma *kernel_gamma *kernel_gamma))
+ ((float)(kernel_gamma * kernel_gamma * kernel_gamma * kernel_gamma))
#define kernel_igamma ((float)(1. / kernel_gamma))
-#define kernel_igamma2 ((float)(kernel_igamma *kernel_igamma))
-#define kernel_igamma3 ((float)(kernel_igamma *kernel_igamma *kernel_igamma))
+#define kernel_igamma2 ((float)(kernel_igamma * kernel_igamma))
+#define kernel_igamma3 ((float)(kernel_igamma * kernel_igamma * kernel_igamma))
#define kernel_igamma4 \
- ((float)(kernel_igamma *kernel_igamma *kernel_igamma *kernel_igamma))
+ ((float)(kernel_igamma * kernel_igamma * kernel_igamma * kernel_igamma))
/* The number of branches */
#define kernel_ivals_f ((float)(kernel_ivals))
/* Kernel self contribution (i.e. W(0,h)) */
#define kernel_root \
- ((float)(kernel_coeffs[kernel_degree]) * kernel_constant *kernel_igamma3)
+ ((float)(kernel_coeffs[kernel_degree]) * kernel_constant * kernel_igamma3)
/**
* @brief Computes the kernel function and its derivative.
@@ -200,8 +200,8 @@ __attribute__((always_inline)) INLINE static void kernel_deval(
* @param u The ratio of the distance to the smoothing length $u = x/h$.
* @param W (return) The value of the kernel function $W(x,h)$.
*/
-__attribute__((always_inline))
- INLINE static void kernel_eval(float u, float *restrict W) {
+__attribute__((always_inline)) INLINE static void kernel_eval(
+ float u, float *restrict W) {
/* Go to the range [0,1[ from [0,H[ */
const float x = u * kernel_igamma;
@@ -246,8 +246,8 @@ static const vector kernel_igamma4_vec = FILL_VEC((float)kernel_igamma4);
* @param w (return) The value of the kernel function $W(x,h)$.
* @param dw_dx (return) The norm of the gradient of $|\\nabla W(x,h)|$.
*/
-__attribute__((always_inline))
- INLINE static void kernel_deval_vec(vector *u, vector *w, vector *dw_dx) {
+__attribute__((always_inline)) INLINE static void kernel_deval_vec(
+ vector *u, vector *w, vector *dw_dx) {
/* Go to the range [0,1[ from [0,H[ */
vector x;
diff --git a/src/kick.h b/src/kick.h
new file mode 100644
index 0000000000000000000000000000000000000000..df3bd4cc6ce9819d0b65640db51d2ed20a7d59fe
--- /dev/null
+++ b/src/kick.h
@@ -0,0 +1,114 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_KICK_H
+#define SWIFT_KICK_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+#include "debug.h"
+
+/**
+ * @brief Perform the 'kick' operation on a #gpart
+ *
+ * @param gp The #gpart to kick.
+ * @param new_dti The (integer) time-step for this kick.
+ * @param timeBase The minimal allowed time-step size.
+ */
+__attribute__((always_inline)) INLINE static void kick_gpart(struct gpart* gp,
+ int new_dti,
+ double timeBase) {
+
+ /* Compute the time step for this kick */
+ const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
+ const int ti_end = gp->ti_end + new_dti / 2;
+ const double dt = (ti_end - ti_start) * timeBase;
+ const double half_dt = (ti_end - gp->ti_end) * timeBase;
+
+ /* Move particle forward in time */
+ gp->ti_begin = gp->ti_end;
+ gp->ti_end = gp->ti_begin + new_dti;
+
+ /* Kick particles in momentum space */
+ gp->v_full[0] += gp->a_grav[0] * dt;
+ gp->v_full[1] += gp->a_grav[1] * dt;
+ gp->v_full[2] += gp->a_grav[2] * dt;
+
+ /* Extra kick work */
+ gravity_kick_extra(gp, dt, half_dt);
+}
+
+/**
+ * @brief Perform the 'kick' operation on a #part
+ *
+ * @param p The #part to kick.
+ * @param xp The #xpart of the particle.
+ * @param new_dti The (integer) time-step for this kick.
+ * @param timeBase The minimal allowed time-step size.
+ */
+__attribute__((always_inline)) INLINE static void kick_part(struct part* p,
+ struct xpart* xp,
+ int new_dti,
+ double timeBase) {
+
+ /* Compute the time step for this kick */
+ const int ti_start = (p->ti_begin + p->ti_end) / 2;
+ const int ti_end = p->ti_end + new_dti / 2;
+ const double dt = (ti_end - ti_start) * timeBase;
+ const double half_dt = (ti_end - p->ti_end) * timeBase;
+
+ /* Move particle forward in time */
+ p->ti_begin = p->ti_end;
+ p->ti_end = p->ti_begin + new_dti;
+ if (p->gpart != NULL) {
+ p->gpart->ti_begin = p->ti_begin;
+ p->gpart->ti_end = p->ti_end;
+ }
+
+ /* Get the acceleration */
+ float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]};
+ if (p->gpart != NULL) {
+ a_tot[0] += p->gpart->a_grav[0];
+ a_tot[1] += p->gpart->a_grav[1];
+ a_tot[1] += p->gpart->a_grav[2];
+ }
+
+ /* Kick particles in momentum space */
+ xp->v_full[0] += a_tot[0] * dt;
+ xp->v_full[1] += a_tot[1] * dt;
+ xp->v_full[2] += a_tot[2] * dt;
+ if (p->gpart != NULL) {
+ p->gpart->v_full[0] = xp->v_full[0];
+ p->gpart->v_full[1] = xp->v_full[1];
+ p->gpart->v_full[2] = xp->v_full[2];
+ }
+
+ /* Go back by half-step for the hydro velocity */
+ p->v[0] = xp->v_full[0] - half_dt * a_tot[0];
+ p->v[1] = xp->v_full[1] - half_dt * a_tot[1];
+ p->v[2] = xp->v_full[2] - half_dt * a_tot[2];
+
+ /* Extra kick work */
+ hydro_kick_extra(p, xp, dt, half_dt);
+ if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt, half_dt);
+}
+
+#endif /* SWIFT_KICK_H */
diff --git a/src/lock.h b/src/lock.h
index 84bfa7bec77a5f23a69a7288c460a94a6857e20a..31f80605efb29bfaaadc65e25ba4f31679f8fbe8 100644
--- a/src/lock.h
+++ b/src/lock.h
@@ -27,7 +27,7 @@
#ifdef PTHREAD_SPINLOCK
#include <pthread.h>
-#define lock_type pthread_spinlock_t
+#define swift_lock_type pthread_spinlock_t
#define lock_init(l) (pthread_spin_init(l, PTHREAD_PROCESS_PRIVATE) != 0)
#define lock_destroy(l) (pthread_spin_destroy(l) != 0)
#define lock_lock(l) (pthread_spin_lock(l) != 0)
@@ -37,7 +37,7 @@
#elif defined(PTHREAD_LOCK)
#include <pthread.h>
-#define lock_type pthread_mutex_t
+#define swift_lock_type pthread_mutex_t
#define lock_init(l) (pthread_mutex_init(l, NULL) != 0)
#define lock_destroy(l) (pthread_mutex_destroy(l) != 0)
#define lock_lock(l) (pthread_mutex_lock(l) != 0)
@@ -46,7 +46,7 @@
#define lock_unlock_blind(l) pthread_mutex_unlock(l)
#else
-#define lock_type volatile int
+#define swift_lock_type volatile int
#define lock_init(l) (*(l) = 0)
#define lock_destroy(l) 0
INLINE static int lock_lock(volatile int *l) {
diff --git a/src/map.c b/src/map.c
index da13fbfb4ac00ed58184f7fe818826c82265a1de..fbe57fde7b1e29c49b0f27d86d177245dd9a27e2 100644
--- a/src/map.c
+++ b/src/map.c
@@ -18,10 +18,10 @@
*
******************************************************************************/
+#include "map.h"
#include <stdio.h>
#include <stdlib.h>
#include "error.h"
-#include "map.h"
/**
* @brief Mapping function to draw a specific cell (gnuplot).
diff --git a/src/map.h b/src/map.h
index 0753c2641af6deb050c1dcef6bcd3ae4621ae6aa..950a5fd96ebdc7177b41912b1565163f33de8701 100644
--- a/src/map.h
+++ b/src/map.h
@@ -22,8 +22,8 @@
#ifndef SWIFT_MAP_H
#define SWIFT_MAP_H
-#include "part.h"
#include "cell.h"
+#include "part.h"
void map_cells_plot(struct cell *c, void *data);
void map_check(struct part *p, struct cell *c, void *data);
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 4579be8f04ae687140413383e061988c9783b570..e779b56a85da3db72ea860fb336fa42037fbcd0e 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -529,7 +529,7 @@ void read_ic_parallel(char* fileName, double dim[3], struct part** parts,
void write_output_parallel(struct engine* e, const char* baseName,
struct UnitSystem* us, int mpi_rank, int mpi_size,
MPI_Comm comm, MPI_Info info) {
- hid_t h_file = 0, h_grp = 0, h_grpsph = 0;
+ hid_t h_file = 0, h_grp = 0;
const size_t Ngas = e->s->nr_parts;
const size_t Ntot = e->s->nr_gparts;
int periodic = e->s->periodic;
@@ -633,10 +633,17 @@ void write_output_parallel(struct engine* e, const char* baseName,
writeCodeDescription(h_file);
/* Print the SPH parameters */
- h_grpsph = H5Gcreate(h_file, "/SPH", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if (h_grpsph < 0) error("Error while creating SPH group");
- writeSPHflavour(h_grpsph);
- H5Gclose(h_grpsph);
+ h_grp = H5Gcreate(h_file, "/SPH", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating SPH group");
+ writeSPHflavour(h_grp);
+ H5Gclose(h_grp);
+
+ /* Print the runtime parameters */
+ h_grp =
+ H5Gcreate(h_file, "/Parameters", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating parameters group");
+ parser_write_params_to_hdf5(e->parameter_file, h_grp);
+ H5Gclose(h_grp);
/* Print the system of Units */
writeUnitSystem(h_file, us);
diff --git a/src/parser.c b/src/parser.c
index 49b035db05722d13d1ce48129b8fae8bfc1ad86b..74c277f036c733b5b1fbff4f2cb477b52169e2c6 100644
--- a/src/parser.c
+++ b/src/parser.c
@@ -22,14 +22,15 @@
/* Some standard headers. */
/* Needs to be included so that strtok returns char * instead of a int *. */
-#include <string.h>
-#include <stdlib.h>
#include <ctype.h>
+#include <stdlib.h>
+#include <string.h>
/* This object's header. */
#include "parser.h"
/* Local headers. */
+#include "common_io.h"
#include "error.h"
#define PARSER_COMMENT_STRING "#"
@@ -255,7 +256,7 @@ static void parse_value(char *line, struct swift_params *params) {
/* Check that it is a parameter inside a section.*/
if (*line == ' ' || *line == '\t') {
parse_section_param(line, &isFirstParam, section, params);
- } else {/*Else it is the start of a new section or standalone parameter. */
+ } else { /*Else it is the start of a new section or standalone parameter. */
/* Take first token as the parameter name. */
token = strtok(line, " :\t");
strcpy(tmpStr, token);
@@ -276,7 +277,13 @@ static void parse_value(char *line, struct swift_params *params) {
find_duplicate_params(params, tmpStr);
strcpy(section, tmpSectionName);
- strcpy(params->section[params->sectionCount++].name, tmpSectionName);
+ strcpy(params->section[params->sectionCount].name, tmpSectionName);
+ if (params->sectionCount == PARSER_MAX_NO_OF_SECTIONS - 1) {
+ error(
+ "Maximal number of sections in parameter file reached. Aborting !");
+ } else {
+ params->sectionCount++;
+ }
inSection = 1;
isFirstParam = 1;
} else {
@@ -294,7 +301,14 @@ static void parse_value(char *line, struct swift_params *params) {
/* Must be a standalone parameter so no need to prefix name with a
* section. */
strcpy(params->data[params->paramCount].name, tmpStr);
- strcpy(params->data[params->paramCount++].value, token);
+ strcpy(params->data[params->paramCount].value, token);
+ if (params->paramCount == PARSER_MAX_NO_OF_PARAMS - 1) {
+ error(
+ "Maximal number of parameters in parameter file reached. Aborting "
+ "!");
+ } else {
+ params->paramCount++;
+ }
inSection = 0;
isFirstParam = 1;
}
@@ -346,7 +360,12 @@ static void parse_section_param(char *line, int *isFirstParam,
find_duplicate_params(params, paramName);
strcpy(params->data[params->paramCount].name, paramName);
- strcpy(params->data[params->paramCount++].value, token);
+ strcpy(params->data[params->paramCount].value, token);
+ if (params->paramCount == PARSER_MAX_NO_OF_PARAMS - 1) {
+ error("Maximal number of parameters in parameter file reached. Aborting !");
+ } else {
+ params->paramCount++;
+ }
}
/**
@@ -376,7 +395,8 @@ int parser_get_param_int(const struct swift_params *params, const char *name) {
}
}
- error("Cannot find '%s' in the structure, in file '%s'.", name,params->fileName);
+ error("Cannot find '%s' in the structure, in file '%s'.", name,
+ params->fileName);
return 0;
}
@@ -408,7 +428,8 @@ char parser_get_param_char(const struct swift_params *params,
}
}
- error("Cannot find '%s' in the structure, in file '%s'.", name,params->fileName);
+ error("Cannot find '%s' in the structure, in file '%s'.", name,
+ params->fileName);
return 0;
}
@@ -440,7 +461,8 @@ float parser_get_param_float(const struct swift_params *params,
}
}
- error("Cannot find '%s' in the structure, in file '%s'.", name,params->fileName);
+ error("Cannot find '%s' in the structure, in file '%s'.", name,
+ params->fileName);
return 0.f;
}
@@ -471,7 +493,8 @@ double parser_get_param_double(const struct swift_params *params,
}
}
- error("Cannot find '%s' in the structure, in file '%s'.", name,params->fileName);
+ error("Cannot find '%s' in the structure, in file '%s'.", name,
+ params->fileName);
return 0.;
}
@@ -558,3 +581,11 @@ void parser_write_params_to_file(const struct swift_params *params,
fclose(file);
}
+
+#if defined(HAVE_HDF5)
+void parser_write_params_to_hdf5(const struct swift_params *params, hid_t grp) {
+
+ for (int i = 0; i < params->paramCount; i++)
+ writeAttribute_s(grp, params->data[i].name, params->data[i].value);
+}
+#endif
diff --git a/src/parser.h b/src/parser.h
index 2d50ddc49bf2231235bbb73040b067f77462d33e..21156a37b3bc76d8c3ba01385cb991b2fe213a85 100644
--- a/src/parser.h
+++ b/src/parser.h
@@ -22,9 +22,14 @@
/* Config parameters. */
#include "../config.h"
+#if defined(HAVE_HDF5)
+#include <hdf5.h>
+#endif
+
/* Some constants. */
#define PARSER_MAX_LINE_SIZE 256
-#define PARSER_MAX_NO_OF_PARAMS 512
+#define PARSER_MAX_NO_OF_PARAMS 256
+#define PARSER_MAX_NO_OF_SECTIONS 64
/* A parameter in the input file */
struct parameter {
@@ -38,7 +43,7 @@ struct section {
/* The array of parameters read from a file */
struct swift_params {
- struct section section[PARSER_MAX_NO_OF_PARAMS];
+ struct section section[PARSER_MAX_NO_OF_SECTIONS];
struct parameter data[PARSER_MAX_NO_OF_PARAMS];
int sectionCount;
int paramCount;
@@ -59,4 +64,8 @@ double parser_get_param_double(const struct swift_params *params,
void parser_get_param_string(const struct swift_params *params,
const char *name, char *retParam);
+#if defined(HAVE_HDF5)
+void parser_write_params_to_hdf5(const struct swift_params *params, hid_t grp);
+#endif
+
#endif /* SWIFT_PARSER_H */
diff --git a/src/partition.c b/src/partition.c
index 70249a679edfe72bda97eddfc918a8526659d995..432943f05f78b41086c7d0d9cc29b1456f094d54 100644
--- a/src/partition.c
+++ b/src/partition.c
@@ -31,11 +31,11 @@
#include "../config.h"
/* Standard headers. */
+#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
-#include <float.h>
/* MPI headers. */
#ifdef WITH_MPI
@@ -59,13 +59,12 @@
/* Simple descriptions of initial partition types for reports. */
const char *initial_partition_name[] = {
- "gridded cells", "vectorized point associated cells",
+ "gridded cells", "vectorized point associated cells",
"METIS particle weighted cells", "METIS unweighted cells"};
/* Simple descriptions of repartition types for reports. */
const char *repartition_name[] = {
- "no",
- "METIS edge and vertex time weighted cells",
+ "no", "METIS edge and vertex time weighted cells",
"METIS particle count vertex weighted cells",
"METIS time edge weighted cells",
"METIS particle count vertex and time edge cells"};
@@ -782,8 +781,9 @@ void partition_initial_partition(struct partition *initial_partition,
struct cell *c;
/* If we've got the wrong number of nodes, fail. */
- if (nr_nodes != initial_partition->grid[0] * initial_partition->grid[1] *
- initial_partition->grid[2])
+ if (nr_nodes !=
+ initial_partition->grid[0] * initial_partition->grid[1] *
+ initial_partition->grid[2])
error("Grid size does not match number of nodes.");
/* Run through the cells and set their nodeID. */
@@ -792,8 +792,9 @@ void partition_initial_partition(struct partition *initial_partition,
c = &s->cells[k];
for (j = 0; j < 3; j++)
ind[j] = c->loc[j] / s->dim[j] * initial_partition->grid[j];
- c->nodeID = ind[0] + initial_partition->grid[0] *
- (ind[1] + initial_partition->grid[1] * ind[2]);
+ c->nodeID = ind[0] +
+ initial_partition->grid[0] *
+ (ind[1] + initial_partition->grid[1] * ind[2]);
// message("cell at [%e,%e,%e]: ind = [%i,%i,%i], nodeID = %i", c->loc[0],
// c->loc[1], c->loc[2], ind[0], ind[1], ind[2], c->nodeID);
}
diff --git a/src/potentials.h b/src/potentials.h
index 50a8903095a8a8218e8d02e3f42af07e40be5b43..c3db5a3d2231c93b1d83aa94bca83f7f7d7106bc 100644
--- a/src/potentials.h
+++ b/src/potentials.h
@@ -30,8 +30,8 @@
/* Local includes. */
#include "const.h"
#include "error.h"
-#include "part.h"
#include "parser.h"
+#include "part.h"
#include "physical_constants.h"
#include "units.h"
@@ -58,11 +58,10 @@ struct external_potential {
* @param phys_const The physical constants in internal units.
* @param g Pointer to the g-particle data.
*/
-__attribute__((always_inline))
- INLINE static float external_gravity_pointmass_timestep(
- const struct external_potential* potential,
- const struct phys_const* const phys_const,
- const struct gpart* const g) {
+__attribute__((always_inline)) INLINE static float
+external_gravity_pointmass_timestep(const struct external_potential* potential,
+ const struct phys_const* const phys_const,
+ const struct gpart* const g) {
const float G_newton = phys_const->const_newton_G;
const float dx = g->x[0] - potential->point_mass.x;
diff --git a/src/proxy.c b/src/proxy.c
index 02263a5653bdcdd2d1bf0a86523ed1a599d4bf21..efe3a3eec108d44d5b9bf8b4718dc025464f8762 100644
--- a/src/proxy.c
+++ b/src/proxy.c
@@ -249,8 +249,8 @@ void proxy_parts_exch2(struct proxy *p) {
} while (p->nr_parts_in > p->size_parts_in);
free(p->parts_in);
free(p->xparts_in);
- if ((p->parts_in = (struct part *)malloc(
- sizeof(struct part) *p->size_parts_in)) == NULL ||
+ if ((p->parts_in = (struct part *)malloc(sizeof(struct part) *
+ p->size_parts_in)) == NULL ||
(p->xparts_in = (struct xpart *)malloc(sizeof(struct xpart) *
p->size_parts_in)) == NULL)
error("Failed to re-allocate parts_in buffers.");
@@ -310,7 +310,7 @@ void proxy_parts_load(struct proxy *p, const struct part *parts,
} while (p->nr_parts_out + N > p->size_parts_out);
struct part *tp;
struct xpart *txp;
- if ((tp = (struct part *)malloc(sizeof(struct part) *p->size_parts_out)) ==
+ if ((tp = (struct part *)malloc(sizeof(struct part) * p->size_parts_out)) ==
NULL ||
(txp = (struct xpart *)malloc(sizeof(struct xpart) *
p->size_parts_out)) == NULL)
@@ -395,8 +395,8 @@ void proxy_init(struct proxy *p, int mynodeID, int nodeID) {
/* Allocate the part send and receive buffers, if needed. */
if (p->parts_in == NULL) {
p->size_parts_in = proxy_buffinit;
- if ((p->parts_in = (struct part *)malloc(
- sizeof(struct part) *p->size_parts_in)) == NULL ||
+ if ((p->parts_in = (struct part *)malloc(sizeof(struct part) *
+ p->size_parts_in)) == NULL ||
(p->xparts_in = (struct xpart *)malloc(sizeof(struct xpart) *
p->size_parts_in)) == NULL)
error("Failed to allocate parts_in buffers.");
@@ -404,8 +404,8 @@ void proxy_init(struct proxy *p, int mynodeID, int nodeID) {
p->nr_parts_in = 0;
if (p->parts_out == NULL) {
p->size_parts_out = proxy_buffinit;
- if ((p->parts_out = (struct part *)malloc(
- sizeof(struct part) *p->size_parts_out)) == NULL ||
+ if ((p->parts_out = (struct part *)malloc(sizeof(struct part) *
+ p->size_parts_out)) == NULL ||
(p->xparts_out = (struct xpart *)malloc(sizeof(struct xpart) *
p->size_parts_out)) == NULL)
error("Failed to allocate parts_out buffers.");
diff --git a/src/queue.c b/src/queue.c
index 6b788d7376ba4bdc95f1b1d918ab52a9514e7b4a..dfdac883f9713eb57bb2dc45eb97774983e3a9b2 100644
--- a/src/queue.c
+++ b/src/queue.c
@@ -133,7 +133,7 @@ void queue_init(struct queue *q, struct task *tasks) {
struct task *queue_gettask(struct queue *q, const struct task *prev,
int blocking) {
- lock_type *qlock = &q->lock;
+ swift_lock_type *qlock = &q->lock;
struct task *res = NULL;
/* Grab the task lock. */
diff --git a/src/queue.h b/src/queue.h
index 9ce52ea5404db727f29d0f1cf898f5f5a4f6d935..b28039a30724fa0c086f9111470af2566933e3ac 100644
--- a/src/queue.h
+++ b/src/queue.h
@@ -41,7 +41,7 @@ extern int queue_counter[queue_counter_count];
struct queue {
/* The lock to access this queue. */
- lock_type lock;
+ swift_lock_type lock;
/* Size, count and next element. */
int size, count;
diff --git a/src/riemann.h b/src/riemann.h
index ad37490d7249bafe776b58a609064cbd0bc23abe..d647b021167317d14f4cd7316d09c247794f3d23 100644
--- a/src/riemann.h
+++ b/src/riemann.h
@@ -21,11 +21,11 @@
/* gives us const_hydro_gamma and tells us which floating point type to use */
#include "const.h"
+#include "error.h"
+#include "float.h"
#include "math.h"
#include "stdio.h"
-#include "float.h"
#include "stdlib.h"
-#include "error.h"
#define HLLC_SOLVER
diff --git a/src/riemann/riemann_exact.h b/src/riemann/riemann_exact.h
index b768cde5f4f5dfd0463cc8a582a1af0a17607bbe..a2f3c30fb1daf5d53bf35abe4ca7e73eafba6018 100644
--- a/src/riemann/riemann_exact.h
+++ b/src/riemann/riemann_exact.h
@@ -78,9 +78,9 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_fb(GFLOAT p,
* @param aL The left sound speed
* @param aR The right sound speed
*/
-__attribute__((always_inline))
- INLINE static GFLOAT riemann_f(GFLOAT p, GFLOAT* WL, GFLOAT* WR, GFLOAT vL,
- GFLOAT vR, GFLOAT aL, GFLOAT aR) {
+__attribute__((always_inline)) INLINE static GFLOAT riemann_f(
+ GFLOAT p, GFLOAT* WL, GFLOAT* WR, GFLOAT vL, GFLOAT vR, GFLOAT aL,
+ GFLOAT aR) {
return riemann_fb(p, WL, aL) + riemann_fb(p, WR, aR) + (vR - vL);
}
diff --git a/src/riemann/riemann_hllc.h b/src/riemann/riemann_hllc.h
index 3fcc8b534ce65d4d364c400dc43e1992f39264b9..6c583f6410f53ed64d630082926d816129768fab 100644
--- a/src/riemann/riemann_hllc.h
+++ b/src/riemann/riemann_hllc.h
@@ -63,13 +63,15 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
all these speeds are along the interface normal, since uL and uR are */
qL = 1.;
if (pstar > WL[4]) {
- qL = sqrtf(1. + 0.5 * (const_hydro_gamma + 1.) / const_hydro_gamma *
- (pstar / WL[4] - 1.));
+ qL = sqrtf(1. +
+ 0.5 * (const_hydro_gamma + 1.) / const_hydro_gamma *
+ (pstar / WL[4] - 1.));
}
qR = 1.;
if (pstar > WR[4]) {
- qR = sqrtf(1. + 0.5 * (const_hydro_gamma + 1.) / const_hydro_gamma *
- (pstar / WR[4] - 1.));
+ qR = sqrtf(1. +
+ 0.5 * (const_hydro_gamma + 1.) / const_hydro_gamma *
+ (pstar / WR[4] - 1.));
}
SL = uL - aL * qL;
SR = uR + aR * qR;
diff --git a/src/runner.c b/src/runner.c
index 7a2bcd8563c62250f44183aef3f4b1457623633e..25bd27f109b09f93fafe78fe32350a2dceea0187 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -42,52 +42,67 @@
#include "atomic.h"
#include "const.h"
#include "debug.h"
+#include "drift.h"
#include "engine.h"
#include "error.h"
#include "gravity.h"
-#include "hydro_properties.h"
#include "hydro.h"
+#include "hydro_properties.h"
+#include "kick.h"
#include "minmax.h"
#include "scheduler.h"
#include "space.h"
#include "task.h"
#include "timers.h"
+#include "timestep.h"
/* Orientation of the cell pairs */
const float runner_shift[13 * 3] = {
- 5.773502691896258e-01, 5.773502691896258e-01, 5.773502691896258e-01,
- 7.071067811865475e-01, 7.071067811865475e-01, 0.0,
- 5.773502691896258e-01, 5.773502691896258e-01, -5.773502691896258e-01,
- 7.071067811865475e-01, 0.0, 7.071067811865475e-01,
- 1.0, 0.0, 0.0,
- 7.071067811865475e-01, 0.0, -7.071067811865475e-01,
- 5.773502691896258e-01, -5.773502691896258e-01, 5.773502691896258e-01,
- 7.071067811865475e-01, -7.071067811865475e-01, 0.0,
- 5.773502691896258e-01, -5.773502691896258e-01, -5.773502691896258e-01,
- 0.0, 7.071067811865475e-01, 7.071067811865475e-01,
- 0.0, 1.0, 0.0,
- 0.0, 7.071067811865475e-01, -7.071067811865475e-01,
- 0.0, 0.0, 1.0, };
+ 5.773502691896258e-01,
+ 5.773502691896258e-01,
+ 5.773502691896258e-01,
+ 7.071067811865475e-01,
+ 7.071067811865475e-01,
+ 0.0,
+ 5.773502691896258e-01,
+ 5.773502691896258e-01,
+ -5.773502691896258e-01,
+ 7.071067811865475e-01,
+ 0.0,
+ 7.071067811865475e-01,
+ 1.0,
+ 0.0,
+ 0.0,
+ 7.071067811865475e-01,
+ 0.0,
+ -7.071067811865475e-01,
+ 5.773502691896258e-01,
+ -5.773502691896258e-01,
+ 5.773502691896258e-01,
+ 7.071067811865475e-01,
+ -7.071067811865475e-01,
+ 0.0,
+ 5.773502691896258e-01,
+ -5.773502691896258e-01,
+ -5.773502691896258e-01,
+ 0.0,
+ 7.071067811865475e-01,
+ 7.071067811865475e-01,
+ 0.0,
+ 1.0,
+ 0.0,
+ 0.0,
+ 7.071067811865475e-01,
+ -7.071067811865475e-01,
+ 0.0,
+ 0.0,
+ 1.0,
+};
/* Does the axis need flipping ? */
const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
-/* #define MY_CELL 428428428 */
-/* #define DX 0.1 */
-/* #define NX 1000 */
-/* #define CELL_ID \ */
-/* ((int)(c->loc[0] / DX) * NX *NX + (int)(c->loc[1] / DX) * NX + \ */
-/* (int)(c->loc[2] / DX)) */
-/* #define OUT \ */
-/* if (CELL_ID == MY_CELL) { \ */
-/* message(" cell= %d gcount=%d time=%f \n", CELL_ID, c->gcount,
- * r->e->time); \ */
-/* } */
-/* #define OUT message(" cell %d %d %f \n",CELL_ID, c->count, r->e->time); */
-/* #define OUT if(CELL_ID == MY_CELL) message("\n cell %f %f %f %d %d */
-/* %f\n",c->loc[0],c->loc[1],c->loc[2], CELL_ID, c->count, r->e->time); */
-
/* Import the density loop functions. */
#define FUNCTION density
#include "runner_doiact.h"
@@ -138,52 +153,6 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
if (g->ti_end <= ti_current) {
external_gravity(potential, constants, g);
-
- /* /\* check for energy and angular momentum conservation - begin by */
- /* * synchronizing velocity*\/ */
- /* const float dx = g->x[0] - r->e->potential->point_mass.x; */
- /* const float dy = g->x[1] - r->e->potential->point_mass.y; */
- /* const float dz = g->x[2] - r->e->potential->point_mass.z; */
- /* const float dr = sqrtf((dx * dx) + (dy * dy) + (dz * dz)); */
- /* const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz); */
-
- /* const int current_dti = g->ti_end - g->ti_begin; */
- /* const float dt = 0.5f * current_dti * r->e->timeBase; */
- /* const float vx = g->v_full[0] + dt * g->a_grav[0]; */
- /* const float vy = g->v_full[1] + dt * g->a_grav[1]; */
- /* const float vz = g->v_full[2] + dt * g->a_grav[2]; */
-
- /* /\* E/L *\/ */
- /* float L[3], E; */
- /* E = 0.5 * ((vx * vx) + (vy * vy) + (vz * vz)) - */
- /* r->e->physical_constants->newton_gravity * */
- /* r->e->potential->point_mass.mass * rinv; */
- /* L[0] = dy * vz - dz * vy; */
- /* L[1] = dz * vx - dx * vz; */
- /* L[2] = dx * vy - dy * vx; */
- /* if (abs(g->id) == 1) { */
- /* float v2 = vx * vx + vy * vy + vz * vz; */
- /* float fg = r->e->physical_constants->newton_gravity * */
- /* r->e->potential->point_mass.mass * rinv; */
- /* float fga = sqrtf((g->a_grav[0] * g->a_grav[0]) + */
- /* (g->a_grav[1] * g->a_grav[1]) + */
- /* (g->a_grav[2] * g->a_grav[2])) * */
- /* dr; */
- /* // message("grav_external time= %f\t V_c^2= %f GM/r= %f E= %f L[2]=
- * %f */
- /* // x= %f y= %f vx= %f vy= %f\n", r->e->time, v2, fg, E, L[2],
- * g->x[0], */
- /* // g->x[1], vx, vy); */
- /* message("%f\t %f %f %f %f %f %f %f %f %f %f %f %f %f\n", r->e->time,
- */
- /* g->tx, g->tv, dt, v2, fg, fga, dr, E, L[2], g->x[0], g->x[1],
- */
- /* vx, vy); */
- /* /\* message(" G=%e M=%e\n", r->e->physical_constants->newton_gravity,
- */
- /* * r->e->potential->point_mass.mass); *\/ */
- /* /\* exit(-1); *\/ */
- /* } */
}
}
if (timer) TIMER_TOC(timer_dograv_external);
@@ -795,7 +764,12 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
struct gpart *const gparts = c->gparts;
float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
+ double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
+ double mom[3] = {0.0, 0.0, 0.0};
+ double ang_mom[3] = {0.0, 0.0, 0.0};
+
TIMER_TIC
+
#ifdef TASK_VERBOSE
OUT;
#endif
@@ -811,14 +785,7 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
struct gpart *const gp = &gparts[k];
/* Drift... */
- gp->x[0] += gp->v_full[0] * dt;
- gp->x[1] += gp->v_full[1] * dt;
- gp->x[2] += gp->v_full[2] * dt;
-
- /* Compute offset since last cell construction */
- gp->x_diff[0] -= gp->v_full[0] * dt;
- gp->x_diff[1] -= gp->v_full[1] * dt;
- gp->x_diff[2] -= gp->v_full[2] * dt;
+ drift_gpart(gp, dt, timeBase, ti_old, ti_current);
/* Compute (square of) motion since last cell construction */
const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
@@ -835,40 +802,8 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
struct part *const p = &parts[k];
struct xpart *const xp = &xparts[k];
- /* Useful quantity */
- const float h_inv = 1.0f / p->h;
-
/* Drift... */
- p->x[0] += xp->v_full[0] * dt;
- p->x[1] += xp->v_full[1] * dt;
- p->x[2] += xp->v_full[2] * dt;
-
- /* Predict velocities (for hydro terms) */
- p->v[0] += p->a_hydro[0] * dt;
- p->v[1] += p->a_hydro[1] * dt;
- p->v[2] += p->a_hydro[2] * dt;
-
- /* Predict smoothing length */
- const float w1 = p->h_dt * h_inv * dt;
- if (fabsf(w1) < 0.2f)
- p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
- else
- p->h *= expf(w1);
-
- /* Predict density */
- const float w2 = -3.0f * p->h_dt * h_inv * dt;
- if (fabsf(w2) < 0.2f)
- p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
- else
- p->rho *= expf(w2);
-
- /* Predict the values of the extra fields */
- hydro_predict_extra(p, xp, ti_old, ti_current, timeBase);
-
- /* Compute offset since last cell construction */
- xp->x_diff[0] -= xp->v_full[0] * dt;
- xp->x_diff[1] -= xp->v_full[1] * dt;
- xp->x_diff[2] -= xp->v_full[2] * dt;
+ drift_part(p, xp, dt, timeBase, ti_old, ti_current);
/* Compute (square of) motion since last cell construction */
const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
@@ -878,6 +813,34 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
/* Maximal smoothing length */
h_max = fmaxf(p->h, h_max);
+
+ /* Now collect quantities for statistics */
+
+ const float half_dt =
+ (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
+ const double x[3] = {p->x[0], p->x[1], p->x[2]};
+ const float v[3] = {xp->v_full[0] + p->a_hydro[0] * half_dt,
+ xp->v_full[1] + p->a_hydro[1] * half_dt,
+ xp->v_full[2] + p->a_hydro[2] * half_dt};
+ const float m = p->mass;
+
+ /* Collect mass */
+ mass += m;
+
+ /* Collect momentum */
+ mom[0] += m * v[0];
+ mom[1] += m * v[1];
+ mom[2] += m * v[2];
+
+ /* Collect angular momentum */
+ ang_mom[0] += m * (x[1] * v[2] - x[2] * v[1]);
+ ang_mom[1] += m * (x[2] * v[0] - x[0] * v[2]);
+ ang_mom[2] += m * (x[0] * v[1] - x[1] * v[0]);
+
+ /* Collect energies. */
+ e_kin += 0.5 * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
+ e_pot += 0.;
+ e_int += m * hydro_get_internal_energy(p, half_dt);
}
/* Now, get the maximal particle motion from its square */
@@ -890,52 +853,64 @@ void runner_do_drift(struct runner *r, struct cell *c, int timer) {
/* Loop over the progeny. */
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
+
+ /* Recurse */
struct cell *cp = c->progeny[k];
runner_do_drift(r, cp, 0);
+ /* Collect */
dx_max = fmaxf(dx_max, cp->dx_max);
h_max = fmaxf(h_max, cp->h_max);
+ mass += cp->mass;
+ e_kin += cp->e_kin;
+ e_int += cp->e_int;
+ e_pot += cp->e_pot;
+ mom[0] += cp->mom[0];
+ mom[1] += cp->mom[1];
+ mom[2] += cp->mom[2];
+ ang_mom[0] += cp->ang_mom[0];
+ ang_mom[1] += cp->ang_mom[1];
+ ang_mom[2] += cp->ang_mom[2];
}
}
/* Store the values */
c->h_max = h_max;
c->dx_max = dx_max;
+ c->mass = mass;
+ c->e_kin = e_kin;
+ c->e_int = e_int;
+ c->e_pot = e_pot;
+ c->mom[0] = mom[0];
+ c->mom[1] = mom[1];
+ c->mom[2] = mom[2];
+ c->ang_mom[0] = ang_mom[0];
+ c->ang_mom[1] = ang_mom[1];
+ c->ang_mom[2] = ang_mom[2];
if (timer) TIMER_TOC(timer_drift);
}
/**
- * @brief Kick particles in momentum space and collect statistics
+ * @brief Kick particles in momentum space and collect statistics (fixed
+ * time-step case)
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
-void runner_do_kick(struct runner *r, struct cell *c, int timer) {
+void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer) {
- const float global_dt_min = r->e->dt_min;
- const float global_dt_max = r->e->dt_max;
- const int ti_current = r->e->ti_current;
+ const double global_dt = r->e->dt_max;
const double timeBase = r->e->timeBase;
- const double timeBase_inv = 1.0 / r->e->timeBase;
const int count = c->count;
const int gcount = c->gcount;
struct part *const parts = c->parts;
struct xpart *const xparts = c->xparts;
struct gpart *const gparts = c->gparts;
- const struct external_potential *potential = r->e->external_potential;
- const struct hydro_props *hydro_properties = r->e->hydro_properties;
- const struct phys_const *constants = r->e->physical_constants;
- const float ln_max_h_change = hydro_properties->log_max_h_change;
- const int is_fixdt =
- (r->e->policy & engine_policy_fixdt) == engine_policy_fixdt;
int updated = 0, g_updated = 0;
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
- double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
- float mom[3] = {0.0f, 0.0f, 0.0f};
- float ang[3] = {0.0f, 0.0f, 0.0f};
TIMER_TIC
@@ -943,86 +918,143 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
OUT;
#endif
+ /* The new time-step */
+ const int new_dti = global_dt / timeBase;
+
/* No children? */
if (!c->split) {
- /* Loop over the g-particles and kick the active ones. */
+ /* Loop over the g-particles and kick everyone. */
for (int k = 0; k < gcount; k++) {
/* Get a handle on the part. */
struct gpart *const gp = &gparts[k];
- /* If the g-particle has no counterpart and needs to be kicked */
+ /* If the g-particle has no counterpart */
if (gp->id < 0) {
- if (is_fixdt || gp->ti_end <= ti_current) {
+ /* First, finish the force calculation */
+ gravity_end_force(gp);
- /* First, finish the force calculation */
- gravity_end_force(gp);
+ /* Kick the g-particle forward */
+ kick_gpart(gp, new_dti, timeBase);
- /* Now we are ready to compute the next time-step size */
- int new_dti;
+ /* Number of updated g-particles */
+ g_updated++;
- if (is_fixdt) {
+ /* Minimal time for next end of time-step */
+ ti_end_min = min(gp->ti_end, ti_end_min);
+ ti_end_max = max(gp->ti_end, ti_end_max);
+ }
+ }
- /* Now we have a time step, proceed with the kick */
- new_dti = global_dt_max * timeBase_inv;
+ /* Now do the hydro ones... */
- } else {
+ /* Loop over the particles and kick everyone. */
+ for (int k = 0; k < count; k++) {
- /* Compute the next timestep (gravity condition) */
- const float new_dt_external =
- gravity_compute_timestep_external(potential, constants, gp);
- const float new_dt_self =
- gravity_compute_timestep_self(constants, gp);
+ /* Get a handle on the part. */
+ struct part *const p = &parts[k];
+ struct xpart *const xp = &xparts[k];
- float new_dt = fminf(new_dt_external, new_dt_self);
+ /* First, finish the force loop */
+ p->h_dt *= p->h * 0.333333333f;
- /* Limit timestep within the allowed range */
- new_dt = fminf(new_dt, global_dt_max);
- new_dt = fmaxf(new_dt, global_dt_min);
+ /* And do the same of the extra variable */
+ hydro_end_force(p);
+ if (p->gpart != NULL) gravity_end_force(p->gpart);
- /* Convert to integer time */
- new_dti = new_dt * timeBase_inv;
+ /* Kick the particle forward */
+ kick_part(p, xp, new_dti, timeBase);
- /* Recover the current timestep */
- const int current_dti = gp->ti_end - gp->ti_begin;
+ /* Number of updated particles */
+ updated++;
+ if (p->gpart != NULL) g_updated++;
- /* Limit timestep increase */
- if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
+ /* Minimal time for next end of time-step */
+ ti_end_min = min(p->ti_end, ti_end_min);
+ ti_end_max = max(p->ti_end, ti_end_max);
+ }
+ }
- /* Put this timestep on the time line */
- int dti_timeline = max_nr_timesteps;
- while (new_dti < dti_timeline) dti_timeline /= 2;
+ /* Otherwise, aggregate data from children. */
+ else {
- new_dti = dti_timeline;
+ /* Loop over the progeny. */
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) {
+ struct cell *const cp = c->progeny[k];
- /* Make sure we are allowed to increase the timestep size */
- if (new_dti > current_dti) {
- if ((max_nr_timesteps - gp->ti_end) % new_dti > 0)
- new_dti = current_dti;
- }
+ /* Recurse */
+ runner_do_kick_fixdt(r, cp, 0);
- /* Now we have a time step, proceed with the kick */
- }
+ /* And aggregate */
+ updated += cp->updated;
+ g_updated += cp->g_updated;
+ ti_end_min = min(cp->ti_end_min, ti_end_min);
+ ti_end_max = max(cp->ti_end_max, ti_end_max);
+ }
+ }
- /* Compute the time step for this kick */
- const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
- const int ti_end = gp->ti_end + new_dti / 2;
- const double dt = (ti_end - ti_start) * timeBase;
- const double half_dt = (ti_end - gp->ti_end) * timeBase;
+ /* Store the values. */
+ c->updated = updated;
+ c->g_updated = g_updated;
+ c->ti_end_min = ti_end_min;
+ c->ti_end_max = ti_end_max;
- /* Move particle forward in time */
- gp->ti_begin = gp->ti_end;
- gp->ti_end = gp->ti_begin + new_dti;
+ if (timer) TIMER_TOC(timer_kick);
+}
- /* Kick particles in momentum space */
- gp->v_full[0] += gp->a_grav[0] * dt;
- gp->v_full[1] += gp->a_grav[1] * dt;
- gp->v_full[2] += gp->a_grav[2] * dt;
+/**
+ * @brief Kick particles in momentum space and collect statistics (floating
+ * time-step case)
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_kick(struct runner *r, struct cell *c, int timer) {
- /* Extra kick work */
- gravity_kick_extra(gp, dt, half_dt);
+ const struct engine *e = r->e;
+ const double timeBase = e->timeBase;
+ const int ti_current = r->e->ti_current;
+ const int count = c->count;
+ const int gcount = c->gcount;
+ struct part *const parts = c->parts;
+ struct xpart *const xparts = c->xparts;
+ struct gpart *const gparts = c->gparts;
+
+ int updated = 0, g_updated = 0;
+ int ti_end_min = max_nr_timesteps, ti_end_max = 0;
+
+ TIMER_TIC
+
+#ifdef TASK_VERBOSE
+ OUT;
+#endif
+
+ /* No children? */
+ if (!c->split) {
+
+ /* Loop over the g-particles and kick the active ones. */
+ for (int k = 0; k < gcount; k++) {
+
+ /* Get a handle on the part. */
+ struct gpart *const gp = &gparts[k];
+
+ /* If the g-particle has no counterpart and needs to be kicked */
+ if (gp->id < 0) {
+
+ if (gp->ti_end <= ti_current) {
+
+ /* First, finish the force calculation */
+ gravity_end_force(gp);
+
+ /* Compute the next timestep */
+ const int new_dti = get_gpart_timestep(gp, e);
+
+ /* Now we have a time step, proceed with the kick */
+ kick_gpart(gp, new_dti, timeBase);
/* Number of updated g-particles */
g_updated++;
@@ -1044,7 +1076,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
struct xpart *const xp = &xparts[k];
/* If particle needs to be kicked */
- if (is_fixdt || p->ti_end <= ti_current) {
+ if (p->ti_end <= ti_current) {
/* First, finish the force loop */
p->h_dt *= p->h * 0.333333333f;
@@ -1053,142 +1085,17 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
hydro_end_force(p);
if (p->gpart != NULL) gravity_end_force(p->gpart);
- /* Now we are ready to compute the next time-step size */
- int new_dti;
-
- if (is_fixdt) {
-
- /* Now we have a time step, proceed with the kick */
- new_dti = global_dt_max * timeBase_inv;
-
- } else {
-
- /* Compute the next timestep (hydro condition) */
- const float new_dt_hydro =
- hydro_compute_timestep(p, xp, hydro_properties);
-
- /* Compute the next timestep (gravity condition) */
- float new_dt_grav = FLT_MAX;
- if (p->gpart != NULL) {
-
- const float new_dt_external = gravity_compute_timestep_external(
- potential, constants, p->gpart);
- const float new_dt_self =
- gravity_compute_timestep_self(constants, p->gpart);
-
- new_dt_grav = fminf(new_dt_external, new_dt_self);
- }
-
- /* Final time-step is minimum of hydro and gravity */
- float new_dt = fminf(new_dt_hydro, new_dt_grav);
-
- /* Limit change in h */
- const float dt_h_change =
- (p->h_dt != 0.0f) ? fabsf(ln_max_h_change * p->h / p->h_dt)
- : FLT_MAX;
-
- new_dt = fminf(new_dt, dt_h_change);
-
- /* Limit timestep within the allowed range */
- new_dt = fminf(new_dt, global_dt_max);
- new_dt = fmaxf(new_dt, global_dt_min);
+ /* Compute the next timestep (hydro condition) */
+ const int new_dti = get_part_timestep(p, xp, e);
- /* Convert to integer time */
- new_dti = new_dt * timeBase_inv;
-
- /* Recover the current timestep */
- const int current_dti = p->ti_end - p->ti_begin;
-
- /* Limit timestep increase */
- if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
-
- /* Put this timestep on the time line */
- int dti_timeline = max_nr_timesteps;
- while (new_dti < dti_timeline) dti_timeline /= 2;
-
- new_dti = dti_timeline;
-
- /* Make sure we are allowed to increase the timestep size */
- if (new_dti > current_dti) {
- if ((max_nr_timesteps - p->ti_end) % new_dti > 0)
- new_dti = current_dti;
- }
-
- /* Now we have a time step, proceed with the kick */
- }
-
- /* Compute the time step for this kick */
- const int ti_start = (p->ti_begin + p->ti_end) / 2;
- const int ti_end = p->ti_end + new_dti / 2;
- const double dt = (ti_end - ti_start) * timeBase;
- const double half_dt = (ti_end - p->ti_end) * timeBase;
-
- /* Move particle forward in time */
- p->ti_begin = p->ti_end;
- p->ti_end = p->ti_begin + new_dti;
- if (p->gpart != NULL) {
- p->gpart->ti_begin = p->ti_begin;
- p->gpart->ti_end = p->ti_end;
- }
-
- /* Get the acceleration */
- float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]};
- if (p->gpart != NULL) {
- a_tot[0] += p->gpart->a_grav[0];
- a_tot[1] += p->gpart->a_grav[1];
- a_tot[1] += p->gpart->a_grav[2];
- }
-
- /* Kick particles in momentum space */
- xp->v_full[0] += a_tot[0] * dt;
- xp->v_full[1] += a_tot[1] * dt;
- xp->v_full[2] += a_tot[2] * dt;
-
- if (p->gpart != NULL) {
- p->gpart->v_full[0] = xp->v_full[0];
- p->gpart->v_full[1] = xp->v_full[1];
- p->gpart->v_full[2] = xp->v_full[2];
- }
-
- /* Go back by half-step for the hydro velocity */
- p->v[0] = xp->v_full[0] - half_dt * a_tot[0];
- p->v[1] = xp->v_full[1] - half_dt * a_tot[1];
- p->v[2] = xp->v_full[2] - half_dt * a_tot[2];
-
- /* Extra kick work */
- hydro_kick_extra(p, xp, dt, half_dt);
- if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt, half_dt);
+ /* Now we have a time step, proceed with the kick */
+ kick_part(p, xp, new_dti, timeBase);
/* Number of updated particles */
updated++;
if (p->gpart != NULL) g_updated++;
}
- /* Now collect quantities for statistics */
-
- const double x[3] = {p->x[0], p->x[1], p->x[2]};
- const float v_full[3] = {xp->v_full[0], xp->v_full[1], xp->v_full[2]};
- const float m = p->mass;
-
- /* Collect mass */
- mass += m;
-
- /* Collect momentum */
- mom[0] += m * v_full[0];
- mom[1] += m * v_full[1];
- mom[2] += m * v_full[2];
-
- /* Collect angular momentum */
- ang[0] += m * (x[1] * v_full[2] - x[2] * v_full[1]);
- ang[1] += m * (x[2] * v_full[0] - x[0] * v_full[2]);
- ang[2] += m * (x[0] * v_full[1] - x[1] * v_full[0]);
-
- /* Collect total energy. */
- e_kin += 0.5 * m * (v_full[0] * v_full[0] + v_full[1] * v_full[1] +
- v_full[2] * v_full[2]);
- e_pot += 0.f; /* No gravitational potential thus far */
- e_int += hydro_get_internal_energy(p);
-
/* Minimal time for next end of time-step */
ti_end_min = min(p->ti_end, ti_end_min);
ti_end_max = max(p->ti_end, ti_end_max);
@@ -1209,16 +1116,6 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
/* And aggregate */
updated += cp->updated;
g_updated += cp->g_updated;
- e_kin += cp->e_kin;
- e_int += cp->e_int;
- e_pot += cp->e_pot;
- mass += cp->mass;
- mom[0] += cp->mom[0];
- mom[1] += cp->mom[1];
- mom[2] += cp->mom[2];
- ang[0] += cp->ang[0];
- ang[1] += cp->ang[1];
- ang[2] += cp->ang[2];
ti_end_min = min(cp->ti_end_min, ti_end_min);
ti_end_max = max(cp->ti_end_max, ti_end_max);
}
@@ -1227,16 +1124,6 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
/* Store the values. */
c->updated = updated;
c->g_updated = g_updated;
- c->e_kin = e_kin;
- c->e_int = e_int;
- c->e_pot = e_pot;
- c->mass = mass;
- c->mom[0] = mom[0];
- c->mom[1] = mom[1];
- c->mom[2] = mom[2];
- c->ang[0] = ang[0];
- c->ang[1] = ang[1];
- c->ang[2] = ang[2];
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
@@ -1372,6 +1259,9 @@ void *runner_main(void *data) {
case task_type_kick:
runner_do_kick(r, ci, 1);
break;
+ case task_type_kick_fixdt:
+ runner_do_kick_fixdt(r, ci, 1);
+ break;
case task_type_send:
break;
case task_type_recv:
diff --git a/src/runner.h b/src/runner.h
index f53822cdbf9608a473fc72e6a2d049cfd6d3c5a2..35e5f56a7b94145eec656b02a1c5568ec39352b6 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -51,6 +51,7 @@ void runner_do_ghost(struct runner *r, struct cell *c);
void runner_do_sort(struct runner *r, struct cell *c, int flag, int clock);
void runner_do_gsort(struct runner *r, struct cell *c, int flag, int clock);
void runner_do_kick(struct runner *r, struct cell *c, int timer);
+void runner_do_kick_fixdt(struct runner *r, struct cell *c, int timer);
void runner_do_drift(struct runner *r, struct cell *c, int timer);
void runner_do_init(struct runner *r, struct cell *c, int timer);
void *runner_main(void *data);
diff --git a/src/scheduler.c b/src/scheduler.c
index 1b545aac5a3e4a153f8bbec5824797f7bb0edf60..196031417be11044e46ecb88d7a8e6c6adb0ab66 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -842,7 +842,8 @@ void scheduler_reset(struct scheduler *s, int size) {
if (s->tasks_ind != NULL) free(s->tasks_ind);
/* Allocate the new lists. */
- if ((s->tasks = (struct task *)malloc(sizeof(struct task) *size)) == NULL ||
+ if ((s->tasks = (struct task *)malloc(sizeof(struct task) * size)) ==
+ NULL ||
(s->tasks_ind = (int *)malloc(sizeof(int) * size)) == NULL)
error("Failed to allocate task lists.");
}
@@ -1328,7 +1329,7 @@ void scheduler_init(struct scheduler *s, struct space *space, int nr_tasks,
/* Init the unlocks. */
if ((s->unlocks = (struct task **)malloc(
- sizeof(struct task *) *scheduler_init_nr_unlocks)) == NULL ||
+ sizeof(struct task *) * scheduler_init_nr_unlocks)) == NULL ||
(s->unlock_ind =
(int *)malloc(sizeof(int) * scheduler_init_nr_unlocks)) == NULL)
error("Failed to allocate unlocks.");
diff --git a/src/scheduler.h b/src/scheduler.h
index 754bfb8c224a8b58b6b16e8b5dc71bd6563136f3..7528e13c2f7901ba9979a3ef7567bef1e23d8ac4 100644
--- a/src/scheduler.h
+++ b/src/scheduler.h
@@ -90,7 +90,7 @@ struct scheduler {
volatile int nr_unlocks, size_unlocks, completed_unlock_writes;
/* Lock for this scheduler. */
- lock_type lock;
+ swift_lock_type lock;
/* Waiting queue. */
pthread_mutex_t sleep_mutex;
diff --git a/src/serial_io.c b/src/serial_io.c
index f2c0dcce1cc258e7c9a43027fbdd45e299b47ae6..5abd3ebc28672d68c4135efe5753dc4713c2d3c6 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -608,7 +608,7 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts,
void write_output_serial(struct engine* e, const char* baseName,
struct UnitSystem* us, int mpi_rank, int mpi_size,
MPI_Comm comm, MPI_Info info) {
- hid_t h_file = 0, h_grp = 0, h_grpsph = 0;
+ hid_t h_file = 0, h_grp = 0;
const size_t Ngas = e->s->nr_parts;
const size_t Ntot = e->s->nr_gparts;
int periodic = e->s->periodic;
@@ -710,10 +710,17 @@ void write_output_serial(struct engine* e, const char* baseName,
writeCodeDescription(h_file);
/* Print the SPH parameters */
- h_grpsph = H5Gcreate(h_file, "/SPH", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
- if (h_grpsph < 0) error("Error while creating SPH group");
- writeSPHflavour(h_grpsph);
- H5Gclose(h_grpsph);
+ h_grp = H5Gcreate(h_file, "/SPH", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating SPH group");
+ writeSPHflavour(h_grp);
+ H5Gclose(h_grp);
+
+ /* Print the runtime parameters */
+ h_grp =
+ H5Gcreate(h_file, "/Parameters", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating parameters group");
+ parser_write_params_to_hdf5(e->parameter_file, h_grp);
+ H5Gclose(h_grp);
/* Print the system of Units */
writeUnitSystem(h_file, us);
diff --git a/src/single_io.c b/src/single_io.c
index d545fb086fa4bc63fe58dee2bfe85d9586997850..fb3bf4368feed9892b098228d85c99f0bc3e724b 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -35,8 +35,8 @@
#include "single_io.h"
/* Local includes. */
-#include "const.h"
#include "common_io.h"
+#include "const.h"
#include "error.h"
/*-----------------------------------------------------------------------------
@@ -562,6 +562,13 @@ void write_output_single(struct engine* e, const char* baseName,
writeSPHflavour(h_grp);
H5Gclose(h_grp);
+ /* Print the runtime parameters */
+ h_grp =
+ H5Gcreate(h_file, "/Parameters", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+ if (h_grp < 0) error("Error while creating parameters group");
+ parser_write_params_to_hdf5(e->parameter_file, h_grp);
+ H5Gclose(h_grp);
+
/* Print the system of Units */
writeUnitSystem(h_file, us);
diff --git a/src/space.c b/src/space.c
index a36906f6818c6bab0dbef717626fe537289f1715..720961673453433f124f0a4088704224d850d39e 100644
--- a/src/space.c
+++ b/src/space.c
@@ -28,8 +28,8 @@
#include <float.h>
#include <limits.h>
#include <math.h>
-#include <string.h>
#include <stdlib.h>
+#include <string.h>
/* MPI headers. */
#ifdef WITH_MPI
diff --git a/src/space.h b/src/space.h
index e70e77f15f0b1ea898fb328746b484a6e61f6738..d448cd1279611c26139672608798089f7c6a7afb 100644
--- a/src/space.h
+++ b/src/space.h
@@ -97,7 +97,7 @@ struct space {
int periodic;
/* General-purpose lock for this space. */
- lock_type lock;
+ swift_lock_type lock;
/* Number of queues in the system. */
int nr_queues;
diff --git a/src/task.c b/src/task.c
index c8f69681c8ac3a142a7fdebd3778eedd5cae546a..4b0b649c63f1af9ceb30a0ff311d55f679a0a503 100644
--- a/src/task.c
+++ b/src/task.c
@@ -51,8 +51,8 @@ const char *taskID_names[task_type_count] = {
"ghost", "drift", "kick", "send", "recv", "grav_pp",
"grav_mm", "grav_up", "grav_down", "grav_external", "comm_root"};
-const char *subtaskID_names[task_type_count] = {"none", "density",
- "force", "grav"};
+const char *subtaskID_names[task_type_count] = {"none", "density", "force",
+ "grav"};
/**
* @brief Computes the overlap between the parts array of two given cells.
diff --git a/src/task.h b/src/task.h
index 8193ce970fd67779aa6c8d62e0be5e0b5ea7cd83..a949e02fce1e9e0fa6b86f27e1c192a45b9073f3 100644
--- a/src/task.h
+++ b/src/task.h
@@ -42,6 +42,7 @@ enum task_types {
task_type_ghost,
task_type_drift,
task_type_kick,
+ task_type_kick_fixdt,
task_type_send,
task_type_recv,
task_type_grav_pp,
@@ -74,7 +75,7 @@ struct task {
char skip, tight, implicit;
int flags, wait, rank, weight;
- lock_type lock;
+ swift_lock_type lock;
struct cell *ci, *cj;
diff --git a/src/timestep.h b/src/timestep.h
new file mode 100644
index 0000000000000000000000000000000000000000..1674958812d1f9e72eb618a9cb026befd610d06d
--- /dev/null
+++ b/src/timestep.h
@@ -0,0 +1,136 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_TIMESTEP_H
+#define SWIFT_TIMESTEP_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+#include "debug.h"
+
+/**
+ * @brief Compute a valid integer time-step form a given time-step
+ *
+ * @param new_dt The time-step to convert.
+ * @param ti_begin The (integer) start of the previous time-step.
+ * @param ti_end The (integer) end of the previous time-step.
+ * @param timeBase_inv The inverse of the system's minimal time-step.
+ */
+__attribute__((always_inline)) INLINE static int get_integer_timestep(
+ float new_dt, int ti_begin, int ti_end, double timeBase_inv) {
+
+ /* Convert to integer time */
+ int new_dti = (int)(new_dt * timeBase_inv);
+
+ /* Recover the current timestep */
+ const int current_dti = ti_end - ti_begin;
+
+ /* Limit timestep increase */
+ if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
+
+ /* Put this timestep on the time line */
+ int dti_timeline = max_nr_timesteps;
+ while (new_dti < dti_timeline) dti_timeline /= 2;
+ new_dti = dti_timeline;
+
+ /* Make sure we are allowed to increase the timestep size */
+ if (new_dti > current_dti) {
+ if ((max_nr_timesteps - ti_end) % new_dti > 0) new_dti = current_dti;
+ }
+
+ return new_dti;
+}
+
+/**
+ * @brief Compute the new (integer) time-step of a given #gpart
+ *
+ * @param gp The #gpart.
+ * @param e The #engine (used to get some constants).
+ */
+__attribute__((always_inline)) INLINE static int get_gpart_timestep(
+ const struct gpart *gp, const struct engine *e) {
+
+ const float new_dt_external = gravity_compute_timestep_external(
+ e->external_potential, e->physical_constants, gp);
+ const float new_dt_self =
+ gravity_compute_timestep_self(e->physical_constants, gp);
+
+ float new_dt = fminf(new_dt_external, new_dt_self);
+
+ /* Limit timestep within the allowed range */
+ new_dt = fminf(new_dt, e->dt_max);
+ new_dt = fmaxf(new_dt, e->dt_min);
+
+ /* Convert to integer time */
+ const int new_dti =
+ get_integer_timestep(new_dt, gp->ti_begin, gp->ti_end, e->timeBase_inv);
+
+ return new_dti;
+}
+
+/**
+ * @brief Compute the new (integer) time-step of a given #part
+ *
+ * @param p The #part.
+ * @param xp The #xpart partner of p.
+ * @param e The #engine (used to get some constants).
+ */
+__attribute__((always_inline)) INLINE static int get_part_timestep(
+ const struct part *p, const struct xpart *xp, const struct engine *e) {
+
+ /* Compute the next timestep (hydro condition) */
+ const float new_dt_hydro = hydro_compute_timestep(p, xp, e->hydro_properties);
+
+ /* Compute the next timestep (gravity condition) */
+ float new_dt_grav = FLT_MAX;
+ if (p->gpart != NULL) {
+
+ const float new_dt_external = gravity_compute_timestep_external(
+ e->external_potential, e->physical_constants, p->gpart);
+ const float new_dt_self =
+ gravity_compute_timestep_self(e->physical_constants, p->gpart);
+
+ new_dt_grav = fminf(new_dt_external, new_dt_self);
+ }
+
+ /* Final time-step is minimum of hydro and gravity */
+ float new_dt = fminf(new_dt_hydro, new_dt_grav);
+
+ /* Limit change in h */
+ const float dt_h_change =
+ (p->h_dt != 0.0f)
+ ? fabsf(e->hydro_properties->log_max_h_change * p->h / p->h_dt)
+ : FLT_MAX;
+
+ new_dt = fminf(new_dt, dt_h_change);
+
+ /* Limit timestep within the allowed range */
+ new_dt = fminf(new_dt, e->dt_max);
+ new_dt = fmaxf(new_dt, e->dt_min);
+
+ /* Convert to integer time */
+ const int new_dti =
+ get_integer_timestep(new_dt, p->ti_begin, p->ti_end, e->timeBase_inv);
+
+ return new_dti;
+}
+
+#endif /* SWIFT_TIMESTEP_H */
diff --git a/src/tools.c b/src/tools.c
index 7b468426b894044dc8e215f5150d9536fa2d4cdd..0363100331ad5b298279e9ebadcf87eab5f9e896 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -20,15 +20,15 @@
******************************************************************************/
#include <math.h>
-#include <stdlib.h>
#include <stddef.h>
#include <stdio.h>
+#include <stdlib.h>
+#include "cell.h"
#include "error.h"
#include "part.h"
-#include "cell.h"
-#include "tools.h"
#include "swift.h"
+#include "tools.h"
/**
* Factorize a given integer, attempts to keep larger pair of factors.
diff --git a/src/tools.h b/src/tools.h
index 8e0212652922fb4afd1ac89a64d4a01c71ecd4a9..5f9f41d033ab03983bde3bb37e87f8a39d2deecd 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -22,8 +22,8 @@
#ifndef SWIFT_TOOL_H
#define SWIFT_TOOL_H
-#include "runner.h"
#include "cell.h"
+#include "runner.h"
void factor(int value, int *f1, int *f2);
void density_dump(int N);
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 504c07b4db391a2bc6fa60b7ec354367008aeec1..3e2f11c7aabac0a7ff2df19f2c48f9f81ea55df5 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -18,9 +18,9 @@
******************************************************************************/
#include <fenv.h>
+#include <stdio.h>
#include <stdlib.h>
#include <string.h>
-#include <stdio.h>
#include <unistd.h>
#include "swift.h"
@@ -301,7 +301,8 @@ int main(int argc, char *argv[]) {
/* Help users... */
message("Smoothing length: h = %f", h * size);
message("Kernel: %s", kernel_name);
- message("Neighbour target: N = %f", h * h * h * 4.0 * M_PI * kernel_gamma3 / 3.0);
+ message("Neighbour target: N = %f",
+ h * h * h * 4.0 * M_PI * kernel_gamma3 / 3.0);
message("Density target: rho = %f", rho);
message("div_v target: div = %f", vel == 2 ? 3.f : 0.f);
message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f);
diff --git a/tests/testKernel.c b/tests/testKernel.c
index 27d45375b91e4143036c32b71a217e80ad630222..182bae5334e1a5061e584212a31186dc4e7f0818 100644
--- a/tests/testKernel.c
+++ b/tests/testKernel.c
@@ -18,9 +18,9 @@
*
******************************************************************************/
-#define NO__SSE2__
-#include "vector.h"
+#define NO__AVX__
#include "kernel_hydro.h"
+#include "vector.h"
#include <stdlib.h>
#include <strings.h>
diff --git a/tests/testPair.c b/tests/testPair.c
index a70cb381fd9a6dec061a100eeba17bfd54b0e973..f9539fc1a444828c65b39e56618eb7bb98bd67de 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -18,9 +18,9 @@
******************************************************************************/
#include <fenv.h>
+#include <stdio.h>
#include <stdlib.h>
#include <string.h>
-#include <stdio.h>
#include <unistd.h>
#include "swift.h"
diff --git a/tests/testParser.c b/tests/testParser.c
index 0b08d20c9e2d48de1858877cf186eaa9d0ac84c0..31979a0a572415af04c11063bad9202a0c313426 100644
--- a/tests/testParser.c
+++ b/tests/testParser.c
@@ -17,11 +17,11 @@
*
******************************************************************************/
-#include "parser.h"
#include <assert.h>
-#include <string.h>
-#include <stdio.h>
#include <math.h>
+#include <stdio.h>
+#include <string.h>
+#include "parser.h"
int main(int argc, char *argv[]) {
const char *input_file = argv[1];
diff --git a/tests/testReading.c b/tests/testReading.c
index 33aeb5095ba499bc0fd18ba15b513e351692432e..2fa88855a70a12265f180cd97528dda855322d1d 100644
--- a/tests/testReading.c
+++ b/tests/testReading.c
@@ -17,9 +17,12 @@
*
******************************************************************************/
-#include "swift.h"
+/* Some standard headers. */
#include <stdlib.h>
+/* Includes. */
+#include "swift.h"
+
int main() {
size_t Ngas = 0, Ngpart = 0;
diff --git a/tests/testSingle.c b/tests/testSingle.c
index eb49a570b93b14734c9e6af37d3d8a2b90d04078..d37f20908a28fb2e1098043dd7fbbcf76bd8c247 100644
--- a/tests/testSingle.c
+++ b/tests/testSingle.c
@@ -22,15 +22,15 @@
#include "../config.h"
/* Some standard headers. */
+#include <fenv.h>
+#include <float.h>
+#include <limits.h>
+#include <math.h>
+#include <pthread.h>
#include <stdio.h>
#include <stdlib.h>
-#include <unistd.h>
#include <string.h>
-#include <pthread.h>
-#include <math.h>
-#include <float.h>
-#include <limits.h>
-#include <fenv.h>
+#include <unistd.h>
/* Conditional headers. */
#ifdef HAVE_LIBZ