diff --git a/doc/RTD/source/SubgridModels/GEAR/index.rst b/doc/RTD/source/SubgridModels/GEAR/index.rst
index 82ba003a0ee2ed0916e115c71b1da6edb3cd4150..c1ac89f70be3bc2a75a99e2fb02ef515d7d67d5d 100644
--- a/doc/RTD/source/SubgridModels/GEAR/index.rst
+++ b/doc/RTD/source/SubgridModels/GEAR/index.rst
@@ -55,13 +55,21 @@ and ``grackle_3`` (the numbers correspond to the value of
methods and UV background (on/off with ``GrackleCooling:with_UV_background``). In order to use the Grackle cooling, you will need
to provide a HDF5 table computed by Cloudy (``GrackleCooling:cloudy_table``).
+
When starting a simulation without providing the different element fractions in the non equilibrium mode, the code supposes an equilibrium and computes them automatically.
The code uses an iterative method in order to find the correct initial composition and this method can be tuned with two parameters. ``GrackleCooling:max_steps`` defines the maximal number of steps to reach the convergence and ``GrackleCooling:convergence_limit`` defines the tolerance in the relative error.
In order to compile SWIFT with Grackle, you need to provide the options ``with-chemistry=GEAR`` and ``with-grackle=$GRACKLE_ROOT``
where ``$GRACKLE_ROOT`` is the root of the install directory (not the ``lib``).
-You will need a Grackle version later than 3.1.1. To compile it, run
+.. warning::
+ The actual Grackle version fully supported by SWIFT is 3.2.1. It can be downloaded from
+ `the official Grackle git repository <https://github.com/grackle-project/grackle/archive/refs/tags/grackle-3.2.1.tar.gz>`_.
+ However, this version still had a bug when using threadsafe functions. Alternately, it is possible to get a fixed version
+ using `the following fork frozen for compatibility with SWIFT <https://github.com/mladenivkovic/grackle-swift>`_.
+
+
+To compile it, run
the following commands from the root directory of Grackle:
``./configure; cd src/clib``.
Update the variables ``LOCAL_HDF5_INSTALL`` and ``MACH_INSTALL_PREFIX`` in
@@ -95,6 +103,10 @@ The self shielding method is defined by ``GrackleCooling:self_shielding_method``
self_shielding_method: -1 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
+.. note::
+ A simple example running SWIFT with Grackle can be find in ``examples/Cooling/CoolingBox``. A more advanced example combining heating and cooling (with heating and ionization sources) is given in ``examples/Cooling/CoolingHeatingBox``.
+
+
.. _gear_star_formation:
Star formation
diff --git a/examples/Cooling/CoolingHeatingBox/README b/examples/Cooling/CoolingHeatingBox/README
new file mode 100644
index 0000000000000000000000000000000000000000..999269bf9d7667f5e1b5d49debf0af227299c5c5
--- /dev/null
+++ b/examples/Cooling/CoolingHeatingBox/README
@@ -0,0 +1,7 @@
+Runs a uniform box exposed to constant heating and ionization rates.
+This tests is the Iliev+2006 (ui.adsabs.harvard.edu/abs/2006MNRAS.369.1625I/)
+Test 0 part 3, however without stopping the heating and ionization rates
+after 0.5 Myr. While the cooling is present, gas heats up and ionize.
+
+To run with GEAR-RT, configure swift with:
+ ./configure --with-cooling=grackle_1 --with-grackle=$GRACKLE_ROOT
diff --git a/examples/Cooling/CoolingHeatingBox/coolingBox.yml b/examples/Cooling/CoolingHeatingBox/coolingBox.yml
new file mode 100644
index 0000000000000000000000000000000000000000..362a3c958061edb59de54ead05c7ab59c7d0fe40
--- /dev/null
+++ b/examples/Cooling/CoolingHeatingBox/coolingBox.yml
@@ -0,0 +1,72 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.98841e43 # 1e10 Msol #2.0e33 # Solar masses
+ UnitLength_in_cgs: 3.08567758e24 # Mpc #3.0857e21 # Kiloparsecs
+ UnitVelocity_in_cgs: 1.0e5 # Kilometers per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 5e-6 # The end time of the simulation (in internal units).
+ dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-8 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: snap/snapshot # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1e-7 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-3 # Time between statistics output
+
+Scheduler:
+ tasks_per_cell: 64
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+ minimal_temperature: 100. # Kelvin
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./coolingBox.hdf5 # The file to read
+ periodic: 1
+
+# Dimensionless pre-factor for the time-step condition
+LambdaCooling:
+ lambda_nH2_cgs: 1e-22 # Cooling rate divided by square Hydrogen number density (in cgs units [erg * s^-1 * cm^3])
+ cooling_tstep_mult: 1.0 # Dimensionless pre-factor for the time-step condition
+
+# Cooling with Grackle 3.0
+GrackleCooling:
+ cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
+ with_UV_background: 0 # Enable or not the UV background
+ redshift: 0 # Redshift to use (-1 means time based redshift)
+ with_metal_cooling: 0 # Enable or not the metal cooling
+ max_steps: 10000 # (optional) Max number of step when computing the initial composition
+ convergence_limit: 1e-2 # (optional) Convergence threshold (relative) for initial composition
+ thermal_time_myr: 5 # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
+ self_shielding_method: 0 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
+ self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
+ HydrogenFractionByMass : 1. # Hydrogen fraction by mass (default is 0.76)
+ use_radiative_transfer : 1 # Arrays of ionization and heating rates are provided
+ RT_heating_rate_cgs : 1.65117e-17 # heating rate in units of / nHI_cgs (corresponding to a flux of 10^12 photons s−1 cm−2 , with a 10^5 K blackbody spectrum)
+ RT_HI_ionization_rate_cgs : 1.63054e-06 # HI ionization rate in cgs [1/s] (corresponding to a flux of 10^12 photons s−1 cm−2 , with a 10^5 K blackbody spectrum)
+ RT_HeI_ionization_rate_cgs : 0 # HeI ionization rate in cgs [1/s]
+ RT_HeII_ionization_rate_cgs: 0 # HeII ionization rate in cgs [1/s]
+ RT_H2_dissociation_rate_cgs: 0 # H2 dissociation rate in cgs [1/s]
+
+
+GEARChemistry:
+ initial_metallicity: 0.01295
+
+
+
+GEARPressureFloor:
+ jeans_factor: 0. # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+
diff --git a/examples/Cooling/CoolingHeatingBox/getGlass.sh b/examples/Cooling/CoolingHeatingBox/getGlass.sh
new file mode 100755
index 0000000000000000000000000000000000000000..ffd92e88deae6e91237059adac2a6c2067caee46
--- /dev/null
+++ b/examples/Cooling/CoolingHeatingBox/getGlass.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/glassCube_32.hdf5
diff --git a/examples/Cooling/CoolingHeatingBox/getResults.sh b/examples/Cooling/CoolingHeatingBox/getResults.sh
new file mode 100755
index 0000000000000000000000000000000000000000..d178c5952f8f63c8d1fd28648e4508728a2ff567
--- /dev/null
+++ b/examples/Cooling/CoolingHeatingBox/getResults.sh
@@ -0,0 +1,2 @@
+#!/bin/bash
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ReferenceSolutions/CoolingHeatingBox_results.txt
diff --git a/examples/Cooling/CoolingHeatingBox/makeIC.py b/examples/Cooling/CoolingHeatingBox/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..8746a96204a0145296b0a5d155b2fd81024cb56c
--- /dev/null
+++ b/examples/Cooling/CoolingHeatingBox/makeIC.py
@@ -0,0 +1,122 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Yves Revaz (yves.revaz@epfl.ch)
+# Loic Hausammann (loic.hausammann@id.ethz.ch)
+#
+# 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/>.
+#
+##############################################################################
+
+import h5py
+import numpy as np
+
+# Generates a SWIFT IC file with a constant density and pressure
+
+# Parameters
+periodic = 1 # 1 For periodic box
+boxSize = 1 # 1 kiloparsec
+# rho = 3.2e3 # Density in code units (3.2e6 is 0.1 hydrogen atoms per cm^3)
+
+# Iliev test values
+rho = (
+ 2471402.49842514
+) # Density in code units (1-hydrogen atoms per cm^3 assuming Hydrogen only)
+T = 100 # Initial Temperature
+
+
+gamma = 5.0 / 3.0 # Gas adiabatic index
+fileName = "coolingBox.hdf5"
+# ---------------------------------------------------
+
+# defines some constants
+# need to be changed in plotResults.py too
+h_frac = 1.0 # 0.76
+mu = 4.0 / (1.0 + 3.0 * h_frac)
+
+m_h_cgs = 1.67262192369e-24
+k_b_cgs = 1.380649e-16
+
+# defines units
+unit_length = 3.08567758e24 # Mpc
+unit_mass = 1.98841e43 # 1e10 solar mass
+unit_velocity = 1e5 # 1e5 km/s
+unit_time = unit_length / unit_velocity
+
+# Read id, position and h from glass
+glass = h5py.File("glassCube_32.hdf5", "r")
+ids = glass["/PartType0/ParticleIDs"][:]
+pos = glass["/PartType0/Coordinates"][:, :] * boxSize
+h = glass["/PartType0/SmoothingLength"][:] * boxSize
+
+# Compute basic properties
+numPart = np.size(pos) // 3
+mass = boxSize ** 3 * rho / numPart
+internalEnergy = k_b_cgs * T * mu / ((gamma - 1.0) * m_h_cgs)
+internalEnergy *= (unit_time / unit_length) ** 2
+
+# File
+f = h5py.File(fileName, "w")
+
+# Header
+grp = f.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+
+# Runtime parameters
+grp = f.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+# Units
+grp = f.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = unit_length
+grp.attrs["Unit mass in cgs (U_M)"] = unit_mass
+grp.attrs["Unit time in cgs (U_t)"] = unit_time
+grp.attrs["Unit current in cgs (U_I)"] = 1.0
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
+
+# Particle group
+grp = f.create_group("/PartType0")
+
+v = np.zeros((numPart, 3))
+ds = grp.create_dataset("Velocities", (numPart, 3), "f")
+ds[()] = v
+
+m = np.full((numPart, 1), mass)
+ds = grp.create_dataset("Masses", (numPart, 1), "f")
+ds[()] = m
+
+h = np.reshape(h, (numPart, 1))
+ds = grp.create_dataset("SmoothingLength", (numPart, 1), "f")
+ds[()] = h
+
+u = np.full((numPart, 1), internalEnergy)
+ds = grp.create_dataset("InternalEnergy", (numPart, 1), "f")
+ds[()] = u
+
+ids = np.reshape(ids, (numPart, 1))
+ds = grp.create_dataset("ParticleIDs", (numPart, 1), "L")
+ds[()] = ids
+
+ds = grp.create_dataset("Coordinates", (numPart, 3), "d")
+ds[()] = pos
+
+f.close()
+
+print("Initial condition generated")
diff --git a/examples/Cooling/CoolingHeatingBox/plotResults.py b/examples/Cooling/CoolingHeatingBox/plotResults.py
new file mode 100644
index 0000000000000000000000000000000000000000..2609ddf45f85c05dcf4447f0e6dd263cf4a21ba8
--- /dev/null
+++ b/examples/Cooling/CoolingHeatingBox/plotResults.py
@@ -0,0 +1,118 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2023 Yves Revaz (yves.revaz@epfl.ch)
+# Loic Hausammann (loic.hausammann@id.ethz.ch)
+#
+# 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/>.
+#
+##############################################################################
+import matplotlib
+
+import matplotlib.pyplot as plt
+import numpy as np
+from glob import glob
+import h5py
+
+plt.style.use("../../../tools/stylesheets/mnras.mplstyle")
+
+
+##########################################
+# read specific energy from the solution
+##########################################
+
+# Read in units.
+resultfile = "CoolingHeatingBox_results.txt"
+
+f = open(resultfile, "r")
+firstline = f.readline()
+massline = f.readline()
+lengthline = f.readline()
+velline = f.readline()
+f.close()
+units = []
+for l in [massline, lengthline, velline]:
+ before, after = l.split("used:")
+ val, unit = after.split("[")
+ val = val.strip()
+ units.append(float(val))
+
+mass_units = units[0]
+length_units = units[1]
+velocity_units = units[2]
+time_units = velocity_units / length_units
+density_units = mass_units / length_units ** 3
+
+# Read in all other data
+data = np.loadtxt(resultfile)
+
+Time = data[:, 1]
+Time_Myr = Time * 1e-6
+Energy = data[:, 12] # internal energy IU
+
+
+##########################################
+# compute specific energy from the models
+##########################################
+
+
+# First snapshot
+snap_filename = "snap/snapshot_0000.hdf5"
+
+# Read the initial state of the gas
+f = h5py.File(snap_filename, "r")
+
+# Read the units parameters from the snapshot
+units = f["InternalCodeUnits"]
+unit_mass = units.attrs["Unit mass in cgs (U_M)"]
+unit_length = units.attrs["Unit length in cgs (U_L)"]
+unit_time = units.attrs["Unit time in cgs (U_t)"]
+
+
+# Read snapshots
+files = glob("snap/snapshot_*.hdf5")
+files.sort()
+N = len(files)
+energy_snap = np.zeros(N)
+time_snap_cgs = np.zeros(N)
+for i in range(N):
+ snap = h5py.File(files[i], "r")
+ masses = snap["/PartType0/Masses"][:]
+ u = snap["/PartType0/InternalEnergies"][:] * masses
+ u = sum(u) / masses.sum()
+ energy_snap[i] = u
+ time_snap_cgs[i] = snap["/Header"].attrs["Time"] * unit_time
+
+
+##########################################
+# plot
+##########################################
+
+plt.figure()
+
+Myr_in_s = 1e6 * 365.25 * 24.0 * 60.0 * 60.0
+yr_in_s = 365.25 * 24.0 * 60.0 * 60.0
+
+plt.plot(Time_Myr, Energy, "r", ms=3, label="Expected solution")
+
+plt.plot(time_snap_cgs / Myr_in_s, energy_snap, "k", ms=2)
+
+
+plt.legend(loc="right", fontsize=8, frameon=False, handlelength=3, ncol=1)
+plt.xlabel("${\\rm{Time~[Myr]}}$", labelpad=0)
+plt.ylabel(r"$\rm{Internal Energy\,\,[IU]}$")
+
+
+plt.tight_layout()
+
+plt.show()
diff --git a/examples/Cooling/CoolingHeatingBox/run.sh b/examples/Cooling/CoolingHeatingBox/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..466c2281b987d5adf76555890383db3548ff1d14
--- /dev/null
+++ b/examples/Cooling/CoolingHeatingBox/run.sh
@@ -0,0 +1,38 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e glassCube_32.hdf5 ]
+then
+ echo "Fetching initial glass file for the cooling box example..."
+ ./getGlass.sh
+fi
+if [ ! -e coolingBox.hdf5 ]
+then
+ echo "Generating initial conditions for the cooling box example..."
+ python3 makeIC.py
+fi
+
+# Get the Grackle cooling table
+if [ ! -e CloudyData_UVB=HM2012.h5 ]
+then
+ echo "Fetching the Cloudy tables required by Grackle..."
+ ../getGrackleCoolingTable.sh
+fi
+
+# Get the results
+if [ ! -e CoolingHeatingBox_results.txt]
+then
+ echo "Fetching the the results..."
+ ./getResults.sh
+fi
+
+
+# Create output directory
+rm -rf snap
+mkdir snap
+
+# Run SWIFT
+../../../swift --hydro --cooling --threads=14 coolingBox.yml
+
+# Check energy conservation and cooling rate
+python3 plotResults.py
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 7b1eefbf50a7f279652b136942ce5b614758d7b1..b10e6d82e79947a4b568ce2bbe68229edcf52700 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -504,6 +504,18 @@ GrackleCooling:
thermal_time_myr: 5 # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
self_shielding_method: -1 # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
+ HydrogenFractionByMass : 1. # Hydrogen fraction by mass (default is 0.76)
+ use_radiative_transfer : 1 # Arrays of ionization and heating rates are provided
+ RT_heating_rate_cgs : 0 # heating rate in units of / nHI_cgs
+ RT_HI_ionization_rate_cgs : 0 # HI ionization rate in cgs [1/s]
+ RT_HeI_ionization_rate_cgs : 0 # HeI ionization rate in cgs [1/s]
+ RT_HeII_ionization_rate_cgs: 0 # HeII ionization rate in cgs [1/s]
+ RT_H2_dissociation_rate_cgs: 0 # H2 dissociation rate in cgs [1/s]
+ volumetric_heating_rates_cgs: 0 # Volumetric heating rate in cgs [erg/s/cm3]
+ specific_heating_rates_cgs: 0 # Specific heating rate in cgs [erg/s/g]
+
+
+
# Parameters related to chemistry models -----------------------------------------------
diff --git a/src/cooling/EAGLE/cooling.c b/src/cooling/EAGLE/cooling.c
index 4799488ff3318c7cdb3e592c4e27440a0815dad6..11ab3d97def84b2fd32f55ce23513595615a4d4e 100644
--- a/src/cooling/EAGLE/cooling.c
+++ b/src/cooling/EAGLE/cooling.c
@@ -126,10 +126,13 @@ __attribute__((always_inline)) INLINE void get_redshift_index(
* @param pressure_floor Properties of the pressure floor.
* @param cooling The #cooling_function_data used in the run.
* @param s The space data, including a pointer to array of particles
+ * @param time The current system time
*/
-void cooling_update(const struct cosmology *cosmo,
+void cooling_update(const struct phys_const *phys_const,
+ const struct cosmology *cosmo,
const struct pressure_floor_props *pressure_floor,
- struct cooling_function_data *cooling, struct space *s) {
+ struct cooling_function_data *cooling, struct space *s,
+ const double time) {
/* Current redshift */
const float redshift = cosmo->z;
@@ -573,6 +576,28 @@ __attribute__((always_inline)) INLINE void cooling_first_init_part(
xp->cooling_data.radiated_energy = 0.f;
}
+/**
+ * @brief Perform additional init on the cooling properties of the
+ * (x-)particles that requires the density to be known.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cosmo The current cosmological model.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE void cooling_post_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 *cooling, const struct part *restrict p,
+ struct xpart *restrict xp) {}
+
/**
* @brief Compute the temperature based on gas properties.
*
@@ -1059,7 +1084,8 @@ void cooling_restore_tables(struct cooling_function_data *cooling,
/* Force a re-read of the cooling tables */
cooling->z_index = -10;
cooling->previous_z_index = eagle_cooling_N_redshifts - 2;
- cooling_update(cosmo, /*pfloor=*/NULL, cooling, /*space=*/NULL);
+ cooling_update(/*phys_const=*/NULL, cosmo, /*pfloor=*/NULL, cooling,
+ /*space=*/NULL, /*time=*/0);
}
/**
diff --git a/src/cooling/EAGLE/cooling.h b/src/cooling/EAGLE/cooling.h
index ce590c2cf72399e86dd73dc19d145f2a38f9def4..a44b219d9e65bc61c184f251d3b4ba2d34181402 100644
--- a/src/cooling/EAGLE/cooling.h
+++ b/src/cooling/EAGLE/cooling.h
@@ -37,9 +37,11 @@ struct pressure_floor_props;
struct space;
struct phys_const;
-void cooling_update(const struct cosmology *cosmo,
+void cooling_update(const struct phys_const *phys_const,
+ const struct cosmology *cosmo,
const struct pressure_floor_props *pressure_floor,
- struct cooling_function_data *cooling, struct space *s);
+ struct cooling_function_data *cooling, struct space *s,
+ const double time);
void cooling_cool_part(const struct phys_const *phys_const,
const struct unit_system *us,
@@ -79,6 +81,27 @@ void cooling_first_init_part(
const struct cooling_function_data *restrict cooling,
const struct part *restrict p, struct xpart *restrict xp);
+/**
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state. The function requires the density to be defined and thus must
+ * be called after its computation.
+ *
+ * @param phys_const The #phys_const.
+ * @param us The #unit_system.
+ * @param hydro_props The #hydro_props.
+ * @param cosmo The #cosmology.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+void cooling_post_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_from_gas(
const struct phys_const *phys_const, const struct cosmology *cosmo,
const struct cooling_function_data *cooling, const float rho_phys,
diff --git a/src/cooling/PS2020/cooling.c b/src/cooling/PS2020/cooling.c
index 3292841019f4328f7c767cfde1a1bca5fa3db716..fb90cbf687ca3657c122a5572bf65fa5f80035d7 100644
--- a/src/cooling/PS2020/cooling.c
+++ b/src/cooling/PS2020/cooling.c
@@ -72,10 +72,13 @@ static const double bracket_factor = 1.5;
* @param pressure_floor Properties of the pressure floor.
* @param cooling The #cooling_function_data used in the run.
* @param s The space data, including a pointer to array of particles
+ * @param time The current system time
*/
-void cooling_update(const struct cosmology *cosmo,
+void cooling_update(const struct phys_const *phys_const,
+ const struct cosmology *cosmo,
const struct pressure_floor_props *pressure_floor,
- struct cooling_function_data *cooling, struct space *s) {
+ struct cooling_function_data *cooling, struct space *s,
+ const double time) {
/* Extra energy for reionization? */
if (!cooling->H_reion_done) {
@@ -950,6 +953,26 @@ __attribute__((always_inline)) INLINE void cooling_first_init_part(
p->cooling_data.subgrid_dens = -1.f;
}
+/**
+ * @brief Perform additional init on the cooling properties of the
+ * (x-)particles that requires the density to be known.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cosmo The current cosmological model.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE void cooling_post_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) {}
+
/**
* @brief Compute the fraction of Hydrogen that is in HI based
* on the pressure of the gas.
@@ -1583,7 +1606,8 @@ void cooling_restore_tables(struct cooling_function_data *cooling,
read_cooling_header(cooling);
read_cooling_tables(cooling);
- cooling_update(cosmo, /*pfloor=*/NULL, cooling, /*space=*/NULL);
+ cooling_update(/*phys_const=*/NULL, cosmo, /*pfloor=*/NULL, cooling,
+ /*space=*/NULL, /*time=*/0);
}
/**
diff --git a/src/cooling/PS2020/cooling.h b/src/cooling/PS2020/cooling.h
index 1bacfaa86a2d81e038c4cf71fa03ee05638e4fa2..3af8492898a1b790fa841def1e0ea8aa53838d41 100644
--- a/src/cooling/PS2020/cooling.h
+++ b/src/cooling/PS2020/cooling.h
@@ -37,9 +37,11 @@ struct pressure_floor_props;
struct feedback_props;
struct space;
-void cooling_update(const struct cosmology *cosmo,
+void cooling_update(const struct phys_const *phys_const,
+ const struct cosmology *cosmo,
const struct pressure_floor_props *pressure_floor,
- struct cooling_function_data *cooling, struct space *s);
+ struct cooling_function_data *cooling, struct space *s,
+ const double time);
void cooling_cool_part(const struct phys_const *phys_const,
const struct unit_system *us,
@@ -65,6 +67,26 @@ void cooling_first_init_part(const struct phys_const *phys_const,
const struct cooling_function_data *cooling,
struct part *p, struct xpart *xp);
+/**
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state. The function requires the density to be defined and thus must
+ * be called after its computation.
+ *
+ * @param phys_const The #phys_const.
+ * @param us The #unit_system.
+ * @param hydro_props The #hydro_props.
+ * @param cosmo The #cosmology.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+void cooling_post_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,
diff --git a/src/cooling/QLA/cooling.c b/src/cooling/QLA/cooling.c
index dfe6e04edf357f76e72a573ffcd89b73c983a289..79167d24f547a1b0885897395bce712a6a367a8f 100644
--- a/src/cooling/QLA/cooling.c
+++ b/src/cooling/QLA/cooling.c
@@ -71,10 +71,13 @@ static const double bracket_factor = 1.5;
* @param pressure_floor Properties of the pressure floor.
* @param cooling The #cooling_function_data used in the run.
* @param s The space data, including a pointer to array of particles
+ * @param time The current system time
*/
-void cooling_update(const struct cosmology *cosmo,
+void cooling_update(const struct phys_const *phys_const,
+ const struct cosmology *cosmo,
const struct pressure_floor_props *pressure_floor,
- struct cooling_function_data *cooling, struct space *s) {
+ struct cooling_function_data *cooling, struct space *s,
+ const double time) {
/* Extra energy for reionization? */
if (!cooling->H_reion_done) {
@@ -762,6 +765,26 @@ __attribute__((always_inline)) INLINE void cooling_first_init_part(
const struct cooling_function_data *cooling, struct part *p,
struct xpart *xp) {}
+/**
+ * @brief Perform additional init on the cooling properties of the
+ * (x-)particles that requires the density to be known.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cosmo The current cosmological model.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE void cooling_post_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) {}
+
/**
* @brief Compute the fraction of Hydrogen that is in HI based
* on the pressure of the gas.
@@ -1137,7 +1160,8 @@ void cooling_restore_tables(struct cooling_function_data *cooling,
read_cooling_header(cooling);
read_cooling_tables(cooling);
- cooling_update(cosmo, /*pfloor=*/NULL, cooling, /*space=*/NULL);
+ cooling_update(/*phys_const=*/NULL, cosmo, /*pfloor=*/NULL, cooling,
+ /*space=*/NULL, /*time=*/0);
}
/**
diff --git a/src/cooling/QLA/cooling.h b/src/cooling/QLA/cooling.h
index dca77090a4bac04bcd6316a9f3b8ca8c0e74af05..8ef8b792404529acf46f352ae05cbaa61e2ba366 100644
--- a/src/cooling/QLA/cooling.h
+++ b/src/cooling/QLA/cooling.h
@@ -38,9 +38,11 @@ struct pressure_floor_props;
struct phys_const;
struct space;
-void cooling_update(const struct cosmology *cosmo,
+void cooling_update(const struct phys_const *phys_const,
+ const struct cosmology *cosmo,
const struct pressure_floor_props *pressure_floor,
- struct cooling_function_data *cooling, struct space *s);
+ struct cooling_function_data *cooling, struct space *s,
+ const double time);
void cooling_cool_part(const struct phys_const *phys_const,
const struct unit_system *us,
@@ -66,6 +68,26 @@ void cooling_first_init_part(const struct phys_const *phys_const,
const struct cooling_function_data *cooling,
struct part *p, struct xpart *xp);
+/**
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state. The function requires the density to be defined and thus must
+ * be called after its computation.
+ *
+ * @param phys_const The #phys_const.
+ * @param us The #unit_system.
+ * @param hydro_props The #hydro_props.
+ * @param cosmo The #cosmology.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+void cooling_post_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,
diff --git a/src/cooling/QLA_EAGLE/cooling.c b/src/cooling/QLA_EAGLE/cooling.c
index 8e7bb884c912ad24466f8241392d236d2c6740f5..811db57538d08d9ee6ca4c085ad57af0b46dbd65 100644
--- a/src/cooling/QLA_EAGLE/cooling.c
+++ b/src/cooling/QLA_EAGLE/cooling.c
@@ -126,10 +126,13 @@ __attribute__((always_inline)) INLINE void get_redshift_index(
* @param pressure_floor The properties of the pressure floor.
* @param cooling The #cooling_function_data used in the run.
* @param s The space data, including a pointer to array of particles
+ * @param time The current system time
*/
-void cooling_update(const struct cosmology *cosmo,
+void cooling_update(const struct phys_const *phys_const,
+ const struct cosmology *cosmo,
const struct pressure_floor_props *pressure_floor,
- struct cooling_function_data *cooling, struct space *s) {
+ struct cooling_function_data *cooling, struct space *s,
+ const double time) {
/* Current redshift */
const float redshift = cosmo->z;
@@ -568,6 +571,28 @@ __attribute__((always_inline)) INLINE void cooling_first_init_part(
const struct cooling_function_data *restrict cooling,
const struct part *restrict p, struct xpart *restrict xp) {}
+/**
+ * @brief Perform additional init on the cooling properties of the
+ * (x-)particles that requires the density to be known.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cosmo The current cosmological model.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE void cooling_post_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) {}
+
/**
* @brief Compute the temperature based on gas properties.
*
@@ -1068,7 +1093,8 @@ void cooling_restore_tables(struct cooling_function_data *cooling,
/* Force a re-read of the cooling tables */
cooling->z_index = -10;
cooling->previous_z_index = qla_eagle_cooling_N_redshifts - 2;
- cooling_update(cosmo, /*pfloor=*/NULL, cooling, /*space=*/NULL);
+ cooling_update(/*phys_const=*/NULL, cosmo, /*pfloor=*/NULL, cooling,
+ /*space=*/NULL, /*time=*/0);
}
/**
diff --git a/src/cooling/QLA_EAGLE/cooling.h b/src/cooling/QLA_EAGLE/cooling.h
index 0abafc6323be12715c1b5f2de453c99b6dccf58a..9ddd756a69174694454b0ab5102c275791120912 100644
--- a/src/cooling/QLA_EAGLE/cooling.h
+++ b/src/cooling/QLA_EAGLE/cooling.h
@@ -37,9 +37,11 @@ struct pressure_floor_props;
struct space;
struct phys_const;
-void cooling_update(const struct cosmology *cosmo,
+void cooling_update(const struct phys_const *phys_const,
+ const struct cosmology *cosmo,
const struct pressure_floor_props *pressure_floor,
- struct cooling_function_data *cooling, struct space *s);
+ struct cooling_function_data *cooling, struct space *s,
+ const double time);
void cooling_cool_part(const struct phys_const *phys_const,
const struct unit_system *us,
@@ -88,6 +90,27 @@ void cooling_first_init_part(
const struct cooling_function_data *restrict cooling,
const struct part *restrict p, struct xpart *restrict xp);
+/**
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state. The function requires the density to be defined and thus must
+ * be called after its computation.
+ *
+ * @param phys_const The #phys_const.
+ * @param us The #unit_system.
+ * @param hydro_props The #hydro_props.
+ * @param cosmo The #cosmology.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+void cooling_post_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_from_gas(
const struct phys_const *phys_const, const struct cosmology *cosmo,
const struct cooling_function_data *cooling, const float rho_phys,
diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h
index 22c90bc679437cee9e7057263b930bdc4b5a6286..0ab86457c639fcdfd94361fc2984699d94a3bcd7 100644
--- a/src/cooling/const_du/cooling.h
+++ b/src/cooling/const_du/cooling.h
@@ -58,11 +58,12 @@
* @param pressure_floor The properties of the pressure floor.
* @param cooling The #cooling_function_data used in the run.
* @param s The #space containing all the particles.
+ * @param time The current system time
*/
INLINE static void cooling_update(
- const struct cosmology* cosmo,
+ const struct phys_const* phys_const, const struct cosmology* cosmo,
const struct pressure_floor_props* pressure_floor,
- struct cooling_function_data* cooling, struct space* s) {
+ struct cooling_function_data* cooling, struct space* s, const double time) {
// Add content if required.
}
@@ -226,6 +227,28 @@ __attribute__((always_inline)) INLINE static void cooling_first_init_part(
xp->cooling_data.radiated_energy = 0.f;
}
+/**
+ * @brief Perform additional init on the cooling properties of the
+ * (x-)particles that requires the density to be known.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cosmo The current cosmological model.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void cooling_post_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* cooling, const struct part* restrict p,
+ struct xpart* restrict xp) {}
+
/**
* @brief Compute the temperature of a #part based on the cooling function.
*
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index c204996e40d69f8837d7fbeddb4dbc125e98718f..bc9d3e44fede03864b90184ee9aa577da3cfd45f 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -54,11 +54,12 @@
* @param pressure_floor The properties of the pressure floor.
* @param cooling The #cooling_function_data used in the run.
* @param s The #space containing all the particles.
+ * @param time The current system time
*/
INLINE static void cooling_update(
- const struct cosmology* cosmo,
+ const struct phys_const* phys_const, const struct cosmology* cosmo,
const struct pressure_floor_props* pressure_floor,
- struct cooling_function_data* cooling, struct space* s) {
+ struct cooling_function_data* cooling, struct space* s, const double time) {
// Add content if required.
}
@@ -314,6 +315,28 @@ __attribute__((always_inline)) INLINE static void cooling_first_init_part(
xp->cooling_data.radiated_energy = 0.f;
}
+/**
+ * @brief Perform additional init on the cooling properties of the
+ * (x-)particles that requires the density to be known.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cosmo The current cosmological model.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void cooling_post_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* cooling, const struct part* restrict p,
+ struct xpart* restrict xp) {}
+
/**
* @brief Compute the temperature of a #part based on the cooling function.
*
diff --git a/src/cooling/grackle/cooling.c b/src/cooling/grackle/cooling.c
index 2c4e093c9b8f986c1a7abec36aa68d17eb8cadc6..215ec4cd9db8570dc1603a9f5f46ca8bb3afc473 100644
--- a/src/cooling/grackle/cooling.c
+++ b/src/cooling/grackle/cooling.c
@@ -68,14 +68,18 @@ gr_float cooling_time(const struct phys_const* phys_const,
* @brief Common operations performed on the cooling function at a
* given time-step or redshift.
*
+ * @param phys_const The #phys_const.
* @param cosmo The current cosmological model.
* @param pressure_floor Properties of the pressure floor.
* @param cooling The #cooling_function_data used in the run.
* @param s The #space containing all the particles.
+ * @param time The current system time
*/
-void cooling_update(const struct cosmology* cosmo,
+void cooling_update(const struct phys_const* phys_const,
+ const struct cosmology* cosmo,
const struct pressure_floor_props* pressure_floor,
- struct cooling_function_data* cooling, struct space* s) {
+ struct cooling_function_data* cooling, struct space* s,
+ const double time) {
/* set current time */
if (cooling->redshift == -1)
cooling->units.a_value = cosmo->a;
@@ -180,7 +184,7 @@ void cooling_compute_equilibrium(const struct phys_const* phys_const,
/* get temporary data */
struct part p_tmp = *p;
struct cooling_function_data cooling_tmp = *cooling;
- cooling_tmp.chemistry.with_radiative_cooling = 0;
+ cooling_tmp.chemistry_data.with_radiative_cooling = 0;
/* need density for computation, therefore quick estimate */
p_tmp.rho = 0.2387 * p_tmp.mass / pow(p_tmp.h, 3);
@@ -241,15 +245,19 @@ void cooling_first_init_part(const struct phys_const* phys_const,
#if COOLING_GRACKLE_MODE >= 1
gr_float zero = 1.e-20;
+ /* NOTE: at this stage, we assume neutral gas
+ * a better determination will be done in cooling_post_init_part */
+
/* primordial chemistry >= 1 */
- xp->cooling_data.HI_frac = zero;
- xp->cooling_data.HII_frac = grackle_data->HydrogenFractionByMass;
+ xp->cooling_data.HI_frac = grackle_data->HydrogenFractionByMass;
+ xp->cooling_data.HII_frac = zero;
xp->cooling_data.HeI_frac = zero;
xp->cooling_data.HeII_frac = zero;
xp->cooling_data.HeIII_frac = 1. - grackle_data->HydrogenFractionByMass;
xp->cooling_data.e_frac = xp->cooling_data.HII_frac +
0.25 * xp->cooling_data.HeII_frac +
0.5 * xp->cooling_data.HeIII_frac;
+
#endif // MODE >= 1
#if COOLING_GRACKLE_MODE >= 2
@@ -266,11 +274,36 @@ void cooling_first_init_part(const struct phys_const* phys_const,
xp->cooling_data.DII_frac = zero;
xp->cooling_data.HDI_frac = zero;
#endif // MODE >= 3
+}
+
+/**
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state. The function requires the density to be defined and thus must
+ * be called after its computation.
+ *
+ * @param phys_const The #phys_const.
+ * @param us The #unit_system.
+ * @param hydro_props The #hydro_props.
+ * @param cosmo The #cosmology.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+void cooling_post_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,
+ const struct part* p, struct xpart* xp) {
+
+ // const float rho = hydro_get_physical_density(p, cosmo);
+ // const float energy = hydro_get_physical_internal_energy(p, xp, cosmo);
+ // message("rho = %g energy = %g",rho,energy);
#if COOLING_GRACKLE_MODE > 0
/* TODO: this can fail spectacularly and needs to be replaced. */
- cooling_compute_equilibrium(phys_const, us, hydro_props, cosmo, cooling, p,
- xp);
+ // cooling_compute_equilibrium(phys_const, us, hydro_props, cosmo, cooling,
+ // p,xp);
#endif
}
@@ -292,21 +325,32 @@ float cooling_get_radiated_energy(const struct xpart* xp) {
void cooling_print_backend(const struct cooling_function_data* cooling) {
message("Cooling function is 'Grackle'.");
- message("Using Grackle = %i", cooling->chemistry.use_grackle);
- message("Chemical network = %i", cooling->chemistry.primordial_chemistry);
+ message("Using Grackle = %i", cooling->chemistry_data.use_grackle);
+ message("Chemical network = %i",
+ cooling->chemistry_data.primordial_chemistry);
message("CloudyTable = %s", cooling->cloudy_table);
message("Redshift = %g", cooling->redshift);
message("UV background = %d", cooling->with_uv_background);
- message("Metal cooling = %i", cooling->chemistry.metal_cooling);
+ message("Metal cooling = %i", cooling->chemistry_data.metal_cooling);
message("Self Shielding = %i", cooling->self_shielding_method);
if (cooling->self_shielding_method == -1) {
message("Self Shelding density = %g", cooling->self_shielding_threshold);
}
message("Thermal time = %g", cooling->thermal_time);
- message("Specific Heating Rates = %i",
- cooling->provide_specific_heating_rates);
- message("Volumetric Heating Rates = %i",
- cooling->provide_volumetric_heating_rates);
+ message("Specific Heating Rates = %g", cooling->specific_heating_rates);
+ message("Volumetric Heating Rates = %g", cooling->volumetric_heating_rates);
+
+ message("grackle_chemistry_data.RT_heating_rate = %g",
+ cooling->RT_heating_rate);
+ message("grackle_chemistry_data.RT_HI_ionization_rate = %g",
+ cooling->RT_HI_ionization_rate);
+ message("grackle_chemistry_data.RT_HeI_ionization_rate = %g",
+ cooling->RT_HeI_ionization_rate);
+ message("grackle_chemistry_data.RT_HeII_ionization_rate = %g",
+ cooling->RT_HeII_ionization_rate);
+ message("grackle_chemistry_data.RT_H2_dissociation_rate = %g",
+ cooling->RT_H2_dissociation_rate);
+
message("Units:");
message("\tComoving = %i", cooling->units.comoving_coordinates);
message("\tLength = %g", cooling->units.length_units);
@@ -314,6 +358,35 @@ void cooling_print_backend(const struct cooling_function_data* cooling) {
message("\tTime = %g", cooling->units.time_units);
message("\tScale Factor = %g (units: %g)", cooling->units.a_value,
cooling->units.a_units);
+
+ message("Grackle parameters:");
+ message("grackle_chemistry_data.use_grackle = %d",
+ cooling->chemistry_data.use_grackle);
+ message("grackle_chemistry_data.with_radiative_cooling %d",
+ cooling->chemistry_data.with_radiative_cooling);
+ message("grackle_chemistry_data.primordial_chemistry = %d",
+ cooling->chemistry_data.primordial_chemistry);
+ message("grackle_chemistry_data.dust_chemistry = %d",
+ cooling->chemistry_data.dust_chemistry);
+ message("grackle_chemistry_data.metal_cooling = %d",
+ cooling->chemistry_data.metal_cooling);
+ message("grackle_chemistry_data.UVbackground = %d",
+ cooling->chemistry_data.UVbackground);
+ message("grackle_chemistry_data.CaseBRecombination = %d",
+ cooling->chemistry_data.CaseBRecombination);
+ message("grackle_chemistry_data.grackle_data_file = %s",
+ cooling->chemistry_data.grackle_data_file);
+ message("grackle_chemistry_data.use_radiative_transfer = %d",
+ cooling->chemistry_data.use_radiative_transfer);
+ message("grackle_chemistry_data.use_volumetric_heating_rate = %d",
+ cooling->chemistry_data.use_volumetric_heating_rate);
+ message("grackle_chemistry_data.use_specific_heating_rate = %d",
+ cooling->chemistry_data.use_specific_heating_rate);
+ message("grackle_chemistry_data.self_shielding_method = %d",
+ cooling->chemistry_data.self_shielding_method);
+ message("grackle_chemistry_data.HydrogenFractionByMass = %.3g",
+ cooling->chemistry_data.HydrogenFractionByMass);
+ message("grackle_chemistry_data.Gamma = %.6g", cooling->chemistry_data.Gamma);
}
/**
@@ -529,19 +602,84 @@ void cooling_copy_from_grackle3(grackle_field_data* data, const struct part* p,
*/
void cooling_copy_to_grackle(grackle_field_data* data, const struct part* p,
struct xpart* xp, gr_float rho,
- gr_float species_densities[12]) {
+ gr_float species_densities[12],
+ const struct cooling_function_data* cooling,
+ const struct phys_const* phys_const) {
+
+ const float time_units = cooling->units.time_units;
cooling_copy_to_grackle1(data, p, xp, rho, species_densities);
cooling_copy_to_grackle2(data, p, xp, rho, species_densities);
cooling_copy_to_grackle3(data, p, xp, rho, species_densities);
- data->volumetric_heating_rate = NULL;
- data->specific_heating_rate = NULL;
- data->RT_heating_rate = NULL;
- data->RT_HI_ionization_rate = NULL;
- data->RT_HeI_ionization_rate = NULL;
- data->RT_HeII_ionization_rate = NULL;
- data->RT_H2_dissociation_rate = NULL;
+ if (cooling->chemistry_data.use_volumetric_heating_rate) {
+ gr_float* volumetric_heating_rate = (gr_float*)malloc(sizeof(gr_float));
+ *volumetric_heating_rate = cooling->volumetric_heating_rates;
+ data->volumetric_heating_rate = volumetric_heating_rate;
+ }
+
+ if (cooling->chemistry_data.use_specific_heating_rate) {
+ gr_float* specific_heating_rate = (gr_float*)malloc(sizeof(gr_float));
+ *specific_heating_rate = cooling->specific_heating_rates;
+ data->specific_heating_rate = specific_heating_rate;
+ }
+
+ if (cooling->chemistry_data.use_radiative_transfer) {
+
+ /* heating rate */
+ gr_float* RT_heating_rate = (gr_float*)malloc(sizeof(gr_float));
+ *RT_heating_rate = cooling->RT_heating_rate;
+ /* Note to self:
+ * If cooling->RT_heating_rate is computed properly, i.e. using
+ * the HI density, and then being HI density dependent, we need
+ * to divide it as follow. If it is assumed to be already normed
+ * as it is so when providing it via some parameters, we keep it
+ * unchanged.
+ */
+ /* Grackle wants heating rate in units of / nHI_cgs */
+ // const double nHI_cgs = species_densities[0]
+ // / phys_const->const_proton_mass
+ // / pow(length_units,3);
+ //*RT_heating_rate /= nHI_cgs;
+ data->RT_heating_rate = RT_heating_rate;
+
+ /* HI ionization rate */
+ gr_float* RT_HI_ionization_rate = (gr_float*)malloc(sizeof(gr_float));
+ *RT_HI_ionization_rate = cooling->RT_HI_ionization_rate;
+ /* Grackle wants it in 1/internal_time_units */
+ *RT_HI_ionization_rate /= (1. / time_units);
+ data->RT_HI_ionization_rate = RT_HI_ionization_rate;
+
+ /* HeI ionization rate */
+ gr_float* RT_HeI_ionization_rate = (gr_float*)malloc(sizeof(gr_float));
+ *RT_HeI_ionization_rate = cooling->RT_HeI_ionization_rate;
+ /* Grackle wants it in 1/internal_time_units */
+ *RT_HeI_ionization_rate /= (1. / time_units);
+ data->RT_HeI_ionization_rate = RT_HeI_ionization_rate;
+
+ /* HeII ionization rate */
+ gr_float* RT_HeII_ionization_rate = (gr_float*)malloc(sizeof(gr_float));
+ *RT_HeII_ionization_rate = cooling->RT_HeII_ionization_rate;
+ /* Grackle wants it in 1/internal_time_units */
+ *RT_HeII_ionization_rate /= (1. / time_units);
+ data->RT_HeII_ionization_rate = RT_HeII_ionization_rate;
+
+ /* H2 ionization rate */
+ gr_float* RT_H2_dissociation_rate = (gr_float*)malloc(sizeof(gr_float));
+ *RT_H2_dissociation_rate = cooling->RT_H2_dissociation_rate;
+ /* Grackle wants it in 1/internal_time_units */
+ *RT_H2_dissociation_rate /= (1. / time_units);
+ data->RT_H2_dissociation_rate = RT_H2_dissociation_rate;
+
+ } else {
+ data->volumetric_heating_rate = NULL;
+ data->specific_heating_rate = NULL;
+ data->RT_heating_rate = NULL;
+ data->RT_HI_ionization_rate = NULL;
+ data->RT_HeI_ionization_rate = NULL;
+ data->RT_HeII_ionization_rate = NULL;
+ data->RT_H2_dissociation_rate = NULL;
+ }
gr_float* metal_density = (gr_float*)malloc(sizeof(gr_float));
*metal_density = chemistry_get_total_metal_mass_fraction_for_cooling(p) * rho;
@@ -560,11 +698,25 @@ void cooling_copy_to_grackle(grackle_field_data* data, const struct part* p,
* @param rho The particle density.
*/
void cooling_copy_from_grackle(grackle_field_data* data, const struct part* p,
- struct xpart* xp, gr_float rho) {
+ struct xpart* xp, gr_float rho,
+ const struct cooling_function_data* cooling) {
cooling_copy_from_grackle1(data, p, xp, rho);
cooling_copy_from_grackle2(data, p, xp, rho);
cooling_copy_from_grackle3(data, p, xp, rho);
+ if (cooling->chemistry_data.use_volumetric_heating_rate)
+ free(data->volumetric_heating_rate);
+ if (cooling->chemistry_data.use_specific_heating_rate)
+ free(data->specific_heating_rate);
+
+ if (cooling->chemistry_data.use_radiative_transfer) {
+ free(data->RT_heating_rate);
+ free(data->RT_HI_ionization_rate);
+ free(data->RT_HeI_ionization_rate);
+ free(data->RT_HeII_ionization_rate);
+ free(data->RT_H2_dissociation_rate);
+ }
+
free(data->metal_density);
}
@@ -622,7 +774,8 @@ gr_float cooling_new_energy(const struct phys_const* phys_const,
/* set current time */
code_units units = cooling->units;
- chemistry_data chemistry_grackle = cooling->chemistry;
+ chemistry_data chemistry_grackle = cooling->chemistry_data;
+ chemistry_data_storage rates_grackle = cooling->chemistry_rates;
/* initialize data */
grackle_field_data data;
@@ -660,19 +813,20 @@ gr_float cooling_new_energy(const struct phys_const* phys_const,
data.z_velocity = NULL;
/* copy to grackle structure */
- cooling_copy_to_grackle(&data, p, xp, density, species_densities);
+ cooling_copy_to_grackle(&data, p, xp, density, species_densities, cooling,
+ phys_const);
/* Apply the self shielding if requested */
cooling_apply_self_shielding(cooling, &chemistry_grackle, p, cosmo);
/* solve chemistry */
- if (local_solve_chemistry(&chemistry_grackle, &grackle_rates, &units, &data,
+ if (local_solve_chemistry(&chemistry_grackle, &rates_grackle, &units, &data,
dt) == 0) {
error("Error in solve_chemistry.");
}
/* copy from grackle data to particle */
- cooling_copy_from_grackle(&data, p, xp, density);
+ cooling_copy_from_grackle(&data, p, xp, density, cooling);
return energy;
}
@@ -702,7 +856,8 @@ gr_float cooling_time(const struct phys_const* phys_const,
/* initialize data */
grackle_field_data data;
- chemistry_data chemistry_grackle = cooling->chemistry;
+ chemistry_data chemistry_grackle = cooling->chemistry_data;
+ chemistry_data_storage rates_grackle = cooling->chemistry_rates;
/* set values */
/* grid */
@@ -733,20 +888,21 @@ gr_float cooling_time(const struct phys_const* phys_const,
gr_float species_densities[12];
/* copy data from particle to grackle data */
- cooling_copy_to_grackle(&data, p, xp, density, species_densities);
+ cooling_copy_to_grackle(&data, p, xp, density, species_densities, cooling,
+ phys_const);
/* Apply the self shielding if requested */
cooling_apply_self_shielding(cooling, &chemistry_grackle, p, cosmo);
/* Compute cooling time */
gr_float cooling_time;
- if (local_calculate_cooling_time(&chemistry_grackle, &grackle_rates, &units,
+ if (local_calculate_cooling_time(&chemistry_grackle, &rates_grackle, &units,
&data, &cooling_time) == 0) {
error("Error in calculate_cooling_time.");
}
/* copy from grackle data to particle */
- cooling_copy_from_grackle(&data, p, xp, density);
+ cooling_copy_from_grackle(&data, p, xp, density, cooling);
/* compute rate */
return cooling_time;
@@ -977,7 +1133,7 @@ void cooling_init_grackle(struct cooling_function_data* cooling) {
grackle_verbose = 1;
#endif
- chemistry_data* chemistry = &cooling->chemistry;
+ chemistry_data* chemistry = &cooling->chemistry_data;
/* Create a chemistry object for parameters and rate data. */
if (set_default_chemistry_parameters(chemistry) == 0) {
@@ -996,25 +1152,33 @@ void cooling_init_grackle(struct cooling_function_data* cooling) {
chemistry->grackle_data_file = cooling->cloudy_table;
/* radiative transfer */
- chemistry->use_radiative_transfer = cooling->provide_specific_heating_rates ||
- cooling->provide_volumetric_heating_rates;
- chemistry->use_volumetric_heating_rate =
- cooling->provide_volumetric_heating_rates;
- chemistry->use_specific_heating_rate =
- cooling->provide_specific_heating_rates;
-
- if (cooling->provide_specific_heating_rates &&
- cooling->provide_volumetric_heating_rates)
- message(
- "WARNING: You should specified either the specific or the volumetric "
+ chemistry->use_radiative_transfer = cooling->use_radiative_transfer;
+
+ if (cooling->volumetric_heating_rates > 0)
+ chemistry->use_volumetric_heating_rate = 1;
+
+ if (cooling->specific_heating_rates > 0)
+ chemistry->use_specific_heating_rate = 1;
+
+ /* hydrogen fraction by mass */
+ chemistry->HydrogenFractionByMass = cooling->HydrogenFractionByMass;
+
+ /* use the Case B recombination rates */
+ chemistry->CaseBRecombination = 1;
+
+ if (cooling->specific_heating_rates > 0 &&
+ cooling->volumetric_heating_rates > 0)
+ error(
+ "You should specified either the specific or the volumetric "
"heating rates, not both");
/* self shielding */
chemistry->self_shielding_method = cooling->self_shielding_method;
- /* Initialize the chemistry object. */
- if (initialize_chemistry_data(&cooling->units) == 0) {
- error("Error in initialize_chemistry_data.");
+ if (local_initialize_chemistry_data(&cooling->chemistry_data,
+ &cooling->chemistry_rates,
+ &cooling->units) == 0) {
+ error("Error in initialize_chemistry_data");
}
}
@@ -1052,7 +1216,9 @@ void cooling_init_backend(struct swift_params* parameter_file,
* @param cooling the cooling data structure.
*/
void cooling_clean(struct cooling_function_data* cooling) {
- _free_chemistry_data(&cooling->chemistry, &grackle_rates);
+ /* Clean up grackle data. This is a call to a grackle function */
+ local_free_chemistry_data(&cooling->chemistry_data,
+ &cooling->chemistry_rates);
}
/**
diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h
index 3584bd7f188c9a6ea5bf0501635f771753318d10..92f6c921634bb62e58626effbcdc1019b2ecc75f 100644
--- a/src/cooling/grackle/cooling.h
+++ b/src/cooling/grackle/cooling.h
@@ -44,9 +44,11 @@ struct swift_params;
#define GRACKLE_NPART 1
#define GRACKLE_RANK 3
-void cooling_update(const struct cosmology* cosmo,
+void cooling_update(const struct phys_const* phys_const,
+ const struct cosmology* cosmo,
const struct pressure_floor_props* pressure_floor,
- struct cooling_function_data* cooling, struct space* s);
+ struct cooling_function_data* cooling, struct space* s,
+ const double time);
void cooling_first_init_part(const struct phys_const* phys_const,
const struct unit_system* us,
@@ -55,6 +57,26 @@ void cooling_first_init_part(const struct phys_const* phys_const,
const struct cooling_function_data* cooling,
const struct part* p, struct xpart* xp);
+/**
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state. The function requires the density to be defined and thus must
+ * be called after its computation.
+ *
+ * @param phys_const The #phys_const.
+ * @param us The #unit_system.
+ * @param hydro_props The #hydro_props.
+ * @param cosmo The #cosmology.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+void cooling_post_init_part(const struct phys_const* phys_const,
+ const struct unit_system* us,
+ const struct hydro_props* hydro_properties,
+ const struct cosmology* cosmo,
+ const struct cooling_function_data* cooling,
+ const struct part* p, struct xpart* xp);
+
/**
* @brief Returns the subgrid temperature of a particle.
*
diff --git a/src/cooling/grackle/cooling_io.h b/src/cooling/grackle/cooling_io.h
index 290b9d894b212b2a70c6f87bc7b3347b18d77087..fb15b730f9c75260df0a56fca0d2e19c60e748df 100644
--- a/src/cooling/grackle/cooling_io.h
+++ b/src/cooling/grackle/cooling_io.h
@@ -170,11 +170,32 @@ __attribute__((always_inline)) INLINE static void cooling_read_parameters(
cooling->with_metal_cooling =
parser_get_param_int(parameter_file, "GrackleCooling:with_metal_cooling");
- cooling->provide_volumetric_heating_rates = parser_get_opt_param_int(
- parameter_file, "GrackleCooling:provide_volumetric_heating_rates", 0);
+ cooling->use_radiative_transfer = parser_get_opt_param_int(
+ parameter_file, "GrackleCooling:use_radiative_transfer", 0);
- cooling->provide_specific_heating_rates = parser_get_opt_param_int(
- parameter_file, "GrackleCooling:provide_specific_heating_rates", 0);
+ cooling->RT_heating_rate = parser_get_opt_param_double(
+ parameter_file, "GrackleCooling:RT_heating_rate_cgs", 0);
+
+ cooling->RT_HI_ionization_rate = parser_get_opt_param_double(
+ parameter_file, "GrackleCooling:RT_HI_ionization_rate_cgs", 0);
+
+ cooling->RT_HeI_ionization_rate = parser_get_opt_param_double(
+ parameter_file, "GrackleCooling:RT_HeI_ionization_rate_cgs", 0);
+
+ cooling->RT_HeII_ionization_rate = parser_get_opt_param_double(
+ parameter_file, "GrackleCooling:RT_HeII_ionization_rate_cgs", 0);
+
+ cooling->RT_H2_dissociation_rate = parser_get_opt_param_double(
+ parameter_file, "GrackleCooling:RT_H2_dissociation_rate_cgs", 0);
+
+ cooling->volumetric_heating_rates = parser_get_opt_param_double(
+ parameter_file, "GrackleCooling:volumetric_heating_rates_cgs", 0);
+
+ cooling->specific_heating_rates = parser_get_opt_param_double(
+ parameter_file, "GrackleCooling:specific_heating_rates_cgs", 0);
+
+ cooling->HydrogenFractionByMass = parser_get_opt_param_double(
+ parameter_file, "GrackleCooling:HydrogenFractionByMass", 0.76);
/* Self shielding */
cooling->self_shielding_method = parser_get_opt_param_int(
diff --git a/src/cooling/grackle/cooling_properties.h b/src/cooling/grackle/cooling_properties.h
index 385ecdf59fec48d94130426389acd732c922083c..93b602cffdb4b6130bb773ffa66b7c0c0a6dd394 100644
--- a/src/cooling/grackle/cooling_properties.h
+++ b/src/cooling/grackle/cooling_properties.h
@@ -19,6 +19,9 @@
#ifndef SWIFT_COOLING_PROPERTIES_GRACKLE_H
#define SWIFT_COOLING_PROPERTIES_GRACKLE_H
+/* skip deprecation warnings. I cleaned old API calls. */
+#define OMIT_LEGACY_INTERNAL_GRACKLE_FUNC
+
/* include grackle */
#include <grackle.h>
@@ -48,16 +51,41 @@ struct cooling_function_data {
code_units units;
/*! grackle chemistry data */
- chemistry_data chemistry;
+ chemistry_data chemistry_data;
+
+ /*! grackle chemistry data storage
+ * (needed for local function calls) */
+ chemistry_data_storage chemistry_rates;
/*! Enable/Disable metal cooling */
int with_metal_cooling;
- /*! User provide volumetric heating rates */
- int provide_volumetric_heating_rates;
+ /*! Arrays of ionization and heating rates are provided */
+ int use_radiative_transfer;
+
+ /*! Grackle RT_heating_rate (in IU) */
+ float RT_heating_rate;
+
+ /*! Grackle RT_HI_ionization_rate (in IU) */
+ float RT_HI_ionization_rate;
+
+ /*! Grackle RT_HeI_ionization_rate (in IU) */
+ float RT_HeI_ionization_rate;
+
+ /*! Grackle RT_HeII_ionization_rate (in IU) */
+ float RT_HeII_ionization_rate;
+
+ /*! Grackle RT_H2_dissociation_rate (in IU) */
+ float RT_H2_dissociation_rate;
+
+ /*! Volumetric heating rates */
+ float volumetric_heating_rates;
+
+ /*! Specific heating rates */
+ float specific_heating_rates;
- /*! User provide specific heating rates */
- int provide_specific_heating_rates;
+ /*! Hydrogen fraction by mass */
+ float HydrogenFractionByMass;
/*! Self shielding method (1 -> 3 for grackle's ones, 0 for none and -1 for
* GEAR) */
diff --git a/src/cooling/none/cooling.h b/src/cooling/none/cooling.h
index 404f5ff48497f4461c1a1ab0cab8a5a7b163c9c8..dbb5e9872bf4e197f372c3f4a121c1beaaa775f9 100644
--- a/src/cooling/none/cooling.h
+++ b/src/cooling/none/cooling.h
@@ -45,11 +45,12 @@
* @param cooling The #cooling_function_data used in the run.
* @param pressure_floor Properties of the pressure floor.
* @param s The #space containing all the particles.
+ * @param time The current system time
*/
INLINE static void cooling_update(
- const struct cosmology* cosmo,
+ const struct phys_const* phys_const, const struct cosmology* cosmo,
const struct pressure_floor_props* pressure_floor,
- struct cooling_function_data* cooling, struct space* s) {
+ struct cooling_function_data* cooling, struct space* s, const double time) {
// Add content if required.
}
@@ -127,6 +128,28 @@ __attribute__((always_inline)) INLINE static void cooling_first_init_part(
const struct cooling_function_data* data, const struct part* restrict p,
struct xpart* restrict xp) {}
+/**
+ * @brief Perform additional init on the cooling properties of the
+ * (x-)particles that requires the density to be known.
+ *
+ * Nothing to do here.
+ *
+ * @param phys_const The physical constant in internal units.
+ * @param us The unit system.
+ * @param hydro_props The properties of the hydro scheme.
+ * @param cosmo The current cosmological model.
+ * @param cooling The properties of the cooling function.
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void cooling_post_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* cooling, const struct part* restrict p,
+ struct xpart* restrict xp) {}
+
/**
* @brief Compute the temperature of a #part based on the cooling function.
*
diff --git a/src/engine.c b/src/engine.c
index a4eabaade5f5edc7ad8f8b082686821fb1c65586..0873b8332277e1e64c19788f62ce1ddf44f8cab2 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1953,8 +1953,8 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
/* Update the cooling function */
if ((e->policy & engine_policy_cooling) ||
(e->policy & engine_policy_temperature))
- cooling_update(e->cosmology, e->pressure_floor_props, e->cooling_func,
- e->s);
+ cooling_update(e->physical_constants, e->cosmology, e->pressure_floor_props,
+ e->cooling_func, e->s, e->time);
#ifdef WITH_CSDS
if (e->policy & engine_policy_csds) {
@@ -2003,6 +2003,9 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
}
}
+ /* Do some post initialisations */
+ space_post_init_parts(e->s, e->verbose);
+
/* Apply some RT conversions (e.g. energy -> energy density) */
if (e->policy & engine_policy_rt)
space_convert_rt_quantities(e->s, e->verbose);
@@ -2374,8 +2377,8 @@ int engine_step(struct engine *e) {
/* Update the cooling function */
if ((e->policy & engine_policy_cooling) ||
(e->policy & engine_policy_temperature))
- cooling_update(e->cosmology, e->pressure_floor_props, e->cooling_func,
- e->s);
+ cooling_update(e->physical_constants, e->cosmology, e->pressure_floor_props,
+ e->cooling_func, e->s, e->time);
/* Update the softening lengths */
if (e->policy & engine_policy_self_gravity)
diff --git a/src/space.c b/src/space.c
index 3402f8f3827a88a4c65aecda27ea096e1d8eb5a7..d63e8a61398ee3cb87cc35f57441f31faff9c9a0 100644
--- a/src/space.c
+++ b/src/space.c
@@ -844,6 +844,58 @@ void space_convert_rt_quantities(struct space *s, int verbose) {
clocks_getunit());
}
+void space_post_init_parts_mapper(void *restrict map_data, int count,
+ void *restrict extra_data) {
+ struct space *s = (struct space *)extra_data;
+ const struct engine *restrict e = s->e;
+
+ const struct hydro_props *restrict hydro_props = e->hydro_properties;
+ const struct phys_const *restrict phys_const = e->physical_constants;
+ const struct unit_system *us = s->e->internal_units;
+ const struct cosmology *restrict cosmo = e->cosmology;
+ const struct cooling_function_data *cool_func = e->cooling_func;
+
+ struct part *restrict p = (struct part *)map_data;
+ const ptrdiff_t delta = p - s->parts;
+ struct xpart *restrict xp = s->xparts + delta;
+
+ /* Loop over all the particles ignoring the extra buffer ones for on-the-fly
+ * creation
+ * Here we can initialize the cooling properties of the (x-)particles
+ * using quantities (like the density) defined only after the neighbour loop.
+ *
+ * */
+
+ for (int k = 0; k < count; k++) {
+ cooling_post_init_part(phys_const, us, hydro_props, cosmo, cool_func, &p[k],
+ &xp[k]);
+ }
+}
+
+/**
+ * @brief Calls the #part post-initialisation function on all particles in the
+ * space.
+ * Here we can initialize the cooling properties of the (x-)particles
+ * using quantities (like the density) defined only after the initial neighbour
+ * loop.
+ *
+ * @param s The #space.
+ * @param verbose Are we talkative?
+ */
+void space_post_init_parts(struct space *s, int verbose) {
+
+ const ticks tic = getticks();
+
+ if (s->nr_parts > 0)
+ threadpool_map(&s->e->threadpool, space_post_init_parts_mapper, s->parts,
+ s->nr_parts, sizeof(struct part), threadpool_auto_chunk_size,
+ s);
+
+ if (verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
+
void space_collect_sum_part_mass(void *restrict map_data, int count,
void *restrict extra_data) {
diff --git a/src/space.h b/src/space.h
index d98fd0f211ece1fe94385b4c1483329fb13eed90..946bf0a8597482ad633a0b8350868c521ba12572 100644
--- a/src/space.h
+++ b/src/space.h
@@ -428,6 +428,7 @@ void space_init_sinks(struct space *s, int verbose);
void space_after_snap_tracer(struct space *s, int verbose);
void space_convert_quantities(struct space *s, int verbose);
void space_convert_rt_quantities(struct space *s, int verbose);
+void space_post_init_parts(struct space *s, int verbose);
void space_link_cleanup(struct space *s);
void space_check_drift_point(struct space *s, integertime_t ti_drift,
int multipole);