/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2022 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_GRACKLE_UTILS_H
#define SWIFT_RT_GRACKLE_UTILS_H
/* skip deprecation warnings. I cleaned old API calls. */
#define OMIT_LEGACY_INTERNAL_GRACKLE_FUNC
/* need hydro gamma */
#include "hydro.h"
#include
/* need to rework (and check) code if changed */
#define FIELD_SIZE 1
/**
* @file src/rt/GEAR/rt_grackle_utils.h
* @brief Utility and helper functions related to using grackle.
*/
/**
* @brief Update grackle units during run
*
* @param grackle_units grackle units struct
* @param cosmo cosmology struct
*
* NOTE: In the current implementation, this function does nothing.
* However, there might be use-cases in the future (e.g. switching
* UV background on or off depending on redshift) that might be
* needed in the future, which can be implemented into this function.
*/
__attribute__((always_inline)) INLINE void update_grackle_units_cosmo(
code_units *grackle_units, const struct unit_system *us,
const struct cosmology *restrict cosmo) {}
/**
* @brief initialize grackle during rt_props_init
*
* @param grackle_units grackle units struct to fill up correctly.
* @param grackle_chemistry_dat grackle chemistry data struct to fill up
*correctly.
* @param hydrogen_mass_fraction global hydrogen mass fraction.
* @param grackle_verb run grackle in verbose mode?
* @param case_B_recombination use grackle with case B recombination?
* @param us #unit_system struct
**/
__attribute__((always_inline)) INLINE static void rt_init_grackle(
code_units *grackle_units, chemistry_data *grackle_chemistry_data,
chemistry_data_storage *grackle_chemistry_rates,
float hydrogen_mass_fraction, const int grackle_verb,
const int case_B_recombination, const struct unit_system *us,
const struct cosmology *restrict cosmo) {
grackle_verbose = grackle_verb;
/* Initialize units */
/* ---------------- */
/* we assume all quantities to be physical, not comoving */
grackle_units->a_units = 1.0;
grackle_units->a_value = 1.0;
grackle_units->comoving_coordinates = 0;
grackle_units->density_units =
units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
grackle_units->length_units =
units_cgs_conversion_factor(us, UNIT_CONV_LENGTH);
grackle_units->time_units = units_cgs_conversion_factor(us, UNIT_CONV_TIME);
/* Set velocity units */
set_velocity_units(grackle_units);
/* Chemistry Parameters */
/* -------------------- */
/* More details on
* https://grackle.readthedocs.io/en/grackle-3.2.0/Integration.html#chemistry-data
*/
if (local_initialize_chemistry_parameters(grackle_chemistry_data) == 0) {
error("Error in set_default_chemistry_parameters.");
}
/* chemistry on */
grackle_chemistry_data->use_grackle = 1;
/* cooling on */
/* NOTE: without cooling on, it also won't heat... */
grackle_chemistry_data->with_radiative_cooling = 1;
/* 6 species atomic H and He */
grackle_chemistry_data->primordial_chemistry = 1;
/* No dust processes */
grackle_chemistry_data->dust_chemistry = 0;
/* No H2 formation on dust */
grackle_chemistry_data->h2_on_dust = 0;
/* metal cooling (uses Cloudy) off (for now) */
grackle_chemistry_data->metal_cooling = 0;
/* no cooling below CMB temperature */
grackle_chemistry_data->cmb_temperature_floor = 1;
/* UV background off */
grackle_chemistry_data->UVbackground = 0;
/* data file - currently not used */
grackle_chemistry_data->grackle_data_file = "";
/* adiabatic index */
grackle_chemistry_data->Gamma = hydro_gamma;
/* we'll provide grackle with ionization and heating rates from RT */
grackle_chemistry_data->use_radiative_transfer = 1;
/* fraction by mass of Hydrogen in the metal-free portion of the gas */
grackle_chemistry_data->HydrogenFractionByMass = hydrogen_mass_fraction;
/* Use case B recombination? (On-the-spot approximation) */
grackle_chemistry_data->CaseBRecombination = case_B_recombination;
if (local_initialize_chemistry_data(grackle_chemistry_data,
grackle_chemistry_rates,
grackle_units) == 0) {
error("Error in initialize_chemistry_data");
}
}
/**
* @brief fill out a grackle field struct with the relevant (gas) data from a
*particle
*
* @param grackle_fields (return) grackle field to copy into
* @param density array of particle density
* @param internal_energy array of particle internal_energy
* @param species_densities array of species densities of particle (HI, HII,
*HeI, HeII, HeIII, e-)
* @param iact_rates array of interaction rates (heating, 3 ioniziation, H2
*dissociation)
*
**/
__attribute__((always_inline)) INLINE static void
rt_get_grackle_particle_fields(grackle_field_data *grackle_fields,
gr_float density, gr_float internal_energy,
gr_float species_densities[6],
gr_float iact_rates[5]) {
int *dimension = malloc(3 * sizeof(int));
int *start = malloc(3 * sizeof(int));
int *end = malloc(3 * sizeof(int));
dimension[0] = FIELD_SIZE;
dimension[1] = 0;
dimension[2] = 0;
start[0] = 0;
start[1] = 0;
start[2] = 0;
end[0] = FIELD_SIZE - 1;
end[1] = 0;
end[2] = 0;
grackle_fields->grid_dx = 0.;
grackle_fields->grid_rank = 3;
grackle_fields->grid_dimension = dimension;
grackle_fields->grid_start = start;
grackle_fields->grid_end = end;
/* Set initial quantities */
grackle_fields->density = malloc(FIELD_SIZE * sizeof(gr_float));
grackle_fields->internal_energy = malloc(FIELD_SIZE * sizeof(gr_float));
grackle_fields->x_velocity = NULL;
grackle_fields->y_velocity = NULL;
grackle_fields->z_velocity = NULL;
/* for primordial_chemistry >= 1 */
grackle_fields->HI_density = malloc(FIELD_SIZE * sizeof(gr_float));
grackle_fields->HII_density = malloc(FIELD_SIZE * sizeof(gr_float));
grackle_fields->HeI_density = malloc(FIELD_SIZE * sizeof(gr_float));
grackle_fields->HeII_density = malloc(FIELD_SIZE * sizeof(gr_float));
grackle_fields->HeIII_density = malloc(FIELD_SIZE * sizeof(gr_float));
grackle_fields->e_density = malloc(FIELD_SIZE * sizeof(gr_float));
/* for primordial_chemistry >= 2 */
grackle_fields->HM_density = NULL;
grackle_fields->H2I_density = NULL;
grackle_fields->H2II_density = NULL;
/* for primordial_chemistry >= 3 */
grackle_fields->DI_density = NULL;
grackle_fields->DII_density = NULL;
grackle_fields->HDI_density = NULL;
/* for metal_cooling = 1 */
grackle_fields->metal_density = NULL;
/* for use_dust_density_field = 1 */
grackle_fields->dust_density = NULL;
/* volumetric heating rate (provide in units [erg s^-1 cm^-3]) */
grackle_fields->volumetric_heating_rate = NULL;
/* specific heating rate (provide in units [egs s^-1 g^-1] */
grackle_fields->specific_heating_rate = NULL;
/* radiative transfer ionization / dissociation rate fields (provide in units
* [1/s]) */
grackle_fields->RT_HI_ionization_rate = malloc(FIELD_SIZE * sizeof(gr_float));
grackle_fields->RT_HeI_ionization_rate =
malloc(FIELD_SIZE * sizeof(gr_float));
grackle_fields->RT_HeII_ionization_rate =
malloc(FIELD_SIZE * sizeof(gr_float));
grackle_fields->RT_H2_dissociation_rate =
malloc(FIELD_SIZE * sizeof(gr_float));
/* radiative transfer heating rate field
* (provide in units [erg s^-1 cm^-3] / nHI in cgs) */
grackle_fields->RT_heating_rate = malloc(FIELD_SIZE * sizeof(gr_float));
grackle_fields->H2_self_shielding_length = NULL;
grackle_fields->H2_custom_shielding_factor = NULL;
grackle_fields->isrf_habing = NULL;
for (int i = 0; i < FIELD_SIZE; i++) {
grackle_fields->density[i] = density;
grackle_fields->internal_energy[i] = internal_energy;
grackle_fields->HI_density[i] = species_densities[0];
grackle_fields->HII_density[i] = species_densities[1];
grackle_fields->HeI_density[i] = species_densities[2];
grackle_fields->HeII_density[i] = species_densities[3];
grackle_fields->HeIII_density[i] = species_densities[4];
/* e_density = electron density*mh/me = n_e * m_h */
grackle_fields->e_density[i] = species_densities[5];
/* grackle_fields->HM_density[i] = species_densities[6]; */
/* grackle_fields->H2I_density[i] = species_densities[7]; */
/* grackle_fields->H2II_density[i] = species_densities[8]; */
/* grackle_fields->DI_density[i] = species_densities[9]; */
/* grackle_fields->DII_density[i] = species_densities[10]; */
/* grackle_fields->HDI_density[i] = species_densities[11]; */
/* grackle_fields->metal_density[i] = 0.0; */
/* solar metallicity */
/* grackle_chemistry_data.SolarMetalFractionByMass *
* grackle_fields->density[i]; */
/* grackle_fields->x_velocity[i] = 0.0; */
/* grackle_fields->y_velocity[i] = 0.0; */
/* grackle_fields->z_velocity[i] = 0.0; */
/* grackle_fields->volumetric_heating_rate[i] = 0.0; */
/* grackle_fields->specific_heating_rate[i] = 0.0; */
grackle_fields->RT_heating_rate[i] = iact_rates[0];
grackle_fields->RT_HI_ionization_rate[i] = iact_rates[1];
grackle_fields->RT_HeI_ionization_rate[i] = iact_rates[2];
grackle_fields->RT_HeII_ionization_rate[i] = iact_rates[3];
grackle_fields->RT_H2_dissociation_rate[i] = iact_rates[4];
}
}
/**
* @brief free arrays allocated in grackle_fields.
*
* @param grackle_fields grackle fields to clean up
*
**/
__attribute__((always_inline)) INLINE static void rt_clean_grackle_fields(
grackle_field_data *grackle_fields) {
free(grackle_fields->grid_dimension);
free(grackle_fields->grid_start);
free(grackle_fields->grid_end);
/* initial quantities */
free(grackle_fields->density);
free(grackle_fields->internal_energy);
/* free(grackle_fields->x_velocity); */
/* free(grackle_fields->y_velocity); */
/* free(grackle_fields->z_velocity); */
/* for primordial_chemistry >= 1 */
free(grackle_fields->HI_density);
free(grackle_fields->HII_density);
free(grackle_fields->HeI_density);
free(grackle_fields->HeII_density);
free(grackle_fields->HeIII_density);
free(grackle_fields->e_density);
/* for primordial_chemistry >= 2 */
/* free(grackle_fields->HM_density); */
/* free(grackle_fields->H2I_density); */
/* free(grackle_fields->H2II_density); */
/* for primordial_chemistry >= 3 */
/* free(grackle_fields->DI_density); */
/* free(grackle_fields->DII_density); */
/* free(grackle_fields->HDI_density); */
/* for metal_cooling = 1 */
/* free(grackle_fields->metal_density); */
/* for use_dust_density_field = 1 */
/* free(grackle_fields->dust_density); */
/* free(grackle_fields->volumetric_heating_rate); */
/* free(grackle_fields->specific_heating_rate); */
free(grackle_fields->RT_HI_ionization_rate);
free(grackle_fields->RT_HeI_ionization_rate);
free(grackle_fields->RT_HeII_ionization_rate);
free(grackle_fields->RT_H2_dissociation_rate);
free(grackle_fields->RT_heating_rate);
/* free(grackle_fields->H2_self_shielding_length); */
/* free(grackle_fields->H2_custom_shielding_factor); */
/* free(grackle_fields->isrf_habing); */
}
/**
* @brief Write out all available grackle field data for a given index
* and setup to a file.
* This function is intended for debugging.
*
* @param fp FILE pointer to write into.
* @param grackle_fields grackle field data
* @param grackle_chemistry_data grackle chemistry data.
* @param grackle_units units used by grackle
* @param field_index grackle field index to print out.
**/
__attribute__((always_inline)) INLINE static void
rt_write_grackle_setup_and_field(FILE *fp, grackle_field_data grackle_fields,
chemistry_data *grackle_chemistry_data,
code_units *grackle_units, int field_index) {
fprintf(fp, "Grackle chemistry parameters:\n");
fprintf(fp, "use_grackle = %d\n",
grackle_chemistry_data->use_grackle);
fprintf(fp, "with_radiative_cooling = %d\n",
grackle_chemistry_data->with_radiative_cooling);
fprintf(fp, "primordial_chemistry = %d\n",
grackle_chemistry_data->primordial_chemistry);
fprintf(fp, "dust_chemistry = %d\n",
grackle_chemistry_data->dust_chemistry);
fprintf(fp, "metal_cooling = %d\n",
grackle_chemistry_data->metal_cooling);
fprintf(fp, "UVbackground = %d\n",
grackle_chemistry_data->UVbackground);
fprintf(fp, "grackle_data_file = %s\n",
grackle_chemistry_data->grackle_data_file);
fprintf(fp, "cmb_temperature_floor = %d\n",
grackle_chemistry_data->cmb_temperature_floor);
fprintf(fp, "Gamma = %g\n",
grackle_chemistry_data->Gamma);
fprintf(fp, "h2_on_dust = %d\n",
grackle_chemistry_data->h2_on_dust);
fprintf(fp, "use_dust_density_field = %d\n",
grackle_chemistry_data->use_dust_density_field);
fprintf(fp, "dust_recombination_cooling = %d\n",
grackle_chemistry_data->dust_recombination_cooling);
fprintf(fp, "photoelectric_heating = %d\n",
grackle_chemistry_data->photoelectric_heating);
fprintf(fp, "photoelectric_heating_rate = %g\n",
grackle_chemistry_data->photoelectric_heating_rate);
fprintf(fp, "use_isrf_field = %d\n",
grackle_chemistry_data->use_isrf_field);
fprintf(fp, "interstellar_radiation_field = %g\n",
grackle_chemistry_data->interstellar_radiation_field);
fprintf(fp, "use_volumetric_heating_rate = %d\n",
grackle_chemistry_data->use_volumetric_heating_rate);
fprintf(fp, "use_specific_heating_rate = %d\n",
grackle_chemistry_data->use_specific_heating_rate);
fprintf(fp, "three_body_rate = %d\n",
grackle_chemistry_data->three_body_rate);
fprintf(fp, "cie_cooling = %d\n",
grackle_chemistry_data->cie_cooling);
fprintf(fp, "h2_optical_depth_approximation = %d\n",
grackle_chemistry_data->h2_optical_depth_approximation);
fprintf(fp, "ih2co = %d\n",
grackle_chemistry_data->ih2co);
fprintf(fp, "ipiht = %d\n",
grackle_chemistry_data->ipiht);
fprintf(fp, "HydrogenFractionByMass = %g\n",
grackle_chemistry_data->HydrogenFractionByMass);
fprintf(fp, "DeuteriumToHydrogenRatio = %g\n",
grackle_chemistry_data->DeuteriumToHydrogenRatio);
fprintf(fp, "SolarMetalFractionByMass = %g\n",
grackle_chemistry_data->SolarMetalFractionByMass);
fprintf(fp, "local_dust_to_gas_ratio = %g\n",
grackle_chemistry_data->local_dust_to_gas_ratio);
fprintf(fp, "NumberOfTemperatureBins = %d\n",
grackle_chemistry_data->NumberOfTemperatureBins);
fprintf(fp, "CaseBRecombination = %d\n",
grackle_chemistry_data->CaseBRecombination);
fprintf(fp, "TemperatureStart = %g\n",
grackle_chemistry_data->TemperatureStart);
fprintf(fp, "TemperatureEnd = %g\n",
grackle_chemistry_data->TemperatureEnd);
fprintf(fp, "NumberOfDustTemperatureBins = %d\n",
grackle_chemistry_data->NumberOfDustTemperatureBins);
fprintf(fp, "DustTemperatureStart = %g\n",
grackle_chemistry_data->DustTemperatureStart);
fprintf(fp, "DustTemperatureEnd = %g\n",
grackle_chemistry_data->DustTemperatureEnd);
fprintf(fp, "Compton_xray_heating = %d\n",
grackle_chemistry_data->Compton_xray_heating);
fprintf(fp, "LWbackground_sawtooth_suppression = %d\n",
grackle_chemistry_data->LWbackground_sawtooth_suppression);
fprintf(fp, "LWbackground_intensity = %g\n",
grackle_chemistry_data->LWbackground_intensity);
fprintf(fp, "UVbackground_redshift_on = %g\n",
grackle_chemistry_data->UVbackground_redshift_on);
fprintf(fp, "UVbackground_redshift_off = %g\n",
grackle_chemistry_data->UVbackground_redshift_off);
fprintf(fp, "UVbackground_redshift_fullon = %g\n",
grackle_chemistry_data->UVbackground_redshift_fullon);
fprintf(fp, "UVbackground_redshift_drop = %g\n",
grackle_chemistry_data->UVbackground_redshift_drop);
fprintf(fp, "cloudy_electron_fraction_factor = %g\n",
grackle_chemistry_data->cloudy_electron_fraction_factor);
fprintf(fp, "use_radiative_transfer = %d\n",
grackle_chemistry_data->use_radiative_transfer);
fprintf(fp, "radiative_transfer_coupled_rate_solver = %d\n",
grackle_chemistry_data->radiative_transfer_coupled_rate_solver);
fprintf(fp, "radiative_transfer_intermediate_step = %d\n",
grackle_chemistry_data->radiative_transfer_intermediate_step);
fprintf(fp, "radiative_transfer_hydrogen_only = %d\n",
grackle_chemistry_data->radiative_transfer_hydrogen_only);
fprintf(fp, "self_shielding_method = %d\n",
grackle_chemistry_data->self_shielding_method);
fprintf(fp, "H2_custom_shielding = %d\n",
grackle_chemistry_data->H2_custom_shielding);
fprintf(fp, "H2_self_shielding = %d\n",
grackle_chemistry_data->H2_self_shielding);
fprintf(fp, "\nUnits:\n");
fprintf(fp, "a_units = %g\n", grackle_units->a_units);
fprintf(fp, "a_value = %g\n", grackle_units->a_value);
fprintf(fp, "comoving_coordinates = %d\n",
grackle_units->comoving_coordinates);
fprintf(fp, "density_units = %g\n", grackle_units->density_units);
fprintf(fp, "length_units = %g\n", grackle_units->length_units);
fprintf(fp, "time_units = %g\n", grackle_units->time_units);
fprintf(fp, "velocity_units = %g\n", grackle_units->velocity_units);
#define rt_print_grackle_field(v) \
if (grackle_fields.v != NULL) \
fprintf(fp, "grackle_fields." #v " = %g\n", grackle_fields.v[field_index])
fprintf(fp, "\nGrackle field data:\n");
rt_print_grackle_field(density);
rt_print_grackle_field(internal_energy);
rt_print_grackle_field(HI_density);
rt_print_grackle_field(HII_density);
rt_print_grackle_field(HeI_density);
rt_print_grackle_field(HeII_density);
rt_print_grackle_field(HeIII_density);
rt_print_grackle_field(e_density);
rt_print_grackle_field(HM_density);
rt_print_grackle_field(H2I_density);
rt_print_grackle_field(H2II_density);
rt_print_grackle_field(DI_density);
rt_print_grackle_field(DII_density);
rt_print_grackle_field(HDI_density);
rt_print_grackle_field(metal_density);
rt_print_grackle_field(x_velocity);
rt_print_grackle_field(y_velocity);
rt_print_grackle_field(z_velocity);
rt_print_grackle_field(volumetric_heating_rate);
rt_print_grackle_field(specific_heating_rate);
rt_print_grackle_field(RT_HI_ionization_rate);
rt_print_grackle_field(RT_HeI_ionization_rate);
rt_print_grackle_field(RT_HeII_ionization_rate);
rt_print_grackle_field(RT_H2_dissociation_rate);
rt_print_grackle_field(RT_heating_rate);
#undef rt_print_grackle_field
}
#endif /* SWIFT_RT_GRACKLE_UTILS_H */