Commit e7a5317e authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Use the P20 cooling tables for the Quick Lyman-alpha model

parent 9cdbe7b0
......@@ -1489,7 +1489,7 @@ AC_SUBST([NUMA_INCS])
# As an example for this, see the call to AC_ARG_WITH for cooling.
AC_ARG_WITH([subgrid],
[AS_HELP_STRING([--with-subgrid=<subgrid>],
[Master switch for subgrid methods. Inexperienced user should start here. Options are: @<:@none, GEAR, QLA, EAGLE default: none@:>@]
[Master switch for subgrid methods. Inexperienced user should start here. Options are: @<:@none, GEAR, QLA, EAGLE, EAGLE-XL default: none@:>@]
)],
[with_subgrid="$withval"],
[with_subgrid=none]
......
......@@ -86,12 +86,13 @@ LineOfSight:
# Quick Lyman-alpha cooling (EAGLE with fixed primoridal Z)
QLACooling:
dir_name: ./coolingtables/
H_reion_z: 7.5 # Planck 2018
dir_name: ./UV_dust1_CR1_G1_shield1.hdf5 # Location of the cooling tables
H_reion_z: 7.5 # Redshift of Hydrogen re-ionization (Planck 2018)
H_reion_eV_p_H: 2.0
He_reion_z_centre: 3.5
He_reion_z_sigma: 0.5
He_reion_eV_p_H: 2.0
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
rapid_cooling_threshold: 0.333333 # Switch to rapid cooling regime for dt / t_cool above this threshold.
# Quick Lyman-alpha star formation parameters
QLAStarFormation:
......
#!/bin/bash
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/EAGLE/coolingtables.tar.gz
tar -xf coolingtables.tar.gz
wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/CoolingTables/COLIBRE/UV_dust1_CR1_G1_shield1.hdf5
......@@ -389,14 +389,26 @@ EAGLECooling:
Ca_over_Si_in_solar: 1. # (Optional) Ratio of Ca/Si to use in units of solar. If set to 1, the code uses [Ca/Si] = 0, i.e. Ca/Si = 0.0941736.
S_over_Si_in_solar: 1. # (Optional) Ratio of S/Si to use in units of solar. If set to 1, the code uses [S/Si] = 0, i.e. S/Si = 0.6054160.
# Quick Lyman-alpha cooling (EAGLE with fixed primoridal Z)
# Quick Lyman-alpha cooling (EAGLE-XL with fixed primoridal Z)
QLACooling:
dir_name: ./coolingtables/ # Location of the Wiersma+08 cooling tables
H_reion_z: 7.5 # Redshift of Hydrogen re-ionization
H_reion_eV_p_H: 2.0 # Energy inject by Hydrogen re-ionization in electron-volt per Hydrogen atom
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
dir_name: ./UV_dust1_CR1_G1_shield1.hdf5 # Location of the Ploeckinger+20 cooling tables
H_reion_z: 7.5 # Redshift of Hydrogen re-ionization (Planck 2018)
H_reion_eV_p_H: 2.0 # Energy inject by Hydrogen re-ionization in electron-volt per Hydrogen atom
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
rapid_cooling_threshold: 0.333333 # Switch to rapid cooling regime for dt / t_cool above this threshold.
# COLIBRE cooling parameters (EAGLE-XL)
COLIBRECooling:
dir_name: ./UV_dust1_CR1_G1_shield1.hdf5 # Location of the cooling tables
H_reion_z: 7.5 # Redshift of Hydrogen re-ionization (Planck 2018)
H_reion_eV_p_H: 2.0 # Energy inject by Hydrogen re-ionization in electron-volt per Hydrogen atom
He_reion_z_centre: 3.5 # Redshift of the centre of the Helium re-ionization Gaussian
He_reion_z_sigma: 0.5 # Spread in redshift of the Helium re-ionization Gaussian
He_reion_eV_p_H: 2.0 # Energy inject by Helium re-ionization in electron-volt per Hydrogen atom
rapid_cooling_threshold: 0.333333 # Switch to rapid cooling regime for dt / t_cool above this threshold.
delta_logTEOS_subgrid_properties: 0.3 # delta log T above the EOS below which the subgrid properties use Teq assumption
# Cooling with Grackle 3.0
GrackleCooling:
......
This diff is collapsed.
......@@ -20,8 +20,9 @@
#define SWIFT_COOLING_QLA_H
/**
* @file src/cooling/EAGLE/cooling.h
* @brief EAGLE cooling function declarations
* @file src/cooling/QLA/cooling.h
* @brief Quick Lyman alpha cooling (COLIBRE with primoridal Z) function
* declarations
*/
/* Local includes. */
......@@ -32,10 +33,8 @@ struct xpart;
struct cosmology;
struct hydro_props;
struct entropy_floor_properties;
struct space;
struct unit_system;
struct phys_const;
struct swift_params;
struct space;
void cooling_update(const struct cosmology *cosmo,
struct cooling_function_data *cooling, struct space *s);
......@@ -46,34 +45,36 @@ void cooling_cool_part(const struct phys_const *phys_const,
const struct hydro_props *hydro_properties,
const struct entropy_floor_properties *floor_props,
const struct cooling_function_data *cooling,
struct part *restrict p, struct xpart *restrict xp,
const float dt, const float dt_therm, const double time);
struct part *p, struct xpart *xp, const float dt,
const float dt_therm, const double time);
float cooling_timestep(const struct cooling_function_data *restrict cooling,
const struct phys_const *restrict phys_const,
const struct cosmology *restrict cosmo,
const struct unit_system *restrict us,
float cooling_timestep(const struct cooling_function_data *cooling,
const struct phys_const *phys_const,
const struct cosmology *cosmo,
const struct unit_system *us,
const struct hydro_props *hydro_props,
const struct part *restrict p,
const struct xpart *restrict xp);
void cooling_first_init_part(
const struct phys_const *restrict phys_const,
const struct unit_system *restrict us,
const struct hydro_props *hydro_props,
const struct cosmology *restrict cosmo,
const struct cooling_function_data *restrict cooling,
const struct part *restrict p, struct xpart *restrict xp);
float cooling_get_temperature(
const struct phys_const *restrict phys_const,
const struct hydro_props *restrict hydro_props,
const struct unit_system *restrict us,
const struct cosmology *restrict cosmo,
const struct cooling_function_data *restrict cooling,
const struct part *restrict p, const struct xpart *restrict xp);
float cooling_get_radiated_energy(const struct xpart *restrict xp);
const struct part *p, const struct xpart *xp);
void cooling_first_init_part(const struct phys_const *phys_const,
const struct unit_system *us,
const struct hydro_props *hydro_props,
const struct cosmology *cosmo,
const struct cooling_function_data *cooling,
struct part *p, struct xpart *xp);
float cooling_get_temperature_from_gas(
const struct phys_const *phys_const, const struct cosmology *cosmo,
const struct cooling_function_data *cooling, const float rho_phys,
const float XH, const float logZZsol, const float u_phys);
float cooling_get_temperature(const struct phys_const *phys_const,
const struct hydro_props *hydro_props,
const struct unit_system *us,
const struct cosmology *cosmo,
const struct cooling_function_data *cooling,
const struct part *p, const struct xpart *xp);
float cooling_get_radiated_energy(const struct xpart *xp);
void cooling_split_part(struct part *p, struct xpart *xp, double n);
......
......@@ -41,7 +41,7 @@ __attribute__((always_inline)) INLINE static void cooling_write_flavour(
const struct cooling_function_data* cooling) {
io_write_attribute_s(h_grp, "Cooling Model",
"Quick Lyman-alpha (EAGLE with primordial Z only)");
"Quick Lyman-alpha (COLIBRE with primordial Z only)");
}
#endif
......@@ -70,6 +70,7 @@ __attribute__((always_inline)) INLINE static int cooling_write_particles(
list[0] = io_make_output_field_convert_part(
"Temperatures", FLOAT, 1, UNIT_CONV_TEMPERATURE, 0.f, parts, xparts,
convert_part_T, "Temperatures of the gas particles");
return 1;
}
......
This diff is collapsed.
......@@ -26,20 +26,50 @@
*/
struct cooling_tables {
/* array of heating rates due to metals */
float *metal_heating;
/* array of all mean particle masses mu (temperature) */
float *Tmu;
/* array of heating rates due to hydrogen and helium */
float *H_plus_He_heating;
/* array of all mean particle masses mu (internal energy) */
float *Umu;
/* array of electron abundances due to hydrogen and helium */
float *H_plus_He_electron_abundance;
/* array of all cooling processes (temperature) */
float *Tcooling;
/* array of temperatures */
float *temperature;
/* array of all cooling processes (internal energy) */
float *Ucooling;
/* array of electron abundances due to metals */
float *electron_abundance;
/* array of all heating processes (temperature) */
float *Theating;
/* array of all heating processes (internal energy) */
float *Uheating;
/* array of all electron abundances (temperature) */
float *Telectron_fraction;
/* array of all electron abundances (internal energy) */
float *Uelectron_fraction;
/* array to get T from U */
float *T_from_U;
/* array to get U from T */
float *U_from_T;
/* array of equilibrium temperatures */
float *logTeq;
/* array of mean particle masses at equilibrium temperatures */
float *meanpartmass_Teq;
/* array of pressures at equilibrium temperatures */
float *logPeq;
/* array of hydrogen fractions at equilibrium temperature */
float *logHfracs_Teq;
/* array of all hydrogen fractions */
float *logHfracs_all;
};
/**
......@@ -59,17 +89,33 @@ struct cooling_function_data {
/*! Temperature bins */
float *Temp;
/*! Helium fraction bins */
float *HeFrac;
/*! Metallicity bins */
float *Metallicity;
/*! Internal energy bins */
float *Therm;
/*! Mass fractions of elements for solar abundances (from the tables) */
float *SolarAbundances;
/*! Abundance ratios for each metallicity bin and for each included element */
float *LogAbundances;
float *Abundances;
float *Abundances_inv;
/*! Atomic masses for all included elements */
float *atomicmass;
float *atomicmass_inv;
/*! Mass fractions of all included elements */
float *LogMassFractions;
float *MassFractions;
/*! Index for solar metallicity in the metallicity dimension */
int indxZsol;
/*! Inverse of the solar mass fractions */
float *SolarAbundances_inv;
/*! Solar metallicity (metal mass fraction) */
float *Zsol;
/*! Inverse of solar metallicity (metal mass fraction) */
float *Zsol_inv;
/*! Filepath to the directory containing the HDF5 cooling tables */
char cooling_table_path[qla_table_path_name_length];
......@@ -80,7 +126,7 @@ struct cooling_function_data {
/*! H reionization energy in CGS units */
float H_reion_heat_cgs;
/*! Have we already done H reioisation? */
/*! Have we already done H reionization? */
int H_reion_done;
/*! Redshift of He reionization */
......@@ -100,26 +146,35 @@ struct cooling_function_data {
*/
double internal_energy_from_cgs;
/*! Pressure conversion from internal units to CGS (for quick access) */
double pressure_to_cgs;
/*! Number density conversion from internal units to CGS (for quick access) */
double number_density_to_cgs;
/*! Number density conversion from CGS to internal units (for quick access) */
double number_density_from_cgs;
/*! Inverse of proton mass in cgs (for quick access) */
double inv_proton_mass_cgs;
/*! Temperature of the CMB at present day (for quick access) */
/*! Logarithm base 10 of the Boltzmann constant in CGS (for quick access) */
double log10_kB_cgs;
/*! Temperatur of the CMB at present day (for quick access) */
double T_CMB_0;
/*! Compton rate in cgs units */
double compton_rate_cgs;
/*! Index of the current redshift along the redshift index of the tables */
int z_index;
/*! Minimal temperature allowed for the gas particles */
double Tmin;
/*! Distance between the current redshift and table[z_index] */
float dz;
/*! Minimal internal energy in cgs allowed for the gas particles */
double umin_cgs;
/*! Index of the previous tables along the redshift index of the tables */
int previous_z_index;
/*! Threshold to switch between rapid and slow cooling regimes. */
double rapid_cooling_threshold;
};
/**
......
This diff is collapsed.
......@@ -20,46 +20,99 @@
#define SWIFT_QLA_COOL_TABLES_H
/**
* @file src/cooling/QLA/cooling.h
* @brief Quick Lyman-alpha cooling function (EAGLE with primordial Z only)
* @file src/cooling/QLA/cooling_tables.h
* @brief QLA cooling function
*/
/* Config parameters. */
#include "../config.h"
#include "config.h"
#include "cooling_struct.h"
/*! Number of different bins along the temperature axis of the tables */
#define qla_cooling_N_temperature 86
/*! Number of different bins along the redhsift axis of the tables */
#define qla_cooling_N_redshifts 49
/*! Number of different bins along the redshift axis of the tables */
#define qla_cooling_N_redshifts 46
/*! Number of redshift bins loaded at any given point int time */
#define qla_cooling_N_loaded_redshifts 2
/*! Number of different bins along the density axis of the tables */
#define qla_cooling_N_density 71
/*! Number of different bins along the temperature axis of the tables */
#define qla_cooling_N_temperature 176
/*! Number of different bins along the metallicity axis of the tables */
#define qla_cooling_N_metallicity 11
/*! Number of different bins along the density axis of the tables */
#define qla_cooling_N_density 41
/*! Number of different bins along the internal energy axis of the tables */
#define qla_cooling_N_internalenergy 191
/*! Number of different bins along the metal axis of the tables */
#define qla_cooling_N_metal 9
/*! Number of different cooling channels in the tables */
#define qla_cooling_N_cooltypes 22
/*! Number of different bins along the metal axis of the tables */
#define qla_cooling_N_He_frac 7
/*! Number of different heating channels in the tables */
#define qla_cooling_N_heattypes 24
/*! Number of different bins along the abundances axis of the tables */
#define qla_cooling_N_abundances 11
/*! Number of different electron fractions (each element - other atoms
* + tot prim + tot metal + tot) in the tables */
#define qla_cooling_N_electrontypes 14
void get_cooling_redshifts(struct cooling_function_data *cooling);
/*! Number of different elements in the tables */
#define qla_cooling_N_elementtypes 12
void read_cooling_header(const char *fname,
struct cooling_function_data *cooling);
/**
* @brief Elements present in the tables
*/
enum qla_cooling_element {
element_H,
element_He,
element_C,
element_N,
element_O,
element_Ne,
element_Mg,
element_Si,
element_S,
element_Ca,
element_Fe,
element_OA
};
void allocate_cooling_tables(struct cooling_function_data *restrict cooling);
/**
* @brief Hydrogen species
*/
enum qla_hydrogen_species { neutral = 0, ionized = 1, molecular = 2 };
void get_redshift_invariant_table(
struct cooling_function_data *restrict cooling, const int photodis);
void get_cooling_table(struct cooling_function_data *restrict cooling,
const int low_z_index, const int high_z_index);
/**
* @brief Cooling channels beyond the metal lines
*/
enum qla_cooling_channels {
cooltype_H2 = element_OA + 1,
cooltype_molecules,
cooltype_HD,
cooltype_NetFFH,
cooltype_NetFFM,
cooltype_eeBrems,
cooltype_Compton,
cooltype_Dust
};
/**
* @brief Heating channels beyond the metal lines
*/
enum qla_heating_channels {
heattype_H2 = element_OA + 1,
heattype_COdiss,
heattype_CosmicRay,
heattype_UTA,
heattype_line,
heattype_Hlin,
heattype_ChaT,
heattype_HFF,
heattype_Compton,
heattype_Dust
};
/* Pre-declaration */
struct cooling_function_data;
void get_cooling_redshifts(struct cooling_function_data *cooling);
void read_cooling_header(struct cooling_function_data *cooling);
void read_cooling_tables(struct cooling_function_data *cooling);
#endif /* SWIFT_QLA_COOL_TABLES_H */
#endif
......@@ -21,7 +21,7 @@
/**
* @file src/cooling/QLA/interpolate.h
* @brief Interpolation functions for Quick Lyman-alpha cooling tables
* @brief Interpolation functions for Quick Lyman-alpha cooling tables (COLIBRE)
*/
/* Config parameters. */
......@@ -30,6 +30,7 @@
/* Local includes. */
#include "align.h"
#include "error.h"
#include "exp10.h"
#include "inline.h"
/**
......@@ -87,6 +88,28 @@ __attribute__((always_inline)) INLINE int row_major_index_4d(
return x * Ny * Nz * Nw + y * Nz * Nw + z * Nw + w;
}
/**
* @brief Returns the 1d index of element with 5d indices x,y,z,w
* from a flattened 5d array in row major order
*
* @param x, y, z, v, w Indices of element of interest
* @param Nx, Ny, Nz, Nv, Nw Sizes of array dimensions
*/
__attribute__((always_inline)) INLINE int row_major_index_5d(
const int x, const int y, const int z, const int w, const int v,
const int Nx, const int Ny, const int Nz, const int Nw, const int Nv) {
#ifdef SWIFT_DEBUG_CHECKS
assert(x < Nx);
assert(y < Ny);
assert(z < Nz);
assert(w < Nw);
assert(v < Nv);
#endif
return x * Ny * Nz * Nw * Nv + y * Nz * Nw * Nv + z * Nw * Nv + w * Nv + v;
}
/**
* @brief Finds the index of a value in a table and compute delta to nearest
* element.
......@@ -119,29 +142,32 @@ __attribute__((always_inline)) INLINE void get_index_1d(
swift_align_information(float, table, SWIFT_STRUCT_ALIGNMENT);
/* Distance between elements in the array */
const float delta = (size - 1) / (table[size - 1] - table[0]);
/* Do not use first or last entry, might be an extra bin with uneven spacing
*/
const float delta = (size - 3) / (table[size - 2] - table[1]);
if (x < table[0] + epsilon) {
/* We are below the first element */
*i = 0;
*dx = 0.f;
} else if (x < table[size - 1] - epsilon) {
/* Normal case */
*i = (x - table[0]) * delta;
/* Check for an extra entry at the beginning (e.g. metallicity) */
int istart = 0;
int iend = size - 1;
#ifdef SWIFT_DEBUG_CHECKS
if (*i > size || *i < 0) {
error(
"trying to get index for value outside table range. Table size: %d, "
"calculated index: %d, value: %.5e, table[0]: %.5e, grid size: %.5e",
size, *i, x, table[0], delta);
}
#endif
if (fabsf(table[1] - table[0]) > delta + epsilon) {
istart = 1;
}
if (fabsf(table[size - 1] - table[size - 2]) > delta + epsilon) {
iend = size - 2;
}
/*extra array at the beginning */
if (x < table[istart] + epsilon) {
/* We are before the first element */
*i = 0;
*dx = 0.f;
} else if (x < table[iend] - epsilon) {
*i = (x - table[1]) * delta + 1;
*dx = (x - table[*i]) * delta;
} else {
/* We are after the last element */
*i = size - 2;
*i = iend - 1;
*dx = 1.f;
}
......@@ -242,6 +268,67 @@ __attribute__((always_inline)) INLINE float interpolation_3d(
return result;
}
/**
* @brief Interpolate a flattened 3D table at a given position but avoid the
* z-dimension.
*
* This function uses linear interpolation along each axis.
* We look at the zi coordoniate but do not interpolate around it. We just
* interpolate the remaining 2 dimensions.
* The function also assumes that the table is aligned on
* SWIFT_STRUCT_ALIGNMENT.
*
* @param table The 3D table to interpolate.
* @param xi, yi, zi Indices of element of interest.
* @param Nx, Ny, Nz Sizes of array dimensions.
* @param dx, dy, dz Distance between the point and the index in units of
* the grid spacing.
*/
__attribute__((always_inline)) INLINE float interpolation_3d_no_z(
const float *table, const int xi, const int yi, const int zi,
const float dx, const float dy, const float dz, const int Nx, const int Ny,
const int Nz) {
#ifdef SWIFT_DEBUG_CHECKS
if (dx < -0.001f || dx > 1.001f) error("Invalid dx=%e", dx);
if (dy < -0.001f || dy > 1.001f) error("Invalid dy=%e", dy);
if (dz != 0.f) error("Attempting to interpolate along z!");
#endif
const float tx = 1.f - dx;
const float ty = 1.f - dy;
const float tz = 1.f;
/* Indicate that the whole array is aligned on page boundaries */
swift_align_information(float, table, SWIFT_STRUCT_ALIGNMENT);
/* Linear interpolation along each axis. We read the table 2^2=4 times */
/* Note that we intentionally kept the table access along the axis where */
/* we do not interpolate as comments in the code to allow readers to */
/* understand what is going on. */
float result = tx * ty * tz *
table[row_major_index_3d(xi + 0, yi + 0, zi + 0, Nx, Ny, Nz)];
/* result += tx * ty * dz *
table[row_major_index_3d(xi + 0, yi + 0, zi + 1, Nx, Ny, Nz)]; */
result += tx * dy * tz *
table[row_major_index_3d(xi + 0, yi + 1, zi + 0, Nx, Ny, Nz)];
result += dx * ty * tz *
table[row_major_index_3d(xi + 1, yi + 0, zi + 0, Nx, Ny, Nz)];
/* result += tx * dy * dz *
table[row_major_index_3d(xi + 0, yi + 1, zi + 1, Nx, Ny, Nz)]; */
/* result += dx * ty * dz *
table[row_major_index_3d(xi + 1, yi + 0, zi + 1, Nx, Ny, Nz)]; */
result += dx * dy * tz *
table[row_major_index_3d(xi + 1, yi + 1, zi + 0, Nx, Ny, Nz)];
/* result += dx * dy * dz *
table[row_major_index_3d(xi + 1, yi + 1, zi + 1, Nx, Ny, Nz)]; */
return result;
}
/**
* @brief Interpolate a flattened 3D table at a given position but avoid the
* x-dimension.
......@@ -392,6 +479,96 @@ __attribute__((always_inline)) INLINE float interpolation_4d(