/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2021 Mladen Ivkovic (mladen.ivkovic@hotmail.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 .
*
******************************************************************************/
#ifndef SWIFT_RT_PROPERTIES_GEAR_H
#define SWIFT_RT_PROPERTIES_GEAR_H
#include "rt_grackle_utils.h"
#include "rt_interaction_rates.h"
#include "rt_parameters.h"
#include "rt_stellar_emission_model.h"
#include
/**
* @file src/rt/GEAR/rt_properties.h
* @brief Main header file for the 'GEAR' radiative transfer scheme
* properties.
*/
#define RT_IMPLEMENTATION "GEAR M1closure"
#if defined(RT_RIEMANN_SOLVER_GLF)
#define RT_RIEMANN_SOLVER_NAME "GLF Riemann Solver"
#elif defined(RT_RIEMANN_SOLVER_HLL)
#define RT_RIEMANN_SOLVER_NAME "HLL Riemann Solver"
#else
#error "No valid choice of RT Riemann solver has been selected"
#endif
/**
* @brief Properties of the 'GEAR' radiative transfer model
*/
struct rt_props {
/* Which stellar emission model to use */
enum rt_stellar_emission_models stellar_emission_model;
/* (Lower) frequency bin edges for photon groups */
float photon_groups[RT_NGROUPS];
/* Global constant stellar emission rates */
double stellar_const_emission_rates[RT_NGROUPS];
/* CFL condition */
float CFL_condition;
/* Factor to limit cooling time by */
float f_limit_cooling_time;
/* do we set initial ionization mass fractions manually? */
int set_initial_ionization_mass_fractions;
int set_equilibrium_initial_ionization_mass_fractions;
/* initial mass fractions for ionization species */
/* the following are required for manually setting exact values */
float mass_fraction_HI_init;
float mass_fraction_HII_init;
float mass_fraction_HeI_init;
float mass_fraction_HeII_init;
float mass_fraction_HeIII_init;
/* float number_density_electrons_init; [> todo: do we need this? <] */
/* Hydrogen and Helium mass fractions of the non-metal portion of the gas */
float hydrogen_mass_fraction;
float helium_mass_fraction;
/* Skip thermochemistry? For testing/debugging only! */
int skip_thermochemistry;
/* Re-do thermochemistry recursively if difference in internal energy is too
* big? */
int max_tchem_recursion;
/* Optionally restrict maximal timestep for stars */
float stars_max_timestep;
/* Which stellar spectrum type to use? */
int stellar_spectrum_type;
/* If constant: get max frequency */
double const_stellar_spectrum_max_frequency;
/* If blackbody: get temperature */
double stellar_spectrum_blackbody_T;
/* Storage for integrated photoionization cross sections */
/* Note: they are always in cgs. */
double** energy_weighted_cross_sections;
double** number_weighted_cross_sections;
/* Mean photon energy in frequency bin for user provided spectrum. In erg.*/
double average_photon_energy[RT_NGROUPS];
/* Integral over photon numbers of user provided spectrum. */
double photon_number_integral[RT_NGROUPS];
/* Grackle Stuff */
/* ------------- */
/*! grackle unit system */
code_units grackle_units;
/*! grackle chemistry data */
chemistry_data grackle_chemistry_data;
/*! grackle chemistry data storage
* (needed for local function calls) */
chemistry_data_storage grackle_chemistry_rates;
/*! use case B recombination? */
int case_B_recombination;
/*! make grackle talkative? */
int grackle_verbose;
#ifdef SWIFT_RT_DEBUG_CHECKS
/* radiation emitted by stars this step. This is not really a property,
* but a placeholder to sum up a global variable. It's being reset
* every timestep. */
unsigned long long debug_radiation_emitted_this_step;
/* total radiation emitted by stars. This is not really a property,
* but a placeholder to sum up a global variable */
unsigned long long debug_radiation_emitted_tot;
/* radiation absorbed by gas this step. This is not really a property,
* but a placeholder to sum up a global variable */
unsigned long long debug_radiation_absorbed_this_step;
/* total radiation absorbed by gas. This is not really a property,
* but a placeholder to sum up a global variable */
unsigned long long debug_radiation_absorbed_tot;
/* Max number of subcycles per hydro step */
int debug_max_nr_subcycles;
#endif
};
/* Some declarations to avoid cyclical inclusions. */
/* ----------------------------------------------- */
/* Keep the declarations for *after* the definition of rt_props struct */
/**
* @brief allocate and pre-compute the averaged cross sections
* for each photon group and ionizing species.
* Declare this here to avoid cyclical inclusions.
*
* @param rt_props RT properties struct
* @param phys_const physical constants struct
* @param us internal units struct
**/
void rt_cross_sections_init(struct rt_props* restrict rt_props,
const struct phys_const* restrict phys_const,
const struct unit_system* restrict us);
/* Now for the good stuff */
/* ------------------------------------- */
/**
* @brief Print the RT model.
*
* @param rtp The #rt_props
*/
__attribute__((always_inline)) INLINE static void rt_props_print(
const struct rt_props* rtp) {
/* Only the master print */
if (engine_rank != 0) return;
message("Radiative transfer scheme: '%s'", RT_IMPLEMENTATION);
message("RT Riemann Solver used: '%s'", RT_RIEMANN_SOLVER_NAME);
char messagestring[200] = "Using photon frequency bins: [ ";
char freqstring[20];
for (int g = 0; g < RT_NGROUPS; g++) {
sprintf(freqstring, "%.3g ", rtp->photon_groups[g]);
strcat(messagestring, freqstring);
}
strcat(messagestring, "]");
message("%s", messagestring);
if (rtp->stellar_emission_model == rt_stellar_emission_model_const) {
strcpy(messagestring, "Using constant stellar emission rates: [ ");
for (int g = 0; g < RT_NGROUPS; g++) {
sprintf(freqstring, "%.3g ", rtp->stellar_const_emission_rates[g]);
strcat(messagestring, freqstring);
}
strcat(messagestring, "]");
message("%s", messagestring);
} else if (rtp->stellar_emission_model ==
rt_stellar_emission_model_IlievTest) {
message("Using Iliev+06 Test 4 stellar emission model.");
} else {
error("Unknown stellar emission model %d", rtp->stellar_emission_model);
}
if (rtp->set_equilibrium_initial_ionization_mass_fractions)
message(
"Setting initial ionization mass fractions "
"assuming ionization equilibrium");
if (rtp->set_initial_ionization_mass_fractions)
message(
"Using initial ionization mass fractions specified in parameter file");
if (rtp->skip_thermochemistry)
message("WARNING: Thermochemistry will be skipped.");
}
/**
* @brief Initialize the global properties of the RT scheme.
*
* @param rtp The #rt_props.
* @param phys_const The physical constants in the internal unit system.
* @param us The internal unit system.
* @param params The parsed parameters.
* @param cosmo The cosmological model.
*/
__attribute__((always_inline)) INLINE static void rt_props_init(
struct rt_props* rtp, const struct phys_const* phys_const,
const struct unit_system* us, struct swift_params* params,
struct cosmology* cosmo) {
/* Make sure we reset debugging counters correctly after
* zeroth step. */
/* Read in photon frequency group properties */
/* ----------------------------------------- */
if (RT_NGROUPS <= 0) {
error(
"You need to run GEAR-RT with at least 1 photon group, "
"you have %d",
RT_NGROUPS);
} else {
/* !! Keep the frequencies in Hz for now. !! */
parser_get_param_float_array(params, "GEARRT:photon_groups_Hz", RT_NGROUPS,
rtp->photon_groups);
}
/* Sanity check: photon group edges must be in increasing order. */
for (int g = 0; g < RT_NGROUPS - 1; g++) {
if (rtp->photon_groups[g + 1] <= rtp->photon_groups[g])
error(
"Photon frequency bin edges need to be in increasing order. "
"Found index %d %.3g <= %.3g",
g + 1, rtp->photon_groups[g + 1], rtp->photon_groups[g]);
}
/* Get stellar emission rate model related parameters */
/* -------------------------------------------------- */
/* First initialize everything */
for (int g = 0; g < RT_NGROUPS; g++) {
rtp->stellar_const_emission_rates[g] = 0.;
}
rtp->stellar_emission_model = rt_stellar_emission_model_none;
char stellar_model_str[80];
parser_get_param_string(params, "GEARRT:stellar_luminosity_model",
stellar_model_str);
if (strcmp(stellar_model_str, "const") == 0) {
rtp->stellar_emission_model = rt_stellar_emission_model_const;
} else if (strcmp(stellar_model_str, "IlievTest4") == 0) {
rtp->stellar_emission_model = rt_stellar_emission_model_IlievTest;
} else {
error("Unknown stellar luminosity model '%s'", stellar_model_str);
}
if (rtp->stellar_emission_model == rt_stellar_emission_model_const) {
/* Read the luminosities from the parameter file */
double emission_rates[RT_NGROUPS];
parser_get_param_double_array(params,
"GEARRT:const_stellar_luminosities_LSol",
RT_NGROUPS, emission_rates);
const double unit_power = units_cgs_conversion_factor(us, UNIT_CONV_POWER);
const double unit_power_inv = 1. / unit_power;
for (int g = 0; g < RT_NGROUPS; g++) {
rtp->stellar_const_emission_rates[g] = emission_rates[g] * unit_power_inv;
}
} else if (rtp->stellar_emission_model ==
rt_stellar_emission_model_IlievTest) {
/* Nothing to do here */
} else {
error("Unknown stellar emission model %d", rtp->stellar_emission_model);
}
/* get reduced speed of light factor */
/* --------------------------------- */
const float f_r = parser_get_param_float(params, "GEARRT:f_reduce_c");
if (f_r <= 0.f)
error("Invalid speed of light reduction factor: %.3e <= 0.", f_r);
rt_params.reduced_speed_of_light = phys_const->const_speed_light_c * f_r;
rt_params.reduced_speed_of_light_inverse =
1.f / rt_params.reduced_speed_of_light;
/* get CFL condition */
/* ----------------- */
const float CFL = parser_get_param_float(params, "GEARRT:CFL_condition");
if (CFL <= 0.f) error("Invalid CFL number: %.3e <= 0.", CFL);
rtp->CFL_condition = CFL;
const float f_limit_cooling_time = parser_get_opt_param_float(
params, "GEARRT:f_limit_cooling_time", /*default=*/0.6);
if (f_limit_cooling_time < 0.f)
error("Invalid cooling time reduction factor: %.3e < 0.",
f_limit_cooling_time);
else if (f_limit_cooling_time == 0.f)
message("Warning: Computation of cooling time will be skipped");
rtp->f_limit_cooling_time = f_limit_cooling_time;
/* Get thermochemistry set-up */
/* -------------------------- */
rtp->hydrogen_mass_fraction =
parser_get_param_float(params, "GEARRT:hydrogen_mass_fraction");
rtp->helium_mass_fraction = 1.f - rtp->hydrogen_mass_fraction;
if (rtp->hydrogen_mass_fraction <= 0.f || rtp->hydrogen_mass_fraction > 1.f)
error("Invalid hydrogen mass fraction: %g", rtp->hydrogen_mass_fraction);
/* Are we manually overwriting initial mass fractions of H and He? */
rtp->set_initial_ionization_mass_fractions = parser_get_opt_param_int(
params, "GEARRT:set_initial_ionization_mass_fractions",
/* default = */ 0);
if (rtp->set_initial_ionization_mass_fractions) {
/* Read in mass fractions */
rtp->mass_fraction_HI_init =
parser_get_param_float(params, "GEARRT:mass_fraction_HI");
rtp->mass_fraction_HII_init =
parser_get_param_float(params, "GEARRT:mass_fraction_HII");
rtp->mass_fraction_HeI_init =
parser_get_param_float(params, "GEARRT:mass_fraction_HeI");
rtp->mass_fraction_HeII_init =
parser_get_param_float(params, "GEARRT:mass_fraction_HeII");
rtp->mass_fraction_HeIII_init =
parser_get_param_float(params, "GEARRT:mass_fraction_HeIII");
/* Temporary check neglecting metals. Make sure we sum up to 1. */
const float h_sum =
rtp->mass_fraction_HI_init + rtp->mass_fraction_HII_init;
if (fabsf(h_sum - rtp->hydrogen_mass_fraction) > 1e-4)
error(
"Inconsistent H mass fractions: XH_tot %.6g != XHI %.6g + XHII %.6g",
rtp->hydrogen_mass_fraction, rtp->mass_fraction_HI_init,
rtp->mass_fraction_HII_init);
const float he_sum = rtp->mass_fraction_HeI_init +
rtp->mass_fraction_HeII_init +
rtp->mass_fraction_HeIII_init;
if (fabsf(he_sum - rtp->helium_mass_fraction) > 1e-4)
error(
"Inconsistent He mass fractions: XHe_tot %.6g != XHeI %.6g + XHeII "
"%.6g + XHeIII %.6g",
rtp->helium_mass_fraction, rtp->mass_fraction_HeI_init,
rtp->mass_fraction_HeII_init, rtp->mass_fraction_HeIII_init);
const float mass_fraction_sum = h_sum + he_sum;
if (fabsf(mass_fraction_sum - 1.f) > 1e-5)
error("Constituent species mass fraction sums up to %.6f, I expect 1.0",
mass_fraction_sum);
} else {
/* Initialize properties to deliberately bogus values */
rtp->mass_fraction_HI_init = -1.f;
rtp->mass_fraction_HII_init = -1.f;
rtp->mass_fraction_HeI_init = -1.f;
rtp->mass_fraction_HeII_init = -1.f;
rtp->mass_fraction_HeIII_init = -1.f;
}
/* Are we setting up initial mass fractions in equilibrium? */
rtp->set_equilibrium_initial_ionization_mass_fractions =
parser_get_opt_param_int(
params, "GEARRT:set_equilibrium_initial_ionization_mass_fractions",
/* default = */ 0);
if (rtp->set_equilibrium_initial_ionization_mass_fractions &&
rtp->set_initial_ionization_mass_fractions)
error(
"Can't use equilibrium initial ionization mass fractions "
"simultaneously with manually set mass fractions. Pick one.");
/* Are we skipping thermochemistry? */
rtp->skip_thermochemistry = parser_get_opt_param_int(
params, "GEARRT:skip_thermochemistry", /* default = */ 0);
/* Are we re-doing thermochemistry? */
rtp->max_tchem_recursion = parser_get_opt_param_int(
params, "GEARRT:max_tchem_recursion", /* default = */ 0);
/* Stellar Spectra */
/* --------------- */
/* Initialize conditional parameters to bogus values */
rtp->const_stellar_spectrum_max_frequency = -1.;
rtp->stellar_spectrum_blackbody_T = -1.;
rtp->stellar_spectrum_type =
parser_get_param_int(params, "GEARRT:stellar_spectrum_type");
if (rtp->stellar_spectrum_type == 0) {
/* Constant spectrum: Read additional parameter */
/* TODO: also translate back to internal units at later. For now, keep it in
* Hz */
rtp->const_stellar_spectrum_max_frequency = parser_get_param_float(
params, "GEARRT:stellar_spectrum_const_max_frequency_Hz");
} else if (rtp->stellar_spectrum_type == 1) {
/* Blackbody spectrum: Read additional parameter */
rtp->stellar_spectrum_blackbody_T = parser_get_param_float(
params, "GEARRT:stellar_spectrum_blackbody_temperature_K");
rtp->stellar_spectrum_blackbody_T /=
units_cgs_conversion_factor(us, UNIT_CONV_TEMPERATURE);
} else {
error("Selected unknown stellar spectrum type %d",
rtp->stellar_spectrum_type);
}
/* Maximal Star Timestep? */
/* ---------------------- */
rtp->stars_max_timestep = parser_get_opt_param_float(
params, "GEARRT:stars_max_timestep", /*default=*/FLT_MAX);
/* Turn off if negative value given */
if (rtp->stars_max_timestep < 0.f) rtp->stars_max_timestep = FLT_MAX;
/* Better safe than sorry */
if (rtp->stars_max_timestep == 0.f)
error("You are restricting star time step to 0. That's a no-no.");
#ifdef SWIFT_RT_DEBUG_CHECKS
rtp->debug_radiation_emitted_tot = 0ULL;
rtp->debug_radiation_emitted_this_step = 0ULL;
rtp->debug_radiation_absorbed_tot = 0ULL;
rtp->debug_radiation_absorbed_this_step = 0ULL;
/* Don't make it an optional parameter here so we crash
* if I forgot to provide it */
rtp->debug_max_nr_subcycles =
parser_get_param_int(params, "TimeIntegration:max_nr_rt_subcycles");
#endif
/* Grackle setup */
/* ------------- */
rtp->grackle_verbose =
parser_get_opt_param_int(params, "GEARRT:grackle_verbose", /*default=*/0);
rtp->case_B_recombination = parser_get_opt_param_int(
params, "GEARRT:case_B_recombination", /*default=*/1);
rt_init_grackle(&rtp->grackle_units, &rtp->grackle_chemistry_data,
&rtp->grackle_chemistry_rates, rtp->hydrogen_mass_fraction,
rtp->grackle_verbose, rtp->case_B_recombination, us, cosmo);
/* Pre-compute interaction rates/cross sections */
/* -------------------------------------------- */
rtp->energy_weighted_cross_sections = NULL;
rtp->number_weighted_cross_sections = NULL;
rt_cross_sections_init(rtp, phys_const, us);
/* Finishers */
/* --------- */
}
__attribute__((always_inline)) INLINE static void rt_props_update(
struct rt_props* rtp, const struct unit_system* us,
struct cosmology* cosmo) {
update_grackle_units_cosmo(&(rtp->grackle_units), us, cosmo);
}
/**
* @brief Write an RT properties struct to the given FILE as a
* stream of bytes.
*
* @param props the struct
* @param stream the file stream
*/
__attribute__((always_inline)) INLINE static void rt_struct_dump(
const struct rt_props* props, FILE* stream) {
restart_write_blocks((void*)props, sizeof(struct rt_props), 1, stream,
"RT props", "RT properties struct");
/* The RT parameters, in particular the reduced speed of light, are
* not defined at compile time. So we need to read them in again. */
restart_write_blocks(&rt_params, sizeof(struct rt_parameters), 1, stream,
"RT global parameters", "RT global parameters struct");
}
/**
* @brief Restore an RT properties struct from the given FILE as
* a stream of bytes.
*
* @param props the struct
* @param stream the file stream
* @param phys_const The physical constants in the internal unit system.
* @param us The internal unit system.
* @param cosmo the #cosmology
*/
__attribute__((always_inline)) INLINE static void rt_struct_restore(
struct rt_props* props, FILE* stream, const struct phys_const* phys_const,
const struct unit_system* us, const struct cosmology* restrict cosmo) {
restart_read_blocks((void*)props, sizeof(struct rt_props), 1, stream, NULL,
"RT properties struct");
/* Set up stuff that needs array allocation */
rt_init_grackle(&props->grackle_units, &props->grackle_chemistry_data,
&props->grackle_chemistry_rates,
props->hydrogen_mass_fraction, props->grackle_verbose,
props->case_B_recombination, us, cosmo);
props->energy_weighted_cross_sections = NULL;
props->number_weighted_cross_sections = NULL;
rt_cross_sections_init(props, phys_const, us);
/* The RT parameters, in particular the reduced speed of light, are
* not defined at compile time. So we need to write them down. */
restart_read_blocks(&rt_params, sizeof(struct rt_parameters), 1, stream, NULL,
"RT global parameters struct");
}
#endif /* SWIFT_RT_PROPERTIES_GEAR_H */