From 54a7cde3fecb7c2c70e7695148f239d34f8181d1 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Tue, 13 Aug 2024 08:59:25 -0400 Subject: [PATCH 1/4] add forgotten call to gr_initialize_field_data into cxx_example.C --- src/example/cxx_example.C | 1 + 1 file changed, 1 insertion(+) diff --git a/src/example/cxx_example.C b/src/example/cxx_example.C index c91f4f20..b7589e8f 100644 --- a/src/example/cxx_example.C +++ b/src/example/cxx_example.C @@ -77,6 +77,7 @@ int main(int argc, char *argv[]) // Create struct for storing grackle field data grackle_field_data my_fields; + gr_initialize_field_data(&my_fields); // Set grid dimension and size. // grid_start and grid_end are used to ignore ghost zones. From d9113f28e98dc356df4aef3b53e0a90b9f382b4f Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 22 Jul 2024 12:29:48 -0400 Subject: [PATCH 2/4] add the current_a_value to the grackle_field_data struct. --- src/clib/set_default_chemistry_parameters.c | 4 ++-- src/include/grackle_fortran_interface.def | 4 +++- src/include/grackle_types.h | 1 + 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/src/clib/set_default_chemistry_parameters.c b/src/clib/set_default_chemistry_parameters.c index 964338c0..d9940394 100644 --- a/src/clib/set_default_chemistry_parameters.c +++ b/src/clib/set_default_chemistry_parameters.c @@ -16,8 +16,7 @@ #include #include "grackle_macros.h" -#include "grackle_types.h" -#include "grackle_chemistry_data.h" +#include "grackle.h" int grackle_verbose = FALSE; @@ -53,6 +52,7 @@ int gr_initialize_field_data(grackle_field_data *my_fields) my_fields->grid_dimension = NULL; my_fields->grid_start = NULL; my_fields->grid_end = NULL; + my_fields->current_a_value = GR_SPECIFY_INITIAL_A_VALUE; my_fields->grid_dx = -1.0; // now, modify all members holding datafields to have values of NULL diff --git a/src/include/grackle_fortran_interface.def b/src/include/grackle_fortran_interface.def index 6629c69e..6f166d24 100644 --- a/src/include/grackle_fortran_interface.def +++ b/src/include/grackle_fortran_interface.def @@ -40,7 +40,9 @@ c This is the fortran definition of grackle_field_data TYPE(C_PTR) :: grid_dimension TYPE(C_PTR) :: grid_start TYPE(C_PTR) :: grid_end + REAL(C_DOUBLE) :: current_a_value REAL(C_DOUBLE) :: grid_dx + TYPE(C_PTR) :: density TYPE(C_PTR) :: HI_density @@ -252,4 +254,4 @@ c The following define the fortran interfaces to the C routines IMPORT TYPE(grackle_field_data), INTENT(INOUT) :: my_fields END FUNCTION gr_initialize_field_data - END INTERFACE \ No newline at end of file + END INTERFACE diff --git a/src/include/grackle_types.h b/src/include/grackle_types.h index 3a5f2d86..f24b71a1 100644 --- a/src/include/grackle_types.h +++ b/src/include/grackle_types.h @@ -50,6 +50,7 @@ typedef struct int *grid_start; int *grid_end; + double current_a_value; gr_float grid_dx; gr_float *density; From 2beb25818370dd1085de82cb65cabc552b8f652f Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 22 Jul 2024 12:52:48 -0400 Subject: [PATCH 3/4] Introduced the unit_handling chemistry parameter. --- src/clib/grackle_chemistry_data_fields.def | 4 ++++ src/clib/initialize_chemistry_data.c | 17 +++++++++++++++-- src/include/grackle.h | 3 +++ src/include/grackle_chemistry_data.h | 3 +++ 4 files changed, 25 insertions(+), 2 deletions(-) diff --git a/src/clib/grackle_chemistry_data_fields.def b/src/clib/grackle_chemistry_data_fields.def index bcc8aa4a..9124b3df 100644 --- a/src/clib/grackle_chemistry_data_fields.def +++ b/src/clib/grackle_chemistry_data_fields.def @@ -208,6 +208,10 @@ ENTRY(collisional_ionisation_rates, INT, 1) //Collisional ionisation ENTRY(recombination_cooling_rates, INT, 1) //Recombination cooling ENTRY(bremsstrahlung_cooling_rates, INT, 1) //Bremsstrahlung cooling +/* Flag to specify unit-handling. Use the GR_UNIT_HANDLING_LEGACY and + GR_UNIT_HANDLING_AUTOMATIC to set the flag values */ +ENTRY(unit_handling, INT, -3) + /* maximum number of subcycle iterations for solve_chemistry */ ENTRY(max_iterations, INT, 10000) diff --git a/src/clib/initialize_chemistry_data.c b/src/clib/initialize_chemistry_data.c index ed51108b..f045f8d7 100644 --- a/src/clib/initialize_chemistry_data.c +++ b/src/clib/initialize_chemistry_data.c @@ -236,6 +236,21 @@ int local_initialize_chemistry_data(chemistry_data *my_chemistry, //omp_set_schedule( omp_sched_auto, chunk_size ); # endif + if (my_chemistry->unit_handling == -3) { + if (grackle_verbose) { + fprintf(stderr, ("WARNING: unit_handling is unset. Defaulting to legacy " + "handling.\n")); + } + my_chemistry->unit_handling = GR_UNIT_HANDLING_LEGACY; + } else if ((my_chemistry->unit_handling != GR_UNIT_HANDLING_LEGACY) && + (my_chemistry->unit_handling != GR_UNIT_HANDLING_AUTOMATIC)) { + fprintf(stderr, "unit_handling has an invalid value\n"); + return GR_FAIL; + } + + /* store a copy of the initial units */ + my_rates->initial_units = *my_units; + /* Only allow a units to be one with proper coordinates. */ if (my_units->comoving_coordinates == FALSE && my_units->a_units != 1.0) { @@ -342,8 +357,6 @@ int local_initialize_chemistry_data(chemistry_data *my_chemistry, return GR_FAIL; } - /* store a copy of the initial units */ - my_rates->initial_units = *my_units; if (grackle_verbose) { time_t timer; diff --git a/src/include/grackle.h b/src/include/grackle.h index 3796e787..14e1e3d5 100644 --- a/src/include/grackle.h +++ b/src/include/grackle.h @@ -26,6 +26,9 @@ extern "C" { #define GR_FAIL 0 #define GR_SPECIFY_INITIAL_A_VALUE -1 +#define GR_UNIT_HANDLING_LEGACY -2 +#define GR_UNIT_HANDLING_AUTOMATIC -1 + extern int grackle_verbose; diff --git a/src/include/grackle_chemistry_data.h b/src/include/grackle_chemistry_data.h index 22ce5007..9a2cc223 100644 --- a/src/include/grackle_chemistry_data.h +++ b/src/include/grackle_chemistry_data.h @@ -176,6 +176,9 @@ typedef struct int recombination_cooling_rates; //Recombination cooling int bremsstrahlung_cooling_rates; //Bremsstrahlung cooling + /* Flag to specify unit-handling */ + int unit_handling; + /* maximum number of subcycle iterations for solve_chemistry */ int max_iterations; From 719059456486162d7a971ddcf0c97127bb424277 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 8 Aug 2024 09:10:39 -0400 Subject: [PATCH 4/4] introduce the proof-of-concept --- src/clib/CMakeLists.txt | 1 + src/clib/calculate_cooling_time.c | 17 +++- src/clib/calculate_dust_temperature.c | 17 +++- src/clib/calculate_gamma.c | 16 +++- src/clib/calculate_pressure.c | 18 +++- src/clib/calculate_temperature.c | 17 +++- src/clib/solve_chemistry.c | 17 +++- src/clib/unit_handling.h | 113 ++++++++++++++++++++++++++ 8 files changed, 194 insertions(+), 22 deletions(-) create mode 100644 src/clib/unit_handling.h diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 4a5772f0..a471ea19 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -129,6 +129,7 @@ add_library(Grackle_Grackle set_default_chemistry_parameters.c solve_chemistry.c update_UVbackground_rates.c + unit_handling.h utils.c # auto-generated C source files diff --git a/src/clib/calculate_cooling_time.c b/src/clib/calculate_cooling_time.c index 7c76d75c..f897c124 100644 --- a/src/clib/calculate_cooling_time.c +++ b/src/clib/calculate_cooling_time.c @@ -14,10 +14,10 @@ #include #include #include +#include "grackle.h" #include "grackle_macros.h" -#include "grackle_types.h" -#include "grackle_chemistry_data.h" #include "phys_constants.h" +#include "unit_handling.h" #include "utils.h" extern chemistry_data *grackle_data; @@ -25,8 +25,6 @@ extern chemistry_data_storage grackle_rates; /* function prototypes */ -double get_temperature_units(code_units *my_units); - int update_UVbackground_rates(chemistry_data *my_chemistry, chemistry_data_storage *my_rates, photo_rate_storage *my_uvb_rates, @@ -90,6 +88,17 @@ int local_calculate_cooling_time(chemistry_data *my_chemistry, if (!my_chemistry->use_grackle) return SUCCESS; + /* do unit-handling */ + code_units units = determine_code_units(my_units, my_rates, + my_fields->current_a_value, + my_chemistry->unit_handling, + "local_calculate_cooling_time"); + if (units.a_units < 0) { + return FAIL; + } else { + my_units = &units; + } + /* Update UV background rates. */ photo_rate_storage my_uvb_rates; diff --git a/src/clib/calculate_dust_temperature.c b/src/clib/calculate_dust_temperature.c index 025da5a6..096fe73b 100644 --- a/src/clib/calculate_dust_temperature.c +++ b/src/clib/calculate_dust_temperature.c @@ -14,10 +14,10 @@ #include #include #include +#include "grackle.h" #include "grackle_macros.h" -#include "grackle_types.h" -#include "grackle_chemistry_data.h" #include "phys_constants.h" +#include "unit_handling.h" #ifdef _OPENMP #include #endif @@ -27,8 +27,6 @@ extern chemistry_data_storage grackle_rates; /* function prototypes */ -double get_temperature_units(code_units *my_units); - int local_calculate_temperature(chemistry_data *my_chemistry, chemistry_data_storage *my_rates, code_units *my_units, @@ -65,6 +63,17 @@ int local_calculate_dust_temperature(chemistry_data *my_chemistry, if (my_chemistry->dust_chemistry < 1 && my_chemistry->h2_on_dust < 1) return SUCCESS; + /* do unit-handling */ + code_units units = determine_code_units(my_units, my_rates, + my_fields->current_a_value, + my_chemistry->unit_handling, + "local_calculate_dust_temperature"); + if (units.a_units < 0) { + return FAIL; + } else { + my_units = &units; + } + double co_length_units, co_density_units; if (my_units->comoving_coordinates == TRUE) { co_length_units = my_units->length_units; diff --git a/src/clib/calculate_gamma.c b/src/clib/calculate_gamma.c index 1b177904..5bda5bd0 100644 --- a/src/clib/calculate_gamma.c +++ b/src/clib/calculate_gamma.c @@ -14,11 +14,11 @@ #include #include #include +#include "grackle.h" #include "grackle_macros.h" -#include "grackle_types.h" -#include "grackle_chemistry_data.h" #include "phys_constants.h" #include "index_helper.h" +#include "unit_handling.h" #ifdef _OPENMP #include #endif @@ -43,6 +43,18 @@ int local_calculate_gamma(chemistry_data *my_chemistry, if (!my_chemistry->use_grackle) return SUCCESS; + + /* do unit-handling */ + code_units units = determine_code_units(my_units, my_rates, + my_fields->current_a_value, + my_chemistry->unit_handling, + "local_calculate_gamma"); + if (units.a_units < 0) { + return FAIL; + } else { + my_units = &units; + } + const grackle_index_helper ind_helper = _build_index_helper(my_fields); int outer_ind, index; diff --git a/src/clib/calculate_pressure.c b/src/clib/calculate_pressure.c index 40f14274..9a2ddd60 100644 --- a/src/clib/calculate_pressure.c +++ b/src/clib/calculate_pressure.c @@ -14,11 +14,11 @@ #include #include #include +#include "grackle.h" #include "grackle_macros.h" -#include "grackle_types.h" -#include "grackle_chemistry_data.h" #include "phys_constants.h" #include "index_helper.h" +#include "unit_handling.h" #ifdef _OPENMP #include #endif @@ -26,8 +26,6 @@ extern chemistry_data *grackle_data; extern chemistry_data_storage grackle_rates; -double get_temperature_units(code_units *my_units); - int local_calculate_pressure(chemistry_data *my_chemistry, chemistry_data_storage *my_rates, code_units *my_units, @@ -38,6 +36,18 @@ int local_calculate_pressure(chemistry_data *my_chemistry, if (!my_chemistry->use_grackle) return SUCCESS; + /* do unit-handling */ + code_units units = determine_code_units(my_units, my_rates, + my_fields->current_a_value, + my_chemistry->unit_handling, + "calculate_pressure"); + if (units.a_units < 0) { + return FAIL; + } else { + my_units = &units; + } + + double tiny_number = 1.e-20; const grackle_index_helper ind_helper = _build_index_helper(my_fields); int outer_ind, index; diff --git a/src/clib/calculate_temperature.c b/src/clib/calculate_temperature.c index 4466db60..61c4abb4 100644 --- a/src/clib/calculate_temperature.c +++ b/src/clib/calculate_temperature.c @@ -14,11 +14,11 @@ #include #include #include +#include "grackle.h" #include "grackle_macros.h" -#include "grackle_types.h" -#include "grackle_chemistry_data.h" #include "phys_constants.h" #include "index_helper.h" +#include "unit_handling.h" #ifdef _OPENMP #include #endif @@ -47,8 +47,6 @@ extern void FORTRAN_NAME(calc_temp_cloudy_g)( double *priPar1, double *priPar2, double *priPar3, long long *priDataSize, double *priMMW); -double get_temperature_units(code_units *my_units); - int local_calculate_pressure(chemistry_data *my_chemistry, chemistry_data_storage *my_rates, code_units *my_units, @@ -71,6 +69,17 @@ int local_calculate_temperature(chemistry_data *my_chemistry, if (!my_chemistry->use_grackle) return SUCCESS; + /* do unit-handling */ + code_units units = determine_code_units(my_units, my_rates, + my_fields->current_a_value, + my_chemistry->unit_handling, + "local_solve_chemistry"); + if (units.a_units < 0) { + return FAIL; + } else { + my_units = &units; + } + /* Compute the pressure first. */ if (my_chemistry->primordial_chemistry > 0) { diff --git a/src/clib/solve_chemistry.c b/src/clib/solve_chemistry.c index d4398af0..23d0dbe5 100644 --- a/src/clib/solve_chemistry.c +++ b/src/clib/solve_chemistry.c @@ -13,10 +13,10 @@ #include #include +#include "grackle.h" #include "grackle_macros.h" -#include "grackle_types.h" -#include "grackle_chemistry_data.h" #include "phys_constants.h" +#include "unit_handling.h" #include "utils.h" extern chemistry_data *grackle_data; @@ -24,8 +24,6 @@ extern chemistry_data_storage grackle_rates; /* function prototypes */ -double get_temperature_units(code_units *my_units); - int update_UVbackground_rates(chemistry_data *my_chemistry, chemistry_data_storage *my_rates, photo_rate_storage *my_uvb_rates, @@ -102,6 +100,17 @@ int local_solve_chemistry(chemistry_data *my_chemistry, if (!my_chemistry->use_grackle) return SUCCESS; + /* do unit-handling */ + code_units units = determine_code_units(my_units, my_rates, + my_fields->current_a_value, + my_chemistry->unit_handling, + "local_solve_chemistry"); + if (units.a_units < 0) { + return FAIL; + } else { + my_units = &units; + } + /* Update UV background rates. */ photo_rate_storage my_uvb_rates; diff --git a/src/clib/unit_handling.h b/src/clib/unit_handling.h new file mode 100644 index 00000000..05fa532d --- /dev/null +++ b/src/clib/unit_handling.h @@ -0,0 +1,113 @@ +/*********************************************************************** +/ +/ Implement and declare automatic unit-handling +/ +/ 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. +************************************************************************/ + +#ifndef UNIT_HANDLING_H +#define UNIT_HANDLING_H + +#include // provides fprintf +#include "grackle.h" + +#ifdef __cplusplus +extern "C" { +#endif /* __cplusplus */ + + + +/// Automatically compute a code_units instance when Grackle is configured +/// with automatic code-units. +/// +/// @param specified_units Pointer to the specified units +/// @param my_rates Pointer to the chemistry_data_storage object that contains +/// the initial rates that were specified when the grackle_solver was +/// initialized +/// @param current_a_value The specified current scale-factor +/// @param unit_handling The value of Grackle's unit handling flag +/// @param func_name Name of the function that is calling the error check +/// +/// @returns a fully initialized code_units struct. If there was any kind of +/// issue, then the a_units member will be negative. +/// +/// @note +/// This is declared inline to reduce as much overhead as possible +/// +/// @note +/// This is currently a proof-of-concept, some improvements are definitely +/// possible. For example, we could: +/// * directly pass in my_rates->initial_units rather than my_rates +/// * refactor gr_query_units so that we could directly call into its +/// internals (the error-checking and string comparisons introduce a bunch +/// of unnecessary overhead) +static inline code_units determine_code_units ( + const code_units* specified_units, + const chemistry_data_storage* my_rates, + double current_a_value, + int unit_handling, + const char* func_name) +{ + code_units out; + out.a_units = -1.0; + + if (unit_handling == GR_UNIT_HANDLING_LEGACY) { + // handle error-cases + if (current_a_value != GR_SPECIFY_INITIAL_A_VALUE) { + fprintf(stderr, + "ERROR in %s: current_a_value shouldn't be specified when using " + "Grackle's legacy unit handling\n", func_name); + out.a_units = -1.0; + return out; + } else if (specified_units == NULL) { + fprintf(stderr, + "ERROR in %s: code_units argument can't be NULL when using " + "Grackle's legacy unit handling\n", func_name); + return out; + } + + // this is the happy path + out = *specified_units; + return out; + + } else if (unit_handling == GR_UNIT_HANDLING_AUTOMATIC) { + // handle error-case + if (specified_units != NULL) { + fprintf(stderr, + "ERROR in %s: code_units argument must be NULL when using Grackle's " + "automatic unit handling\n", func_name); + return out; + } + + // this is the happy path + out = my_rates->initial_units; + out.a_value = current_a_value; + out.density_units = gr_query_units(my_rates, "density_units", + current_a_value); + out.length_units = gr_query_units(my_rates, "length_units", + current_a_value); + out.time_units = gr_query_units(my_rates, "time_units", + current_a_value); + set_velocity_units(&out); + return out; + + } else { + fprintf(stderr, + "ERROR in %s: Grackle's unit_handling parameter has an unrecognized " + "value\n", func_name); + return out; + } + +} + +#ifdef __cplusplus +} /* extern "C" */ +#endif /* __cplusplus */ + +#endif /* UNIT_HANDLING_H */ +