diff --git a/configure.ac b/configure.ac index c0ce630351f471de337d73a0f401227a91a38dc2..1130f5d875ccbaa6450371640f0156cefd9a56ca 100644 --- a/configure.ac +++ b/configure.ac @@ -1371,7 +1371,7 @@ esac # Hydro scheme. AC_ARG_WITH([hydro], [AS_HELP_STRING([--with-hydro=], - [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, anarchy-pu debug default: gadget2@:>@] + [Hydro dynamics to use @<:@gadget2, minimal, pressure-entropy, pressure-energy, pressure-energy-monaghan, default, gizmo-mfv, gizmo-mfm, shadowfax, planetary, anarchy-pu default: gadget2@:>@] )], [with_hydro="$withval"], [with_hydro="gadget2"] diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in index 8e5fa52b236d91cabf51ce0c4f481a7334fd315d..a3da1b4f8edc8c1af225ac941441a92b14a13f99 100644 --- a/doc/Doxyfile.in +++ b/doc/Doxyfile.in @@ -811,7 +811,7 @@ RECURSIVE = NO # Note that relative paths are relative to the directory from which doxygen is # run. -EXCLUDE = +EXCLUDE = @top_srcdir@/src/cooling/EAGLE/newton_cooling.c # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or # directories that are symbolic links (a Unix file system feature) are excluded diff --git a/examples/Cooling/CoolingRates/cooling_rates.c b/examples/Cooling/CoolingRates/cooling_rates.c index 5ffdc7f8d56bfea68119d275b8bb574658be9d91..07d82b3056223b0a36498c4815e08db7c764b911 100644 --- a/examples/Cooling/CoolingRates/cooling_rates.c +++ b/examples/Cooling/CoolingRates/cooling_rates.c @@ -104,9 +104,9 @@ INLINE static double eagle_print_metal_cooling_rate( /* calculate cooling rates */ for (int j = 0; j < eagle_cooling_N_metal + 2; j++) element_lambda[j] = 0.0; - lambda_net = eagle_metal_cooling_rate( - log10(u), cosmo->z, n_h, abundance_ratio, n_h_i, d_n_h, He_i, d_He, - cooling, /*dLambdaNet_du=*/NULL, element_lambda); + lambda_net = + eagle_metal_cooling_rate(log10(u), cosmo->z, n_h, abundance_ratio, n_h_i, + d_n_h, He_i, d_He, cooling, element_lambda); /* write cooling rate contributions to their own files. */ for (int j = 0; j < eagle_cooling_N_metal + 2; j++) { @@ -235,18 +235,26 @@ int main(int argc, char **argv) { // Init cooling cooling_init(params, &us, &internal_const, &hydro_properties, &cooling); + cooling.H_reion_done = 1; cooling_print(&cooling); cooling_update(&cosmo, &cooling, &s); + // Copy over the raw metals into the smoothed metals + memcpy(&p.chemistry_data.smoothed_metal_mass_fraction, + &p.chemistry_data.metal_mass_fraction, + chemistry_element_count * sizeof(float)); + p.chemistry_data.smoothed_metal_mass_fraction_total = + p.chemistry_data.metal_mass_fraction_total; + // Calculate abundance ratios float abundance_ratio[(chemistry_element_count + 2)]; abundance_ratio_to_solar(&p, &cooling, abundance_ratio); // 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]); + float XH = p.chemistry_data.smoothed_metal_mass_fraction[chemistry_element_H]; + float XHe = + p.chemistry_data.smoothed_metal_mass_fraction[chemistry_element_He]; + float HeFrac = XHe / (XH + XHe); int He_i, n_h_i; float d_He, d_n_h; get_index_1d(cooling.HeFrac, eagle_cooling_N_He_frac, HeFrac, &He_i, &d_He); @@ -286,7 +294,7 @@ int main(int argc, char **argv) { // calculate cooling rates const double temperature = eagle_convert_u_to_temp( - log10(u), cosmo.z, 0, NULL, n_h_i, He_i, d_n_h, d_He, &cooling); + log10(u), cosmo.z, n_h_i, He_i, d_n_h, d_He, &cooling); const double cooling_du_dt = eagle_print_metal_cooling_rate( n_h_i, d_n_h, He_i, d_He, &p, &xp, &cooling, &cosmo, &internal_const, diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml index 13ed73554936abc824ee990322cfd49464897c5b..821dca7180275f813fcc8a1271c9062a69c9ff87 100644 --- a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_high_z.yml @@ -44,7 +44,20 @@ InitialConditions: periodic: 1 # Parameters for the EAGLE chemistry -EAGLEChemistry: # Solar abundances +# EAGLEChemistry: # Solar abundances +# init_abundance_metal: 0.014 +# init_abundance_Hydrogen: 0.70649785 +# init_abundance_Helium: 0.28055534 +# init_abundance_Carbon: 2.0665436e-3 +# init_abundance_Nitrogen: 8.3562563e-4 +# init_abundance_Oxygen: 5.4926244e-3 +# init_abundance_Neon: 1.4144605e-3 +# init_abundance_Magnesium: 5.907064e-4 +# init_abundance_Silicon: 6.825874e-4 +# init_abundance_Iron: 1.1032152e-3 + +# Parameters for the EAGLE chemistry +EAGLEChemistry: # Primoridal abundances init_abundance_metal: 0. init_abundance_Hydrogen: 0.752 init_abundance_Helium: 0.248 @@ -82,4 +95,4 @@ LambdaCooling: ConstCooling: cooling_rate: 1. # Cooling rate (du/dt) (internal units) min_energy: 1. # Minimal internal energy per unit mass (internal units) - cooling_tstep_mult: 1. # Dimensionless pre-factor for the time-step condition \ No newline at end of file + cooling_tstep_mult: 1. # Dimensionless pre-factor for the time-step condition diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml index 5b336bc4bf86d2596ec52e84e24da2710d96cb8a..b90e020803be0e58c4e3539ee0bdfb8dbe540ace 100644 --- a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_low_z.yml @@ -44,7 +44,19 @@ InitialConditions: periodic: 1 # Parameters for the EAGLE chemistry -EAGLEChemistry: # Primordial +# EAGLEChemistry: # Solar abundances +# init_abundance_metal: 0.014 +# init_abundance_Hydrogen: 0.70649785 +# init_abundance_Helium: 0.28055534 +# init_abundance_Carbon: 2.0665436e-3 +# init_abundance_Nitrogen: 8.3562563e-4 +# init_abundance_Oxygen: 5.4926244e-3 +# init_abundance_Neon: 1.4144605e-3 +# init_abundance_Magnesium: 5.907064e-4 +# init_abundance_Silicon: 6.825874e-4 +# init_abundance_Iron: 1.1032152e-3 + +EAGLEChemistry: # Primordial abundances init_abundance_metal: 0. init_abundance_Hydrogen: 0.752 init_abundance_Helium: 0.248 diff --git a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml index 8bce6498432b3782225c7e53397868fc1356a26e..bb42bf5cb2cd0b8f7be9c87083262f0cfe5c9db6 100644 --- a/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml +++ b/examples/Cooling/CoolingRedshiftDependence/cooling_redshift_dependence_no_z.yml @@ -35,7 +35,20 @@ InitialConditions: periodic: 1 # Parameters for the EAGLE chemistry -EAGLEChemistry: # Solar abundances +# EAGLEChemistry: # Solar abundances +# init_abundance_metal: 0.014 +# init_abundance_Hydrogen: 0.70649785 +# init_abundance_Helium: 0.28055534 +# init_abundance_Carbon: 2.0665436e-3 +# init_abundance_Nitrogen: 8.3562563e-4 +# init_abundance_Oxygen: 5.4926244e-3 +# init_abundance_Neon: 1.4144605e-3 +# init_abundance_Magnesium: 5.907064e-4 +# init_abundance_Silicon: 6.825874e-4 +# init_abundance_Iron: 1.1032152e-3 + +# Parameters for the EAGLE chemistry +EAGLEChemistry: # Primoridal abundances init_abundance_metal: 0. init_abundance_Hydrogen: 0.752 init_abundance_Helium: 0.248 @@ -73,4 +86,4 @@ LambdaCooling: ConstCooling: cooling_rate: 1. # Cooling rate (du/dt) (internal units) min_energy: 1. # Minimal internal energy per unit mass (internal units) - cooling_tstep_mult: 1. # Dimensionless pre-factor for the time-step condition \ No newline at end of file + cooling_tstep_mult: 1. # Dimensionless pre-factor for the time-step condition diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml index b7e29d8df7b3510df97079d4a4ac85081498b3eb..e4e084f30020872befa7641f59c876dbcf2ef669 100644 --- a/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml +++ b/examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml @@ -10,8 +10,8 @@ InternalUnitSystem: Gravity: eta: 0.025 # Constant dimensionless multiplier for time integration. theta: 0.7 # Opening angle (Multipole acceptance criterion). - comoving_softening: 0.01 # Comoving softening length (in internal units). - max_physical_softening: 0.01 # Physical softening length (in internal units). + comoving_softening: 0.1 # Comoving softening length (in internal units). + max_physical_softening: 0.1 # Physical softening length (in internal units). # Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.) TimeIntegration: @@ -69,16 +69,16 @@ EAGLEChemistry: # Hernquist potential parameters HernquistPotential: - useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions + useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units) - idealizeddisk: 1 # Run with an idealized galaxy disk + idealizeddisk: 1 # Run with an idealized galaxy disk M200: 137.0 # M200 of the galaxy disk - h: 0.704 # reduced Hubble constant (value does not specify the used units!) - concentration: 9.0 # concentration of the Halo - diskfraction: 0.040 # Disk mass fraction - bulgefraction: 0.014 # Bulge mass fraction - timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration - epsilon: 0.01 # Softening size (internal units) + h: 0.704 # reduced Hubble constant (value does not specify the used units!) + concentration: 9.0 # concentration of the Halo + diskfraction: 0.040 # Disk mass fraction + bulgefraction: 0.014 # Bulge mass fraction + timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration + epsilon: 0.1 # Softening size (internal units) # EAGLE star formation parameters EAGLEStarFormation: diff --git a/src/cell.c b/src/cell.c index c1ee647717c4793508c726ca6a0ac5834388e8de..b52b51fea1749124f36ec4fc36f0d7d0351dbe2d 100644 --- a/src/cell.c +++ b/src/cell.c @@ -53,6 +53,7 @@ #include "chemistry.h" #include "drift.h" #include "engine.h" +#include "entropy_floor.h" #include "error.h" #include "feedback.h" #include "gravity.h" @@ -4182,7 +4183,8 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) { /* Drift... */ drift_part(p, xp, dt_drift, dt_kick_hydro, dt_kick_grav, dt_therm, - ti_old_part, ti_current); + ti_old_part, ti_current, e->cosmology, e->hydro_properties, + e->entropy_floor); /* Update the tracers properties */ tracers_after_drift(p, xp, e->internal_units, e->physical_constants, diff --git a/src/chemistry/EAGLE/chemistry.h b/src/chemistry/EAGLE/chemistry.h index 8c0061eef618fa7c01c48ae0ace763a20d85e762..a04d7f94e7340a68d25914828de40bb2e1e8020c 100644 --- a/src/chemistry/EAGLE/chemistry.h +++ b/src/chemistry/EAGLE/chemistry.h @@ -90,8 +90,9 @@ __attribute__((always_inline)) INLINE static void chemistry_end_density( /* Some smoothing length multiples. */ const float h = p->h; - const float h_inv = 1.0f / h; /* 1/h */ - const float factor = pow_dimension(h_inv) / p->rho; /* 1 / h^d * rho */ + const float h_inv = 1.0f / h; /* 1/h */ + const float rho = hydro_get_comoving_density(p); + const float factor = pow_dimension(h_inv) / rho; /* 1 / (h^d * rho) */ const float m = hydro_get_mass(p); struct chemistry_part_data* cpd = &p->chemistry_data; diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c index af7a4c90eee1e581f24e6b303ee3c99c721b530b..0b5fcde02e4005c1d8c6a3bc4ff830e435b35ad6 100644 --- a/src/cooling/EAGLE/cooling.c +++ b/src/cooling/EAGLE/cooling.c @@ -48,18 +48,13 @@ #include "space.h" #include "units.h" -/* Maximum number of iterations for newton - * and bisection integration schemes */ -static const int newton_max_iterations = 15; +/* Maximum number of iterations for bisection scheme */ static const int bisection_max_iterations = 150; /* Tolerances for termination criteria. */ static const float explicit_tolerance = 0.05; -static const float newton_tolerance = 1.0e-4; static const float bisection_tolerance = 1.0e-6; -static const float rounding_tolerance = 1.0e-4; static const double bracket_factor = 1.5; -static const double newton_log_u_guess_cgs = 12; /** * @brief Find the index of the current redshift along the redshift dimension @@ -189,102 +184,6 @@ void cooling_update(const struct cosmology *cosmo, cooling->z_index = z_index; } -/** - * @brief Newton Raphson integration scheme to calculate particle cooling over - * timestep. This replaces bisection scheme used in EAGLE to minimize the - * number of array accesses. Integration defaults to bisection scheme (see - * function bisection_iter) if this function does not converge within a - * specified number of steps - * - * @param logu_init Initial guess for log(internal energy) - * @param u_ini Internal energy at beginning of hydro step - * @param n_H_index Particle hydrogen number density index - * @param d_n_H Particle hydrogen number density offset - * @param He_index Particle helium fraction index - * @param d_He Particle helium fraction offset - * @param He_reion_heat Heating due to helium reionization - * (only depends on redshift, so passed as parameter) - * @param p #part structure - * @param cosmo #cosmology structure - * @param cooling #cooling_function_data structure - * @param phys_const #phys_const data structure - * @param abundance_ratio Array of ratios of metal abundance to solar - * @param dt timestep - * @param bisection_flag Flag to identify if scheme failed to converge - */ -INLINE static float newton_iter( - float logu_init, double u_ini, int n_H_index, float d_n_H, int He_index, - float d_He, float He_reion_heat, struct part *restrict p, - const struct cosmology *restrict cosmo, - const struct cooling_function_data *restrict cooling, - const struct phys_const *restrict phys_const, - const float abundance_ratio[chemistry_element_count + 2], float dt, - int *bisection_flag) { - - double logu, logu_old; - double dLambdaNet_du = 0.0, LambdaNet; - - /* table bounds */ - const float log_table_bound_high = - (cooling->Therm[eagle_cooling_N_temperature - 1] - 0.05) / M_LOG10E; - const float log_table_bound_low = (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 double n_H = - hydro_get_physical_density(p, cosmo) * XH / phys_const->const_proton_mass; - const double n_H_cgs = n_H * cooling->number_density_to_cgs; - - /* compute ratefact = n_H * n_H / rho; Might lead to round-off error: - * replaced by equivalent expression below */ - const double ratefact_cgs = n_H_cgs * XH * cooling->inv_proton_mass_cgs; - - logu_old = logu_init; - logu = logu_old; - int i = 0; - - float LambdaNet_old = 0; - LambdaNet = 0; - do /* iterate to convergence */ - { - logu_old = logu; - LambdaNet_old = LambdaNet; - LambdaNet = (He_reion_heat / (dt * ratefact_cgs)) + - eagle_cooling_rate(logu_old, cosmo->z, n_H_cgs, abundance_ratio, - n_H_index, d_n_H, He_index, d_He, cooling, - &dLambdaNet_du); - - /* Newton iteration. For details on how the cooling equation is integrated - * see documentation in theory/Cooling/ */ - logu = logu_old - (1.0 - u_ini * exp(-logu_old) - - LambdaNet * ratefact_cgs * dt * exp(-logu_old)) / - (1.0 - dLambdaNet_du * ratefact_cgs * dt); - /* Check if first step passes over equilibrium solution, if it does adjust - * next guess */ - if (i == 1 && LambdaNet_old * LambdaNet < 0) logu = newton_log_u_guess_cgs; - - /* check whether iterations go within about 10% of the table bounds, - * if they do default to bisection method */ - if (logu > log_table_bound_high) { - i = newton_max_iterations; - break; - } else if (logu < log_table_bound_low) { - i = newton_max_iterations; - break; - } - - i++; - } while (fabs(logu - logu_old) > newton_tolerance && - i < newton_max_iterations); - if (i >= newton_max_iterations) { - /* flag to trigger bisection scheme */ - *bisection_flag = 1; - } - - return logu; -} - /** * @brief Bisection integration scheme * @@ -304,11 +203,12 @@ INLINE static float newton_iter( */ INLINE static double bisection_iter( const double u_ini_cgs, const double n_H_cgs, const double redshift, - int n_H_index, float d_n_H, int He_index, float d_He, - double Lambda_He_reion_cgs, double ratefact_cgs, + const int n_H_index, const float d_n_H, const int He_index, + const float d_He, const double Lambda_He_reion_cgs, + const double ratefact_cgs, const struct cooling_function_data *restrict cooling, - const float abundance_ratio[chemistry_element_count + 2], double dt_cgs, - long long ID) { + const float abundance_ratio[chemistry_element_count + 2], + const double dt_cgs, const long long ID) { /* Bracketing */ double u_lower_cgs = u_ini_cgs; @@ -320,9 +220,8 @@ INLINE static double bisection_iter( double LambdaNet_cgs = Lambda_He_reion_cgs + - eagle_cooling_rate(log(u_ini_cgs), redshift, n_H_cgs, abundance_ratio, - n_H_index, d_n_H, He_index, d_He, cooling, - /*dLambdaNet_du=*/NULL); + eagle_cooling_rate(log10(u_ini_cgs), redshift, n_H_cgs, abundance_ratio, + n_H_index, d_n_H, He_index, d_He, cooling); /*************************************/ /* Let's try to bracket the solution */ @@ -335,11 +234,10 @@ INLINE static double bisection_iter( u_upper_cgs *= bracket_factor; /* Compute a new rate */ - LambdaNet_cgs = - Lambda_He_reion_cgs + - eagle_cooling_rate(log(u_lower_cgs), redshift, n_H_cgs, abundance_ratio, - n_H_index, d_n_H, He_index, d_He, cooling, - /*dLambdaNet_du=*/NULL); + LambdaNet_cgs = Lambda_He_reion_cgs + + eagle_cooling_rate(log10(u_lower_cgs), redshift, n_H_cgs, + abundance_ratio, n_H_index, d_n_H, + He_index, d_He, cooling); int i = 0; while (u_lower_cgs - u_ini_cgs - LambdaNet_cgs * ratefact_cgs * dt_cgs > @@ -351,10 +249,9 @@ INLINE static double bisection_iter( /* Compute a new rate */ LambdaNet_cgs = Lambda_He_reion_cgs + - eagle_cooling_rate(log(u_lower_cgs), redshift, n_H_cgs, + eagle_cooling_rate(log10(u_lower_cgs), redshift, n_H_cgs, abundance_ratio, n_H_index, d_n_H, - He_index, d_He, cooling, - /*dLambdaNet_du=*/NULL); + He_index, d_He, cooling); i++; } @@ -371,11 +268,10 @@ INLINE static double bisection_iter( u_upper_cgs *= bracket_factor; /* Compute a new rate */ - LambdaNet_cgs = - Lambda_He_reion_cgs + - eagle_cooling_rate(log(u_upper_cgs), redshift, n_H_cgs, abundance_ratio, - n_H_index, d_n_H, He_index, d_He, cooling, - /*dLambdaNet_du=*/NULL); + LambdaNet_cgs = Lambda_He_reion_cgs + + eagle_cooling_rate(log10(u_upper_cgs), redshift, n_H_cgs, + abundance_ratio, n_H_index, d_n_H, + He_index, d_He, cooling); int i = 0; while (u_upper_cgs - u_ini_cgs - LambdaNet_cgs * ratefact_cgs * dt_cgs < @@ -387,10 +283,9 @@ INLINE static double bisection_iter( /* Compute a new rate */ LambdaNet_cgs = Lambda_He_reion_cgs + - eagle_cooling_rate(log(u_upper_cgs), redshift, n_H_cgs, + eagle_cooling_rate(log10(u_upper_cgs), redshift, n_H_cgs, abundance_ratio, n_H_index, d_n_H, - He_index, d_He, cooling, - /*dLambdaNet_du=*/NULL); + He_index, d_He, cooling); i++; } @@ -417,11 +312,10 @@ INLINE static double bisection_iter( u_next_cgs = 0.5 * (u_lower_cgs + u_upper_cgs); /* New rate */ - LambdaNet_cgs = - Lambda_He_reion_cgs + - eagle_cooling_rate(log(u_next_cgs), redshift, n_H_cgs, abundance_ratio, - n_H_index, d_n_H, He_index, d_He, cooling, - /*dLambdaNet_du=*/NULL); + LambdaNet_cgs = Lambda_He_reion_cgs + + eagle_cooling_rate(log10(u_next_cgs), redshift, n_H_cgs, + abundance_ratio, n_H_index, d_n_H, + He_index, d_He, cooling); // Debugging if (u_next_cgs <= 0) @@ -503,7 +397,7 @@ void cooling_cool_part(const struct phys_const *phys_const, /* Get the change in internal energy due to hydro forces */ const float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo); - /* Get internal energy at the end of the next kick step (assuming dt does not + /* Get internal energy at the end of the step (assuming dt does not * increase) */ double u_0 = (u_start + hydro_du_dt * dt_therm); @@ -511,7 +405,6 @@ void cooling_cool_part(const struct phys_const *phys_const, u_0 = max(u_0, hydro_properties->minimal_internal_energy); /* Convert to CGS units */ - const double u_start_cgs = u_start * cooling->internal_energy_to_cgs; const double u_0_cgs = u_0 * cooling->internal_energy_to_cgs; const double dt_cgs = dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME); @@ -536,7 +429,7 @@ void cooling_cool_part(const struct phys_const *phys_const, * metal-free Helium mass fraction as per the Wiersma+08 definition */ const float HeFrac = XHe / (XH + XHe); - /* convert Hydrogen mass fraction into Hydrogen number density */ + /* convert Hydrogen mass fraction into physical Hydrogen number density */ const double n_H = hydro_get_physical_density(p, cosmo) * XH / phys_const->const_proton_mass; const double n_H_cgs = n_H * cooling->number_density_to_cgs; @@ -573,10 +466,9 @@ void cooling_cool_part(const struct phys_const *phys_const, /* First try an explicit integration (note we ignore the derivative) */ const double LambdaNet_cgs = - Lambda_He_reion_cgs + eagle_cooling_rate(log(u_0_cgs), cosmo->z, n_H_cgs, - abundance_ratio, n_H_index, - d_n_H, He_index, d_He, cooling, - /*dLambdaNet_du=*/NULL); + Lambda_He_reion_cgs + + eagle_cooling_rate(log10(u_0_cgs), cosmo->z, n_H_cgs, abundance_ratio, + n_H_index, d_n_H, He_index, d_He, cooling); /* if cooling rate is small, take the explicit solution */ if (fabs(ratefact_cgs * LambdaNet_cgs * dt_cgs) < @@ -586,76 +478,32 @@ void cooling_cool_part(const struct phys_const *phys_const, } else { - int bisection_flag = 1; - - // MATTHIEU: TO DO restore the Newton-Raphson scheme - if (0 && cooling->newton_flag) { - - /* Ok, try a Newton-Raphson scheme instead */ - double log_u_final_cgs = - newton_iter(log(u_0_cgs), u_0_cgs, n_H_index, d_n_H, He_index, d_He, - Lambda_He_reion_cgs, p, cosmo, cooling, phys_const, - abundance_ratio, dt_cgs, &bisection_flag); - - /* Check if newton scheme sent us to a higher energy despite being in - a cooling regime If it did try newton scheme with a better guess. - (Guess internal energy near equilibrium solution). */ - if (LambdaNet_cgs < 0 && log_u_final_cgs > log(u_0_cgs)) { - bisection_flag = 0; - log_u_final_cgs = - newton_iter(newton_log_u_guess_cgs, u_0_cgs, n_H_index, d_n_H, - He_index, d_He, Lambda_He_reion_cgs, p, cosmo, cooling, - phys_const, abundance_ratio, dt_cgs, &bisection_flag); - } - - u_final_cgs = exp(log_u_final_cgs); - } - - /* Alright, all else failed, let's bisect */ - if (bisection_flag || !(cooling->newton_flag)) { - u_final_cgs = - bisection_iter(u_0_cgs, n_H_cgs, cosmo->z, n_H_index, d_n_H, He_index, - d_He, Lambda_He_reion_cgs, ratefact_cgs, cooling, - abundance_ratio, dt_cgs, p->id); - } + /* Otherwise, go the bisection route. */ + u_final_cgs = + bisection_iter(u_0_cgs, n_H_cgs, cosmo->z, n_H_index, d_n_H, He_index, + d_He, Lambda_He_reion_cgs, ratefact_cgs, cooling, + abundance_ratio, dt_cgs, p->id); } - /* 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; - /* Convert back to internal units */ - double delta_u = delta_u_cgs * cooling->internal_energy_from_cgs; + double u_final = u_final_cgs * cooling->internal_energy_from_cgs; /* We now need to check that we are not going to go below any of the limits */ - /* Limit imposed by the entropy floor */ - const double A_floor = entropy_floor(p, cosmo, floor_props); - const double rho = hydro_get_physical_density(p, cosmo); - const double u_floor = gas_internal_energy_from_entropy(rho, A_floor); - /* Absolute minimum */ const double u_minimal = hydro_properties->minimal_internal_energy; + u_final = max(u_final, u_minimal); - /* Largest of both limits */ - const double u_limit = max(u_minimal, u_floor); - - /* First, check whether we may end up below the minimal energy after - * this step 1/2 kick + another 1/2 kick that could potentially be for - * a time-step twice as big. We hence check for 1.5 delta_u. */ - if (u_start + 1.5 * delta_u < u_limit) { - delta_u = (u_limit - u_start) / 1.5; - } + /* Limit imposed by the entropy floor */ + const double A_floor = entropy_floor(p, cosmo, floor_props); + const double rho_physical = hydro_get_physical_density(p, cosmo); + const double u_floor = + gas_internal_energy_from_entropy(rho_physical, A_floor); + u_final = max(u_final, u_floor); - /* Second, check whether the energy used in the prediction could get negative. - * We need to check for the 1/2 dt kick followed by a full time-step drift - * 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.) { - delta_u = -u_start / (2.5 + rounding_tolerance); - } + /* Expected change in energy over the next kick step + (assuming no change in dt) */ + const double delta_u = u_final - max(u_start, u_floor); /* Turn this into a rate of change (including cosmology term) */ const float cooling_du_dt = delta_u / dt_therm; @@ -774,8 +622,7 @@ float cooling_get_temperature( /* Compute the log10 of the temperature by interpolating the table */ const double log_10_T = eagle_convert_u_to_temp( - log10(u_cgs), cosmo->z, /*compute_dT_du=*/0, /*dT_du=*/NULL, n_H_index, - He_index, d_n_H, d_He, cooling); + log10(u_cgs), cosmo->z, n_H_index, He_index, d_n_H, d_He, cooling); /* Undo the log! */ return exp10(log_10_T); @@ -940,10 +787,6 @@ void cooling_init_backend(struct swift_params *parameter_file, /* set previous_z_index and to last value of redshift table*/ cooling->previous_z_index = eagle_cooling_N_redshifts - 2; - - /* Check if we are running with the newton scheme */ - cooling->newton_flag = parser_get_opt_param_int( - parameter_file, "EAGLECooling:newton_integration", 0); } /** diff --git a/src/cooling/EAGLE/cooling_rates.h b/src/cooling/EAGLE/cooling_rates.h index 35a53b43fb26466686fa87704ec174a849667a78..6707a9a9d70b58eb218cfc8ebeca9b1f08fdcb12 100644 --- a/src/cooling/EAGLE/cooling_rates.h +++ b/src/cooling/EAGLE/cooling_rates.h @@ -99,7 +99,7 @@ __attribute__((always_inline)) INLINE void abundance_ratio_to_solar( cooling->Ca_over_Si_ratio_in_solar; ratio_solar[10] = - p->chemistry_data.metal_mass_fraction[chemistry_element_Fe] * + p->chemistry_data.smoothed_metal_mass_fraction[chemistry_element_Fe] * cooling->SolarAbundances_inv[10 /* Fe */]; } @@ -118,7 +118,8 @@ __attribute__((always_inline)) INLINE void abundance_ratio_to_solar( */ __attribute__((always_inline)) INLINE double eagle_helium_reionization_extraheat( - double z, double delta_z, const struct cooling_function_data *cooling) { + const double z, const double delta_z, + const struct cooling_function_data *cooling) { #ifdef SWIFT_DEBUG_CHECKS if (delta_z > 0.f) error("Invalid value for delta_z. Should be negative."); @@ -169,15 +170,12 @@ eagle_helium_reionization_extraheat( * @param d_He Offset between helium fraction and table[He_index]. * @param cooling #cooling_function_data structure. * - * @param compute_dT_du Do we want to compute dT/du ? - * @param dT_du (return) The value of dT/du - * * @return log_10 of the temperature. */ __attribute__((always_inline)) INLINE double eagle_convert_u_to_temp( - const double log_10_u_cgs, const float redshift, const int compute_dT_du, - float *dT_du, int n_H_index, int He_index, float d_n_H, float d_He, - const struct cooling_function_data *restrict cooling) { + const double log_10_u_cgs, const float redshift, const int n_H_index, + const int He_index, const float d_n_H, const float d_He, + const struct cooling_function_data *cooling) { /* Get index of u along the internal energy axis */ int u_index; @@ -209,67 +207,6 @@ __attribute__((always_inline)) INLINE double eagle_convert_u_to_temp( eagle_cooling_N_temperature); /* */ } - if (compute_dT_du) { - - float log_10_T_high, log_10_T_low; - - /* Interpolate temperature table to return temperature for internal energy - * at grid point above current internal energy for computing dT_du used for - * calculation of dlambda_du in cooling.c (use 3D interpolation for high - * redshift table, otherwise 4D) */ - if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) { - - log_10_T_high = interpolation_3d(cooling->table.temperature, /* */ - n_H_index, He_index, u_index, /* */ - d_n_H, d_He, /*delta_u=*/1.f, /* */ - eagle_cooling_N_density, /* */ - eagle_cooling_N_He_frac, /* */ - eagle_cooling_N_temperature); /* */ - - } else { - - log_10_T_high = - interpolation_4d(cooling->table.temperature, /* */ - /*z_index=*/0, n_H_index, He_index, u_index, /* */ - cooling->dz, d_n_H, d_He, /*delta_u=*/1.f, /* */ - eagle_cooling_N_loaded_redshifts, /* */ - eagle_cooling_N_density, /* */ - eagle_cooling_N_He_frac, /* */ - eagle_cooling_N_temperature); /* */ - } - - /* Interpolate temperature table to return temperature for internal energy - * at grid point below current internal energy for computing dT_du used for - * calculation of dlambda_du in cooling.c (use 3D interpolation for high - * redshift table, otherwise 4D) */ - if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) { - - log_10_T_low = interpolation_3d(cooling->table.temperature, /* */ - n_H_index, He_index, u_index, /* */ - d_n_H, d_He, /*delta_u=*/0.f, /* */ - eagle_cooling_N_density, /* */ - eagle_cooling_N_He_frac, /* */ - eagle_cooling_N_temperature); /* */ - - } else { - - log_10_T_low = - interpolation_4d(cooling->table.temperature, /* */ - /*z_index=*/0, n_H_index, He_index, u_index, /* */ - cooling->dz, d_n_H, d_He, /*delta_u=*/0.f, /* */ - eagle_cooling_N_loaded_redshifts, /* */ - eagle_cooling_N_density, /* */ - eagle_cooling_N_He_frac, /* */ - eagle_cooling_N_temperature); /* */ - } - - /* Calculate dT/du */ - const float delta_u = exp(cooling->Therm[u_index + 1] * M_LN10) - - exp(cooling->Therm[u_index] * M_LN10); - *dT_du = - (exp(M_LN10 * log_10_T_high) - exp(M_LN10 * log_10_T_low)) / delta_u; - } - /* Special case for temperatures below the start of the table */ if (u_index == 0 && d_u == 0.f) { @@ -355,10 +292,6 @@ __attribute__((always_inline)) INLINE double eagle_Compton_cooling_rate( * the redshift, hydrogen number density and helium fraction indices * and offsets passed in. * - * If the arguement dlambda_du is non-NULL, the routine also - * calculates derivative of cooling rate with respect to internal - * energy. - * * If the argument element_lambda is non-NULL, the routine also * returns the cooling rate per element in the array. * @@ -374,36 +307,19 @@ __attribute__((always_inline)) INLINE double eagle_Compton_cooling_rate( * @param d_He Particle helium fraction offset * @param cooling Cooling data structure * - * @param dlambda_du (return) Derivative of the cooling rate with respect to u. * @param element_lambda (return) Cooling rate from each element * * @return The cooling rate */ INLINE static double eagle_metal_cooling_rate( - double log10_u_cgs, double redshift, double n_H_cgs, - const float solar_ratio[chemistry_element_count + 2], int n_H_index, - float d_n_H, int He_index, float d_He, - const struct cooling_function_data *restrict cooling, double *dlambda_du, - double *element_lambda) { - -#ifdef TO_BE_DONE - /* used for calculating dlambda_du */ - double temp_lambda_high = 0, temp_lambda_low = 0; - double h_plus_he_electron_abundance_high = 0; - double h_plus_he_electron_abundance_low = 0; - double solar_electron_abundance_high = 0; - double solar_electron_abundance_low = 0; - double elem_cool_low = 0, elem_cool_high = 0; -#endif - - /* We only need dT_du if dLambda_du is non-NULL */ - const int compute_dT_du = (dlambda_du != NULL) ? 1 : 0; + const double log10_u_cgs, const double redshift, const double n_H_cgs, + const float solar_ratio[chemistry_element_count + 2], const int n_H_index, + const float d_n_H, const int He_index, const float d_He, + const struct cooling_function_data *cooling, double *element_lambda) { /* Temperature */ - float dT_du = -1.f; - const double log_10_T = - eagle_convert_u_to_temp(log10_u_cgs, redshift, compute_dT_du, &dT_du, - n_H_index, He_index, d_n_H, d_He, cooling); + const double log_10_T = eagle_convert_u_to_temp( + log10_u_cgs, redshift, n_H_index, He_index, d_n_H, d_He, cooling); /* Get index along temperature dimension of the tables */ int T_index; @@ -411,12 +327,6 @@ INLINE static double eagle_metal_cooling_rate( get_index_1d(cooling->Temp, eagle_cooling_N_temperature, log_10_T, &T_index, &d_T); -#ifdef TO_BE_DONE - /* Difference between entries on the temperature table around u */ - const float delta_T = exp(M_LN10 * cooling->Temp[T_index + 1]) - - exp(M_LN10 * cooling->Temp[T_index]); -#endif - /**********************/ /* Metal-free cooling */ /**********************/ @@ -433,21 +343,6 @@ INLINE static double eagle_metal_cooling_rate( eagle_cooling_N_density, /* */ eagle_cooling_N_He_frac, /* */ eagle_cooling_N_temperature); /* */ - -#ifdef TO_BE_DONE - /* compute values at temperature gridpoints above and below input - * temperature for calculation of dlambda_du. Pass in NULL pointer for - * dlambda_du in order to skip */ - if (dlambda_du != NULL) { - temp_lambda_high = interpolation_3d( - cooling->table.H_plus_He_heating, n_H_index, He_index, T_index, d_n_h, - d_He, 1.f, cooling->N_nH, cooling->N_He, cooling->N_Temp); - temp_lambda_low = interpolation_3d( - cooling->table.H_plus_He_heating, n_H_index, He_index, T_index, d_n_h, - d_He, 0.f, cooling->N_nH, cooling->N_He, cooling->N_Temp); - } -#endif - } else { /* Using normal tables, have to interpolate in redshift */ @@ -459,29 +354,8 @@ INLINE static double eagle_metal_cooling_rate( eagle_cooling_N_density, /* */ eagle_cooling_N_He_frac, /* */ eagle_cooling_N_temperature); /* */ - -#ifdef TO_BE_DONE - /* compute values at temperature gridpoints above and below input - * temperature for calculation of dlambda_du */ - if (dlambda_du != NULL) { - temp_lambda_high = - interpolation_4d(cooling->table.H_plus_He_heating, 0, n_H_index, - He_index, T_index, cooling->dz, d_n_h, d_He, 1.f, 2, - cooling->N_nH, cooling->N_He, cooling->N_Temp); - temp_lambda_low = - interpolation_4d(cooling->table.H_plus_He_heating, 0, n_H_index, - He_index, T_index, cooling->dz, d_n_h, d_He, 0.f, 2, - cooling->N_nH, cooling->N_He, cooling->N_Temp); - } -#endif } -#ifdef TO_BE_DONE - if (dlambda_du != NULL) { - *dlambda_du += (temp_lambda_high - temp_lambda_low) / delta_T * dT_du; - } -#endif - /* If we're testing cooling rate contributions write to array */ if (element_lambda != NULL) { element_lambda[0] = Lambda_free; @@ -502,21 +376,6 @@ INLINE static double eagle_metal_cooling_rate( eagle_cooling_N_density, /* */ eagle_cooling_N_He_frac, /* */ eagle_cooling_N_temperature); /* */ -#ifdef TO_BE_DONE - /* compute values at temperature gridpoints above and below input - * temperature for calculation of dlambda_du. Pass in NULL pointer for - * dlambda_du in order to skip */ - - h_plus_he_electron_abundance_high = - interpolation_3d(cooling->table.H_plus_He_electron_abundance, n_H_index, - He_index, T_index, d_n_h, d_He, 1.f, cooling->N_nH, - cooling->N_He, cooling->N_Temp); - h_plus_he_electron_abundance_low = - interpolation_3d(cooling->table.H_plus_He_electron_abundance, n_H_index, - He_index, T_index, d_n_h, d_He, 0.f, cooling->N_nH, - cooling->N_He, cooling->N_Temp); - -#endif } else { @@ -528,19 +387,6 @@ INLINE static double eagle_metal_cooling_rate( eagle_cooling_N_density, /* */ eagle_cooling_N_He_frac, /* */ eagle_cooling_N_temperature); /* */ - -#ifdef TO_BE_DONE - /* compute values at temperature gridpoints above and below input - * temperature for calculation of dlambda_du */ - h_plus_he_electron_abundance_high = - interpolation_4d(cooling->table.H_plus_He_electron_abundance, 0, - n_H_index, He_index, T_index, cooling->dz, d_n_h, d_He, - 1.f, 2, cooling->N_nH, cooling->N_He, cooling->N_Temp); - h_plus_he_electron_abundance_low = - interpolation_4d(cooling->table.H_plus_He_electron_abundance, 0, - n_H_index, He_index, T_index, cooling->dz, d_n_h, d_He, - 0.f, 2, cooling->N_nH, cooling->N_He, cooling->N_Temp); -#endif } /**********************/ @@ -583,19 +429,6 @@ INLINE static double eagle_metal_cooling_rate( eagle_cooling_N_density, /* */ eagle_cooling_N_temperature); /* */ -#ifdef TO_BE_DONE - /* compute values at temperature gridpoints above and below input - * temperature for calculation of dlambda_du */ - if (dlambda_du != NULL) { - solar_electron_abundance_high = - interpolation_2d(cooling->table.electron_abundance, n_H_index, - T_index, d_n_h, 1.f, cooling->N_nH, cooling->N_Temp); - solar_electron_abundance_low = - interpolation_2d(cooling->table.electron_abundance, n_H_index, - T_index, d_n_h, 0.f, cooling->N_nH, cooling->N_Temp); - } -#endif - } else { /* Using normal tables, have to interpolate in redshift */ @@ -606,19 +439,6 @@ INLINE static double eagle_metal_cooling_rate( eagle_cooling_N_loaded_redshifts, /* */ eagle_cooling_N_density, /* */ eagle_cooling_N_temperature); /* */ - -#ifdef TO_BE_DONE - /* compute values at temperature gridpoints above and below input - * temperature for calculation of dlambda_du */ - if (dlambda_du != NULL) { - solar_electron_abundance_high = interpolation_3d( - cooling->table.electron_abundance, 0, n_H_index, T_index, cooling->dz, - d_n_h, 1.f, 2, cooling->N_nH, cooling->N_Temp); - solar_electron_abundance_low = interpolation_3d( - cooling->table.electron_abundance, 0, n_H_index, T_index, cooling->dz, - d_n_h, 0.f, 2, cooling->N_nH, cooling->N_Temp); - } -#endif } const double electron_abundance_ratio = @@ -654,26 +474,6 @@ INLINE static double eagle_metal_cooling_rate( lambda_metal[elem] *= electron_abundance_ratio; lambda_metal[elem] *= solar_ratio[elem]; } - -#ifdef TO_BE_DONE - /* compute values at temperature gridpoints above and below input - * temperature for calculation of dlambda_du */ - if (dlambda_du != NULL) { - elem_cool_high = interpolation_3d_no_x( - cooling->table.metal_heating, elem, n_H_index, T_index, 0.f, d_n_h, - 1.f, cooling->N_Elements, cooling->N_nH, cooling->N_Temp); - - elem_cool_low = interpolation_3d_no_x( - cooling->table.metal_heating, elem, n_H_index, T_index, 0.f, d_n_h, - 0.f, cooling->N_nH, cooling->N_Temp, cooling->N_Elements); - - *dlambda_du += (elem_cool_high * h_plus_he_electron_abundance_high / - solar_electron_abundance_high - - elem_cool_low * h_plus_he_electron_abundance_low / - solar_electron_abundance_low) / - delta_T * dT_du * solar_ratio[elem + 2]; - } -#endif } } else { @@ -696,29 +496,9 @@ INLINE static double eagle_metal_cooling_rate( lambda_metal[elem] *= electron_abundance_ratio; lambda_metal[elem] *= solar_ratio[elem]; - } -#ifdef TO_BE_DONE - /* compute values at temperature gridpoints above and below input - * temperature for calculation of dlambda_du */ - if (dlambda_du != NULL) { - elem_cool_high = interpolation_4d_no_x( - cooling->table.metal_heating, elem, 0, n_H_index, T_index, 0., - cooling->dz, d_n_h, 1.f, cooling->N_Elements, 2, cooling->N_nH, - cooling->N_Temp); - - elem_cool_low = interpolation_4d_no_x( - cooling->table.metal_heating, elem, 0, n_H_index, T_index, 0., - cooling->dz, d_n_h, 0.f, cooling->N_Elements, 2, cooling->N_nH, - cooling->N_Temp); - - *dlambda_du += (elem_cool_high * h_plus_he_electron_abundance_high / - solar_electron_abundance_high - - elem_cool_low * h_plus_he_electron_abundance_low / - solar_electron_abundance_low) / - delta_T * dT_du * solar_ratio[elem + 2]; + // message("lambda[%d]=%e", elem, lambda_metal[elem]); } -#endif } } @@ -738,11 +518,11 @@ INLINE static double eagle_metal_cooling_rate( } /** - * @brief Wrapper function used to calculate cooling rate and dLambda_du. + * @brief Wrapper function used to calculate cooling rate. * Table indices and offsets for redshift, hydrogen number density and * helium fraction are passed it so as to compute them only once per particle. * - * @param log_u_cgs Natural log of internal energy per unit mass in CGS units. + * @param log10_u_cgs Log base 10 of internal energy per unit mass in CGS units. * @param redshift The current redshift. * @param n_H_cgs Hydrogen number density in CGS units. * @param abundance_ratio Ratio of element abundance to solar. @@ -753,22 +533,17 @@ INLINE static double eagle_metal_cooling_rate( * @param d_He Particle helium fraction offset * @param cooling #cooling_function_data structure * - * @param dLambdaNet_du (return) Derivative of the cooling rate with respect to - * u. - * * @return The cooling rate */ INLINE static double eagle_cooling_rate( - double log_u_cgs, double redshift, double n_H_cgs, - const float abundance_ratio[chemistry_element_count + 2], int n_H_index, - float d_n_H, int He_index, float d_He, - const struct cooling_function_data *restrict cooling, - double *dLambdaNet_du) { + const double log10_u_cgs, const double redshift, const double n_H_cgs, + const float abundance_ratio[chemistry_element_count + 2], + const int n_H_index, const float d_n_H, const int He_index, + const float d_He, const struct cooling_function_data *cooling) { - return eagle_metal_cooling_rate(log_u_cgs / M_LN10, redshift, n_H_cgs, + return eagle_metal_cooling_rate(log10_u_cgs, redshift, n_H_cgs, abundance_ratio, n_H_index, d_n_H, He_index, - d_He, cooling, dLambdaNet_du, - /*element_lambda=*/NULL); + d_He, cooling, /* element_lambda=*/NULL); } #endif /* SWIFT_EAGLE_COOLING_RATES_H */ diff --git a/src/cooling/EAGLE/cooling_struct.h b/src/cooling/EAGLE/cooling_struct.h index 88e49036f62962d35025f7fa81d5703a1245e3ce..8c7c44e6bb79a1f9cce32967cb808d149dcb3740 100644 --- a/src/cooling/EAGLE/cooling_struct.h +++ b/src/cooling/EAGLE/cooling_struct.h @@ -112,7 +112,7 @@ struct cooling_function_data { /*! Inverse of proton mass in cgs (for quick access) */ double inv_proton_mass_cgs; - /*! Temperatur of the CMB at present day (for quick access) */ + /*! Temperature of the CMB at present day (for quick access) */ double T_CMB_0; /*! Compton rate in cgs units */ @@ -126,9 +126,6 @@ struct cooling_function_data { /*! Index of the previous tables along the redshift index of the tables */ int previous_z_index; - - /*! Are we doing Newton-Raphson iterations? */ - int newton_flag; }; /** diff --git a/src/cooling/EAGLE/cooling_tables.c b/src/cooling/EAGLE/cooling_tables.c index cd6678c70614aa62c8145afbecc2c8f673dbf0de..cddc8d50e0c2f00c30846ebaf65b2bf250e9cc8a 100644 --- a/src/cooling/EAGLE/cooling_tables.c +++ b/src/cooling/EAGLE/cooling_tables.c @@ -296,6 +296,7 @@ void read_cooling_header(const char *fname, for (int i = 0; i < N_nH; i++) { cooling->nH[i] = log10(cooling->nH[i]); } + /* Compute inverse of solar mass fractions */ for (int i = 0; i < N_SolarAbundances; ++i) { cooling->SolarAbundances_inv[i] = 1.f / cooling->SolarAbundances[i]; diff --git a/src/cooling/EAGLE/newton_cooling.c b/src/cooling/EAGLE/newton_cooling.c new file mode 100644 index 0000000000000000000000000000000000000000..c68a77614425f85da011ab45d2ec488710e068dc --- /dev/null +++ b/src/cooling/EAGLE/newton_cooling.c @@ -0,0 +1,678 @@ +/** + * @file Backup file containing the now-defunct Newton integration scheme + * for the cooling. + */ + +/** + * @brief Computes the log_10 of the temperature corresponding to a given + * internal energy, hydrogen number density, Helium fraction and redshift. + * + * Note that the redshift is implicitly passed in via the currently loaded + * tables in the #cooling_function_data. + * + * For the low-z case, we interpolate the flattened 4D table 'u_to_temp' that + * is arranged in the following way: + * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts + * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density + * - 3rd dim: Helium fraction, length = eagle_cooling_N_He_frac + * - 4th dim: Internal energy, length = eagle_cooling_N_temperature + * + * For the high-z case, we interpolate the flattened 3D table 'u_to_temp' that + * is arranged in the following way: + * - 1st dim: Hydrogen density, length = eagle_cooling_N_density + * - 2nd dim: Helium fraction, length = eagle_cooling_N_He_frac + * - 3rd dim: Internal energy, length = eagle_cooling_N_temperature + * + * @param log_10_u_cgs Log base 10 of internal energy in cgs. + * @param redshift Current redshift. + * @param n_H_index Index along the Hydrogen density dimension. + * @param He_index Index along the Helium fraction dimension. + * @param d_n_H Offset between Hydrogen density and table[n_H_index]. + * @param d_He Offset between helium fraction and table[He_index]. + * @param cooling #cooling_function_data structure. + * + * @param compute_dT_du Do we want to compute dT/du ? + * @param dT_du (return) The value of dT/du + * + * @return log_10 of the temperature. + */ +__attribute__((always_inline)) INLINE double eagle_convert_u_to_temp( + const double log_10_u_cgs, const float redshift, const int compute_dT_du, + float *dT_du, int n_H_index, int He_index, float d_n_H, float d_He, + const struct cooling_function_data *restrict cooling) { + + /* Get index of u along the internal energy axis */ + int u_index; + float d_u; + get_index_1d(cooling->Therm, eagle_cooling_N_temperature, log_10_u_cgs, + &u_index, &d_u); + + /* Interpolate temperature table to return temperature for current + * internal energy (use 3D interpolation for high redshift table, + * otherwise 4D) */ + float log_10_T; + if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) { + + log_10_T = interpolation_3d(cooling->table.temperature, /* */ + n_H_index, He_index, u_index, /* */ + d_n_H, d_He, d_u, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_He_frac, /* */ + eagle_cooling_N_temperature); /* */ + } else { + + log_10_T = + interpolation_4d(cooling->table.temperature, /* */ + /*z_index=*/0, n_H_index, He_index, u_index, /* */ + cooling->dz, d_n_H, d_He, d_u, /* */ + eagle_cooling_N_loaded_redshifts, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_He_frac, /* */ + eagle_cooling_N_temperature); /* */ + } + + if (compute_dT_du) { + + float log_10_T_high, log_10_T_low; + + /* Interpolate temperature table to return temperature for internal energy + * at grid point above current internal energy for computing dT_du used for + * calculation of dlambda_du in cooling.c (use 3D interpolation for high + * redshift table, otherwise 4D) */ + if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) { + + log_10_T_high = interpolation_3d(cooling->table.temperature, /* */ + n_H_index, He_index, u_index, /* */ + d_n_H, d_He, /*delta_u=*/1.f, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_He_frac, /* */ + eagle_cooling_N_temperature); /* */ + + } else { + + log_10_T_high = + interpolation_4d(cooling->table.temperature, /* */ + /*z_index=*/0, n_H_index, He_index, u_index, /* */ + cooling->dz, d_n_H, d_He, /*delta_u=*/1.f, /* */ + eagle_cooling_N_loaded_redshifts, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_He_frac, /* */ + eagle_cooling_N_temperature); /* */ + } + + /* Interpolate temperature table to return temperature for internal energy + * at grid point below current internal energy for computing dT_du used for + * calculation of dlambda_du in cooling.c (use 3D interpolation for high + * redshift table, otherwise 4D) */ + if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) { + + log_10_T_low = interpolation_3d(cooling->table.temperature, /* */ + n_H_index, He_index, u_index, /* */ + d_n_H, d_He, /*delta_u=*/0.f, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_He_frac, /* */ + eagle_cooling_N_temperature); /* */ + + } else { + + log_10_T_low = + interpolation_4d(cooling->table.temperature, /* */ + /*z_index=*/0, n_H_index, He_index, u_index, /* */ + cooling->dz, d_n_H, d_He, /*delta_u=*/0.f, /* */ + eagle_cooling_N_loaded_redshifts, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_He_frac, /* */ + eagle_cooling_N_temperature); /* */ + } + + /* Calculate dT/du */ + const float delta_u = exp(cooling->Therm[u_index + 1] * M_LN10) - + exp(cooling->Therm[u_index] * M_LN10); + *dT_du = + (exp(M_LN10 * log_10_T_high) - exp(M_LN10 * log_10_T_low)) / delta_u; + } + + /* Special case for temperatures below the start of the table */ + if (u_index == 0 && d_u == 0.f) { + + /* The temperature is multiplied by u / 10^T[0] + * where T[0] is the first entry in the table */ + log_10_T += log_10_u_cgs - cooling->Temp[0]; + } + + return log_10_T; +} + +/** + * @brief Computes the cooling rate corresponding to a given internal energy, + * hydrogen number density, Helium fraction, redshift and metallicity from + * all the possible channels. + * + * 1) Metal-free cooling: + * We interpolate the flattened 4D table 'H_and_He_net_heating' that is + * arranged in the following way: + * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts + * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density + * - 3rd dim: Helium fraction, length = eagle_cooling_N_He_frac + * - 4th dim: Internal energy, length = eagle_cooling_N_temperature + * + * 2) Electron abundance + * We compute the electron abundance by interpolating the flattened 4d table + * 'H_and_He_electron_abundance' that is arranged in the following way: + * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts + * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density + * - 3rd dim: Helium fraction, length = eagle_cooling_N_He_frac + * - 4th dim: Internal energy, length = eagle_cooling_N_temperature + * + * 3) Compton cooling is applied via the analytic formula. + * + * 4) Solar electron abudance + * We compute the solar electron abundance by interpolating the flattened 3d + * table 'solar_electron_abundance' that is arranged in the following way: + * - 1st dim: redshift, length = eagle_cooling_N_loaded_redshifts + * - 2nd dim: Hydrogen density, length = eagle_cooling_N_density + * - 3rd dim: Internal energy, length = eagle_cooling_N_temperature + * + * 5) Metal-line cooling + * For each tracked element we interpolate the flattened 4D table + * 'table_metals_net_heating' that is arrange in the following way: + * - 1st dim: element, length = eagle_cooling_N_metal + * - 2nd dim: redshift, length = eagle_cooling_N_loaded_redshifts + * - 3rd dim: Hydrogen density, length = eagle_cooling_N_density + * - 4th dim: Internal energy, length = eagle_cooling_N_temperature + * + * Note that this is a fake 4D interpolation as we do not interpolate + * along the 1st dimension. We just do this once per element. + * + * Since only the temperature changes when cooling a given particle, + * the redshift, hydrogen number density and helium fraction indices + * and offsets passed in. + * + * If the arguement dlambda_du is non-NULL, the routine also + * calculates derivative of cooling rate with respect to internal + * energy. + * + * If the argument element_lambda is non-NULL, the routine also + * returns the cooling rate per element in the array. + * + * @param log10_u_cgs Log base 10 of internal energy per unit mass in CGS units. + * @param redshift The current redshift + * @param n_H_cgs The Hydrogen number density in CGS units. + * @param solar_ratio Array of ratios of particle metal abundances + * to solar metal abundances + * + * @param n_H_index Particle hydrogen number density index + * @param d_n_H Particle hydrogen number density offset + * @param He_index Particle helium fraction index + * @param d_He Particle helium fraction offset + * @param cooling Cooling data structure + * + * @param dlambda_du (return) Derivative of the cooling rate with respect to u. + * @param element_lambda (return) Cooling rate from each element + * + * @return The cooling rate + */ +INLINE static double eagle_metal_cooling_rate( + double log10_u_cgs, double redshift, double n_H_cgs, + const float solar_ratio[chemistry_element_count + 2], int n_H_index, + float d_n_H, int He_index, float d_He, + const struct cooling_function_data *restrict cooling, double *dlambda_du, + double *element_lambda) { + + /* used for calculating dlambda_du */ + double temp_lambda_high = 0, temp_lambda_low = 0; + double h_plus_he_electron_abundance_high = 0; + double h_plus_he_electron_abundance_low = 0; + double solar_electron_abundance_high = 0; + double solar_electron_abundance_low = 0; + double elem_cool_low = 0, elem_cool_high = 0; + + /* We only need dT_du if dLambda_du is non-NULL */ + const int compute_dT_du = (dlambda_du != NULL) ? 1 : 0; + + /* Temperature */ + float dT_du = -1.f; + const double log_10_T = + eagle_convert_u_to_temp(log10_u_cgs, redshift, compute_dT_du, &dT_du, + n_H_index, He_index, d_n_H, d_He, cooling); + + /* Get index along temperature dimension of the tables */ + int T_index; + float d_T; + get_index_1d(cooling->Temp, eagle_cooling_N_temperature, log_10_T, &T_index, + &d_T); + + /* Difference between entries on the temperature table around u */ + const float delta_T = exp(M_LN10 * cooling->Temp[T_index + 1]) - + exp(M_LN10 * cooling->Temp[T_index]); + + /**********************/ + /* Metal-free cooling */ + /**********************/ + + double Lambda_free; + + if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) { + + /* If we're using the high redshift tables then we don't interpolate + * in redshift */ + Lambda_free = interpolation_3d(cooling->table.H_plus_He_heating, /* */ + n_H_index, He_index, T_index, /* */ + d_n_H, d_He, d_T, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_He_frac, /* */ + eagle_cooling_N_temperature); /* */ + + /* compute values at temperature gridpoints above and below input + * temperature for calculation of dlambda_du. Pass in NULL pointer for + * dlambda_du in order to skip */ + if (dlambda_du != NULL) { + temp_lambda_high = interpolation_3d( + cooling->table.H_plus_He_heating, n_H_index, He_index, T_index, d_n_h, + d_He, 1.f, cooling->N_nH, cooling->N_He, cooling->N_Temp); + temp_lambda_low = interpolation_3d( + cooling->table.H_plus_He_heating, n_H_index, He_index, T_index, d_n_h, + d_He, 0.f, cooling->N_nH, cooling->N_He, cooling->N_Temp); + } + + } else { + + /* Using normal tables, have to interpolate in redshift */ + Lambda_free = + interpolation_4d(cooling->table.H_plus_He_heating, /* */ + /*z_index=*/0, n_H_index, He_index, T_index, /* */ + cooling->dz, d_n_H, d_He, d_T, /* */ + eagle_cooling_N_loaded_redshifts, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_He_frac, /* */ + eagle_cooling_N_temperature); /* */ + + /* compute values at temperature gridpoints above and below input + * temperature for calculation of dlambda_du */ + if (dlambda_du != NULL) { + temp_lambda_high = + interpolation_4d(cooling->table.H_plus_He_heating, 0, n_H_index, + He_index, T_index, cooling->dz, d_n_h, d_He, 1.f, 2, + cooling->N_nH, cooling->N_He, cooling->N_Temp); + temp_lambda_low = + interpolation_4d(cooling->table.H_plus_He_heating, 0, n_H_index, + He_index, T_index, cooling->dz, d_n_h, d_He, 0.f, 2, + cooling->N_nH, cooling->N_He, cooling->N_Temp); + } + } + + if (dlambda_du != NULL) { + *dlambda_du += (temp_lambda_high - temp_lambda_low) / delta_T * dT_du; + } + + /* If we're testing cooling rate contributions write to array */ + if (element_lambda != NULL) { + element_lambda[0] = Lambda_free; + } + + /**********************/ + /* Electron abundance */ + /**********************/ + + double H_plus_He_electron_abundance; + + if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) { + + H_plus_He_electron_abundance = + interpolation_3d(cooling->table.H_plus_He_electron_abundance, /* */ + n_H_index, He_index, T_index, /* */ + d_n_H, d_He, d_T, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_He_frac, /* */ + eagle_cooling_N_temperature); /* */ + /* compute values at temperature gridpoints above and below input + * temperature for calculation of dlambda_du. Pass in NULL pointer for + * dlambda_du in order to skip */ + + h_plus_he_electron_abundance_high = + interpolation_3d(cooling->table.H_plus_He_electron_abundance, n_H_index, + He_index, T_index, d_n_h, d_He, 1.f, cooling->N_nH, + cooling->N_He, cooling->N_Temp); + h_plus_he_electron_abundance_low = + interpolation_3d(cooling->table.H_plus_He_electron_abundance, n_H_index, + He_index, T_index, d_n_h, d_He, 0.f, cooling->N_nH, + cooling->N_He, cooling->N_Temp); + + } else { + + H_plus_He_electron_abundance = + interpolation_4d(cooling->table.H_plus_He_electron_abundance, /* */ + /*z_index=*/0, n_H_index, He_index, T_index, /* */ + cooling->dz, d_n_H, d_He, d_T, /* */ + eagle_cooling_N_loaded_redshifts, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_He_frac, /* */ + eagle_cooling_N_temperature); /* */ + + /* compute values at temperature gridpoints above and below input + * temperature for calculation of dlambda_du */ + h_plus_he_electron_abundance_high = + interpolation_4d(cooling->table.H_plus_He_electron_abundance, 0, + n_H_index, He_index, T_index, cooling->dz, d_n_h, d_He, + 1.f, 2, cooling->N_nH, cooling->N_He, cooling->N_Temp); + h_plus_he_electron_abundance_low = + interpolation_4d(cooling->table.H_plus_He_electron_abundance, 0, + n_H_index, He_index, T_index, cooling->dz, d_n_h, d_He, + 0.f, 2, cooling->N_nH, cooling->N_He, cooling->N_Temp); + } + + /**********************/ + /* Compton cooling */ + /**********************/ + + double Lambda_Compton = 0.; + + /* Do we need to add the inverse Compton cooling? */ + /* It is *not* stored in the tables before re-ionisation */ + if ((redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) || + (redshift > cooling->H_reion_z)) { + + const double T = exp10(log_10_T); + + /* Note the minus sign */ + Lambda_Compton -= eagle_Compton_cooling_rate(cooling, redshift, n_H_cgs, T, + H_plus_He_electron_abundance); + } + + /* If we're testing cooling rate contributions write to array */ + if (element_lambda != NULL) { + element_lambda[1] = Lambda_Compton; + } + + /*******************************/ + /* Solar electron abundance */ + /*******************************/ + + double solar_electron_abundance; + + if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) { + + /* If we're using the high redshift tables then we don't interpolate + * in redshift */ + solar_electron_abundance = + interpolation_2d(cooling->table.electron_abundance, /* */ + n_H_index, T_index, /* */ + d_n_H, d_T, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_temperature); /* */ + + /* compute values at temperature gridpoints above and below input + * temperature for calculation of dlambda_du */ + if (dlambda_du != NULL) { + solar_electron_abundance_high = + interpolation_2d(cooling->table.electron_abundance, n_H_index, + T_index, d_n_h, 1.f, cooling->N_nH, cooling->N_Temp); + solar_electron_abundance_low = + interpolation_2d(cooling->table.electron_abundance, n_H_index, + T_index, d_n_h, 0.f, cooling->N_nH, cooling->N_Temp); + } + + } else { + + /* Using normal tables, have to interpolate in redshift */ + solar_electron_abundance = + interpolation_3d(cooling->table.electron_abundance, /* */ + /*z_index=*/0, n_H_index, T_index, /* */ + cooling->dz, d_n_H, d_T, /* */ + eagle_cooling_N_loaded_redshifts, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_temperature); /* */ + + /* compute values at temperature gridpoints above and below input + * temperature for calculation of dlambda_du */ + if (dlambda_du != NULL) { + solar_electron_abundance_high = interpolation_3d( + cooling->table.electron_abundance, 0, n_H_index, T_index, cooling->dz, + d_n_h, 1.f, 2, cooling->N_nH, cooling->N_Temp); + solar_electron_abundance_low = interpolation_3d( + cooling->table.electron_abundance, 0, n_H_index, T_index, cooling->dz, + d_n_h, 0.f, 2, cooling->N_nH, cooling->N_Temp); + } + + const double electron_abundance_ratio = + H_plus_He_electron_abundance / solar_electron_abundance; + + /**********************/ + /* Metal-line cooling */ + /**********************/ + + /* for each element the cooling rate is multiplied by the ratio of H, He + * electron abundance to solar electron abundance then by the ratio of the + * particle metal abundance to solar metal abundance. */ + + double lambda_metal[eagle_cooling_N_metal + 2] = {0.}; + + if (redshift > cooling->Redshifts[eagle_cooling_N_redshifts - 1]) { + + /* Loop over the metals (ignore H and He) */ + for (int elem = 2; elem < eagle_cooling_N_metal + 2; elem++) { + + if (solar_ratio[elem] > 0.) { + + /* Note that we do not interpolate along the x-axis + * (element dimension) */ + lambda_metal[elem] = + interpolation_3d_no_x(cooling->table.metal_heating, /* */ + elem - 2, n_H_index, T_index, /* */ + /*delta_elem=*/0.f, d_n_H, d_T, /* */ + eagle_cooling_N_metal, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_temperature); /* */ + + lambda_metal[elem] *= electron_abundance_ratio; + lambda_metal[elem] *= solar_ratio[elem]; + } + + /* compute values at temperature gridpoints above and below input + * temperature for calculation of dlambda_du */ + if (dlambda_du != NULL) { + elem_cool_high = interpolation_3d_no_x( + cooling->table.metal_heating, elem, n_H_index, T_index, 0.f, + d_n_h, 1.f, cooling->N_Elements, cooling->N_nH, cooling->N_Temp); + + elem_cool_low = interpolation_3d_no_x( + cooling->table.metal_heating, elem, n_H_index, T_index, 0.f, + d_n_h, 0.f, cooling->N_nH, cooling->N_Temp, cooling->N_Elements); + + *dlambda_du += (elem_cool_high * h_plus_he_electron_abundance_high / + solar_electron_abundance_high - + elem_cool_low * h_plus_he_electron_abundance_low / + solar_electron_abundance_low) / + delta_T * dT_du * solar_ratio[elem + 2]; + } + } + + } else { + + /* Loop over the metals (ignore H and He) */ + for (int elem = 2; elem < eagle_cooling_N_metal + 2; elem++) { + + if (solar_ratio[elem] > 0.) { + + /* Note that we do not interpolate along the x-axis + * (element dimension) */ + lambda_metal[elem] = interpolation_4d_no_x( + cooling->table.metal_heating, /* */ + elem - 2, /*z_index=*/0, n_H_index, T_index, /* */ + /*delta_elem=*/0.f, cooling->dz, d_n_H, d_T, /* */ + eagle_cooling_N_metal, /* */ + eagle_cooling_N_loaded_redshifts, /* */ + eagle_cooling_N_density, /* */ + eagle_cooling_N_temperature); /* */ + + lambda_metal[elem] *= electron_abundance_ratio; + lambda_metal[elem] *= solar_ratio[elem]; + } + + /* compute values at temperature gridpoints above and below input + * temperature for calculation of dlambda_du */ + if (dlambda_du != NULL) { + elem_cool_high = interpolation_4d_no_x( + cooling->table.metal_heating, elem, 0, n_H_index, T_index, 0., + cooling->dz, d_n_h, 1.f, cooling->N_Elements, 2, cooling->N_nH, + cooling->N_Temp); + + elem_cool_low = interpolation_4d_no_x( + cooling->table.metal_heating, elem, 0, n_H_index, T_index, 0., + cooling->dz, d_n_h, 0.f, cooling->N_Elements, 2, cooling->N_nH, + cooling->N_Temp); + + *dlambda_du += (elem_cool_high * h_plus_he_electron_abundance_high / + solar_electron_abundance_high - + elem_cool_low * h_plus_he_electron_abundance_low / + solar_electron_abundance_low) / + delta_T * dT_du * solar_ratio[elem + 2]; + } + } + } + + if (element_lambda != NULL) { + for (int elem = 2; elem < eagle_cooling_N_metal + 2; ++elem) { + element_lambda[elem] = lambda_metal[elem]; + } + } + + /* Sum up all the contributions */ + double Lambda_net = Lambda_free + Lambda_Compton; + for (int elem = 2; elem < eagle_cooling_N_metal + 2; ++elem) { + Lambda_net += lambda_metal[elem]; + } + + return Lambda_net; + } + + /** + * @brief Wrapper function used to calculate cooling rate and dLambda_du. + * Table indices and offsets for redshift, hydrogen number density and + * helium fraction are passed it so as to compute them only once per particle. + * + * @param log_u_cgs Natural log of internal energy per unit mass in CGS units. + * @param redshift The current redshift. + * @param n_H_cgs Hydrogen number density in CGS units. + * @param abundance_ratio Ratio of element abundance to solar. + * + * @param n_H_index Particle hydrogen number density index + * @param d_n_H Particle hydrogen number density offset + * @param He_index Particle helium fraction index + * @param d_He Particle helium fraction offset + * @param cooling #cooling_function_data structure + * + * @param dLambdaNet_du (return) Derivative of the cooling rate with respect + * to u. + * + * @return The cooling rate + */ + INLINE static double eagle_cooling_rate( + double log_u_cgs, double redshift, double n_H_cgs, + const float abundance_ratio[chemistry_element_count + 2], int n_H_index, + float d_n_H, int He_index, float d_He, + const struct cooling_function_data *restrict cooling, + double *dLambdaNet_du) { + + return eagle_metal_cooling_rate(log_u_cgs / M_LN10, redshift, n_H_cgs, + abundance_ratio, n_H_index, d_n_H, He_index, + d_He, cooling, dLambdaNet_du, + /*element_lambda=*/NULL); + } + + /** + * @brief Newton Raphson integration scheme to calculate particle cooling over + * timestep. This replaces bisection scheme used in EAGLE to minimize the + * number of array accesses. Integration defaults to bisection scheme (see + * function bisection_iter) if this function does not converge within a + * specified number of steps + * + * @param logu_init Initial guess for log(internal energy) + * @param u_ini Internal energy at beginning of hydro step + * @param n_H_index Particle hydrogen number density index + * @param d_n_H Particle hydrogen number density offset + * @param He_index Particle helium fraction index + * @param d_He Particle helium fraction offset + * @param He_reion_heat Heating due to helium reionization + * (only depends on redshift, so passed as parameter) + * @param p #part structure + * @param cosmo #cosmology structure + * @param cooling #cooling_function_data structure + * @param phys_const #phys_const data structure + * @param abundance_ratio Array of ratios of metal abundance to solar + * @param dt timestep + * @param bisection_flag Flag to identify if scheme failed to converge + */ + float newton_iter(float logu_init, double u_ini, int n_H_index, float d_n_H, + int He_index, float d_He, float He_reion_heat, + struct part *restrict p, + const struct cosmology *restrict cosmo, + const struct cooling_function_data *restrict cooling, + const struct phys_const *restrict phys_const, + const float abundance_ratio[chemistry_element_count + 2], + float dt, int *bisection_flag) { + + double logu, logu_old; + double dLambdaNet_du = 0.0, LambdaNet; + + /* table bounds */ + const float log_table_bound_high = + (cooling->Therm[eagle_cooling_N_temperature - 1] - 0.05) / M_LOG10E; + const float log_table_bound_low = (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 double n_H = hydro_get_physical_density(p, cosmo) * XH / + phys_const->const_proton_mass; + const double n_H_cgs = n_H * cooling->number_density_to_cgs; + + /* compute ratefact = n_H * n_H / rho; Might lead to round-off error: + * replaced by equivalent expression below */ + const double ratefact_cgs = n_H_cgs * XH * cooling->inv_proton_mass_cgs; + + logu_old = logu_init; + logu = logu_old; + int i = 0; + + float LambdaNet_old = 0; + LambdaNet = 0; + do /* iterate to convergence */ + { + logu_old = logu; + LambdaNet_old = LambdaNet; + LambdaNet = (He_reion_heat / (dt * ratefact_cgs)) + + eagle_cooling_rate(logu_old, cosmo->z, n_H_cgs, + abundance_ratio, n_H_index, d_n_H, + He_index, d_He, cooling, &dLambdaNet_du); + + /* Newton iteration. For details on how the cooling equation is integrated + * see documentation in theory/Cooling/ */ + logu = logu_old - (1.0 - u_ini * exp(-logu_old) - + LambdaNet * ratefact_cgs * dt * exp(-logu_old)) / + (1.0 - dLambdaNet_du * ratefact_cgs * dt); + /* Check if first step passes over equilibrium solution, if it does adjust + * next guess */ + if (i == 1 && LambdaNet_old * LambdaNet < 0) + logu = newton_log_u_guess_cgs; + + /* check whether iterations go within about 10% of the table bounds, + * if they do default to bisection method */ + if (logu > log_table_bound_high) { + i = newton_max_iterations; + break; + } else if (logu < log_table_bound_low) { + i = newton_max_iterations; + break; + } + + i++; + } while (fabs(logu - logu_old) > newton_tolerance && + i < newton_max_iterations); + if (i >= newton_max_iterations) { + /* flag to trigger bisection scheme */ + *bisection_flag = 1; + } + + return logu; + } diff --git a/src/drift.h b/src/drift.h index faa8743a86e47ab53d9880e086d5acf454ca7c38..b0b3e972995eefb57f3807034c26e2e1a6dbc200 100644 --- a/src/drift.h +++ b/src/drift.h @@ -27,7 +27,9 @@ #include "const.h" #include "debug.h" #include "dimension.h" +#include "entropy_floor.h" #include "hydro.h" +#include "hydro_properties.h" #include "part.h" #include "stars.h" @@ -71,11 +73,16 @@ __attribute__((always_inline)) INLINE static void drift_gpart( * @param dt_therm The drift time-step for thermodynamic quantities. * @param ti_old Integer start of time-step (for debugging checks). * @param ti_current Integer end of time-step (for debugging checks). + * @param cosmo The cosmological model. + * @param hydro_props The properties of the hydro scheme. + * @param floor The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void drift_part( struct part *restrict p, struct xpart *restrict xp, double dt_drift, double dt_kick_hydro, double dt_kick_grav, double dt_therm, - integertime_t ti_old, integertime_t ti_current) { + integertime_t ti_old, integertime_t ti_current, + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor) { #ifdef SWIFT_DEBUG_CHECKS if (p->ti_drift != ti_old) @@ -102,7 +109,7 @@ __attribute__((always_inline)) INLINE static void drift_part( p->v[2] += xp->a_grav[2] * dt_kick_grav; /* Predict the values of the extra fields */ - hydro_predict_extra(p, xp, dt_drift, dt_therm); + hydro_predict_extra(p, xp, dt_drift, dt_therm, cosmo, hydro_props, floor); /* Compute offsets since last cell construction */ for (int k = 0; k < 3; k++) { diff --git a/src/entropy_floor.h b/src/entropy_floor.h index 9f2e97ccc815bee1087884d060492ff3715f1c6f..771bd098bc0bb8390da64bd7331404ad982f7d3b 100644 --- a/src/entropy_floor.h +++ b/src/entropy_floor.h @@ -27,10 +27,16 @@ /* Config parameters. */ #include "../config.h" +/* Local includes */ #include "common_io.h" #include "error.h" #include "inline.h" +/* Pre-declarations to avoid cyclic inclusions */ +static INLINE float hydro_get_physical_density(const struct part *restrict p, + const struct cosmology *cosmo); +static INLINE float hydro_get_comoving_density(const struct part *restrict p); + /* Import the right entropy floor definition */ #if defined(ENTROPY_FLOOR_NONE) #include "./entropy_floor/none/entropy_floor.h" diff --git a/src/entropy_floor/EAGLE/entropy_floor.h b/src/entropy_floor/EAGLE/entropy_floor.h index 5a53b3b39c9df3aa624b3fd93e82bed9d221204c..c52f8d111302718c8c7fd7e410d97fe9ed3031de 100644 --- a/src/entropy_floor/EAGLE/entropy_floor.h +++ b/src/entropy_floor/EAGLE/entropy_floor.h @@ -21,6 +21,7 @@ #include "adiabatic_index.h" #include "cosmology.h" +#include "equation_of_state.h" #include "hydro.h" #include "hydro_properties.h" #include "parser.h" diff --git a/src/entropy_floor/none/entropy_floor.h b/src/entropy_floor/none/entropy_floor.h index 6ce5319c8c76cf8f9e6dc1716f383452d7cce7d6..8c89e08ab2ecf2af5eaf0d6a4a8527c63d700e6b 100644 --- a/src/entropy_floor/none/entropy_floor.h +++ b/src/entropy_floor/none/entropy_floor.h @@ -64,6 +64,7 @@ static INLINE float entropy_floor( static INLINE float entropy_floor_temperature( const struct part *p, const struct cosmology *cosmo, const struct entropy_floor_properties *props) { + return 0.f; } diff --git a/src/hydro/AnarchyPU/hydro.h b/src/hydro/AnarchyPU/hydro.h index 0ac52165ec4b6dc5c193cf3d22d20458d2e643d3..6d340f2e976fd9b79fc968fc3d65ae26abd731c8 100644 --- a/src/hydro/AnarchyPU/hydro.h +++ b/src/hydro/AnarchyPU/hydro.h @@ -40,6 +40,7 @@ #include "approx_math.h" #include "cosmology.h" #include "dimension.h" +#include "entropy_floor.h" #include "equation_of_state.h" #include "hydro_properties.h" #include "hydro_space.h" @@ -604,6 +605,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient( */ __attribute__((always_inline)) INLINE static void hydro_reset_gradient( struct part *restrict p) { + p->viscosity.v_sig = 2.f * p->force.soundspeed; } @@ -628,7 +630,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient( /* Include the extra factors in the del^2 u */ - p->diffusion.laplace_u *= 2 * h_inv_dim_plus_one; + p->diffusion.laplace_u *= 2.f * h_inv_dim_plus_one; } /** @@ -696,16 +698,17 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* We use in this function that h is the radius of support */ const float kernel_support_physical = p->h * cosmo->a * kernel_gamma; + const float kernel_support_physical_inv = 1.f / kernel_support_physical; const float v_sig_physical = p->viscosity.v_sig * cosmo->a_factor_sound_speed; const float soundspeed_physical = hydro_get_physical_soundspeed(p, cosmo); const float sound_crossing_time_inverse = - soundspeed_physical / kernel_support_physical; + soundspeed_physical * kernel_support_physical_inv; /* Construct time differential of div.v implicitly following the ANARCHY spec */ - float div_v_dt = + const float div_v_dt = dt_alpha == 0.f ? 0.f : (p->viscosity.div_v - p->viscosity.div_v_previous_step) / dt_alpha; @@ -734,9 +737,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( hydro_props->viscosity.length); } - if (p->viscosity.alpha < hydro_props->viscosity.alpha_min) { - p->viscosity.alpha = hydro_props->viscosity.alpha_min; - } + /* Check that we did not hit the minimum */ + p->viscosity.alpha = + max(p->viscosity.alpha, hydro_props->viscosity.alpha_min); /* Set our old div_v to the one for the next loop */ p->viscosity.div_v_previous_step = p->viscosity.div_v; @@ -744,18 +747,20 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( /* Now for the diffusive alpha */ const float diffusion_timescale_physical_inverse = - v_sig_physical / kernel_support_physical; + v_sig_physical * kernel_support_physical_inv; + + const float sqrt_u_inv = 1.f / sqrtf(p->u); - const float sqrt_u = sqrtf(p->u); /* Calculate initial value of alpha dt before bounding */ /* Evolution term: following Schaller+ 2015. This is made up of several - cosmology factors: physical h, sound speed from u / sqrt(u), and the - 1 / a^2 coming from the laplace operator. */ + cosmology factors: physical smoothing length, sound speed from laplace(u) / + sqrt(u), and the 1 / a^2 coming from the laplace operator. */ float alpha_diff_dt = hydro_props->diffusion.beta * kernel_support_physical * - p->diffusion.laplace_u * cosmo->a_factor_sound_speed / - (sqrt_u * cosmo->a * cosmo->a); + p->diffusion.laplace_u * cosmo->a_factor_sound_speed * + sqrt_u_inv * cosmo->a2_inv; + /* Decay term: not documented in Schaller+ 2015 but was present - * in the original EAGLE code and in Schaye+ 2015 */ + * in the original EAGLE code and in appendix of Schaye+ 2015 */ alpha_diff_dt -= (p->diffusion.alpha - hydro_props->diffusion.alpha_min) * diffusion_timescale_physical_inverse; @@ -763,11 +768,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force( new_diffusion_alpha += alpha_diff_dt * dt_alpha; /* Consistency checks to ensure min < alpha < max */ - if (new_diffusion_alpha > hydro_props->diffusion.alpha_max) { - new_diffusion_alpha = hydro_props->diffusion.alpha_max; - } else if (new_diffusion_alpha < hydro_props->diffusion.alpha_min) { - new_diffusion_alpha = hydro_props->diffusion.alpha_min; - } + new_diffusion_alpha = + min(new_diffusion_alpha, hydro_props->diffusion.alpha_max); + new_diffusion_alpha = + max(new_diffusion_alpha, hydro_props->diffusion.alpha_min); p->diffusion.alpha = new_diffusion_alpha; } @@ -810,6 +814,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( /* Re-set the entropy */ p->u = xp->u_full; + + /* Compute the sound speed */ + const float soundspeed = hydro_get_comoving_soundspeed(p); + + p->force.soundspeed = soundspeed; } /** @@ -825,10 +834,29 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param xp The extended data of the particle. * @param dt_drift The drift time-step for positions. * @param dt_therm The drift time-step for thermal quantities. + * @param cosmo The cosmological model. + * @param hydro_props The properties of the hydro scheme. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, - float dt_therm) { + float dt_therm, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Predict the internal energy */ + p->u += p->u_dt * dt_therm; + + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = + gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, floor_A); + + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; + + p->u = max(p->u, floor_u); + p->u = max(p->u, min_u); const float h_inv = 1.f / p->h; @@ -852,9 +880,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( p->pressure_bar *= expf_exact; } - /* Predict the internal energy */ - p->u += p->u_dt * dt_therm; - /* Compute the new sound speed */ const float soundspeed = hydro_get_comoving_soundspeed(p); @@ -893,30 +918,35 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @param dt_kick_corr The time-step for this kick (for gravity corrections). * @param cosmo The cosmological model. * @param hydro_props The constants used in the scheme + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( struct part *restrict p, struct xpart *restrict xp, float dt_therm, float dt_grav, float dt_hydro, float dt_kick_corr, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Integrate the internal energy forward in time */ + const float delta_u = p->u_dt * dt_therm; /* 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) { - p->u_dt = -0.5f * xp->u_full / dt_therm; - } - xp->u_full += p->u_dt * dt_therm; + xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full); - /* 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; - } + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = + gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, floor_A); - /* Compute the sound speed */ - const float soundspeed = hydro_get_comoving_soundspeed(p); + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; - p->force.soundspeed = soundspeed; + /* Take highest of both limits */ + const float energy_min = max(min_u, floor_u); + + if (xp->u_full < energy_min) { + xp->u_full = energy_min; + p->u_dt = 0.f; + } } /** diff --git a/src/hydro/AnarchyPU/hydro_iact.h b/src/hydro/AnarchyPU/hydro_iact.h index d091ebebcf4f165db006ca58667e943ddbaf8599..f3a8734fbf037b43d509750e018d51c9b01933b4 100644 --- a/src/hydro/AnarchyPU/hydro_iact.h +++ b/src/hydro/AnarchyPU/hydro_iact.h @@ -364,7 +364,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force( /* Includes the hubble flow term; not used for du/dt */ const float dvdr_Hubble = dvdr + a2_Hubble * r2; - /* Are the part*icles moving towards each others ? */ + /* Are the particles moving towards each others ? */ const float omega_ij = min(dvdr_Hubble, 0.f); const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ @@ -495,7 +495,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force( /* Includes the hubble flow term; not used for du/dt */ const float dvdr_Hubble = dvdr + a2_Hubble * r2; - /* Are the part*icles moving towards each others ? */ + /* Are the particles moving towards each others ? */ const float omega_ij = min(dvdr_Hubble, 0.f); const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */ @@ -600,4 +600,4 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_limiter( } } -#endif /* SWIFT_MINIMAL_HYDRO_IACT_H */ +#endif /* SWIFT_ANARCHY_PU_HYDRO_IACT_H */ diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h index 175265dd4707be19e995e2ca6912629876ca2ca2..c575fc54f0214c85a37df85586cc97574c2b0e3d 100644 --- a/src/hydro/Default/hydro.h +++ b/src/hydro/Default/hydro.h @@ -23,6 +23,7 @@ #include "approx_math.h" #include "cosmology.h" #include "dimension.h" +#include "entropy_floor.h" #include "equation_of_state.h" #include "hydro_properties.h" #include "hydro_space.h" @@ -617,11 +618,29 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param xp The extended data of the particle * @param dt_drift The drift time-step for positions. * @param dt_therm The drift time-step for thermal quantities. + * @param cosmo The cosmological model. + * @param hydro_props The properties of the hydro scheme. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, struct xpart *restrict xp, float dt_drift, - float dt_therm) { - float u, w; + float dt_therm, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Predict the internal energy */ + p->u += p->force.u_dt * dt_therm; + + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = + gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, floor_A); + + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; + + p->u = max(p->u, floor_u); + p->u = max(p->u, min_u); const float h_inv = 1.f / p->h; @@ -639,15 +658,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( else p->rho *= expf(w2); - /* Predict internal energy */ - w = p->force.u_dt / p->u * dt_therm; - if (fabsf(w) < 0.2f) - u = p->u *= approx_expf(w); - else - u = p->u *= expf(w); - /* Predict gradient term */ - p->force.P_over_rho2 = u * hydro_gamma_minus_one / (p->rho * xp->omega); + p->force.P_over_rho2 = p->u * hydro_gamma_minus_one / (p->rho * xp->omega); } /** @@ -673,7 +685,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( __attribute__((always_inline)) INLINE static void hydro_kick_extra( struct part *restrict p, struct xpart *restrict xp, float dt_therm, float dt_grav, float dt_hydro, float dt_kick_corr, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) {} + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) {} /** * @brief Converts hydro quantity of a particle at the start of a run diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h index e6df246dd10e84822b5c1cb806c6ad9840089f20..e89396182c2acfe6fd080435b27cdfddf62d43fb 100644 --- a/src/hydro/Gadget2/hydro.h +++ b/src/hydro/Gadget2/hydro.h @@ -35,6 +35,7 @@ #include "approx_math.h" #include "cosmology.h" #include "dimension.h" +#include "entropy_floor.h" #include "equation_of_state.h" #include "hydro_properties.h" #include "hydro_space.h" @@ -638,6 +639,20 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( /* Re-set the entropy */ p->entropy = xp->entropy_full; + + /* Re-compute the pressure */ + const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); + + /* Compute the new sound speed */ + const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + + /* Divide the pressure by the density squared to get the SPH term */ + const float rho_inv = 1.f / p->rho; + const float P_over_rho2 = pressure * rho_inv * rho_inv; + + /* Update variables */ + p->force.soundspeed = soundspeed; + p->force.P_over_rho2 = P_over_rho2; } /** @@ -647,10 +662,29 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param xp The extended data of the particle * @param dt_drift The drift time-step for positions. * @param dt_therm The drift time-step for thermal quantities. + * @param cosmo The cosmological model. + * @param hydro_props The properties of the hydro scheme. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, - float dt_therm) { + float dt_therm, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Predict the entropy */ + p->entropy += p->entropy_dt * dt_therm; + + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; + const float min_A = + gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, min_u); + + p->entropy = max(p->entropy, floor_A); + p->entropy = max(p->entropy, min_A); const float h_inv = 1.f / p->h; @@ -668,16 +702,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( else p->rho *= expf(w2); - /* Predict the entropy */ -#ifdef SWIFT_DEBUG_CHECKS - if (p->entropy + p->entropy_dt * dt_therm <= 0) - error( - "Negative entropy for particle id %llu old entropy %.e d_entropy %.e " - "entropy_dt %.e dt therm %.e", - p->id, p->entropy, p->entropy_dt * dt_therm, p->entropy_dt, dt_therm); -#endif - p->entropy += p->entropy_dt * dt_therm; - /* Re-compute the pressure */ const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); @@ -720,42 +744,37 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @param dt_hydro The time-step for this kick (for hydro forces) * @param dt_kick_corr The time-step for this kick (for correction of the kick) * @param cosmo The cosmological model. - * @param hydro_props The constants used in the scheme + * @param hydro_props The constants used in the scheme. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( struct part *restrict p, struct xpart *restrict xp, float dt_therm, float dt_grav, float dt_hydro, float dt_kick_corr, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_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) { - p->entropy_dt = -0.5f * xp->entropy_full / dt_therm; - } - xp->entropy_full += p->entropy_dt * dt_therm; + /* Integrate the entropy forward in time */ + const float delta_entropy = p->entropy_dt * dt_therm; - /* Apply the minimal energy limit */ - const float physical_density = p->rho * cosmo->a3_inv; - const float min_physical_energy = hydro_props->minimal_internal_energy; - const float min_physical_entropy = - gas_entropy_from_internal_energy(physical_density, min_physical_energy); - const float min_comoving_entropy = min_physical_entropy; /* A' = A */ - if (xp->entropy_full < min_comoving_entropy) { - xp->entropy_full = min_comoving_entropy; - p->entropy_dt = 0.f; - } + /* Do not decrease the entropy by more than a factor of 2 */ + xp->entropy_full = + max(xp->entropy_full + delta_entropy, 0.5f * xp->entropy_full); - /* Compute the pressure */ - const float pressure = gas_pressure_from_entropy(p->rho, xp->entropy_full); + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); - /* Compute the new sound speed */ - const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; + const float min_A = + gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, min_u); - /* Divide the pressure by the density squared to get the SPH term */ - const float rho_inv = 1.f / p->rho; - const float P_over_rho2 = pressure * rho_inv * rho_inv; + /* Take highest of both limits */ + const float entropy_min = max(min_A, floor_A); - p->force.soundspeed = soundspeed; - p->force.P_over_rho2 = P_over_rho2; + if (xp->entropy_full < entropy_min) { + xp->entropy_full = entropy_min; + p->entropy_dt = 0.f; + } } /** diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h index 3456e6f9e0a4c1715951030a8b0a10c081d746bf..c744a189a4a41720d7f0ee1f856fafd43e469033 100644 --- a/src/hydro/GizmoMFM/hydro.h +++ b/src/hydro/GizmoMFM/hydro.h @@ -23,6 +23,7 @@ #include "adiabatic_index.h" #include "approx_math.h" #include "cosmology.h" +#include "entropy_floor.h" #include "equation_of_state.h" #include "hydro_gradients.h" #include "hydro_properties.h" @@ -498,7 +499,9 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * @param xp The extended data of this particle. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part* restrict p, const struct xpart* restrict xp) {} + struct part* restrict p, const struct xpart* restrict xp) { + // MATTHIEU: Do we need something here? +} /** * @brief Converts the hydrodynamic variables from the initial condition file to @@ -524,9 +527,14 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( * @param xp The extended particle data to act upon. * @param dt_drift The drift time-step for positions. * @param dt_therm The drift time-step for thermal quantities. + * @param cosmo The cosmological model. + * @param hydro_props The properties of the hydro scheme. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( - struct part* p, struct xpart* xp, float dt_drift, float dt_therm) { + struct part* p, struct xpart* xp, float dt_drift, float dt_therm, + const struct cosmology* cosmo, const struct hydro_props* hydro_props, + const struct entropy_floor_properties* floor_props) { const float h_inv = 1.0f / p->h; @@ -565,6 +573,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( #endif } + // MATTHIEU: Apply the entropy floor here. + #ifdef SWIFT_DEBUG_CHECKS if (p->h <= 0.0f) { error("Zero or negative smoothing length (%g)!", p->h); @@ -611,11 +621,13 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @param dt_kick_corr Gravity correction time-step @f$adt@f$. * @param cosmo Cosmology. * @param hydro_props Additional hydro properties. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( struct part* p, struct xpart* xp, float dt_therm, float dt_grav, float dt_hydro, float dt_kick_corr, const struct cosmology* cosmo, - const struct hydro_props* hydro_props) { + const struct hydro_props* hydro_props, + const struct entropy_floor_properties* floor_props) { float a_grav[3]; @@ -646,6 +658,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( p->flux.energy = 0.0f; } + // MATTHIEU: Apply the entropy floor here. + gizmo_check_physical_quantities( "mass", "energy", p->conserved.mass, p->conserved.momentum[0], p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.energy); diff --git a/src/hydro/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h index 526d84289554f7e2f4cd183bbf82484f0804be28..1f26f29fa5f3f9d496c73116c7608f788f93d605 100644 --- a/src/hydro/GizmoMFV/hydro.h +++ b/src/hydro/GizmoMFV/hydro.h @@ -23,6 +23,7 @@ #include "adiabatic_index.h" #include "approx_math.h" #include "cosmology.h" +#include "entropy_floor.h" #include "equation_of_state.h" #include "hydro_gradients.h" #include "hydro_properties.h" @@ -540,7 +541,9 @@ __attribute__((always_inline)) INLINE static void hydro_reset_gradient( * @param xp The extended data of this particle. */ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( - struct part* restrict p, const struct xpart* restrict xp) {} + struct part* restrict p, const struct xpart* restrict xp) { + // MATTHIEU: Apply the entropy floor here. +} /** * @brief Converts the hydrodynamic variables from the initial condition file to @@ -568,7 +571,9 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( * @param dt_therm The drift time-step for thermal quantities. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( - struct part* p, struct xpart* xp, float dt_drift, float dt_therm) { + struct part* p, struct xpart* xp, float dt_drift, float dt_therm, + const struct cosmology* cosmo, const struct hydro_props* hydro_props, + const struct entropy_floor_properties* floor_props) { #ifdef GIZMO_LLOYD_ITERATION return; @@ -620,6 +625,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( #endif } + // MATTHIEU: Apply the entropy floor here. + /* we use a sneaky way to get the gravitational contribution to the velocity update */ p->primitives.v[0] += p->v[0] - xp->v_full[0]; @@ -671,11 +678,13 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @param dt_kick_corr Gravity correction time-step @f$adt@f$. * @param cosmo Cosmology. * @param hydro_props Additional hydro properties. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( struct part* p, struct xpart* xp, float dt_therm, float dt_grav, float dt_hydro, float dt_kick_corr, const struct cosmology* cosmo, - const struct hydro_props* hydro_props) { + const struct hydro_props* hydro_props, + const struct entropy_floor_properties* floor_props) { float a_grav[3]; @@ -741,6 +750,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( p->conserved.flux.energy = 0.f; } + // MATTHIEU: Apply the entropy floor here. + gizmo_check_physical_quantities( "mass", "energy", p->conserved.mass, p->conserved.momentum[0], p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.energy); diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h index e2fd0069524de32c49cbc3cb46553b18928d14bf..6169bcdfef72e261ec93fe14ca05d2c307fcfb0e 100644 --- a/src/hydro/Minimal/hydro.h +++ b/src/hydro/Minimal/hydro.h @@ -37,6 +37,7 @@ #include "approx_math.h" #include "cosmology.h" #include "dimension.h" +#include "entropy_floor.h" #include "equation_of_state.h" #include "hydro_properties.h" #include "hydro_space.h" @@ -628,6 +629,16 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( /* Re-set the entropy */ p->u = xp->u_full; + + /* Re-compute the pressure */ + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + + /* Compute the new sound speed */ + const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + + /* Update variables */ + p->force.pressure = pressure; + p->force.soundspeed = soundspeed; } /** @@ -643,10 +654,29 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param xp The extended data of the particle. * @param dt_drift The drift time-step for positions. * @param dt_therm The drift time-step for thermal quantities. + * @param cosmo The cosmological model. + * @param hydro_props The properties of the hydro scheme. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, - float dt_therm) { + float dt_therm, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Predict the internal energy */ + p->u += p->u_dt * dt_therm; + + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = + gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, floor_A); + + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; + + p->u = max(p->u, floor_u); + p->u = max(p->u, min_u); const float h_inv = 1.f / p->h; @@ -664,9 +694,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( else p->rho *= expf(w2); - /* Predict the internal energy */ - p->u += p->u_dt * dt_therm; - /* Compute the new pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); @@ -708,36 +735,36 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @param dt_hydro The time-step for this kick (for hydro quantities). * @param dt_kick_corr The time-step for this kick (for gravity corrections). * @param cosmo The cosmological model. - * @param hydro_props The constants used in the scheme + * @param hydro_props The constants used in the scheme. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( struct part *restrict p, struct xpart *restrict xp, float dt_therm, float dt_grav, float dt_hydro, float dt_kick_corr, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Integrate the internal energy forward in time */ + const float delta_u = p->u_dt * dt_therm; /* 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) { - p->u_dt = -0.5f * xp->u_full / dt_therm; - } - xp->u_full += p->u_dt * dt_therm; + xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full); - /* Apply the minimal energy limit */ - const float min_comoving_energy = - hydro_props->minimal_internal_energy / cosmo->a_factor_internal_energy; - if (xp->u_full < min_comoving_energy) { - xp->u_full = min_comoving_energy; - p->u_dt = 0.f; - } + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = + gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, floor_A); - /* Compute the pressure */ - const float pressure = gas_pressure_from_internal_energy(p->rho, xp->u_full); + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; - /* Compute the sound speed */ - const float soundspeed = - gas_soundspeed_from_internal_energy(p->rho, xp->u_full); + /* Take highest of both limits */ + const float energy_min = max(min_u, floor_u); - p->force.pressure = pressure; - p->force.soundspeed = soundspeed; + if (xp->u_full < energy_min) { + xp->u_full = energy_min; + p->u_dt = 0.f; + } } /** diff --git a/src/hydro/Planetary/hydro.h b/src/hydro/Planetary/hydro.h index e03b307bca754eb3a809613b20ebf3b1aae67a9f..ba9228289fd43e6cf130b58e968de63e4fee3147 100644 --- a/src/hydro/Planetary/hydro.h +++ b/src/hydro/Planetary/hydro.h @@ -39,6 +39,7 @@ #include "cosmology.h" #include "debug.h" #include "dimension.h" +#include "entropy_floor.h" #include "equation_of_state.h" #include "hydro_properties.h" #include "hydro_space.h" @@ -638,6 +639,17 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( /* Re-set the entropy */ p->u = xp->u_full; + + /* Compute the pressure */ + const float pressure = + gas_pressure_from_internal_energy(p->rho, xp->u_full, p->mat_id); + + /* Compute the sound speed */ + const float soundspeed = + gas_soundspeed_from_internal_energy(p->rho, p->u, p->mat_id); + + p->force.pressure = pressure; + p->force.soundspeed = soundspeed; } /** @@ -653,10 +665,23 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param xp The extended data of the particle. * @param dt_drift The drift time-step for positions. * @param dt_therm The drift time-step for thermal quantities. + * @param cosmo The cosmological model. + * @param hydro_props The constants used in the scheme + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, - float dt_therm) { + float dt_therm, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Predict the internal energy */ + p->u += p->u_dt * dt_therm; + + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; + + p->u = max(p->u, min_u); const float h_inv = 1.f / p->h; @@ -674,9 +699,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( else p->rho *= expf(w2); - /* Predict the internal energy */ - p->u += p->u_dt * dt_therm; - /* Compute the new pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u, p->mat_id); @@ -718,36 +740,27 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @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 + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( struct part *restrict p, struct xpart *restrict xp, float dt_therm, float dt_grav, float dt_hydro, float dt_kick_corr, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Integrate the internal energy forward in time */ + const float delta_u = p->u_dt * dt_therm; /* 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) { - p->u_dt = -0.5f * xp->u_full / dt_therm; - } - xp->u_full += p->u_dt * dt_therm; + xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full); + + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; - /* 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; + if (xp->u_full < min_u) { + xp->u_full = min_u; p->u_dt = 0.f; } - - /* Compute the pressure */ - const float pressure = - gas_pressure_from_internal_energy(p->rho, xp->u_full, p->mat_id); - - /* Compute the sound speed */ - const float soundspeed = - gas_soundspeed_from_internal_energy(p->rho, p->u, p->mat_id); - - p->force.pressure = pressure; - p->force.soundspeed = soundspeed; } /** diff --git a/src/hydro/PressureEnergy/hydro.h b/src/hydro/PressureEnergy/hydro.h index 3fe3c805fd4df4b495b102f9061b1e0a507e84e9..b342120f16e2fb50012e211dd49eacfa63690f45 100644 --- a/src/hydro/PressureEnergy/hydro.h +++ b/src/hydro/PressureEnergy/hydro.h @@ -40,6 +40,7 @@ #include "approx_math.h" #include "cosmology.h" #include "dimension.h" +#include "entropy_floor.h" #include "equation_of_state.h" #include "hydro_properties.h" #include "hydro_space.h" @@ -660,6 +661,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( /* Re-set the entropy */ p->u = xp->u_full; + + /* Compute the sound speed */ + const float soundspeed = hydro_get_comoving_soundspeed(p); + + p->force.soundspeed = soundspeed; } /** @@ -675,10 +681,29 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param xp The extended data of the particle. * @param dt_drift The drift time-step for positions. * @param dt_therm The drift time-step for thermal quantities. + * @param cosmo The cosmological model. + * @param hydro_props The properties of the hydro scheme. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, - float dt_therm) { + float dt_therm, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Predict the internal energy */ + p->u += p->u_dt * dt_therm; + + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = + gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, floor_A); + + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; + + p->u = max(p->u, floor_u); + p->u = max(p->u, min_u); const float h_inv = 1.f / p->h; @@ -702,9 +727,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( p->pressure_bar *= expf_exact; } - /* Predict the internal energy */ - p->u += p->u_dt * dt_therm; - /* Compute the new sound speed */ const float soundspeed = hydro_get_comoving_soundspeed(p); @@ -743,30 +765,35 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @param dt_kick_corr The time-step for this kick (for gravity corrections). * @param cosmo The cosmological model. * @param hydro_props The constants used in the scheme + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( struct part *restrict p, struct xpart *restrict xp, float dt_therm, float dt_grav, float dt_hydro, float dt_kick_corr, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Integrate the internal energy forward in time */ + const float delta_u = p->u_dt * dt_therm; /* 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) { - p->u_dt = -0.5f * xp->u_full / dt_therm; - } - xp->u_full += p->u_dt * dt_therm; + xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full); - /* 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; - } + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = + gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, floor_A); - /* Compute the sound speed */ - const float soundspeed = hydro_get_comoving_soundspeed(p); + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; - p->force.soundspeed = soundspeed; + /* Take highest of both limits */ + const float energy_min = max(min_u, floor_u); + + if (xp->u_full < energy_min) { + xp->u_full = energy_min; + p->u_dt = 0.f; + } } /** diff --git a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h index be6a789ebdcd664a99a95b83f70167da04bd7534..de5ebafb2f51cc3739a806bbbbbfc9e20e205bcc 100644 --- a/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h +++ b/src/hydro/PressureEnergyMorrisMonaghanAV/hydro.h @@ -41,6 +41,7 @@ #include "approx_math.h" #include "cosmology.h" #include "dimension.h" +#include "entropy_floor.h" #include "equation_of_state.h" #include "hydro_properties.h" #include "hydro_space.h" @@ -674,6 +675,11 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( /* Re-set the entropy */ p->u = xp->u_full; + + /* Compute the sound speed */ + const float soundspeed = hydro_get_comoving_soundspeed(p); + + p->force.soundspeed = soundspeed; } /** @@ -689,10 +695,29 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param xp The extended data of the particle. * @param dt_drift The drift time-step for positions. * @param dt_therm The drift time-step for thermal quantities. + * @param cosmo The cosmological model. + * @param hydro_props The properties of the hydro scheme. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, - float dt_therm) { + float dt_therm, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Predict the internal energy */ + p->u += p->u_dt * dt_therm; + + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = + gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, floor_A); + + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; + + p->u = max(p->u, floor_u); + p->u = max(p->u, min_u); const float h_inv = 1.f / p->h; @@ -716,9 +741,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( p->pressure_bar *= expf_exact; } - /* Predict the internal energy */ - p->u += p->u_dt * dt_therm; - /* Compute the new sound speed */ const float soundspeed = hydro_get_comoving_soundspeed(p); @@ -757,30 +779,35 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @param dt_kick_corr The time-step for this kick (for gravity corrections). * @param cosmo The cosmological model. * @param hydro_props The constants used in the scheme + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( struct part *restrict p, struct xpart *restrict xp, float dt_therm, float dt_grav, float dt_hydro, float dt_kick_corr, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Integrate the internal energy forward in time */ + const float delta_u = p->u_dt * dt_therm; /* 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) { - p->u_dt = -0.5f * xp->u_full / dt_therm; - } - xp->u_full += p->u_dt * dt_therm; + xp->u_full = max(xp->u_full + delta_u, 0.5f * xp->u_full); - /* 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; - } + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + const float floor_u = + gas_internal_energy_from_entropy(p->rho * cosmo->a3_inv, floor_A); - /* Compute the sound speed */ - const float soundspeed = hydro_get_comoving_soundspeed(p); + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; - p->force.soundspeed = soundspeed; + /* Take highest of both limits */ + const float energy_min = max(min_u, floor_u); + + if (xp->u_full < energy_min) { + xp->u_full = energy_min; + p->u_dt = 0.f; + } } /** diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h index 578bebf635149fe29299800c2a4854c0158d3de9..8fdd201cf29ce43e5c2207dd0b219d08e9e34515 100644 --- a/src/hydro/PressureEntropy/hydro.h +++ b/src/hydro/PressureEntropy/hydro.h @@ -35,6 +35,7 @@ #include "approx_math.h" #include "cosmology.h" #include "dimension.h" +#include "entropy_floor.h" #include "equation_of_state.h" #include "hydro_properties.h" #include "hydro_space.h" @@ -614,6 +615,21 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( /* Re-set the entropy */ p->entropy = xp->entropy_full; + + /* Re-compute the pressure */ + const float pressure = gas_pressure_from_entropy(p->rho, p->entropy); + + /* Compute the new sound speed */ + const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); + + /* Divide the pressure by the density squared to get the SPH term */ + const float rho_inv = 1.f / p->rho; + const float P_over_rho2 = pressure * rho_inv * rho_inv; + + /* Update variables */ + p->force.soundspeed = soundspeed; + p->force.P_over_rho2 = P_over_rho2; + p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy); } /** @@ -623,10 +639,29 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( * @param xp The extended data of the particle * @param dt_drift The drift time-step for positions. * @param dt_therm The drift time-step for thermal quantities. + * @param cosmo The cosmological model. + * @param hydro_props The properties of the hydro scheme. + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( struct part *restrict p, const struct xpart *restrict xp, float dt_drift, - float dt_therm) { + float dt_therm, const struct cosmology *cosmo, + const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_props) { + + /* Predict the entropy */ + p->entropy += p->entropy_dt * dt_therm; + + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); + + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; + const float min_A = + gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, min_u); + + p->entropy = max(p->entropy, floor_A); + p->entropy = max(p->entropy, min_A); const float h_inv = 1.f / p->h; @@ -647,9 +682,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( p->rho_bar *= expf(w2); } - /* Predict the entropy */ - p->entropy += p->entropy_dt * dt_therm; - /* Compute the pressure */ const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy); @@ -691,42 +723,36 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @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 + * @param floor_props The properties of the entropy floor. */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( struct part *restrict p, struct xpart *restrict xp, float dt_therm, float dt_grav, float dt_hydro, float dt_kick_corr, - const struct cosmology *cosmo, const struct hydro_props *hydro_props) { + const struct cosmology *cosmo, const struct hydro_props *hydro_props, + const struct entropy_floor_properties *floor_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) { - 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_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; - } + /* Integrate the entropy forward in time */ + const float delta_entropy = p->entropy_dt * dt_therm; - /* Compute the pressure */ - const float pressure = - gas_pressure_from_entropy(p->rho_bar, xp->entropy_full); + /* Do not decrease the entropy by more than a factor of 2 */ + xp->entropy_full = + max(xp->entropy_full + delta_entropy, 0.5f * xp->entropy_full); - /* Compute the new sound speed */ - const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure); + /* Check against entropy floor */ + const float floor_A = entropy_floor(p, cosmo, floor_props); - /* Divide the pressure by the density squared to get the SPH term */ - const float rho_bar_inv = 1.f / p->rho_bar; - const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv; + /* Check against absolute minimum */ + const float min_u = hydro_props->minimal_internal_energy; + const float min_A = + gas_entropy_from_internal_energy(p->rho * cosmo->a3_inv, min_u); - p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy); - p->force.soundspeed = soundspeed; - p->force.P_over_rho2 = P_over_rho2; + /* Take highest of both limits */ + const float entropy_min = max(min_A, floor_A); + + if (xp->entropy_full < entropy_min) { + xp->entropy_full = entropy_min; + p->entropy_dt = 0.f; + } } /** diff --git a/src/hydro/Shadowswift/hydro.h b/src/hydro/Shadowswift/hydro.h index 9ac89fa02a3663e521f54b4df0f15884b041fad6..ed1b000199f3641ff458c05e6e5268a8b183c468 100644 --- a/src/hydro/Shadowswift/hydro.h +++ b/src/hydro/Shadowswift/hydro.h @@ -23,6 +23,7 @@ #include "adiabatic_index.h" #include "approx_math.h" #include "cosmology.h" +#include "entropy_floor.h" #include "equation_of_state.h" #include "hydro_gradients.h" #include "hydro_properties.h" @@ -431,7 +432,9 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( * @param dt The drift time-step. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( - struct part* p, struct xpart* xp, float dt_drift, float dt_therm) {} + struct part* p, struct xpart* xp, float dt_drift, float dt_therm, + const struct cosmology* cosmo, const struct hydro_props* hydro_props, + const struct entropy_floor_properties* floor_props) {} /** * @brief Set the particle acceleration after the flux loop. @@ -453,7 +456,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( __attribute__((always_inline)) INLINE static void hydro_kick_extra( struct part* p, struct xpart* xp, float dt, float dt_grav, float dt_hydro, float dt_kick_corr, const struct cosmology* cosmo, - const struct hydro_props* hydro_props) { + const struct hydro_props* hydro_props, + const struct entropy_floor_properties* floor_props) { /* Update the conserved variables. We do this here and not in the kick, since we need the updated variables below. */ diff --git a/src/kick.h b/src/kick.h index f2085bf1f427cf5f15ed0e8791ad1923f0b22bed..404f7f31170dcec4728b8fbc52717dbf4a35bb0c 100644 --- a/src/kick.h +++ b/src/kick.h @@ -71,7 +71,7 @@ __attribute__((always_inline)) INLINE static void kick_gpart( * @param dt_kick_corr The kick time-step for the gizmo-mfv gravity correction. * @param cosmo The cosmological model. * @param hydro_props The constants used in the scheme. - * @param entropy_floor_props Properties of the entropy floor. + * @param floor_props Properties of the entropy floor. * @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). @@ -80,8 +80,8 @@ __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, double dt_kick_corr, const struct cosmology *cosmo, const struct hydro_props *hydro_props, - const struct entropy_floor_properties *entropy_floor_props, - integertime_t ti_start, integertime_t ti_end) { + const struct entropy_floor_properties *floor_props, integertime_t ti_start, + integertime_t ti_end) { #ifdef SWIFT_DEBUG_CHECKS if (p->ti_kick != ti_start) @@ -114,15 +114,8 @@ __attribute__((always_inline)) INLINE static void kick_part( /* Extra kick work */ hydro_kick_extra(p, xp, dt_kick_therm, dt_kick_grav, dt_kick_hydro, - dt_kick_corr, cosmo, hydro_props); + dt_kick_corr, cosmo, hydro_props, floor_props); if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt_kick_grav); - - /* Verify that the particle is not below the entropy floor */ - const float floor = entropy_floor(p, cosmo, entropy_floor_props); - if (hydro_get_physical_entropy(p, xp, cosmo) < floor) { - hydro_set_physical_entropy(p, xp, cosmo, floor); - hydro_set_physical_internal_energy_dt(p, cosmo, 0.f); - } } /** diff --git a/src/star_formation.h b/src/star_formation.h index 5f873d2da142ec6fda90e986b523f60f7ef0d4ef..3fb6a93249aea13b35e44f6c6265a9d5dc105967 100644 --- a/src/star_formation.h +++ b/src/star_formation.h @@ -27,6 +27,9 @@ /* Config parameters. */ #include "../config.h" +/* Local includes */ +#include "hydro.h" + /* Import the right star formation law definition */ #if defined(STAR_FORMATION_NONE) #include "./star_formation/none/star_formation.h"