/******************************************************************************* * 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_COLIBRE_COOLING_SUBGRID_H #define SWIFT_COLIBRE_COOLING_SUBGRID_H /* Config parameters. */ #include /* Local includes. */ #include "colibre_tables_restrict.h" #include "cooling.h" #include "equation_of_state.h" /** * @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 = subgrid_density: the physical subgrid density in internal units * isub = subgrid_temperature: the subgrid temperature temperature * isub = subgrid_HI_fraction: the subgrid HI fraction (n_HI / n_H) * isub = subgrid_HII_fraction: the subgrid HII fraction (n_HII / n_H) * isub = 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 colibre_subgrid_properties isub) { /* 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 == colibre_compute_subgrid_density) { standard_return = rho_phys; } else if (isub == colibre_compute_subgrid_temperature) { standard_return = exp10(log10_T); } else { /* This should never happen! */ standard_return = -1.; error("Unknown subgrid property to calculate"); } /* 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 */ if (HII_region) { if (isub == colibre_compute_subgrid_density) { const double mu_HII = get_mu_HII_region(cooling, XH); const double P_cgs = (double)P_phys * cooling->colibre_table.pressure_to_cgs; const double rho_cgs = P_cgs * cooling->colibre_table.proton_mass_cgs / (exp10(cooling->colibre_table.log10_kB_cgs) * mu_HII * XH * cooling->HIIregion_temp); return rho_cgs * cooling->colibre_table.density_from_cgs; } else if (isub == colibre_compute_subgrid_temperature) { return cooling->HIIregion_temp; } else { error("Unknown subgrid property to calculate"); return 0.; } } /* 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->colibre_table.number_density_to_cgs; /* Get index along the different table axis for this particle */ int imet, iden, item; float dmet, dden, dtem; cooling_get_index_1d(cooling->colibre_table.Metallicity, colibre_cooling_N_metallicity, logZZsol, &imet, &dmet); cooling_get_index_1d(cooling->colibre_table.nH, colibre_cooling_N_density, log10f(n_H_cgs), &iden, &dden); cooling_get_index_1d(cooling->colibre_table.Temp, colibre_cooling_N_temperature, log10_T, &item, &dtem); float log10_Teq = interpolation_3d(cooling->colibre_table.logTeq, /* */ 0 /*ired*/, imet, iden, /* */ cooling->colibre_table.dz, dmet, dden, /* */ colibre_cooling_N_loaded_redshifts, /* */ colibre_cooling_N_metallicity, /* */ colibre_cooling_N_density); /* Return case 3: * The tabulated thermal equilibrium temperature for the SPH property * is larger than the SPH/EOS gas temperature. This is possible for low * densities and high metallicities, but also if non-equilibrium effects * increase the cooling rate compared to the equilibrium cooling rate. * The subgrid properties cannot be found by projecting the SPH properties * on the thermal equilibrium function and we return the SPH properties */ if (log10_Teq >= log10_T) { return standard_return; } /* Particle is on the EOS and the SPH/EOS temperature is above the thermal * equilibrium temperature. * The subgrid property is calculated by projecting the hydro properties * on the thermal equilibrium temperature where the pressure equals * the hydro pressure */ const float P_cgs = P_phys * cooling->colibre_table.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->colibre_table.logPeq, /* */ 0 /*ired*/, imet, colibre_cooling_N_density - 1, /* */ cooling->colibre_table.dz, dmet, 0., /* */ colibre_cooling_N_loaded_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 == colibre_compute_subgrid_density) { const double rho_cgs = exp10(cooling->colibre_table.nH[colibre_cooling_N_density - 1]) * cooling->colibre_table.proton_mass_cgs / XH; return rho_cgs * cooling->colibre_table.density_from_cgs; } else if (isub == colibre_compute_subgrid_temperature) { const float log10_T_at_Peq = interpolation_3d_no_z( cooling->colibre_table.logTeq, 0 /*ired*/, imet, colibre_cooling_N_density - 1, cooling->colibre_table.dz, dmet, 0., colibre_cooling_N_loaded_redshifts, colibre_cooling_N_metallicity, colibre_cooling_N_density); return exp10(log10_T_at_Peq); } 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_n_at_Peq; float log10_Peq_current = interpolation_3d_no_z( cooling->colibre_table.logPeq, 0 /*ired*/, imet, i, cooling->colibre_table.dz, dmet, 0., colibre_cooling_N_loaded_redshifts, colibre_cooling_N_metallicity, colibre_cooling_N_density); 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->colibre_table.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->colibre_table.nH[i - 1] + delta_P_eq * (cooling->colibre_table.nH[i] - cooling->colibre_table.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 == colibre_compute_subgrid_density) { const double rho_cgs = exp10(log10_n_at_Peq) * cooling->colibre_table.proton_mass_cgs / XH; /* only project to higher densities */ if (rho_cgs * cooling->colibre_table.density_from_cgs < rho_phys) { return standard_return; } else { return rho_cgs * cooling->colibre_table.density_from_cgs; } } else if (isub == colibre_compute_subgrid_temperature) { /* We found a valid density: Get the index along the density * axis of the tables. */ int iden_eq; float dden_eq; cooling_get_index_1d(cooling->colibre_table.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->colibre_table.logTeq, /* */ 0 /*ired*/, imet, iden_eq, /* */ cooling->colibre_table.dz, dmet, dden_eq, /* */ colibre_cooling_N_loaded_redshifts, /* */ colibre_cooling_N_metallicity, /* */ colibre_cooling_N_density); /* only project to lower temperatures */ if (log10_T_at_Peq > log10_T) { return standard_return; } else { return exp10(log10_T_at_Peq); } } else { error("Unknown subgrid property to calculate"); return 0.; } } 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_COLIBRE_COOLING_SUBGRID_H */