Skip to content
Snippets Groups Projects
Commit 7f01b842 authored by Loic Hausammann's avatar Loic Hausammann
Browse files

Cleaning grackle cooling

parent 7d8ac575
No related branches found
No related tags found
1 merge request!740Grackle cosmo
Hide whitespace changes
Inline Side-by-side
src/cooling/grackle/cooling.h
+ 75
− 52
View file @ 7f01b842
......@@ -24,9 +24,6 @@
* @brief Cooling using the GRACKLE 3.0 library.
*/
/* Config parameters. */
#include "../config.h"
/* Some standard headers. */
#include <float.h>
#include <math.h>
......@@ -45,10 +42,26 @@
#include "physical_constants.h"
#include "units.h"
static const float rounding_tolerance = 1.0e-4;
/* need to rework (and check) code if changed */
#define GRACKLE_NPART 1
#define GRACKLE_RANK 3
/* prototype */
static gr_float cooling_time(
const struct phys_const* restrict phys_const,
const struct unit_system* restrict us,
const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling,
const struct part* restrict p, const struct xpart* restrict xp);
static gr_float cooling_new_energy(
const struct phys_const* restrict phys_const,
const struct unit_system* restrict us,
const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling,
const struct part* restrict p, struct xpart* restrict xp, double dt);
/**
* @brief Common operations performed on the cooling function at a
* given time-step or redshift.
......@@ -66,21 +79,6 @@ INLINE static void cooling_update(const struct cosmology* cosmo,
}
/* prototypes */
static gr_float cooling_time(
const struct phys_const* restrict phys_const,
const struct unit_system* restrict us,
const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling,
const struct part* restrict p, struct xpart* restrict xp);
static double cooling_rate(const struct phys_const* restrict phys_const,
const struct unit_system* restrict us,
const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling,
const struct part* restrict p,
struct xpart* restrict xp, double dt);
/**
* @brief Print the chemical network
*
......@@ -182,7 +180,7 @@ __attribute__((always_inline)) INLINE static void cooling_compute_equilibrium(
const double alpha = 0.01;
double dt =
fabs(cooling_time(phys_const, us, cosmo, &cooling_tmp, &p_tmp, xp));
cooling_rate(phys_const, us, cosmo, &cooling_tmp, &p_tmp, xp, dt);
cooling_new_energy(phys_const, us, cosmo, &cooling_tmp, &p_tmp, xp, dt);
dt = alpha *
fabs(cooling_time(phys_const, us, cosmo, &cooling_tmp, &p_tmp, xp));
......@@ -198,7 +196,7 @@ __attribute__((always_inline)) INLINE static void cooling_compute_equilibrium(
old = *xp;
/* update chemistry */
cooling_rate(phys_const, us, cosmo, &cooling_tmp, &p_tmp, xp, dt);
cooling_new_energy(phys_const, us, cosmo, &cooling_tmp, &p_tmp, xp, dt);
} while (step < max_step && !cooling_converged(xp, &old, conv_limit));
if (step == max_step)
......@@ -295,7 +293,8 @@ __attribute__((always_inline)) INLINE static void cooling_print_backend(
message("\tLength = %g", cooling->units.length_units);
message("\tDensity = %g", cooling->units.density_units);
message("\tTime = %g", cooling->units.time_units);
message("\tScale Factor = %g", cooling->units.a_units);
message("\tScale Factor = %g (units: %g)",
cooling->units.a_value, cooling->units.a_units);
}
/**
......@@ -492,7 +491,7 @@ __attribute__((always_inline)) INLINE static void cooling_print_backend(
cooling_copy_from_grackle3(data, p, xp, rho);
/**
* @brief Compute the cooling rate and update the particle chemistry data
* @brief Compute the energy of a particle after dt and update the particle chemistry data
*
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
......@@ -503,7 +502,7 @@ __attribute__((always_inline)) INLINE static void cooling_print_backend(
*
* @return du / dt
*/
__attribute__((always_inline)) INLINE static gr_float cooling_rate(
__attribute__((always_inline)) INLINE static gr_float cooling_new_energy(
const struct phys_const* restrict phys_const,
const struct unit_system* restrict us,
const struct cosmology* restrict cosmo,
......@@ -530,8 +529,8 @@ __attribute__((always_inline)) INLINE static gr_float cooling_rate(
/* general particle data */
gr_float density = hydro_get_physical_density(p, cosmo);
const double energy_before =
hydro_get_drifted_physical_internal_energy(p, cosmo);
const float energy_before =
hydro_get_physical_internal_energy(p, xp, cosmo);
gr_float energy = energy_before;
/* initialize density */
......@@ -549,18 +548,16 @@ __attribute__((always_inline)) INLINE static gr_float cooling_rate(
cooling_copy_to_grackle(data, p, xp, density);
/* solve chemistry */
chemistry_data chemistry_grackle = cooling->chemistry;
chemistry_data_storage chemistry_rates = grackle_rates;
if (local_solve_chemistry(&chemistry_grackle, &chemistry_rates, &units, &data,
dt) == 0) {
chemistry_data chemistry_grackle = cooling->chemistry;
if (local_solve_chemistry(&chemistry_grackle, &grackle_rates,
&units, &data, dt) == 0) {
error("Error in solve_chemistry.");
}
/* copy from grackle data to particle */
cooling_copy_from_grackle(data, p, xp, density);
/* compute rate */
return (energy - energy_before) / dt;
return energy;
}
/**
......@@ -577,18 +574,14 @@ __attribute__((always_inline)) INLINE static gr_float cooling_time(
const struct unit_system* restrict us,
const struct cosmology* restrict cosmo,
const struct cooling_function_data* restrict cooling,
const struct part* restrict p, struct xpart* restrict xp) {
const struct part* restrict p, const struct xpart* restrict xp) {
/* set current time */
code_units units = cooling->units;
if (cooling->redshift == -1)
error("TODO time dependant redshift");
else
units.a_value = 1. / (1. + cooling->redshift);
/* initialize data */
grackle_field_data data;
/* set values */
/* grid */
int grid_dimension[GRACKLE_RANK] = {GRACKLE_NPART, 1, 1};
......@@ -602,7 +595,7 @@ __attribute__((always_inline)) INLINE static gr_float cooling_time(
/* general particle data */
const gr_float energy_before =
hydro_get_drifted_physical_internal_energy(p, cosmo);
hydro_get_physical_internal_energy(p, xp, cosmo);
gr_float density = hydro_get_physical_density(p, cosmo);
gr_float energy = energy_before;
......@@ -659,19 +652,47 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
struct part* restrict p, struct xpart* restrict xp, double dt,
double dt_therm) {
/* Nothing to do here? */
if (dt == 0.) return;
/* Current du_dt */
/* Current energy */
const float u_old = hydro_get_physical_internal_energy(p, xp, cosmo);
/* Current du_dt in physical coordinates (internal units) */
const float hydro_du_dt = hydro_get_physical_internal_energy_dt(p, cosmo);
/* compute cooling rate */
const float du_dt = cooling_rate(phys_const, us, cosmo, cooling, p, xp, dt);
/* Calculate energy after dt */
gr_float u_new = cooling_new_energy(phys_const, us, cosmo, cooling, p, xp, dt);
float delta_u = u_new - u_old + hydro_du_dt * dt_therm;
/* We now need to check that we are not going to go below any of the limits */
/* First, check whether we may end up below the minimal energy after
* this step 1/2 kick + another 1/2 kick that could potentially be for
* a time-step twice as big. We hence check for 1.5 delta_u. */
if (u_old + 1.5 * delta_u < hydro_props->minimal_internal_energy) {
delta_u = (hydro_props->minimal_internal_energy - u_old) / 1.5;
}
/* Second, check whether the energy used in the prediction could get negative.
* We need to check for the 1/2 dt kick followed by a full time-step drift
* that could potentially be for a time-step twice as big. We hence check
* for 2.5 delta_u but this time against 0 energy not the minimum.
* To avoid numerical rounding bringing us below 0., we add a tiny tolerance.
*/
if (u_old + 2.5 * delta_u < 0.) {
delta_u = -u_old / (2.5 + rounding_tolerance);
}
/* Turn this into a rate of change (including cosmology term) */
const float cooling_du_dt = delta_u / dt_therm;
/* record energy lost */
xp->cooling_data.radiated_energy += -du_dt * dt * hydro_get_mass(p);
/* Update the internal energy time derivative */
hydro_set_physical_internal_energy_dt(p, cosmo, cooling_du_dt);
/* Update the internal energy */
hydro_set_physical_internal_energy_dt(p, cosmo, hydro_du_dt + du_dt);
/* Store the radiated energy */
xp->cooling_data.radiated_energy -= hydro_get_mass(p) * cooling_du_dt * dt;
}
static INLINE float cooling_get_temperature(
......@@ -720,7 +741,7 @@ __attribute__((always_inline)) INLINE static void cooling_init_units(
/* These are conversions from code units to cgs. */
/* first cosmo */
cooling->units.a_units = 1.0; // units for the expansion factor (1/1+zi)
cooling->units.a_units = 1.0; // units for the expansion factor
cooling->units.a_value = 1.0;
/* We assume here all physical quantities to
......@@ -729,12 +750,13 @@ __attribute__((always_inline)) INLINE static void cooling_init_units(
/* then units */
cooling->units.density_units =
us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
cooling->units.length_units = us->UnitLength_in_cgs;
cooling->units.time_units = us->UnitTime_in_cgs;
cooling->units.velocity_units = cooling->units.a_units *
cooling->units.length_units /
cooling->units.time_units;
units_cgs_conversion_factor(us, UNIT_CONV_DENSITY);
cooling->units.length_units =
units_cgs_conversion_factor(us, UNIT_CONV_LENGTH);
cooling->units.time_units =
units_cgs_conversion_factor(us, UNIT_CONV_TIME);
cooling->units.velocity_units =
units_cgs_conversion_factor(us, UNIT_CONV_VELOCITY);
}
/**
......@@ -813,6 +835,7 @@ __attribute__((always_inline)) INLINE static void cooling_init_backend(
/* Set up the units system. */
cooling_init_units(us, cooling);
/* Set up grackle */
cooling_init_grackle(cooling);
}
......
src/cooling/grackle/cooling_io.h
+ 0
− 3
View file @ 7f01b842
......@@ -19,9 +19,6 @@
#ifndef SWIFT_COOLING_GRACKLE_IO_H
#define SWIFT_COOLING_GRACKLE_IO_H
/* Config parameters. */
#include "../config.h"
/* Local includes */
#include "cooling_struct.h"
#include "io_properties.h"
......
src/cooling/grackle/cooling_struct.h
+ 0
− 2
View file @ 7f01b842
......@@ -19,8 +19,6 @@
#ifndef SWIFT_COOLING_STRUCT_GRACKLE_H
#define SWIFT_COOLING_STRUCT_GRACKLE_H
#include "../config.h"
/* include grackle */
#include <grackle.h>
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment