diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index 4703c7091550c8c496952d9e96b623e180c78a69..7226f8838542abf580d5ab26aada7420179ec0d5 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -760,8 +760,10 @@ WARN_LOGFILE =
# Note: If this tag is empty the current directory is searched.
INPUT = @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
-INPUT += @top_srcdir@/src/hydro/Minimal @top_srcdir@/src/gravity/Default
-INPUT += @top_srcdir@/src/riemann
+INPUT += @top_srcdir@/src/hydro/Minimal
+INPUT += @top_srcdir@/src/gravity/Default
+INPUT += @top_srcdir@/src/riemann
+INPUT += @top_srcdir@/src/cooling/const
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
diff --git a/src/Makefile.am b/src/Makefile.am
index a69d59ba4f1fbd541c335e9f208609582c50251b..df1491f32e70bf298ea2b18539f6f16e16ed6871 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -51,7 +51,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
units.c common_io.c single_io.c multipole.c version.c map.c \
kernel_hydro.c tools.c part.c partition.c clocks.c parser.c \
physical_constants.c potentials.c hydro_properties.c \
- runner_doiact_fft.c threadpool.c cooling.c
+ runner_doiact_fft.c threadpool.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
@@ -71,7 +71,8 @@ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_h
hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h hydro/Gizmo/hydro_io.h \
hydro/Gizmo/hydro_debug.h hydro/Gizmo/hydro_part.h \
riemann.h \
- riemann/riemann_hllc.h riemann/riemann_trrs.h riemann/riemann_exact.h
+ riemann/riemann_hllc.h riemann/riemann_trrs.h riemann/riemann_exact.h \
+ cooling/const/cooling.h
# Sources and flags for regular library
libswiftsim_la_SOURCES = $(AM_SOURCES)
diff --git a/src/const.h b/src/const.h
index 3c2ef65219d57559284ba11b189daef18de0958f..e83d71d8733bdc43d1d8735de0d3cf0112019a61 100644
--- a/src/const.h
+++ b/src/const.h
@@ -95,8 +95,9 @@
//#define EXTERNAL_POTENTIAL_DISK_PATCH
/* Cooling properties */
-//#define CONST_COOLING
-#define CREASEY_COOLING
+#define COOLING_CONST_COOLING
+//#define COOLING_CREASEY_COOLING
+//#define COOLING_GRACKLE_COOLING
/* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS
diff --git a/src/cooling.c b/src/cooling.c
deleted file mode 100644
index 4ee03aa9bcbe0fb616fe90c1844da98a3c1ccf79..0000000000000000000000000000000000000000
--- a/src/cooling.c
+++ /dev/null
@@ -1,175 +0,0 @@
-
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- * Richard Bower (r.g.bower@durham.ac.uk)
- * Stefan Arridge (stefan.arridge@durham.ac.uk)
- *
- * 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/>.
- *
- ******************************************************************************/
-
-/* Config parameters. */
-#include "../config.h"
-
-/* This object's header. */
-#include "adiabatic_index.h"
-#include "cooling.h"
-#include "hydro.h"
-
-/**
- * @brief Initialises the cooling properties in the internal system
- * of units.
- *
- * @param parameter_file The parsed parameter file
- * @param us The current internal system of units
- * @param cooling The cooling properties to initialize
- */
-void cooling_init(const struct swift_params* parameter_file,
- struct UnitSystem* us,
- const struct phys_const* const phys_const,
- struct cooling_data* cooling) {
-
-#ifdef CONST_COOLING
- cooling->const_cooling.lambda =
- parser_get_param_double(parameter_file, "Cooling:lambda");
- cooling->const_cooling.min_energy =
- parser_get_param_double(parameter_file, "Cooling:min_energy");
- cooling->const_cooling.cooling_tstep_mult =
- parser_get_param_double(parameter_file, "Cooling:cooling_tstep_mult");
-#endif /* CONST_COOLING */
-
-#ifdef CREASEY_COOLING
- cooling->creasey_cooling.lambda =
- parser_get_param_double(parameter_file, "CreaseyCooling:Lambda");
- cooling->creasey_cooling.min_temperature = parser_get_param_double(
- parameter_file, "CreaseyCooling:minimum_temperature");
- cooling->creasey_cooling.mean_molecular_weight = parser_get_param_double(
- parameter_file, "CreaseyCooling:mean_molecular_weight");
- cooling->creasey_cooling.hydrogen_mass_abundance = parser_get_param_double(
- parameter_file, "CreaseyCooling:hydrogen_mass_abundance");
- cooling->creasey_cooling.cooling_tstep_mult = parser_get_param_double(
- parameter_file, "CreaseyCooling:cooling_tstep_mult");
-
- /*convert minimum temperature into minimum internal energy*/
- float u_floor =
- phys_const->const_boltzmann_k * cooling->creasey_cooling.min_temperature /
- (hydro_gamma_minus_one * cooling->creasey_cooling.mean_molecular_weight *
- phys_const->const_proton_mass);
- float u_floor_cgs =
- u_floor * units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
-
- cooling->creasey_cooling.min_internal_energy = u_floor;
- cooling->creasey_cooling.min_internal_energy_cgs = u_floor_cgs;
-#endif /* CREASEY_COOLING */
-}
-
-/**
- * @brief Prints the properties of the cooling model to stdout.
- *
- * @param cooling The cooling properties.
- */
-void cooling_print(const struct cooling_data* cooling) {
-
-#ifdef CONST_COOLING
- message(
- "Cooling properties are (lambda, min_energy, tstep multiplier) %g %g %g ",
- cooling->const_cooling.lambda, cooling->const_cooling.min_energy,
- cooling->const_cooling.cooling_tstep_mult);
-#endif /* CONST_COOLING */
-
-#ifdef CREASEY_COOLING
- message(
- "Cooling properties for Creasey cooling are (lambda, min_temperature, "
- "hydrogen_mass_abundance, mean_molecular_weight, tstep multiplier) %g %g "
- "%g %g %g",
- cooling->creasey_cooling.lambda, cooling->creasey_cooling.min_temperature,
- cooling->creasey_cooling.hydrogen_mass_abundance,
- cooling->creasey_cooling.mean_molecular_weight,
- cooling->creasey_cooling.cooling_tstep_mult);
-#endif /* CREASEY_COOLING */
-}
-
-void update_entropy(const struct phys_const* const phys_const,
- const struct UnitSystem* us,
- const struct cooling_data* cooling, struct part* p,
- float dt) {
-
- /*updates the entropy of a particle after integrating the cooling equation*/
- float u_old;
- float u_new;
- float new_entropy;
- // float old_entropy = p->entropy;
- float rho = p->rho;
-
- // u_old = old_entropy/(GAMMA_MINUS1) * pow(rho,GAMMA_MINUS1);
- u_old =
- hydro_get_internal_energy(p, 0); // dt = 0 because using current entropy
- u_new = calculate_new_thermal_energy(u_old, rho, dt, cooling, phys_const, us);
- new_entropy = u_new * pow_minus_gamma_minus_one(rho) * hydro_gamma_minus_one;
- p->entropy = new_entropy;
-}
-
-/*This function integrates the cooling equation, given the initial
- thermal energy, density and the timestep dt. Returns the final internal
- energy*/
-
-float calculate_new_thermal_energy(float u_old, float rho, float dt,
- const struct cooling_data* cooling,
- const struct phys_const* const phys_const,
- const struct UnitSystem* us) {
-#ifdef CONST_COOLING
- /*du/dt = -lambda, independent of density*/
- float du_dt = -cooling->const_cooling.lambda;
- float u_floor = cooling->const_cooling.min_energy;
- float u_new;
- if (u_old - du_dt * dt > u_floor) {
- u_new = u_old + du_dt * dt;
- } else {
- u_new = u_floor;
- }
-#endif /*CONST_COOLING*/
-
-#ifdef CREASEY_COOLING
- /* rho*du/dt = -lambda*n_H^2 */
- float u_new;
- float X_H = cooling->creasey_cooling.hydrogen_mass_abundance;
- float lambda_cgs = cooling->creasey_cooling.lambda; // this is always in cgs
- float u_floor_cgs = cooling->creasey_cooling.min_internal_energy_cgs;
-
- /*convert from internal code units to cgs*/
- float dt_cgs = dt * units_cgs_conversion_factor(us, UNIT_CONV_TIME);
- float rho_cgs = rho * units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
- float m_p_cgs = phys_const->const_proton_mass *
- units_cgs_conversion_factor(us, UNIT_CONV_MASS);
- float n_H_cgs = X_H * rho_cgs / m_p_cgs;
- float u_old_cgs =
- u_old * units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
- float du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs;
- float u_new_cgs;
-
- if (u_old_cgs + du_dt_cgs * dt_cgs > u_floor_cgs) {
- u_new_cgs = u_old_cgs + du_dt_cgs * dt_cgs;
- } else {
- u_new_cgs = u_floor_cgs;
- }
- /*convert back to internal code units when returning new internal energy*/
-
- u_new = u_new_cgs /
- units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
-
-#endif /*CREASEY_COOLING*/
- return u_new;
-}
diff --git a/src/cooling.h b/src/cooling.h
deleted file mode 100644
index 7493e73dd3888983af7bd094e4e3b660053e378d..0000000000000000000000000000000000000000
--- a/src/cooling.h
+++ /dev/null
@@ -1,103 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- * Richard Bower (r.g.bower@durham.ac.uk)
- * Stefan Arridge (stefan.arridge@durham.ac.uk)
- *
- * 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_H
-#define SWIFT_COOLING_H
-
-/* Config parameters. */
-#include "../config.h"
-
-/* Some standard headers. */
-#include <math.h>
-
-/* Local includes. */
-#include "const.h"
-#include "error.h"
-#include "hydro.h"
-#include "parser.h"
-#include "part.h"
-#include "physical_constants.h"
-#include "units.h"
-
-/* Cooling Properties */
-struct cooling_data {
-
-#ifdef CONST_COOLING
- struct {
- float lambda;
- float min_energy;
- float cooling_tstep_mult;
- } const_cooling;
-#endif
-
-#ifdef CREASEY_COOLING
- struct {
- float lambda;
- float min_temperature;
- float hydrogen_mass_abundance;
- float mean_molecular_weight;
- float min_internal_energy;
- float min_internal_energy_cgs;
- float cooling_tstep_mult;
- } creasey_cooling;
-#endif
-};
-
-/* Include Cooling */
-#ifdef CONST_COOLING
-
-/**
- * @brief Computes the time-step due to cooling
- *
- * @param cooling The #cooling_data used in the run.
- * @param phys_const The physical constants in internal units.
- * @param Pointer to the particle data.
- */
-__attribute__((always_inline)) INLINE static double cooling_timestep(
- const struct cooling_data* cooling,
- const struct phys_const* const phys_const, const struct part* const p) {
-
- const double cooling_rate = cooling->const_cooling.lambda;
- const double internal_energy =
- hydro_get_internal_energy(p, 0); // dt = 0 because using current entropy
- return cooling->const_cooling.cooling_tstep_mult * internal_energy /
- cooling_rate;
-}
-
-#endif /* CONST_COOLING */
-
-/* Now, some generic functions, defined in the source file */
-void cooling_init(const struct swift_params* parameter_file,
- struct UnitSystem* us,
- const struct phys_const* const phys_const,
- struct cooling_data* cooling);
-
-void cooling_print(const struct cooling_data* cooling);
-void update_entropy(const struct phys_const* const phys_const,
- const struct UnitSystem* us,
- const struct cooling_data* cooling, struct part* p,
- float dt);
-float calculate_new_thermal_energy(float u_old, float rho, float dt,
- const struct cooling_data* cooling,
- const struct phys_const* const phys_const,
- const struct UnitSystem* us);
-#endif /* SWIFT_COOLING_H */
diff --git a/src/cooling/const/cooling.h b/src/cooling/const/cooling.h
new file mode 100644
index 0000000000000000000000000000000000000000..7e71d58f78bb8954536f6eae74eaf9a757d42025
--- /dev/null
+++ b/src/cooling/const/cooling.h
@@ -0,0 +1,140 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Richard Bower (r.g.bower@durham.ac.uk)
+ * Stefan Arridge (stefan.arridge@durham.ac.uk)
+ *
+ * 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_CONST_COOLING_H
+#define SWIFT_CONST_COOLING_H
+
+/**
+ * @file src/cooling/const/cooling.h
+ * @brief Routines related to the "constant cooling" cooling function.
+ *
+ * This is the simplest possible cooling function. A constant cooling rate with
+ * a minimal energy floor is applied.
+ */
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local includes. */
+#include "const.h"
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/**
+ * @brief Properties of the cooling function.
+ */
+struct cooling_data {
+
+ /*! Cooling rate in internal units */
+ float lambda;
+
+ /*! Minimal internal energy of the particles */
+ float min_energy;
+
+ /*! Constant multiplication factor for time-step criterion */
+ float cooling_tstep_mult;
+};
+
+/**
+ * @brief Apply the cooling function to a particle.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_data used in the run.
+ * @param p Pointer to the particle data.
+ * @param dt The time-step of this particle.
+ */
+__attribute__((always_inline)) INLINE static void cooling_cool_part(
+ const struct phys_const* const phys_const, const struct UnitSystem* us,
+ const struct cooling_data* cooling, struct part* p, float dt) {
+
+ /* Get current internal energy (dt=0) */
+ const float u_old = hydro_get_internal_energy(p, 0.f);
+
+ /* Get cooling function properties */
+ const float du_dt = -cooling->lambda;
+ const float u_floor = cooling->min_energy;
+
+ /* Constant cooling with a minimal floor */
+ float u_new;
+ if (u_old - du_dt * dt > u_floor) {
+ u_new = u_old + du_dt * dt;
+ } else {
+ u_new = u_floor;
+ }
+
+ /* Update the internal energy */
+ hydro_set_internal_energy(p, u_new);
+}
+
+/**
+ * @brief Computes the cooling time-step.
+ *
+ * @param cooling The #cooling_data used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static double cooling_timestep(
+ const struct cooling_data* cooling,
+ const struct phys_const* const phys_const, const struct part* const p) {
+
+ const float cooling_rate = cooling->lambda;
+ const float internal_energy =
+ hydro_get_internal_energy(p, 0); // dt = 0 because using current entropy
+ return cooling->cooling_tstep_mult * internal_energy / cooling_rate;
+}
+
+/**
+ * @brief Initialises the cooling properties in the internal system
+ * of units.
+ *
+ * @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 cooling The cooling properties to initialize
+ */
+inline inline void cooling_init(const struct swift_params* parameter_file,
+ const struct UnitSystem* us,
+ const struct phys_const* phys_const,
+ struct cooling_data* cooling) {
+
+ cooling->lambda = parser_get_param_double(parameter_file, "Cooling:lambda");
+ cooling->min_energy =
+ parser_get_param_double(parameter_file, "Cooling:min_energy");
+ cooling->cooling_tstep_mult =
+ parser_get_param_double(parameter_file, "Cooling:cooling_tstep_mult");
+}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * @param cooling The properties of the cooling function.
+ */
+static inline void cooling_print(const struct cooling_data* cooling) {
+
+ message("Cooling function is 'Constant cooling' with rate %f and floor %f",
+ cooling->lambda, cooling->min_energy);
+}
+
+#endif /* SWIFT_CONST_COOLING_H */
diff --git a/src/engine.c b/src/engine.c
index b0626f7d734457bbb0b521913e75ae40031d58a8..f88005b37d99e8faa4aea11f120629748fc42204 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3105,6 +3105,7 @@ void engine_unpin() {
* @param physical_constants The #phys_const used for this run.
* @param hydro The #hydro_props used for this run.
* @param potential The properties of the external potential.
+ * @param cooling The properties of the cooling function.
*/
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
diff --git a/src/runner.c b/src/runner.c
index 8d8b4760ee0664095f907f0071729ceaac2ed1fd..004ec2e316fbe1e08a1df7a14d3b979df5d7387f 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -163,7 +163,6 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
const struct phys_const *constants = r->e->physical_constants;
const struct UnitSystem *us = r->e->internalUnits;
const double timeBase = r->e->timeBase;
- double dt;
TIMER_TIC;
@@ -186,8 +185,10 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
/* Kick has already updated ti_end, so need to check ti_begin */
if (p->ti_begin == ti_current) {
- dt = (p->ti_end - p->ti_begin) * timeBase;
- update_entropy(constants, us, cooling, p, dt);
+
+ const double dt = (p->ti_end - p->ti_begin) * timeBase;
+
+ cooling_cool_part(constants, us, cooling, p, dt);
}
}