Skip to content

GitLab

Commit c92464e2 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'grackle_rework' into 'master' 

Rework of Grackle

See merge request !484
parents 19ca1efd 7bca4ea7
Hide whitespace changes
Inline Side-by-side
Showing with 231 additions and 391 deletions
configure.ac
View file @ c92464e2
......@@ -478,8 +478,8 @@ if test "x$with_grackle" != "xno"; then
AC_CHECK_LIB(
[grackle],
[initialize_grackle],
[AC_DEFINE([HAVE_GRACKLE],1,[The GRACKLE library appears to be present.]),
[initialize_chemistry_data],
[AC_DEFINE([HAVE_GRACKLE],1,[The GRACKLE library appears to be present.])
AC_DEFINE([CONFIG_BFLOAT_8],1,[Use doubles in grackle])
],
[AC_MSG_ERROR(Cannot find grackle library!)],
......
src/Makefile.am
View file @ c92464e2
......@@ -48,11 +48,6 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h
GRACKLE_SRC =
if HAVEGRACKLE
GRACKLE_SRC += cooling/grackle/grackle_wrapper.c
endif
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
......@@ -62,8 +57,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
part_type.c xmf.c gravity_properties.c gravity.c \
collectgroup.c hydro_space.c equation_of_state.c \
$(GRACKLE_SRC)
collectgroup.c hydro_space.c equation_of_state.c
# Include files for distribution, not installation.
......@@ -123,7 +117,6 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
cooling/const_du/cooling.h cooling/const_du/cooling_struct.h \
cooling/const_lambda/cooling.h cooling/const_lambda/cooling_struct.h \
cooling/grackle/cooling.h cooling/grackle/cooling_struct.h \
cooling/grackle/grackle_wrapper.h \
memswap.h dump.h logger.h
......
src/cooling/grackle/cooling.h
View file @ c92464e2
......@@ -27,6 +27,7 @@
/* Some standard headers. */
#include <float.h>
#include <math.h>
#include <grackle.h>
/* Local includes. */
#include "error.h"
......@@ -36,13 +37,67 @@
#include "physical_constants.h"
#include "units.h"
/* include the grackle wrapper */
#include "grackle_wrapper.h"
/* need to rework (and check) code if changed */
#define GRACKLE_NPART 1
#define GRACKLE_RANK 3
/**
* @brief Compute the cooling rate
* @brief Sets the cooling properties of the (x-)particles to a valid start
* state.
*
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
*/
__attribute__((always_inline)) INLINE static void cooling_init_part(
const struct part* restrict p, struct xpart* restrict xp) {
xp->cooling_data.radiated_energy = 0.f;
}
/**
* @brief Returns the total radiated energy by this particle.
*
* @param xp The extended particle data
*/
__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
const struct xpart* restrict xp) {
return xp->cooling_data.radiated_energy;
}
/**
* @brief Prints the properties of the cooling model to stdout.
*
* We do nothing.
* @param cooling The properties of the cooling function.
*/
__attribute__((always_inline))INLINE static void cooling_print_backend(
const struct cooling_function_data* cooling) {
message("Cooling function is 'Grackle'.");
message("Using Grackle = %i", cooling->chemistry.use_grackle);
message("Chemical network = %i", cooling->chemistry.primordial_chemistry);
message("Radiative cooling = %i", cooling->chemistry.with_radiative_cooling);
message("Metal cooling = %i", cooling->chemistry.metal_cooling);
message("CloudyTable = %s",
cooling->cloudy_table);
message("UVbackground = %d", cooling->uv_background);
message("Redshift = %g", cooling->redshift);
message("Density Self Shielding = %g",
cooling->density_self_shielding);
message("Units:");
message("\tComoving = %i", cooling->units.comoving_coordinates);
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);
}
/**
* @brief Compute the cooling rate
*
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
......@@ -58,33 +113,103 @@ __attribute__((always_inline)) INLINE static double cooling_rate(
const struct cooling_function_data* restrict cooling,
struct part* restrict p, float dt) {
if (cooling->GrackleRedshift == -1) error("TODO time dependant redshift");
/* Get current internal energy (dt=0) */
double u_old = hydro_get_internal_energy(p);
double u_new = u_old;
/* Get current density */
const float rho = hydro_get_density(p);
/* Actual scaling fractor */
const float a_now = 1. / (1. + cooling->GrackleRedshift);
/* 0.02041 (= 1 Zsun in Grackle v2.0, but = 1.5761 Zsun in
Grackle v2.1) */
double Z = 0.02041;
if (wrap_do_cooling(rho, &u_new, dt, Z, a_now) == 0) {
error("Error in do_cooling.\n");
return 0;
if (cooling->chemistry.primordial_chemistry > 1)
error("Not implemented");
/* 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};
int grid_start[GRACKLE_RANK] = {0, 0, 0};
int grid_end[GRACKLE_RANK] = {GRACKLE_NPART - 1, 0, 0};
data.grid_rank = GRACKLE_RANK;
data.grid_dimension = grid_dimension;
data.grid_start = grid_start;
data.grid_end = grid_end;
/* general particle data */
gr_float density = hydro_get_density(p);
const double energy_before = hydro_get_internal_energy(p);
gr_float energy = energy_before;
gr_float vx = 0;
gr_float vy = 0;
gr_float vz = 0;
data.density = &density;
data.internal_energy = &energy;
data.x_velocity = &vx;
data.y_velocity = &vy;
data.z_velocity = &vz;
/* /\* primordial chemistry >= 1 *\/ */
/* gr_float HI_density = density; */
/* gr_float HII_density = 0.; */
/* gr_float HeI_density = 0.; */
/* gr_float HeII_density = 0.; */
/* gr_float HeIII_density = 0.; */
/* gr_float e_density = 0.; */
/* data.HI_density = &HI_density; */
/* data.HII_density = &HII_density; */
/* data.HeI_density = &HeI_density; */
/* data.HeII_density = &HeII_density; */
/* data.HeIII_density = &HeIII_density; */
/* data.e_density = &e_density; */
/* /\* primordial chemistry >= 2 *\/ */
/* gr_float HM_density = 0.; */
/* gr_float H2I_density = 0.; */
/* gr_float H2II_density = 0.; */
/* data.HM_density = &HM_density; */
/* data.H2I_density = &H2I_density; */
/* data.H2II_density = &H2II_density; */
/* /\* primordial chemistry >= 3 *\/ */
/* gr_float DI_density = 0.; */
/* gr_float DII_density = 0.; */
/* gr_float HDI_density = 0.; */
/* data.DI_density = &DI_density; */
/* data.DII_density = &DII_density; */
/* data.HDI_density = &HDI_density; */
/* metal cooling = 1 */
gr_float metal_density = density * grackle_data->SolarMetalFractionByMass;
data.metal_density = &metal_density;
/* /\* volumetric heating rate *\/ */
/* gr_float volumetric_heating_rate = 0.; */
/* data.volumetric_heating_rate = &volumetric_heating_rate; */
/* /\* specific heating rate *\/ */
/* gr_float specific_heating_rate = 0.; */
/* data.specific_heating_rate = &specific_heating_rate; */
/* solve chemistry with table */
if (solve_chemistry(&units, &data, dt) == 0) {
error("Error in solve_chemistry.");
}
return (u_new - u_old) / dt;
return (energy - energy_before) / dt;
}
/**
* @brief Apply the cooling function to a particle.
*
* We do nothing.
*
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
* @param cooling The #cooling_function_data used in the run.
......@@ -102,15 +227,11 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
/* Current du_dt */
const float hydro_du_dt = hydro_get_internal_energy_dt(p);
float du_dt;
float delta_u;
du_dt = cooling_rate(phys_const, us, cooling, p, dt);
delta_u = du_dt * dt;
/* compute cooling rate */
const float du_dt = cooling_rate(phys_const, us, cooling, p, dt);
/* record energy lost */
xp->cooling_data.radiated_energy += -delta_u * hydro_get_mass(p);
xp->cooling_data.radiated_energy += - du_dt * dt * hydro_get_mass(p);
/* Update the internal energy */
hydro_set_internal_energy_dt(p, hydro_du_dt + du_dt);
......@@ -134,30 +255,6 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
return FLT_MAX;
}
/**
* @brief Sets the cooling properties of the (x-)particles to a valid start
* state.
*
* @param p Pointer to the particle data.
* @param xp Pointer to the extended particle data.
*/
__attribute__((always_inline)) INLINE static void cooling_init_part(
const struct part* restrict p, struct xpart* restrict xp) {
xp->cooling_data.radiated_energy = 0.f;
}
/**
* @brief Returns the total radiated energy by this particle.
*
* @param xp The extended particle data
*/
__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
const struct xpart* restrict xp) {
return xp->cooling_data.radiated_energy;
}
/**
* @brief Initialises the cooling properties.
*
......@@ -166,33 +263,76 @@ __attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
* @param phys_const The physical constants in internal units.
* @param cooling The cooling properties to initialize
*/
static INLINE void cooling_init_backend(
__attribute__((always_inline))INLINE static void cooling_init_backend(
const struct swift_params* parameter_file, const struct unit_system* us,
const struct phys_const* phys_const,
struct cooling_function_data* cooling) {
double units_density, units_length, units_time;
int grackle_chemistry;
int UVbackground;
/* read parameters */
parser_get_param_string(parameter_file, "GrackleCooling:GrackleCloudyTable",
cooling->GrackleCloudyTable);
cooling->UVbackground =
parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");
cooling->GrackleRedshift =
parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");
cooling->GrackleHSShieldingDensityThreshold = parser_get_param_double(
cooling->cloudy_table);
cooling->uv_background =
parser_get_param_int(parameter_file, "GrackleCooling:UVbackground");
cooling->redshift =
parser_get_param_double(parameter_file, "GrackleCooling:GrackleRedshift");
cooling->density_self_shielding = parser_get_param_double(
parameter_file, "GrackleCooling:GrackleHSShieldingDensityThreshold");
#ifdef SWIFT_DEBUG_CHECKS
/* enable verbose for grackle */
grackle_verbose = 1;
#endif
UVbackground = cooling->UVbackground;
grackle_chemistry = 0; /* forced to be zero : read table */
/* Set up the units system.
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_value = 1.0;
/* We assume here all physical quantities to
be in proper coordinate (not comobile) */
cooling->units.comoving_coordinates = 0;
/* 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;
chemistry_data *chemistry = &cooling->chemistry;
/* Create a chemistry object for parameters and rate data. */
if (set_default_chemistry_parameters(chemistry) == 0) {
error("Error in set_default_chemistry_parameters.");
}
// Set parameter values for chemistry.
chemistry->use_grackle = 1;
chemistry->with_radiative_cooling = 1;
/* molecular network with H, He, D
From Cloudy table */
chemistry->primordial_chemistry = 0;
chemistry->metal_cooling = 1; // metal cooling on
chemistry->UVbackground = cooling->uv_background;
chemistry->grackle_data_file = cooling->cloudy_table;
chemistry->use_radiative_transfer = 0;
chemistry->use_volumetric_heating_rate = 0;
chemistry->use_specific_heating_rate = 0;
/* Initialize the chemistry object. */
if (initialize_chemistry_data(&cooling->units) == 0) {
error("Error in initialize_chemistry_data.");
}
units_density = us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
units_length = us->UnitLength_in_cgs;
units_time = us->UnitTime_in_cgs;
#ifdef SWIFT_DEBUG_CHECKS
float threshold = cooling->GrackleHSShieldingDensityThreshold;
if (GRACKLE_NPART != 1)
error("Grackle with multiple particles not implemented");
float threshold = cooling->density_self_shielding;
threshold /= phys_const->const_proton_mass;
threshold /= pow(us->UnitLength_in_cgs, 3);
......@@ -200,41 +340,14 @@ static INLINE void cooling_init_backend(
message("***************************************");
message("initializing grackle cooling function");
message("");
message("CloudyTable = %s",
cooling->GrackleCloudyTable);
message("UVbackground = %d", UVbackground);
message("GrackleRedshift = %g", cooling->GrackleRedshift);
message("GrackleHSShieldingDensityThreshold = %g atom/cm3", threshold);
#endif
cooling_print_backend(cooling);
message("Density Self Shielding = %g atom/cm3", threshold);
if (wrap_init_cooling(cooling->GrackleCloudyTable, UVbackground,
units_density, units_length, units_time,
grackle_chemistry) != 1) {
error("Error in initialize_chemistry_data.");
}
#ifdef SWIFT_DEBUG_CHECKS
grackle_print_data();
message("");
message("***************************************");
#endif
}
/**
* @brief Prints the properties of the cooling model to stdout.
*
* @param cooling The properties of the cooling function.
*/
static INLINE void cooling_print_backend(
const struct cooling_function_data* cooling) {
message("Cooling function is 'Grackle'.");
message("CloudyTable = %s",
cooling->GrackleCloudyTable);
message("UVbackground = %d", cooling->UVbackground);
message("GrackleRedshift = %g", cooling->GrackleRedshift);
message("GrackleHSShieldingDensityThreshold = %g atom/cm3",
cooling->GrackleHSShieldingDensityThreshold);
}
#endif /* SWIFT_COOLING_GRACKLE_H */
src/cooling/grackle/cooling_struct.h
View file @ c92464e2
......@@ -19,6 +19,9 @@
#ifndef SWIFT_COOLING_STRUCT_NONE_H
#define SWIFT_COOLING_STRUCT_NONE_H
/* include grackle */
#include <grackle.h>
/**
* @file src/cooling/none/cooling_struct.h
* @brief Empty infrastructure for the cases without cooling function
......@@ -30,16 +33,22 @@
struct cooling_function_data {
/* Filename of the Cloudy Table */
char GrackleCloudyTable[200];
char cloudy_table[200];
/* Enable/Disable UV backgroud */
int UVbackground;
int uv_background;
/* Redshift to use for the UV backgroud (-1 to use cosmological one) */
double GrackleRedshift;
double redshift;
/* Density Threshold for the shielding */
double GrackleHSShieldingDensityThreshold;
double density_self_shielding;
/* unit system */
code_units units;
/* grackle chemistry data */
chemistry_data chemistry;
};
/**
......
src/cooling/grackle/grackle_wrapper.c deleted 100644 → 0
View file @ 19ca1efd
/***********************************************************************
/
/ Grackle c wrapper
/
/
/ Copyright (c) 2013, Enzo/Grackle Development Team.
/
/ Distributed under the terms of the Enzo Public Licence.
/
/ The full license is in the file LICENSE, distributed with this
/ software.
************************************************************************/
#include "grackle_wrapper.h"
#ifdef SWIFT_DEBUG_CHECKS
#include <assert.h>
#define GRACKLE_ASSERT(X) assert((X))
#else
#define GRACKLE_ASSERT(X)
#endif
code_units my_units;
// arrays passed to grackle as input and to be filled
#define FIELD_SIZE 1
gr_float HDI_density[FIELD_SIZE];
// Set grid dimension and size.
// grid_start and grid_end are used to ignore ghost zones.
const int grid_rank = 1;
int wrap_init_cooling(char *CloudyTable, int UVbackground, double udensity,
double ulength, double utime, int grackle_chemistry) {
#ifdef SWIFT_DEBUG_CHECKS
grackle_verbose = 1;
#endif
double velocity_units;
// First, set up the units system.
// These are conversions from code units to cgs.
my_units.a_units = 1.0; // units for the expansion factor (1/1+zi)
my_units.comoving_coordinates =
0; /* so, according to the doc, we assume here all physical quantities to
be in proper coordinate (not comobile) */
my_units.density_units = udensity;
my_units.length_units = ulength;
my_units.time_units = utime;
velocity_units =
my_units.a_units * my_units.length_units / my_units.time_units;
my_units.velocity_units = velocity_units;
// Second, create a chemistry object for parameters and rate data.
if (set_default_chemistry_parameters() == 0) {
error("Error in set_default_chemistry_parameters.");
}
// Set parameter values for chemistry.
grackle_data.use_grackle = 1;
grackle_data.with_radiative_cooling = 1;
grackle_data.primordial_chemistry =
grackle_chemistry; // molecular network with H, He, D
grackle_data.metal_cooling = 1; // metal cooling on
grackle_data.UVbackground = UVbackground;
grackle_data.grackle_data_file = CloudyTable;
// Finally, initialize the chemistry object.
// snl: a_value is not the true initial ae!! This should get set during
// update_UVbackground
double initial_redshift = 0.;
double a_value = 1. / (1. + initial_redshift);
// Finally, initialize the chemistry object.
if (initialize_chemistry_data(&my_units, a_value) == 0) {
error("Error in initialize_chemistry_data.");
}
return 1;
}
int wrap_set_UVbackground_on() {
// The UV background rates is enabled
grackle_data.UVbackground = 1;
return 1;
}
int wrap_set_UVbackground_off() {
// The UV background rates is disabled
grackle_data.UVbackground = 0;
return 1;
}
int wrap_get_cooling_time(double rho, double u, double Z, double a_now,
double *coolingtime) {
gr_float den_factor = 1.0;
gr_float u_factor = 1.0;
gr_float x_velocity[FIELD_SIZE] = {0.0};
gr_float y_velocity[FIELD_SIZE] = {0.0};
gr_float z_velocity[FIELD_SIZE] = {0.0};
gr_float density[FIELD_SIZE] = {rho * den_factor};
gr_float metal_density[FIELD_SIZE] = {Z * density[0]};
gr_float energy[FIELD_SIZE] = {u * u_factor};
gr_float cooling_time[FIELD_SIZE] = {0.0};
int grid_dimension[3] = {1, 0, 0};
int grid_start[3] = {0, 0, 0};
int grid_end[3] = {0, 0, 0};
if (FIELD_SIZE != 1) {
error("field_size must currently be set to 1.");
}
if (calculate_cooling_time_table(&my_units, a_now, grid_rank, grid_dimension,
grid_start, grid_end, density, energy,
x_velocity, y_velocity, z_velocity,
metal_density, cooling_time) == 0) {
error("Error in calculate_cooling_time.");
}
// return updated chemistry and energy
for (int i = 0; i < FIELD_SIZE; i++) {
coolingtime[i] = cooling_time[i];
}
return 1;
}
int wrap_do_cooling(double rho, double *u, double dt, double Z, double a_now) {
GRACKLE_ASSERT(FIELD_SIZE == 1);
gr_float den_factor = 1.0;
gr_float u_factor = 1.0;
gr_float x_velocity[FIELD_SIZE] = {0.0};
gr_float y_velocity[FIELD_SIZE] = {0.0};
gr_float z_velocity[FIELD_SIZE] = {0.0};
gr_float density[FIELD_SIZE] = {rho * den_factor};
gr_float metal_density[FIELD_SIZE] = {Z * density[0]};
gr_float energy[FIELD_SIZE] = {(*u) * u_factor};
int grid_dimension[3] = {1, 0, 0};
int grid_start[3] = {0, 0, 0};