Skip to content

SWIFTsim
SWIFTsim
Commit a26725c1 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files
Options

Merge branch 'QLA_P20' into 'master' 

Use the P20 cooling tables for the Quick Lyman-alpha model

See merge request !1138
parents 9cdbe7b0 e7a5317e
Hide whitespace changes
Inline Side-by-side
Showing with 1867 additions and 1375 deletions
configure.ac
View file @ a26725c1
......@@ -1489,7 +1489,7 @@ AC_SUBST([NUMA_INCS])
# As an example for this, see the call to AC_ARG_WITH for cooling.
AC_ARG_WITH([subgrid],
[AS_HELP_STRING([--with-subgrid=<subgrid>],
[Master switch for subgrid methods. Inexperienced user should start here. Options are: @<:@none, GEAR, QLA, EAGLE default: none@:>@]
[Master switch for subgrid methods. Inexperienced user should start here. Options are: @<:@none, GEAR, QLA, EAGLE, EAGLE-XL default: none@:>@]
)],
[with_subgrid="$withval"],
[with_subgrid=none]
......
examples/QuickLymanAlpha/L050N0752/qla_50.yml
View file @ a26725c1
......@@ -86,12 +86,13 @@ LineOfSight:
# Quick Lyman-alpha cooling (EAGLE with fixed primoridal Z)
QLACooling:
dir_name: ./coolingtables/
H_reion_z: 7.5 # Planck 2018
dir_name: ./UV_dust1_CR1_G1_shield1.hdf5 # Location of the cooling tables
H_reion_z: 7.5 # Redshift of Hydrogen re-ionization (Planck 2018)
H_reion_eV_p_H: 2.0
He_reion_z_centre: 3.5
He_reion_z_sigma: 0.5
He_reion_eV_p_H: 2.0
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
rapid_cooling_threshold: 0.333333 # Switch to rapid cooling regime for dt / t_cool above this threshold.
# Quick Lyman-alpha star formation parameters
QLAStarFormation:
......
examples/QuickLymanAlpha/getQLACoolingTable.sh
View file @ a26725c1
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/EAGLE/coolingtables.tar.gz
tar -xf coolingtables.tar.gz
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/COLIBRE/UV_dust1_CR1_G1_shield1.hdf5
examples/parameter_example.yml
View file @ a26725c1
......@@ -389,14 +389,26 @@ EAGLECooling:
Ca_over_Si_in_solar: 1. # (Optional) Ratio of Ca/Si to use in units of solar. If set to 1, the code uses [Ca/Si] = 0, i.e. Ca/Si = 0.0941736.
S_over_Si_in_solar: 1. # (Optional) Ratio of S/Si to use in units of solar. If set to 1, the code uses [S/Si] = 0, i.e. S/Si = 0.6054160.
# Quick Lyman-alpha cooling (EAGLE with fixed primoridal Z)
# Quick Lyman-alpha cooling (EAGLE-XL with fixed primoridal Z)
QLACooling:
dir_name: ./coolingtables/ # Location of the Wiersma+08 cooling tables
H_reion_z: 7.5 # Redshift of Hydrogen re-ionization
H_reion_eV_p_H: 2.0 # Energy inject by Hydrogen re-ionization in electron-volt per Hydrogen atom
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
dir_name: ./UV_dust1_CR1_G1_shield1.hdf5 # Location of the Ploeckinger+20 cooling tables
H_reion_z: 7.5 # Redshift of Hydrogen re-ionization (Planck 2018)
H_reion_eV_p_H: 2.0 # Energy inject by Hydrogen re-ionization in electron-volt per Hydrogen atom
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
rapid_cooling_threshold: 0.333333 # Switch to rapid cooling regime for dt / t_cool above this threshold.
# COLIBRE cooling parameters (EAGLE-XL)
COLIBRECooling:
dir_name: ./UV_dust1_CR1_G1_shield1.hdf5 # Location of the cooling tables
H_reion_z: 7.5 # Redshift of Hydrogen re-ionization (Planck 2018)
H_reion_eV_p_H: 2.0 # Energy inject by Hydrogen re-ionization in electron-volt per Hydrogen atom
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
rapid_cooling_threshold: 0.333333 # Switch to rapid cooling regime for dt / t_cool above this threshold.
delta_logTEOS_subgrid_properties: 0.3 # delta log T above the EOS below which the subgrid properties use Teq assumption
# Cooling with Grackle 3.0
GrackleCooling:
......
src/cooling/QLA/cooling.c
View file @ a26725c1
......@@ -18,7 +18,7 @@
******************************************************************************/
/**
* @file src/cooling/QLA/cooling.c
* @brief Quick Lyman-alpha cooling functions
* @brief QLA cooling functions
*/
/* Config parameters. */
......@@ -31,6 +31,7 @@
#include <time.h>
/* Local includes. */
#include "adiabatic_index.h"
#include "chemistry.h"
#include "cooling.h"
#include "cooling_rates.h"
......@@ -48,7 +49,8 @@
#include "space.h"
#include "units.h"
/* Maximum number of iterations for bisection scheme */
/* Maximum number of iterations for
* bisection integration schemes */
static const int bisection_max_iterations = 150;
/* Tolerances for termination criteria. */
......@@ -56,64 +58,6 @@ static const float explicit_tolerance = 0.05;
static const float bisection_tolerance = 1.0e-6;
static const double bracket_factor = 1.5;
/**
* @brief Find the index of the current redshift along the redshift dimension
* of the cooling tables.
*
* Since the redshift table is not evenly spaced, compare z with each
* table value in decreasing order starting with the previous redshift index
*
* The returned difference is expressed in units of the table separation. This
* means dx = (x - table[i]) / (table[i+1] - table[i]). It is always between
* 0 and 1.
*
* @param z Redshift we are searching for.
* @param z_index (return) Index of the redshift in the table.
* @param dz (return) Difference in redshift between z and table[z_index].
* @param cooling #cooling_function_data structure containing redshift table.
*/
__attribute__((always_inline)) INLINE void get_redshift_index(
const float z, int *z_index, float *dz,
struct cooling_function_data *restrict cooling) {
/* Before the earliest redshift or before hydrogen reionization, flag for
* collisional cooling */
if (z > cooling->H_reion_z) {
*z_index = qla_cooling_N_redshifts;
*dz = 0.0;
}
/* From reionization use the cooling tables */
else if (z > cooling->Redshifts[qla_cooling_N_redshifts - 1] &&
z <= cooling->H_reion_z) {
*z_index = qla_cooling_N_redshifts + 1;
*dz = 0.0;
}
/* At the end, just use the last value */
else if (z <= cooling->Redshifts[0]) {
*z_index = 0;
*dz = 0.0;
}
/* Normal case: search... */
else {
/* start at the previous index and search */
for (int i = cooling->previous_z_index; i >= 0; i--) {
if (z > cooling->Redshifts[i]) {
*z_index = i;
cooling->previous_z_index = i;
*dz = (z - cooling->Redshifts[i]) /
(cooling->Redshifts[i + 1] - cooling->Redshifts[i]);
break;
}
}
}
}
/**
* @brief Common operations performed on the cooling function at a
* given time-step or redshift. Predominantly used to read cooling tables
......@@ -128,15 +72,6 @@ __attribute__((always_inline)) INLINE void get_redshift_index(
void cooling_update(const struct cosmology *cosmo,
struct cooling_function_data *cooling, struct space *s) {
/* Current redshift */
const float redshift = cosmo->z;
/* What is the current table index along the redshift axis? */
int z_index = -1;
float dz = 0.f;
get_redshift_index(redshift, &z_index, &dz, cooling);
cooling->dz = dz;
/* Extra energy for reionization? */
if (!cooling->H_reion_done) {
......@@ -153,35 +88,164 @@ void cooling_update(const struct cosmology *cosmo,
cooling->H_reion_done = 1;
}
}
}
/* Do we already have the correct tables loaded? */
if (cooling->z_index == z_index) return;
/**
* @brief Compute the internal energy of a #part based on the cooling function
* but for a given temperature.
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param p #part data.
* @param xp Pointer to the #xpart data.
* @param T temperature of the gas (internal units).
*/
float cooling_get_internalenergy_for_temperature(
const struct phys_const *phys_const, const struct hydro_props *hydro_props,
const struct unit_system *us, const struct cosmology *cosmo,
const struct cooling_function_data *cooling, const struct part *p,
const struct xpart *xp, const float T) {
/* Which table should we load ? */
if (z_index >= qla_cooling_N_redshifts) {
#ifdef SWIFT_DEBUG_CHECKS
if (cooling->Redshifts == NULL)
error(
"Cooling function has not been initialised. Did you forget the "
"--temperature runtime flag?");
#endif
if (z_index == qla_cooling_N_redshifts + 1) {
/* Get the Hydrogen mass fraction */
const float XH = 1. - phys_const->const_primordial_He_fraction;
/* Bewtween re-ionization and first table */
get_redshift_invariant_table(cooling, /* photodis=*/0);
/* Convert Hydrogen mass fraction into Hydrogen number density */
const float rho = hydro_get_physical_density(p, cosmo);
const double n_H = rho * XH / phys_const->const_proton_mass;
const double n_H_cgs = n_H * cooling->number_density_to_cgs;
} else {
/* Get this particle's metallicity ratio to solar.
*
* Note that we do not need the individual element's ratios that
* the function also computes. */
float dummy[qla_cooling_N_elementtypes];
const float logZZsol =
abundance_ratio_to_solar(p, cooling, phys_const, dummy);
/* Above re-ionization */
get_redshift_invariant_table(cooling, /* photodis=*/1);
}
/* compute hydrogen number density, metallicity and redshift indices and
* offsets */
} else {
float d_red, d_met, d_n_H;
int red_index, met_index, n_H_index;
/* Normal case: two tables bracketing the current z */
const int low_z_index = z_index;
const int high_z_index = z_index + 1;
get_index_1d(cooling->Redshifts, qla_cooling_N_redshifts, cosmo->z,
&red_index, &d_red);
get_index_1d(cooling->Metallicity, qla_cooling_N_metallicity, logZZsol,
&met_index, &d_met);
get_index_1d(cooling->nH, qla_cooling_N_density, log10(n_H_cgs), &n_H_index,
&d_n_H);
get_cooling_table(cooling, low_z_index, high_z_index);
}
/* Compute the log10 of the temperature by interpolating the table */
const double log_10_U =
qla_convert_temp_to_u(log10(T), cosmo->z, n_H_index, d_n_H, met_index,
d_met, red_index, d_red, cooling);
/* Undo the log! */
return exp10(log_10_U);
}
/**
* @brief Compute the temperature based on gas properties.
*
* The temperature returned is consistent with the cooling rates.
*
* @param phys_const #phys_const data structure.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param rho_phys Density of the gas in internal physical units.
* @param logZZsol Logarithm base 10 of the gas' metallicity in units of solar
* metallicity.
* @param XH The Hydrogen abundance of the gas.
* @param u_phys Internal energy of the gas in internal physical units.
*/
float cooling_get_temperature_from_gas(
const struct phys_const *phys_const, const struct cosmology *cosmo,
const struct cooling_function_data *cooling, const float rho_phys,
const float logZZsol, const float XH, const float u_phys) {
/* Convert to CGS */
const double u_cgs = u_phys * cooling->internal_energy_to_cgs;
const double n_H = rho_phys * XH / phys_const->const_proton_mass;
const double n_H_cgs = n_H * cooling->number_density_to_cgs;
/* compute hydrogen number density, metallicity and redshift indices and
* offsets */
float d_red, d_met, d_n_H;
int red_index, met_index, n_H_index;
get_index_1d(cooling->Redshifts, qla_cooling_N_redshifts, cosmo->z,
&red_index, &d_red);
get_index_1d(cooling->Metallicity, qla_cooling_N_metallicity, logZZsol,
&met_index, &d_met);
get_index_1d(cooling->nH, qla_cooling_N_density, log10(n_H_cgs), &n_H_index,
&d_n_H);
/* Compute the log10 of the temperature by interpolating the table */
const double log_10_T =
qla_convert_u_to_temp(log10(u_cgs), cosmo->z, n_H_index, d_n_H, met_index,
d_met, red_index, d_red, cooling);
/* Undo the log! */
return exp10(log_10_T);
}
/**
* @brief Compute the temperature of a #part based on the cooling function.
*
* The temperature returned is consistent with the cooling rates.
*
* @param phys_const #phys_const data structure.
* @param hydro_props The properties of the hydro scheme.
* @param us The internal system of units.
* @param cosmo #cosmology data structure.
* @param cooling #cooling_function_data struct.
* @param p #part data.
* @param xp Pointer to the #xpart data.
*/
float cooling_get_temperature(const struct phys_const *phys_const,
const struct hydro_props *hydro_props,
const struct unit_system *us,
const struct cosmology *cosmo,
const struct cooling_function_data *cooling,
const struct part *p, const struct xpart *xp) {
#ifdef SWIFT_DEBUG_CHECKS
if (cooling->Redshifts == NULL)
error(
"Cooling function has not been initialised. Did you forget the "
"--temperature runtime flag?");
#endif
/* Get physical internal energy */
const float u_phys = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Get the Hydrogen mass fraction */
const float XH = 1. - phys_const->const_primordial_He_fraction;
/* Store the currently loaded index */
cooling->z_index = z_index;
/* Convert Hydrogen mass fraction into Hydrogen number density */
const float rho_phys = hydro_get_physical_density(p, cosmo);
/* Get this particle's metallicity ratio to solar.
*
* Note that we do not need the individual element's ratios that
* the function also computes. */
float dummy[qla_cooling_N_elementtypes];
const float logZZsol =
abundance_ratio_to_solar(p, cooling, phys_const, dummy);
return cooling_get_temperature_from_gas(phys_const, cosmo, cooling, rho_phys,
logZZsol, XH, u_phys);
}
/**
......@@ -192,8 +256,10 @@ void cooling_update(const struct cosmology *cosmo,
* @param redshift Current redshift.
* @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 met_index Particle metallicity index.
* @param d_met Particle metallicity offset.
* @param red_index Redshift index.
* @param d_red Redshift offset.
* @param Lambda_He_reion_cgs Cooling rate coming from He reionization.
* @param ratefact_cgs Multiplication factor to get a cooling rate.
* @param cooling #cooling_function_data structure.
......@@ -201,18 +267,17 @@ void cooling_update(const struct cosmology *cosmo,
* @param dt_cgs timestep in CGS.
* @param ID ID of the particle (for debugging).
*/
INLINE static double bisection_iter(
static INLINE double bisection_iter(
const double u_ini_cgs, const double n_H_cgs, const double redshift,
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[qla_cooling_N_abundances], const double dt_cgs,
const long long ID) {
int n_H_index, float d_n_H, int met_index, float d_met, int red_index,
float d_red, double Lambda_He_reion_cgs, double ratefact_cgs,
const struct cooling_function_data *cooling,
const float abundance_ratio[qla_cooling_N_elementtypes], double dt_cgs,
long long ID) {
/* Bracketing */
double u_lower_cgs = u_ini_cgs;
double u_upper_cgs = u_ini_cgs;
double u_lower_cgs = max(u_ini_cgs, cooling->umin_cgs);
double u_upper_cgs = max(u_ini_cgs, cooling->umin_cgs);
/*************************************/
/* Let's get a first guess */
......@@ -221,7 +286,8 @@ INLINE static double bisection_iter(
double LambdaNet_cgs =
Lambda_He_reion_cgs +
qla_cooling_rate(log10(u_ini_cgs), redshift, n_H_cgs, abundance_ratio,
n_H_index, d_n_H, He_index, d_He, cooling);
n_H_index, d_n_H, met_index, d_met, red_index, d_red,
cooling, 0, 0, 0, 0);
/*************************************/
/* Let's try to bracket the solution */
......@@ -230,36 +296,52 @@ INLINE static double bisection_iter(
if (LambdaNet_cgs < 0) {
/* we're cooling! */
u_lower_cgs /= bracket_factor;
u_upper_cgs *= bracket_factor;
u_lower_cgs = max(u_lower_cgs / bracket_factor, cooling->umin_cgs);
u_upper_cgs = max(u_upper_cgs * bracket_factor, cooling->umin_cgs);
/* Compute a new rate */
LambdaNet_cgs =
Lambda_He_reion_cgs +
qla_cooling_rate(log10(u_lower_cgs), redshift, n_H_cgs, abundance_ratio,
n_H_index, d_n_H, He_index, d_He, cooling);
n_H_index, d_n_H, met_index, d_met, red_index, d_red,
cooling, 0, 0, 0, 0);
int i = 0;
while (u_lower_cgs - u_ini_cgs - LambdaNet_cgs * ratefact_cgs * dt_cgs >
0 &&
i < bisection_max_iterations) {
u_lower_cgs /= bracket_factor;
u_upper_cgs /= bracket_factor;
u_lower_cgs = max(u_lower_cgs / bracket_factor, cooling->umin_cgs);
u_upper_cgs = max(u_upper_cgs / bracket_factor, cooling->umin_cgs);
/* Compute a new rate */
LambdaNet_cgs = Lambda_He_reion_cgs +
qla_cooling_rate(log10(u_lower_cgs), redshift, n_H_cgs,
abundance_ratio, n_H_index, d_n_H,
He_index, d_He, cooling);
LambdaNet_cgs =
Lambda_He_reion_cgs +
qla_cooling_rate(log10(u_lower_cgs), redshift, n_H_cgs,
abundance_ratio, n_H_index, d_n_H, met_index, d_met,
red_index, d_red, cooling, 0, 0, 0, 0);
/* If the energy is below or equal the minimum energy and we are still
* cooling, return the minimum energy */
if ((u_lower_cgs <= cooling->umin_cgs) && (LambdaNet_cgs < 0.))
return cooling->umin_cgs;
i++;
}
if (i >= bisection_max_iterations) {
error(
"particle %llu exceeded max iterations searching for bounds when "
"cooling, u_ini_cgs %.5e n_H_cgs %.5e",
ID, u_ini_cgs, n_H_cgs);
"cooling \n more info: n_H_cgs = %.4e, u_ini_cgs = %.4e, redshift = "
"%.4f\n"
"n_H_index = %i, d_n_H = %.4f\n"
"met_index = %i, d_met = %.4f, red_index = %i, d_red = %.4f, initial "
"Lambda = %.4e",
ID, n_H_cgs, u_ini_cgs, redshift, n_H_index, d_n_H, met_index, d_met,
red_index, d_red,
qla_cooling_rate(log10(u_ini_cgs), redshift, n_H_cgs, abundance_ratio,
n_H_index, d_n_H, met_index, d_met, red_index, d_red,
cooling, 0, 0, 0, 0));
}
} else {
......@@ -271,7 +353,8 @@ INLINE static double bisection_iter(
LambdaNet_cgs =
Lambda_He_reion_cgs +
qla_cooling_rate(log10(u_upper_cgs), redshift, n_H_cgs, abundance_ratio,
n_H_index, d_n_H, He_index, d_He, cooling);
n_H_index, d_n_H, met_index, d_met, red_index, d_red,
cooling, 0, 0, 0, 0);
int i = 0;
while (u_upper_cgs - u_ini_cgs - LambdaNet_cgs * ratefact_cgs * dt_cgs <
......@@ -282,18 +365,33 @@ INLINE static double bisection_iter(
u_upper_cgs *= bracket_factor;
/* Compute a new rate */
LambdaNet_cgs = Lambda_He_reion_cgs +
qla_cooling_rate(log10(u_upper_cgs), redshift, n_H_cgs,
abundance_ratio, n_H_index, d_n_H,
He_index, d_He, cooling);
LambdaNet_cgs =
Lambda_He_reion_cgs +
qla_cooling_rate(log10(u_upper_cgs), redshift, n_H_cgs,
abundance_ratio, n_H_index, d_n_H, met_index, d_met,
red_index, d_red, cooling, 0, 0, 0, 0);
i++;
}
if (i >= bisection_max_iterations) {
message("Aborting...");
message("particle %llu", ID);
message("n_H_cgs = %.4e", n_H_cgs);
message("u_ini_cgs = %.4e", u_ini_cgs);
message("redshift = %.4f", redshift);
message("indices nH, met, red = %i, %i, %i", n_H_index, met_index,
red_index);
message("index weights nH, met, red = %.4e, %.4e, %.4e", d_n_H, d_met,
d_red);
fflush(stdout);
message("cooling rate = %.4e",
qla_cooling_rate(log10(u_ini_cgs), redshift, n_H_cgs,
abundance_ratio, n_H_index, d_n_H, met_index,
d_met, red_index, d_red, cooling, 0, 0, 0, 0));
error(
"particle %llu exceeded max iterations searching for bounds when "
"heating, u_ini_cgs %.5e n_H_cgs %.5e",
ID, u_ini_cgs, n_H_cgs);
"cooling",
ID);
}
}
......@@ -315,14 +413,8 @@ INLINE static double bisection_iter(
LambdaNet_cgs =
Lambda_He_reion_cgs +
qla_cooling_rate(log10(u_next_cgs), redshift, n_H_cgs, abundance_ratio,
n_H_index, d_n_H, He_index, d_He, cooling);
#ifdef SWIFT_DEBUG_CHECKS
if (u_next_cgs <= 0)
error(
"Got negative energy! u_next_cgs=%.5e u_upper=%.5e u_lower=%.5e "
"Lambda=%.5e",
u_next_cgs, u_upper_cgs, u_lower_cgs, LambdaNet_cgs);
#endif
n_H_index, d_n_H, met_index, d_met, red_index, d_red,
cooling, 0, 0, 0, 0);
/* Where do we go next? */
if (u_next_cgs - u_ini_cgs - LambdaNet_cgs * ratefact_cgs * dt_cgs > 0.0) {
......@@ -356,9 +448,7 @@ INLINE static double bisection_iter(
*
* f(u_new) = u_new - u_old - dt * du/dt(u_new, X) = 0
*
* We first try a few Newton-Raphson iteration if it does not converge, we
* revert to a bisection scheme.
*
* A bisection scheme is used.
* This is done by first bracketing the solution and then iterating
* towards the solution by reducing the window down to a certain tolerance.
* Note there is always at least one solution since
......@@ -374,8 +464,7 @@ INLINE static double bisection_iter(
* @param xp Pointer to the extended particle data.
* @param dt The cooling time-step of this particle.
* @param dt_therm The hydro time-step of this particle.
* @param time The current time (since the Big Bang or start of the run) in
* internal units.
* @param time Time since Big Bang
*/
void cooling_cool_part(const struct phys_const *phys_const,
const struct unit_system *us,
......@@ -383,12 +472,13 @@ void cooling_cool_part(const struct phys_const *phys_const,
const struct hydro_props *hydro_properties,
const struct entropy_floor_properties *floor_props,
const struct cooling_function_data *cooling,
struct part *restrict p, struct xpart *restrict xp,
const float dt, const float dt_therm,
const double time) {
struct part *p, struct xpart *xp, const float dt,
const float dt_therm, const double time) {
/* No cooling happens over zero time */
if (dt == 0.) return;
if (dt == 0.) {
return;
}
#ifdef SWIFT_DEBUG_CHECKS
if (cooling->Redshifts == NULL)
......@@ -403,7 +493,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 step (assuming dt does not
/* Get internal energy at the end of the next kick step (assuming dt does not
* increase) */
double u_0 = (u_start + hydro_du_dt * dt_therm);
......@@ -421,19 +511,15 @@ void cooling_cool_part(const struct phys_const *phys_const,
/* Get this particle's abundance ratios compared to solar
* Note that we need to add S and Ca that are in the tables but not tracked
* by the particles themselves.
* The order is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe] */
float abundance_ratio[qla_cooling_N_abundances];
abundance_ratio_to_solar(p, cooling, phys_const, abundance_ratio);
* The order is [H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe, OA] */
float abundance_ratio[qla_cooling_N_elementtypes];
float logZZsol =
abundance_ratio_to_solar(p, cooling, phys_const, abundance_ratio);
/* Get the Hydrogen and Helium mass fractions */
const float XH = 1. - phys_const->const_primordial_He_fraction;
const float XHe = phys_const->const_primordial_He_fraction;
/* Get the Helium mass fraction. Note that this is He / (H + He), i.e. a
* metal-free Helium mass fraction as per the Wiersma+08 definition */
const float HeFrac = XHe / (XH + XHe);
/* convert Hydrogen mass fraction into physical Hydrogen number density */
/* convert Hydrogen mass fraction into 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;
......@@ -442,13 +528,17 @@ void cooling_cool_part(const struct phys_const *phys_const,
* equivalent expression below */
const double ratefact_cgs = n_H_cgs * (XH * cooling->inv_proton_mass_cgs);
/* compute hydrogen number density and helium fraction table indices and
/* compute hydrogen number density, metallicity and redshift indices and
* offsets (These are fixed for any value of u, so no need to recompute them)
*/
int He_index, n_H_index;
float d_He, d_n_H;
get_index_1d(cooling->HeFrac, qla_cooling_N_He_frac, HeFrac, &He_index,
&d_He);
float d_red, d_met, d_n_H;
int red_index, met_index, n_H_index;
get_index_1d(cooling->Redshifts, qla_cooling_N_redshifts, cosmo->z,
&red_index, &d_red);
get_index_1d(cooling->Metallicity, qla_cooling_N_metallicity, logZZsol,
&met_index, &d_met);
get_index_1d(cooling->nH, qla_cooling_N_density, log10(n_H_cgs), &n_H_index,
&d_n_H);
......@@ -457,22 +547,21 @@ void cooling_cool_part(const struct phys_const *phys_const,
/* Get helium and hydrogen reheating term */
const double Helium_reion_heat_cgs =
qla_helium_reionization_extraheat(cosmo->z, delta_redshift, cooling);
eagle_helium_reionization_extraheat(cosmo->z, delta_redshift, cooling);
/* Convert this into a rate */
const double Lambda_He_reion_cgs =
Helium_reion_heat_cgs / (dt_cgs * ratefact_cgs);
/* Let's compute the internal energy at the end of the step */
/* Initialise to the initial energy to appease compiler; this will never not
be overwritten. */
double u_final_cgs = u_0_cgs;
double u_final_cgs;
/* First try an explicit integration (note we ignore the derivative) */
const double LambdaNet_cgs =
Lambda_He_reion_cgs + qla_cooling_rate(log10(u_0_cgs), cosmo->z, n_H_cgs,
abundance_ratio, n_H_index, d_n_H,
He_index, d_He, cooling);
met_index, d_met, red_index, d_red,
cooling, 0, 0, 0, 0);
/* if cooling rate is small, take the explicit solution */
if (fabs(ratefact_cgs * LambdaNet_cgs * dt_cgs) <
......@@ -482,11 +571,10 @@ void cooling_cool_part(const struct phys_const *phys_const,
} else {
/* 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);
bisection_iter(u_0_cgs, n_H_cgs, cosmo->z, n_H_index, d_n_H, met_index,
d_met, red_index, d_red, Lambda_He_reion_cgs,
ratefact_cgs, cooling, abundance_ratio, dt_cgs, p->id);
}
/* Convert back to internal units */
......@@ -509,11 +597,34 @@ void cooling_cool_part(const struct phys_const *phys_const,
(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;
/* Determine if we are in the slow- or rapid-cooling regime,
* by comparing dt / t_cool to the rapid_cooling_threshold.
*
* Note that dt / t_cool = fabs(delta_u) / u_start. */
const double dt_over_t_cool = fabs(delta_u) / max(u_start, u_floor);
/* If rapid_cooling_threshold < 0, always use the slow-cooling
* regime. */
if ((cooling->rapid_cooling_threshold >= 0.0) &&
(dt_over_t_cool >= cooling->rapid_cooling_threshold)) {
/* Update the internal energy time derivative */
hydro_set_physical_internal_energy_dt(p, cosmo, cooling_du_dt);
/* Rapid-cooling regime. */
/* Update the particle's u and du/dt */
hydro_set_physical_internal_energy(p, xp, cosmo, u_final);
hydro_set_drifted_physical_internal_energy(p, cosmo, u_final);
hydro_set_physical_internal_energy_dt(p, cosmo, 0.);
} else {
/* Slow-cooling regime. */
/* Update du/dt so that we can subsequently drift internal energy. */
const float cooling_du_dt = delta_u / dt_therm;
/* Update the internal energy time derivative */
hydro_set_physical_internal_energy_dt(p, cosmo, cooling_du_dt);
}
}
/**
......@@ -530,12 +641,10 @@ void cooling_cool_part(const struct phys_const *phys_const,
* @param xp extended particle data.
*/
__attribute__((always_inline)) INLINE float cooling_timestep(
const struct cooling_function_data *restrict cooling,
const struct phys_const *restrict phys_const,
const struct cosmology *restrict cosmo,
const struct unit_system *restrict us,
const struct