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

Removed unnecessary files. Corrected header guards. Renamed grackle-cooling

parent 3fc24757
No related branches found
No related tags found
1 merge request!242Add the cooling infrastructure and the const_du and const_lambda
Hide whitespace changes
Inline Side-by-side
src/const.h
+ 1
− 1
View file @ d3ab083d
...@@ -97,7 +97,7 @@ ...@@ -97,7 +97,7 @@
/* Cooling properties */ /* Cooling properties */
//#define COOLING_CONST_DU //#define COOLING_CONST_DU
#define COOLING_CONST_LAMBDA #define COOLING_CONST_LAMBDA
//#define COOLING_GRACKLE_COOLING //#define COOLING_GRACKLE
/* Are we debugging ? */ /* Are we debugging ? */
//#define SWIFT_DEBUG_CHECKS //#define SWIFT_DEBUG_CHECKS
......
src/cooling.h
+ 1
− 1
View file @ d3ab083d
...@@ -35,7 +35,7 @@ ...@@ -35,7 +35,7 @@
#include "./cooling/const_du/cooling.h" #include "./cooling/const_du/cooling.h"
#elif defined(COOLING_CONST_LAMBDA) #elif defined(COOLING_CONST_LAMBDA)
#include "./cooling/const_lambda/cooling.h" #include "./cooling/const_lambda/cooling.h"
#elif defined(COOLING_GRACKLE_COOLING) #elif defined(COOLING_GRACKLE)
#include "./cooling/grackle/cooling.h" #include "./cooling/grackle/cooling.h"
#else #else
#error "Invalid choice of cooling function." #error "Invalid choice of cooling function."
......
src/cooling/const_du/cooling.h
+ 5
− 4
View file @ d3ab083d
...@@ -18,8 +18,8 @@ ...@@ -18,8 +18,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. * along with this program. If not, see <http://www.gnu.org/licenses/>.
* *
******************************************************************************/ ******************************************************************************/
#ifndef SWIFT_COOLING_CONST_DU #ifndef SWIFT_COOLING_CONST_DU_H
#define SWIFT_COOLING_CONST_DU #define SWIFT_COOLING_CONST_DU_H
/** /**
* @file src/cooling/const/cooling.h * @file src/cooling/const/cooling.h
...@@ -118,7 +118,8 @@ INLINE void cooling_init(const struct swift_params* parameter_file, ...@@ -118,7 +118,8 @@ INLINE void cooling_init(const struct swift_params* parameter_file,
const struct phys_const* phys_const, const struct phys_const* phys_const,
struct cooling_data* cooling) { 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 = cooling->min_energy =
parser_get_param_double(parameter_file, "Cooling:min_energy"); parser_get_param_double(parameter_file, "Cooling:min_energy");
cooling->cooling_tstep_mult = cooling->cooling_tstep_mult =
...@@ -136,4 +137,4 @@ INLINE void cooling_print(const struct cooling_data* cooling) { ...@@ -136,4 +137,4 @@ INLINE void cooling_print(const struct cooling_data* cooling) {
cooling->cooling_rate, cooling->min_energy); cooling->cooling_rate, cooling->min_energy);
} }
#endif /* SWIFT_CONST_COOLING_H */ #endif /* SWIFT_COOLING_CONST_DU_H */
src/cooling/const_lambda/cooling.c deleted 100644 → 0
+ 0
− 150
View file @ 3fc24757
/*******************************************************************************
* 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;
}
src/cooling/const_lambda/cooling.h
+ 39
− 28
View file @ d3ab083d
...@@ -20,8 +20,8 @@ ...@@ -20,8 +20,8 @@
* *
******************************************************************************/ ******************************************************************************/
#ifndef SWIFT_COOLING_CONST_LAMBDA #ifndef SWIFT_COOLING_CONST_LAMBDA_H
#define SWIFT_COOLING_CONST_LAMBDA #define SWIFT_COOLING_CONST_LAMBDA_H
/* Some standard headers. */ /* Some standard headers. */
#include <math.h> #include <math.h>
...@@ -37,7 +37,7 @@ ...@@ -37,7 +37,7 @@
/* Cooling Properties */ /* Cooling Properties */
struct cooling_data { struct cooling_data {
/*! Cooling rate in cgs units. Defined by 'rho * du/dt = -lambda * n_H^2'*/ /*! Cooling rate in cgs units. Defined by 'rho * du/dt = -lambda * n_H^2'*/
float lambda; float lambda;
...@@ -79,17 +79,19 @@ __attribute__((always_inline)) INLINE static float cooling_rate( ...@@ -79,17 +79,19 @@ __attribute__((always_inline)) INLINE static float cooling_rate(
const float lambda_cgs = cooling->lambda; const float lambda_cgs = cooling->lambda;
/*convert from internal code units to cgs*/ /*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 * 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; const float n_H_cgs = X_H * rho_cgs / m_p_cgs;
/* Calculate du_dt */ /* Calculate du_dt */
const float du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs; const float du_dt_cgs = -lambda_cgs * n_H_cgs * n_H_cgs / rho_cgs;
/* Convert du/dt back to internal code units */ /* Convert du/dt back to internal code units */
const float du_dt = du_dt_cgs * units_cgs_conversion_factor(us, UNIT_CONV_TIME)/ const float du_dt =
units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS); 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; return du_dt;
} }
...@@ -107,7 +109,6 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( ...@@ -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 phys_const* const phys_const, const struct UnitSystem* us,
const struct cooling_data* cooling, struct part* p, float dt) { const struct cooling_data* cooling, struct part* p, float dt) {
/* Get current internal energy (dt=0) */ /* Get current internal energy (dt=0) */
const float u_old = hydro_get_internal_energy(p, 0.f); const float u_old = hydro_get_internal_energy(p, 0.f);
...@@ -115,7 +116,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( ...@@ -115,7 +116,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
const float u_floor = cooling->min_energy; const float u_floor = cooling->min_energy;
/* Calculate du_dt */ /* 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 */ /* Intergrate cooling equation, but enforce energy floor */
float u_new; float u_new;
if (u_old + du_dt * dt > u_floor) { if (u_old + du_dt * dt > u_floor) {
...@@ -123,17 +124,18 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( ...@@ -123,17 +124,18 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
} else { } else {
u_new = u_floor; u_new = u_floor;
} }
/* Update the internal energy */ /* Update the internal energy */
hydro_set_internal_energy(p, u_new); 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){ */ /* 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); */ /* exit(-1); */
/* } */ /* } */
} }
/** /**
* @brief Computes the time-step due to cooling * @brief Computes the time-step due to cooling
* *
...@@ -143,16 +145,17 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part( ...@@ -143,16 +145,17 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
* @param Pointer to the particle data. * @param Pointer to the particle data.
*/ */
__attribute__((always_inline)) INLINE static float cooling_timestep( __attribute__((always_inline)) INLINE static float cooling_timestep(
const struct cooling_data* cooling, const struct phys_const* const phys_const, const struct cooling_data* cooling,
const struct UnitSystem* us, const struct part* const p) { const struct phys_const* const phys_const, const struct UnitSystem* us,
const struct part* const p) {
/* Get du_dt */ /* 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) */ /* Get current internal energy (dt=0) */
const float u = hydro_get_internal_energy(p, 0.f); 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, ...@@ -169,15 +172,20 @@ INLINE void cooling_init(const struct swift_params* parameter_file,
struct cooling_data* cooling) { struct cooling_data* cooling) {
cooling->lambda = parser_get_param_double(parameter_file, "Cooling:lambda"); cooling->lambda = parser_get_param_double(parameter_file, "Cooling:lambda");
cooling->min_temperature = parser_get_param_double(parameter_file, "Cooling:minimum_temperature"); cooling->min_temperature =
cooling->hydrogen_mass_abundance = parser_get_param_double(parameter_file, "Cooling:hydrogen_mass_abundance"); 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(
cooling->cooling_tstep_mult = parser_get_param_double(parameter_file, "Cooling:cooling_tstep_mult"); 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*/ /*convert minimum temperature into minimum internal energy*/
const float u_floor = const float u_floor =
phys_const->const_boltzmann_k * cooling->min_temperature / phys_const->const_boltzmann_k * cooling->min_temperature /
(hydro_gamma_minus_one * cooling->mean_molecular_weight * phys_const->const_proton_mass); (hydro_gamma_minus_one * cooling->mean_molecular_weight *
phys_const->const_proton_mass);
const float u_floor_cgs = const float u_floor_cgs =
u_floor * units_cgs_conversion_factor(us, UNIT_CONV_ENERGY_PER_UNIT_MASS); 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, ...@@ -192,9 +200,12 @@ INLINE void cooling_init(const struct swift_params* parameter_file,
*/ */
INLINE void cooling_print(const struct cooling_data* cooling) { 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)", message(
cooling->lambda, cooling->min_temperature,cooling->hydrogen_mass_abundance,cooling->mean_molecular_weight); "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 */
src/runner.c
+ 5
− 4
View file @ d3ab083d
...@@ -160,7 +160,8 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) { ...@@ -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 r runner task
* @param c cell * @param c cell
...@@ -1330,9 +1331,9 @@ void *runner_main(void *data) { ...@@ -1330,9 +1331,9 @@ void *runner_main(void *data) {
case task_type_grav_external: case task_type_grav_external:
runner_do_grav_external(r, t->ci, 1); runner_do_grav_external(r, t->ci, 1);
break; break;
case task_type_cooling: case task_type_cooling:
runner_do_cooling(r, t->ci, 1); runner_do_cooling(r, t->ci, 1);
break; break;
default: default:
error("Unknown task type."); error("Unknown task type.");
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment