diff --git a/README.md b/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..8ef8a65e890507c0467a965405e488aeff8ecf36
--- /dev/null
+++ b/README.md
@@ -0,0 +1,93 @@
+SWIFT: SPH WIth Fine-grained inter-dependent Tasking
+====================================================
+
+[](https://gitlab.cosma.dur.ac.uk/jenkins/job/GNU%20SWIFT%20build/)
+
+SWIFT is a gravity and SPH solver designed to run cosmological simulations
+on peta-scale machines, scaling well up to 10's of thousands of compute
+node.
+
+More general information about SWIFT is available on the project
+[webpages](http://www.swiftsim.com).
+
+For information on how to _run_ SWIFT, please consult the onboarding guide
+available [here](http://www.swiftsim.com/onboarding.pdf). This includes
+dependencies, and a few examples to get you going.
+
+We suggest that you use the latest release branch of SWIFT, rather than the
+current master branch as this will change rapidly. We do, however, like to
+ensure that the master branch will build and run.
+
+This GitHub repository is designed to be an issue tracker, and a space for
+the public to submit patches through pull requests. It is synchronised with
+the main development repository that is available on the
+[ICC](http://icc.dur.ac.uk)'s GitLab server which is available
+[here](https://gitlab.cosma.dur.ac.uk/swift/swiftsim).
+
+Please feel free to submit issues to this repository, or even pull
+requests. We will try to deal with them as soon as possible, but as the
+core development team is quite small this could take some time.
+
+Contribution Guidelines
+-----------------------
+
+The SWIFT source code uses a variation of the 'Google' formatting 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 pull request.
+
+Please check that the test suite still runs with your changes applied before
+submitting a pull request and add relevant unit tests probing the correctness
+of new modules. An example of how to add a test to the suite can be found by
+considering the tests/testGreeting case.
+
+Any contributions that fail any of the automated tests will not be accepted.
+Contributions that include tests of the proposed modules (or any current ones!)
+are highly encouraged.
+
+```
+ Welcome to the cosmological hydrodynamical code
+ ______ _________________
+ / ___/ | / / _/ ___/_ __/
+ \__ \| | /| / // // /_ / /
+ ___/ /| |/ |/ // // __/ / /
+ /____/ |__/|__/___/_/ /_/
+ SPH With Inter-dependent Fine-grained Tasking
+
+ Website: www.swiftsim.com
+ Twitter: @SwiftSimulation
+
+See INSTALL.swift for install instructions.
+
+Usage: swift [OPTION]... PARAMFILE
+ swift_mpi [OPTION]... PARAMFILE
+
+Valid options are:
+ -a Pin runners using processor affinity.
+ -c Run with cosmological time integration.
+ -C Run with cooling.
+ -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.
+ Allows user to check validity of parameter and IC files as well as memory limits.
+ -D Always drift all particles even the ones far from active particles. This emulates
+ Gadget-[23] and GIZMO's default behaviours.
+ -e Enable floating-point exceptions (debugging mode).
+ -f {int} Overwrite the CPU frequency (Hz) to be used for time measurements.
+ -g Run with an external gravitational potential.
+ -G Run with self-gravity.
+ -M Reconstruct the multipoles every time-step.
+ -n {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop.
+ -P {sec:par:val} Set parameter value and overwrites values read from the parameters file. Can be used more than once.
+ -s Run with hydrodynamics.
+ -S Run with stars.
+ -t {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
+ -T Print timers every time-step.
+ -v [12] Increase the level of verbosity:
+ 1: MPI-rank 0 writes,
+ 2: All MPI-ranks write.
+ -y {int} Time-step frequency at which task graphs are dumped.
+ -Y {int} Time-step frequency at which threadpool tasks are dumped.
+ -h Print this help message and exit.
+
+See the file examples/parameter_example.yml for an example of parameter file.
+```
diff --git a/examples/EAGLE_100/eagle_100.yml b/examples/EAGLE_100/eagle_100.yml
index 9cf1c9432ed6c79e4b91262cfdd30bad14833558..b3cbbccc7a89131fedc27aa1024ed8f460494582 100644
--- a/examples/EAGLE_100/eagle_100.yml
+++ b/examples/EAGLE_100/eagle_100.yml
@@ -37,15 +37,18 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- epsilon: 0.001 # Softening length (in internal units).
- theta: 0.85 # Opening angle (Multipole acceptance criterion)
-
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.85 # Opening angle (Multipole acceptance criterion)
+ comoving_softening: 0.0026994 # Comoving softening length (in internal units).
+ max_physical_softening: 0.0007 # Physical softening length (in internal units).
+
# 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).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100 # (internal units)
# Parameters related to the initial conditions
InitialConditions:
file_name: ./EAGLE_ICs_100.hdf5 # The file to read
+ cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
diff --git a/examples/EAGLE_100/run.sh b/examples/EAGLE_100/run.sh
index 6ef47d5d98172cc8a318242923ede37332bd5590..9310fc370acb51256b46e6ba0fb739fb2728caf9 100755
--- a/examples/EAGLE_100/run.sh
+++ b/examples/EAGLE_100/run.sh
@@ -7,5 +7,5 @@ then
./getIC.sh
fi
-../swift -s -t 16 eagle_100.yml 2>&1 | tee output.log
+../swift -c -s -t 16 eagle_100.yml 2>&1 | tee output.log
diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index d7aedf591a6988f89ad700f0a6729b3682bb6881..933f085ccd0bcbaf80a6ca83b7cb06c5763d4a20 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -38,16 +38,18 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- epsilon: 0.001 # Softening length (in internal units).
- theta: 0.85 # Opening angle (Multipole acceptance criterion)
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.85 # Opening angle (Multipole acceptance criterion)
+ comoving_softening: 0.0026994 # Comoving softening length (in internal units).
+ max_physical_softening: 0.0007 # Physical softening length (in internal units).
# 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).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100 # (internal units)
# Parameters related to the initial conditions
InitialConditions:
file_name: ./EAGLE_ICs_12.hdf5 # The file to read
-
+ cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
diff --git a/examples/EAGLE_12/run.sh b/examples/EAGLE_12/run.sh
index 21b965050ab1c786fd08c0dfd28b1612ac09e9e5..82d0e3e4d4418f0d00851f1279587a1ed45144df 100755
--- a/examples/EAGLE_12/run.sh
+++ b/examples/EAGLE_12/run.sh
@@ -7,5 +7,5 @@ then
./getIC.sh
fi
-../swift -s -t 16 eagle_12.yml 2>&1 | tee output.log
+../swift -c -s -t 16 eagle_12.yml 2>&1 | tee output.log
diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml
index 24daa13e06c2058d8a93243e9838b6f6b2cfe9bb..37ad91a29442825fcc843084cbbc4224203cab87 100644
--- a/examples/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_25/eagle_25.yml
@@ -38,16 +38,19 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- epsilon: 0.001 # Softening length (in internal units).
- theta: 0.85 # Opening angle (Multipole acceptance criterion)
-
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.85 # Opening angle (Multipole acceptance criterion)
+ comoving_softening: 0.0026994 # Comoving softening length (in internal units).
+ max_physical_softening: 0.0007 # Physical softening length (in internal units).
+
# 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).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100 # (internal units)
# Parameters related to the initial conditions
InitialConditions:
file_name: ./EAGLE_ICs_25.hdf5 # The file to read
+ cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
diff --git a/examples/EAGLE_25/run.sh b/examples/EAGLE_25/run.sh
index 21d0c7baf4482bdc6df83c552253a02d47432ba2..ae36e0d8e6ae6377682f259f506f720fbd585b29 100755
--- a/examples/EAGLE_25/run.sh
+++ b/examples/EAGLE_25/run.sh
@@ -7,5 +7,5 @@ then
./getIC.sh
fi
-../swift -s -t 16 eagle_25.yml 2>&1 | tee output.log
+../swift -c -s -t 16 eagle_25.yml 2>&1 | tee output.log
diff --git a/examples/EAGLE_50/eagle_50.yml b/examples/EAGLE_50/eagle_50.yml
index adf64f6963e2a0bc545637757f5826efd0c027a9..021e6b8ab347d38b7ecaeddd5a198324ab57f31b 100644
--- a/examples/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_50/eagle_50.yml
@@ -37,16 +37,19 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- epsilon: 0.001 # Softening length (in internal units).
- theta: 0.85 # Opening angle (Multipole acceptance criterion)
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.85 # Opening angle (Multipole acceptance criterion)
+ comoving_softening: 0.0026994 # Comoving softening length (in internal units).
+ max_physical_softening: 0.0007 # Physical softening length (in internal units).
# 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).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100 # (internal units)
# Parameters related to the initial conditions
InitialConditions:
file_name: ./EAGLE_ICs_50.hdf5 # The file to read
+ cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
diff --git a/examples/EAGLE_50/run.sh b/examples/EAGLE_50/run.sh
index 253e7e6eb11a04e1731c605e9ea067ba2bd13221..c3b93e84192ff58e445d4e164a4d6f66d19f85fe 100755
--- a/examples/EAGLE_50/run.sh
+++ b/examples/EAGLE_50/run.sh
@@ -7,5 +7,5 @@ then
./getIC.sh
fi
-../swift -s -t 16 eagle_50.yml 2>&1 | tee output.log
+../swift -c -s -t 16 eagle_50.yml 2>&1 | tee output.log
diff --git a/examples/EAGLE_6/eagle_6.yml b/examples/EAGLE_6/eagle_6.yml
index d0055607d3cd7c1453e3b93ed8b21cd5ff2673d5..839ccf7e429248f7b7eb81635e20e43df7d1af9d 100644
--- a/examples/EAGLE_6/eagle_6.yml
+++ b/examples/EAGLE_6/eagle_6.yml
@@ -28,7 +28,7 @@ TimeIntegration:
# Parameters governing the snapshots
Snapshots:
basename: eagle # Common part of the name of output files
- time_first: 0. # Time of the first output (in internal units)
+ time_first: 1. # Time of the first output (in internal units)
delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
compression: 4
@@ -38,16 +38,19 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- epsilon: 0.001 # Softening length (in internal units).
- theta: 0.85 # Opening angle (Multipole acceptance criterion)
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.85 # Opening angle (Multipole acceptance criterion)
+ comoving_softening: 0.0026994 # Comoving softening length (in internal units).
+ max_physical_softening: 0.0007 # Physical softening length (in internal units).
# 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).
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100 # (internal units)
# Parameters related to the initial conditions
InitialConditions:
file_name: ./EAGLE_ICs_6.hdf5 # The file to read
+ cleanup_h_factors: 1 # Remove the h-factors inherited from Gadget
diff --git a/examples/EAGLE_6/run.sh b/examples/EAGLE_6/run.sh
index d8e5592467a115460bb455ab31bb5e1f4017a948..3e1d21954d8402c9e1929119db50708859332b6c 100755
--- a/examples/EAGLE_6/run.sh
+++ b/examples/EAGLE_6/run.sh
@@ -7,5 +7,5 @@ then
./getIC.sh
fi
-../swift -s -t 16 eagle_6.yml 2>&1 | tee output.log
+../swift -c -s -t 16 eagle_6.yml 2>&1 | tee output.log
diff --git a/examples/EvrardCollapse_3D/evrard.yml b/examples/EvrardCollapse_3D/evrard.yml
index 94cd6e7a95d2671f7021ae3cd2ba40477cd8e7fd..a2b39fb1e7fa36098abe720945ece4b611eb4fe8 100644
--- a/examples/EvrardCollapse_3D/evrard.yml
+++ b/examples/EvrardCollapse_3D/evrard.yml
@@ -31,11 +31,10 @@ SPH:
# Parameters for the self-gravity scheme
Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- epsilon: 0.001 # Softening length (in internal units).
- theta: 0.9
- a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
- r_cut: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.9
+ comoving_softening: 0.001 # Comoving softening length (in internal units).
+ max_physical_softening: 0.001 # Physical softening length (in internal units).
# Parameters related to the initial conditions
InitialConditions:
diff --git a/examples/EvrardCollapse_3D/makeIC.py b/examples/EvrardCollapse_3D/makeIC.py
index b68f0da0e1869ce54ae7ab2ad3b33b8c3deb8b51..f4d3c4c5bf7f91e5f79cfcd4e9ae23388932144e 100644
--- a/examples/EvrardCollapse_3D/makeIC.py
+++ b/examples/EvrardCollapse_3D/makeIC.py
@@ -18,9 +18,24 @@
################################################################################
import h5py
+import argparse as ap
from numpy import *
-# Generates a swift IC file for the Evrard collapse
+parser = ap.ArgumentParser(
+ description="Generates a swift IC file for the Evrard collapse"
+)
+
+parser.add_argument(
+ "-n",
+ "--nparts",
+ help="""
+ Number of particles to be used in the Evrard collapse.
+ """,
+ required=False,
+ default=100000
+)
+
+args = vars(parser.parse_args())
# Parameters
gamma = 5. / 3. # Gas adiabatic index
@@ -28,7 +43,7 @@ M = 1. # total mass of the sphere
R = 1. # radius of the sphere
u0 = 0.05 / M # initial thermal energy
fileName = "evrard.hdf5"
-numPart = 100000
+numPart = int(args["nparts"])
r = R * sqrt(random.random(numPart))
phi = 2. * pi * random.random(numPart)
diff --git a/examples/SedovBlast_3D/sedov.yml b/examples/SedovBlast_3D/sedov.yml
index 057e5a18acba2a1e8211d40f15ffde35b2019a9b..2df2c432cef8afec49c687b643a872b06f9abb60 100644
--- a/examples/SedovBlast_3D/sedov.yml
+++ b/examples/SedovBlast_3D/sedov.yml
@@ -31,6 +31,6 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
- file_name: ./sedov.hdf5 # The file to read
- h_scaling: 3.33
+ file_name: ./sedov.hdf5
+ smoothing_length_scaling: 3.33
diff --git a/examples/SquareTest_2D/run.sh b/examples/SquareTest_2D/run.sh
index 242f1b8b729979b026cfe31002e84cd9ef741129..7d77e9c5bd89732970b47feb3a297ef92b345a01 100755
--- a/examples/SquareTest_2D/run.sh
+++ b/examples/SquareTest_2D/run.sh
@@ -11,4 +11,4 @@ fi
../swift -s -t 1 square.yml 2>&1 | tee output.log
# Plot the solution
-python plotSolution.py 40
+python plotSolution.py 5
diff --git a/examples/analyse_dump_cells.py b/examples/analyse_dump_cells.py
index 3140e799566c75fe494a75895db0f4a8dcff4e57..2adfaf319e9c0da33f86a6158da68e6620c47361 100755
--- a/examples/analyse_dump_cells.py
+++ b/examples/analyse_dump_cells.py
@@ -39,9 +39,11 @@ zcol = 2
xwcol = 3
ywcol = 4
zwcol = 5
-localcol = 18
supercol = 15
-activecol = 16
+topcol = 16
+activecol = 17
+localcol = 19
+mpicol = 20
# Command-line arguments.
if len(sys.argv) < 5:
@@ -59,11 +61,12 @@ for i in range(4, len(sys.argv)):
# Read the file.
data = pl.loadtxt(sys.argv[i])
- #print data
+ if len(data) == 0 or len(data) == 20:
+ continue
- # Select cells that are on the current rank and are super cells.
+ # Select cells that are on the current rank and are top-level cells.
rdata = data[data[:,localcol] == 1]
- tdata = rdata[rdata[:,supercol] == 1]
+ tdata = rdata[rdata[:,topcol] == 1]
# Separation of the cells is in data.
xwidth = tdata[0,xwcol]
diff --git a/examples/main.c b/examples/main.c
index 80e4dfb479172d9e2bca6fb19ee676371d130afe..ca91c9e718cca9f7f8d82fce2f5d237f9f5b98e4 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -484,7 +484,7 @@ int main(int argc, char *argv[]) {
#endif
/* Common variables for restart and IC sections. */
- int clean_h_values = 0;
+ int clean_smoothing_length_values = 0;
int flag_entropy_ICs = 0;
/* Work out where we will read and write restart files. */
@@ -609,20 +609,32 @@ int main(int argc, char *argv[]) {
if (myrank == 0 && with_cosmology) cosmology_print(&cosmo);
/* Initialise the hydro properties */
- if (with_hydro) hydro_props_init(&hydro_properties, params);
+ if (with_hydro)
+ hydro_props_init(&hydro_properties, &prog_const, &us, params);
if (with_hydro) eos_init(&eos, params);
/* Initialise the gravity properties */
- if (with_self_gravity) gravity_props_init(&gravity_properties, params);
+ if (with_self_gravity)
+ gravity_props_init(&gravity_properties, params, &cosmo);
/* Read particles and space information from (GADGET) ICs */
char ICfileName[200] = "";
parser_get_param_string(params, "InitialConditions:file_name", ICfileName);
const int replicate =
parser_get_opt_param_int(params, "InitialConditions:replicate", 1);
- clean_h_values =
- parser_get_opt_param_int(params, "InitialConditions:cleanup_h", 0);
+ clean_smoothing_length_values = parser_get_opt_param_int(
+ params, "InitialConditions:cleanup_smoothing_lengths", 0);
+ const int cleanup_h = parser_get_opt_param_int(
+ params, "InitialConditions:cleanup_h_factors", 0);
+ const int generate_gas_in_ics = parser_get_opt_param_int(
+ params, "InitialConditions:generate_gas_in_ics", 0);
+ if (generate_gas_in_ics && flag_entropy_ICs)
+ error("Can't generate gas if the entropy flag is set in the ICs.");
+ if (generate_gas_in_ics && !with_cosmology)
+ error("Can't generate gas if the run is not cosmological.");
if (myrank == 0) message("Reading ICs from file '%s'", ICfileName);
+ if (myrank == 0 && cleanup_h)
+ message("Cleaning up h-factors (h=%f)", cosmo.h);
fflush(stdout);
/* Get ready to read particles of all kinds */
@@ -630,25 +642,27 @@ int main(int argc, char *argv[]) {
double dim[3] = {0., 0., 0.};
int periodic = 0;
if (myrank == 0) clocks_gettime(&tic);
+#if defined(HAVE_HDF5)
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
read_ic_parallel(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
&Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
- myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL,
- nr_threads, dry_run);
+ cleanup_h, cosmo.h, myrank, nr_nodes, MPI_COMM_WORLD,
+ MPI_INFO_NULL, nr_threads, dry_run);
#else
read_ic_serial(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
&Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
- myrank, nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL, nr_threads,
- dry_run);
+ cleanup_h, cosmo.h, myrank, nr_nodes, MPI_COMM_WORLD,
+ MPI_INFO_NULL, nr_threads, dry_run);
#endif
#else
read_ic_single(ICfileName, &us, dim, &parts, &gparts, &sparts, &Ngas,
&Ngpart, &Nspart, &periodic, &flag_entropy_ICs, with_hydro,
(with_external_gravity || with_self_gravity), with_stars,
- nr_threads, dry_run);
+ cleanup_h, cosmo.h, nr_threads, dry_run);
+#endif
#endif
if (myrank == 0) {
clocks_gettime(&toc);
@@ -659,14 +673,18 @@ int main(int argc, char *argv[]) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check once and for all that we don't have unwanted links */
- if (!with_stars) {
+ if (!with_stars && !dry_run) {
for (size_t k = 0; k < Ngpart; ++k)
if (gparts[k].type == swift_type_star) error("Linking problem");
}
- if (!with_hydro) {
+ if (!with_hydro && !dry_run) {
for (size_t k = 0; k < Ngpart; ++k)
if (gparts[k].type == swift_type_gas) error("Linking problem");
}
+
+ /* Check that the other links are correctly set */
+ if (!dry_run)
+ part_verify_links(parts, gparts, sparts, Ngas, Ngpart, Nspart, 1);
#endif
/* Get the total number of particles across all nodes. */
@@ -690,8 +708,9 @@ int main(int argc, char *argv[]) {
/* Initialize the space with these data. */
if (myrank == 0) clocks_gettime(&tic);
- space_init(&s, params, dim, parts, gparts, sparts, Ngas, Ngpart, Nspart,
- periodic, replicate, with_self_gravity, talking, dry_run);
+ space_init(&s, params, &cosmo, dim, parts, gparts, sparts, Ngas, Ngpart,
+ Nspart, periodic, replicate, generate_gas_in_ics,
+ with_self_gravity, talking, dry_run);
if (myrank == 0) {
clocks_gettime(&toc);
@@ -700,6 +719,11 @@ int main(int argc, char *argv[]) {
fflush(stdout);
}
+ /* Check that the matter content matches the cosmology given in the
+ * parameter file. */
+ if (with_cosmology && with_self_gravity && !dry_run)
+ space_check_cosmology(&s, &cosmo, myrank);
+
/* Also update the total counts (in case of changes due to replication) */
#if defined(WITH_MPI)
N_long[0] = s.nr_parts;
@@ -819,7 +843,7 @@ int main(int argc, char *argv[]) {
/* Initialise the table of Ewald corrections for the gravity checks */
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
- if (periodic) gravity_exact_force_ewald_init(e.s->dim[0]);
+ if (s.periodic) gravity_exact_force_ewald_init(e.s->dim[0]);
#endif
/* Init the runner history. */
@@ -836,7 +860,7 @@ int main(int argc, char *argv[]) {
#endif
/* Initialise the particles */
- engine_init_particles(&e, flag_entropy_ICs, clean_h_values);
+ engine_init_particles(&e, flag_entropy_ICs, clean_smoothing_length_values);
/* Write the state of the system before starting time integration. */
engine_dump_snapshot(&e);
@@ -845,10 +869,10 @@ int main(int argc, char *argv[]) {
/* Legend */
if (myrank == 0)
- printf("# %6s %14s %14s %14s %9s %12s %12s %12s %16s [%s] %6s\n", "Step",
- "Time", "Scale-factor", "Time-step", "Time-bins", "Updates",
- "g-Updates", "s-Updates", "Wall-clock time", clocks_getunit(),
- "Props");
+ printf("# %6s %14s %14s %10s %14s %9s %12s %12s %12s %16s [%s] %6s\n",
+ "Step", "Time", "Scale-factor", "Redshift", "Time-step", "Time-bins",
+ "Updates", "g-Updates", "s-Updates", "Wall-clock time",
+ clocks_getunit(), "Props");
/* File for the timers */
if (with_verbose_timers) timers_open_file(myrank);
@@ -1000,6 +1024,24 @@ int main(int argc, char *argv[]) {
(double)runner_hist_bins[k]);
#endif
+ /* Write final time information */
+ if (myrank == 0) {
+
+ /* Print some information to the screen */
+ printf(" %6d %14e %14e %10.5f %14e %4d %4d %12zu %12zu %12zu %21.3f %6d\n",
+ e.step, e.time, e.cosmology->a, e.cosmology->z, e.time_step,
+ e.min_active_bin, e.max_active_bin, e.updates, e.g_updates,
+ e.s_updates, e.wallclock_time, e.step_props);
+ fflush(stdout);
+
+ fprintf(e.file_timesteps,
+ " %6d %14e %14e %14e %4d %4d %12zu %12zu %12zu %21.3f %6d\n",
+ e.step, e.time, e.cosmology->a, e.time_step, e.min_active_bin,
+ e.max_active_bin, e.updates, e.g_updates, e.s_updates,
+ e.wallclock_time, e.step_props);
+ fflush(e.file_timesteps);
+ }
+
/* Write final output. */
engine_drift_all(&e);
engine_print_stats(&e);
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 4784a877a21e698af28f72aab473ac9c3cad605e..0a81716a3e34bbb305916b1148a9d8199b84fa5a 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -30,15 +30,19 @@ SPH:
h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step.
max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
-
+ initial_temperature: 0 # (Optional) Initial temperature (in internal units) to set the gas particles at start-up. Value is ignored if set to 0.
+ minimal_temperature: 0 # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
+ H_mass_fraction: 0.76 # (Optional) Hydrogen mass fraction used for initial conversion from temp to internal energy.
+
# Parameters for the self-gravity scheme
Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 0.7 # Opening angle (Multipole acceptance criterion)
- epsilon: 0.1 # Softening length (in internal units).
- a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
- r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
- r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
+ comoving_softening: 0.0026994 # Comoving softening length (in internal units).
+ max_physical_softening: 0.0007 # Physical softening length (in internal units).
+ a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
+ r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+ r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
# Parameters for the task scheduling
Scheduler:
@@ -81,12 +85,14 @@ Statistics:
# Parameters related to the initial conditions
InitialConditions:
file_name: SedovBlast/sedov.hdf5 # The file to read
- cleanup_h: 0 # (Optional) Clean the values of h that are read in. Set to 1 to activate.
- h_scaling: 1. # (Optional) A scaling factor to apply to all smoothing lengths in the ICs.
+ generate_gas_in_ics: 0 # (Optional) Generate gas particles from the DM-only ICs (e.g. from panphasia).
+ cleanup_h_factors: 0 # (Optional) Clean up the h-factors used in the ICs (e.g. in Gadget files).
+ cleanup_smoothing_lengths: 0 # (Optional) Clean the values of the smoothing lengths that are read in to remove stupid values. Set to 1 to activate.
+ smoothing_length_scaling: 1. # (Optional) A scaling factor to apply to all smoothing lengths in the ICs.
shift_x: 0. # (Optional) A shift to apply to all particles read from the ICs (in internal units).
shift_y: 0.
shift_z: 0.
- replicate: 2 # (Optional) Replicate all particles along each axis a given number of times. Default 1.
+ replicate: 2 # (Optional) Replicate all particles along each axis a given integer number of times. Default 1.
# Parameters controlling restarts
Restarts:
diff --git a/examples/process_cells b/examples/process_cells
index eead4572387f2732cf0b86265f14d20039f73ae1..b57ed9e73c76807707d158485e31a5b643d8dec0 100755
--- a/examples/process_cells
+++ b/examples/process_cells
@@ -48,7 +48,7 @@ echo "Number of steps = $steps"
# And process them,
echo "Processing cell dumps files..."
-#echo $files | xargs -P $NPROCS -n 4 /bin/bash -c "${SCRIPTHOME}/process_cells_helper $NX $NY $NZ \$0 \$1 \$2 \$3"
+echo $files | xargs -P $NPROCS -n 4 /bin/bash -c "${SCRIPTHOME}/process_cells_helper $NX $NY $NZ \$0 \$1 \$2 \$3"
# Create summary.
grep "top cells" step*-active-cells.dat | sort -h > active_cells.log
diff --git a/src/Makefile.am b/src/Makefile.am
index 6877cbe71f45d3efc5b77c72b70d8d0722bf3e13..a7bc4774ea3fb162a2b3af1ea825607659e34c6a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -65,11 +65,13 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h runner_doiact_fft.h \
runner_doiact_nosort.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
- dimension.h equation_of_state.h part_type.h periodic.h memswap.h dump.h logger.h \
+ dimension.h part_type.h periodic.h memswap.h dump.h logger.h \
gravity.h gravity_io.h gravity_cache.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 \
sourceterms.h \
+ equation_of_state.h \
+ equation_of_state/ideal_gas/equation_of_state.h equation_of_state/isothermal/equation_of_state.h \
hydro.h hydro_io.h \
hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
diff --git a/src/approx_math.h b/src/approx_math.h
index 48319ddfd7a86c132a1cd18b4a08fa849a36a15a..aa9ed2b4efa6e0542e2eb2432132f4b0232f7403 100644
--- a/src/approx_math.h
+++ b/src/approx_math.h
@@ -49,4 +49,16 @@ __attribute__((always_inline)) INLINE static float good_approx_expf(float x) {
x * ((1.f / 120.f) + (1.f / 720.f) * x)))));
}
+/**
+ * @brief Approximate version of exp(x) using a 6th order Taylor expansion
+ */
+__attribute__((always_inline)) INLINE static double good_approx_exp(double x) {
+ return 1. +
+ x * (1. +
+ x * (0.5 +
+ x * ((1. / 6.) +
+ x * ((1. / 24.) +
+ x * ((1. / 120.) + (1. / 720.) * x)))));
+}
+
#endif /* SWIFT_APPROX_MATH_H */
diff --git a/src/common_io.c b/src/common_io.c
index 67326817397118568e060f5dff49421bb32fba4d..c235060b1008aef0aa157bfe6e41dea76357b7b5 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -21,21 +21,6 @@
/* Config parameters. */
#include "../config.h"
-#if defined(HAVE_HDF5)
-
-/* Some standard headers. */
-#include <hdf5.h>
-#include <math.h>
-#include <stddef.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-/* MPI headers. */
-#ifdef WITH_MPI
-#include <mpi.h>
-#endif
-
/* This object's header. */
#include "common_io.h"
@@ -50,6 +35,22 @@
#include "units.h"
#include "version.h"
+/* Some standard headers. */
+#include <math.h>
+#include <stddef.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+#if defined(HAVE_HDF5)
+
+#include <hdf5.h>
+
+/* MPI headers. */
+#ifdef WITH_MPI
+#include <mpi.h>
+#endif
+
/**
* @brief Converts a C data type to the HDF5 equivalent.
*
@@ -84,36 +85,6 @@ hid_t io_hdf5_type(enum IO_DATA_TYPE type) {
}
}
-/**
- * @brief Returns the memory size of the data type
- */
-size_t io_sizeof_type(enum IO_DATA_TYPE type) {
-
- switch (type) {
- case INT:
- return sizeof(int);
- case UINT:
- return sizeof(unsigned int);
- case LONG:
- return sizeof(long);
- case ULONG:
- return sizeof(unsigned long);
- case LONGLONG:
- return sizeof(long long);
- case ULONGLONG:
- return sizeof(unsigned long long);
- case FLOAT:
- return sizeof(float);
- case DOUBLE:
- return sizeof(double);
- case CHAR:
- return sizeof(char);
- default:
- error("Unknown type");
- return 0;
- }
-}
-
/**
* @brief Return 1 if the type has double precision
*
@@ -429,6 +400,36 @@ void io_write_engine_policy(hid_t h_file, const struct engine* e) {
#endif /* HAVE_HDF5 */
+/**
+ * @brief Returns the memory size of the data type
+ */
+size_t io_sizeof_type(enum IO_DATA_TYPE type) {
+
+ switch (type) {
+ case INT:
+ return sizeof(int);
+ case UINT:
+ return sizeof(unsigned int);
+ case LONG:
+ return sizeof(long);
+ case ULONG:
+ return sizeof(unsigned long);
+ case LONGLONG:
+ return sizeof(long long);
+ case ULONGLONG:
+ return sizeof(unsigned long long);
+ case FLOAT:
+ return sizeof(float);
+ case DOUBLE:
+ return sizeof(double);
+ case CHAR:
+ return sizeof(char);
+ default:
+ error("Unknown type");
+ return 0;
+ }
+}
+
/**
* @brief Mapper function to copy #part or #gpart fields into a buffer.
*/
@@ -457,6 +458,7 @@ void io_convert_part_f_mapper(void* restrict temp, int N,
const struct io_props props = *((const struct io_props*)extra_data);
const struct part* restrict parts = props.parts;
+ const struct xpart* restrict xparts = props.xparts;
const struct engine* e = props.e;
const size_t dim = props.dimension;
@@ -465,7 +467,8 @@ void io_convert_part_f_mapper(void* restrict temp, int N,
const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
for (int i = 0; i < N; i++)
- props.convert_part_f(e, parts + delta + i, &temp_f[i * dim]);
+ props.convert_part_f(e, parts + delta + i, xparts + delta + i,
+ &temp_f[i * dim]);
}
/**
@@ -477,6 +480,7 @@ void io_convert_part_d_mapper(void* restrict temp, int N,
const struct io_props props = *((const struct io_props*)extra_data);
const struct part* restrict parts = props.parts;
+ const struct xpart* restrict xparts = props.xparts;
const struct engine* e = props.e;
const size_t dim = props.dimension;
@@ -485,7 +489,8 @@ void io_convert_part_d_mapper(void* restrict temp, int N,
const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
for (int i = 0; i < N; i++)
- props.convert_part_d(e, parts + delta + i, &temp_d[i * dim]);
+ props.convert_part_d(e, parts + delta + i, xparts + delta + i,
+ &temp_d[i * dim]);
}
/**
@@ -704,7 +709,7 @@ void io_duplicate_hydro_gparts_mapper(void* restrict data, int Ngas,
gparts[i + Ndm].type = swift_type_gas;
/* Link the particles */
- gparts[i + Ndm].id_or_neg_offset = -i;
+ gparts[i + Ndm].id_or_neg_offset = -(long long)(offset + i);
parts[i].gpart = &gparts[i + Ndm];
}
}
@@ -760,7 +765,7 @@ void io_duplicate_hydro_sparts_mapper(void* restrict data, int Nstars,
gparts[i + Ndm].type = swift_type_star;
/* Link the particles */
- gparts[i + Ndm].id_or_neg_offset = -i;
+ gparts[i + Ndm].id_or_neg_offset = -(long long)(offset + i);
sparts[i].gpart = &gparts[i + Ndm];
}
}
diff --git a/src/common_io.h b/src/common_io.h
index 49358d8a374553d803de144f1135df6eee80cf15..4f6d75fad0482f62b8d5c3f80a632d7eb6f54d60 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -37,8 +37,6 @@ struct io_props;
struct engine;
struct threadpool;
-#if defined(HAVE_HDF5)
-
/**
* @brief The different types of data used in the GADGET IC files.
*
@@ -56,9 +54,9 @@ enum IO_DATA_TYPE {
CHAR
};
+#if defined(HAVE_HDF5)
+
hid_t io_hdf5_type(enum IO_DATA_TYPE type);
-size_t io_sizeof_type(enum IO_DATA_TYPE type);
-int io_is_double_precision(enum IO_DATA_TYPE type);
void io_read_attribute(hid_t grp, const char* name, enum IO_DATA_TYPE type,
void* data);
@@ -86,6 +84,9 @@ void io_copy_temp_buffer(void* temp, const struct engine* e,
#endif /* defined HDF5 */
+size_t io_sizeof_type(enum IO_DATA_TYPE type);
+int io_is_double_precision(enum IO_DATA_TYPE type);
+
void io_collect_dm_gparts(const struct gpart* const gparts, size_t Ntot,
struct gpart* const dmparts, size_t Ndm);
void io_prepare_dm_gparts(struct threadpool* tp, struct gpart* const gparts,
diff --git a/src/cosmology.c b/src/cosmology.c
index b1bf71ad25451d2185d911320fb757193e990945..a5034240d07b9b9bf401d470c0bacbf6962a8006 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -130,9 +130,11 @@ double cosmology_get_time_since_big_bang(const struct cosmology *c, double a) {
* @brief Update the cosmological parameters to the current simulation time.
*
* @param c The #cosmology struct.
+ * @param phys_const The physical constants in the internal units.
* @param ti_current The current (integer) time.
*/
-void cosmology_update(struct cosmology *c, integertime_t ti_current) {
+void cosmology_update(struct cosmology *c, const struct phys_const *phys_const,
+ integertime_t ti_current) {
/* Get scale factor and powers of it */
const double a = c->a_begin * exp(ti_current * c->time_base);
@@ -173,6 +175,10 @@ void cosmology_update(struct cosmology *c, integertime_t ti_current) {
/* Expansion rate */
c->a_dot = c->H * c->a;
+ /* Critical density */
+ c->critical_density =
+ 3. * c->H * c->H / (8. * M_PI * phys_const->const_newton_G);
+
/* Time-step conversion factor */
c->time_step_factor = c->H;
@@ -284,8 +290,6 @@ double time_integrand(double a, void *param) {
*/
void cosmology_init_tables(struct cosmology *c) {
- const ticks tic = getticks();
-
#ifdef HAVE_LIBGSL
/* Retrieve some constants */
@@ -373,9 +377,6 @@ void cosmology_init_tables(struct cosmology *c) {
error("Code not compiled with GSL. Can't compute cosmology integrals.");
#endif
-
- message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
- clocks_getunit());
}
/**
@@ -431,7 +432,7 @@ void cosmology_init(const struct swift_params *params,
cosmology_init_tables(c);
/* Set remaining variables to alid values */
- cosmology_update(c, 0);
+ cosmology_update(c, phys_const, 0);
/* Update the times */
c->time_begin = cosmology_get_time_since_big_bang(c, c->a_begin);
@@ -454,7 +455,7 @@ void cosmology_init_no_cosmo(struct cosmology *c) {
c->Omega_b = 0.;
c->w_0 = 0.;
c->w_a = 0.;
- c->h = 0.;
+ c->h = 1.;
c->w = 0.;
c->a_begin = 1.;
@@ -474,6 +475,8 @@ void cosmology_init_no_cosmo(struct cosmology *c) {
c->a_factor_hydro_accel = 1.;
c->a_factor_grav_accel = 1.;
+ c->critical_density = 0.;
+
c->time_step_factor = 1.;
c->a_dot = 0.;
diff --git a/src/cosmology.h b/src/cosmology.h
index b541e20fa1dd4f55422f9ff1d0952783d8de0ad7..109b80a57d8dbc4eb942dd4ecbbc0db84198100b 100644
--- a/src/cosmology.h
+++ b/src/cosmology.h
@@ -70,6 +70,9 @@ struct cosmology {
/*! Hubble constant at the current redshift (in internal units) */
double H;
+ /*! The critical density at the current redshift (in internal units) */
+ double critical_density;
+
/*! Conversion factor from internal time-step size to cosmological step */
double time_step_factor;
@@ -160,7 +163,8 @@ struct cosmology {
double universe_age_at_present_day;
};
-void cosmology_update(struct cosmology *c, integertime_t ti_current);
+void cosmology_update(struct cosmology *c, const struct phys_const *phys_const,
+ integertime_t ti_current);
double cosmology_get_drift_factor(const struct cosmology *cosmo,
integertime_t ti_start, integertime_t ti_end);
diff --git a/src/debug.c b/src/debug.c
index 77d791a31787c7997a9a56e3204dac1160178976..b1e2cb08bc7fa99330da3d9c9382dbef81b3215a 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -271,17 +271,17 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
}
/**
- * @brief map function for dumping cells. In MPI mode locally active cells
- * only.
+ * @brief map function for dumping cells.
*/
static void dumpCells_map(struct cell *c, void *data) {
size_t *ldata = (size_t *)data;
FILE *file = (FILE *)ldata[0];
struct engine *e = (struct engine *)ldata[1];
float ntasks = c->nr_tasks;
- int active = (int)ldata[2];
- int mpiactive = (int)ldata[3];
- int pactive = (int)ldata[4];
+ int super = (int)ldata[2];
+ int active = (int)ldata[3];
+ int mpiactive = (int)ldata[4];
+ int pactive = (int)ldata[5];
#if SWIFT_DEBUG_CHECKS
/* The c->nr_tasks field does not include all the tasks. So let's check this
@@ -305,11 +305,13 @@ static void dumpCells_map(struct cell *c, void *data) {
/* Only cells with particles are dumped. */
if (c->count > 0 || c->gcount > 0 || c->scount > 0) {
-/* In MPI mode we may only output cells with foreign partners.
- * These define the edges of the partitions. */
+ /* In MPI mode we may only output cells with foreign partners.
+ * These define the edges of the partitions. */
+ int ismpiactive = 0;
#if WITH_MPI
+ ismpiactive = (c->send_xv != NULL);
if (mpiactive)
- mpiactive = (c->send_xv != NULL);
+ mpiactive = ismpiactive;
else
mpiactive = 1;
#else
@@ -322,9 +324,10 @@ static void dumpCells_map(struct cell *c, void *data) {
else
active = 1;
- /* So output local super cells that are active and have MPI tasks as
- * requested. */
- if (c->nodeID == e->nodeID && (c->super == c) && active && mpiactive) {
+ /* So output local super cells that are active and have MPI
+ * tasks as requested. */
+ if (c->nodeID == e->nodeID && (!super || (super && c->super == c)) &&
+ active && mpiactive) {
/* If requested we work out how many particles are active in this cell. */
int pactcount = 0;
@@ -342,12 +345,13 @@ static void dumpCells_map(struct cell *c, void *data) {
fprintf(file,
" %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6d %6d %6d %6d %6d %6d "
- "%6.1f %20lld %6d %6d %6d %6d %6d\n",
+ "%6.1f %20lld %6d %6d %6d %6d %6d %6d %6d\n",
c->loc[0], c->loc[1], c->loc[2], c->width[0], c->width[1],
c->width[2], e->step, c->count, c->gcount, c->scount, pactcount,
c->depth, ntasks, c->ti_hydro_end_min,
get_time_bin(c->ti_hydro_end_min), (c->super == c),
- cell_is_active_hydro(c, e), c->nodeID, c->nodeID == e->nodeID);
+ (c->parent == NULL), cell_is_active_hydro(c, e), c->nodeID,
+ c->nodeID == e->nodeID, ismpiactive);
}
}
}
@@ -359,6 +363,7 @@ static void dumpCells_map(struct cell *c, void *data) {
*
* @param prefix base output filename, result is written to
* %prefix%_%rank%_%step%.dat
+ * @param super just output the super cells.
* @param active just output active cells.
* @param mpiactive just output MPI active cells, i.e. those with foreign cells.
* @param pactive also output a count of active particles.
@@ -366,8 +371,8 @@ static void dumpCells_map(struct cell *c, void *data) {
* @param rank node ID of MPI rank, or 0 if not relevant.
* @param step the current engine step, or some unique integer.
*/
-void dumpCells(const char *prefix, int active, int mpiactive, int pactive,
- struct space *s, int rank, int step) {
+void dumpCells(const char *prefix, int super, int active, int mpiactive,
+ int pactive, struct space *s, int rank, int step) {
FILE *file = NULL;
@@ -379,17 +384,18 @@ void dumpCells(const char *prefix, int active, int mpiactive, int pactive,
/* Header. */
fprintf(file,
"# %6s %6s %6s %6s %6s %6s %6s %6s %6s %6s %6s %6s %6s "
- "%20s %6s %6s %6s %6s %6s\n",
+ "%20s %6s %6s %6s %6s %6s %6s %6s\n",
"x", "y", "z", "xw", "yw", "zw", "step", "count", "gcount", "scount",
"actcount", "depth", "tasks", "ti_end_min", "timebin", "issuper",
- "active", "rank", "local");
+ "istop", "active", "rank", "local", "mpiactive");
- size_t data[5];
+ size_t data[6];
data[0] = (size_t)file;
data[1] = (size_t)s->e;
- data[2] = (size_t)active;
- data[3] = (size_t)mpiactive;
- data[4] = (size_t)pactive;
+ data[2] = (size_t)super;
+ data[3] = (size_t)active;
+ data[4] = (size_t)mpiactive;
+ data[5] = (size_t)pactive;
space_map_cells_pre(s, 1, dumpCells_map, &data);
fclose(file);
}
diff --git a/src/debug.h b/src/debug.h
index 5a646948204cc809c66430bb01c9ee90e21d720a..1e482c05c5af2dfebff1a254018fb1802df6cc5d 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -36,8 +36,8 @@ void printgParticle_single(struct gpart *gp);
int checkSpacehmax(struct space *s);
int checkCellhdxmax(const struct cell *c, int *depth);
-void dumpCells(const char *prefix, int active, int mpiactive, int pactive,
- struct space *s, int rank, int step);
+void dumpCells(const char *prefix, int super, int active, int mpiactive,
+ int pactive, struct space *s, int rank, int step);
#if defined(WITH_MPI) && defined(HAVE_METIS)
#include "metis.h"
diff --git a/src/engine.c b/src/engine.c
index 3e63166d3465c42ff5db21dac7308deb9cf96935..6f8c0e9937ee40912964b88a3be69cacf060a0ae 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3654,7 +3654,7 @@ int engine_estimate_nr_tasks(struct engine *e) {
#endif
double ntasks = n1 * ntop + n2 * (ncells - ntop);
- tasks_per_cell = ceil(ntasks / ncells);
+ if (ncells > 0) tasks_per_cell = ceil(ntasks / ncells);
if (tasks_per_cell < 1.0) tasks_per_cell = 1.0;
if (e->verbose)
@@ -3668,10 +3668,10 @@ int engine_estimate_nr_tasks(struct engine *e) {
* @brief Rebuild the space and tasks.
*
* @param e The #engine.
- * @param clean_h_values Are we cleaning up the values of h before building
- * the tasks ?
+ * @param clean_smoothing_length_values Are we cleaning up the values of
+ * the smoothing lengths before building the tasks ?
*/
-void engine_rebuild(struct engine *e, int clean_h_values) {
+void engine_rebuild(struct engine *e, int clean_smoothing_length_values) {
const ticks tic = getticks();
@@ -3687,7 +3687,7 @@ void engine_rebuild(struct engine *e, int clean_h_values) {
#endif
/* Initial cleaning up session ? */
- if (clean_h_values) space_sanitize(e->s);
+ if (clean_smoothing_length_values) space_sanitize(e->s);
/* If in parallel, exchange the cell structure, top-level and neighbouring
* multipoles. */
@@ -4379,10 +4379,10 @@ void engine_step(struct engine *e) {
if (e->nodeID == 0) {
/* Print some information to the screen */
- printf(" %6d %14e %14e %14e %4d %4d %12zu %12zu %12zu %21.3f %6d\n",
- e->step, e->time, e->cosmology->a, e->time_step, e->min_active_bin,
- e->max_active_bin, e->updates, e->g_updates, e->s_updates,
- e->wallclock_time, e->step_props);
+ printf(" %6d %14e %14e %10.5f %14e %4d %4d %12zu %12zu %12zu %21.3f %6d\n",
+ e->step, e->time, e->cosmology->a, e->cosmology->z, e->time_step,
+ e->min_active_bin, e->max_active_bin, e->updates, e->g_updates,
+ e->s_updates, e->wallclock_time, e->step_props);
fflush(stdout);
fprintf(e->file_timesteps,
@@ -4403,7 +4403,7 @@ void engine_step(struct engine *e) {
if (e->policy & engine_policy_cosmology) {
e->time_old = e->time;
- cosmology_update(e->cosmology, e->ti_current);
+ cosmology_update(e->cosmology, e->physical_constants, e->ti_current);
e->time = e->cosmology->time;
e->time_step = e->time - e->time_old;
} else {
@@ -4412,6 +4412,10 @@ void engine_step(struct engine *e) {
e->time_step = (e->ti_current - e->ti_old) * e->time_base;
}
+ /* Update the softening lengths */
+ if (e->policy & engine_policy_self_gravity)
+ gravity_update(e->gravity_properties, e->cosmology);
+
/* Prepare the tasks to be launched, rebuild or repartition if needed. */
engine_prepare(e);
@@ -4436,7 +4440,7 @@ void engine_step(struct engine *e) {
if (e->verbose) engine_print_task_counts(e);
/* Dump local cells and active particle counts. */
-/* dumpCells("cells", 0, 0, 1, e->s, e->nodeID, e->step); */
+/* dumpCells("cells", 0, 0, 0, 0, e->s, e->nodeID, e->step); */
#ifdef SWIFT_DEBUG_CHECKS
/* Check that we have the correct total mass in the top-level multipoles */
@@ -5093,6 +5097,7 @@ void engine_dump_snapshot(struct engine *e) {
#endif
/* Dump... */
+#if defined(HAVE_HDF5)
#if defined(WITH_MPI)
#if defined(HAVE_PARALLEL_HDF5)
write_output_parallel(e, e->snapshotBaseName, e->internal_units,
@@ -5106,6 +5111,7 @@ void engine_dump_snapshot(struct engine *e) {
#else
write_output_single(e, e->snapshotBaseName, e->internal_units,
e->snapshotUnits);
+#endif
#endif
e->dump_snapshot = 0;
@@ -5201,7 +5207,7 @@ void engine_init(
long long Ngas, long long Ndm, int policy, int verbose,
struct repartition *reparttype, const struct unit_system *internal_units,
const struct phys_const *physical_constants, struct cosmology *cosmo,
- const struct hydro_props *hydro, const struct gravity_props *gravity,
+ const struct hydro_props *hydro, struct gravity_props *gravity,
const struct external_potential *potential,
const struct cooling_function_data *cooling_func,
const struct chemistry_data *chemistry, struct sourceterms *sourceterms) {
diff --git a/src/engine.h b/src/engine.h
index 0fa8ca93b84760d0d92d9494f6e8e6224e3d1d9a..4c37aab593fdffc5f68ea16cb83b29ad782debf0 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -278,7 +278,7 @@ struct engine {
const struct hydro_props *hydro_properties;
/* Properties of the self-gravity scheme */
- const struct gravity_props *gravity_properties;
+ struct gravity_props *gravity_properties;
/* Properties of external gravitational potential */
const struct external_potential *external_potential;
@@ -335,7 +335,7 @@ void engine_init(
long long Ngas, long long Ndm, int policy, int verbose,
struct repartition *reparttype, const struct unit_system *internal_units,
const struct phys_const *physical_constants, struct cosmology *cosmo,
- const struct hydro_props *hydro, const struct gravity_props *gravity,
+ const struct hydro_props *hydro, struct gravity_props *gravity,
const struct external_potential *potential,
const struct cooling_function_data *cooling_func,
const struct chemistry_data *chemistry, struct sourceterms *sourceterms);
diff --git a/src/equation_of_state.c b/src/equation_of_state.c
index a792bc58cf13fa5f710dc1cd80713ff90e700b18..0953ee7114314f3baf6e86a14a6ddb0e98442dd9 100644
--- a/src/equation_of_state.c
+++ b/src/equation_of_state.c
@@ -20,48 +20,6 @@
/* This object's header. */
#include "equation_of_state.h"
-/* local headers */
-#include "common_io.h"
-
/* Equation of state for the physics model
* (temporary ugly solution as a global variable) */
struct eos_parameters eos;
-
-void eos_init(struct eos_parameters *e, const struct swift_params *params) {
-
-#if defined(EOS_IDEAL_GAS)
-/* nothing to do here */
-#elif defined(EOS_ISOTHERMAL_GAS)
- e->isothermal_internal_energy =
- parser_get_param_float(params, "EoS:isothermal_internal_energy");
-#endif
-}
-
-void eos_print(const struct eos_parameters *e) {
-
-#if defined(EOS_IDEAL_GAS)
- message("Equation of state: Ideal gas.");
-#elif defined(EOS_ISOTHERMAL_GAS)
- message(
- "Equation of state: Isothermal with internal energy "
- "per unit mass set to %f.",
- e->isothermal_internal_energy);
-#endif
-
- message("Adiabatic index gamma: %f.", hydro_gamma);
-}
-
-#if defined(HAVE_HDF5)
-void eos_print_snapshot(hid_t h_grpsph, const struct eos_parameters *e) {
-
- io_write_attribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
-
-#if defined(EOS_IDEAL_GAS)
- io_write_attribute_s(h_grpsph, "Equation of state", "Ideal gas");
-#elif defined(EOS_ISOTHERMAL_GAS)
- io_write_attribute_s(h_grpsph, "Equation of state", "Isothermal gas");
- io_write_attribute_f(h_grpsph, "Thermal energy per unit mass",
- e->isothermal_internal_energy);
-#endif
-}
-#endif
diff --git a/src/equation_of_state.h b/src/equation_of_state.h
index 76b8bb922133c9c4f29453a14068fe3f9044d66f..195b52514f2acc0c40959e09c088a06f0a411869 100644
--- a/src/equation_of_state.h
+++ b/src/equation_of_state.h
@@ -29,315 +29,13 @@
/* Config parameters. */
#include "../config.h"
-/* Some standard headers. */
-#include <math.h>
-
-/* Local headers. */
-#include "adiabatic_index.h"
-#include "debug.h"
-#include "inline.h"
-
-extern struct eos_parameters eos;
-
-/* ------------------------------------------------------------------------- */
+/* Import the right functions */
#if defined(EOS_IDEAL_GAS)
-
-struct eos_parameters {};
-
-/**
- * @brief Returns the internal energy given density and entropy
- *
- * Computes \f$u = \frac{S\rho^{\gamma-1} }{\gamma - 1}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_entropy(float density, float entropy) {
-
- return entropy * pow_gamma_minus_one(density) *
- hydro_one_over_gamma_minus_one;
-}
-
-/**
- * @brief Returns the pressure given density and entropy
- *
- * Computes \f$P = S\rho^\gamma\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
- float density, float entropy) {
-
- return entropy * pow_gamma(density);
-}
-
-/**
- * @brief Returns the entropy given density and pressure.
- *
- * Computes \f$A = \frac{P}{\rho^\gamma}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param pressure The pressure \f$P\f$.
- * @return The entropy \f$A\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
- float density, float pressure) {
-
- return pressure * pow_minus_gamma(density);
-}
-
-/**
- * @brief Returns the sound speed given density and entropy
- *
- * Computes \f$c = \sqrt{\gamma S \rho^{\gamma-1}}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
- float density, float entropy) {
-
- return sqrtf(hydro_gamma * pow_gamma_minus_one(density) * entropy);
-}
-
-/**
- * @brief Returns the entropy given density and internal energy
- *
- * Computes \f$S = \frac{(\gamma - 1)u}{\rho^{\gamma-1}}\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_entropy_from_internal_energy(float density, float u) {
-
- return hydro_gamma_minus_one * u * pow_minus_gamma_minus_one(density);
-}
-
-/**
- * @brief Returns the pressure given density and internal energy
- *
- * Computes \f$P = (\gamma - 1)u\rho\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_pressure_from_internal_energy(float density, float u) {
-
- return hydro_gamma_minus_one * u * density;
-}
-
-/**
- * @brief Returns the internal energy given density and pressure.
- *
- * Computes \f$u = \frac{1}{\gamma - 1}\frac{P}{\rho}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param pressure The pressure \f$P\f$.
- * @return The internal energy \f$u\f$.
- */
-__attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_pressure(float density, float pressure) {
- return hydro_one_over_gamma_minus_one * pressure / density;
-}
-
-/**
- * @brief Returns the sound speed given density and internal energy
- *
- * Computes \f$c = \sqrt{\gamma (\gamma - 1) u }\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_soundspeed_from_internal_energy(float density, float u) {
-
- return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
-}
-
-/**
- * @brief Returns the sound speed given density and pressure
- *
- * Computes \f$c = \sqrt{\frac{\gamma P}{\rho} }\f$.
- *
- * @param density The density \f$\rho\f$
- * @param P The pressure \f$P\f$
- */
-__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
- float density, float P) {
-
- const float density_inv = 1.f / density;
- return sqrtf(hydro_gamma * P * density_inv);
-}
-
-/* ------------------------------------------------------------------------- */
+#include "./equation_of_state/ideal_gas/equation_of_state.h"
#elif defined(EOS_ISOTHERMAL_GAS)
-
-struct eos_parameters {
-
- /*! Thermal energy per unit mass */
- float isothermal_internal_energy;
-};
-
-/**
- * @brief Returns the internal energy given density and entropy
- *
- * Since we are using an isothermal EoS, the entropy and density values are
- * ignored.
- * Computes \f$u = u_{cst}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_entropy(float density, float entropy) {
-
- return eos.isothermal_internal_energy;
-}
-
-/**
- * @brief Returns the pressure given density and entropy
- *
- * Since we are using an isothermal EoS, the entropy value is ignored.
- * Computes \f$P = (\gamma - 1)u_{cst}\rho\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
- float density, float entropy) {
-
- return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
-}
-
-/**
- * @brief Returns the entropy given density and pressure.
- *
- * Since we are using an isothermal EoS, the pressure value is ignored.
- * Computes \f$A = (\gamma - 1)u_{cst}\rho^{-(\gamma-1)}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param pressure The pressure \f$P\f$ (ignored).
- * @return The entropy \f$A\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
- float density, float pressure) {
-
- return hydro_gamma_minus_one * eos.isothermal_internal_energy *
- pow_minus_gamma_minus_one(density);
-}
-
-/**
- * @brief Returns the sound speed given density and entropy
- *
- * Since we are using an isothermal EoS, the entropy and density values are
- * ignored.
- * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
- *
- * @param density The density \f$\rho\f$.
- * @param entropy The entropy \f$S\f$.
- */
-__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
- float density, float entropy) {
-
- return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
- hydro_gamma_minus_one);
-}
-
-/**
- * @brief Returns the entropy given density and internal energy
- *
- * Since we are using an isothermal EoS, the energy value is ignored.
- * Computes \f$S = \frac{(\gamma - 1)u_{cst}}{\rho^{\gamma-1}}\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_entropy_from_internal_energy(float density, float u) {
-
- return hydro_gamma_minus_one * eos.isothermal_internal_energy *
- pow_minus_gamma_minus_one(density);
-}
-
-/**
- * @brief Returns the pressure given density and internal energy
- *
- * Since we are using an isothermal EoS, the energy value is ignored.
- * Computes \f$P = (\gamma - 1)u_{cst}\rho\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_pressure_from_internal_energy(float density, float u) {
-
- return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
-}
-
-/**
- * @brief Returns the internal energy given density and pressure.
- *
- * Just returns the constant internal energy.
- *
- * @param density The density \f$\rho\f$ (ignored).
- * @param pressure The pressure \f$P\f$ (ignored).
- * @return The internal energy \f$u\f$ (which is constant).
- */
-__attribute__((always_inline)) INLINE static float
-gas_internal_energy_from_pressure(float density, float pressure) {
- return eos.isothermal_internal_energy;
-}
-
-/**
- * @brief Returns the sound speed given density and internal energy
- *
- * Since we are using an isothermal EoS, the energy and density values are
- * ignored.
- * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
- *
- * @param density The density \f$\rho\f$
- * @param u The internal energy \f$u\f$
- */
-__attribute__((always_inline)) INLINE static float
-gas_soundspeed_from_internal_energy(float density, float u) {
-
- return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
- hydro_gamma_minus_one);
-}
-
-/**
- * @brief Returns the sound speed given density and pressure
- *
- * Since we are using an isothermal EoS, the pressure and density values are
- * ignored.
- * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
- *
- * @param density The density \f$\rho\f$
- * @param P The pressure \f$P\f$
- */
-__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
- float density, float P) {
-
- return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
- hydro_gamma_minus_one);
-}
-
-/* ------------------------------------------------------------------------- */
+#include "./equation_of_state/isothermal/equation_of_state.h"
#else
-
-#error "An Equation of state needs to be chosen in const.h !"
-
-#endif
-
-void eos_init(struct eos_parameters *e, const struct swift_params *params);
-void eos_print(const struct eos_parameters *e);
-
-#if defined(HAVE_HDF5)
-void eos_print_snapshot(hid_t h_grpsph, const struct eos_parameters *e);
+#error "Invalid choice of equation of state"
#endif
#endif /* SWIFT_EQUATION_OF_STATE_H */
diff --git a/src/equation_of_state/ideal_gas/equation_of_state.h b/src/equation_of_state/ideal_gas/equation_of_state.h
new file mode 100644
index 0000000000000000000000000000000000000000..42314e0fd87cb5ee2b81bc8cf29029ab4951c869
--- /dev/null
+++ b/src/equation_of_state/ideal_gas/equation_of_state.h
@@ -0,0 +1,203 @@
+/*******************************************************************************
+ * 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_IDEAL_GAS_EQUATION_OF_STATE_H
+#define SWIFT_IDEAL_GAS_EQUATION_OF_STATE_H
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "common_io.h"
+#include "debug.h"
+#include "inline.h"
+
+extern struct eos_parameters eos;
+/* ------------------------------------------------------------------------- */
+
+struct eos_parameters {};
+
+/**
+ * @brief Returns the internal energy given density and entropy
+ *
+ * Computes \f$u = \frac{S\rho^{\gamma-1} }{\gamma - 1}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_entropy(float density, float entropy) {
+
+ return entropy * pow_gamma_minus_one(density) *
+ hydro_one_over_gamma_minus_one;
+}
+
+/**
+ * @brief Returns the pressure given density and entropy
+ *
+ * Computes \f$P = S\rho^\gamma\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
+ float density, float entropy) {
+
+ return entropy * pow_gamma(density);
+}
+
+/**
+ * @brief Returns the entropy given density and pressure.
+ *
+ * Computes \f$A = \frac{P}{\rho^\gamma}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param pressure The pressure \f$P\f$.
+ * @return The entropy \f$A\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
+ float density, float pressure) {
+
+ return pressure * pow_minus_gamma(density);
+}
+
+/**
+ * @brief Returns the sound speed given density and entropy
+ *
+ * Computes \f$c = \sqrt{\gamma S \rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
+ float density, float entropy) {
+
+ return sqrtf(hydro_gamma * pow_gamma_minus_one(density) * entropy);
+}
+
+/**
+ * @brief Returns the entropy given density and internal energy
+ *
+ * Computes \f$S = \frac{(\gamma - 1)u}{\rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_entropy_from_internal_energy(float density, float u) {
+
+ return hydro_gamma_minus_one * u * pow_minus_gamma_minus_one(density);
+}
+
+/**
+ * @brief Returns the pressure given density and internal energy
+ *
+ * Computes \f$P = (\gamma - 1)u\rho\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_pressure_from_internal_energy(float density, float u) {
+
+ return hydro_gamma_minus_one * u * density;
+}
+
+/**
+ * @brief Returns the internal energy given density and pressure.
+ *
+ * Computes \f$u = \frac{1}{\gamma - 1}\frac{P}{\rho}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param pressure The pressure \f$P\f$.
+ * @return The internal energy \f$u\f$.
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_pressure(float density, float pressure) {
+ return hydro_one_over_gamma_minus_one * pressure / density;
+}
+
+/**
+ * @brief Returns the sound speed given density and internal energy
+ *
+ * Computes \f$c = \sqrt{\gamma (\gamma - 1) u }\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_soundspeed_from_internal_energy(float density, float u) {
+
+ return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
+}
+
+/**
+ * @brief Returns the sound speed given density and pressure
+ *
+ * Computes \f$c = \sqrt{\frac{\gamma P}{\rho} }\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param P The pressure \f$P\f$
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
+ float density, float P) {
+
+ const float density_inv = 1.f / density;
+ return sqrtf(hydro_gamma * P * density_inv);
+}
+
+/**
+ * @brief Initialize the eos parameters
+ *
+ * @param e The #eos_paramters
+ * @param params The parsed parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_init(
+ struct eos_parameters *e, const struct swift_params *params) {}
+
+/**
+ * @brief Print the equation of state
+ *
+ * @param e The #eos_parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_print(
+ const struct eos_parameters *e) {
+
+ message("Equation of state: Ideal gas.");
+
+ message("Adiabatic index gamma: %f.", hydro_gamma);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Write equation of state information to the snapshot
+ *
+ * @param h_grpsph The HDF5 group in which to write
+ * @param e The #eos_parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_print_snapshot(
+ hid_t h_grpsph, const struct eos_parameters *e) {
+
+ io_write_attribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
+
+ io_write_attribute_s(h_grpsph, "Equation of state", "Ideal gas");
+}
+#endif
+
+#endif /* SWIFT_IDEAL_GAS_EQUATION_OF_STATE_H */
diff --git a/src/equation_of_state/isothermal/equation_of_state.h b/src/equation_of_state/isothermal/equation_of_state.h
new file mode 100644
index 0000000000000000000000000000000000000000..71890b4df656cb5f44d3cb0fbb3bd6005bd6ab6a
--- /dev/null
+++ b/src/equation_of_state/isothermal/equation_of_state.h
@@ -0,0 +1,233 @@
+/*******************************************************************************
+ * 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_ISOTHERMAL_EQUATION_OF_STATE_H
+#define SWIFT_ISOTHERMAL_EQUATION_OF_STATE_H
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "common_io.h"
+#include "debug.h"
+#include "inline.h"
+
+extern struct eos_parameters eos;
+/* ------------------------------------------------------------------------- */
+
+struct eos_parameters {
+
+ /*! Thermal energy per unit mass */
+ float isothermal_internal_energy;
+};
+
+/**
+ * @brief Returns the internal energy given density and entropy
+ *
+ * Since we are using an isothermal EoS, the entropy and density values are
+ * ignored.
+ * Computes \f$u = u_{cst}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_entropy(float density, float entropy) {
+
+ return eos.isothermal_internal_energy;
+}
+
+/**
+ * @brief Returns the pressure given density and entropy
+ *
+ * Since we are using an isothermal EoS, the entropy value is ignored.
+ * Computes \f$P = (\gamma - 1)u_{cst}\rho\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
+ float density, float entropy) {
+
+ return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
+}
+
+/**
+ * @brief Returns the entropy given density and pressure.
+ *
+ * Since we are using an isothermal EoS, the pressure value is ignored.
+ * Computes \f$A = (\gamma - 1)u_{cst}\rho^{-(\gamma-1)}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param pressure The pressure \f$P\f$ (ignored).
+ * @return The entropy \f$A\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_entropy_from_pressure(
+ float density, float pressure) {
+
+ return hydro_gamma_minus_one * eos.isothermal_internal_energy *
+ pow_minus_gamma_minus_one(density);
+}
+
+/**
+ * @brief Returns the sound speed given density and entropy
+ *
+ * Since we are using an isothermal EoS, the entropy and density values are
+ * ignored.
+ * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
+ *
+ * @param density The density \f$\rho\f$.
+ * @param entropy The entropy \f$S\f$.
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_entropy(
+ float density, float entropy) {
+
+ return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
+ hydro_gamma_minus_one);
+}
+
+/**
+ * @brief Returns the entropy given density and internal energy
+ *
+ * Since we are using an isothermal EoS, the energy value is ignored.
+ * Computes \f$S = \frac{(\gamma - 1)u_{cst}}{\rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_entropy_from_internal_energy(float density, float u) {
+
+ return hydro_gamma_minus_one * eos.isothermal_internal_energy *
+ pow_minus_gamma_minus_one(density);
+}
+
+/**
+ * @brief Returns the pressure given density and internal energy
+ *
+ * Since we are using an isothermal EoS, the energy value is ignored.
+ * Computes \f$P = (\gamma - 1)u_{cst}\rho\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_pressure_from_internal_energy(float density, float u) {
+
+ return hydro_gamma_minus_one * eos.isothermal_internal_energy * density;
+}
+
+/**
+ * @brief Returns the internal energy given density and pressure.
+ *
+ * Just returns the constant internal energy.
+ *
+ * @param density The density \f$\rho\f$ (ignored).
+ * @param pressure The pressure \f$P\f$ (ignored).
+ * @return The internal energy \f$u\f$ (which is constant).
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_pressure(float density, float pressure) {
+ return eos.isothermal_internal_energy;
+}
+
+/**
+ * @brief Returns the sound speed given density and internal energy
+ *
+ * Since we are using an isothermal EoS, the energy and density values are
+ * ignored.
+ * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param u The internal energy \f$u\f$
+ */
+__attribute__((always_inline)) INLINE static float
+gas_soundspeed_from_internal_energy(float density, float u) {
+
+ return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
+ hydro_gamma_minus_one);
+}
+
+/**
+ * @brief Returns the sound speed given density and pressure
+ *
+ * Since we are using an isothermal EoS, the pressure and density values are
+ * ignored.
+ * Computes \f$c = \sqrt{u_{cst} \gamma (\gamma-1)}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param P The pressure \f$P\f$
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
+ float density, float P) {
+
+ return sqrtf(eos.isothermal_internal_energy * hydro_gamma *
+ hydro_gamma_minus_one);
+}
+
+/* ------------------------------------------------------------------------- */
+
+/**
+ * @brief Initialize the eos parameters
+ *
+ * @param e The #eos_paramters
+ * @param params The parsed parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_init(
+ struct eos_parameters *e, const struct swift_params *params) {
+
+ e->isothermal_internal_energy =
+ parser_get_param_float(params, "EoS:isothermal_internal_energy");
+}
+
+/**
+ * @brief Print the equation of state
+ *
+ * @param e The #eos_parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_print(
+ const struct eos_parameters *e) {
+
+ message(
+ "Equation of state: Isothermal with internal energy "
+ "per unit mass set to %f.",
+ e->isothermal_internal_energy);
+
+ message("Adiabatic index gamma: %f.", hydro_gamma);
+}
+
+#if defined(HAVE_HDF5)
+/**
+ * @brief Write equation of state information to the snapshot
+ *
+ * @param h_grpsph The HDF5 group in which to write
+ * @param e The #eos_parameters
+ */
+__attribute__((always_inline)) INLINE static void eos_print_snapshot(
+ hid_t h_grpsph, const struct eos_parameters *e) {
+
+ io_write_attribute_f(h_grpsph, "Adiabatic index", hydro_gamma);
+
+ io_write_attribute_s(h_grpsph, "Equation of state", "Isothermal gas");
+ io_write_attribute_f(h_grpsph, "Thermal energy per unit mass",
+ e->isothermal_internal_energy);
+}
+#endif
+
+#endif /* SWIFT_ISOTHERMAL_EQUATION_OF_STATE_H */
diff --git a/src/gravity.c b/src/gravity.c
index 0d951b4a9f8babb68d413e1183b1e2b52129abaf..1642d37f7a6e543c6e2f263253aebc213ad4b19f 100644
--- a/src/gravity.c
+++ b/src/gravity.c
@@ -372,7 +372,9 @@ int gravity_exact_force_file_exits(const struct engine *e) {
/* Check whether it matches the current parameters */
if (N == SWIFT_GRAVITY_FORCE_CHECKS && periodic == e->s->periodic &&
- (fabs(epsilon - e->gravity_properties->epsilon) / epsilon < 1e-5) &&
+ (fabs(epsilon - gravity_get_softening(0, e->gravity_properties)) /
+ epsilon <
+ 1e-5) &&
(fabs(newton_G - e->physical_constants->const_newton_G) / newton_G <
1e-5)) {
return 1;
@@ -396,6 +398,8 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
struct gpart *restrict gparts = (struct gpart *)map_data;
struct exact_force_data *data = (struct exact_force_data *)extra_data;
const struct space *s = data->s;
+ const struct part *parts = s->parts;
+ const struct spart *sparts = s->sparts;
const struct engine *e = data->e;
const int periodic = s->periodic;
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
@@ -406,13 +410,22 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
struct gpart *gpi = &gparts[i];
+ long long id = 0;
+ if (gpi->type == swift_type_gas)
+ id = parts[-gpi->id_or_neg_offset].id;
+ else if (gpi->type == swift_type_star)
+ id = sparts[-gpi->id_or_neg_offset].id;
+ else if (gpi->type == swift_type_black_hole)
+ error("Unexisting type");
+ else
+ id = gpi->id_or_neg_offset;
+
/* Is the particle active and part of the subset to be tested ? */
- if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
- gpart_is_active(gpi, e)) {
+ if (id % SWIFT_GRAVITY_FORCE_CHECKS == 0 && gpart_is_active(gpi, e)) {
/* Get some information about the particle */
const double pix[3] = {gpi->x[0], gpi->x[1], gpi->x[2]};
- const double hi = gpi->epsilon;
+ const double hi = gravity_get_softening(gpi, e->gravity_properties);
const double hi_inv = 1. / hi;
const double hi_inv3 = hi_inv * hi_inv * hi_inv;
@@ -556,17 +569,21 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ const struct part *parts = s->parts;
+ const struct spart *sparts = s->sparts;
+
/* File name */
char file_name_swift[100];
sprintf(file_name_swift, "gravity_checks_swift_step%d_order%d.dat", e->step,
SELF_GRAVITY_MULTIPOLE_ORDER);
/* Creare files and write header */
+ const double epsilon = gravity_get_softening(0, e->gravity_properties);
FILE *file_swift = fopen(file_name_swift, "w");
fprintf(file_swift, "# Gravity accuracy test - SWIFT FORCES\n");
fprintf(file_swift, "# G= %16.8e\n", e->physical_constants->const_newton_G);
fprintf(file_swift, "# N= %d\n", SWIFT_GRAVITY_FORCE_CHECKS);
- fprintf(file_swift, "# epsilon= %16.8e\n", e->gravity_properties->epsilon);
+ fprintf(file_swift, "# epsilon= %16.8e\n", epsilon);
fprintf(file_swift, "# periodic= %d\n", s->periodic);
fprintf(file_swift, "# theta= %16.8e\n", e->gravity_properties->theta_crit);
fprintf(file_swift, "# Git Branch: %s\n", git_branch());
@@ -580,14 +597,23 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
struct gpart *gpi = &s->gparts[i];
+ long long id = 0;
+ if (gpi->type == swift_type_gas)
+ id = parts[-gpi->id_or_neg_offset].id;
+ else if (gpi->type == swift_type_star)
+ id = sparts[-gpi->id_or_neg_offset].id;
+ else if (gpi->type == swift_type_black_hole)
+ error("Unexisting type");
+ else
+ id = gpi->id_or_neg_offset;
+
/* Is the particle was active and part of the subset to be tested ? */
- if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
- gpart_is_starting(gpi, e)) {
+ if (id % SWIFT_GRAVITY_FORCE_CHECKS == 0 && gpart_is_starting(gpi, e)) {
fprintf(file_swift,
- "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n",
- gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
- gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2], gpi->potential);
+ "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n", id,
+ gpi->x[0], gpi->x[1], gpi->x[2], gpi->a_grav[0], gpi->a_grav[1],
+ gpi->a_grav[2], gpi->potential);
}
}
@@ -609,7 +635,7 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
fprintf(file_exact, "# Gravity accuracy test - EXACT FORCES\n");
fprintf(file_exact, "# G= %16.8e\n", e->physical_constants->const_newton_G);
fprintf(file_exact, "# N= %d\n", SWIFT_GRAVITY_FORCE_CHECKS);
- fprintf(file_exact, "# epsilon=%16.8e\n", e->gravity_properties->epsilon);
+ fprintf(file_exact, "# epsilon=%16.8e\n", epsilon);
fprintf(file_exact, "# periodic= %d\n", s->periodic);
fprintf(file_exact, "# Git Branch: %s\n", git_branch());
fprintf(file_exact, "# Git Revision: %s\n", git_revision());
@@ -622,15 +648,24 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
struct gpart *gpi = &s->gparts[i];
+ long long id = 0;
+ if (gpi->type == swift_type_gas)
+ id = parts[-gpi->id_or_neg_offset].id;
+ else if (gpi->type == swift_type_star)
+ id = sparts[-gpi->id_or_neg_offset].id;
+ else if (gpi->type == swift_type_black_hole)
+ error("Unexisting type");
+ else
+ id = gpi->id_or_neg_offset;
+
/* Is the particle was active and part of the subset to be tested ? */
- if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
- gpart_is_starting(gpi, e)) {
+ if (id % SWIFT_GRAVITY_FORCE_CHECKS == 0 && gpart_is_starting(gpi, e)) {
fprintf(file_exact,
- "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n",
- gpi->id_or_neg_offset, gpi->x[0], gpi->x[1], gpi->x[2],
- gpi->a_grav_exact[0], gpi->a_grav_exact[1],
- gpi->a_grav_exact[2], gpi->potential_exact);
+ "%18lld %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n", id,
+ gpi->x[0], gpi->x[1], gpi->x[2], gpi->a_grav_exact[0],
+ gpi->a_grav_exact[1], gpi->a_grav_exact[2],
+ gpi->potential_exact);
}
}
diff --git a/src/gravity.h b/src/gravity.h
index 85e42370bc456dceb577c42ee609e3f0724e14ea..33bef8ba329aa1264334e3319b6250c01b974338 100644
--- a/src/gravity.h
+++ b/src/gravity.h
@@ -24,16 +24,17 @@
/* Local headers. */
#include "const.h"
-#include "engine.h"
#include "inline.h"
#include "part.h"
-#include "space.h"
/* So far only one model here */
/* Straight-forward import */
#include "./gravity/Default/gravity.h"
#include "./gravity/Default/gravity_iact.h"
+struct engine;
+struct space;
+
void gravity_exact_force_ewald_init(double boxSize);
void gravity_exact_force_ewald_free();
void gravity_exact_force_compute(struct space *s, const struct engine *e);
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 099cbcd158436b414b9bc6a912727116b763cff9..eabc95fc7eb1e1d4039dc9cea6611f1d17451a37 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -21,8 +21,10 @@
#define SWIFT_DEFAULT_GRAVITY_H
#include <float.h>
+
#include "cosmology.h"
#include "gravity_properties.h"
+#include "kernel_gravity.h"
#include "minmax.h"
/**
@@ -40,11 +42,12 @@ __attribute__((always_inline)) INLINE static float gravity_get_mass(
* @brief Returns the softening of a particle
*
* @param gp The particle of interest
+ * @param grav_props The global gravity properties.
*/
__attribute__((always_inline)) INLINE static float gravity_get_softening(
- const struct gpart* restrict gp) {
+ const struct gpart* gp, const struct gravity_props* restrict grav_props) {
- return gp->epsilon;
+ return grav_props->epsilon_cur;
}
/**
@@ -102,12 +105,10 @@ gravity_compute_timestep_self(const struct gpart* const gp,
const float ac_inv = (ac2 > 0.f) ? 1.f / sqrtf(ac2) : FLT_MAX;
- // MATTHIEU cosmological evolution of the softening?
- const float epsilon = gravity_get_softening(gp);
+ const float epsilon = gravity_get_softening(gp, grav_props);
- /* Note that 0.66666667 = 2. (from Gadget) / 3. (Plummer softening) */
- const float dt =
- sqrtf(0.66666667f * cosmo->a * grav_props->eta * epsilon * ac_inv);
+ const float dt = sqrtf(2. * kernel_gravity_softening_plummer_equivalent_inv *
+ cosmo->a * grav_props->eta * epsilon * ac_inv);
return dt;
}
@@ -183,7 +184,6 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
struct gpart* gp, const struct gravity_props* grav_props) {
gp->time_bin = 0;
- gp->epsilon = grav_props->epsilon;
gravity_init_gpart(gp);
}
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
index b83f8c73b4678145f2ecd6d94b136b5aadffccbd..dce038c58e1769446861bdf6c9a2a44415642c68 100644
--- a/src/gravity/Default/gravity_debug.h
+++ b/src/gravity/Default/gravity_debug.h
@@ -22,9 +22,9 @@
__attribute__((always_inline)) INLINE static void gravity_debug_particle(
const struct gpart* p) {
printf(
- "mass=%.3e epsilon=%.5e time_bin=%d\n"
+ "mass=%.3e time_bin=%d\n"
"x=[%.5e,%.5e,%.5e], v_full=[%.5e,%.5e,%.5e], a=[%.5e,%.5e,%.5e]\n",
- p->mass, p->epsilon, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
+ p->mass, p->time_bin, p->x[0], p->x[1], p->x[2], p->v_full[0],
p->v_full[1], p->v_full[2], p->a_grav[0], p->a_grav[1], p->a_grav[2]);
#ifdef SWIFT_DEBUG_CHECKS
printf("num_interacted=%lld ti_drift=%lld ti_kick=%lld\n", p->num_interacted,
diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h
index 26f7a6d15e11cf2e2d9ccd9978a1c67bab17d852..7b8ec2c8fef04cc4a4fc6836d8bb895b24d3c41f 100644
--- a/src/gravity/Default/gravity_io.h
+++ b/src/gravity/Default/gravity_io.h
@@ -35,6 +35,38 @@ void convert_gpart_pos(const struct engine* e, const struct gpart* gp,
}
}
+void convert_gpart_vel(const struct engine* e, const struct gpart* gp,
+ float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, gp->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ ret[0] = gp->v_full[0] + gp->a_grav[0] * dt_kick_grav;
+ ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav;
+ ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav;
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a2_inv;
+ ret[1] *= cosmo->a2_inv;
+ ret[2] *= cosmo->a2_inv;
+}
+
/**
* @brief Specifies which g-particle fields to read from a dataset
*
@@ -55,7 +87,7 @@ void darkmatter_read_particles(struct gpart* gparts, struct io_props* list,
UNIT_CONV_SPEED, gparts, v_full);
list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
gparts, mass);
- list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY,
+ list[3] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset);
}
@@ -69,20 +101,20 @@ void darkmatter_read_particles(struct gpart* gparts, struct io_props* list,
void darkmatter_write_particles(const struct gpart* gparts,
struct io_props* list, int* num_fields) {
- /* Say how much we want to read */
+ /* Say how much we want to write */
*num_fields = 5;
- /* List what we want to read */
+ /* List what we want to write */
list[0] = io_make_output_field_convert_gpart(
"Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, gparts, convert_gpart_pos);
- list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED,
- gparts, v_full);
+ list[1] = io_make_output_field_convert_gpart(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, gparts, convert_gpart_vel);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, gparts, mass);
list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset);
- list[4] = io_make_output_field("Acceleration", FLOAT, 3,
- UNIT_CONV_ACCELERATION, gparts, a_grav);
+ list[4] = io_make_output_field("Potential", FLOAT, 1, UNIT_CONV_POTENTIAL,
+ gparts, potential);
}
#endif /* SWIFT_DEFAULT_GRAVITY_IO_H */
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 4c7a7b5d4af1c26ee96ac797963f5588d3303d77..06b9c52e14a655aa97b6c36ff9dca1f1cdd6e9a0 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -22,35 +22,32 @@
/* Gravity particle. */
struct gpart {
- /* Particle ID. If negative, it is the negative offset of the #part with
+ /*! Particle ID. If negative, it is the negative offset of the #part with
which this gpart is linked. */
long long id_or_neg_offset;
- /* Particle position. */
+ /*! Particle position. */
double x[3];
- /* Offset between current position and position at last tree rebuild. */
+ /*! Offset between current position and position at last tree rebuild. */
float x_diff[3];
- /* Particle velocity. */
+ /*! Particle velocity. */
float v_full[3];
- /* Particle acceleration. */
+ /*! Particle acceleration. */
float a_grav[3];
- /* Particle mass. */
+ /*! Particle mass. */
float mass;
- /* Gravitational potential */
+ /*! Gravitational potential */
float potential;
- /* Softening length */
- float epsilon;
-
- /* Time-step length */
+ /*! Time-step length */
timebin_t time_bin;
- /* Type of the #gpart (DM, gas, star, ...) */
+ /*! Type of the #gpart (DM, gas, star, ...) */
enum part_type type;
#ifdef SWIFT_DEBUG_CHECKS
diff --git a/src/gravity_cache.h b/src/gravity_cache.h
index be349cfd4a4583a493219eb020757341c98f540c..0012be193c657abccd5d83ed872e4d21764065a0 100644
--- a/src/gravity_cache.h
+++ b/src/gravity_cache.h
@@ -149,14 +149,16 @@ static INLINE void gravity_cache_init(struct gravity_cache *c, int count) {
* @param shift A shift to apply to all the particles.
* @param CoM The position of the multipole.
* @param r_max2 The square of the multipole radius.
- * @param theta_crit2 The square of the opening angle.
* @param cell The cell we play with (to get reasonable padding positions).
+ * @param grav_props The global gravity properties.
*/
__attribute__((always_inline)) INLINE static void gravity_cache_populate(
timebin_t max_active_bin, struct gravity_cache *c,
const struct gpart *restrict gparts, int gcount, int gcount_padded,
- const double shift[3], const float CoM[3], float r_max2, float theta_crit2,
- const struct cell *cell) {
+ const double shift[3], const float CoM[3], float r_max2,
+ const struct cell *cell, const struct gravity_props *grav_props) {
+
+ const float theta_crit2 = grav_props->theta_crit2;
/* Make the compiler understand we are in happy vectorization land */
swift_declare_aligned_ptr(float, x, c->x, SWIFT_CACHE_ALIGNMENT);
@@ -174,7 +176,7 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
x[i] = (float)(gparts[i].x[0] - shift[0]);
y[i] = (float)(gparts[i].x[1] - shift[1]);
z[i] = (float)(gparts[i].x[2] - shift[2]);
- epsilon[i] = gparts[i].epsilon;
+ epsilon[i] = gravity_get_softening(&gparts[i], grav_props);
m[i] = gparts[i].mass;
active[i] = (int)(gparts[i].time_bin <= max_active_bin);
@@ -220,13 +222,15 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
* multiple of the vector length.
* @param shift A shift to apply to all the particles.
* @param cell The cell we play with (to get reasonable padding positions).
+ * @param grav_props The global gravity properties.
*/
__attribute__((always_inline)) INLINE static void
gravity_cache_populate_no_mpole(timebin_t max_active_bin,
struct gravity_cache *c,
const struct gpart *restrict gparts, int gcount,
int gcount_padded, const double shift[3],
- const struct cell *cell) {
+ const struct cell *cell,
+ const struct gravity_props *grav_props) {
/* Make the compiler understand we are in happy vectorization land */
swift_declare_aligned_ptr(float, x, c->x, SWIFT_CACHE_ALIGNMENT);
@@ -242,7 +246,7 @@ gravity_cache_populate_no_mpole(timebin_t max_active_bin,
x[i] = (float)(gparts[i].x[0] - shift[0]);
y[i] = (float)(gparts[i].x[1] - shift[1]);
z[i] = (float)(gparts[i].x[2] - shift[2]);
- epsilon[i] = gparts[i].epsilon;
+ epsilon[i] = gravity_get_softening(&gparts[i], grav_props);
m[i] = gparts[i].mass;
active[i] = (int)(gparts[i].time_bin <= max_active_bin);
}
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 1df421a8b21424dc99e52bc9dbc59560e54d14a5..79be11a8bab2389c19c58a59851379ea23d5fb3f 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -37,7 +37,8 @@
#define gravity_props_default_r_cut_min 0.1f
void gravity_props_init(struct gravity_props *p,
- const struct swift_params *params) {
+ const struct swift_params *params,
+ const struct cosmology *cosmo) {
/* Tree-PM parameters */
p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
@@ -56,11 +57,35 @@ void gravity_props_init(struct gravity_props *p,
p->theta_crit2 = p->theta_crit * p->theta_crit;
p->theta_crit_inv = 1. / p->theta_crit;
- /* Softening lengths */
- p->epsilon = 3. * parser_get_param_double(params, "Gravity:epsilon");
- p->epsilon2 = p->epsilon * p->epsilon;
- p->epsilon_inv = 1.f / p->epsilon;
- p->epsilon_inv3 = p->epsilon_inv * p->epsilon_inv * p->epsilon_inv;
+ /* Softening parameters */
+ p->epsilon_comoving =
+ parser_get_param_double(params, "Gravity:comoving_softening");
+ p->epsilon_max_physical =
+ parser_get_param_double(params, "Gravity:max_physical_softening");
+
+ /* Set the softening to the current time */
+ gravity_update(p, cosmo);
+}
+
+void gravity_update(struct gravity_props *p, const struct cosmology *cosmo) {
+
+ /* Current softening lengths */
+ double softening;
+ if (p->epsilon_comoving * cosmo->a > p->epsilon_max_physical)
+ softening = p->epsilon_max_physical / cosmo->a;
+ else
+ softening = p->epsilon_comoving;
+
+ /* Plummer equivalent -> internal */
+ softening *= kernel_gravity_softening_plummer_equivalent;
+
+ /* Store things */
+ p->epsilon_cur = softening;
+
+ /* Other factors */
+ p->epsilon_cur2 = softening * softening;
+ p->epsilon_cur_inv = 1. / softening;
+ p->epsilon_cur_inv3 = 1. / (softening * softening * softening);
}
void gravity_props_print(const struct gravity_props *p) {
@@ -72,8 +97,17 @@ void gravity_props_print(const struct gravity_props *p) {
message("Self-gravity opening angle: theta=%.4f", p->theta_crit);
- message("Self-gravity softening: epsilon=%.4f (Plummer equivalent: %.4f)",
- p->epsilon, p->epsilon / 3.);
+ message(
+ "Self-gravity comoving softening: epsilon=%.4f (Plummer equivalent: "
+ "%.4f)",
+ p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent,
+ p->epsilon_comoving);
+
+ message(
+ "Self-gravity maximal physical softening: epsilon=%.4f (Plummer "
+ "equivalent: %.4f)",
+ p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent,
+ p->epsilon_max_physical);
message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth);
@@ -86,9 +120,18 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
const struct gravity_props *p) {
io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
- io_write_attribute_f(h_grpgrav, "Softening length", p->epsilon);
- io_write_attribute_f(h_grpgrav, "Softening length (Plummer equivalent)",
- p->epsilon / 3.);
+ io_write_attribute_f(
+ h_grpgrav, "Comoving softening length",
+ p->epsilon_comoving * kernel_gravity_softening_plummer_equivalent);
+ io_write_attribute_f(h_grpgrav,
+ "Comoving Softening length (Plummer equivalent)",
+ p->epsilon_comoving);
+ io_write_attribute_f(
+ h_grpgrav, "Maximal physical softening length",
+ p->epsilon_max_physical * kernel_gravity_softening_plummer_equivalent);
+ io_write_attribute_f(h_grpgrav,
+ "Maximal physical softening length (Plummer equivalent)",
+ p->epsilon_max_physical);
io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
diff --git a/src/gravity_properties.h b/src/gravity_properties.h
index 826ffd0de05d376b930523c8a5d937e457c6d795..f36a39a6e187b8b397a5f597692f5a37c982aa7a 100644
--- a/src/gravity_properties.h
+++ b/src/gravity_properties.h
@@ -27,6 +27,7 @@
#endif
/* Local includes. */
+#include "cosmology.h"
#include "parser.h"
#include "restart.h"
@@ -58,22 +59,30 @@ struct gravity_props {
/*! Inverse of opening angle */
double theta_crit_inv;
- /*! Softening length */
- float epsilon;
+ /*! Comoving softening */
+ double epsilon_comoving;
- /*! Square of softening length */
- float epsilon2;
+ /*! Maxium physical softening */
+ double epsilon_max_physical;
- /*! Inverse of softening length */
- float epsilon_inv;
+ /*! Current sftening length */
+ float epsilon_cur;
- /*! Cube of the inverse of softening length */
- float epsilon_inv3;
+ /*! Square of current softening length */
+ float epsilon_cur2;
+
+ /*! Inverse of current softening length */
+ float epsilon_cur_inv;
+
+ /*! Cube of the inverse of current oftening length */
+ float epsilon_cur_inv3;
};
void gravity_props_print(const struct gravity_props *p);
void gravity_props_init(struct gravity_props *p,
- const struct swift_params *params);
+ const struct swift_params *params,
+ const struct cosmology *cosmo);
+void gravity_update(struct gravity_props *p, const struct cosmology *cosmo);
#if defined(HAVE_HDF5)
void gravity_props_print_snapshot(hid_t h_grpsph,
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index 1228aa69abfe812e5b2e73f066bb13be3292aa20..7da005456233b800e05e75389922395ea3c1702c 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -160,6 +160,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->mass;
}
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
@@ -492,7 +504,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param dt_therm The time-step for this kick (for thermodynamic quantities)
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
- struct part *restrict p, struct xpart *restrict xp, float dt_therm) {}
+ struct part *restrict p, struct xpart *restrict xp, float dt_therm,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props) {}
/**
* @brief Converts hydro quantity of a particle at the start of a run
@@ -500,9 +513,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* Requires the density to be known
*
* @param p The particle to act upon
+ * @param xp The extended data.
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part *restrict p, struct xpart *restrict xp) {}
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {}
/**
* @brief Initialises the particles for the first time
@@ -529,4 +545,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
hydro_init_part(p, NULL);
}
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+ p->u = u_init;
+}
+
#endif /* SWIFT_DEFAULT_HYDRO_H */
diff --git a/src/hydro/Default/hydro_io.h b/src/hydro/Default/hydro_io.h
index bfff08dcdf061d86986ffdeabee6e7a84cba1c7b..e70a26f9edc3f7497dd4c548efd464726efe0e39 100644
--- a/src/hydro/Default/hydro_io.h
+++ b/src/hydro/Default/hydro_io.h
@@ -54,7 +54,7 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
}
void convert_part_pos(const struct engine* e, const struct part* p,
- double* ret) {
+ const struct xpart* xp, double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
@@ -67,6 +67,49 @@ void convert_part_pos(const struct engine* e, const struct part* p,
}
}
+void convert_part_vel(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a2_inv;
+ ret[1] *= cosmo->a2_inv;
+ ret[2] *= cosmo->a2_inv;
+}
+
+void convert_part_potential(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ if (p->gpart != NULL)
+ ret[0] = gravity_get_comoving_potential(p->gpart);
+ else
+ ret[0] = 0.f;
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -74,16 +117,17 @@ void convert_part_pos(const struct engine* e, const struct part* p,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
-void hydro_write_particles(struct part* parts, struct io_props* list,
- int* num_fields) {
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
- *num_fields = 7;
+ *num_fields = 8;
/* List what we want to write */
- list[0] = io_make_output_field_convert_part(
- "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
- list[1] =
- io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
@@ -94,6 +138,9 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
UNIT_CONV_NO_UNITS, parts, id);
list[6] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
+ list[7] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+ UNIT_CONV_POTENTIAL, parts,
+ xparts, convert_part_potential);
}
/**
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 185a98b6d207bb851512a686be51414341fe7740..bc06a24e2a8245556a1042f2459273b8d750489e 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -170,6 +170,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->mass;
}
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
@@ -513,19 +525,29 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt_therm The time-step for this kick (for thermodynamic quantities)
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
- struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
+ struct part *restrict p, struct xpart *restrict xp, float dt_therm,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Do not decrease the entropy by more than a factor of 2 */
if (dt_therm > 0. && p->entropy_dt * dt_therm < -0.5f * xp->entropy_full) {
- /* message("Warning! Limiting entropy_dt. Possible cooling error.\n
- * entropy_full = %g \n entropy_dt * dt =%g \n", */
- /* xp->entropy_full,p->entropy_dt * dt); */
p->entropy_dt = -0.5f * xp->entropy_full / dt_therm;
}
xp->entropy_full += p->entropy_dt * dt_therm;
+ /* Apply the minimal energy limit */
+ const float density = p->rho * cosmo->a3_inv;
+ const float min_energy = hydro_props->minimal_internal_energy;
+ const float min_entropy =
+ gas_entropy_from_internal_energy(density, min_energy);
+ if (xp->entropy_full < min_entropy) {
+ xp->entropy_full = min_entropy;
+ p->entropy_dt = 0.f;
+ }
+
/* Compute the pressure */
const float pressure = gas_pressure_from_entropy(p->rho, xp->entropy_full);
@@ -545,13 +567,18 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*
* Requires the density to be known
*
- * @param p The particle to act upon
+ * @param p The particle to act upon.
+ * @param xp The extended data.
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part *restrict p, struct xpart *restrict xp) {
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
- /* We read u in the entropy field. We now get S from u */
- xp->entropy_full = gas_entropy_from_internal_energy(p->rho, p->entropy);
+ /* We read u in the entropy field. We now get (comoving) S from (physical) u
+ * and (physical) rho. Note that comoving S == physical S */
+ xp->entropy_full =
+ gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, p->entropy);
p->entropy = xp->entropy_full;
/* Compute the pressure */
@@ -593,4 +620,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
hydro_init_part(p, NULL);
}
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+ p->entropy = u_init;
+}
+
#endif /* SWIFT_GADGET2_HYDRO_H */
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index b7fa7dd4aa00249d009037961888c89abc0d0313..b2af8909bed1780586a5130370222c9b8157d724 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -59,9 +59,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
const float mi = pi->mass;
const float mj = pj->mass;
- /* Get r and r inverse. */
- const float r = sqrtf(r2);
- const float r_inv = 1.0f / r;
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Compute the kernel function for pi */
const float hi_inv = 1.f / hi;
@@ -148,9 +148,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Get the masses. */
const float mj = pj->mass;
- /* Get r and r inverse. */
- const float r = sqrtf(r2);
- const float r_inv = 1.0f / r;
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Compute the kernel function */
const float hi_inv = 1.0f / hi;
@@ -440,8 +440,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
- const float r = sqrtf(r2);
- const float r_inv = 1.0f / r;
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Get some values in local variables. */
const float mi = pi->mass;
@@ -559,8 +560,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
- const float r = sqrtf(r2);
- const float r_inv = 1.0f / r;
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Get some values in local variables. */
// const float mi = pi->mass;
diff --git a/src/hydro/Gadget2/hydro_io.h b/src/hydro/Gadget2/hydro_io.h
index 5c329e9d7297109c378e04f0e4931a0717e07278..28c0eea4772f51fab35a08d43c0564472694eeeb 100644
--- a/src/hydro/Gadget2/hydro_io.h
+++ b/src/hydro/Gadget2/hydro_io.h
@@ -55,18 +55,20 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho);
}
-void convert_u(const struct engine* e, const struct part* p, float* ret) {
+void convert_part_u(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_internal_energy(p);
}
-void convert_P(const struct engine* e, const struct part* p, float* ret) {
+void convert_part_P(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_pressure(p);
}
void convert_part_pos(const struct engine* e, const struct part* p,
- double* ret) {
+ const struct xpart* xp, double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
@@ -79,6 +81,49 @@ void convert_part_pos(const struct engine* e, const struct part* p,
}
}
+void convert_part_vel(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a2_inv;
+ ret[1] *= cosmo->a2_inv;
+ ret[2] *= cosmo->a2_inv;
+}
+
+void convert_part_potential(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ if (p->gpart != NULL)
+ ret[0] = gravity_get_comoving_potential(p->gpart);
+ else
+ ret[0] = 0.f;
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -86,16 +131,17 @@ void convert_part_pos(const struct engine* e, const struct part* p,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
-void hydro_write_particles(const struct part* parts, struct io_props* list,
- int* num_fields) {
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
- *num_fields = 9;
+ *num_fields = 10;
/* List what we want to write */
- list[0] = io_make_output_field_convert_part(
- "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
- list[1] =
- io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
@@ -108,9 +154,13 @@ void hydro_write_particles(const struct part* parts, struct io_props* list,
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
list[7] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
- parts, convert_u);
+ parts, xparts, convert_part_u);
list[8] = io_make_output_field_convert_part(
- "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
+ "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_part_P);
+
+ list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+ UNIT_CONV_POTENTIAL, parts,
+ xparts, convert_part_potential);
#ifdef DEBUG_INTERACTIONS_SPH
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 0faa4bf39ae533ee9942898bdc08118ee4078868..bc90a6790fa8123b19badf63761a6dd1378197f6 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -517,7 +517,7 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part* p, struct xpart* xp) {}
+ struct part* p, struct xpart* xp, const struct cosmology* cosmo) {}
/**
* @brief Extra operations to be done during the drift
@@ -611,7 +611,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param half_dt Half the physical time step.
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
- struct part* p, struct xpart* xp, float dt) {
+ struct part* p, struct xpart* xp, float dt, const struct cosmology* cosmo,
+ const struct hydro_props* hydro_props) {
float a_grav[3];
@@ -628,6 +629,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.energy += p->conserved.flux.energy * dt;
#endif
+ /* Apply the minimal energy limit */
+ const float min_energy =
+ hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
+ if (p->conserved.energy < min_energy) {
+ p->conserved.energy = min_energy;
+ p->conserved.flux.energy = 0.f;
+ }
+
gizmo_check_physical_quantity("mass", p->conserved.mass);
gizmo_check_physical_quantity("energy", p->conserved.energy);
@@ -812,6 +821,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->conserved.mass;
}
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part* restrict p, float m) {
+
+ p->conserved.mass = m;
+}
+
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
@@ -917,4 +938,29 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
p->primitives.P = S * pow_gamma(p->primitives.rho);
}
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part* p, float u_init) {
+
+ p->conserved.energy = u_init * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+ /* add the kinetic energy */
+ p->conserved.energy += 0.5f * p->conserved.mass *
+ (p->conserved.momentum[0] * p->primitives.v[0] +
+ p->conserved.momentum[1] * p->primitives.v[1] +
+ p->conserved.momentum[2] * p->primitives.v[2]);
+#endif
+ p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u_init;
+}
+
#endif /* SWIFT_GIZMO_HYDRO_H */
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index b8e2e30646c3aafb73d387010ab7ea0b0a663a77..7aa56b3bde5247f4839b8f4f60cdb3de67253392 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -72,7 +72,8 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
* @param p Particle.
* @param ret (return) Internal energy of the particle
*/
-void convert_u(const struct engine* e, const struct part* p, float* ret) {
+void convert_u(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_internal_energy(p);
}
@@ -84,7 +85,8 @@ void convert_u(const struct engine* e, const struct part* p, float* ret) {
* @param p Particle.
* @param ret (return) Entropic function of the particle
*/
-void convert_A(const struct engine* e, const struct part* p, float* ret) {
+void convert_A(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_entropy(p);
}
@@ -95,7 +97,8 @@ void convert_A(const struct engine* e, const struct part* p, float* ret) {
* @param p Particle.
* @return Total energy of the particle
*/
-void convert_Etot(const struct engine* e, const struct part* p, float* ret) {
+void convert_Etot(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
#ifdef GIZMO_TOTAL_ENERGY
ret[0] = p->conserved.energy;
#else
@@ -110,7 +113,7 @@ void convert_Etot(const struct engine* e, const struct part* p, float* ret) {
}
void convert_part_pos(const struct engine* e, const struct part* p,
- double* ret) {
+ const struct xpart* xp, double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
@@ -123,6 +126,49 @@ void convert_part_pos(const struct engine* e, const struct part* p,
}
}
+void convert_part_vel(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a2_inv;
+ ret[1] *= cosmo->a2_inv;
+ ret[2] *= cosmo->a2_inv;
+}
+
+void convert_part_potential(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ if (p->gpart != NULL)
+ ret[0] = gravity_get_comoving_potential(p->gpart);
+ else
+ ret[0] = 0.f;
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -130,33 +176,39 @@ void convert_part_pos(const struct engine* e, const struct part* p,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
-void hydro_write_particles(struct part* parts, struct io_props* list,
- int* num_fields) {
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
- *num_fields = 10;
+ *num_fields = 11;
/* List what we want to write */
- list[0] = io_make_output_field_convert_part(
- "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
- list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts,
- primitives.v);
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
+
list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
conserved.mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
parts, h);
list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
- parts, convert_u);
+ parts, xparts, convert_u);
list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
UNIT_CONV_NO_UNITS, parts, id);
list[6] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
primitives.rho);
list[7] = io_make_output_field_convert_part(
- "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, convert_A);
+ "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY, parts, xparts, convert_A);
list[8] = io_make_output_field("Pressure", FLOAT, 1, UNIT_CONV_PRESSURE,
parts, primitives.P);
list[9] = io_make_output_field_convert_part(
- "TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, convert_Etot);
+ "TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY, parts, xparts, convert_Etot);
+
+ list[10] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+ UNIT_CONV_POTENTIAL, parts,
+ xparts, convert_part_potential);
}
/**
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 97cd5ef577dcfd7b07c5183a92349313096a6412..3f9d99683bde4ed6db64d8aaa5b111e2f67f0969 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -197,6 +197,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->mass;
}
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
@@ -510,9 +522,12 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p The particle to act upon.
* @param xp The particle extended data to act upon.
* @param dt_therm The time-step for this kick (for thermodynamic quantities).
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
- struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
+ struct part *restrict p, struct xpart *restrict xp, float dt_therm,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Do not decrease the energy by more than a factor of 2*/
if (dt_therm > 0. && p->u_dt * dt_therm < -0.5f * xp->u_full) {
@@ -520,6 +535,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
}
xp->u_full += p->u_dt * dt_therm;
+ /* Apply the minimal energy limit */
+ const float min_energy =
+ hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
+ if (xp->u_full < min_energy) {
+ xp->u_full = min_energy;
+ p->u_dt = 0.f;
+ }
+
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, xp->u_full);
@@ -540,9 +563,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
*
* @param p The particle to act upon
* @param xp The extended particle to act upon
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part *restrict p, struct xpart *restrict xp) {
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
@@ -580,4 +605,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
hydro_init_part(p, NULL);
}
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+ p->u = u_init;
+}
+
#endif /* SWIFT_MINIMAL_HYDRO_H */
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index bfd8f1dff5febbd14b8a1e3a1cb8c01d49bc85ba..42fd93d6062cbfcea5cf5297eeda0bb6525f3cad 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -53,7 +53,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
float wi, wj, wi_dx, wj_dx;
- const float r = sqrtf(r2);
+ /* Get r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Get the masses. */
const float mi = pi->mass;
@@ -101,8 +103,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Get the masses. */
const float mj = pj->mass;
- /* Get r and r inverse. */
- const float r = sqrtf(r2);
+ /* Get r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
const float h_inv = 1.f / hi;
const float ui = r * h_inv;
@@ -134,8 +137,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
- const float r = sqrtf(r2);
- const float r_inv = 1.0f / r;
+ /* Get r and r inverse. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Recover some data */
const float mi = pi->mass;
@@ -246,8 +250,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
- const float r = sqrtf(r2);
- const float r_inv = 1.0f / r;
+ /* Get r and r inverse. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Recover some data */
// const float mi = pi->mass;
diff --git a/src/hydro/Minimal/hydro_io.h b/src/hydro/Minimal/hydro_io.h
index 2ff70fd2f7e95344c17ba0f3237c3b2ee131752c..380d6120e05acf2e015ece6f133df02fad3b761d 100644
--- a/src/hydro/Minimal/hydro_io.h
+++ b/src/hydro/Minimal/hydro_io.h
@@ -69,18 +69,20 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho);
}
-void convert_S(const struct engine* e, const struct part* p, float* ret) {
+void convert_S(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_entropy(p);
}
-void convert_P(const struct engine* e, const struct part* p, float* ret) {
+void convert_P(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_pressure(p);
}
void convert_part_pos(const struct engine* e, const struct part* p,
- double* ret) {
+ const struct xpart* xp, double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
@@ -93,23 +95,68 @@ void convert_part_pos(const struct engine* e, const struct part* p,
}
}
+void convert_part_vel(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a2_inv;
+ ret[1] *= cosmo->a2_inv;
+ ret[2] *= cosmo->a2_inv;
+}
+
+void convert_part_potential(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ if (p->gpart != NULL)
+ ret[0] = gravity_get_comoving_potential(p->gpart);
+ else
+ ret[0] = 0.f;
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
* @param parts The particle array.
+ * @param xparts The extended particle array.
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
-void hydro_write_particles(struct part* parts, struct io_props* list,
- int* num_fields) {
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
- *num_fields = 9;
+ *num_fields = 10;
/* List what we want to write */
- list[0] = io_make_output_field_convert_part(
- "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
- list[1] =
- io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
@@ -120,10 +167,15 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
UNIT_CONV_NO_UNITS, parts, id);
list[6] =
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
- list[7] = io_make_output_field_convert_part(
- "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, convert_S);
+ list[7] = io_make_output_field_convert_part("Entropy", FLOAT, 1,
+ UNIT_CONV_ENTROPY_PER_UNIT_MASS,
+ parts, xparts, convert_S);
list[8] = io_make_output_field_convert_part(
- "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
+ "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_P);
+
+ list[9] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+ UNIT_CONV_POTENTIAL, parts,
+ xparts, convert_part_potential);
}
/**
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index 37de855b375a5ac11f8a33f3993f7d2101c58522..87d46c6d43f0d4f6de6d18f5400b38f0fc4d0f55 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -170,6 +170,18 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
return p->mass;
}
+/**
+ * @brief Sets the mass of a particle
+ *
+ * @param p The particle of interest
+ * @param m The mass to set.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_mass(
+ struct part *restrict p, float m) {
+
+ p->mass = m;
+}
+
/**
* @brief Returns the velocities drifted to the current time of a particle.
*
@@ -524,9 +536,12 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
* @param p The particle to act upon
* @param xp The particle extended data to act upon
* @param dt_therm The time-step for this kick (for thermodynamic quantities)
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme
*/
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
- struct part *restrict p, struct xpart *restrict xp, float dt_therm) {
+ struct part *restrict p, struct xpart *restrict xp, float dt_therm,
+ const struct cosmology *cosmo, const struct hydro_props *hydro_props) {
/* Do not decrease the entropy (temperature) by more than a factor of 2*/
if (dt_therm > 0. && p->entropy_dt * dt_therm < -0.5f * xp->entropy_full) {
@@ -534,6 +549,16 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
}
xp->entropy_full += p->entropy_dt * dt_therm;
+ /* Apply the minimal energy limit */
+ const float density = p->rho_bar * cosmo->a3_inv;
+ const float min_energy = hydro_props->minimal_internal_energy;
+ const float min_entropy =
+ gas_entropy_from_internal_energy(density, min_energy);
+ if (xp->entropy_full < min_entropy) {
+ xp->entropy_full = min_entropy;
+ p->entropy_dt = 0.f;
+ }
+
/* Compute the pressure */
const float pressure =
gas_pressure_from_entropy(p->rho_bar, xp->entropy_full);
@@ -556,12 +581,16 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* Requires the density to be known
*
* @param p The particle to act upon
+ * @param xp The extended data.
+ * @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part *restrict p, struct xpart *restrict xp) {
+ struct part *restrict p, struct xpart *restrict xp,
+ const struct cosmology *cosmo) {
/* We read u in the entropy field. We now get S from u */
- xp->entropy_full = gas_entropy_from_internal_energy(p->rho_bar, p->entropy);
+ xp->entropy_full =
+ gas_entropy_from_internal_energy(p->rho_bar * cosmo->a3_inv, p->entropy);
p->entropy = xp->entropy_full;
p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
@@ -605,4 +634,21 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
hydro_init_part(p, NULL);
}
+/**
+ * @brief Overwrite the initial internal energy of a particle.
+ *
+ * Note that in the cases where the thermodynamic variable is not
+ * internal energy but gets converted later, we must overwrite that
+ * field. The conversion to the actual variable happens later after
+ * the initial fake time-step.
+ *
+ * @param p The #part to write to.
+ * @param u_init The new initial internal energy.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_set_init_internal_energy(struct part *p, float u_init) {
+
+ p->entropy = u_init;
+}
+
#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */
diff --git a/src/hydro/PressureEntropy/hydro_iact.h b/src/hydro/PressureEntropy/hydro_iact.h
index c9e3c4ec7728df8a8e7d59b8c21f22bd613463d8..b8f8c1983a3b1fb67781f7228194deb770273988 100644
--- a/src/hydro/PressureEntropy/hydro_iact.h
+++ b/src/hydro/PressureEntropy/hydro_iact.h
@@ -54,9 +54,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
const float mi = pi->mass;
const float mj = pj->mass;
- /* Get r and r inverse. */
- const float r = sqrtf(r2);
- const float r_inv = 1.0f / r;
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Compute the kernel function for pi */
const float hi_inv = 1.f / hi;
@@ -142,9 +142,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
/* Get the masses. */
const float mj = pj->mass;
- /* Get r and r inverse. */
- const float r = sqrtf(r2);
- const float ri = 1.0f / r;
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Compute the kernel function */
const float h_inv = 1.0f / hi;
@@ -164,7 +164,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
pi->density.pressure_dh -=
mj * pj->entropy_one_over_gamma * (hydro_dimension * wi + ui * wi_dx);
- const float fac = mj * wi_dx * ri;
+ const float fac = mj * wi_dx * r_inv;
/* Compute dv dot r */
dv[0] = pi->v[0] - pj->v[0];
@@ -205,8 +205,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
- const float r = sqrtf(r2);
- const float r_inv = 1.0f / r;
+ /* Get r and 1/r . */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Get some values in local variables. */
const float mi = pi->mass;
@@ -318,8 +319,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float fac_mu = pow_three_gamma_minus_five_over_two(a);
const float a2_Hubble = a * a * H;
- const float r = sqrtf(r2);
- const float r_inv = 1.0f / r;
+ /* Get r and 1/r. */
+ const float r_inv = 1.0f / sqrtf(r2);
+ const float r = r2 * r_inv;
/* Get some values in local variables. */
// const float mi = pi->mass;
diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h
index ed8e0e45d240ea37ebfb7df6339afd454dc0e4c4..78371c1eb21fafed56e89d46690d0cf1e0f2a0f0 100644
--- a/src/hydro/PressureEntropy/hydro_io.h
+++ b/src/hydro/PressureEntropy/hydro_io.h
@@ -67,18 +67,20 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
UNIT_CONV_DENSITY, parts, rho);
}
-void convert_u(const struct engine* e, const struct part* p, float* ret) {
+void convert_u(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_internal_energy(p);
}
-void convert_P(const struct engine* e, const struct part* p, float* ret) {
+void convert_P(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
ret[0] = hydro_get_comoving_pressure(p);
}
void convert_part_pos(const struct engine* e, const struct part* p,
- double* ret) {
+ const struct xpart* xp, double* ret) {
if (e->s->periodic) {
ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
@@ -91,6 +93,49 @@ void convert_part_pos(const struct engine* e, const struct part* p,
}
}
+void convert_part_vel(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav, dt_kick_hydro;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_hydro -=
+ cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a2_inv;
+ ret[1] *= cosmo->a2_inv;
+ ret[2] *= cosmo->a2_inv;
+}
+
+void convert_part_potential(const struct engine* e, const struct part* p,
+ const struct xpart* xp, float* ret) {
+
+ if (p->gpart != NULL)
+ ret[0] = gravity_get_comoving_potential(p->gpart);
+ else
+ ret[0] = 0.f;
+}
+
/**
* @brief Specifies which particle fields to write to a dataset
*
@@ -98,16 +143,17 @@ void convert_part_pos(const struct engine* e, const struct part* p,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
-void hydro_write_particles(struct part* parts, struct io_props* list,
- int* num_fields) {
+void hydro_write_particles(const struct part* parts, const struct xpart* xparts,
+ struct io_props* list, int* num_fields) {
- *num_fields = 10;
+ *num_fields = 11;
/* List what we want to write */
- list[0] = io_make_output_field_convert_part(
- "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
- list[1] =
- io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
+ list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+ UNIT_CONV_LENGTH, parts, xparts,
+ convert_part_pos);
+ list[1] = io_make_output_field_convert_part(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
list[2] =
io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
@@ -120,11 +166,14 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
list[7] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
- parts, convert_u);
+ parts, xparts, convert_u);
list[8] = io_make_output_field_convert_part(
- "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
+ "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_P);
list[9] = io_make_output_field("WeightedDensity", FLOAT, 1, UNIT_CONV_DENSITY,
parts, rho_bar);
+ list[10] = io_make_output_field_convert_part("Potential", FLOAT, 1,
+ UNIT_CONV_POTENTIAL, parts,
+ xparts, convert_part_potential);
}
/**
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 31e18913b7414110db68c6202c512de5fb90a3c1..bd8328946d4bfd078bd10b439e96b35dce457c49 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -36,8 +36,13 @@
#define hydro_props_default_volume_change 1.4f
#define hydro_props_default_h_max FLT_MAX
#define hydro_props_default_h_tolerance 1e-4
+#define hydro_props_default_init_temp 0.f
+#define hydro_props_default_min_temp 0.f
+#define hydro_props_default_H_fraction 0.76
void hydro_props_init(struct hydro_props *p,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
const struct swift_params *params) {
/* Kernel properties */
@@ -73,6 +78,53 @@ void hydro_props_init(struct hydro_props *p,
const float max_volume_change = parser_get_opt_param_float(
params, "SPH:max_volume_change", hydro_props_default_volume_change);
p->log_max_h_change = logf(powf(max_volume_change, hydro_dimension_inv));
+
+ /* Initial temperature */
+ p->initial_temperature = parser_get_opt_param_float(
+ params, "SPH:initial_temperature", hydro_props_default_init_temp);
+
+ /* Initial temperature */
+ p->minimal_temperature = parser_get_opt_param_float(
+ params, "SPH:minimal_temperature", hydro_props_default_min_temp);
+
+ if ((p->initial_temperature != 0.) &&
+ (p->initial_temperature < p->minimal_temperature))
+ error("Initial temperature lower than minimal allowed temperature!");
+
+ /* Hydrogen mass fraction */
+ p->hydrogen_mass_fraction = parser_get_opt_param_double(
+ params, "SPH:H_mass_fraction", hydro_props_default_H_fraction);
+
+ /* Compute the initial energy (Note the temp. read is in internal units) */
+ double u_init = phys_const->const_boltzmann_k / phys_const->const_proton_mass;
+ u_init *= p->initial_temperature;
+ u_init *= hydro_one_over_gamma_minus_one;
+
+ /* Correct for hydrogen mass fraction */
+ double mean_molecular_weight;
+ if (p->initial_temperature *
+ units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE) >
+ 1e4)
+ mean_molecular_weight = 4. / (8. - 5. * (1. - p->hydrogen_mass_fraction));
+ else
+ mean_molecular_weight = 4. / (1. + 3. * p->hydrogen_mass_fraction);
+
+ p->initial_internal_energy = u_init / mean_molecular_weight;
+
+ /* Compute the minimal energy (Note the temp. read is in internal units) */
+ double u_min = phys_const->const_boltzmann_k / phys_const->const_proton_mass;
+ u_min *= p->initial_temperature;
+ u_min *= hydro_one_over_gamma_minus_one;
+
+ /* Correct for hydrogen mass fraction */
+ if (p->minimal_temperature *
+ units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE) >
+ 1e4)
+ mean_molecular_weight = 4. / (8. - 5. * (1. - p->hydrogen_mass_fraction));
+ else
+ mean_molecular_weight = 4. / (1. + 3. * p->hydrogen_mass_fraction);
+
+ p->minimal_internal_energy = u_min / mean_molecular_weight;
}
void hydro_props_print(const struct hydro_props *p) {
@@ -103,6 +155,12 @@ void hydro_props_print(const struct hydro_props *p) {
if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
message("Maximal iterations in ghost task set to %d (default is %d)",
p->max_smoothing_iterations, hydro_props_default_max_iterations);
+
+ if (p->initial_temperature != hydro_props_default_init_temp)
+ message("Initial gas temperature set to %f", p->initial_temperature);
+
+ if (p->minimal_temperature != hydro_props_default_min_temp)
+ message("Minimal gas temperature set to %f", p->minimal_temperature);
}
#if defined(HAVE_HDF5)
@@ -125,6 +183,12 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
pow_dimension(expf(p->log_max_h_change)));
io_write_attribute_i(h_grpsph, "Max ghost iterations",
p->max_smoothing_iterations);
+ io_write_attribute_f(h_grpsph, "Minimal temperature", p->minimal_temperature);
+ io_write_attribute_f(h_grpsph, "Initial temperature", p->initial_temperature);
+ io_write_attribute_f(h_grpsph, "Initial energy per unit mass",
+ p->initial_internal_energy);
+ io_write_attribute_f(h_grpsph, "Hydrogen mass fraction",
+ p->hydrogen_mass_fraction);
}
#endif
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
index 6d325c35d6ae6a9bcddbfff9ccc4601e026fc4e6..2799f6c86eec7f0ebf140643c5ab7fa9b60e6273 100644
--- a/src/hydro_properties.h
+++ b/src/hydro_properties.h
@@ -33,7 +33,9 @@
/* Local includes. */
#include "parser.h"
+#include "physical_constants.h"
#include "restart.h"
+#include "units.h"
/**
* @brief Contains all the constants and parameters of the hydro scheme
@@ -63,10 +65,28 @@ struct hydro_props {
/*! Maximal change of h over one time-step */
float log_max_h_change;
+
+ /*! Minimal temperature allowed */
+ float minimal_temperature;
+
+ /*! Minimal internal energy per unit mass */
+ float minimal_internal_energy;
+
+ /*! Initial temperature */
+ float initial_temperature;
+
+ /*! Initial internal energy per unit mass */
+ float initial_internal_energy;
+
+ /*! Primoridal hydrogen mass fraction for initial energy conversion */
+ float hydrogen_mass_fraction;
};
void hydro_props_print(const struct hydro_props *p);
-void hydro_props_init(struct hydro_props *p, const struct swift_params *params);
+void hydro_props_init(struct hydro_props *p,
+ const struct phys_const *phys_const,
+ const struct unit_system *us,
+ const struct swift_params *params);
#if defined(HAVE_HDF5)
void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p);
diff --git a/src/io_properties.h b/src/io_properties.h
index 17db123ba8325ede0c7af4e4b005e359cace7a7f..55e128e890b75db2d01911bd3c3b46388b0c6597 100644
--- a/src/io_properties.h
+++ b/src/io_properties.h
@@ -34,9 +34,11 @@ enum DATA_IMPORTANCE { COMPULSORY = 1, OPTIONAL = 0, UNUSED = -1 };
/* Helper typedefs */
typedef void (*conversion_func_part_float)(const struct engine*,
- const struct part*, float*);
+ const struct part*,
+ const struct xpart*, float*);
typedef void (*conversion_func_part_double)(const struct engine*,
- const struct part*, double*);
+ const struct part*,
+ const struct xpart*, double*);
typedef void (*conversion_func_gpart_float)(const struct engine*,
const struct gpart*, float*);
typedef void (*conversion_func_gpart_double)(const struct engine*,
@@ -78,6 +80,7 @@ struct io_props {
/* The particle arrays */
const struct part* parts;
+ const struct xpart* xparts;
const struct gpart* gparts;
/* Are we converting? */
@@ -125,6 +128,7 @@ INLINE static struct io_props io_make_input_field_(
r.field = field;
r.partSize = partSize;
r.parts = NULL;
+ r.xparts = NULL;
r.gparts = NULL;
r.conversion = 0;
r.convert_part_f = NULL;
@@ -179,10 +183,10 @@ INLINE static struct io_props io_make_output_field_(
/**
* @brief Constructs an #io_props (with conversion) from its parameters
*/
-#define io_make_output_field_convert_part(name, type, dim, units, part, \
- convert) \
- io_make_output_field_convert_part_##type(name, type, dim, units, \
- sizeof(part[0]), part, convert)
+#define io_make_output_field_convert_part(name, type, dim, units, part, xpart, \
+ convert) \
+ io_make_output_field_convert_part_##type( \
+ name, type, dim, units, sizeof(part[0]), part, xpart, convert)
/**
* @brief Construct an #io_props from its parameters
@@ -193,6 +197,7 @@ INLINE static struct io_props io_make_output_field_(
* @param units The units of the dataset
* @param partSize The size in byte of the particle
* @param parts The particle array
+ * @param xparts The xparticle array
* @param functionPtr The function used to convert a particle to a float
*
* Do not call this function directly. Use the macro defined above.
@@ -200,7 +205,8 @@ INLINE static struct io_props io_make_output_field_(
INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
enum unit_conversion_factor units, size_t partSize,
- const struct part* parts, conversion_func_part_float functionPtr) {
+ const struct part* parts, const struct xpart* xparts,
+ conversion_func_part_float functionPtr) {
struct io_props r;
strcpy(r.name, name);
@@ -211,6 +217,7 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
r.field = NULL;
r.partSize = partSize;
r.parts = parts;
+ r.xparts = xparts;
r.gparts = NULL;
r.conversion = 1;
r.convert_part_f = functionPtr;
@@ -230,6 +237,7 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
* @param units The units of the dataset
* @param partSize The size in byte of the particle
* @param parts The particle array
+ * @param xparts The xparticle array
* @param functionPtr The function used to convert a particle to a float
*
* Do not call this function directly. Use the macro defined above.
@@ -237,7 +245,8 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
INLINE static struct io_props io_make_output_field_convert_part_DOUBLE(
const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
enum unit_conversion_factor units, size_t partSize,
- const struct part* parts, conversion_func_part_double functionPtr) {
+ const struct part* parts, const struct xpart* xparts,
+ conversion_func_part_double functionPtr) {
struct io_props r;
strcpy(r.name, name);
@@ -248,6 +257,7 @@ INLINE static struct io_props io_make_output_field_convert_part_DOUBLE(
r.field = NULL;
r.partSize = partSize;
r.parts = parts;
+ r.xparts = xparts;
r.gparts = NULL;
r.conversion = 1;
r.convert_part_f = NULL;
@@ -293,6 +303,7 @@ INLINE static struct io_props io_make_output_field_convert_gpart_FLOAT(
r.field = NULL;
r.partSize = gpartSize;
r.parts = NULL;
+ r.xparts = NULL;
r.gparts = gparts;
r.conversion = 1;
r.convert_part_f = NULL;
@@ -330,6 +341,7 @@ INLINE static struct io_props io_make_output_field_convert_gpart_DOUBLE(
r.field = NULL;
r.partSize = gpartSize;
r.parts = NULL;
+ r.xparts = NULL;
r.gparts = gparts;
r.conversion = 1;
r.convert_part_f = NULL;
diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h
index 24d30136c8213dab6db424f2f459766c45ec6b1e..dc9db63f9d59ad0b3e72792b33b38a961ea6c30f 100644
--- a/src/kernel_gravity.h
+++ b/src/kernel_gravity.h
@@ -26,6 +26,10 @@
#include "inline.h"
#include "minmax.h"
+/*! Conversion factor between Plummer softening and internal softening */
+#define kernel_gravity_softening_plummer_equivalent 3.
+#define kernel_gravity_softening_plummer_equivalent_inv (1. / 3.)
+
/**
* @brief Computes the gravity softening function for potential.
*
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index 77bc7c17e2ab8038a4860c74daebd19a2d785070..578d7c2154f6801e21c004c01cade8bf6bd728ae 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -49,7 +49,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_pot_eval(
#else
const float x = 2.f * u;
- const float exp_x = expf(x);
+ const float exp_x = good_approx_expf(x);
const float alpha = 1.f / (1.f + exp_x);
/* We want 2 - 2 exp(x) * alpha */
@@ -83,7 +83,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
#else
const float arg = 2.f * u;
- const float exp_arg = expf(arg);
+ const float exp_arg = good_approx_expf(arg);
const float term = 1.f / (1.f + exp_arg);
*W = arg * exp_arg * term * term - exp_arg * term + 1.f;
@@ -111,4 +111,4 @@ __attribute__((always_inline)) INLINE static void fourier_kernel_long_grav_eval(
#endif
}
-#endif // SWIFT_KERNEL_LONG_GRAVITY_H
+#endif /* SWIFT_KERNEL_LONG_GRAVITY_H */
diff --git a/src/kick.h b/src/kick.h
index 9b1f4f18112a1cb2affcff70e7773ba4c48681b5..9d10f1e78d3934c4277c14217cbbc46514e87033 100644
--- a/src/kick.h
+++ b/src/kick.h
@@ -68,13 +68,16 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
* @param dt_kick_hydro The kick time-step for hydro accelerations.
* @param dt_kick_grav The kick time-step for gravity accelerations.
* @param dt_kick_therm The kick time-step for changes in thermal state.
+ * @param cosmo The cosmological model.
+ * @param hydro_props The constants used in the scheme
* @param ti_start The starting (integer) time of the kick (for debugging
* checks).
* @param ti_end The ending (integer) time of the kick (for debugging checks).
*/
__attribute__((always_inline)) INLINE static void kick_part(
struct part *restrict p, struct xpart *restrict xp, double dt_kick_hydro,
- double dt_kick_grav, double dt_kick_therm, integertime_t ti_start,
+ double dt_kick_grav, double dt_kick_therm, const struct cosmology *cosmo,
+ const struct hydro_props *hydro_props, integertime_t ti_start,
integertime_t ti_end) {
#ifdef SWIFT_DEBUG_CHECKS
@@ -107,7 +110,7 @@ __attribute__((always_inline)) INLINE static void kick_part(
}
/* Extra kick work */
- hydro_kick_extra(p, xp, dt_kick_therm);
+ hydro_kick_extra(p, xp, dt_kick_therm, cosmo, hydro_props);
if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt_kick_grav);
}
diff --git a/src/logger.c b/src/logger.c
index 38cb220b94dd14680584e9edca75883e1df16907..5fd4145aa1b042ed806dd3fe5487d094600b66c4 100644
--- a/src/logger.c
+++ b/src/logger.c
@@ -218,12 +218,6 @@ void logger_log_gpart(struct gpart *p, unsigned int mask, size_t *offset,
buff += 3 * sizeof(float);
}
- /* Particle smoothing length as a single float. */
- if (mask & logger_mask_h) {
- memcpy(buff, &p->epsilon, sizeof(float));
- buff += sizeof(float);
- }
-
/* Particle constants, which is a bit more complicated. */
if (mask & logger_mask_rho) {
memcpy(buff, &p->mass, sizeof(float));
@@ -385,12 +379,6 @@ int logger_read_gpart(struct gpart *p, size_t *offset, const char *buff) {
buff += 3 * sizeof(float);
}
- /* Particle smoothing length as a single float. */
- if (mask & logger_mask_h) {
- memcpy(&p->epsilon, buff, sizeof(float));
- buff += sizeof(float);
- }
-
/* Particle constants, which is a bit more complicated. */
if (mask & logger_mask_rho) {
memcpy(&p->mass, buff, sizeof(float));
diff --git a/src/multipole.h b/src/multipole.h
index 852976ad552b3ccf16f3e084612166c9212df2a6..dd74c4a40d79cc393e252bfe061591dd868f2a55 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1524,8 +1524,8 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
const double dim[3]) {
/* Recover some constants */
- const float eps = props->epsilon;
- const float eps_inv = props->epsilon_inv;
+ const float eps = props->epsilon_cur;
+ const float eps_inv = props->epsilon_cur_inv;
/* Compute distance vector */
float dx = (float)(pos_b[0] - pos_a[0]);
@@ -2380,6 +2380,7 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
gp->a_grav[2] += a_grav[2];
gp->potential += pot;
}
+
/**
* @brief Checks whether a cell-cell interaction can be appromixated by a M-M
* interaction using the distance and cell radius.
diff --git a/src/parallel_io.c b/src/parallel_io.c
index d3f5f8dc36dd63ad1b7cf4d65311751cd15f887b..51c88e95396eb6af1ab7f7eb2d1ed48eada59e0d 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -69,11 +69,14 @@
* @param offset Offset in the array where this mpi task starts writing.
* @param internal_units The #unit_system used internally.
* @param ic_units The #unit_system used in the snapshots.
+ * @param cleanup_h Are we removing h-factors from the ICs?
+ * @param h The value of the reduced Hubble constant to use for cleaning.
*/
void readArray_chunk(hid_t h_data, hid_t h_plist_id,
const struct io_props props, size_t N, long long offset,
const struct unit_system* internal_units,
- const struct unit_system* ic_units) {
+ const struct unit_system* ic_units, int cleanup_h,
+ double h) {
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
@@ -136,6 +139,23 @@ void readArray_chunk(hid_t h_data, hid_t h_plist_id,
}
}
+ /* Clean-up h if necessary */
+ const float h_factor_exp = units_h_factor(internal_units, props.units);
+ if (cleanup_h && h_factor_exp != 0.f) {
+ const double h_factor = pow(h, h_factor_exp);
+
+ /* message("Multipltying '%s' by h^%f=%f", props.name, h_factor_exp,
+ * h_factor); */
+
+ if (io_is_double_precision(props.type)) {
+ double* temp_d = (double*)temp;
+ for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= h_factor;
+ } else {
+ float* temp_f = (float*)temp;
+ for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= h_factor;
+ }
+ }
+
/* Copy temporary buffer to particle data */
char* temp_c = (char*)temp;
for (size_t i = 0; i < N; ++i)
@@ -158,11 +178,13 @@ void readArray_chunk(hid_t h_data, hid_t h_plist_id,
* @param offset The offset in the array on disk for this rank.
* @param internal_units The #unit_system used internally.
* @param ic_units The #unit_system used in the ICs.
+ * @param cleanup_h Are we removing h-factors from the ICs?
+ * @param h The value of the reduced Hubble constant to use for cleaning.
*/
void readArray(hid_t grp, struct io_props props, size_t N, long long N_total,
int mpi_rank, long long offset,
const struct unit_system* internal_units,
- const struct unit_system* ic_units) {
+ const struct unit_system* ic_units, int cleanup_h, double h) {
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
@@ -207,7 +229,7 @@ void readArray(hid_t grp, struct io_props props, size_t N, long long N_total,
/* Write the first chunk */
const size_t this_chunk = (N > max_chunk_size) ? max_chunk_size : N;
readArray_chunk(h_data, h_plist_id, props, this_chunk, offset,
- internal_units, ic_units);
+ internal_units, ic_units, cleanup_h, h);
/* Compute how many items are left */
if (N > max_chunk_size) {
@@ -309,8 +331,7 @@ void prepareArray(struct engine* e, hid_t grp, char* fileName, FILE* xmfFile,
io_write_attribute_d(
h_data, "CGS conversion factor",
units_cgs_conversion_factor(snapshot_units, props.units));
- io_write_attribute_f(h_data, "h-scale exponent",
- units_h_factor(snapshot_units, props.units));
+ io_write_attribute_f(h_data, "h-scale exponent", 0);
io_write_attribute_f(h_data, "a-scale exponent",
units_a_factor(snapshot_units, props.units));
io_write_attribute_s(h_data, "Conversion factor", buffer);
@@ -541,6 +562,8 @@ void writeArray(struct engine* e, hid_t grp, char* fileName,
* @param with_hydro Are we running with hydro ?
* @param with_gravity Are we running with gravity ?
* @param with_stars Are we running with stars ?
+ * @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
+ * @param h The value of the reduced Hubble constant to use for correction.
* @param mpi_rank The MPI rank of this node
* @param mpi_size The number of MPI ranks
* @param comm The MPI communicator
@@ -554,8 +577,9 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nstars, int* periodic, int* flag_entropy,
int with_hydro, int with_gravity, int with_stars,
- int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
- int n_threads, int dry_run) {
+ int cleanup_h, double h, int mpi_rank, int mpi_size,
+ MPI_Comm comm, MPI_Info info, int n_threads,
+ int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
@@ -626,6 +650,13 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
else if (hydro_dimension == 1)
dim[2] = dim[1] = dim[0];
+ /* Convert the box size if we want to clean-up h-factors */
+ if (cleanup_h) {
+ dim[0] /= h;
+ dim[1] /= h;
+ dim[2] /= h;
+ }
+
/* message("Found %lld particles in a %speriodic box of size [%f %f %f].", */
/* N_total[0], (periodic ? "": "non-"), dim[0], */
/* dim[1], dim[2]); */
@@ -771,7 +802,7 @@ void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
if (!dry_run)
for (int i = 0; i < num_fields; ++i)
readArray(h_grp, list[i], Nparticles, N_total[ptype], mpi_rank,
- offset[ptype], internal_units, ic_units);
+ offset[ptype], internal_units, ic_units, cleanup_h, h);
/* Close particle group */
H5Gclose(h_grp);
@@ -821,9 +852,10 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
const struct unit_system* internal_units,
const struct unit_system* snapshot_units) {
- struct part* parts = e->s->parts;
- struct gpart* gparts = e->s->gparts;
- struct spart* sparts = e->s->sparts;
+ const struct part* parts = e->s->parts;
+ const struct xpart* xparts = e->s->xparts;
+ const struct gpart* gparts = e->s->gparts;
+ const struct spart* sparts = e->s->sparts;
FILE* xmfFile = 0;
int periodic = e->s->periodic;
int numFiles = 1;
@@ -980,7 +1012,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
switch (ptype) {
case swift_type_gas:
- hydro_write_particles(parts, list, &num_fields);
+ hydro_write_particles(parts, xparts, list, &num_fields);
num_fields += chemistry_write_particles(parts, list + num_fields);
break;
@@ -1045,10 +1077,11 @@ void write_output_parallel(struct engine* e, const char* baseName,
const size_t Ngas = e->s->nr_parts;
const size_t Nstars = e->s->nr_sparts;
const size_t Ntot = e->s->nr_gparts;
- struct part* parts = e->s->parts;
- struct gpart* gparts = e->s->gparts;
+ const struct part* parts = e->s->parts;
+ const struct xpart* xparts = e->s->xparts;
+ const struct gpart* gparts = e->s->gparts;
struct gpart* dmparts = NULL;
- struct spart* sparts = e->s->sparts;
+ const struct spart* sparts = e->s->sparts;
/* Number of unassociated gparts */
const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
@@ -1208,7 +1241,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
case swift_type_gas:
Nparticles = Ngas;
- hydro_write_particles(parts, list, &num_fields);
+ hydro_write_particles(parts, xparts, list, &num_fields);
num_fields += chemistry_write_particles(parts, list + num_fields);
break;
diff --git a/src/parallel_io.h b/src/parallel_io.h
index 207e478c341488deb056761b37fb9c7370744b38..5ad3b34cdc4320d6fe3eb860615db33958ad612f 100644
--- a/src/parallel_io.h
+++ b/src/parallel_io.h
@@ -22,6 +22,8 @@
/* Config parameters. */
#include "../config.h"
+#if defined(HAVE_HDF5) && defined(WITH_MPI) && defined(HAVE_PARALLEL_HDF5)
+
/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
@@ -32,15 +34,14 @@
#include "part.h"
#include "units.h"
-#if defined(HAVE_HDF5) && defined(WITH_MPI) && defined(HAVE_PARALLEL_HDF5)
-
void read_ic_parallel(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nsparts, int* periodic, int* flag_entropy,
int with_hydro, int with_gravity, int with_stars,
- int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
- int nr_threads, int dry_run);
+ int cleanup_h, double h, int mpi_rank, int mpi_size,
+ MPI_Comm comm, MPI_Info info, int nr_threads,
+ int dry_run);
void write_output_parallel(struct engine* e, const char* baseName,
const struct unit_system* internal_units,
diff --git a/src/parser.h b/src/parser.h
index 4e61b16ab53b3688e60f85a42e786c44b095120a..58f4a53ea2114dbed29f49e2b5ccca69239b45e3 100644
--- a/src/parser.h
+++ b/src/parser.h
@@ -23,6 +23,9 @@
/* Config parameters. */
#include "../config.h"
+/* Standard headers */
+#include <stdio.h>
+
#if defined(HAVE_HDF5)
#include <hdf5.h>
#endif
diff --git a/src/part.c b/src/part.c
index b4c3a11bb0f7fbce4fbfbbbee198d1f88efdf6a3..1b696a8cbc135fd2c128b5ad705a0e6e24a2d5c8 100644
--- a/src/part.c
+++ b/src/part.c
@@ -139,8 +139,8 @@ void part_verify_links(struct part *parts, struct gpart *gparts,
if (gparts[k].type == swift_type_dark_matter) {
/* Check that it's not linked */
- if (gparts[k].id_or_neg_offset < 0)
- error("DM gpart particle linked to something!");
+ if (gparts[k].id_or_neg_offset <= 0)
+ error("DM gpart particle linked to something !");
}
/* We have a gas particle */
diff --git a/src/partition.c b/src/partition.c
index 8203c255d484672254b249042db22f5952579b4e..8feae3b7a8dce7c78310a3b35762d31b439c69e3 100644
--- a/src/partition.c
+++ b/src/partition.c
@@ -48,6 +48,7 @@
/* Local headers. */
#include "debug.h"
+#include "engine.h"
#include "error.h"
#include "partition.h"
#include "restart.h"
diff --git a/src/runner.c b/src/runner.c
index 1b8f9a2fc84b08c590396f1251a6bf3c5cc3189e..35996a3ba9a760d5f019ede35abc207bbdcbe1b3 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -981,6 +981,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
+ const struct hydro_props *hydro_props = e->hydro_properties;
const int with_cosmology = (e->policy & engine_policy_cosmology);
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
@@ -1044,8 +1045,8 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
}
/* do the kick */
- kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, ti_begin,
- ti_begin + ti_step / 2);
+ kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, cosmo,
+ hydro_props, ti_begin, ti_begin + ti_step / 2);
/* Update the accelerations to be used in the drift for hydro */
if (p->gpart != NULL) {
@@ -1150,6 +1151,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const struct cosmology *cosmo = e->cosmology;
+ const struct hydro_props *hydro_props = e->hydro_properties;
const int with_cosmology = (e->policy & engine_policy_cosmology);
const int count = c->count;
const int gcount = c->gcount;
@@ -1209,8 +1211,8 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
}
/* Finish the time-step with a second half-kick */
- kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm,
- ti_begin + ti_step / 2, ti_begin + ti_step);
+ kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, cosmo,
+ hydro_props, ti_begin + ti_step / 2, ti_begin + ti_step);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index cec450d84449d2b8cbb6790fcf9a596bb24649be..4f353a317e9f312f5ae6627cfaf0c06d176e8fad 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -466,7 +466,6 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
struct space *s = e->s;
const int periodic = s->periodic;
const double cell_width = s->width[0];
- const float theta_crit2 = e->gravity_properties->theta_crit2;
const double a_smooth = e->gravity_properties->a_smooth;
const double r_cut_min = e->gravity_properties->r_cut_min;
const double rlr = cell_width * a_smooth;
@@ -532,11 +531,11 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
/* Fill the caches */
gravity_cache_populate(e->max_active_bin, ci_cache, ci->gparts, gcount_i,
- gcount_padded_i, shift_i, CoM_j, rmax2_j, theta_crit2,
- ci);
+ gcount_padded_i, shift_i, CoM_j, rmax2_j, ci,
+ e->gravity_properties);
gravity_cache_populate(e->max_active_bin, cj_cache, cj->gparts, gcount_j,
- gcount_padded_j, shift_j, CoM_i, rmax2_i, theta_crit2,
- cj);
+ gcount_padded_j, shift_j, CoM_i, rmax2_i, cj,
+ e->gravity_properties);
/* Can we use the Newtonian version or do we need the truncated one ? */
if (!periodic) {
@@ -675,7 +674,7 @@ static INLINE void runner_doself_grav_pp_full(struct runner *r,
const int gcount_padded = gcount - (gcount % VEC_SIZE) + VEC_SIZE;
gravity_cache_populate_no_mpole(e->max_active_bin, ci_cache, gparts, gcount,
- gcount_padded, loc, c);
+ gcount_padded, loc, c, e->gravity_properties);
/* Ok... Here we go ! */
@@ -805,7 +804,7 @@ static INLINE void runner_doself_grav_pp_truncated(struct runner *r,
const int gcount_padded = gcount - (gcount % VEC_SIZE) + VEC_SIZE;
gravity_cache_populate_no_mpole(e->max_active_bin, ci_cache, gparts, gcount,
- gcount_padded, loc, c);
+ gcount_padded, loc, c, e->gravity_properties);
/* Ok... Here we go ! */
diff --git a/src/runner_doiact_vec.c b/src/runner_doiact_vec.c
index 0daea0cf8a8d868ccb14f7d9b27f04aca5a677c8..2e86280d64491ee1750f41c2cd22ab01c08e30b8 100644
--- a/src/runner_doiact_vec.c
+++ b/src/runner_doiact_vec.c
@@ -1689,7 +1689,7 @@ void runner_dopair_subset_density_vec(struct runner *r,
struct cell *restrict cj, const int sid,
const int flipped, const double *shift) {
-#ifdef WITH_VECTORIZATION
+#if defined(WITH_VECTORIZATION) && defined(GADGET2_SPH)
TIMER_TIC;
diff --git a/src/serial_io.c b/src/serial_io.c
index aa8d126e4c6f3e19aed7d205c7ced39f7eded6f8..dd623f946ce6ec32415586e5048979de3adb58fa 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -67,6 +67,8 @@
* @param offset The offset position where this rank starts reading.
* @param internal_units The #unit_system used internally
* @param ic_units The #unit_system used in the ICs
+ * @param cleanup_h Are we removing h-factors from the ICs?
+ * @param h The value of the reduced Hubble constant to use for cleaning.
*
* @todo A better version using HDF5 hyper-slabs to read the file directly into
* the part array will be written once the structures have been stabilized.
@@ -74,7 +76,7 @@
void readArray(hid_t grp, const struct io_props props, size_t N,
long long N_total, long long offset,
const struct unit_system* internal_units,
- const struct unit_system* ic_units) {
+ const struct unit_system* ic_units, int cleanup_h, double h) {
const size_t typeSize = io_sizeof_type(props.type);
const size_t copySize = typeSize * props.dimension;
@@ -161,6 +163,23 @@ void readArray(hid_t grp, const struct io_props props, size_t N,
}
}
+ /* Clean-up h if necessary */
+ const float h_factor_exp = units_h_factor(internal_units, props.units);
+ if (cleanup_h && h_factor_exp != 0.f && exist != 0) {
+ const double h_factor = pow(h, h_factor_exp);
+
+ /* message("Multipltying '%s' by h^%f=%f", props.name, h_factor_exp,
+ * h_factor); */
+
+ if (io_is_double_precision(props.type)) {
+ double* temp_d = (double*)temp;
+ for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= h_factor;
+ } else {
+ float* temp_f = (float*)temp;
+ for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= h_factor;
+ }
+ }
+
/* Copy temporary buffer to particle data */
char* temp_c = temp;
for (size_t i = 0; i < N; ++i)
@@ -256,8 +275,7 @@ void prepareArray(const struct engine* e, hid_t grp, char* fileName,
io_write_attribute_d(
h_data, "CGS conversion factor",
units_cgs_conversion_factor(snapshot_units, props.units));
- io_write_attribute_f(h_data, "h-scale exponent",
- units_h_factor(snapshot_units, props.units));
+ io_write_attribute_f(h_data, "h-scale exponent", 0);
io_write_attribute_f(h_data, "a-scale exponent",
units_a_factor(snapshot_units, props.units));
io_write_attribute_s(h_data, "Conversion factor", buffer);
@@ -382,6 +400,8 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
* @param with_hydro Are we reading gas particles ?
* @param with_gravity Are we reading/creating #gpart arrays ?
* @param with_stars Are we reading star particles ?
+ * @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
+ * @param h The value of the reduced Hubble constant to use for correction.
* @param mpi_rank The MPI rank of this node
* @param mpi_size The number of MPI ranks
* @param comm The MPI communicator
@@ -402,8 +422,8 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nstars, int* periodic, int* flag_entropy,
int with_hydro, int with_gravity, int with_stars,
- int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
- int n_threads, int dry_run) {
+ int cleanup_h, double h, int mpi_rank, int mpi_size,
+ MPI_Comm comm, MPI_Info info, int n_threads, int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
@@ -476,6 +496,13 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
else if (hydro_dimension == 1)
dim[2] = dim[1] = dim[0];
+ /* Convert the box size if we want to clean-up h-factors */
+ if (cleanup_h) {
+ dim[0] /= h;
+ dim[1] /= h;
+ dim[2] /= h;
+ }
+
/* message("Found %lld particles in a %speriodic box of size [%f %f %f].",
*/
/* N_total, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
@@ -641,7 +668,7 @@ void read_ic_serial(char* fileName, const struct unit_system* internal_units,
if (!dry_run)
for (int i = 0; i < num_fields; ++i)
readArray(h_grp, list[i], Nparticles, N_total[ptype], offset[ptype],
- internal_units, ic_units);
+ internal_units, ic_units, cleanup_h, h);
/* Close particle group */
H5Gclose(h_grp);
@@ -712,10 +739,11 @@ void write_output_serial(struct engine* e, const char* baseName,
const size_t Ntot = e->s->nr_gparts;
int periodic = e->s->periodic;
int numFiles = 1;
- struct part* parts = e->s->parts;
- struct gpart* gparts = e->s->gparts;
+ const struct part* parts = e->s->parts;
+ const struct xpart* xparts = e->s->xparts;
+ const struct gpart* gparts = e->s->gparts;
struct gpart* dmparts = NULL;
- struct spart* sparts = e->s->sparts;
+ const struct spart* sparts = e->s->sparts;
FILE* xmfFile = 0;
/* Number of unassociated gparts */
@@ -962,7 +990,7 @@ void write_output_serial(struct engine* e, const char* baseName,
case swift_type_gas:
Nparticles = Ngas;
- hydro_write_particles(parts, list, &num_fields);
+ hydro_write_particles(parts, xparts, list, &num_fields);
num_fields += chemistry_write_particles(parts, list + num_fields);
break;
diff --git a/src/serial_io.h b/src/serial_io.h
index 025faffc0087a8b26d4149914e621a978e63c9be..cad428916f400bc9b144dcf4c23f9d38c75c1e9d 100644
--- a/src/serial_io.h
+++ b/src/serial_io.h
@@ -22,6 +22,8 @@
/* Config parameters. */
#include "../config.h"
+#if defined(HAVE_HDF5) && defined(WITH_MPI) && !defined(HAVE_PARALLEL_HDF5)
+
/* MPI headers. */
#ifdef WITH_MPI
#include <mpi.h>
@@ -32,15 +34,13 @@
#include "part.h"
#include "units.h"
-#if defined(HAVE_HDF5) && defined(WITH_MPI) && !defined(HAVE_PARALLEL_HDF5)
-
void read_ic_serial(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nstars, int* periodic, int* flag_entropy,
int with_hydro, int with_gravity, int with_stars,
- int mpi_rank, int mpi_size, MPI_Comm comm, MPI_Info info,
- int nr_threads, int dry_run);
+ int cleanup_h, double h, int mpi_rank, int mpi_size,
+ MPI_Comm comm, MPI_Info info, int nr_threads, int dry_run);
void write_output_serial(struct engine* e, const char* baseName,
const struct unit_system* internal_units,
diff --git a/src/single_io.c b/src/single_io.c
index 0fecd724908f7ce2779526792ce8512b9236d58e..2170bcffc4ce3ab21f1edd168d1dc37b2b4af963 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -64,13 +64,15 @@
* @param N The number of particles.
* @param internal_units The #unit_system used internally
* @param ic_units The #unit_system used in the ICs
+ * @param cleanup_h Are we removing h-factors from the ICs?
+ * @param The value of the reduced Hubble constant.
*
* @todo A better version using HDF5 hyper-slabs to read the file directly into
* the part array will be written once the structures have been stabilized.
*/
void readArray(hid_t h_grp, const struct io_props prop, size_t N,
const struct unit_system* internal_units,
- const struct unit_system* ic_units) {
+ const struct unit_system* ic_units, int cleanup_h, double h) {
const size_t typeSize = io_sizeof_type(prop.type);
const size_t copySize = typeSize * prop.dimension;
@@ -124,18 +126,35 @@ void readArray(hid_t h_grp, const struct io_props prop, size_t N,
}
/* Unit conversion if necessary */
- const double factor =
+ const double unit_factor =
units_conversion_factor(ic_units, internal_units, prop.units);
- if (factor != 1. && exist != 0) {
+ if (unit_factor != 1. && exist != 0) {
/* message("Converting ! factor=%e", factor); */
if (io_is_double_precision(prop.type)) {
double* temp_d = (double*)temp;
- for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= factor;
+ for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= unit_factor;
} else {
float* temp_f = (float*)temp;
- for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= factor;
+ for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= unit_factor;
+ }
+ }
+
+ /* Clean-up h if necessary */
+ const float h_factor_exp = units_h_factor(internal_units, prop.units);
+ if (cleanup_h && h_factor_exp != 0.f && exist != 0) {
+ const double h_factor = pow(h, h_factor_exp);
+
+ /* message("Multipltying '%s' by h^%f=%f", prop.name, h_factor_exp,
+ * h_factor); */
+
+ if (io_is_double_precision(prop.type)) {
+ double* temp_d = (double*)temp;
+ for (size_t i = 0; i < num_elements; ++i) temp_d[i] *= h_factor;
+ } else {
+ float* temp_f = (float*)temp;
+ for (size_t i = 0; i < num_elements; ++i) temp_f[i] *= h_factor;
}
}
@@ -270,8 +289,7 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
io_write_attribute_d(
h_data, "CGS conversion factor",
units_cgs_conversion_factor(snapshot_units, props.units));
- io_write_attribute_f(h_data, "h-scale exponent",
- units_h_factor(snapshot_units, props.units));
+ io_write_attribute_f(h_data, "h-scale exponent", 0);
io_write_attribute_f(h_data, "a-scale exponent",
units_a_factor(snapshot_units, props.units));
io_write_attribute_s(h_data, "Conversion factor", buffer);
@@ -301,6 +319,9 @@ void writeArray(const struct engine* e, hid_t grp, char* fileName,
* @param with_hydro Are we reading gas particles ?
* @param with_gravity Are we reading/creating #gpart arrays ?
* @param with_stars Are we reading star particles ?
+ * @param cleanup_h Are we cleaning-up h-factors from the quantities we read?
+ * @param h The value of the reduced Hubble constant to use for correction.
+ * @prarm n_threads The number of threads to use for the temporary threadpool.
* @param dry_run If 1, don't read the particle. Only allocates the arrays.
*
* Opens the HDF5 file fileName and reads the particles contained
@@ -316,7 +337,7 @@ void read_ic_single(char* fileName, const struct unit_system* internal_units,
struct spart** sparts, size_t* Ngas, size_t* Ngparts,
size_t* Nstars, int* periodic, int* flag_entropy,
int with_hydro, int with_gravity, int with_stars,
- int n_threads, int dry_run) {
+ int cleanup_h, double h, int n_threads, int dry_run) {
hid_t h_file = 0, h_grp = 0;
/* GADGET has only cubic boxes (in cosmological mode) */
@@ -384,6 +405,13 @@ void read_ic_single(char* fileName, const struct unit_system* internal_units,
else if (hydro_dimension == 1)
dim[2] = dim[1] = dim[0];
+ /* Convert the box size if we want to clean-up h-factors */
+ if (cleanup_h) {
+ dim[0] /= h;
+ dim[1] /= h;
+ dim[2] /= h;
+ }
+
/* message("Found %d particles in a %speriodic box of size [%f %f %f].", */
/* *N, (periodic ? "": "non-"), dim[0], dim[1], dim[2]); */
@@ -517,7 +545,8 @@ void read_ic_single(char* fileName, const struct unit_system* internal_units,
/* Read everything */
if (!dry_run)
for (int i = 0; i < num_fields; ++i)
- readArray(h_grp, list[i], Nparticles, internal_units, ic_units);
+ readArray(h_grp, list[i], Nparticles, internal_units, ic_units,
+ cleanup_h, h);
/* Close particle group */
H5Gclose(h_grp);
@@ -578,10 +607,11 @@ void write_output_single(struct engine* e, const char* baseName,
const size_t Ntot = e->s->nr_gparts;
int periodic = e->s->periodic;
int numFiles = 1;
- struct part* parts = e->s->parts;
- struct gpart* gparts = e->s->gparts;
+ const struct part* parts = e->s->parts;
+ const struct xpart* xparts = e->s->xparts;
+ const struct gpart* gparts = e->s->gparts;
struct gpart* dmparts = NULL;
- struct spart* sparts = e->s->sparts;
+ const struct spart* sparts = e->s->sparts;
/* Number of unassociated gparts */
const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
@@ -779,7 +809,7 @@ void write_output_single(struct engine* e, const char* baseName,
case swift_type_gas:
N = Ngas;
- hydro_write_particles(parts, list, &num_fields);
+ hydro_write_particles(parts, xparts, list, &num_fields);
num_fields += chemistry_write_particles(parts, list + num_fields);
break;
diff --git a/src/single_io.h b/src/single_io.h
index efed4ef27ade61f7726777d71a2fdc586ea26030..26b849716e3e018d9a10c5c5c513ad26c7ccb274 100644
--- a/src/single_io.h
+++ b/src/single_io.h
@@ -22,24 +22,24 @@
/* Config parameters. */
#include "../config.h"
+#if defined(HAVE_HDF5) && !defined(WITH_MPI)
+
/* Includes. */
#include "engine.h"
#include "part.h"
#include "units.h"
-#if defined(HAVE_HDF5) && !defined(WITH_MPI)
-
void read_ic_single(char* fileName, const struct unit_system* internal_units,
double dim[3], struct part** parts, struct gpart** gparts,
struct spart** sparts, size_t* Ngas, size_t* Ndm,
size_t* Nstars, int* periodic, int* flag_entropy,
int with_hydro, int with_gravity, int with_stars,
- int nr_threads, int dry_run);
+ int cleanup_h, double h, int nr_threads, int dry_run);
void write_output_single(struct engine* e, const char* baseName,
const struct unit_system* internal_units,
const struct unit_system* snapshot_units);
-#endif
+#endif /* HAVE_HDF5 && !WITH_MPI */
#endif /* SWIFT_SINGLE_IO_H */
diff --git a/src/space.c b/src/space.c
index 31e4250c83d86314d38c3dcf739ef179ef813235..bac6240512edbf2ad31a3216bcd518d56e3573a4 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2321,8 +2321,20 @@ void space_first_init_parts(struct space *s,
struct part *restrict p = s->parts;
struct xpart *restrict xp = s->xparts;
+ const struct cosmology *cosmo = s->e->cosmology;
+ const float a_factor_vel = cosmo->a * cosmo->a;
+
+ const struct hydro_props *hydro_props = s->e->hydro_properties;
+ const float u_init = hydro_props->initial_internal_energy;
+ const float u_min = hydro_props->minimal_internal_energy;
+
for (size_t i = 0; i < nr_parts; ++i) {
+ /* Convert velocities to internal units */
+ p[i].v[0] *= a_factor_vel;
+ p[i].v[1] *= a_factor_vel;
+ p[i].v[2] *= a_factor_vel;
+
#ifdef HYDRO_DIMENSION_2D
p[i].x[2] = 0.f;
p[i].v[2] = 0.f;
@@ -2335,6 +2347,10 @@ void space_first_init_parts(struct space *s,
hydro_first_init_part(&p[i], &xp[i]);
+ /* Overwrite the internal energy? */
+ if (u_init > 0.f) hydro_set_init_internal_energy(&p[i], u_init);
+ if (u_min > 0.f) hydro_set_init_internal_energy(&p[i], u_min);
+
/* Also initialise the chemistry */
chemistry_first_init_part(&p[i], &xp[i], chemistry);
@@ -2361,9 +2377,16 @@ void space_first_init_gparts(struct space *s,
const size_t nr_gparts = s->nr_gparts;
struct gpart *restrict gp = s->gparts;
+ const struct cosmology *cosmo = s->e->cosmology;
+ const float a_factor_vel = cosmo->a * cosmo->a;
for (size_t i = 0; i < nr_gparts; ++i) {
+ /* Convert velocities to internal units */
+ gp[i].v_full[0] *= a_factor_vel;
+ gp[i].v_full[1] *= a_factor_vel;
+ gp[i].v_full[2] *= a_factor_vel;
+
#ifdef HYDRO_DIMENSION_2D
gp[i].x[2] = 0.f;
gp[i].v_full[2] = 0.f;
@@ -2392,9 +2415,16 @@ void space_first_init_sparts(struct space *s) {
const size_t nr_sparts = s->nr_sparts;
struct spart *restrict sp = s->sparts;
+ const struct cosmology *cosmo = s->e->cosmology;
+ const float a_factor_vel = cosmo->a * cosmo->a;
for (size_t i = 0; i < nr_sparts; ++i) {
+ /* Convert velocities to internal units */
+ sp[i].v[0] *= a_factor_vel;
+ sp[i].v[1] *= a_factor_vel;
+ sp[i].v[2] *= a_factor_vel;
+
#ifdef HYDRO_DIMENSION_2D
sp[i].x[2] = 0.f;
sp[i].v[2] = 0.f;
@@ -2472,11 +2502,12 @@ void space_init_gparts(struct space *s, int verbose) {
void space_convert_quantities_mapper(void *restrict map_data, int count,
void *restrict extra_data) {
struct space *s = (struct space *)extra_data;
+ const struct cosmology *cosmo = s->e->cosmology;
struct part *restrict parts = (struct part *)map_data;
const ptrdiff_t index = parts - s->parts;
struct xpart *restrict xparts = s->xparts + index;
for (int k = 0; k < count; k++)
- hydro_convert_quantities(&parts[k], &xparts[k]);
+ hydro_convert_quantities(&parts[k], &xparts[k], cosmo);
}
/**
@@ -2504,6 +2535,7 @@ void space_convert_quantities(struct space *s, int verbose) {
*
* @param s The #space to initialize.
* @param params The parsed parameter file.
+ * @param cosmo The current cosmological model.
* @param dim Spatial dimensions of the domain.
* @param parts Array of Gas particles.
* @param gparts Array of Gravity particles.
@@ -2512,8 +2544,9 @@ void space_convert_quantities(struct space *s, int verbose) {
* @param Ngpart The number of Gravity particles in the space.
* @param Nspart The number of star particles in the space.
* @param periodic flag whether the domain is periodic or not.
- * @param replicate How many replications along each direction do we want ?
- * @param gravity flag whether we are doing gravity or not.
+ * @param replicate How many replications along each direction do we want?
+ * @param generate_gas_in_ics Are we generating gas particles from the gparts?
+ * @param self_gravity flag whether we are doing gravity or not?
* @param verbose Print messages to stdout or not.
* @param dry_run If 1, just initialise stuff, don't do anything with the parts.
*
@@ -2523,9 +2556,10 @@ void space_convert_quantities(struct space *s, int verbose) {
* recursively.
*/
void space_init(struct space *s, const struct swift_params *params,
- double dim[3], struct part *parts, struct gpart *gparts,
- struct spart *sparts, size_t Npart, size_t Ngpart,
- size_t Nspart, int periodic, int replicate, int gravity,
+ const struct cosmology *cosmo, double dim[3],
+ struct part *parts, struct gpart *gparts, struct spart *sparts,
+ size_t Npart, size_t Ngpart, size_t Nspart, int periodic,
+ int replicate, int generate_gas_in_ics, int self_gravity,
int verbose, int dry_run) {
/* Clean-up everything */
@@ -2536,7 +2570,7 @@ void space_init(struct space *s, const struct swift_params *params,
s->dim[1] = dim[1];
s->dim[2] = dim[2];
s->periodic = periodic;
- s->gravity = gravity;
+ s->gravity = self_gravity;
s->nr_parts = Npart;
s->size_parts = Npart;
s->parts = parts;
@@ -2548,6 +2582,19 @@ void space_init(struct space *s, const struct swift_params *params,
s->sparts = sparts;
s->nr_queues = 1; /* Temporary value until engine construction */
+ /* Are we generating gas from the DM-only ICs? */
+ if (generate_gas_in_ics) {
+ space_generate_gas(s, cosmo, verbose);
+ parts = s->parts;
+ gparts = s->gparts;
+ Npart = s->nr_parts;
+ Ngpart = s->nr_gparts;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ part_verify_links(parts, gparts, sparts, Npart, Ngpart, Nspart, 1);
+#endif
+ }
+
/* Are we replicating the space ? */
if (replicate < 1)
error("Value of 'InitialConditions:replicate' (%d) is too small",
@@ -2560,6 +2607,10 @@ void space_init(struct space *s, const struct swift_params *params,
Npart = s->nr_parts;
Ngpart = s->nr_gparts;
Nspart = s->nr_sparts;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ part_verify_links(parts, gparts, sparts, Npart, Ngpart, Nspart, 1);
+#endif
}
/* Decide on the minimal top-level cell size */
@@ -2606,8 +2657,8 @@ void space_init(struct space *s, const struct swift_params *params,
}
/* Apply h scaling */
- const double scaling =
- parser_get_opt_param_double(params, "InitialConditions:h_scaling", 1.0);
+ const double scaling = parser_get_opt_param_double(
+ params, "InitialConditions:smoothing_length_scaling", 1.0);
if (scaling != 1.0 && !dry_run) {
message("Re-scaling smoothing lengths by a factor %e", scaling);
for (size_t k = 0; k < Npart; k++) parts[k].h *= scaling;
@@ -2824,6 +2875,152 @@ void space_replicate(struct space *s, int replicate, int verbose) {
#endif
}
+void space_generate_gas(struct space *s, const struct cosmology *cosmo,
+ int verbose) {
+
+ if (verbose) message("Generating gas particles from gparts");
+
+ /* Store the current values */
+ const size_t nr_parts = s->nr_parts;
+ const size_t nr_gparts = s->nr_gparts;
+
+ if (nr_parts != 0)
+ error("Generating gas particles from DM but gas already exists!");
+
+ s->size_parts = s->nr_parts = nr_gparts;
+ s->size_gparts = s->nr_gparts = 2 * nr_gparts;
+
+ /* Allocate space for new particles */
+ struct part *parts = NULL;
+ struct gpart *gparts = NULL;
+
+ if (posix_memalign((void **)&parts, part_align,
+ s->nr_parts * sizeof(struct part)) != 0)
+ error("Failed to allocate new part array.");
+
+ if (posix_memalign((void **)&gparts, gpart_align,
+ s->nr_gparts * sizeof(struct gpart)) != 0)
+ error("Failed to allocate new gpart array.");
+
+ /* Start by copying the gparts */
+ memcpy(gparts, s->gparts, nr_gparts * sizeof(struct gpart));
+ memcpy(gparts + nr_gparts, s->gparts, nr_gparts * sizeof(struct gpart));
+
+ /* And zero the parts */
+ bzero(parts, s->nr_parts * sizeof(struct part));
+
+ /* Compute some constants */
+ const double mass_ratio = cosmo->Omega_b / cosmo->Omega_m;
+ const double bg_density = cosmo->Omega_m * cosmo->critical_density;
+ const double bg_density_inv = 1. / bg_density;
+
+ /* Update the particle properties */
+ for (size_t i = 0; i < nr_gparts; ++i) {
+
+ struct part *p = &parts[i];
+ struct gpart *gp_gas = &gparts[i];
+ struct gpart *gp_dm = &gparts[nr_gparts + i];
+
+ /* Set the IDs */
+ p->id = gp_gas->id_or_neg_offset * 2 + 1;
+ gp_dm->id_or_neg_offset *= 2;
+
+ if (gp_dm->id_or_neg_offset <= 0) error("DM particle ID overflowd");
+
+ if (p->id <= 0) error("gas particle ID overflowd");
+
+ /* Set the links correctly */
+ p->gpart = gp_gas;
+ gp_gas->id_or_neg_offset = -i;
+ gp_gas->type = swift_type_gas;
+
+ /* Compute positions shift */
+ const double d = cbrt(gp_dm->mass * bg_density_inv);
+ const double shift_dm = d * mass_ratio;
+ const double shift_gas = d * (1. - mass_ratio);
+
+ /* Set the masses */
+ gp_dm->mass *= (1. - mass_ratio);
+ gp_gas->mass *= mass_ratio;
+ hydro_set_mass(p, gp_gas->mass);
+
+ /* Set the new positions */
+ gp_dm->x[0] += shift_dm;
+ gp_dm->x[1] += shift_dm;
+ gp_dm->x[2] += shift_dm;
+ gp_gas->x[0] += shift_gas;
+ gp_gas->x[1] += shift_gas;
+ gp_gas->x[2] += shift_gas;
+ p->x[0] = gp_gas->x[0];
+ p->x[1] = gp_gas->x[1];
+ p->x[2] = gp_gas->x[2];
+
+ /* Also copy the velocities */
+ p->v[0] = gp_gas->v_full[0];
+ p->v[1] = gp_gas->v_full[1];
+ p->v[2] = gp_gas->v_full[2];
+
+ /* Set the smoothing length to the mean inter-particle separation */
+ p->h = 30. * d;
+ }
+
+ /* Replace the content of the space */
+ free(s->gparts);
+ s->parts = parts;
+ s->gparts = gparts;
+}
+
+/**
+ * @brief Verify that the matter content matches the cosmology model.
+ *
+ * @param s The #space.
+ * @param cosmo The current cosmology model.
+ * @param rank The MPI rank of this #space.
+ */
+void space_check_cosmology(struct space *s, const struct cosmology *cosmo,
+ int rank) {
+
+ struct gpart *gparts = s->gparts;
+ const size_t nr_gparts = s->nr_gparts;
+
+ /* Sum up the mass in this space */
+ double mass = 0.;
+ for (size_t i = 0; i < nr_gparts; ++i) {
+ mass += gparts[i].mass;
+ }
+
+/* Reduce the total mass */
+#ifdef WITH_MPI
+ double total_mass;
+ MPI_Reduce(&mass, &total_mass, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
+#else
+ double total_mass = mass;
+#endif
+
+ if (rank == 0) {
+
+ const double volume = s->dim[0] * s->dim[1] * s->dim[2];
+
+ /* Current Hubble constant */
+ const double H = cosmo->H;
+
+ /* z=0 Hubble parameter */
+ const double H0 = cosmo->H0;
+
+ /* Critical density at z=0 */
+ const double rho_crit0 = cosmo->critical_density * H0 * H0 / (H * H);
+
+ /* Compute the mass density */
+ const double Omega_m = (total_mass / volume) / rho_crit0;
+
+ if (fabs(Omega_m - cosmo->Omega_m) > 1e-3)
+ error(
+ "The matter content of the simulation does not match the cosmology "
+ "in the parameter file comso.Omega_m=%e Omega_m=%e",
+ cosmo->Omega_m, Omega_m);
+ }
+}
+
/**
* @brief Cleans-up all the cell links in the space
*
diff --git a/src/space.h b/src/space.h
index c8277d51ed70204efbd52b8688888c83b05a6979..571b595eda8758c510a014bc784ae4bde85abaea 100644
--- a/src/space.h
+++ b/src/space.h
@@ -38,6 +38,7 @@
/* Avoid cyclic inclusions */
struct cell;
+struct cosmology;
/* Some constants. */
#define space_cellallocchunk 1000
@@ -184,9 +185,10 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells);
int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
double *shift);
void space_init(struct space *s, const struct swift_params *params,
- double dim[3], struct part *parts, struct gpart *gparts,
- struct spart *sparts, size_t Npart, size_t Ngpart,
- size_t Nspart, int periodic, int replicate, int gravity,
+ const struct cosmology *cosmo, double dim[3],
+ struct part *parts, struct gpart *gparts, struct spart *sparts,
+ size_t Npart, size_t Ngpart, size_t Nspart, int periodic,
+ int replicate, int generate_gas_in_ics, int gravity,
int verbose, int dry_run);
void space_sanitize(struct space *s);
void space_map_cells_pre(struct space *s, int full,
@@ -235,6 +237,10 @@ void space_check_top_multipoles_drift_point(struct space *s,
integertime_t ti_drift);
void space_check_timesteps(struct space *s);
void space_replicate(struct space *s, int replicate, int verbose);
+void space_generate_gas(struct space *s, const struct cosmology *cosmo,
+ int verbose);
+void space_check_cosmology(struct space *s, const struct cosmology *cosmo,
+ int rank);
void space_reset_task_counters(struct space *s);
void space_clean(struct space *s);
void space_free_cells(struct space *s);
diff --git a/src/stars/Default/star_io.h b/src/stars/Default/star_io.h
index 96bbdce6d83dc241d05e7dd1754f476dc0b8e5f9..c3dc31096383533e1e15fa65615d2c9aac0f43e3 100644
--- a/src/stars/Default/star_io.h
+++ b/src/stars/Default/star_io.h
@@ -52,7 +52,7 @@ void star_read_particles(struct spart* sparts, struct io_props* list,
* @param list The list of i/o properties to write.
* @param num_fields The number of i/o fields to write.
*/
-void star_write_particles(struct spart* sparts, struct io_props* list,
+void star_write_particles(const struct spart* sparts, struct io_props* list,
int* num_fields) {
/* Say how much we want to read */
diff --git a/src/units.c b/src/units.c
index 5c50bb063b9411f47357407678458e6094f69e2b..4b632e735b7c6e1c12afe8aebc16aa44abc5b597 100644
--- a/src/units.c
+++ b/src/units.c
@@ -220,6 +220,11 @@ void units_get_base_unit_exponants_array(float baseUnitsExp[5],
baseUnitsExp[UNIT_TIME] = -2.f;
break;
+ case UNIT_CONV_POTENTIAL:
+ baseUnitsExp[UNIT_LENGTH] = 2.f;
+ baseUnitsExp[UNIT_TIME] = -2.f;
+ break;
+
case UNIT_CONV_FORCE:
baseUnitsExp[UNIT_MASS] = 1.f;
baseUnitsExp[UNIT_LENGTH] = 1.f;
diff --git a/src/units.h b/src/units.h
index 5ac70a909a77146ba5f7a441d7747acfc80c3dfa..87c44cc6eb4934980027a60642dd135f03029f7c 100644
--- a/src/units.h
+++ b/src/units.h
@@ -73,6 +73,7 @@ enum unit_conversion_factor {
UNIT_CONV_DENSITY,
UNIT_CONV_SPEED,
UNIT_CONV_ACCELERATION,
+ UNIT_CONV_POTENTIAL,
UNIT_CONV_FORCE,
UNIT_CONV_ENERGY,
UNIT_CONV_ENERGY_PER_UNIT_MASS,
diff --git a/tests/fft_params.yml b/tests/fft_params.yml
index 05d6d8f0b0578d11645fc1d78c1a6322160ae87a..6938e3658148874e58f50ab768a5e1fbc41d9573 100644
--- a/tests/fft_params.yml
+++ b/tests/fft_params.yml
@@ -3,8 +3,9 @@ Scheduler:
# Parameters for the self-gravity scheme
Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 0.7 # Opening angle (Multipole acceptance criterion)
- epsilon: 0.00001 # Softening length (in internal units).
- a_smooth: 0.
- r_cut: 0.
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
+ comoving_softening: 0.00001 # Comoving softening length (in internal units).
+ max_physical_softening: 0.00001 # Physical softening length (in internal units).
+ a_smooth: 0.
+ r_cut: 0.
diff --git a/tests/testFFT.c b/tests/testFFT.c
index ba737935afc4568a1509d2c43a3534b60a6ef867..7b67181ebd3e29bffbf564d00f702e6c15669fab 100644
--- a/tests/testFFT.c
+++ b/tests/testFFT.c
@@ -57,15 +57,18 @@ int main() {
struct swift_params *params = malloc(sizeof(struct swift_params));
parser_read_file("fft_params.yml", params);
+ struct cosmology cosmo;
+ cosmology_init_no_cosmo(&cosmo);
+
/* Initialise the gravity properties */
struct gravity_props gravity_properties;
- gravity_props_init(&gravity_properties, params);
+ gravity_props_init(&gravity_properties, params, &cosmo);
/* Build the infrastructure */
struct space space;
double dim[3] = {1., 1., 1.};
- space_init(&space, params, dim, NULL, gparts, NULL, 0, nr_gparts, 0, 1, 1, 1,
- 0, 0);
+ space_init(&space, params, &cosmo, dim, NULL, gparts, NULL, 0, nr_gparts, 0,
+ 1, 1, 0, 1, 0, 0);
struct engine engine;
engine.s = &space;
diff --git a/tests/testPotentialPair.c b/tests/testPotentialPair.c
index 5fbfb04ecd0b4a8c05c07d622dd77bd38af970ae..53fc54ccdd63a9a9150b6701c1a76ac20af91d4c 100644
--- a/tests/testPotentialPair.c
+++ b/tests/testPotentialPair.c
@@ -106,6 +106,7 @@ int main() {
props.a_smooth = 1.25;
props.r_cut_min = 0.;
props.theta_crit2 = 0.;
+ props.epsilon_cur = eps;
e.gravity_properties = &props;
struct runner r;
@@ -178,7 +179,6 @@ int main() {
gp->x[1] = 0.5;
gp->x[2] = 0.5;
gp->mass = 0.;
- gp->epsilon = eps;
gp->time_bin = 1;
gp->type = swift_type_dark_matter;
gp->id_or_neg_offset = n + 1;
@@ -196,7 +196,6 @@ int main() {
ci.gparts[0].x[1] = 0.5;
ci.gparts[0].x[2] = 0.5;
ci.gparts[0].mass = 1.;
- ci.gparts[0].epsilon = eps;
ci.gparts[0].time_bin = 1;
ci.gparts[0].type = swift_type_dark_matter;
ci.gparts[0].id_or_neg_offset = 1;
@@ -211,11 +210,12 @@ int main() {
for (int n = 0; n < num_tests; ++n) {
const struct gpart *gp = &cj.gparts[n];
const struct gpart *gp2 = &ci.gparts[0];
+ const double epsilon = gravity_get_softening(gp, &props);
double pot_true =
- potential(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr);
+ potential(ci.gparts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr);
double acc_true =
- acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr);
+ acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr);
message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - gp2->x[0],
gp->a_grav[0], acc_true, gp->potential, pot_true);
@@ -252,12 +252,12 @@ int main() {
for (int n = 0; n < num_tests; ++n) {
const struct gpart *gp = &cj.gparts[n];
const struct gravity_tensors *mpole = ci.multipole;
+ const double epsilon = gravity_get_softening(gp, &props);
double pot_true = potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0],
- gp->epsilon, rlr * FLT_MAX);
- double acc_true =
- acceleration(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], gp->epsilon,
- rlr * FLT_MAX);
+ epsilon, rlr * FLT_MAX);
+ double acc_true = acceleration(
+ mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr * FLT_MAX);
message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - mpole->CoM[0],
gp->a_grav[0], acc_true, gp->potential, pot_true);
@@ -304,7 +304,6 @@ int main() {
ci.gparts[n].mass = 1. / 8.;
- ci.gparts[n].epsilon = eps;
ci.gparts[n].time_bin = 1;
ci.gparts[n].type = swift_type_dark_matter;
ci.gparts[n].id_or_neg_offset = 1;
@@ -333,18 +332,19 @@ int main() {
for (int i = 0; i < 8; ++i) {
const struct gpart *gp2 = &ci.gparts[i];
+ const double epsilon = gravity_get_softening(gp, &props);
const double dx[3] = {gp2->x[0] - gp->x[0], gp2->x[1] - gp->x[1],
gp2->x[2] - gp->x[2]};
const double d = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]);
- pot_true += potential(gp2->mass, d, gp->epsilon, rlr * FLT_MAX);
+ pot_true += potential(gp2->mass, d, epsilon, rlr * FLT_MAX);
acc_true[0] -=
- acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[0] / d;
+ acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[0] / d;
acc_true[1] -=
- acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[1] / d;
+ acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[1] / d;
acc_true[2] -=
- acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[2] / d;
+ acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[2] / d;
}
message("x=%e f=%e f_true=%e pot=%e pot_true=%e %e %e",
diff --git a/tests/testPotentialSelf.c b/tests/testPotentialSelf.c
index 33fe574e34bcc49ff8bac6f19f22b55fda19f186..6d31f079fa79f7463637ec71dc2c75f37a10b129 100644
--- a/tests/testPotentialSelf.c
+++ b/tests/testPotentialSelf.c
@@ -41,15 +41,18 @@ const double eps = 0.02;
* @param s String used to identify this check in messages
*/
void check_value(double a, double b, const char *s) {
- if (fabs(a - b) / fabs(a + b) > 2e-6 && fabs(a - b) > 1.e-6)
+ if (fabs(a - b) / fabs(a + b) > 1e-6 && fabs(a - b) > 1.e-6)
error("Values are inconsistent: %12.15e %12.15e (%s)!", a, b, s);
}
/* Definitions of the potential and force that match
exactly the theory document */
-double S(double x) { return exp(x) / (1. + exp(x)); }
+double S(double x) { return good_approx_exp(x) / (1. + good_approx_exp(x)); }
-double S_prime(double x) { return exp(x) / ((1. + exp(x)) * (1. + exp(x))); }
+double S_prime(double x) {
+ return good_approx_exp(x) /
+ ((1. + good_approx_exp(x)) * (1. + good_approx_exp(x)));
+}
double potential(double mass, double r, double H, double rlr) {
@@ -103,6 +106,7 @@ int main() {
struct gravity_props props;
props.a_smooth = 1.25;
+ props.epsilon_cur = eps;
e.gravity_properties = &props;
struct runner r;
@@ -139,7 +143,6 @@ int main() {
c.gparts[0].x[1] = 0.5;
c.gparts[0].x[2] = 0.5;
c.gparts[0].mass = 1.;
- c.gparts[0].epsilon = eps;
c.gparts[0].time_bin = 1;
c.gparts[0].type = swift_type_dark_matter;
c.gparts[0].id_or_neg_offset = 1;
@@ -156,7 +159,6 @@ int main() {
gp->x[1] = 0.5;
gp->x[2] = 0.5;
gp->mass = 0.;
- gp->epsilon = eps;
gp->time_bin = 1;
gp->type = swift_type_dark_matter;
gp->id_or_neg_offset = n + 1;
@@ -172,9 +174,10 @@ int main() {
for (int n = 1; n < num_tests + 1; ++n) {
const struct gpart *gp = &c.gparts[n];
- double pot_true = potential(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
- double acc_true =
- acceleration(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr);
+ const double epsilon = gravity_get_softening(gp, &props);
+
+ double pot_true = potential(c.gparts[0].mass, gp->x[0], epsilon, rlr);
+ double acc_true = acceleration(c.gparts[0].mass, gp->x[0], epsilon, rlr);
// message("x=%e f=%e f_true=%e", gp->x[0], gp->a_grav[0], acc_true);
diff --git a/tests/testReading.c b/tests/testReading.c
index f5e2757c2fe3bd507e877f134e8c768aba0ae645..ca1e0ef69078c5e384a9cd4eab1098923ce9f279 100644
--- a/tests/testReading.c
+++ b/tests/testReading.c
@@ -48,7 +48,8 @@ int main() {
/* Read data */
read_ic_single("input.hdf5", &us, dim, &parts, &gparts, &sparts, &Ngas,
- &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, 1, 1, 0, 1, 0);
+ &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, 1, 1, 0, 0, 1.,
+ 1, 0);
/* Check global properties read are correct */
assert(dim[0] == boxSize);
diff --git a/tests/testSymmetry.c b/tests/testSymmetry.c
index 7a62f735b57b3a8fe5b91308b4dfe8ce00706e0e..6c6cd291dc4b919650d568354d7949f9ed4e0155 100644
--- a/tests/testSymmetry.c
+++ b/tests/testSymmetry.c
@@ -59,6 +59,8 @@ void test() {
pj.h = 2.f;
pi.id = 1ll;
pj.id = 2ll;
+ pi.time_bin = 1;
+ pj.time_bin = 1;
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
/* Give the primitive variables sensible values, since the Riemann solver does
@@ -157,8 +159,20 @@ void test() {
i_not_ok = memcmp(&pi, &pi2, sizeof(struct part));
j_not_ok = memcmp(&pj, &pj2, sizeof(struct part));
- if (i_not_ok) error("Particles 'pi' do not match after density");
- if (j_not_ok) error("Particles 'pj' do not match after density");
+ if (i_not_ok) {
+ printParticle_single(&pi, &xpi);
+ printParticle_single(&pi2, &xpi);
+ print_bytes(&pj, sizeof(struct part));
+ print_bytes(&pj2, sizeof(struct part));
+ error("Particles 'pi' do not match after force (byte = %d)", i_not_ok);
+ }
+ if (j_not_ok) {
+ printParticle_single(&pj, &xpj);
+ printParticle_single(&pj2, &xpj);
+ print_bytes(&pj, sizeof(struct part));
+ print_bytes(&pj2, sizeof(struct part));
+ error("Particles 'pj' do not match after force (byte = %d)", j_not_ok);
+ }
/* --- Test the force loop --- */
diff --git a/tests/tolerance_125_normal.dat b/tests/tolerance_125_normal.dat
index 0f11d03507b23c76b5703e118eede1359fe2afba..ff3d42e22b5b91378114294c4c81405492eafad2 100644
--- a/tests/tolerance_125_normal.dat
+++ b/tests/tolerance_125_normal.dat
@@ -1,4 +1,4 @@
# ID pos_x pos_y pos_z v_x v_y v_z h rho div_v S u P c a_x a_y a_z h_dt v_sig dS/dt du/dt
0 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4
0 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4 1e-4
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5 1e-5
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-5 2e-5 2e-5 1e-5 1e-5 1e-5 1e-5
diff --git a/tests/tolerance_27_normal.dat b/tests/tolerance_27_normal.dat
index 0fe55e84a42e7541068744e1e554afff1731ed3f..713e6c20c39ae9be0a4069a17c1ce45728fb5a34 100644
--- a/tests/tolerance_27_normal.dat
+++ b/tests/tolerance_27_normal.dat
@@ -1,4 +1,4 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 4e-5 4e-4 1e-2 1e-5 6e-6 6e-6 6e-6
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.2e-4 1e-4 2e-4 2e-4 1e-4 1e-4 1e-4
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 6.2e-5 4e-4 1.2e-2 1e-5 6e-6 6e-6 8e-6
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-4 1e-4 2e-4 2e-4 1e-4 1e-4 1e-4
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6
diff --git a/tests/tolerance_27_perturbed.dat b/tests/tolerance_27_perturbed.dat
index e1fd7518e55811386967938dfff0d62f04325cb0..4b7534f4ba43e7117452af819937f2d2b0451352 100644
--- a/tests/tolerance_27_perturbed.dat
+++ b/tests/tolerance_27_perturbed.dat
@@ -1,4 +1,4 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 1e-4 2e-4 1e-2 1e-5 3e-6 3e-6 7e-6
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.3e-3 1e-5 2e-3 6e-5 2e-3 2e-3 2e-3
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.3e-3 1e-5 2e-3 6e-5 3e-3 2e-3 2e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-3 1e-6 1e0 1e-6 2e-6 2e-6 2e-6
diff --git a/tests/tolerance_27_perturbed_h.dat b/tests/tolerance_27_perturbed_h.dat
index 5d6710c3946c0515f94c58eb3e0ecbfe14135b71..a1703f5b296d16f045fd28462fc9c61fc9358585 100644
--- a/tests/tolerance_27_perturbed_h.dat
+++ b/tests/tolerance_27_perturbed_h.dat
@@ -1,4 +1,4 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 3e-6 1e-4 5e-4 1.4e-2 1.1e-5 3e-6 3e-6 8e-6
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.5e-6 1.4e-2 1e-5 2e-3 2.5e-4 3e-3 3e-3 3e-3
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e0 1e-6 4e-6 4e-6 4e-6
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 3e-6 1e-4 5e-4 1.4e-2 1.1e-5 3e-6 3e-6 8e-6
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.5e-6 1.4e-2 1e-5 2e-3 2.5e-4 3e-3 3e-3 3e-3
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e0 1e-6 4e-6 4e-6 4e-6
diff --git a/tests/tolerance_27_perturbed_h2.dat b/tests/tolerance_27_perturbed_h2.dat
index 882554741734aeb270427b62580d6077907056ad..0b3921f0294f03c199e2045364144756dc5b066e 100644
--- a/tests/tolerance_27_perturbed_h2.dat
+++ b/tests/tolerance_27_perturbed_h2.dat
@@ -1,4 +1,4 @@
-# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 3e-6 1e-4 5e-4 1.5e-2 1.4e-5 3e-6 3e-6 1e-5
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.5e-6 2.5e-2 1e-5 5.86e-3 1.17e-3 3e-3 5.65e-3 3e-3
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e0 1e-6 4e-6 4e-6 4e-6
+# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 3e-6 1e-4 5e-4 1.5e-2 1.4e-5 3e-6 3e-6 1e-5
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.5e-6 2.5e-2 1e-5 5.86e-3 1.17e-3 3e-3 8e-3 3e-3
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e0 1e-6 4e-6 4e-6 4e-6
diff --git a/tests/tolerance_pair_force_active.dat b/tests/tolerance_pair_force_active.dat
index c5b5620338b848f1b8f4658049e53a6424e93a60..f4394edf35e7287968673dfb949c726462b5d8ad 100644
--- a/tests/tolerance_pair_force_active.dat
+++ b/tests/tolerance_pair_force_active.dat
@@ -1,4 +1,4 @@
# ID wcount h_dt
0 1 2.4e-3
0 1 2.4e-3
- 0 1 1e-5
+ 0 1 3e-5
diff --git a/tests/tolerance_periodic_BC_normal.dat b/tests/tolerance_periodic_BC_normal.dat
index 823e4af488b343f57e3c90e89ee2d4f13d3ca94b..fdbcd53a1f0233a87d657e5bc3751e685631a4c7 100644
--- a/tests/tolerance_periodic_BC_normal.dat
+++ b/tests/tolerance_periodic_BC_normal.dat
@@ -1,4 +1,4 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-6 4e-5 1e-3 1e-2 2e-4 2e-4 2e-4 2e-4
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 4e-6 7e-5 1e-3 1.3e-2 2e-4 2e-4 2e-4 2e-4
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 2e-4 1e-4 2e-4 6e-4 2e-3 2e-3 2e-3
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-4 1e-6 1e-4 5e-4 2e-4 2e-4 2e-4
diff --git a/theory/Multipoles/potential_derivatives.tex b/theory/Multipoles/potential_derivatives.tex
index 9dc9244ca56eaa56047f561aa221dd6a761e7a5a..c71cc849d5c5e0f7fbddb6b43e745838addd08cb 100644
--- a/theory/Multipoles/potential_derivatives.tex
+++ b/theory/Multipoles/potential_derivatives.tex
@@ -18,8 +18,10 @@ by constructing derivatives of the truncated potentials:
\chi^{(4)}(r, r_s) &= \frac{16}{r_s^4} \left(48\alpha(x)^5 - 120\alpha(x)^4 + 100\alpha(x)^3 -30 \alpha(x)^2 + 2\alpha(x)\right) \nonumber \\
\chi^{(5)}(r, r_s) &= \frac{32}{r_s^5} \left(240\alpha(x)^6 - 720\alpha(x)^5 + 780\alpha(x)^4 - 360\alpha(x)^3 + 62\alpha(x)^2 - 2\alpha(x) \right) \nonumber
\end{align}
-We can now construct common quantities that appear in derivatives of
-multiple orders:
+In the Newtonian limit ($r_s\rightarrow\infty$) the first expression
+reduces to $\chi(r,r_s) = 1$ whilst all higher-order derivatives
+vanish. We can now construct common quantities that appear in
+derivatives of multiple orders:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{align}