From 6762eafc5ef4ac90c9da08cac7f5ffb11d4642ae Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 8 Aug 2024 09:10:39 -0400 Subject: [PATCH] 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 */ +