diff --git a/src/cooling/grackle/cooling.h b/src/cooling/grackle/cooling.h
index 50700071904be662609e513b34a0215d0ef61d7f..925244b668517ff80bc850d695cfaec60be6dc27 100644
--- a/src/cooling/grackle/cooling.h
+++ b/src/cooling/grackle/cooling.h
@@ -27,6 +27,7 @@
/* Some standard headers. */
#include <float.h>
#include <math.h>
+#include <grackle.h>
/* Local includes. */
#include "error.h"
@@ -36,13 +37,37 @@
#include "physical_constants.h"
#include "units.h"
-/* include the grackle wrapper */
-#include "grackle_wrapper.h"
+/* need to rework code if changed */
+#define GRACKLE_NPART 1
+
/**
- * @brief Compute the cooling rate
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state.
*
- * We do nothing.
+ * @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 Compute the cooling rate
*
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
@@ -134,30 +159,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.
*
@@ -171,27 +172,67 @@ static INLINE void cooling_init_backend(
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)
+
+ /* We assume here all physical quantities to
+ be in proper coordinate (not comobile) */
+ cooling->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;
+
+ /* Create a chemistry object for parameters and rate data. */
+ if (set_default_chemistry_parameters() == 0) {
+ error("Error in set_default_chemistry_parameters.");
+ }
- units_density = us->UnitMass_in_cgs / pow(us->UnitLength_in_cgs, 3);
- units_length = us->UnitLength_in_cgs;
- units_time = us->UnitTime_in_cgs;
+ // Set parameter values for chemistry.
+ grackle_data.use_grackle = 1;
+ grackle_data.with_radiative_cooling = 1;
+ /* molecular network with H, He, D
+ From Cloudy table */
+ grackle_data.primordial_chemistry = 0;
+ grackle_data.metal_cooling = 1; // metal cooling on
+ grackle_data.UVbackground = cooling->uv_background;
+ grackle_data.grackle_data_file = cooling->cloudy_table;
+
+ /* Initialize the chemistry object.
+ a_value is not the true initial a
+ This should get set before any computation */
+ double a_value = 1.;
+
+ if (initialize_chemistry_data(&cooling->units, a_value) == 0) {
+ error("Error in initialize_chemistry_data.");
+ }
#ifdef SWIFT_DEBUG_CHECKS
+ if (GRACKLE_NPART != 1)
+ error("Grackle with multiple particles not implemented");
float threshold = cooling->GrackleHSShieldingDensityThreshold;
threshold /= phys_const->const_proton_mass;
@@ -200,21 +241,11 @@ 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();
+ //grackle_print_data();
message("");
message("***************************************");
#endif
@@ -229,12 +260,90 @@ 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);
+ 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 = %g", 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 print data in grackle struct
+ *
+ * Should only be used for debugging
+ */
+void grackle_print_data() {
+ message("Grackle Data:");
+ message("\t Data file: %s", grackle_data.grackle_data_file);
+ message("\t With grackle: %i", grackle_data.use_grackle);
+ message("\t With radiative cooling: %i", grackle_data.with_radiative_cooling);
+ message("\t With UV background: %i", grackle_data.UVbackground);
+ message("\t With primordial chemistry: %i",
+ grackle_data.primordial_chemistry);
+ message("\t Number temperature bins: %i",
+ grackle_data.NumberOfTemperatureBins);
+ message("\t T = (%g, ..., %g)", grackle_data.TemperatureStart,
+ grackle_data.TemperatureEnd);
+
+ message("Primordial Cloudy");
+ cloudy_print_data(grackle_data.cloudy_primordial, 1);
+ if (grackle_data.metal_cooling) {
+ message("Metal Cooling");
+ cloudy_print_data(grackle_data.cloudy_metal, 0);
+ }
+
+ message("\t Gamma: %g", grackle_data.Gamma);
+
+ /* UVB */
+ if (grackle_data.UVbackground && grackle_data.primordial_chemistry != 0) {
+ struct UVBtable uvb = grackle_data.UVbackground_table;
+ long long N = uvb.Nz;
+ message("\t UV Background");
+ message("\t\t Redshift from %g to %g with %lli steps", uvb.zmin, uvb.zmax,
+ N);
+ message("\t\t z = (%g, ..., %g)", uvb.z[0], uvb.z[N - 1]);
+ }
+}
+
+/**
+ * @brief print data in cloudy struct
+ *
+ * Should only be used for debugging
+ */
+void cloudy_print_data(const cloudy_data c, const int print_mmw) {
+ long long N = c.data_size;
+ message("\t Data size: %lli", N);
+ message("\t Grid rank: %lli", c.grid_rank);
+
+ char msg[200] = "\t Dimension: (";
+ for (long long i = 0; i < c.grid_rank; i++) {
+ char tmp[200] = "%lli%s";
+ if (i == c.grid_rank - 1)
+ sprintf(tmp, tmp, c.grid_dimension[i], ")");
+ else
+ sprintf(tmp, tmp, c.grid_dimension[i], ", ");
+
+ strcat(msg, tmp);
+ }
+ message("%s", msg);
+
+ if (c.heating_data)
+ message("\t Heating: (%g, ..., %g)", c.heating_data[0],
+ c.heating_data[N - 1]);
+ if (c.cooling_data)
+ message("\t Cooling: (%g, ..., %g)", c.cooling_data[0],
+ c.cooling_data[N - 1]);
+ if (c.mmw_data && print_mmw)
+ message("\t Mean molecular weigth: (%g, ..., %g)", c.mmw_data[0],
+ c.mmw_data[N - 1]);
}
#endif /* SWIFT_COOLING_GRACKLE_H */
diff --git a/src/cooling/grackle/cooling_struct.h b/src/cooling/grackle/cooling_struct.h
index b59af375cba0ed35fcb0159cb63605af8d94f030..4dcb7c415c8bbf04faffeb130d981f98547c62ec 100644
--- a/src/cooling/grackle/cooling_struct.h
+++ b/src/cooling/grackle/cooling_struct.h
@@ -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,20 @@
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;
+
};
/**
diff --git a/src/cooling/grackle/grackle_wrapper.c b/src/cooling/grackle/grackle_wrapper.c
index 9974717b670c937811cea54349124d9a38bed476..f1aa79dd32e0789c1e61ed16c03eddff0679860b 100644
--- a/src/cooling/grackle/grackle_wrapper.c
+++ b/src/cooling/grackle/grackle_wrapper.c
@@ -82,18 +82,6 @@ int wrap_init_cooling(char *CloudyTable, int UVbackground, double udensity,
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;
@@ -165,63 +153,3 @@ int wrap_do_cooling(double rho, double *u, double dt, double Z, double a_now) {
return 1;
}
-void grackle_print_data() {
- message("Grackle Data:");
- message("\t Data file: %s", grackle_data.grackle_data_file);
- message("\t With grackle: %i", grackle_data.use_grackle);
- message("\t With radiative cooling: %i", grackle_data.with_radiative_cooling);
- message("\t With UV background: %i", grackle_data.UVbackground);
- message("\t With primordial chemistry: %i",
- grackle_data.primordial_chemistry);
- message("\t Number temperature bins: %i",
- grackle_data.NumberOfTemperatureBins);
- message("\t T = (%g, ..., %g)", grackle_data.TemperatureStart,
- grackle_data.TemperatureEnd);
-
- message("Primordial Cloudy");
- cloudy_print_data(grackle_data.cloudy_primordial, 1);
- if (grackle_data.metal_cooling) {
- message("Metal Cooling");
- cloudy_print_data(grackle_data.cloudy_metal, 0);
- }
-
- message("\t Gamma: %g", grackle_data.Gamma);
-
- /* UVB */
- if (grackle_data.UVbackground && grackle_data.primordial_chemistry != 0) {
- struct UVBtable uvb = grackle_data.UVbackground_table;
- long long N = uvb.Nz;
- message("\t UV Background");
- message("\t\t Redshift from %g to %g with %lli steps", uvb.zmin, uvb.zmax,
- N);
- message("\t\t z = (%g, ..., %g)", uvb.z[0], uvb.z[N - 1]);
- }
-}
-
-void cloudy_print_data(const cloudy_data c, const int print_mmw) {
- long long N = c.data_size;
- message("\t Data size: %lli", N);
- message("\t Grid rank: %lli", c.grid_rank);
-
- char msg[200] = "\t Dimension: (";
- for (long long i = 0; i < c.grid_rank; i++) {
- char tmp[200] = "%lli%s";
- if (i == c.grid_rank - 1)
- sprintf(tmp, tmp, c.grid_dimension[i], ")");
- else
- sprintf(tmp, tmp, c.grid_dimension[i], ", ");
-
- strcat(msg, tmp);
- }
- message("%s", msg);
-
- if (c.heating_data)
- message("\t Heating: (%g, ..., %g)", c.heating_data[0],
- c.heating_data[N - 1]);
- if (c.cooling_data)
- message("\t Cooling: (%g, ..., %g)", c.cooling_data[0],
- c.cooling_data[N - 1]);
- if (c.mmw_data && print_mmw)
- message("\t Mean molecular weigth: (%g, ..., %g)", c.mmw_data[0],
- c.mmw_data[N - 1]);
-}