/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2021 Tsang Keung Chan (chantsangkeung@gmail.com)
*
* 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 .
*
******************************************************************************/
/**
* @file src/rt/SPHM1RT/rt_cooling.c
* @brief SPHM1RT cooling functions
*/
/* Some standard headers. */
#include
#include /* access to CVDls interface */
#include
#include
#include
/* Local includes. */
#include "rt_cooling.h"
#include "rt_cooling_rates.h"
#include "rt_getters.h"
#include "rt_setters.h"
/**
* @brief Main function for the thermochemistry step.
*
* @param p Particle to work on.
* @param xp Pointer to the particle' extended data.
* @param rt_props RT properties struct
* @param cosmo The current cosmological model.
* @param hydro_props The #hydro_props.
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
* @param dt The time-step of this particle.
*/
void rt_do_thermochemistry(struct part* restrict p, struct xpart* restrict xp,
struct rt_props* rt_props,
const struct cosmology* restrict cosmo,
const struct hydro_props* hydro_props,
const struct phys_const* restrict phys_const,
const struct unit_system* restrict us,
const double dt) {
/* Nothing to do here? */
if (rt_props->skip_thermochemistry == 1) return;
if (dt == 0.0) return;
rt_check_unphysical_elem_spec(p, rt_props);
struct rt_part_data* rpd = &p->rt_data;
struct RTUserData data; /* data for CVODE */
const double dt_cgs = dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME);
const double conv_factor_internal_energy_to_cgs =
units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
const double conv_factor_opacity_from_cgs =
units_cgs_conversion_factor(us, UNIT_CONV_MASS) /
units_cgs_conversion_factor(us, UNIT_CONV_LENGTH) /
units_cgs_conversion_factor(us, UNIT_CONV_LENGTH);
// TK remark: here we only consider the on-the-spot approximation
/**************************/
/* INITIALIZATION */
/**************************/
int useparams = rt_props->useparams;
/* adopt on the spot approximation (OTSA) by default */
/* TODO: currently there is no non-OTSA implementation; to do in the future */
int onthespot = rt_props->onthespot;
data.onthespot = onthespot;
int coolingon = rt_props->coolingon;
data.coolingon = coolingon;
int fixphotondensity = rt_props->fixphotondensity;
data.fixphotondensity = fixphotondensity;
double metal_mass_fraction[rt_chemistry_element_count];
for (int elem = 0; elem < rt_chemistry_element_count; elem++) {
metal_mass_fraction[elem] = (double)(rpd->tchem.metal_mass_fraction[elem]);
data.metal_mass_fraction[elem] = metal_mass_fraction[elem];
}
const double X_H = metal_mass_fraction[rt_chemistry_element_H];
const double m_H_cgs = phys_const->const_proton_mass *
units_cgs_conversion_factor(us, UNIT_CONV_MASS);
const double proton_mass_cgs_inv = 1.0 / m_H_cgs;
data.m_H_cgs = m_H_cgs;
const double k_B_cgs = phys_const->const_boltzmann_k *
units_cgs_conversion_factor(us, UNIT_CONV_ENERGY) /
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
data.k_B_cgs = k_B_cgs;
const double cred_phys = rt_get_physical_cred(p, cosmo->a);
const double cred_cgs =
cred_phys * units_cgs_conversion_factor(us, UNIT_CONV_VELOCITY);
data.cred_cgs = cred_cgs;
/* Get particle density [ and convert to g * cm^-3] */
const double rho = hydro_get_physical_density(p, cosmo);
double rho_cgs = rho * units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
data.rho_cgs = rho_cgs;
/* Hydrogen number density (X_H * rho / m_p) [cm^-3] */
const double n_H_cgs = X_H * rho_cgs * proton_mass_cgs_inv;
data.n_H_cgs = n_H_cgs;
/* Current energy (in internal units) */
float urad[RT_NGROUPS];
rt_get_physical_urad_multifrequency(p, cosmo, urad);
/* need to convert to cgs */
double ngamma_cgs[3];
/* for now, the 0th bin for urad is 0-HI, so we ignore it */
for (int g = 0; g < 3; g++) {
ngamma_cgs[g] =
(double)(rho_cgs * urad[g + 1] * conv_factor_internal_energy_to_cgs /
rt_props->ionizing_photon_energy_cgs[g]);
data.ngamma_cgs[g] = ngamma_cgs[g];
}
/* overwrite the photon density if we choose to fix it */
for (int i = 0; i < 3; i++) {
if ((rt_props->Fgamma_fixed_cgs[i] > 0.0) &&
(rt_props->fixphotondensity == 1)) {
ngamma_cgs[i] = rt_props->Fgamma_fixed_cgs[i] / cred_cgs;
data.ngamma_cgs[i] = ngamma_cgs[i];
urad[i + 1] = (float)(data.ngamma_cgs[i] / rho_cgs /
conv_factor_internal_energy_to_cgs *
rt_props->ionizing_photon_energy_cgs[i]);
}
}
double abundances[rt_species_count];
for (int spec = 0; spec < rt_species_count; spec++) {
abundances[spec] = (double)(rpd->tchem.abundances[spec]);
data.abundances[spec] = abundances[spec];
}
const double u = hydro_get_physical_internal_energy(p, xp, cosmo);
double u_cgs = u * conv_factor_internal_energy_to_cgs;
double T_cgs = rt_convert_u_to_temp(k_B_cgs, m_H_cgs, X_H, u_cgs, abundances);
double T_min_cgs = hydro_props->minimal_temperature;
double u_min_cgs =
rt_convert_temp_to_u(k_B_cgs, m_H_cgs, T_min_cgs, X_H, abundances);
u_cgs = fmax(u_cgs, u_min_cgs);
data.u_cgs = u_cgs;
data.u_min_cgs = u_min_cgs;
/**************************/
/* GET RATE COEFFICIENTS */
/**************************/
int aindex[3];
rt_get_index_to_species(aindex);
for (int i = 0; i < 3; i++) {
data.aindex[i] = aindex[i];
}
double alphalist[rt_species_count], betalist[rt_species_count],
Gammalist[rt_species_count], sigmalist[3][3], epsilonlist[3][3];
if (useparams == 1) {
betalist[rt_sp_elec] = 0.0;
betalist[rt_sp_HI] = rt_props->beta_cgs_H;
betalist[rt_sp_HII] = 0.0;
betalist[rt_sp_HeI] = 0.0;
betalist[rt_sp_HeII] = 0.0;
betalist[rt_sp_HeIII] = 0.0;
alphalist[rt_sp_elec] = 0.0;
alphalist[rt_sp_HI] = 0.0;
alphalist[rt_sp_HeI] = 0.0;
if (onthespot == 1) {
alphalist[rt_sp_HII] = rt_props->alphaB_cgs_H;
alphalist[rt_sp_HeII] = 0.0;
alphalist[rt_sp_HeIII] = 0.0;
} else {
alphalist[rt_sp_HII] = rt_props->alphaA_cgs_H;
alphalist[rt_sp_HeII] = 0.0;
alphalist[rt_sp_HeIII] = 0.0;
}
sigmalist[0][0] = rt_props->sigma_cross_cgs_H[0];
sigmalist[1][0] = rt_props->sigma_cross_cgs_H[1];
sigmalist[2][0] = rt_props->sigma_cross_cgs_H[2];
sigmalist[0][1] = 0.0;
sigmalist[1][1] = 0.0;
sigmalist[2][1] = 0.0;
sigmalist[0][2] = 0.0;
sigmalist[1][2] = 0.0;
sigmalist[2][2] = 0.0;
data.alphaA_cgs_H = rt_props->alphaA_cgs_H;
data.alphaB_cgs_H = rt_props->alphaB_cgs_H;
data.beta_cgs_H = rt_props->beta_cgs_H;
data.sigma_cross_cgs_H[0] = rt_props->sigma_cross_cgs_H[0];
data.sigma_cross_cgs_H[1] = rt_props->sigma_cross_cgs_H[1];
data.sigma_cross_cgs_H[2] = rt_props->sigma_cross_cgs_H[2];
for (int spec = 0; spec < rt_species_count; spec++) {
Gammalist[spec] = 0.0;
}
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
epsilonlist[i][j] = 0.0;
}
}
} else {
rt_compute_rate_coefficients(T_cgs, onthespot, alphalist, betalist,
Gammalist, sigmalist, epsilonlist);
}
data.useparams = rt_props->useparams;
/**************************/
/* SOLVING RATE EQUTAION */
/**************************/
/* Try explicit solution */
double new_abundances[rt_species_count], finish_abundances[rt_species_count],
max_relative_change, new_ngamma_cgs[3], u_new_cgs;
max_relative_change = 0.0;
/* compute net changes and cooling and heating for explicit solution */
rt_compute_explicit_thermochemistry_solution(
n_H_cgs, cred_cgs, dt_cgs, rho_cgs, u_cgs, u_min_cgs, abundances,
ngamma_cgs, alphalist, betalist, Gammalist, sigmalist, epsilonlist,
aindex, &u_new_cgs, new_abundances, new_ngamma_cgs, &max_relative_change);
/* check whether xHI bigger than one */
int errorHI = 0;
if (new_abundances[rt_sp_HI] > 1.01) {
errorHI = 1;
} else {
rt_enforce_constraint_equations(new_abundances, metal_mass_fraction,
finish_abundances);
}
if ((max_relative_change < rt_props->explicitRelTolerance) &&
(errorHI == 0)) {
for (int spec = 0; spec < rt_species_count; spec++) {
if (finish_abundances[spec] > 0.f) {
if (finish_abundances[spec] < FLT_MAX) {
rpd->tchem.abundances[spec] = (float)(finish_abundances[spec]);
} else {
error("finish_abundances larger than FLT_MAX");
}
} else {
rpd->tchem.abundances[spec] = 0.f;
}
}
if (coolingon == 1) {
float u_new = 0.0f;
if (u_new_cgs / conv_factor_internal_energy_to_cgs > 0.f) {
if (u_new_cgs / conv_factor_internal_energy_to_cgs < FLT_MAX) {
u_new = (float)(u_new_cgs / conv_factor_internal_energy_to_cgs);
}
}
hydro_set_physical_internal_energy(p, xp, cosmo, u_new);
}
/* set radiation energy */
float urad_new[RT_NGROUPS];
urad_new[0] = 0.f;
if (fixphotondensity == 0) {
for (int i = 0; i < 3; i++) {
urad_new[i + 1] = 0.f;
if (new_ngamma_cgs[i] / rho_cgs / conv_factor_internal_energy_to_cgs *
rt_props->ionizing_photon_energy_cgs[i] >
0.f) {
if (new_ngamma_cgs[i] / rho_cgs / conv_factor_internal_energy_to_cgs *
rt_props->ionizing_photon_energy_cgs[i] <
FLT_MAX) {
urad_new[i + 1] = (float)(new_ngamma_cgs[i] / rho_cgs /
conv_factor_internal_energy_to_cgs *
rt_props->ionizing_photon_energy_cgs[i]);
}
}
}
} else {
for (int i = 0; i < 3; i++) {
urad_new[i + 1] = urad[i + 1];
}
}
rt_set_physical_urad_multifrequency(p, cosmo, urad_new);
/* chi is in physical unit (L^2/M) */
float chi_new[RT_NGROUPS];
for (int i = 0; i < RT_NGROUPS; i++) {
chi_new[i] = 0.0f;
}
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
if (finish_abundances[aindex[j]] * n_H_cgs / rho_cgs * sigmalist[i][j] *
conv_factor_opacity_from_cgs >
0.f) {
if (finish_abundances[aindex[j]] * n_H_cgs / rho_cgs *
sigmalist[i][j] * conv_factor_opacity_from_cgs <
FLT_MAX) {
chi_new[i + 1] +=
(float)(finish_abundances[aindex[j]] * n_H_cgs / rho_cgs *
sigmalist[i][j] * conv_factor_opacity_from_cgs);
}
}
}
}
rt_set_physical_radiation_opacity(p, cosmo, chi_new);
rt_check_unphysical_elem_spec(p, rt_props);
return;
} else {
/**************************************
* Explicit solution is insufficient. *
* Use implicit solver. *
**************************************/
realtype reltol, t;
N_Vector abstol_vector, y;
int maxsteps = 100000;
int network_size, icount = 0;
/* 3 for species; */
network_size = 3;
/* 1 for thermal energy; */
if (coolingon == 1) {
network_size += 1;
}
/* 3 for radiation bins */
if (fixphotondensity == 0) {
network_size += 3;
}
y = N_VNew_Serial(network_size);
abstol_vector = N_VNew_Serial(network_size);
for (int i = 0; i < 3; i++) {
NV_Ith_S(y, icount) = (realtype)data.abundances[aindex[i]];
NV_Ith_S(abstol_vector, icount) = (realtype)(rt_props->absoluteTolerance);
icount += 1;
}
if (coolingon == 1) {
NV_Ith_S(y, icount) = (realtype)u_cgs;
NV_Ith_S(abstol_vector, icount) = (realtype)rt_props->absoluteTolerance;
icount += 1;
}
if (fixphotondensity == 0) {
for (int i = 0; i < 3; i++) {
NV_Ith_S(y, icount) = (realtype)data.ngamma_cgs[i];
NV_Ith_S(abstol_vector, icount) = (realtype)rt_props->absoluteTolerance;
icount += 1;
}
}
/* Set up the solver */
/* Set the tolerances*/
reltol = (realtype)rt_props->relativeTolerance;
/* Use CVodeCreate to create the solver
* memory and specify the Backward Differentiation
* Formula. Note that CVODE now uses Newton iteration
* iteration by default, so no need to specify this. */
void* cvode_mem;
cvode_mem = CVodeCreate(CV_BDF);
/* Set the user data for CVode */
CVodeSetUserData(cvode_mem, &data);
/* Use CVodeSetMaxNumSteps to set the maximum number
* of steps CVode takes. */
CVodeSetMaxNumSteps(cvode_mem, maxsteps);
/* Use CVodeInit to initialise the integrator
* memory and specify the right hand side
* function in y' = f(t,y) (i.e. the rate
* equations), the initial time 0.0 and the
* initial conditions, in y. */
CVodeInit(cvode_mem, rt_frateeq, 0.0f, y);
/* Use CVodeSVtolerances to specify the scalar
* relative and absolute tolerances. */
CVodeSVtolerances(cvode_mem, reltol, abstol_vector);
/* Create a dense SUNMatrix to use in the
* linear solver. */
SUNMatrix A_sun;
A_sun = SUNDenseMatrix(network_size, network_size);
/* Create a denst SUNLinearSolver object
* to use in CVode. */
SUNLinearSolver LS_sun;
LS_sun = SUNDenseLinearSolver(y, A_sun);
/* Attach the matrix and linear
* solver to CVode. */
CVDlsSetLinearSolver(cvode_mem, LS_sun, A_sun);
/* Specify the maximum number of convergence
* test failures. */
CVodeSetMaxConvFails(cvode_mem, 5000);
/* Call CVode() to integrate the chemistry. */
CVode(cvode_mem, (realtype)dt_cgs, y, &t, CV_NORMAL);
/* Write the output abundances to the gas cell
* Note that species not included in the reduced
* network are kept constant in the GasVars struct. */
icount = 0;
for (int i = 0; i < 3; i++) {
new_abundances[aindex[i]] = (double)NV_Ith_S(y, icount);
icount += 1;
}
if (coolingon == 1) {
u_cgs = (double)NV_Ith_S(y, icount);
icount += 1;
}
if (fixphotondensity == 0) {
for (int i = 0; i < 3; i++) {
new_ngamma_cgs[i] = (double)NV_Ith_S(y, icount);
icount += 1;
}
}
if (new_abundances[rt_sp_HI] > 1.01)
error("HI fraction bigger than one after the CVODE solver");
rt_enforce_constraint_equations(new_abundances, metal_mass_fraction,
finish_abundances);
for (int spec = 0; spec < rt_species_count; spec++) {
if (finish_abundances[spec] > 0.f) {
if (finish_abundances[spec] < FLT_MAX) {
rpd->tchem.abundances[spec] = (float)(finish_abundances[spec]);
} else {
error("finish_abundances larger than FLT_MAX");
}
} else {
rpd->tchem.abundances[spec] = 0.f;
}
}
if (coolingon == 1) {
float u_new = 0.0f;
if (u_new_cgs / conv_factor_internal_energy_to_cgs > 0.f) {
if (u_new_cgs / conv_factor_internal_energy_to_cgs < FLT_MAX) {
u_new = (float)(u_new_cgs / conv_factor_internal_energy_to_cgs);
}
}
hydro_set_physical_internal_energy(p, xp, cosmo, u_new);
}
/* set radiation energy */
float urad_new[RT_NGROUPS];
urad_new[0] = 0.f;
if (fixphotondensity == 0) {
for (int i = 0; i < 3; i++) {
urad_new[i + 1] = 0.f;
if (new_ngamma_cgs[i] / rho_cgs / conv_factor_internal_energy_to_cgs *
rt_props->ionizing_photon_energy_cgs[i] >
0.f) {
if (new_ngamma_cgs[i] / rho_cgs / conv_factor_internal_energy_to_cgs *
rt_props->ionizing_photon_energy_cgs[i] <
FLT_MAX) {
urad_new[i + 1] = (float)(new_ngamma_cgs[i] / rho_cgs /
conv_factor_internal_energy_to_cgs *
rt_props->ionizing_photon_energy_cgs[i]);
}
}
}
} else {
for (int i = 0; i < 3; i++) {
urad_new[i + 1] = urad[i + 1];
}
}
rt_set_physical_urad_multifrequency(p, cosmo, urad_new);
/* chi is in physical unit (L^2/M) */
float chi_new[RT_NGROUPS];
for (int i = 0; i < RT_NGROUPS; i++) {
chi_new[i] = 0.0f;
}
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
if (finish_abundances[aindex[j]] * n_H_cgs / rho_cgs * sigmalist[i][j] *
conv_factor_opacity_from_cgs >
0.f) {
if (finish_abundances[aindex[j]] * n_H_cgs / rho_cgs *
sigmalist[i][j] * conv_factor_opacity_from_cgs <
FLT_MAX) {
chi_new[i + 1] +=
(float)(finish_abundances[aindex[j]] * n_H_cgs / rho_cgs *
sigmalist[i][j] * conv_factor_opacity_from_cgs);
}
}
}
}
rt_set_physical_radiation_opacity(p, cosmo, chi_new);
SUNLinSolFree(LS_sun);
SUNMatDestroy(A_sun);
N_VDestroy_Serial(y);
N_VDestroy_Serial(abstol_vector);
CVodeFree(&cvode_mem);
rt_check_unphysical_elem_spec(p, rt_props);
}
}
/**
* @brief Do the thermochemistry on a particle.
*
* @param p Particle to work on.
* @param xp Pointer to the particle' extended data.
* @param rt_props RT properties struct
* @param cosmo The current cosmological model.
* @param hydro_props The #hydro_props.
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
* @param dt The time-step of this particle.
*/
void rt_tchem(struct part* restrict p, struct xpart* restrict xp,
struct rt_props* rt_props, const struct cosmology* restrict cosmo,
const struct hydro_props* hydro_props,
const struct phys_const* restrict phys_const,
const struct unit_system* restrict us, const double dt) {
rt_do_thermochemistry(p, xp, rt_props, cosmo, hydro_props, phys_const, us,
dt);
}