Skip to content
Snippets Groups Projects
Commit fad25194 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added Quick Lyman-alpha cooling and chemistry model

parent 3896d2bc
No related branches found
No related tags found
1 merge request!1036Quick lyman alpha module
Showing
with 3656 additions and 1 deletion
...@@ -1443,7 +1443,7 @@ AC_SUBST([NUMA_INCS]) ...@@ -1443,7 +1443,7 @@ AC_SUBST([NUMA_INCS])
# As an example for this, see the call to AC_ARG_WITH for cooling. # As an example for this, see the call to AC_ARG_WITH for cooling.
AC_ARG_WITH([subgrid], AC_ARG_WITH([subgrid],
[AS_HELP_STRING([--with-subgrid=<subgrid>], [AS_HELP_STRING([--with-subgrid=<subgrid>],
[Master switch for subgrid methods. Inexperienced user should start from here @<:@none, GEAR, EAGLE default: none@:>@] [Master switch for subgrid methods. Inexperienced user should start here. Options are: @<:@none, GEAR, QLA, EAGLE default: none@:>@]
)], )],
[with_subgrid="$withval"], [with_subgrid="$withval"],
[with_subgrid=none] [with_subgrid=none]
...@@ -1478,6 +1478,18 @@ case "$with_subgrid" in ...@@ -1478,6 +1478,18 @@ case "$with_subgrid" in
with_subgrid_task_order=none with_subgrid_task_order=none
enable_fof=no enable_fof=no
;; ;;
QLA)
with_subgrid_cooling=QLA
with_subgrid_chemistry=QLA
with_subgrid_tracers=none
with_subgrid_entropy_floor=QLA
with_subgrid_stars=none
with_subgrid_star_formation=QLA
with_subgrid_feedback=none
with_subgrid_black_holes=none
with_subgrid_task_order=default
enable_fof=no
;;
EAGLE) EAGLE)
with_subgrid_cooling=EAGLE with_subgrid_cooling=EAGLE
with_subgrid_chemistry=EAGLE with_subgrid_chemistry=EAGLE
...@@ -1810,6 +1822,9 @@ case "$with_cooling" in ...@@ -1810,6 +1822,9 @@ case "$with_cooling" in
primordial_chemistry=${with_cooling:8} primordial_chemistry=${with_cooling:8}
AC_DEFINE_UNQUOTED([COOLING_GRACKLE_MODE], [$primordial_chemistry], [Grackle chemistry network]) AC_DEFINE_UNQUOTED([COOLING_GRACKLE_MODE], [$primordial_chemistry], [Grackle chemistry network])
;; ;;
QLA)
AC_DEFINE([COOLING_QLA], [1], [Cooling following the Quick-Lyman-alpha model])
;;
EAGLE) EAGLE)
AC_DEFINE([COOLING_EAGLE], [1], [Cooling following the EAGLE model]) AC_DEFINE([COOLING_EAGLE], [1], [Cooling following the EAGLE model])
;; ;;
...@@ -1845,6 +1860,9 @@ case "$with_chemistry" in ...@@ -1845,6 +1860,9 @@ case "$with_chemistry" in
number_element=${with_chemistry:5} number_element=${with_chemistry:5}
AC_DEFINE_UNQUOTED([GEAR_CHEMISTRY_ELEMENT_COUNT], [$number_element], [Number of element to follow]) AC_DEFINE_UNQUOTED([GEAR_CHEMISTRY_ELEMENT_COUNT], [$number_element], [Number of element to follow])
;; ;;
QLA)
AC_DEFINE([CHEMISTRY_QLA], [1], [Chemistry taken from the Quick-Lyman-alpha model])
;;
EAGLE) EAGLE)
AC_DEFINE([CHEMISTRY_EAGLE], [1], [Chemistry taken from the EAGLE model]) AC_DEFINE([CHEMISTRY_EAGLE], [1], [Chemistry taken from the EAGLE model])
;; ;;
...@@ -2083,6 +2101,9 @@ case "$with_entropy_floor" in ...@@ -2083,6 +2101,9 @@ case "$with_entropy_floor" in
none) none)
AC_DEFINE([ENTROPY_FLOOR_NONE], [1], [No entropy floor]) AC_DEFINE([ENTROPY_FLOOR_NONE], [1], [No entropy floor])
;; ;;
QLA)
AC_DEFINE([ENTROPY_FLOOR_QLA], [1], [Quick Lyman-alpha entropy floor])
;;
EAGLE) EAGLE)
AC_DEFINE([ENTROPY_FLOOR_EAGLE], [1], [EAGLE entropy floor]) AC_DEFINE([ENTROPY_FLOOR_EAGLE], [1], [EAGLE entropy floor])
;; ;;
...@@ -2140,6 +2161,9 @@ case "$with_star_formation" in ...@@ -2140,6 +2161,9 @@ case "$with_star_formation" in
none) none)
AC_DEFINE([STAR_FORMATION_NONE], [1], [No star formation]) AC_DEFINE([STAR_FORMATION_NONE], [1], [No star formation])
;; ;;
QLA)
AC_DEFINE([STAR_FORMATION_QLA], [1], [Quick Lyman-alpha star formation model)])
;;
EAGLE) EAGLE)
AC_DEFINE([STAR_FORMATION_EAGLE], [1], [EAGLE star formation model (Schaye and Dalla Vecchia (2008))]) AC_DEFINE([STAR_FORMATION_EAGLE], [1], [EAGLE star formation model (Schaye and Dalla Vecchia (2008))])
;; ;;
...@@ -2180,6 +2204,9 @@ AC_SUBST([GIT_CMD]) ...@@ -2180,6 +2204,9 @@ AC_SUBST([GIT_CMD])
DX_INIT_DOXYGEN(libswift,doc/Doxyfile,doc/) DX_INIT_DOXYGEN(libswift,doc/Doxyfile,doc/)
AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""]) AM_CONDITIONAL([HAVE_DOXYGEN], [test "$ac_cv_path_ac_pt_DX_DOXYGEN" != ""])
# Check if using QLA cooling
AM_CONDITIONAL([HAVEQLACOOLING], [test $with_cooling = "QLA"])
# Check if using EAGLE cooling # Check if using EAGLE cooling
AM_CONDITIONAL([HAVEEAGLECOOLING], [test $with_cooling = "EAGLE"]) AM_CONDITIONAL([HAVEEAGLECOOLING], [test $with_cooling = "EAGLE"])
......
...@@ -59,6 +59,12 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \ ...@@ -59,6 +59,12 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
black_holes.h black_holes_io.h black_holes_properties.h black_holes_struct.h \ black_holes.h black_holes_io.h black_holes_properties.h black_holes_struct.h \
feedback.h feedback_struct.h feedback_properties.h task_order.h feedback.h feedback_struct.h feedback_properties.h task_order.h
# source files for EAGLE cooling
QLA_COOLING_SOURCES =
if HAVEQLACOOLING
QLA_COOLING_SOURCES += cooling/QLA/cooling.c cooling/QLA/cooling_tables.c
endif
# source files for EAGLE cooling # source files for EAGLE cooling
EAGLE_COOLING_SOURCES = EAGLE_COOLING_SOURCES =
if HAVEEAGLECOOLING if HAVEEAGLECOOLING
...@@ -102,6 +108,7 @@ AM_SOURCES = space.c runner_main.c runner_doiact_hydro.c runner_doiact_limiter.c ...@@ -102,6 +108,7 @@ AM_SOURCES = space.c runner_main.c runner_doiact_hydro.c runner_doiact_limiter.c
chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \ chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \
outputlist.c velociraptor_dummy.c logger_io.c memuse.c mpiuse.c memuse_rnodes.c fof.c \ outputlist.c velociraptor_dummy.c logger_io.c memuse.c mpiuse.c memuse_rnodes.c fof.c \
hashmap.c pressure_floor.c \ hashmap.c pressure_floor.c \
$(QLA_COOLING_SOURCES) \
$(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES) \ $(EAGLE_COOLING_SOURCES) $(EAGLE_FEEDBACK_SOURCES) \
$(GRACKLE_COOLING_SOURCES) $(GEAR_FEEDBACK_SOURCES) $(GRACKLE_COOLING_SOURCES) $(GEAR_FEEDBACK_SOURCES)
...@@ -222,6 +229,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h ...@@ -222,6 +229,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
cooling/grackle/cooling_io.h \ cooling/grackle/cooling_io.h \
cooling/EAGLE/cooling.h cooling/EAGLE/cooling_struct.h cooling/EAGLE/cooling_tables.h \ cooling/EAGLE/cooling.h cooling/EAGLE/cooling_struct.h cooling/EAGLE/cooling_tables.h \
cooling/EAGLE/cooling_io.h cooling/EAGLE/interpolate.h cooling/EAGLE/cooling_rates.h \ cooling/EAGLE/cooling_io.h cooling/EAGLE/interpolate.h cooling/EAGLE/cooling_rates.h \
cooling/QLA/cooling.h cooling/QLA/cooling_struct.h cooling/QLA/cooling_tables.h \
cooling/QLA/cooling_io.h cooling/QLA/interpolate.h cooling/QLA/cooling_rates.h \
chemistry/none/chemistry.h \ chemistry/none/chemistry.h \
chemistry/none/chemistry_io.h \ chemistry/none/chemistry_io.h \
chemistry/none/chemistry_struct.h \ chemistry/none/chemistry_struct.h \
...@@ -234,6 +243,10 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h ...@@ -234,6 +243,10 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
chemistry/EAGLE/chemistry_io.h \ chemistry/EAGLE/chemistry_io.h \
chemistry/EAGLE/chemistry_struct.h\ chemistry/EAGLE/chemistry_struct.h\
chemistry/EAGLE/chemistry_iact.h \ chemistry/EAGLE/chemistry_iact.h \
chemistry/QLA/chemistry.h \
chemistry/QLA/chemistry_io.h \
chemistry/QLA/chemistry_struct.h\
chemistry/QLA/chemistry_iact.h \
entropy_floor/none/entropy_floor.h \ entropy_floor/none/entropy_floor.h \
entropy_floor/EAGLE/entropy_floor.h \ entropy_floor/EAGLE/entropy_floor.h \
tracers/none/tracers.h tracers/none/tracers_struct.h \ tracers/none/tracers.h tracers/none/tracers_struct.h \
......
...@@ -35,6 +35,9 @@ ...@@ -35,6 +35,9 @@
#elif defined(CHEMISTRY_GEAR) #elif defined(CHEMISTRY_GEAR)
#include "./chemistry/GEAR/chemistry.h" #include "./chemistry/GEAR/chemistry.h"
#include "./chemistry/GEAR/chemistry_iact.h" #include "./chemistry/GEAR/chemistry_iact.h"
#elif defined(CHEMISTRY_QLA)
#include "./chemistry/QLA/chemistry.h"
#include "./chemistry/QLA/chemistry_iact.h"
#elif defined(CHEMISTRY_EAGLE) #elif defined(CHEMISTRY_EAGLE)
#include "./chemistry/EAGLE/chemistry.h" #include "./chemistry/EAGLE/chemistry.h"
#include "./chemistry/EAGLE/chemistry_iact.h" #include "./chemistry/EAGLE/chemistry_iact.h"
......
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2020 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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_CHEMISTRY_QLA_H
#define SWIFT_CHEMISTRY_QLA_H
/**
* @file src/chemistry/none/chemistry.h
* @brief Empty infrastructure for the cases without chemistry function
*/
/* Some standard headers. */
#include <float.h>
#include <math.h>
/* Local includes. */
#include "chemistry_struct.h"
#include "cosmology.h"
#include "error.h"
#include "hydro.h"
#include "parser.h"
#include "part.h"
#include "physical_constants.h"
#include "units.h"
/**
* @brief Return a string containing the name of a given #chemistry_element.
*/
__attribute__((always_inline)) INLINE static const char*
chemistry_get_element_name(enum chemistry_element elem) {
static const char* chemistry_element_names[chemistry_element_count] = {};
return chemistry_element_names[elem];
}
/**
* @brief Initialises the chemistry properties.
*
* Nothing to do here.
*
* @param parameter_file The parsed parameter file.
* @param us The current internal system of units.
* @param phys_const The physical constants in internal units.
* @param data The global chemistry information (to be filled).
*/
static INLINE void chemistry_init_backend(struct swift_params* parameter_file,
const struct unit_system* us,
const struct phys_const* phys_const,
struct chemistry_global_data* data) {}
/**
* @brief Prints the properties of the chemistry model to stdout.
*
* Nothing to do here.
*
* @brief The #chemistry_global_data containing information about the current
* model.
*/
static INLINE void chemistry_print_backend(
const struct chemistry_global_data* data) {
message("Chemistry function is 'Quick Lyman-alpha'.");
}
/**
* @brief Finishes the density calculation.
*
* Nothing to do here.
*
* @param p The particle to act upon
* @param cd The global chemistry information.
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void chemistry_end_density(
struct part* restrict p, const struct chemistry_global_data* cd,
const struct cosmology* cosmo) {}
/**
* @brief Updates to the chemistry data after the hydro force loop.
*
* @param p The particle to act upon.
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void chemistry_end_force(
struct part* restrict p, const struct cosmology* cosmo) {}
/**
* @brief Computes the chemistry-related time-step constraint.
*
* @param phys_const The physical constants in internal units.
* @param cosmo The current cosmological model.
* @param us The internal system of units.
* @param hydro_props The properties of the hydro scheme.
* @param cd The global properties of the chemistry scheme.
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float chemistry_timestep(
const struct phys_const* restrict phys_const,
const struct cosmology* restrict cosmo,
const struct unit_system* restrict us,
const struct hydro_props* hydro_props,
const struct chemistry_global_data* cd, const struct part* restrict p) {
return FLT_MAX;
}
/**
* @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
*
* Nothing to do here.
*
* @param p The particle to act upon
* @param xp The extended particle data to act upon
* @param cd #chemistry_global_data containing chemistry informations.
* @param cosmo The current cosmological model.
*/
__attribute__((always_inline)) INLINE static void
chemistry_part_has_no_neighbours(struct part* restrict p,
struct xpart* restrict xp,
const struct chemistry_global_data* cd,
const struct cosmology* cosmo) {}
/**
* @brief Sets the chemistry properties of the (x-)particles to a valid start
* state.
*
* Nothing to do here.
*
* @param phys_const The physical constant in internal units.
* @param us The unit system.
* @param cosmo The current cosmological model.
* @param data The global chemistry information used for this run.
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
*/
__attribute__((always_inline)) INLINE static void chemistry_first_init_part(
const struct phys_const* restrict phys_const,
const struct unit_system* restrict us,
const struct cosmology* restrict cosmo,
const struct chemistry_global_data* data, const struct part* restrict p,
struct xpart* restrict xp) {}
/**
* @brief Sets the chemistry properties of the (x-)particles to a valid start
* state.
*
* Nothing to do here.
*
* @param p Pointer to the particle data.
* @param data The global chemistry information.
*/
__attribute__((always_inline)) INLINE static void chemistry_init_part(
struct part* restrict p, const struct chemistry_global_data* data) {}
/**
* @brief Sets the chemistry properties of the sparticles to a valid start
* state.
*
* Nothing to do here.
*
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
* @param cosmo The current cosmological model.
* @param data The global chemistry information.
* @param sp Pointer to the sparticle data.
* @param xp Pointer to the extended particle data.
*/
__attribute__((always_inline)) INLINE static void chemistry_first_init_spart(
const struct chemistry_global_data* data, struct spart* restrict sp) {}
/**
* @brief Initialise the chemistry properties of a black hole with
* the chemistry properties of the gas it is born from.
*
* Nothing to do here.
*
* @param bp_data The black hole data to initialise.
* @param p_data The gas data to use.
* @param gas_mass The mass of the gas particle.
*/
__attribute__((always_inline)) INLINE static void chemistry_bpart_from_part(
struct chemistry_bpart_data* bp_data,
const struct chemistry_part_data* p_data, const double gas_mass) {}
/**
* @brief Add the chemistry data of a gas particle to a black hole.
*
* Nothing to do here.
*
* @param bp_data The black hole data to add to.
* @param p_data The gas data to use.
* @param gas_mass The mass of the gas particle.
*/
__attribute__((always_inline)) INLINE static void chemistry_add_part_to_bpart(
struct chemistry_bpart_data* bp_data,
const struct chemistry_part_data* p_data, const double gas_mass) {}
/**
* @brief Add the chemistry data of a black hole to another one.
*
* Nothing to do here.
*
* @param bp_data The black hole data to add to.
* @param swallowed_data The black hole data to use.
*/
__attribute__((always_inline)) INLINE static void chemistry_add_bpart_to_bpart(
struct chemistry_bpart_data* bp_data,
const struct chemistry_bpart_data* swallowed_data) {}
/**
* @brief Split the metal content of a particle into n pieces
*
* Nothing to do here.
*
* @param p The #part.
* @param n The number of pieces to split into.
*/
__attribute__((always_inline)) INLINE static void chemistry_split_part(
struct part* p, const double n) {}
/**
* @brief Returns the total metallicity (metal mass fraction) of the
* star particle to be used in feedback/enrichment related routines.
*
* No metallicity treatment here -> return 0.
*
* @param sp Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float
chemistry_get_total_metal_mass_fraction_for_feedback(const struct spart* sp) {
return 0.f;
}
/**
* @brief Returns the abundance array (metal mass fractions) of the
* star particle to be used in feedback/enrichment related routines.
*
* No metallicity treatment here -> return NULL array.
*
* @param sp Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float const*
chemistry_get_metal_mass_fraction_for_feedback(const struct spart* sp) {
return NULL;
}
/**
* @brief Returns the total metallicity (metal mass fraction) of the
* gas particle to be used in cooling related routines.
*
* No metallicity treatment here -> return 0.
*
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float
chemistry_get_total_metal_mass_fraction_for_cooling(const struct part* p) {
return 0.f;
}
/**
* @brief Returns the abundance array (metal mass fractions) of the
* gas particle to be used in cooling related routines.
*
* No metallicity treatment here -> return NULL array.
*
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float const*
chemistry_get_metal_mass_fraction_for_cooling(const struct part* p) {
return NULL;
}
/**
* @brief Returns the total metallicity (metal mass fraction) of the
* gas particle to be used in star formation related routines.
*
* No metallicity treatment here -> return 0.
*
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float
chemistry_get_total_metal_mass_fraction_for_star_formation(
const struct part* p) {
return 0.f;
}
/**
* @brief Returns the abundance array (metal mass fractions) of the
* gas particle to be used in star formation related routines.
*
* No metallicity treatment here -> return NULL array.
*
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float const*
chemistry_get_metal_mass_fraction_for_star_formation(const struct part* p) {
return NULL;
}
#endif /* SWIFT_CHEMISTRY_QLA_H */
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2020 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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_QLA_CHEMISTRY_IACT_H
#define SWIFT_QLA_CHEMISTRY_IACT_H
/**
* @file none/chemistry_iact.h
* @brief Density computation
*/
/**
* @brief do chemistry computation after the runner_iact_density (symmetric
* version)
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi First particle.
* @param pj Second particle.
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_chemistry(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
struct part *restrict pj, float a, float H) {}
/**
* @brief do chemistry computation after the runner_iact_density (non symmetric
* version)
*
* @param r2 Comoving square distance between the two particles.
* @param dx Comoving vector separating both particles (pi - pj).
* @param hi Comoving smoothing-length of particle i.
* @param hj Comoving smoothing-length of particle j.
* @param pi First particle.
* @param pj Second particle (not updated).
* @param a Current scale factor.
* @param H Current Hubble parameter.
*/
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_chemistry(
float r2, const float *dx, float hi, float hj, struct part *restrict pi,
const struct part *restrict pj, float a, float H) {}
#endif /* SWIFT_QLA_CHEMISTRY_IACT_H */
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2020 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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_CHEMISTRY_IO_QLA_H
#define SWIFT_CHEMISTRY_IO_QLA_H
#include "io_properties.h"
/**
* @brief Specifies which particle fields to read from a dataset
*
* @param parts The particle array.
* @param list The list of i/o properties to read.
*
* @return Returns the number of fields to write.
*/
INLINE static int chemistry_read_particles(struct part* parts,
struct io_props* list) {
/* update list according to hydro_io */
/* Return the number of fields to read */
return 0;
}
/**
* @brief Specifies which particle fields to write to a dataset
*
* @param parts The particle array.
* @param list The list of i/o properties to write.
*
* @return Returns the number of fields to write.
*/
INLINE static int chemistry_write_particles(const struct part* parts,
struct io_props* list) {
/* update list according to hydro_io */
/* Return the number of fields to write */
return 0;
}
/**
* @brief Specifies which sparticle fields to write to a dataset
*
* @param sparts The sparticle array.
* @param list The list of i/o properties to write.
*
* @return Returns the number of fields to write.
*/
INLINE static int chemistry_write_sparticles(const struct spart* sparts,
struct io_props* list) {
/* update list according to hydro_io */
/* Return the number of fields to write */
return 0;
}
/**
* @brief Specifies which bparticle fields to write to a dataset
*
* @param bparts The bparticle array.
* @param list The list of i/o properties to write.
*
* @return Returns the number of fields to write.
*/
INLINE static int chemistry_write_bparticles(const struct bpart* bparts,
struct io_props* list) {
/* update list according to hydro_io */
/* Return the number of fields to write */
return 0;
}
#ifdef HAVE_HDF5
/**
* @brief Writes the current model of chemistry to the file
* @param h_grp The HDF5 group in which to write
*/
INLINE static void chemistry_write_flavour(hid_t h_grp) {
io_write_attribute_s(h_grp, "Chemistry Model", "Quick Lyman-alpha");
io_write_attribute_d(h_grp, "Chemistry element count", 0);
}
#endif
#endif /* SWIFT_CHEMISTRY_IO_QLA_H */
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2020 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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_CHEMISTRY_STRUCT_QLA_H
#define SWIFT_CHEMISTRY_STRUCT_QLA_H
/**
* @file src/chemistry/none/chemistry_struct.h
* @brief Empty infrastructure for the cases without chemistry function
*/
/**
* @brief The individual elements traced in the model.
*/
enum chemistry_element { chemistry_element_count = 0 };
/**
* @brief Global chemical abundance information.
*
* Nothing here.
*/
struct chemistry_global_data {};
/**
* @brief Chemistry properties carried by the #part.
*
* Nothing here.
*/
struct chemistry_part_data {};
/**
* @brief Chemistry properties carried by the #bpart.
*
* Nothing here.
*/
struct chemistry_bpart_data {};
#endif /* SWIFT_CHEMISTRY_STRUCT_QLA_H */
...@@ -27,6 +27,8 @@ ...@@ -27,6 +27,8 @@
#include "./chemistry/none/chemistry_io.h" #include "./chemistry/none/chemistry_io.h"
#elif defined(CHEMISTRY_GEAR) #elif defined(CHEMISTRY_GEAR)
#include "./chemistry/GEAR/chemistry_io.h" #include "./chemistry/GEAR/chemistry_io.h"
#elif defined(CHEMISTRY_QLA)
#include "./chemistry/QLA/chemistry_io.h"
#elif defined(CHEMISTRY_EAGLE) #elif defined(CHEMISTRY_EAGLE)
#include "./chemistry/EAGLE/chemistry_io.h" #include "./chemistry/EAGLE/chemistry_io.h"
#else #else
......
...@@ -32,6 +32,8 @@ ...@@ -32,6 +32,8 @@
#include "./chemistry/none/chemistry_struct.h" #include "./chemistry/none/chemistry_struct.h"
#elif defined(CHEMISTRY_GEAR) #elif defined(CHEMISTRY_GEAR)
#include "./chemistry/GEAR/chemistry_struct.h" #include "./chemistry/GEAR/chemistry_struct.h"
#elif defined(CHEMISTRY_QLA)
#include "./chemistry/QLA/chemistry_struct.h"
#elif defined(CHEMISTRY_EAGLE) #elif defined(CHEMISTRY_EAGLE)
#include "./chemistry/EAGLE/chemistry_struct.h" #include "./chemistry/EAGLE/chemistry_struct.h"
#else #else
......
...@@ -45,6 +45,8 @@ ...@@ -45,6 +45,8 @@
#include "./cooling/Compton/cooling.h" #include "./cooling/Compton/cooling.h"
#elif defined(COOLING_GRACKLE) #elif defined(COOLING_GRACKLE)
#include "./cooling/grackle/cooling.h" #include "./cooling/grackle/cooling.h"
#elif defined(COOLING_QLA)
#include "./cooling/QLA/cooling.h"
#elif defined(COOLING_EAGLE) #elif defined(COOLING_EAGLE)
#include "./cooling/EAGLE/cooling.h" #include "./cooling/EAGLE/cooling.h"
#else #else
......
This diff is collapsed.
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2020 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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_COOLING_QLA_H
#define SWIFT_COOLING_QLA_H
/**
* @file src/cooling/EAGLE/cooling.h
* @brief EAGLE cooling function declarations
*/
/* Local includes. */
#include "cooling_struct.h"
struct part;
struct xpart;
struct cosmology;
struct hydro_props;
struct entropy_floor_properties;
struct space;
void cooling_update(const struct cosmology *cosmo,
struct cooling_function_data *cooling, struct space *s);
void cooling_cool_part(const struct phys_const *phys_const,
const struct unit_system *us,
const struct cosmology *cosmo,
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);
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,
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);
void cooling_split_part(struct part *p, struct xpart *xp, double n);
void cooling_Hydrogen_reionization(const struct cooling_function_data *cooling,
const struct cosmology *cosmo,
struct space *s);
void cooling_init_backend(struct swift_params *parameter_file,
const struct unit_system *us,
const struct phys_const *phys_const,
const struct hydro_props *hydro_props,
struct cooling_function_data *cooling);
void cooling_print_backend(const struct cooling_function_data *cooling);
void cooling_clean(struct cooling_function_data *data);
#endif /* SWIFT_COOLING_QLA_H */
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2020 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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_COOLING_QLA_IO_H
#define SWIFT_COOLING_QLA_IO_H
/* Config parameters. */
#include "../config.h"
/* Local includes */
#include "cooling.h"
#include "io_properties.h"
#ifdef HAVE_HDF5
/**
* @brief Writes the current model of cooling to the file.
*
* @param h_grp The HDF5 group in which to write
* @param cooling The #cooling_function_data
*/
__attribute__((always_inline)) INLINE static void cooling_write_flavour(
hid_t h_grp, const struct cooling_function_data* cooling) {
io_write_attribute_s(h_grp, "Cooling Model",
"Quick Lyman-alpha (EAGLE with primordial Z only)");
}
#endif
INLINE static void convert_part_T(const struct engine* e, const struct part* p,
const struct xpart* xp, float* ret) {
ret[0] = cooling_get_temperature(e->physical_constants, e->hydro_properties,
e->internal_units, e->cosmology,
e->cooling_func, p, xp);
}
/**
* @brief Specifies which particle fields to write to a dataset
*
* @param parts The particle array.
* @param xparts The extended data particle array.
* @param list The list of i/o properties to write.
* @param cooling The #cooling_function_data
*
* @return Returns the number of fields to write.
*/
__attribute__((always_inline)) INLINE static int cooling_write_particles(
const struct part* parts, const struct xpart* xparts, struct io_props* list,
const struct cooling_function_data* cooling) {
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;
}
#endif /* SWIFT_COOLING_QLA_IO_H */
This diff is collapsed.
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2020 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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_COOLING_STRUCT_QLA_H
#define SWIFT_COOLING_STRUCT_QLA_H
#define qla_table_path_name_length 500
/**
* @brief struct containing cooling tables
*/
struct cooling_tables {
/* array of heating rates due to metals */
float *metal_heating;
/* array of heating rates due to hydrogen and helium */
float *H_plus_He_heating;
/* array of electron abundances due to hydrogen and helium */
float *H_plus_He_electron_abundance;
/* array of temperatures */
float *temperature;
/* array of electron abundances due to metals */
float *electron_abundance;
};
/**
* @brief Properties of the cooling function.
*/
struct cooling_function_data {
/*! Cooling tables */
struct cooling_tables table;
/*! Redshift bins */
float *Redshifts;
/*! Hydrogen number density bins */
float *nH;
/*! Temperature bins */
float *Temp;
/*! Helium fraction bins */
float *HeFrac;
/*! Internal energy bins */
float *Therm;
/*! Mass fractions of elements for solar abundances (from the tables) */
float *SolarAbundances;
/*! Inverse of the solar mass fractions */
float *SolarAbundances_inv;
/*! Filepath to the directory containing the HDF5 cooling tables */
char cooling_table_path[qla_table_path_name_length];
/*! Redshift of H reionization */
float H_reion_z;
/*! H reionization energy in CGS units */
float H_reion_heat_cgs;
/*! Have we already done H reioisation? */
int H_reion_done;
/*! Ca over Si abundance divided by the solar ratio for these elements */
float Ca_over_Si_ratio_in_solar;
/*! S over Si abundance divided by the solar ratio for these elements */
float S_over_Si_ratio_in_solar;
/*! Redshift of He reionization */
float He_reion_z_centre;
/*! Spread of the He reionization */
float He_reion_z_sigma;
/*! He reionization energy in CGS units */
float He_reion_heat_cgs;
/*! Internal energy conversion from internal units to CGS (for quick access)
*/
double internal_energy_to_cgs;
/*! Internal energy conversion from CGS to internal units (for quick access)
*/
double internal_energy_from_cgs;
/*! Number density conversion from internal units to CGS (for quick access) */
double number_density_to_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) */
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;
/*! Distance between the current redshift and table[z_index] */
float dz;
/*! Index of the previous tables along the redshift index of the tables */
int previous_z_index;
};
/**
* @brief Properties of the cooling stored in the extended particle data.
*/
struct cooling_xpart_data {};
#endif /* SWIFT_COOLING_STRUCT_QLA_H */
This diff is collapsed.
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2020 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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_QLA_COOL_TABLES_H
#define SWIFT_QLA_COOL_TABLES_H
/**
* @file src/cooling/QLA/cooling.h
* @brief Quick Lyman-alpha cooling function (EAGLE with primordial Z only)
*/
/* Config parameters. */
#include "../config.h"
#include "cooling_struct.h"
/*! Number of different bins along the redhsift axis of the tables */
#define qla_cooling_N_redshifts 49
/*! Number of redshift bins loaded at any given point int time */
#define qla_cooling_N_loaded_redshifts 2
/*! Number of different bins along the temperature axis of the tables */
#define qla_cooling_N_temperature 176
/*! Number of different bins along the density axis of the tables */
#define qla_cooling_N_density 41
/*! Number of different bins along the metal axis of the tables */
#define qla_cooling_N_metal 9
/*! Number of different bins along the metal axis of the tables */
#define qla_cooling_N_He_frac 7
/*! Number of different bins along the abundances axis of the tables */
#define qla_cooling_N_abundances 11
void get_cooling_redshifts(struct cooling_function_data *cooling);
void read_cooling_header(const char *fname,
struct cooling_function_data *cooling);
void allocate_cooling_tables(struct cooling_function_data *restrict cooling);
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);
#endif /* SWIFT_QLA_COOL_TABLES_H */
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2020 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 <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_INTERPOL_QLA_H
#define SWIFT_INTERPOL_QLA_H
/**
* @file src/cooling/QLA/interpolate.h
* @brief Interpolation functions for Quick Lyman-alpha cooling tables
*/
/* Config parameters. */
#include "../config.h"
/* Local includes. */
#include "align.h"
#include "error.h"
#include "inline.h"
/**
* @brief Returns the 1d index of element with 2d indices x,y
* from a flattened 2d array in row major order
*
* @param x, y Indices of element of interest
* @param Nx, Ny Sizes of array dimensions
*/
__attribute__((always_inline)) INLINE int row_major_index_2d(const int x,
const int y,
const int Nx,
const int Ny) {
#ifdef SWIFT_DEBUG_CHECKS
assert(x < Nx);
assert(y < Ny);
#endif
return x * Ny + y;
}
/**
* @brief Returns the 1d index of element with 3d indices x,y,z
* from a flattened 3d array in row major order
*
* @param x, y, z Indices of element of interest
* @param Nx, Ny, Nz Sizes of array dimensions
*/
__attribute__((always_inline)) INLINE int row_major_index_3d(
const int x, const int y, const int z, const int Nx, const int Ny,
const int Nz) {
#ifdef SWIFT_DEBUG_CHECKS
assert(x < Nx);
assert(y < Ny);
assert(z < Nz);
#endif
return x * Ny * Nz + y * Nz + z;
}
/**
* @brief Returns the 1d index of element with 4d indices x,y,z,w
* from a flattened 4d array in row major order
*
* @param x, y, z, w Indices of element of interest
* @param Nx, Ny, Nz, Nw Sizes of array dimensions
*/
__attribute__((always_inline)) INLINE int row_major_index_4d(
const int x, const int y, const int z, const int w, const int Nx,
const int Ny, const int Nz, const int Nw) {
#ifdef SWIFT_DEBUG_CHECKS
assert(x < Nx);
assert(y < Ny);
assert(z < Nz);
assert(w < Nw);
#endif
return x * Ny * Nz * Nw + y * Nz * Nw + z * Nw + w;
}
/**
* @brief Finds the index of a value in a table and compute delta to nearest
* element.
*
* This function assumes the table is monotonically increasing with a constant
* difference between adjacent values.
*
* The returned difference is expressed in units of the table separation. This
* means dx = (x - table[i]) / (table[i+1] - table[i]). It is always between
* 0 and 1.
*
* We use a small epsilon of 1e-4 to avoid out-of-range accesses due to
* rounding errors.
*
* @param table The table to search in.
* @param size The number of elements in the table.
* @param x The value to search for.
* @param i (return) The index in the table of the element.
* @param *dx (return) The difference between x and table[i]
*/
__attribute__((always_inline)) INLINE void get_index_1d(
const float *restrict table, const int size, const float x, int *i,
float *restrict dx) {
/* Small epsilon to avoid rounding issues leading to out-of-bound
* access when using the indices later to read data from the tables. */
const float epsilon = 1e-4f;
/* Indicate that the whole array is aligned on boundaries */
swift_align_information(float, table, SWIFT_STRUCT_ALIGNMENT);
/* Distance between elements in the array */
const float delta = (size - 1) / (table[size - 1] - table[0]);
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;
#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
*dx = (x - table[*i]) * delta;
} else {
/* We are after the last element */
*i = size - 2;
*dx = 1.f;
}
#ifdef SWIFT_DEBUG_CHECKS
if (*dx < -0.001f || *dx > 1.001f) error("Invalid distance found dx=%e", *dx);
#endif
}
/**
* @brief Interpolate a flattened 2D table at a given position.
*
* This function uses linear interpolation along each axis. It also
* assumes that the table is aligned on SWIFT_STRUCT_ALIGNMENT.
*
* @param table The 2D table to interpolate.
* @param xi, yi Indices of element of interest.
* @param Nx, Ny Sizes of array dimensions.
* @param dx, dy Distance between the point and the index in units of
* the grid spacing.
*/
__attribute__((always_inline)) INLINE float interpolation_2d(
const float *table, const int xi, const int yi, const float dx,
const float dy, const int Nx, const int Ny) {
#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);
#endif
const float tx = 1.f - dx;
const float ty = 1.f - dy;
/* Indicate that the whole array is aligned on boundaries */
swift_align_information(float, table, SWIFT_STRUCT_ALIGNMENT);
/* Linear interpolation along each axis. We read the table 2^2=4 times */
float result = tx * ty * table[row_major_index_2d(xi + 0, yi + 0, Nx, Ny)];
result += tx * dy * table[row_major_index_2d(xi + 0, yi + 1, Nx, Ny)];
result += dx * ty * table[row_major_index_2d(xi + 1, yi + 0, Nx, Ny)];
result += dx * dy * table[row_major_index_2d(xi + 1, yi + 1, Nx, Ny)];
return result;
}
/**
* @brief Interpolate a flattened 3D table at a given position.
*
* This function uses linear interpolation along each axis. It 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(
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.001f || dz > 1.001f) error("Invalid dz=%e", dz);
#endif
const float tx = 1.f - dx;
const float ty = 1.f - dy;
const float tz = 1.f - dz;
/* 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^3=8 times */
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.
*
* This function uses linear interpolation along each axis.
* We look at the xi 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_x(
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.f) error("Attempting to interpolate along x!");
if (dy < -0.001f || dy > 1.001f) error("Invalid dy=%e", dy);
if (dz < -0.001f || dz > 1.001f) error("Invalid dz=%e", dz);
#endif
const float tx = 1.f;
const float ty = 1.f - dy;
const float tz = 1.f - dz;
/* 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 4D table at a given position.
*
* This function uses linear interpolation along each axis. It also
* assumes that the table is aligned on SWIFT_STRUCT_ALIGNMENT.
*
* @param table The 4D table to interpolate.
* @param xi, yi, zi, wi Indices of element of interest.
* @param Nx, Ny, Nz, Nw Sizes of array dimensions.
* @param dx, dy, dz, dw Distance between the point and the index in units of
* the grid spacing.
*/
__attribute__((always_inline)) INLINE float interpolation_4d(
const float *table, const int xi, const int yi, const int zi, const int wi,
const float dx, const float dy, const float dz, const float dw,
const int Nx, const int Ny, const int Nz, const int Nw) {
#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.001f || dz > 1.001f) error("Invalid dz=%e", dz);
if (dw < -0.001f || dw > 1.001f) error("Invalid dw=%e", dw);
#endif
const float tx = 1.f - dx;
const float ty = 1.f - dy;
const float tz = 1.f - dz;
const float tw = 1.f - dw;
/* 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^4=16 times */
float result =
tx * ty * tz * tw *
table[row_major_index_4d(xi + 0, yi + 0, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
result +=
tx * ty * tz * dw *
table[row_major_index_4d(xi + 0, yi + 0, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
result +=
tx * ty * dz * tw *
table[row_major_index_4d(xi + 0, yi + 0, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
result +=
tx * dy * tz * tw *
table[row_major_index_4d(xi + 0, yi + 1, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
result +=
dx * ty * tz * tw *
table[row_major_index_4d(xi + 1, yi + 0, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
result +=
tx * ty * dz * dw *
table[row_major_index_4d(xi + 0, yi + 0, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
result +=
tx * dy * tz * dw *
table[row_major_index_4d(xi + 0, yi + 1, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
result +=
dx * ty * tz * dw *
table[row_major_index_4d(xi + 1, yi + 0, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
result +=
tx * dy * dz * tw *
table[row_major_index_4d(xi + 0, yi + 1, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
result +=
dx * ty * dz * tw *
table[row_major_index_4d(xi + 1, yi + 0, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
result +=
dx * dy * tz * tw *
table[row_major_index_4d(xi + 1, yi + 1, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
result +=
dx * dy * dz * tw *
table[row_major_index_4d(xi + 1, yi + 1, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
result +=
dx * dy * tz * dw *
table[row_major_index_4d(xi + 1, yi + 1, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
result +=
dx * ty * dz * dw *
table[row_major_index_4d(xi + 1, yi + 0, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
result +=
tx * dy * dz * dw *
table[row_major_index_4d(xi + 0, yi + 1, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
result +=
dx * dy * dz * dw *
table[row_major_index_4d(xi + 1, yi + 1, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
return result;
}
/**
* @brief Interpolate a flattened 4D table at a given position but avoid the
* x-dimension.
*
* This function uses linear interpolation along each axis.
* We look at the xi coordoniate but do not interpolate around it. We just
* interpolate the remaining 3 dimensions.
* The function also assumes that the table is aligned on
* SWIFT_STRUCT_ALIGNMENT.
*
* @param table The 4D table to interpolate.
* @param xi, yi, zi, wi Indices of element of interest.
* @param Nx, Ny, Nz, Nw Sizes of array dimensions.
* @param dx, dy, dz, dw Distance between the point and the index in units of
* the grid spacing.
*/
__attribute__((always_inline)) INLINE float interpolation_4d_no_x(
const float *table, const int xi, const int yi, const int zi, const int wi,
const float dx, const float dy, const float dz, const float dw,
const int Nx, const int Ny, const int Nz, const int Nw) {
#ifdef SWIFT_DEBUG_CHECKS
if (dx != 0.f) error("Attempting to interpolate along x!");
if (dy < -0.001f || dy > 1.001f) error("Invalid dy=%e", dy);
if (dz < -0.001f || dz > 1.001f) error("Invalid dz=%e", dz);
if (dw < -0.001f || dw > 1.001f) error("Invalid dw=%e", dw);
#endif
const float tx = 1.f;
const float ty = 1.f - dy;
const float tz = 1.f - dz;
const float tw = 1.f - dw;
/* Indicate that the whole array is aligned on boundaries */
swift_align_information(float, table, SWIFT_STRUCT_ALIGNMENT);
/* Linear interpolation along each axis. We read the table 2^3=8 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 * tw *
table[row_major_index_4d(xi + 0, yi + 0, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
result +=
tx * ty * tz * dw *
table[row_major_index_4d(xi + 0, yi + 0, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
result +=
tx * ty * dz * tw *
table[row_major_index_4d(xi + 0, yi + 0, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
result +=
tx * dy * tz * tw *
table[row_major_index_4d(xi + 0, yi + 1, zi + 0, wi + 0, Nx, Ny, Nz, Nw)];
/* result += */
/* dx * ty * tz * tw * */
/* table[row_major_index_4d(xi + 1, yi + 0, zi + 0, wi + 0, Nx, Ny, Nz,
* Nw)]; */
result +=
tx * ty * dz * dw *
table[row_major_index_4d(xi + 0, yi + 0, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
result +=
tx * dy * tz * dw *
table[row_major_index_4d(xi + 0, yi + 1, zi + 0, wi + 1, Nx, Ny, Nz, Nw)];
/* result += */
/* dx * ty * tz * dw * */
/* table[row_major_index_4d(xi + 1, yi + 0, zi + 0, wi + 1, Nx, Ny, Nz,
* Nw)]; */
result +=
tx * dy * dz * tw *
table[row_major_index_4d(xi + 0, yi + 1, zi + 1, wi + 0, Nx, Ny, Nz, Nw)];
/* result += */
/* dx * ty * dz * tw * */
/* table[row_major_index_4d(xi + 1, yi + 0, zi + 1, wi + 0, Nx, Ny, Nz,
* Nw)]; */
/* result += */
/* dx * dy * tz * tw * */
/* table[row_major_index_4d(xi + 1, yi + 1, zi + 0, wi + 0, Nx, Ny, Nz, */
/* Nw)]; */
/* result += */
/* dx * dy * dz * tw * */
/* table[row_major_index_4d(xi + 1, yi + 1, zi + 1, wi + 0, Nx, Ny, Nz, */
/* Nw)]; */
/* result += */
/* dx * dy * tz * dw * */
/* table[row_major_index_4d(xi + 1, yi + 1, zi + 0, wi + 1, Nx, Ny, Nz, */
/* Nw)]; */
/* result += */
/* dx * ty * dz * dw * */
/* table[row_major_index_4d(xi + 1, yi + 0, zi + 1, wi + 1, Nx, Ny, Nz,
* Nw)]; */
result +=
tx * dy * dz * dw *
table[row_major_index_4d(xi + 0, yi + 1, zi + 1, wi + 1, Nx, Ny, Nz, Nw)];
/* result += */
/* dx * dy * dz * dw * */
/* table[row_major_index_4d(xi + 1, yi + 1, zi + 1, wi + 1, Nx, Ny, Nz, */
/* Nw)]; */
return result;
}
#endif
...@@ -33,6 +33,8 @@ ...@@ -33,6 +33,8 @@
#include "./cooling/Compton/cooling_io.h" #include "./cooling/Compton/cooling_io.h"
#elif defined(COOLING_GRACKLE) #elif defined(COOLING_GRACKLE)
#include "./cooling/grackle/cooling_io.h" #include "./cooling/grackle/cooling_io.h"
#elif defined(COOLING_QLA)
#include "./cooling/QLA/cooling_io.h"
#elif defined(COOLING_EAGLE) #elif defined(COOLING_EAGLE)
#include "./cooling/EAGLE/cooling_io.h" #include "./cooling/EAGLE/cooling_io.h"
#else #else
......
...@@ -38,6 +38,8 @@ ...@@ -38,6 +38,8 @@
#include "./cooling/Compton/cooling_struct.h" #include "./cooling/Compton/cooling_struct.h"
#elif defined(COOLING_GRACKLE) #elif defined(COOLING_GRACKLE)
#include "./cooling/grackle/cooling_struct.h" #include "./cooling/grackle/cooling_struct.h"
#elif defined(COOLING_QLA)
#include "./cooling/QLA/cooling_struct.h"
#elif defined(COOLING_EAGLE) #elif defined(COOLING_EAGLE)
#include "./cooling/EAGLE/cooling_struct.h" #include "./cooling/EAGLE/cooling_struct.h"
#else #else
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment