/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2019 Matthieu Schaller (schaller@strw.leidenuniv.nl)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see .
*
******************************************************************************/
#ifndef SWIFT_PS2020_COOLING_SUBGRID_H
#define SWIFT_PS2020_COOLING_SUBGRID_H
/* Config parameters. */
#include
/* Local includes. */
#include "cooling.h"
#include "equation_of_state.h"
/**
* @brief Computes the gas pressure at thermal equilibrium (cooling = heating)
* for gas with at a given redshift, abundance ratio, metallicity, and density
*
* @return The logarithm base 10 of the thermal pressure at thermal equilibrium
* in cgs
*
* @param cooling The properties of the cooling scheme.
* @param abundance_ratio element abundance ratio of gas particle
* @param log_u_cgs The logarithm base 10 of the internal energy in cgs
* @param log10nH_local The logarithm base 10 of the hydrogen particle density
* in cgs
* @param rho_cgs Physical density of the gas in cgs
* @param redshift The redshift
* @param ired Index of redshift within the redshift dimension of the cooling
* tables
* @param imet Index of metallicity within the metallicity dimension of the
* cooling tables
* @param dred delta redshift between redshift bins, for interpolation
* @param dmet delta metallicity between metallicity, for interpolation
*/
double get_thermal_equilibrium_pressure(
const struct cooling_function_data *cooling,
const float abundance_ratio[colibre_cooling_N_elementtypes],
const double log_u_cgs, const float log10nH_local, const double rho_cgs,
const double redshift, const int ired, const int imet, const float dred,
const float dmet) {
/* The thermal equilibrium might have multiple solutions for a given density
* (WNM, CNM), so it is dangerous to do a bisection; instead the internal
* energy is decreased in steps of dlogu until heating > cooling (lambda > 0),
* then interpolate between u_previous and u_current*/
const double du = exp10(0.3f);
/* Get H number density */
double nH_cgs = exp10(log10nH_local);
/* Get index along the density index of the table */
int iden;
float dden;
get_index_1d(cooling->nH, colibre_cooling_N_density, log10nH_local, &iden,
&dden);
/* Prepare for iterations */
double u_current = exp10(log_u_cgs);
double u_prev = u_current;
double Lambda_current = colibre_cooling_rate(
log10f(u_current), redshift, nH_cgs, abundance_ratio, iden, dden, imet,
dmet, ired, dred, cooling, /*onlyicool=*/0,
/*onlyiheat=*/0, /*icool=*/0, /*iheat=*/0);
/* Are we in the case where we actually heat?
* This can happen but should not */
if (Lambda_current > 0.f) {
/* Increase u_current by a factor of 10 */
u_current *= 10.f;
Lambda_current = colibre_cooling_rate(
log10f(u_current), redshift, nH_cgs, abundance_ratio, iden, dden, imet,
dmet, ired, dred, cooling, /*onlyicool=*/0,
/*onlyiheat=*/0, /*icool=*/0, /*iheat=*/0);
/* If we are still heating, print a message */
if (Lambda_current > 0.f) message("Should not be here!");
}
/* The equilibirum energy we are looking for */
double u_eq_cgs = -1.;
/* Main iteration */
do {
Lambda_current = colibre_cooling_rate(
log10f(u_current), redshift, nH_cgs, abundance_ratio, iden, dden, imet,
dmet, ired, dred, cooling, 0, 0, 0, 0);
/* still cooling, need to decrease u */
if (Lambda_current < 0) {
u_current = u_current / du;
/* heating: interpolate to find eq */
} else {
u_prev = u_current * du;
double Lambda_prev = colibre_cooling_rate(
log10f(u_prev), redshift, nH_cgs, abundance_ratio, iden, dden, imet,
dmet, ired, dred, cooling, /*onlyicool=*/0,
/*onlyiheat=*/0, /*icool=*/0, /*iheat=*/0);
if (Lambda_prev * Lambda_current > 0.) {
message("Something wrong! %.4e, %.4e", Lambda_prev, Lambda_current);
}
/* Interpolate to get the equilibrium internal energy */
u_eq_cgs = u_prev - Lambda_prev / (Lambda_current - Lambda_prev) *
(u_current - u_prev);
/* And we are done here */
break;
}
} while (u_current >= cooling->umin_cgs);
/* Equilibrium internal energy is below the minimum internal energy
* set u_eq_cgs to minimum internal energy */
if (u_eq_cgs < 0.f) {
u_eq_cgs = cooling->umin_cgs;
}
/* Finish by computing the pressure from the density
* and equilibirum internal energy */
const float P_eq_cgs = gas_pressure_from_internal_energy(rho_cgs, u_eq_cgs);
return (double)log10f(P_eq_cgs);
}
/**
* @brief Compute the subgrid properties of the gas.
* For particles on the entropy floor, we use pressure equilibrium to
* infer the properties of the particle.
* Note that we return the density in physical coordinates.
*
* @return depends on the input for isub:
* isub = colibre_compute_subgrid_density: the physical subgrid density in
* internal units isub = colibre_compute_subgrid_temperature: the subgrid
* temperature temperature isub = colibre_compute_subgrid_HI_fraction: the
* subgrid HI fraction (n_HI / n_H) isub =
* colibre_compute_subgrid_HII_fraction: the subgrid HII fraction (n_HII /
* n_H) isub = colibre_compute_subgrid_H2_fraction: the subgrid HII fraction
* (n_H2 / n_H)
*
* @param cooling The properties of the cooling scheme.
* @param phys_const The physical constants.
* @param floor_props The properties of the entropy floor.
* @param cosmo The cosmological model.
* @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 P_phys Pressure of the gas in internal physical units.
* @param log10_T The logarithm base 10 of the temperature of the gas.
* @param log10_u_EOS_max_cgs The logarithm base 10 of the maximal energy in cgs
* to be considered on the EOS at the density of the gas.
* @param HII_region Is this patch of gas in an HII region?
* @param abundance_ratio element abundance ratio of gas particle
* @param log_u_cgs The logarithm base 10 of the internal energy in cgs
* @param isub which subgrid property to calculate
*/
double compute_subgrid_property(
const struct cooling_function_data *cooling,
const struct phys_const *phys_const,
const struct entropy_floor_properties *floor_props,
const struct cosmology *cosmo, const float rho_phys, const float logZZsol,
const float XH, const float P_phys, const float log10_T,
const float log10_u_EOS_max_cgs, const int HII_region,
const float abundance_ratio[colibre_cooling_N_elementtypes],
const double log_u_cgs, const enum cooling_subgrid_properties isub) {
if (HII_region)
error("HII regions are not implemented in the EAGLE-XL flavour");
const float weights[3] = {1.0, 1.0, 1.0};
/* Set the value for returning the values based on the hydro properties,
* e.g. for particles above the EOS, or if no EOS is present */
double standard_return;
/* Start by computing the subgrid property as if it were not subgrid */
if (isub == cooling_compute_subgrid_density) {
standard_return = rho_phys;
} else if (isub == cooling_compute_subgrid_temperature) {
standard_return = exp10(log10_T);
} else {
/* Get the Hydrogen number density */
const double n_H = XH * rho_phys / phys_const->const_proton_mass;
const double n_H_cgs = n_H * cooling->number_density_to_cgs;
/* Get index along the different table axis for this particle */
int ired, imet, iden, item;
float dred, dmet, dden, dtem;
get_index_1d(cooling->Redshifts, colibre_cooling_N_redshifts, cosmo->z,
&ired, &dred);
get_index_1d(cooling->Metallicity, colibre_cooling_N_metallicity, logZZsol,
&imet, &dmet);
get_index_1d(cooling->nH, colibre_cooling_N_density, log10f(n_H_cgs), &iden,
&dden);
get_index_1d(cooling->Temp, colibre_cooling_N_temperature, log10_T, &item,
&dtem);
if (isub == cooling_compute_subgrid_HI_fraction) {
const float nHI_over_nH =
interpolation4d_plus_summation(cooling->table.logHfracs_all, /* */
weights, neutral, neutral, /* */
ired, item, imet, iden, /* */
dred, dtem, dmet, dden, /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_temperature, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density, /* */
3);
standard_return = nHI_over_nH;
} else if (isub == cooling_compute_subgrid_HII_fraction) {
const float nHII_over_nH =
interpolation4d_plus_summation(cooling->table.logHfracs_all, /* */
weights, ionized, ionized, /* */
ired, item, imet, iden, /* */
dred, dtem, dmet, dden, /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_temperature, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density, /* */
3);
standard_return = nHII_over_nH;
} else if (isub == cooling_compute_subgrid_H2_fraction) {
const float mH2_over_mH =
interpolation4d_plus_summation(cooling->table.logHfracs_all, /* */
weights, molecular, molecular, /* */
ired, item, imet, iden, /* */
dred, dtem, dmet, dden, /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_temperature, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density, /* */
3);
/* molecular fraction are stored as mass fractions in the table (2nH2/nH)
* but we want particle fractions (nH2/nH) as output */
standard_return = 0.5 * mH2_over_mH;
} else {
/* This should never happen! */
standard_return = -1.;
error("Unknown subgrid property to calculate");
}
}
/* We can now look at whether we need to correct that first guess */
/* Return case 1:
* Particle is above the equation of state, nothing to be done */
if (log_u_cgs >= log10_u_EOS_max_cgs) {
return standard_return;
}
/* Return case 2:
* Particle is in an HII region
* get density by mapping particle onto T_HII and mu_HII at equal pressure */
/* NOTE: HII regions DO NOT exit in the EAGLE-XL flavour
* --> Nothing to do here, move straight to case 3 */
double Lambda;
/* Get the Hydrogen number density */
const double n_H = XH * rho_phys / phys_const->const_proton_mass;
const double n_H_cgs = n_H * cooling->number_density_to_cgs;
/* Get index along the different table axis for this particle */
int ired, imet, iden, item;
float dred, dmet, dden, dtem;
get_index_1d(cooling->Redshifts, colibre_cooling_N_redshifts, cosmo->z, &ired,
&dred);
get_index_1d(cooling->Metallicity, colibre_cooling_N_metallicity, logZZsol,
&imet, &dmet);
get_index_1d(cooling->nH, colibre_cooling_N_density, log10f(n_H_cgs), &iden,
&dden);
get_index_1d(cooling->Temp, colibre_cooling_N_temperature, log10_T, &item,
&dtem);
/* Get the cooling rate for the particle properties */
Lambda =
colibre_cooling_rate(log_u_cgs, cosmo->z, n_H_cgs, abundance_ratio, iden,
dden, imet, dmet, ired, dred, cooling, 0, 0, 0, 0);
/* Return case 3:
* Particle is on the EOS, but heating dominates cooling, so the
* particle is below or at its thermal equilibrium temperature;
* This might happen for a high EOS and low densities, where
* radiative cooling does not dominate the total (net) cooling rate;
* In this case we return the hydro properties */
if (Lambda >= 0.) {
return standard_return;
}
/* Particle is on the EOS and cooling dominates heating, so the
* subgrid property is calculated by projecting the hydro properties
* on the thermal equilibrium temperature where the pressure equals
* the hydro pressure */
if (Lambda < 0) {
float log10_Teq = interpolation_3d(cooling->table.logTeq, /* */
ired, imet, iden, /* */
dred, dmet, dden, /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density);
const float P_cgs = P_phys * cooling->pressure_to_cgs;
const float log10_P_cgs = log10f(P_cgs);
/* Maximal equilibrium pressure (n = n_max, and Teq = Teq(n_max)) from the
* table at the current redshift and metallicity */
const float log10_Peq_max_cgs =
interpolation_3d_no_z(cooling->table.logPeq, /* */
ired, imet, colibre_cooling_N_density - 1, /* */
dred, dmet, 0., /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density);
/* Return case 4:
* EOS pressure (logP) is larger than maximum Peq (possible for very
* steep EOS); use the equilibrium temperature from the highest density
* bin */
if (log10_P_cgs > log10_Peq_max_cgs) {
if (isub == cooling_compute_subgrid_density) {
const double rho_cgs =
exp10(cooling->nH[colibre_cooling_N_density - 1]) *
cooling->proton_mass_cgs / XH;
return rho_cgs * cooling->density_from_cgs;
} else if (isub == cooling_compute_subgrid_temperature) {
const float log10_T_at_Peq = interpolation_3d_no_z(
cooling->table.logTeq, ired, imet, colibre_cooling_N_density - 1,
dred, dmet, 0., colibre_cooling_N_redshifts,
colibre_cooling_N_metallicity, colibre_cooling_N_density);
return exp10(log10_T_at_Peq);
} else if (isub == cooling_compute_subgrid_HI_fraction) {
const float log10_HI_over_nH = interpolation_4d_no_z_no_w(
cooling->table.logHfracs_Teq, /* */
ired, imet, colibre_cooling_N_density - 1, neutral, /* */
dred, dmet, 0., 0., /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density, /* */
3);
return exp10(log10_HI_over_nH);
} else if (isub == cooling_compute_subgrid_HII_fraction) {
const float log10_HII_over_nH = interpolation_4d_no_z_no_w(
cooling->table.logHfracs_Teq, /* */
ired, imet, colibre_cooling_N_density - 1, ionized, /* */
dred, dmet, 0., 0., /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density, /* */
3);
return exp10(log10_HII_over_nH);
} else if (isub == cooling_compute_subgrid_H2_fraction) {
const float log10_H2_over_nH = interpolation_4d_no_z_no_w(
cooling->table.logHfracs_Teq, /* */
ired, imet, colibre_cooling_N_density - 1, molecular, /* */
dred, dmet, 0., 0., /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density, /* */
3);
return 0.5 * exp10(log10_H2_over_nH);
} else {
error("Unknown subgrid property to calculate");
return 0.;
}
}
/* Find equilibrium property by increasing the density
* Note that the logPeq table is neither equally spaced nor
* necessarilly monotically increasing.
* We hence loop over densities and pick the first one where
* log10_P < log10_Peq. We start with the next resolved density (index
* iden+1), as we require the subgrid density to be larger than the
* physical density */
float log10_Peq_prev = log10_P_cgs;
for (int i = iden + 1; i < colibre_cooling_N_density; i++) {
float log10_Peq_current, log10_n_at_Peq;
if ((log10_Teq <= log10_T) &&
(logZZsol <=
cooling->Metallicity[colibre_cooling_N_metallicity - 1])) {
/* Standard case: the thermal equilibrium temperature is below the
* EOS temperature and the metallicity is within the table range */
log10_Peq_current = interpolation_3d_no_z(
cooling->table.logPeq, ired, imet, i, dred, dmet, 0.,
colibre_cooling_N_redshifts, colibre_cooling_N_metallicity,
colibre_cooling_N_density);
} else {
/* The thermal equilibrium tables cannot be used, need to calculate
* the thermal equilibrium temperatures; this only occurs when
* the metallicity is higher than the maximum metallicity in the
* tables */
const double rho_cgs =
exp10(cooling->nH[i]) * cooling->proton_mass_cgs / XH;
log10_Peq_current = get_thermal_equilibrium_pressure(
cooling, abundance_ratio, log_u_cgs, cooling->nH[i], rho_cgs,
cosmo->z, ired, imet, dred, dmet);
}
if (log10_Peq_current >= log10_P_cgs) {
if (log10_Peq_current == log10_Peq_prev) {
/* the interpolated pressure exactly equals the input pressure
* on the first density iteration (where log10_Peq_prev =
* log10_P_cgs), no interpolation necessary */
log10_n_at_Peq = cooling->nH[i];
} else {
/* first density bin where log P_eq >= log P_SPH
* the solution is therefore between log10_Peq_prev and
* log10_Peq_current */
const float delta_P_eq = (log10_P_cgs - log10_Peq_prev) /
(log10_Peq_current - log10_Peq_prev);
/* Interpolate to get the density at equilibrium */
log10_n_at_Peq = cooling->nH[i - 1] +
delta_P_eq * (cooling->nH[i] - cooling->nH[i - 1]);
}
/* Return case 5:
* SPH density and temperatures are projected onto the thermal
* equilibrium temperature function Teq(density,
* metallicity/abundances, redshift) for equal pressure */
if (isub == cooling_compute_subgrid_density) {
const double rho_cgs =
exp10(log10_n_at_Peq) * cooling->proton_mass_cgs / XH;
return rho_cgs * cooling->density_from_cgs;
} else {
/* We found a valid density: Get the index along the density
* axis of the tables. */
int iden_eq;
float dden_eq;
get_index_1d(cooling->nH, colibre_cooling_N_density, log10_n_at_Peq,
&iden_eq, &dden_eq);
/* Finish by interpolating the tables to get the temperature */
const float log10_T_at_Peq =
interpolation_3d(cooling->table.logTeq, /* */
ired, imet, iden_eq, /* */
dred, dmet, dden_eq, /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density);
if (isub == cooling_compute_subgrid_temperature) {
return exp10(log10_T_at_Peq);
} else {
int item_eq;
float dtem_eq;
get_index_1d(cooling->Temp, colibre_cooling_N_temperature,
log10_T_at_Peq, &item_eq, &dtem_eq);
if (isub == cooling_compute_subgrid_HI_fraction) {
const float nHI_over_nH_eq = interpolation4d_plus_summation(
cooling->table.logHfracs_all, /* */
weights, neutral, neutral, /* */
ired, item_eq, imet, iden_eq, /* */
dred, dtem_eq, dmet, dden_eq, /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_temperature, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density, /* */
3);
return nHI_over_nH_eq;
} else if (isub == cooling_compute_subgrid_HII_fraction) {
const float nHII_over_nH_eq = interpolation4d_plus_summation(
cooling->table.logHfracs_all, /* */
weights, ionized, ionized, /* */
ired, item_eq, imet, iden_eq, /* */
dred, dtem_eq, dmet, dden_eq, /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_temperature, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density, /* */
3);
return nHII_over_nH_eq;
} else if (isub == cooling_compute_subgrid_H2_fraction) {
const float mH2_over_mH_eq = interpolation4d_plus_summation(
cooling->table.logHfracs_all, /* */
weights, molecular, molecular, /* */
ired, item_eq, imet, iden_eq, /* */
dred, dtem_eq, dmet, dden_eq, /* */
colibre_cooling_N_redshifts, /* */
colibre_cooling_N_temperature, /* */
colibre_cooling_N_metallicity, /* */
colibre_cooling_N_density, /* */
3);
return 0.5 * mH2_over_mH_eq;
} else {
error("Unknown subgrid property to calculate");
return 0.;
}
}
}
}
/* Move to the next iteration */
log10_Peq_prev = log10_Peq_current;
}
}
/* Return case 6:
* Nothing worked, return the property at the SPH density and temperature
* The code should not reach this point during the simulation, but still
* does before and at the first timestep */
return standard_return;
}
#endif /* SWIFT_PS2020_COOLING_SUBGRID_H */