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/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/main.c b/examples/main.c
index 80e4dfb479172d9e2bca6fb19ee676371d130afe..f490465822859bc991eb31b9a857753c85934d6a 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 */
@@ -635,20 +647,20 @@ int main(int argc, char *argv[]) {
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
if (myrank == 0) {
clocks_gettime(&toc);
@@ -667,6 +679,9 @@ int main(int argc, char *argv[]) {
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 */
+ part_verify_links(parts, gparts, sparts, Ngas, Ngpart, Nspart, 1);
#endif
/* Get the total number of particles across all nodes. */
@@ -690,8 +705,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);
@@ -836,7 +852,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);
@@ -1000,6 +1016,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 %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(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/src/cell.c b/src/cell.c
index b415cc077b569c84c85462634458e68aa5e4833a..ea036c97081e625c4e4ca2bb1d06cd9ea8251bd3 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1081,23 +1081,6 @@ void cell_sanitize(struct cell *c, int treated) {
c->h_max = h_max;
}
-/**
- * @brief Converts hydro quantities to a valid state after the initial density
- * calculation
- *
- * @param c Cell to act upon
- * @param data Unused parameter
- */
-void cell_convert_hydro(struct cell *c, void *data) {
-
- struct part *p = c->parts;
- struct xpart *xp = c->xparts;
-
- for (int i = 0; i < c->count; ++i) {
- hydro_convert_quantities(&p[i], &xp[i]);
- }
-}
-
/**
* @brief Cleans the links in a given cell.
*
diff --git a/src/cell.h b/src/cell.h
index 91663d45bf048fff6e4bbf02d404ec62d216365b..87af742baf935668b0f065e70f065b978b6c9b9f 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -492,7 +492,6 @@ int cell_getsize(struct cell *c);
int cell_link_parts(struct cell *c, struct part *parts);
int cell_link_gparts(struct cell *c, struct gpart *gparts);
int cell_link_sparts(struct cell *c, struct spart *sparts);
-void cell_convert_hydro(struct cell *c, void *data);
void cell_clean_links(struct cell *c, void *data);
void cell_make_multipoles(struct cell *c, integertime_t ti_current);
void cell_check_multipole(struct cell *c, void *data);
diff --git a/src/common_io.c b/src/common_io.c
index 67326817397118568e060f5dff49421bb32fba4d..19adc57fbd51972336840669dab473f846d23dee 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -457,6 +457,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 +466,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 +479,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 +488,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 +708,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 +764,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/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/engine.c b/src/engine.c
index 18a3d4d39d390e435d1a6a5c4ee8b9fbdbc1506b..b66ceb4037f4de87252ce3203d0f7c2933d9ffb9 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3665,10 +3665,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();
@@ -3679,7 +3679,7 @@ void engine_rebuild(struct engine *e, int clean_h_values) {
space_rebuild(e->s, e->verbose);
/* 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. */
@@ -4395,7 +4395,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 {
@@ -4404,6 +4404,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);
@@ -5193,7 +5197,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/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..2508f69a636b93c218cdcb0770cf00ef969389c2 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
*
@@ -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_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_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_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/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/part.c b/src/part.c
index 712a13d42a629355290bcf89f31b94fb9670a215..c43ffa4820504282b4fb02e63b5cfc4196f3be77 100644
--- a/src/part.c
+++ b/src/part.c
@@ -119,7 +119,7 @@ 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)
+ if (gparts[k].id_or_neg_offset <= 0)
error("DM gpart particle linked to something !");
}
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/scheduler.c b/src/scheduler.c
index 072235be1807db972208514739409cdb1e76c0d6..151304293749a29abe9bd9680d8f8f81bc845884 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1435,7 +1435,8 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
switch (t->type) {
case task_type_self:
case task_type_sub_self:
- if (t->subtype == task_subtype_grav || t->subtype == task_subtype_external_grav)
+ if (t->subtype == task_subtype_grav ||
+ t->subtype == task_subtype_external_grav)
qid = t->ci->super_gravity->owner;
else
qid = t->ci->super_hydro->owner;
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 329526fe14e38220750903058bc5201ef7486f63..273e01858f443174d5fb90477cc38e19d4375d67 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2647,8 +2647,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;
@@ -2661,6 +2673,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);
@@ -2684,9 +2700,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;
@@ -2715,9 +2738,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;
@@ -2792,11 +2822,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);
}
/**
@@ -2824,6 +2855,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.
@@ -2832,8 +2864,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.
*
@@ -2843,9 +2876,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 */
@@ -2856,7 +2890,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;
@@ -2868,6 +2902,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",
@@ -2880,6 +2927,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 */
@@ -2926,8 +2977,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;
@@ -3144,6 +3195,97 @@ 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;
+
+ /* 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 = d;
+ }
+
+ /* Replace the content of the space */
+ free(s->gparts);
+ s->parts = parts;
+ s->gparts = gparts;
+}
+
/**
* @brief Cleans-up all the cell links in the space
*
diff --git a/src/space.h b/src/space.h
index 11cbaabdc9fc3ae17024c042ae868d464954d501..76ff2ea9908195d618663eb1874eb48bec463bb7 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
@@ -182,9 +183,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,
@@ -239,6 +241,8 @@ 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_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..246139ba243a37dbc5f6cafbc488cb1bf6cdd8eb 100644
--- a/tests/testPotentialSelf.c
+++ b/tests/testPotentialSelf.c
@@ -103,6 +103,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 +140,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 +156,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 +171,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);