diff --git a/src/const.h b/src/const.h
index cbd814f4a9646e0709b6c4ebd9b4ac3c19be2f5b..552d49c8f1e4dd9f4fa2855e5fed548acd5c0b3f 100644
--- a/src/const.h
+++ b/src/const.h
@@ -97,7 +97,7 @@
/* Cooling properties */
//#define COOLING_CONST_DU
#define COOLING_CONST_LAMBDA
-//#define COOLING_GRACKLE_COOLING
+//#define COOLING_GRACKLE
/* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS
diff --git a/src/cooling.h b/src/cooling.h
index 691b2fd7008a6cdd8c6688bd5fe644d5b48ac897..d2dac6a804675dc7aeec0cfde6884414aa41ad31 100644
--- a/src/cooling.h
+++ b/src/cooling.h
@@ -35,7 +35,7 @@
#include "./cooling/const_du/cooling.h"
#elif defined(COOLING_CONST_LAMBDA)
#include "./cooling/const_lambda/cooling.h"
-#elif defined(COOLING_GRACKLE_COOLING)
+#elif defined(COOLING_GRACKLE)
#include "./cooling/grackle/cooling.h"
#else
#error "Invalid choice of cooling function."
diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h
index a3b4714fe597b3f4098e6f335a0a43026f514ae5..382352a31217c343810c66467747f86bc2cab84b 100644
--- a/src/cooling/const_du/cooling.h
+++ b/src/cooling/const_du/cooling.h
@@ -18,8 +18,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_COOLING_CONST_DU
-#define SWIFT_COOLING_CONST_DU
+#ifndef SWIFT_COOLING_CONST_DU_H
+#define SWIFT_COOLING_CONST_DU_H
/**
* @file src/cooling/const/cooling.h
@@ -118,7 +118,8 @@ INLINE void cooling_init(const struct swift_params* parameter_file,
const struct phys_const* phys_const,
struct cooling_data* cooling) {
- cooling->cooling_rate = parser_get_param_double(parameter_file, "Cooling:cooling_rate");
+ cooling->cooling_rate =
+ parser_get_param_double(parameter_file, "Cooling:cooling_rate");
cooling->min_energy =
parser_get_param_double(parameter_file, "Cooling:min_energy");
cooling->cooling_tstep_mult =
@@ -136,4 +137,4 @@ INLINE void cooling_print(const struct cooling_data* cooling) {
cooling->cooling_rate, cooling->min_energy);
}
-#endif /* SWIFT_CONST_COOLING_H */
+#endif /* SWIFT_COOLING_CONST_DU_H */
diff --git a/src/cooling/const_lambda/cooling.c b/src/cooling/const_lambda/cooling.c
deleted file mode 100644
index 1daa5d52d5794583d1818ffa0d1f7acec83998b5..0000000000000000000000000000000000000000
--- a/src/cooling/const_lambda/cooling.c
+++ /dev/null
@@ -1,150 +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) {
-
- cooling->lambda = parser_get_param_double(parameter_file, "Cooling:Lambda");
- cooling->min_temperature = parser_get_param_double(parameter_file, "Cooling:minimum_temperature");
- cooling->mean_molecular_weight = parser_get_param_double(parameter_file, "Cooling:mean_molecular_weight");
- cooling->hydrogen_mass_abundance = parser_get_param_double(parameter_file, "Cooling:hydrogen_mass_abundance");
- cooling->cooling_tstep_mult = parser_get_param_double(parameter_file, "Cooling:cooling_tstep_mult");
-
- /*convert minimum temperature into minimum internal energy*/
- float u_floor =
- phys_const->const_boltzmann_k * cooling->min_temperature /
- (hydro_gamma_minus_one * 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->min_internal_energy = u_floor;
- cooling->min_internal_energy_cgs = u_floor_cgs;
-}
-
-/**
- * @brief Prints the properties of the cooling model to stdout.
- *
- * @param cooling The cooling properties.
- */
-void cooling_print(const struct cooling_data* cooling) {
-
- message(
- "Cooling properties are (lambda, min_temperature, "
- "hydrogen_mass_abundance, mean_molecular_weight, tstep multiplier) %g %g "
- "%g %g %g",
- cooling->lambda,
- cooling->min_temperature,
- cooling->hydrogen_mass_abundance,
- cooling->mean_molecular_weight,
- cooling->cooling_tstep_mult);
-}
-
-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->hydrogen_mass_abundance;
- float lambda_cgs = cooling->lambda; // this is always in cgs
- float u_floor_cgs = 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/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index d5d7fe136605bb2acd44ea162388afab81f568e1..88ac8977422ca37522670b91a7c01bf4ff87a1f9 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -20,8 +20,8 @@
*
******************************************************************************/
-#ifndef SWIFT_COOLING_CONST_LAMBDA
-#define SWIFT_COOLING_CONST_LAMBDA
+#ifndef SWIFT_COOLING_CONST_LAMBDA_H
+#define SWIFT_COOLING_CONST_LAMBDA_H
/* Some standard headers. */
#include <math.h>
@@ -37,7 +37,7 @@
/* Cooling Properties */
struct cooling_data {
-
+
/*! Cooling rate in cgs units. Defined by 'rho * du/dt = -lambda * n_H^2'*/
float lambda;
@@ -79,17 +79,19 @@ __attribute__((always_inline)) INLINE static float cooling_rate(
const float lambda_cgs = cooling->lambda;
/*convert from internal code units to cgs*/
- const float rho_cgs = rho * units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
+ const float rho_cgs =
+ rho * units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
const float m_p_cgs = phys_const->const_proton_mass *
- units_cgs_conversion_factor(us, UNIT_CONV_MASS);
+ units_cgs_conversion_factor(us, UNIT_CONV_MASS);
const float n_H_cgs = X_H * rho_cgs / m_p_cgs;
/* Calculate du_dt */
const float du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs;
/* Convert du/dt back to internal code units */
- const float du_dt = du_dt_cgs * units_cgs_conversion_factor(us, UNIT_CONV_TIME)/
- units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
+ const float du_dt =
+ du_dt_cgs * units_cgs_conversion_factor(us, UNIT_CONV_TIME) /
+ units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
return du_dt;
}
@@ -107,7 +109,6 @@ __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);
@@ -115,7 +116,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
const float u_floor = cooling->min_energy;
/* Calculate du_dt */
- const float du_dt = cooling_rate(phys_const,us,cooling,p);
+ const float du_dt = cooling_rate(phys_const, us, cooling, p);
/* Intergrate cooling equation, but enforce energy floor */
float u_new;
if (u_old + du_dt * dt > u_floor) {
@@ -123,17 +124,18 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
} else {
u_new = u_floor;
}
-
+
/* Update the internal energy */
hydro_set_internal_energy(p, u_new);
- //const float u_new_test = hydro_get_internal_energy(p, 0.f);
+ // const float u_new_test = hydro_get_internal_energy(p, 0.f);
/* if (-(u_new_test - u_old)/u_old > 1.0e-6){ */
- /* printf("Particle has successfully cooled: u_old = %g , du_dt = %g , dt = %g , du_dt*dt = %g, u_old + du_dt*dt = %g, u_new = %g\n",u_old,du_dt,dt,du_dt*dt,u_new,u_new_test); */
+ /* printf("Particle has successfully cooled: u_old = %g , du_dt = %g , dt = %g
+ * , du_dt*dt = %g, u_old + du_dt*dt = %g, u_new =
+ * %g\n",u_old,du_dt,dt,du_dt*dt,u_new,u_new_test); */
/* exit(-1); */
/* } */
}
-
/**
* @brief Computes the time-step due to cooling
*
@@ -143,16 +145,17 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
* @param Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float cooling_timestep(
- const struct cooling_data* cooling, const struct phys_const* const phys_const,
- const struct UnitSystem* us, const struct part* const p) {
+ const struct cooling_data* cooling,
+ const struct phys_const* const phys_const, const struct UnitSystem* us,
+ const struct part* const p) {
/* Get du_dt */
- const float du_dt = cooling_rate(phys_const,us,cooling,p);
-
+ const float du_dt = cooling_rate(phys_const, us, cooling, p);
+
/* Get current internal energy (dt=0) */
const float u = hydro_get_internal_energy(p, 0.f);
- return u / du_dt;
+ return u / du_dt;
}
/**
@@ -169,15 +172,20 @@ INLINE void cooling_init(const struct swift_params* parameter_file,
struct cooling_data* cooling) {
cooling->lambda = parser_get_param_double(parameter_file, "Cooling:lambda");
- cooling->min_temperature = parser_get_param_double(parameter_file, "Cooling:minimum_temperature");
- cooling->hydrogen_mass_abundance = parser_get_param_double(parameter_file, "Cooling:hydrogen_mass_abundance");
- cooling->mean_molecular_weight = parser_get_param_double(parameter_file, "Cooling:mean_molecular_weight");
- cooling->cooling_tstep_mult = parser_get_param_double(parameter_file, "Cooling:cooling_tstep_mult");
+ cooling->min_temperature =
+ parser_get_param_double(parameter_file, "Cooling:minimum_temperature");
+ cooling->hydrogen_mass_abundance = parser_get_param_double(
+ parameter_file, "Cooling:hydrogen_mass_abundance");
+ cooling->mean_molecular_weight =
+ parser_get_param_double(parameter_file, "Cooling:mean_molecular_weight");
+ cooling->cooling_tstep_mult =
+ parser_get_param_double(parameter_file, "Cooling:cooling_tstep_mult");
/*convert minimum temperature into minimum internal energy*/
const float u_floor =
- phys_const->const_boltzmann_k * cooling->min_temperature /
- (hydro_gamma_minus_one * cooling->mean_molecular_weight * phys_const->const_proton_mass);
+ phys_const->const_boltzmann_k * cooling->min_temperature /
+ (hydro_gamma_minus_one * cooling->mean_molecular_weight *
+ phys_const->const_proton_mass);
const float u_floor_cgs =
u_floor * units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
@@ -192,9 +200,12 @@ INLINE void cooling_init(const struct swift_params* parameter_file,
*/
INLINE void cooling_print(const struct cooling_data* cooling) {
- message("Cooling function is 'Constant lambda' with (lambda,min_temperature,hydrogen_mass_abundance,mean_molecular_weight) = (%g,%g,%g,%g)",
- cooling->lambda, cooling->min_temperature,cooling->hydrogen_mass_abundance,cooling->mean_molecular_weight);
+ message(
+ "Cooling function is 'Constant lambda' with "
+ "(lambda,min_temperature,hydrogen_mass_abundance,mean_molecular_weight) "
+ "= (%g,%g,%g,%g)",
+ cooling->lambda, cooling->min_temperature,
+ cooling->hydrogen_mass_abundance, cooling->mean_molecular_weight);
}
-#endif /* SWIFT_COOLING_CONST_LAMBDA */
-
+#endif /* SWIFT_COOLING_CONST_LAMBDA_H */
diff --git a/src/runner.c b/src/runner.c
index 44ab5d5d3f31b6a9b6a1d1fdbb9bfea17419536a..352a64664460343dd9004d60a9102b4875fce600 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -160,7 +160,8 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
}
/**
- * @brief Calculate change in entropy from cooling
+ * @brief Calculate change in thermal state of particles induced
+ * by radiative cooling and heating.
*
* @param r runner task
* @param c cell
@@ -1330,9 +1331,9 @@ void *runner_main(void *data) {
case task_type_grav_external:
runner_do_grav_external(r, t->ci, 1);
break;
- case task_type_cooling:
- runner_do_cooling(r, t->ci, 1);
- break;
+ case task_type_cooling:
+ runner_do_cooling(r, t->ci, 1);
+ break;
default:
error("Unknown task type.");
}