diff --git a/examples/CoolingBox/coolingBox_solar.yml b/examples/CoolingBox/coolingBox_solar.yml
index 1dd5ee6a3cbf9efea6609a8c2cdc3dff62bd3867..ccdb85eb802d6ee8e1368529839e058a47d165a3 100644
--- a/examples/CoolingBox/coolingBox_solar.yml
+++ b/examples/CoolingBox/coolingBox_solar.yml
@@ -9,8 +9,8 @@ InternalUnitSystem:
# Cosmological parameters
Cosmology:
h: 0.6777 # Reduced Hubble constant
- a_begin: 0.142857142857 # Initial scale-factor of the simulation
- a_end: 0.144285714286 # Final scale factor of the simulation
+ a_begin: 0.1 # Initial scale-factor of the simulation
+ a_end: 1.0 # Final scale factor of the simulation
Omega_m: 0.307 # Matter density parameter
Omega_lambda: 0.693 # Dark-energy density parameter
Omega_b: 0.0455 # Baryon density parameter
@@ -20,7 +20,7 @@ TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 1e-2 # The end time of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
- dt_max: 1e-6 # The maximal time-step size of the simulation (in internal units).
+ dt_max: 1e-7 # The maximal time-step size of the simulation (in internal units).
Scheduler:
max_top_level_cells: 15
@@ -55,6 +55,7 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: ./coolingBox.hdf5 # The file to read
+ periodic: 1
# Dimensionless pre-factor for the time-step condition
LambdaCooling:
diff --git a/examples/CoolingBox/testCooling.c b/examples/CoolingBox/testCooling.c
index fca5429bd1437514c43b7eb705533ccd4fb00ac9..148c36cee0fb185c9a770245b62021974b85f8e7 100644
--- a/examples/CoolingBox/testCooling.c
+++ b/examples/CoolingBox/testCooling.c
@@ -40,29 +40,25 @@ void set_quantities(struct part *restrict p,
struct xpart *restrict xp,
const struct unit_system *restrict us,
const struct cooling_function_data *restrict cooling,
- const struct cosmology *restrict cosmo,
+ struct cosmology *restrict cosmo,
const struct phys_const *restrict phys_const, float nh_cgs,
- double u) {
+ double u,
+ integertime_t ti_current) {
- float scale_factor = 1.0 / (1.0 + cosmo->z);
- double hydrogen_number_density =
- nh_cgs * units_cgs_conversion_factor(us, UNIT_CONV_LENGTH) *
- units_cgs_conversion_factor(us, UNIT_CONV_LENGTH) *
- units_cgs_conversion_factor(us, UNIT_CONV_LENGTH);
+ /* Update cosmology quantities */
+ cosmology_update(cosmo,phys_const,ti_current);
+
+ /* calculate density */
+ double hydrogen_number_density = nh_cgs / cooling->number_density_scale;
p->rho = hydrogen_number_density * phys_const->const_proton_mass /
p->chemistry_data.metal_mass_fraction[chemistry_element_H] * (cosmo->a * cosmo->a * cosmo->a);
- float pressure = (u * scale_factor * scale_factor) / cooling->internal_energy_scale *
+ /* update entropy based on internal energy */
+ float pressure = (u * cosmo->a * cosmo->a) / cooling->internal_energy_scale *
p->rho * (hydro_gamma_minus_one);
p->entropy = pressure * (pow(p->rho, -hydro_gamma));
xp->entropy_full = p->entropy;
- // Using hydro_set_init_internal_energy seems to work better for higher z for
- // setting the internal energy correctly However, with Gadget2 this just sets
- // the entropy to the internal energy, which needs to be converted somehow
- if (cosmo->z >= 1)
- hydro_set_init_internal_energy(
- p, (u * scale_factor * scale_factor) / cooling->internal_energy_scale);
}
/*
@@ -81,132 +77,106 @@ int main(int argc, char **argv) {
struct phys_const phys_const;
struct cooling_function_data cooling;
struct cosmology cosmo;
- char *parametersFileName = "./coolingBox.yml";
+ char *parametersFileName = "./coolingBox_solar.yml";
float nh; // hydrogen number density
double u; // internal energy
- const int npts = 250; // number of values for the internal energy at which
- // cooling rate is evaluated
-
- // Read options
- int param;
- float redshift = -1.0,
- log_10_nh =
- 100; // unreasonable values will be changed if not set in options
- while ((param = getopt(argc, argv, "z:d:m:t")) != -1) switch (param) {
- case 'z':
- // read redshift
- redshift = atof(optarg);
- break;
- case 'd':
- // read log10 of hydrogen number density
- log_10_nh = atof(optarg);
- break;
- case 'm':
- // read which yml file we need to use
- parametersFileName = optarg;
- break;
- case '?':
- if (optopt == 'z')
- printf("Option -%c requires an argument.\n", optopt);
- else
- printf("Unknown option character `\\x%x'.\n", optopt);
- error("invalid option(s) to testCooling");
- }
- // Read the parameter file
+ /* Number of values to test for in redshift,
+ * hydrogen number density and internal energy */
+ const int n_z = 50;
+ const int n_nh = 50;
+ const int n_u = 50;
+
+ /* Number of subcycles and tolerance used to compare
+ * subcycled and implicit solution. Note, high value
+ * of tolerance due to mismatch between explicit and
+ * implicit solutioin for large timesteps */
+ const int n_subcycle = 1000;
+ const float integration_tolerance = 0.2;
+
+ /* Set dt */
+ const float dt_cool = 1.0e-5;
+ const float dt_therm = 1.0e-5;
+
+ /* Read the parameter file */
if (params == NULL) error("Error allocating memory for the parameter file.");
message("Reading runtime parameters from file '%s'", parametersFileName);
parser_read_file(parametersFileName, params);
- // Init units
+ /* Init units */
units_init_from_params(&us, params, "InternalUnitSystem");
phys_const_init(&us, params, &phys_const);
- // Init chemistry
+ /* Init chemistry */
chemistry_init(params, &us, &phys_const, &chem_data);
chemistry_first_init_part(&phys_const, &us, &cosmo, &chem_data, &p, &xp);
chemistry_print(&chem_data);
- // Init cosmology
+ /* Init cosmology */
cosmology_init(params, &us, &phys_const, &cosmo);
cosmology_print(&cosmo);
- if (redshift == -1.0) {
- cosmo.z = 7.0;
- } else {
- cosmo.z = redshift;
- }
- message("redshift %.5e", cosmo.z);
- // Init cooling
+ /* Init cooling */
cooling_init(params, &us, &phys_const, &cooling);
cooling_print(&cooling);
cooling_update(&cosmo, &cooling, 0);
- // Calculate abundance ratios
+ /* Calculate abundance ratios */
float *abundance_ratio;
abundance_ratio = malloc((chemistry_element_count + 2) * sizeof(float));
abundance_ratio_to_solar(&p, &cooling, abundance_ratio);
- // extract mass fractions, calculate table indices and offsets
+ /* extract mass fractions, calculate table indices and offsets */
float XH = p.chemistry_data.metal_mass_fraction[chemistry_element_H];
float HeFrac =
p.chemistry_data.metal_mass_fraction[chemistry_element_He] /
(XH + p.chemistry_data.metal_mass_fraction[chemistry_element_He]);
- int He_i, n_h_i;
- float d_He, d_n_h;
+ int He_i;
+ float d_He;
get_index_1d(cooling.HeFrac, cooling.N_He, HeFrac, &He_i, &d_He);
- // Calculate contributions from metals to cooling rate
- // open file
- const char *output_filename = "cooling_output.dat";
- FILE *output_file = fopen(output_filename, "w");
- if (output_file == NULL) {
- printf("Error opening file!\n");
- exit(1);
- }
-
- // set hydrogen number density
- if (log_10_nh == 100) {
- // hydrogen number density not specified in options
- nh = 1.0e-1;
- } else {
- nh = exp(M_LN10 * log_10_nh);
- }
- nh = 5.6e-2;
-
- // set internal energy to dummy value, will get reset when looping over
- // internal energies
- u = 1.0e14;
- set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh, u);
- float inn_h = hydro_get_physical_density(&p, &cosmo) * XH / phys_const.const_proton_mass
- *cooling.number_density_scale;
- get_index_1d(cooling.nH, cooling.N_nH, log10(inn_h), &n_h_i, &d_n_h);
- message("inn_h %.5e nh %.5e", inn_h, nh);
-
- // Loop over internal energy
- float du;
- double temperature;
- for (int j = 0; j < npts; j++) {
- set_quantities(&p, &xp, &us, &cooling, &cosmo, &phys_const, nh,
- exp(M_LN10 * (10.0 + j * 8.0 / npts)));
- u = hydro_get_physical_internal_energy(&p, &xp, &cosmo) *
- cooling.internal_energy_scale;
- float cooling_du_dt;
-
- // calculate cooling rates
- double dLambdaNet_du;
- temperature = eagle_convert_u_to_temp(log10(u), &du, n_h_i, He_i, d_n_h, d_He, &cooling, &cosmo);
- //cooling_du_dt = eagle_print_metal_cooling_rate(
- // n_h_i, d_n_h, He_i, d_He, &p, &xp, &cooling, &cosmo, &phys_const,
- // abundance_ratio);
- cooling_du_dt = eagle_cooling_rate(
- log(u), &dLambdaNet_du, n_h_i, d_n_h, He_i, d_He, &p, &cooling, &cosmo, &phys_const,
- abundance_ratio);
- fprintf(output_file, "%.5e %.5e\n", u, cooling_du_dt);
+ /* Cooling function needs to know the minimal energy. Set it to the lowest
+ * internal energy in the cooling table. */
+ struct hydro_props hydro_properties;
+ hydro_properties.minimal_internal_energy = exp(M_LN10 * cooling.Therm[0]) / cooling.internal_energy_scale;
+
+ /* calculate spacing in nh and u */
+ const float delta_nh = (cooling.nH[cooling.N_nH - 1] - cooling.nH[0])/n_nh;
+ const float delta_u = (cooling.Therm[cooling.N_Temp - 1] - cooling.Therm[0])/n_u;
+
+ for(int z_i = 0; z_i < n_z; z_i++) {
+ integertime_t ti_current = max_nr_timesteps/n_z*z_i;
+ for(int nh_i = 0; nh_i < n_nh; nh_i++) {
+ nh = exp(M_LN10 * cooling.nH[0] + delta_nh*nh_i);
+ for (int u_i = 0; u_i < n_u; u_i++) {
+ message("z_i %d n_h_i %d u_i %d", z_i, nh_i, u_i);
+ u = exp(M_LN10 * cooling.Therm[0] + delta_u*u_i);
+
+ /* update nh, u, z */
+ set_quantities(&p,&xp,&us,&cooling,&cosmo,&phys_const, nh, u, ti_current);
+
+ /* calculate subcycled solution */
+ for (int t_subcycle = 0; t_subcycle < n_subcycle; t_subcycle++) {
+ cooling_cool_part(&phys_const,&us,&cosmo,&hydro_properties,&cooling,&p,&xp,dt_cool/n_subcycle,dt_therm/n_subcycle);
+ xp.entropy_full += p.entropy_dt*dt_therm/n_subcycle;
+ }
+ double u_subcycled = hydro_get_physical_internal_energy(&p,&xp,&cosmo)*cooling.internal_energy_scale;
+
+ /* reset quantities to nh, u, and z that we want to test */
+ set_quantities(&p,&xp,&us,&cooling,&cosmo,&phys_const, nh, u, ti_current);
+
+ /* compute implicit solution */
+ cooling_cool_part(&phys_const,&us,&cosmo,&hydro_properties,&cooling,&p,&xp,dt_cool,dt_therm);
+ double u_implicit = hydro_get_physical_internal_energy(&p,&xp,&cosmo)*cooling.internal_energy_scale;
+
+ /* check if the two solutions are consistent */
+ if (fabs((u_implicit-u_subcycled)/u_subcycled) > integration_tolerance) message("implicit and subcycled solutions do not match. z_i %d nh_i %d u_i %d implicit %.5e subcycled %.5e error %.5e", z_i, nh_i, u_i, u_implicit, u_subcycled, fabs((u_implicit-u_subcycled)/u_subcycled));
+ }
+ }
}
- fclose(output_file);
- message("done cooling rates test");
+
+ message("done test");
free(params);
return 0;
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index d6de3bf7788817c12875f5ec359cb0cb89d2038e..72a5efeed28e03622279dcbc83c691efe334e579 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -132,6 +132,7 @@ void cooling_cool_part(const struct phys_const *restrict phys_const,
/* Get internal energy at the end of the next kick step (assuming dt does not increase) */
double u_0 = (u_start + hydro_du_dt * dt_therm);
+ // for unit test
u_0 = u_start;
/* Check for minimal energy */
@@ -147,9 +148,9 @@ void cooling_cool_part(const struct phys_const *restrict phys_const,
abundance_ratio_to_solar(p, cooling, abundance_ratio);
/* Get the H and He mass fractions */
- const float XH = p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
- const float HeFrac = p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_He] /
- (XH + p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_He]);
+ const float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
+ const float HeFrac = p->chemistry_data.metal_mass_fraction[chemistry_element_He] /
+ (XH + p->chemistry_data.metal_mass_fraction[chemistry_element_He]);
/* convert Hydrogen mass fraction in Hydrogen number density */
const double n_h = hydro_get_physical_density(p, cosmo) * XH / phys_const->const_proton_mass
@@ -213,13 +214,12 @@ void cooling_cool_part(const struct phys_const *restrict phys_const,
abundance_ratio, dt_cgs);
if (log_u_final_cgs == -100) {
message("particle %llu failed to converge with bisection method, assuming no cooling.", p->id);
- log_u_final_cgs = log(u_0_cgs);
+ log_u_final_cgs = log(u_0_cgs);
}
}
u_final_cgs = exp(log_u_final_cgs);
}
- if (p->id == 1) message("id %llu u_final %.5e u_0 %.5e ratefact %.5e Lambda %.5e dt %.5e rho %.5e nh %.5e rho internal comoving %.5e", p->id, u_final_cgs, u_0_cgs, ratefact, LambdaNet, dt_cgs, hydro_get_physical_density(p, cosmo)*units_cgs_conversion_factor(us,UNIT_CONV_DENSITY), n_h, p->rho);
/* Expected change in energy over the next kick step (assuming no change in dt) */
const double delta_u_cgs = u_final_cgs - u_start_cgs;
@@ -241,8 +241,9 @@ void cooling_cool_part(const struct phys_const *restrict phys_const,
* that could potentially be for a time-step twice as big. We hence check
* for 2.5 delta_u but this time against 0 energy not the minimum.
* To avoid numerical rounding bringing us below 0., we add a tiny tolerance. */
- if (u_start + 2.5 * delta_u < 0.)
+ if (u_start + 2.5 * delta_u < 0.) {
delta_u = -u_start / (2.5 + rounding_tolerance);
+ }
/* Turn this into a rate of change */
const float cooling_du_dt = delta_u / dt_therm ;
@@ -264,16 +265,18 @@ __attribute__((always_inline)) INLINE double
eagle_helium_reionization_extraheat(
double z, double dz, const struct cooling_function_data *restrict cooling) {
- //double he_reion_erg_pG = cooling->he_reion_ev_pH / cooling->proton_mass_cgs;
- //double extra_heating = he_reion_erg_pG *
- // (erf((z - dz - cooling->he_reion_z_center) /
- // (M_SQRT2 * cooling->he_reion_z_sigma)) -
- // erf((z - cooling->he_reion_z_center) /
- // (M_SQRT2 * cooling->he_reion_z_sigma))) /
- // 2.0;
+ double he_reion_erg_pG = cooling->he_reion_ev_pH / cooling->proton_mass_cgs;
+ double extra_heating = he_reion_erg_pG *
+ (erf((z - dz - cooling->he_reion_z_center) /
+ (M_SQRT2 * cooling->he_reion_z_sigma)) -
+ erf((z - cooling->he_reion_z_center) /
+ (M_SQRT2 * cooling->he_reion_z_sigma))) /
+ 2.0;
+
+ return extra_heating;
- //return extra_heating;
- return 0.0;
+ // for unit test
+ //return 0.0;
}
/*
@@ -320,7 +323,7 @@ __attribute__((always_inline)) INLINE double eagle_metal_cooling_rate(
double solar_electron_abundance;
/* convert Hydrogen mass fraction in Hydrogen number density */
- const float XH = p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
+ const float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
const double n_h = hydro_get_physical_density(p, cosmo) * XH / phys_const->const_proton_mass
*cooling->number_density_scale;
@@ -677,10 +680,10 @@ __attribute__((always_inline)) INLINE void abundance_ratio_to_solar(
if (elem == chemistry_element_Fe) {
/* NOTE: solar abundances have iron last with calcium and sulphur directly
* before, hence +2 */
- ratio_solar[elem] = p->chemistry_data.smoothed_metal_mass_fraction[elem] /
+ ratio_solar[elem] = p->chemistry_data.metal_mass_fraction[elem] /
cooling->SolarAbundances[elem + 2];
} else {
- ratio_solar[elem] = p->chemistry_data.smoothed_metal_mass_fraction[elem] /
+ ratio_solar[elem] = p->chemistry_data.metal_mass_fraction[elem] /
cooling->SolarAbundances[elem];
}
}
@@ -688,11 +691,11 @@ __attribute__((always_inline)) INLINE void abundance_ratio_to_solar(
/* assign ratios for Ca and S, note positions of these elements occur before
* Fe */
ratio_solar[chemistry_element_count] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Si] *
+ p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
cooling->sulphur_over_silicon_ratio /
cooling->SolarAbundances[chemistry_element_count - 1];
ratio_solar[chemistry_element_count + 1] =
- p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Si] *
+ p->chemistry_data.metal_mass_fraction[chemistry_element_Si] *
cooling->calcium_over_silicon_ratio /
cooling->SolarAbundances[chemistry_element_count];
}
@@ -738,7 +741,7 @@ __attribute__((always_inline)) INLINE float newton_iter(
(cooling->Therm[0] + 0.05) / M_LOG10E;
/* convert Hydrogen mass fraction in Hydrogen number density */
- const float XH = p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
+ const float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
const double n_h = hydro_get_physical_density(p, cosmo) * XH / phys_const->const_proton_mass
*cooling->number_density_scale;
@@ -820,15 +823,29 @@ __attribute__((always_inline)) INLINE float bisection_iter(
int i = 0;
double u_init = exp(logu_init);
+ FILE *output_file = fopen("bisection_output.dat","w");
+ FILE *bracketing = fopen("bracketing_output.dat","w");
+ FILE *bisection_function = fopen("bisection_function.dat","w");
+
/* convert Hydrogen mass fraction in Hydrogen number density */
- const float XH = p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H];
+ const float XH = p->chemistry_data.metal_mass_fraction[chemistry_element_H];
const double n_h = hydro_get_physical_density(p, cosmo) * XH / phys_const->const_proton_mass
*cooling->number_density_scale;
-
+
/* compute ratefact = n_h * n_h / rho; Might lead to round-off error:
* replaced by equivalent expression below */
const double ratefact = n_h * (XH / cooling->proton_mass_cgs);
+ // Debugging
+ for(int j = 0; j < cooling->N_Temp; j++) {
+ float u_print = exp(M_LN10*cooling->Therm[j]);
+ LambdaNet = (He_reion_heat / (dt * ratefact)) +
+ eagle_cooling_rate(log(u_print), dLambdaNet_du, n_h_i, d_n_h,
+ He_i, d_He, p, cooling, cosmo, phys_const,
+ abundance_ratio);
+ fprintf(bisection_function, "%.5e %.5e\n", u_print, u_print - u_ini - LambdaNet * ratefact * dt);
+ }
+ fclose(bisection_function);
/* Bracketing */
u_lower = u_init;
u_upper = u_init;
@@ -842,6 +859,7 @@ __attribute__((always_inline)) INLINE float bisection_iter(
/* we're cooling */
u_lower /= bracket_factor;
u_upper *= bracket_factor;
+
LambdaNet = (He_reion_heat / (dt * ratefact)) +
eagle_cooling_rate(log(u_lower), dLambdaNet_du, n_h_i, d_n_h,
@@ -854,6 +872,7 @@ __attribute__((always_inline)) INLINE float bisection_iter(
eagle_cooling_rate(log(u_lower), dLambdaNet_du, n_h_i, d_n_h,
He_i, d_He, p, cooling, cosmo, phys_const,
abundance_ratio);
+ fprintf(bracketing, "%.5e %.5e %.5e %.5e\n", u_lower, u_upper, LambdaNet, u_lower - u_ini - LambdaNet * ratefact * dt);
i++;
}
if (i >= bisection_max_iterations) {
@@ -876,6 +895,7 @@ __attribute__((always_inline)) INLINE float bisection_iter(
eagle_cooling_rate(log(u_upper), dLambdaNet_du, n_h_i, d_n_h,
He_i, d_He, p, cooling, cosmo, phys_const,
abundance_ratio);
+ fprintf(bracketing, "%.5e %.5e %.5e %.5e\n", u_lower, u_upper, LambdaNet, u_upper - u_ini - LambdaNet * ratefact * dt);
i++;
}
if (i >= bisection_max_iterations) {
@@ -883,6 +903,7 @@ __attribute__((always_inline)) INLINE float bisection_iter(
return -100;
}
}
+ fclose(bracketing);
/* bisection iteration */
i = 0;
@@ -897,6 +918,7 @@ __attribute__((always_inline)) INLINE float bisection_iter(
} else {
u_lower = u_next;
}
+ fprintf(output_file,"%.5e\n", u_next);
i++;
} while (fabs(u_upper - u_lower) / u_next > bisection_tolerance &&
@@ -906,6 +928,7 @@ __attribute__((always_inline)) INLINE float bisection_iter(
return -100;
message("Particle id %llu failed to converge", p->id); // WARNING: Do we want this to throw an error? In EAGLE calculation continued.
}
+ fclose(output_file);
return log(u_upper);
}
@@ -1032,6 +1055,7 @@ void cooling_init_backend(struct swift_params *parameter_file,
const double compton_coefficient_cgs =
compton_coefficient *
units_general_cgs_conversion_factor(us, dimension_coefficient);
+
#ifdef SWIFT_DEBUG_CHECKS
const double expected_compton_coefficient_cgs = 1.0178085e-37;
if (fabs(compton_coefficient_cgs - expected_compton_coefficient_cgs) /
@@ -1064,10 +1088,16 @@ void cooling_init_backend(struct swift_params *parameter_file,
void cooling_restore_tables(struct cooling_function_data* cooling,
const struct cosmology* cosmo){
+ /* Read redshifts */
get_cooling_redshifts(cooling);
+
+ /* Read cooling header */
char fname[eagle_table_path_name_length + 12];
sprintf(fname, "%sz_0.000.hdf5", cooling->cooling_table_path);
read_cooling_header(fname, cooling);
+
+ /* Read relevant cooling tables.
+ * Third variable in cooling_update flag to mark restart*/
allocate_cooling_tables(cooling);
cooling_update(cosmo, cooling, 1);
}