Commit 33632ee6 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Moved the cooling functions to a separate directory

parent 47883b9f
......@@ -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
......
......@@ -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)
......
......@@ -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
......
/*******************************************************************************
* 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;
}
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
* Matthieu Schaller (matthieu.schaller@durham.ac.uk)
* 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)
*
......@@ -19,12 +18,16 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_CONST_COOLING_H
#define SWIFT_CONST_COOLING_H
#ifndef SWIFT_COOLING_H
#define SWIFT_COOLING_H
/* Config parameters. */
#include "../config.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>
......@@ -38,66 +41,100 @@
#include "physical_constants.h"
#include "units.h"
/* Cooling Properties */
/**
* @brief Properties of the cooling function.
*/
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
/*! 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;
};
/* Include Cooling */
#ifdef CONST_COOLING
/**
* @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 time-step due to cooling
* @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 Pointer to the particle data.
* @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 double cooling_rate = cooling->const_cooling.lambda;
const double internal_energy =
const float cooling_rate = cooling->lambda;
const float 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;
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 /* 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 */
#endif /* SWIFT_CONST_COOLING_H */
......@@ -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,
......
......@@ -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);
}
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment