From 12868c64a1672d715fd63a12f0e1536f3700a7f0 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sun, 14 Sep 2025 11:05:23 -0400 Subject: [PATCH 01/11] introduce GrainSpeciesInfo --- src/clib/CMakeLists.txt | 1 + src/clib/LUT.hpp | 4 +- src/clib/dust/grain_species_info.hpp | 349 +++++++++++++++++++++++++ src/clib/dust_props.hpp | 2 + src/clib/initialize_chemistry_data.cpp | 2 + src/clib/initialize_rates.cpp | 13 +- src/clib/lookup_cool_rates1d.hpp | 1 + src/clib/opaque_storage.hpp | 11 + tests/unit/CMakeLists.txt | 4 + tests/unit/test_grain_species_info.cpp | 259 ++++++++++++++++++ 10 files changed, 643 insertions(+), 3 deletions(-) create mode 100644 src/clib/dust/grain_species_info.hpp create mode 100644 tests/unit/test_grain_species_info.cpp diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 7f071d562..e3ea52b95 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -110,6 +110,7 @@ add_library(Grackle_Grackle collisional_rate_props.cpp collisional_rate_props.hpp cool_multi_time_g.cpp cool_multi_time_g.h dust_props.hpp + dust/grain_species_info.hpp init_misc_species_cool_rates.cpp init_misc_species_cool_rates.hpp initialize_chemistry_data.cpp initialize_dust_yields.cpp initialize_dust_yields.hpp diff --git a/src/clib/LUT.hpp b/src/clib/LUT.hpp index 91be1f901..04e5706b3 100644 --- a/src/clib/LUT.hpp +++ b/src/clib/LUT.hpp @@ -99,12 +99,12 @@ struct OnlyGrainSpLUT { // XMacros provided in grackle_field_data_fdatamembers.def (or we may need to // slightly revise the system?) enum { + MgSiO3_dust, + AC_dust, SiM_dust, FeM_dust, Mg2SiO4_dust, - MgSiO3_dust, Fe3O4_dust, - AC_dust, SiO2_dust, MgO_dust, FeS_dust, diff --git a/src/clib/dust/grain_species_info.hpp b/src/clib/dust/grain_species_info.hpp new file mode 100644 index 000000000..d6eb14ff7 --- /dev/null +++ b/src/clib/dust/grain_species_info.hpp @@ -0,0 +1,349 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Defines details about the dust grain species +/// +//===----------------------------------------------------------------------===// + +#ifndef GRAIN_SPECIES_INFO_HPP +#define GRAIN_SPECIES_INFO_HPP + +#include "LUT.hpp" +#include "grackle_macros.h" +#include // memcpy + +namespace grackle::impl { + +/// holds information about a single gas species that is an ingredient for +/// grain growth +struct GrainGrowthIngredient { + /// the positive integer coefficient associated with the ingredient + int coef; + + /// species index of the ingredient + int species_idx; + + /// the particle mass of the ingredient species (in atomic mass units) + /// + /// @note + /// It's plausible that this should actually be the particle mass in hydrogen + /// masses, which is slightly different from atomic mass units. + /// Unfortunately, the original code didn't specify constants with enough + /// precision to tell the difference + double mparticle_amu; +}; + +/// summarizes details about a single grain species +struct GrainSpeciesInfoEntry { + /// the species index of the grain in the #GrainSpLUT lookup table + int species_idx; + + /// the species index of the grain in the #OnlyGrainSpLUT lookup table + /// + /// @note + /// It's frankly a little redundant to track this information (since an + /// instance of this struct is found at this index of an out.species_infoay) + int onlygrainsp_idx; + + /// name of the dust species + /// + /// @note + /// This primarily exists for debuging purposes + const char* name; + + /// indicates whether to use the carbonaceous or silicate coefficient table + /// to computing contributions of the grain species to the total h2dust rate + /// (or the rate of H2 formation) + /// + /// @note + /// The choice to encode this information within each GrainSpeciesInfoEntry + /// may be somewhat inefficient since (at the time of writing) Amorphous + /// Carbon grains are the **ONLY** grain species that use the carbonaceous + /// coefficient table & all other species use the silicate coefficient table + /// (to be honest, the correctness of this seems questionable, at best). As + /// soon as more than one grain uses the carbonaceous coefficient table (or + /// some other non sillicate table), it will be clearly be much more optimal + /// to track the information in this manner. + bool h2dust_uses_carbonaceous_table; + + /// The sublimation temperature + double sublimation_temperature; + + /// The number of growth ingredients + int n_growth_ingredients; + + /// The array of growth ingredients + const GrainGrowthIngredient* growth_ingredients; +}; + +/// tracks for details about every relevant grain species for a given Grackle +/// configuration. This includes information necessary for computing growth & +/// destruction rates +/// +/// Relationship with OnlyGrainSpLUT +/// -------------------------------- +/// In the short term, the index of each species in the +/// @ref GrainSpeciesInfo::species_info out.species_infoay is dictated by the +/// order of enumerators in the OnlyGrainSpLUT enumeration. +/// +/// In the medium term, we plan to entirely eliminate the OnlyGrainSpLUT +/// enumeration because all of the grain species can be treated very uniformly +/// uniformly. At the time of writing, just about every place where we would +/// use OnlyGrainSpLUT corresponds to a location where would enumerate every +/// possible grain species and perform nearly identical operations on each +/// species. In each case, it is straight-forward to replace these blocks of +/// logic with for-loops (we just need to encode species-specific variations in +/// the calculations in out.species_infoays that have the same ordering as the +/// species). To phrase it another way, in nearly all of the places where we +/// would use OnlyGrainSpLUT, we don't need to know the grain species identity. +/// +/// The exception to this is when we compute the h2dust rate. In this case +/// we need to identify AC_dust since we need to do something **slightly** +/// different from the other grains, but this is easy to work around +struct GrainSpeciesInfo { + /// number of grain species considered for the current Grackle configuration + int n_species; + + /// an out.species_infoay of length of length @ref n_species where each entry + /// holds info about a separate grain species + GrainSpeciesInfoEntry* species_info; +}; + +/// return the number of grain species +/// +/// @param[in] dust_species_parameter The parameter tracked by #chemistry_data +/// @return The number of grain species. A negative value indicates that the +/// @p dust_species_parameter has an invalid value +inline int get_n_grain_species(int dust_species_parameter) { + switch (dust_species_parameter) { + case 0: + return 0; + case 1: + return 2; + case 2: + return 10; + case 3: + return 13; + default: + return -1; + } +} + +/// Species the maximum number of ingredients that a grain species has +/// +/// @note +/// The correctness of this constant is explicitly checked in a unit test +inline constexpr int max_ingredients_per_grain_species = 3; + +// define a few macros to help implement new_GrainSpeciesInfo +// - we undef all of these macros after they're used for cleanliness +// - we may want to consider an alternative that avoids macros in the future +#define GRIMPL_INGREDIENT_LIST_SENTINEL {-1, -1, -1.0} + +// does the heavy-lifting for implementing GR_MK_GRAIN_INFO_ENTRY +inline GrainSpeciesInfoEntry mk_gsp_info_entry_helper_( + int species_idx, int onlygrainsp_idx, const char* name, + bool h2dust_uses_carbonaceous_table, double sublimation_temperature, + const GrainGrowthIngredient* growth_ingredients) { + GrainGrowthIngredient* out_ingredient_ptr = nullptr; + int n_ingredients = 0; + + if (growth_ingredients != nullptr) { + n_ingredients = 0; + // increment the counter until we reach the sentinel + while (growth_ingredients[n_ingredients].coef != -1) { + n_ingredients++; + } + // allocate and initialize the pointer + out_ingredient_ptr = new GrainGrowthIngredient[n_ingredients]; + std::memcpy(out_ingredient_ptr, growth_ingredients, + (std::size_t)n_ingredients); + } + + return GrainSpeciesInfoEntry{species_idx, + onlygrainsp_idx, + name, + h2dust_uses_carbonaceous_table, + sublimation_temperature, + n_ingredients, + out_ingredient_ptr}; +} + +/// fills the specified out.species_infoay with the number of grain species. +/// +/// The out.species_infoay is expected to have space for the number of entries +/// returned by get_n_grain_species(dust_species_paramter) +/// +/// @param[in] dust_species_parameter The parameter tracked by #chemistry_data +/// @returns A fully initialized GrainSpeciesInfo instance +inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { + GrainSpeciesInfo out{-1, nullptr}; // indicates an error + out.n_species = get_n_grain_species(dust_species_parameter); + + if (out.n_species <= 0) { + return out; + } + + out.species_info = new GrainSpeciesInfoEntry[out.n_species]; + + // At the time of writing: + // - we **only** use h2rate_carbonaceous_coef_table for the AC_dust + // species (amorphous carbon grains) + // - we use h2rate_silicate_coef_table for **ALL** other grain species + // (including species that are obviously NOT silicates) + // + // It is not obvious at all why this is the case... + // + // TODO: we should add documentation clearly addressing each of the + // following questions and replace this todo-item with a reference to the + // part of the documentation that discusses this! The documentation + // should make sure to address: + // - if this choice physically motivated or a common convention? (or if + // it was a quick & dirty choice because Gen didn't know what to do) + // - do we have a sense for how "wrong" the choice is? (e.g. will it + // generally overpredict/underpredict?) + // - why don't we use the generic my_rates->h2dusta rate as a generic + // fallback for non-silicate grains instead of using h2dustS? I realize + // that h2dusta has very different units... + + if (dust_species_parameter > 0) { + const GrainGrowthIngredient MgSiO3_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::Mg, 24.}, + {1, SpLUT::SiOI, 44.}, + {2, SpLUT::H2O, 18.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[0] = mk_gsp_info_entry_helper_( + SpLUT::MgSiO3_dust, OnlyGrainSpLUT::MgSiO3_dust, "MgSiO3_dust", false, + 1222.0, MgSiO3_dust_ingred); + + const GrainGrowthIngredient AC_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::CI, 12.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[1] = + mk_gsp_info_entry_helper_(SpLUT::AC_dust, OnlyGrainSpLUT::AC_dust, + "AC_dust", true, 1800.0, AC_dust_ingred); + } + + if (dust_species_parameter > 1) { + const GrainGrowthIngredient SiM_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::SiI, 28.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[2] = + mk_gsp_info_entry_helper_(SpLUT::SiM_dust, OnlyGrainSpLUT::SiM_dust, + "SiM_dust", false, 1500.0, SiM_dust_ingred); + + const GrainGrowthIngredient FeM_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::Fe, 56.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[3] = + mk_gsp_info_entry_helper_(SpLUT::FeM_dust, OnlyGrainSpLUT::FeM_dust, + "FeM_dust", false, 1500.0, FeM_dust_ingred); + + const GrainGrowthIngredient Mg2SiO4_dust_ingred[] = { + // {coef, species_idx, particle mass} + {2, SpLUT::Mg, 24.}, + {1, SpLUT::SiOI, 44.}, + {3, SpLUT::H2O, 18.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[4] = mk_gsp_info_entry_helper_( + SpLUT::Mg2SiO4_dust, OnlyGrainSpLUT::Mg2SiO4_dust, "Mg2SiO4_dust", + false, 1277.0, Mg2SiO4_dust_ingred); + + const GrainGrowthIngredient Fe3O4_dust_ingred[] = { + // {coef, species_idx, particle mass} + {3, SpLUT::Fe, 56.}, + {4, SpLUT::H2O, 18.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[5] = mk_gsp_info_entry_helper_( + SpLUT::Fe3O4_dust, OnlyGrainSpLUT::Fe3O4_dust, "Fe3O4_dust", false, + 1500.0, Fe3O4_dust_ingred); + + const GrainGrowthIngredient SiO2_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::SiO2I, 60.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[6] = + mk_gsp_info_entry_helper_(SpLUT::SiO2_dust, OnlyGrainSpLUT::SiO2_dust, + "SiO2_dust", false, 1500.0, SiO2_dust_ingred); + + const GrainGrowthIngredient MgO_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::Mg, 24.}, + {1, SpLUT::H2O, 18.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[7] = + mk_gsp_info_entry_helper_(SpLUT::MgO_dust, OnlyGrainSpLUT::MgO_dust, + "MgO_dust", false, 1500.0, MgO_dust_ingred); + + const GrainGrowthIngredient FeS_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::Fe, 56.}, + {1, SpLUT::S, 32.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[8] = + mk_gsp_info_entry_helper_(SpLUT::FeS_dust, OnlyGrainSpLUT::FeS_dust, + "FeS_dust", false, 680.0, FeS_dust_ingred); + + const GrainGrowthIngredient Al2O3_dust_ingred[] = { + // {coef, species_idx, particle mass} + {2, SpLUT::Al, 27.}, + {3, SpLUT::H2O, 18.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[9] = mk_gsp_info_entry_helper_( + SpLUT::Al2O3_dust, OnlyGrainSpLUT::Al2O3_dust, "Al2O3_dust", false, + 1500.0, Al2O3_dust_ingred); + } + + if (dust_species_parameter > 2) { + // We do not consider the growth of refractory organics, volatile + // organics, and water ice because their sublimation temperatures + // are low (100-600 K). They sublimate before the growth occurs. + + out.species_info[10] = mk_gsp_info_entry_helper_( + SpLUT::ref_org_dust, OnlyGrainSpLUT::ref_org_dust, "ref_org_dust", + false, 575.0, nullptr); + out.species_info[11] = mk_gsp_info_entry_helper_( + SpLUT::vol_org_dust, OnlyGrainSpLUT::vol_org_dust, "vol_org_dust", + false, 375.0, nullptr); + out.species_info[12] = mk_gsp_info_entry_helper_( + SpLUT::H2O_ice_dust, OnlyGrainSpLUT::H2O_ice_dust, "H2O_ice_dust", + false, 153.0, nullptr); + } + + return out; +} + +#undef GRIMPL_INGREDIENT_LIST_SENTINEL + +/// performs cleanup of the contents of GrainSpeciesInfo +/// +/// This effectively invokes the destructor +inline void drop_GrainSpeciesInfo(GrainSpeciesInfo* ptr) { + if (ptr->n_species == 0) { + return; // avoids double-free + } + + for (int gsp_idx = 0; gsp_idx < ptr->n_species; gsp_idx++) { + if ((ptr->species_info[gsp_idx].growth_ingredients) != nullptr) { + delete[] ptr->species_info[gsp_idx].growth_ingredients; + } + } + delete[] ptr->species_info; + // the following 2 lines are not strictly necessary, but they may help us + // avoid a double-free and a dangling pointer + ptr->n_species = 0; + ptr->species_info = nullptr; +} + +} // namespace grackle::impl + +#endif /* GRAIN_SPECIES_INFO_HPP */ diff --git a/src/clib/dust_props.hpp b/src/clib/dust_props.hpp index 331cd9f4e..c1daa9c4b 100644 --- a/src/clib/dust_props.hpp +++ b/src/clib/dust_props.hpp @@ -19,6 +19,7 @@ #include // malloc, free #include "internal_types.hpp" +#include "LUT.hpp" namespace grackle::impl { @@ -125,6 +126,7 @@ inline void drop_InternalDustPropBuf(InternalDustPropBuf* ptr) { std::free(ptr->dyntab_kappa_tot); } + } // namespace grackle::impl #endif // DUST_PROPS_HPP diff --git a/src/clib/initialize_chemistry_data.cpp b/src/clib/initialize_chemistry_data.cpp index 13d5a0ad3..38a15963f 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -403,6 +403,7 @@ extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry, my_rates->opaque_storage->n_kcol_rate_indices = 0; init_empty_interp_grid_props_( &my_rates->opaque_storage->h2dust_grain_interp_props); + my_rates->opaque_storage->grain_species_info = NULL; double co_length_units, co_density_units; if (my_units->comoving_coordinates == TRUE) { @@ -636,6 +637,7 @@ extern "C" int local_free_chemistry_data(chemistry_data *my_chemistry, } GRACKLE_FREE(my_rates->opaque_storage->used_kcol_rate_indices); free_interp_grid_props_(&my_rates->opaque_storage->h2dust_grain_interp_props); + GRACKLE_FREE(my_rates->opaque_storage->grain_species_info); GRACKLE_FREE(my_rates->opaque_storage); diff --git a/src/clib/initialize_rates.cpp b/src/clib/initialize_rates.cpp index 95b9a594a..8c9f38c0f 100644 --- a/src/clib/initialize_rates.cpp +++ b/src/clib/initialize_rates.cpp @@ -45,6 +45,7 @@ #include "grackle_macros.h" #include "grackle_rate_functions.h" #include "collisional_rate_props.hpp" // init_extra_collisional_rates +#include "dust/grain_species_info.hpp" #include "init_misc_species_cool_rates.hpp" // init_misc_species_cool_rates #include "initialize_dust_yields.hpp" // initialize_dust_yields #include "initialize_rates.hpp" @@ -645,7 +646,17 @@ int grackle::impl::initialize_rates( fprintf(stderr, "Error in initialize_metal_chemistry_rates.\n"); return GR_FAIL; } - /* Dust rates */ + + // Dust Grain Species Information + // (it may make sense want to handle more of the dust separately) + if (my_chemistry->dust_species > 0) { + my_rates->opaque_storage->grain_species_info = + (grackle::impl::GrainSpeciesInfo*)malloc(sizeof(grackle::impl::GrainSpeciesInfo)); + *(my_rates->opaque_storage->grain_species_info) = + grackle::impl::new_GrainSpeciesInfo(my_chemistry->dust_species); + } + + // Dust rates if (grackle::impl::initialize_dust_yields(my_chemistry, my_rates, my_units) == FAIL) { fprintf(stderr, "Error in initialize_dust_yields.\n"); return FAIL; diff --git a/src/clib/lookup_cool_rates1d.hpp b/src/clib/lookup_cool_rates1d.hpp index 237a90b6f..299f8e1f9 100644 --- a/src/clib/lookup_cool_rates1d.hpp +++ b/src/clib/lookup_cool_rates1d.hpp @@ -20,6 +20,7 @@ #include "grackle.h" #include "dust_props.hpp" +#include "dust/grain_species_info.hpp" #include "fortran_func_decls.h" #include "fortran_func_wrappers.hpp" #include "internal_types.hpp" diff --git a/src/clib/opaque_storage.hpp b/src/clib/opaque_storage.hpp index 883bc3b34..a118d42c4 100644 --- a/src/clib/opaque_storage.hpp +++ b/src/clib/opaque_storage.hpp @@ -14,6 +14,7 @@ #define OPAQUE_STORAGE_HPP #include "grackle.h" +#include "dust/grain_species_info.hpp" #include "internal_types.hpp" /// a struct that used to wrap some private storage details @@ -84,6 +85,16 @@ struct gr_opaque_storage { /// explicitly tracking a grid of ln(Tgas) values and a grid of ln(Tdust) /// values (for debugging purposes). gr_interp_grid_props h2dust_grain_interp_props; + + /// Tracks basic information about each relevant grain species + /// + /// > [!note] + /// > A case could be made that we may not want to directly store a + /// > GrainSpeciesInfo instances directly inside of this struct (since it + /// > contains some extra information that is unnecessary during the + /// > calculations). An alternative would be to briefly initialize an + /// > instance during setup and then repack the data. + grackle::impl::GrainSpeciesInfo* grain_species_info; }; #endif /* OPAQUE_STORAGE_HPP */ diff --git a/tests/unit/CMakeLists.txt b/tests/unit/CMakeLists.txt index 9e05050ba..f1018e3ca 100644 --- a/tests/unit/CMakeLists.txt +++ b/tests/unit/CMakeLists.txt @@ -57,6 +57,10 @@ add_executable(runLinAlgTests test_linalg.cpp) target_link_libraries(runLinAlgTests testdeps) gtest_discover_tests(runLinAlgTests) +add_executable(runGrainSpeciesIntoTests test_grain_species_info.cpp) +target_link_libraries(runGrainSpeciesIntoTests testdeps) +gtest_discover_tests(runGrainSpeciesIntoTests) + add_executable(runStatusReporting test_status_reporting.cpp) target_link_libraries(runStatusReporting grtest_utils) gtest_discover_tests(runStatusReporting) diff --git a/tests/unit/test_grain_species_info.cpp b/tests/unit/test_grain_species_info.cpp new file mode 100644 index 000000000..19bd87f4e --- /dev/null +++ b/tests/unit/test_grain_species_info.cpp @@ -0,0 +1,259 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Tests grackle::impl::GrainSpeciesInfo +//===----------------------------------------------------------------------===// + +#include +#include +#include +#include +#include + +#include "LUT.hpp" +#include "dust/grain_species_info.hpp" + +namespace { // stuff in an anonymous namespace is local to this file + +/// checks that the string, @p s, ends with the suffix, @p suffix +/// +/// @param[in] s the string to check +/// @param[in] suffix the desired behavior +bool str_ends_with(std::string_view s, std::string_view suffix) { +#if __cplusplus >= 202002L + return s.ends_with(suffix); +#else + // according to cppreference, this is the behavior of the ends_with method + // introduced in C++20 + std::size_t s_len = s.size(); + std::size_t suf_len = suffix.size(); + return ((s_len >= suf_len) && + (s.compare(s_len - suf_len, std::string_view::npos, suffix) == 0)); +#endif +} + +/// Represents a species name and its location within the global species list +struct SpNameIndexPair { + std::string name; + int index; +}; + +/// Returns name-index pairs for all grains from the list of Grackle Species +/// +/// The index corresponds to the index in global species list +std::vector grain_sp_info_from_species_list() { + std::vector out; + int count = 0; + +#define STRINGIFY_(NAME) #NAME +#define ENTRY(NAME) \ + if (str_ends_with(STRINGIFY_(NAME), "_dust")) { \ + out.push_back(SpNameIndexPair{STRINGIFY_(NAME), count}); \ + } \ + count++; +#include "field_data_evolved_species.def" +#undef ENTRY +#undef STRINGIFY_ + + return out; +} + +/// returns the number of grain species known to grackle +int number_known_grain_species() { + // this is a very silly implementation + return (int)(grain_sp_info_from_species_list().size()); +} + +/// tracks the max value allowed by Grackle's dust_species runtime parameter +constexpr int MAX_dust_species_VAL = 3; + +// down below, we define some machinery to make sure we properly deallocate +// memory when a test fails. This would not be an issue if we were willing to +// make GrainSpeciesInfoDeleter a C++ class with a proper constructor & +// destructor +struct GrainSpeciesInfoDeleter { + void operator()(grackle::impl::GrainSpeciesInfo* ptr) const { + grackle::impl::drop_GrainSpeciesInfo(ptr); + delete ptr; + } +}; + +using unique_GrainSpeciesInfo_ptr = + std::unique_ptr; + +/// creates a unique_ptr holding a grackle::impl::GrainSpeciesInfo +/// +/// This is useful for preventing memory leaks when tests fail +unique_GrainSpeciesInfo_ptr make_unique_GrainSpeciesInfo( + int dust_species_param) { + grackle::impl::GrainSpeciesInfo* ptr = new grackle::impl::GrainSpeciesInfo; + (*ptr) = grackle::impl::new_GrainSpeciesInfo(dust_species_param); + return unique_GrainSpeciesInfo_ptr(ptr, GrainSpeciesInfoDeleter()); +} + +} // anonymous namespace + +// check that OnlyGrainSpLUT contains the entries for every known grain species +TEST(OnlyGrainSpLUTTest, CheckNumEntries) { + ASSERT_EQ(OnlyGrainSpLUT::NUM_ENTRIES, number_known_grain_species()); +} + +/// Define a fixture for running parameterized tests of the GrainSpeciesInfo +/// machinery. The tests are parameterized by the dust_species parameter. +class GrainSpeciesInfoTest : public testing::TestWithParam { +protected: + /// set up the grain_species_info pointer + /// + /// @note + /// We **ONLY** perform setup in this method (rather than in a default + /// constructor) because we want to perform some basic sanity checks + void SetUp() override { + int dust_species = GetParam(); + grain_species_info_ = make_unique_GrainSpeciesInfo(dust_species); + + // perform 3 sanity checks! + ASSERT_GE(grain_species_info_->n_species, 1) + << "GrainSpeciesInfo::n_species should be positive"; + ASSERT_LE(grain_species_info_->n_species, number_known_grain_species()) + << "GrainSpeciesInfo::n_species can't exceed the number of known grain " + << "species"; + ASSERT_NE(grain_species_info_->species_info, nullptr) + << "GrainSpeciesInfo::species_info can't be a nullptr when " + << "GrainSpeciesInfo::n_species is positive"; + } + + unique_GrainSpeciesInfo_ptr grain_species_info_; +}; + +TEST_P(GrainSpeciesInfoTest, CheckOnlyGrainSpeciesLUTConsistency) { + const int n_species = grain_species_info_->n_species; + for (int i = 0; i < n_species; i++) { + // sanity check! + ASSERT_NE(grain_species_info_->species_info[i].name, nullptr); + // actual check! + EXPECT_EQ(i, grain_species_info_->species_info[i].onlygrainsp_idx) + << "element " << i << " of the GrainSpeciesInfo::species_info array " + << "doesn't seem to be synchronized with the OnlyGrainSpeciesLUT " + << "enumeration. At face value (there aren't other related bugs), " + << "OnlyGrainSpeciesLUT::" << grain_species_info_->species_info[i].name + << " seems to have a value of " + << grain_species_info_->species_info[i].onlygrainsp_idx; + } +} + +TEST_P(GrainSpeciesInfoTest, SublimationTemperature) { + const int n_species = grain_species_info_->n_species; + for (int i = 0; i < n_species; i++) { + // sanity check! + ASSERT_NE(grain_species_info_->species_info[i].name, nullptr); + // actual check! + EXPECT_GT(grain_species_info_->species_info[i].sublimation_temperature, 0) + << "element " << i << " of the GrainSpeciesInfo::species_info array " + << "holds an invalid sublimation temperature"; + } +} + +TEST_P(GrainSpeciesInfoTest, SpeciesLUTCompare) { + // this is the reference list of all grain species constructed from the same + // XMacro as that the SpeciesLUT enumeration is built from + std::vector ref_list = grain_sp_info_from_species_list(); + + const int n_species = grain_species_info_->n_species; + for (int i = 0; i < n_species; i++) { + // sanity check! + ASSERT_NE(grain_species_info_->species_info[i].name, nullptr); + // actual check! + std::string actual_name(grain_species_info_->species_info[i].name); + EXPECT_EQ(actual_name, ref_list[i].name) + << "according to the reference list, the grain species at index " << i + << " should be `" << ref_list[i].name << "` (not `" << actual_name + << "`)"; + EXPECT_EQ(grain_species_info_->species_info[i].species_idx, + ref_list[i].index) + << "the value @ index " << i << " of GrainSpeciesInfo::species_info " + << "seems to indicate that the SpLUT::" << actual_name << " enumerator " + << "has a different value compared to the reference list"; + } +} + +INSTANTIATE_TEST_SUITE_P( + , // <- this comma is meaningful + GrainSpeciesInfoTest, testing::Range(1, MAX_dust_species_VAL + 1), + // adjust how the value is formatted in the test name + [](const testing::TestParamInfo& info) { + int val = info.param; + std::string name = "DustSpeciesEq" + std::to_string(val); + return name; + }); + +// check the GrainSpeciesInfo object when constructed from dust_species +// parameters that hold extreme values +TEST(GrainSpeciesInfoTestMisc, DustSpeciesExtremeValues) { + { + unique_GrainSpeciesInfo_ptr ptr = make_unique_GrainSpeciesInfo(0); + EXPECT_EQ(ptr->n_species, 0) + << "GrainSpeciesInfo::n_species should be 0 when the dust_species " + << "parameter is 0."; + EXPECT_EQ(ptr->species_info, nullptr) + << "GrainSpeciesInfo::species_info should be a nullptr when " + << "dust_species parameter is 0."; + } + + int invalid_dust_species_values[2] = {-42423, MAX_dust_species_VAL + 1}; + for (auto dust_species_param : invalid_dust_species_values) { + unique_GrainSpeciesInfo_ptr ptr = + make_unique_GrainSpeciesInfo(dust_species_param); + EXPECT_LE(ptr->n_species, -1) + << "GrainSpeciesInfo::n_species should be negative when the " + << "dust_species parameter is " << dust_species_param << " (i.e. an " + << "invalid value)."; + EXPECT_EQ(ptr->species_info, nullptr) + << "GrainSpeciesInfo::species_info member should be a nullptr when the " + << "dust_species parameter is " << dust_species_param << " (i.e. an " + << "invalid value)."; + } + + int n_known_grain_species = number_known_grain_species(); + { + unique_GrainSpeciesInfo_ptr ptr = + make_unique_GrainSpeciesInfo(MAX_dust_species_VAL); + EXPECT_EQ(ptr->n_species, n_known_grain_species) + << "GrainSpeciesInfo::n_species should specify the number of grain " + << "species known to grackle, " << n_known_grain_species + << ", when the " + << "dust_species parameter holds the max known value, " + << MAX_dust_species_VAL << '.'; + EXPECT_NE(ptr->species_info, nullptr) + << "GrainSpeciesInfo::species_info can't be a nullptr when " + << "GrainSpeciesInfo::n_species is positive"; + } +} + +// check that the grackle::impl::max_ingredients_per_grain_species constant +// indeed holds the maximum number of ingredients that a grain species can +// have. +TEST(GrainSpeciesInfoTestMisc, MaxIngredientsPerGrainSpecies) { + // construct a GrainSpeciesInfo instance that holds values for every dust + // grain species + unique_GrainSpeciesInfo_ptr gsp_info = + make_unique_GrainSpeciesInfo(MAX_dust_species_VAL); + + // go through and compute the number of ingredients for each species + int n_grain_species = gsp_info->n_species; + int max_ingredient_count = 0; + for (int gsp_idx = 0; gsp_idx < n_grain_species; gsp_idx++) { + max_ingredient_count = + std::max(max_ingredient_count, + gsp_info->species_info[gsp_idx].n_growth_ingredients); + } + + EXPECT_EQ(max_ingredient_count, + grackle::impl::max_ingredients_per_grain_species) + << "it appears that grackle::impl::max_ingredients_per_grain_species " + << "should have a value of " << max_ingredient_count; +} From 7a780047d51231e4414c7d526990705146b8c5bc Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Fri, 26 Sep 2025 11:32:15 -0400 Subject: [PATCH 02/11] add some commentary --- .clang-format-ignore | 1 - src/clib/dust_props.hpp | 61 +++++++++++++++++++------------------ src/clib/internal_types.hpp | 8 +++++ 3 files changed, 40 insertions(+), 30 deletions(-) diff --git a/.clang-format-ignore b/.clang-format-ignore index 36fa238b6..ab4b6b6f1 100644 --- a/.clang-format-ignore +++ b/.clang-format-ignore @@ -21,7 +21,6 @@ src/clib/cie_thin_cooling_rate_tables.h src/clib/cool1d_multi_g-cpp.C src/clib/cool1d_multi_g-cpp.h src/clib/cool_multi_time_g.cpp -src/clib/dust_props.hpp src/clib/dynamic_api.c src/clib/fortran_func_decls.h src/clib/fortran_func_wrappers.hpp diff --git a/src/clib/dust_props.hpp b/src/clib/dust_props.hpp index c1daa9c4b..8639b3833 100644 --- a/src/clib/dust_props.hpp +++ b/src/clib/dust_props.hpp @@ -1,7 +1,12 @@ -// See LICENSE file for license and copyright information - -/// @file dust_props.hpp -/// @brief Declarations related to computing dust properties. +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Declarations related to computing dust properties. /// /// In slightly more detail, these is a series of dust-related calculations /// that are largely independent of the gas chemistry. These calculations all @@ -9,22 +14,27 @@ /// filled by calc_grain_size_increment_1d. They are used to compute quantities /// that is important for the gas /// -> by calc_all_tdust_gasgr_1d_g -/// -> by some logic lookup_cool_rates1d_g (this logic should be excised) +/// -> by some logic lookup_cool_rates1d (this logic should be excised) /// /// The premise is that we should try to keep these properties as separated as /// possible from the gas-physics +/// +//===----------------------------------------------------------------------===// #ifndef DUST_PROPS_HPP #define DUST_PROPS_HPP -#include // malloc, free +// In the near future (after GH PR #405 is merged), we plan to: +// 1. rename InternalDustPropBuf to InternalGrainPropBuf so that type's role in +// the multi-grain-species dust model is more obvious +// 2. move this file into the dust subdirectory + #include "internal_types.hpp" #include "LUT.hpp" namespace grackle::impl { struct InternalDustPropBuf { - // a decent argument could be made that this should track grain temperature /// geometric cross-section per unit gas mass of each grain species @@ -39,7 +49,7 @@ struct InternalDustPropBuf { grackle::impl::GrainSpeciesCollection grain_sigma_per_gas_mass; /// buffer closely related to grain_sigma_per_gas_mass - double * sigma_per_gas_mass_tot; + double* sigma_per_gas_mass_tot; /// holds 1D opacity tables at each zone for each grain species zone /// - the values are tabulated with respect to dust temperature. After @@ -85,8 +95,7 @@ struct InternalDustPropBuf { grackle::impl::GrainSpeciesCollection grain_dyntab_kappa; /// buffer closely related to grain_dyntab_kappa - double * dyntab_kappa_tot; - + double* dyntab_kappa_tot; }; /// allocates the contents of a new InternalDustPropBuf @@ -94,21 +103,18 @@ struct InternalDustPropBuf { /// @param nelem The number of elements in each buffer /// @param opacity_table_size The number of elements in a dynamically computed /// opacity table -inline InternalDustPropBuf new_InternalDustPropBuf( - int nelem, int opacity_table_size -) { +inline InternalDustPropBuf new_InternalDustPropBuf(int nelem, + int opacity_table_size) { InternalDustPropBuf out; - out.grain_sigma_per_gas_mass = grackle::impl::new_GrainSpeciesCollection( - nelem - ); - out.sigma_per_gas_mass_tot = (double*)std::malloc(sizeof(double)*nelem); + out.grain_sigma_per_gas_mass = + grackle::impl::new_GrainSpeciesCollection(nelem); + out.sigma_per_gas_mass_tot = new double[nelem]; // this is a little bit of a hack! int table_len = (opacity_table_size < 1) ? 1 : opacity_table_size; - out.grain_dyntab_kappa = grackle::impl::new_GrainSpeciesCollection( - nelem * table_len - ); - out.dyntab_kappa_tot = (double*)std::malloc(sizeof(double)*nelem*table_len); + out.grain_dyntab_kappa = + grackle::impl::new_GrainSpeciesCollection(nelem * table_len); + out.dyntab_kappa_tot = new double[nelem * table_len]; return out; } @@ -118,15 +124,12 @@ inline InternalDustPropBuf new_InternalDustPropBuf( /// This effectively invokes the destructor inline void drop_InternalDustPropBuf(InternalDustPropBuf* ptr) { grackle::impl::drop_GrainSpeciesCollection(&ptr->grain_sigma_per_gas_mass); - std::free(ptr->sigma_per_gas_mass_tot); + delete[] ptr->sigma_per_gas_mass_tot; - grackle::impl::drop_GrainSpeciesCollection( - &ptr->grain_dyntab_kappa - ); - std::free(ptr->dyntab_kappa_tot); + grackle::impl::drop_GrainSpeciesCollection(&ptr->grain_dyntab_kappa); + delete[] ptr->dyntab_kappa_tot; } +} // namespace grackle::impl -} // namespace grackle::impl - -#endif // DUST_PROPS_HPP +#endif // DUST_PROPS_HPP diff --git a/src/clib/internal_types.hpp b/src/clib/internal_types.hpp index 76d70d859..7aaa30fa7 100644 --- a/src/clib/internal_types.hpp +++ b/src/clib/internal_types.hpp @@ -440,6 +440,14 @@ void drop_SpeciesCollection(SpeciesCollection*); /// @note /// This is something we may want to reuse. If we are willing to embrace C++, /// then we may want to use templates +/// +/// @important +/// In the near-term to mid-term, we want to refactor to move away from using +/// @ref OnlyGrainSpLUT and instead start relying upon the order of grain +/// species tracked by @ref grackle::impl::GrainSpeciesInfo (at the moment +/// grains have the same order). This is possible because every block of logic +/// using this struct applies the same logic to each grain species (in other +/// words, the logic can be replaced with a for-loop over the grain species) struct GrainSpeciesCollection { double* data[OnlyGrainSpLUT::NUM_ENTRIES]; }; From cd486d866970a5b5e604ced866b8a74cadad7432 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Thu, 9 Oct 2025 12:33:46 -0400 Subject: [PATCH 03/11] address a memory leak and perform some light refactoring --- src/clib/initialize_chemistry_data.cpp | 26 ++++++++++++++++++++++++-- src/clib/initialize_rates.cpp | 2 +- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/src/clib/initialize_chemistry_data.cpp b/src/clib/initialize_chemistry_data.cpp index ef75340db..b8b39fabf 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -630,13 +630,35 @@ extern "C" int local_free_chemistry_data(chemistry_data *my_chemistry, return FAIL; } + // start freeing memory associated with opaque storage + // --------------------------------------------------- if (my_rates->opaque_storage->kcol_rate_tables != nullptr) { + // delete contents of kcol_rate_tables drop_CollisionalRxnRateCollection(my_rates->opaque_storage->kcol_rate_tables); + // delete kcol_rate_tables, itself delete my_rates->opaque_storage->kcol_rate_tables; } - delete[] my_rates->opaque_storage->used_kcol_rate_indices; + + if (my_rates->opaque_storage->used_kcol_rate_indices !=nullptr) { + // since used_kcol_rate_indices are just integers, we can directly + // deallocate them + delete[] my_rates->opaque_storage->used_kcol_rate_indices; + } + + // delete contents of h2dust_grain_interp_props (automatically handles the + // case where we didn't allocate anything) free_interp_grid_props_(&my_rates->opaque_storage->h2dust_grain_interp_props); - GRACKLE_FREE(my_rates->opaque_storage->grain_species_info); + // since h2dust_grain_interp_props isn't a pointer, there is nothing more to + // allocate right here + + if (my_rates->opaque_storage->grain_species_info != nullptr) { + // delete contents of grain_species_info + grackle::impl::drop_GrainSpeciesInfo( + my_rates->opaque_storage->grain_species_info); + // delete kcol_rate_tables, itself + delete my_rates->opaque_storage->grain_species_info; + } + delete my_rates->opaque_storage; my_rates->opaque_storage = nullptr; diff --git a/src/clib/initialize_rates.cpp b/src/clib/initialize_rates.cpp index 9451ef00f..652aff070 100644 --- a/src/clib/initialize_rates.cpp +++ b/src/clib/initialize_rates.cpp @@ -650,7 +650,7 @@ int grackle::impl::initialize_rates( // (it may make sense want to handle more of the dust separately) if (my_chemistry->dust_species > 0) { my_rates->opaque_storage->grain_species_info = - (grackle::impl::GrainSpeciesInfo*)malloc(sizeof(grackle::impl::GrainSpeciesInfo)); + new grackle::impl::GrainSpeciesInfo; *(my_rates->opaque_storage->grain_species_info) = grackle::impl::new_GrainSpeciesInfo(my_chemistry->dust_species); } From 8122ee4a2e37a511d5a68d4fbc73d96f4bbf910a Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sat, 6 Dec 2025 16:58:22 -0500 Subject: [PATCH 04/11] light reformatting --- src/clib/dust/grain_species_info.hpp | 40 ++++++++++++++-------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/src/clib/dust/grain_species_info.hpp b/src/clib/dust/grain_species_info.hpp index d6eb14ff7..dc55037a2 100644 --- a/src/clib/dust/grain_species_info.hpp +++ b/src/clib/dust/grain_species_info.hpp @@ -226,9 +226,9 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { // {coef, species_idx, particle mass} {1, SpLUT::CI, 12.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[1] = - mk_gsp_info_entry_helper_(SpLUT::AC_dust, OnlyGrainSpLUT::AC_dust, - "AC_dust", true, 1800.0, AC_dust_ingred); + out.species_info[1] = mk_gsp_info_entry_helper_( + SpLUT::AC_dust, OnlyGrainSpLUT::AC_dust, "AC_dust", true, + 1800.0, AC_dust_ingred); } if (dust_species_parameter > 1) { @@ -236,17 +236,17 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { // {coef, species_idx, particle mass} {1, SpLUT::SiI, 28.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[2] = - mk_gsp_info_entry_helper_(SpLUT::SiM_dust, OnlyGrainSpLUT::SiM_dust, - "SiM_dust", false, 1500.0, SiM_dust_ingred); + out.species_info[2] = mk_gsp_info_entry_helper_( + SpLUT::SiM_dust, OnlyGrainSpLUT::SiM_dust, "SiM_dust", false, + 1500.0, SiM_dust_ingred); const GrainGrowthIngredient FeM_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::Fe, 56.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[3] = - mk_gsp_info_entry_helper_(SpLUT::FeM_dust, OnlyGrainSpLUT::FeM_dust, - "FeM_dust", false, 1500.0, FeM_dust_ingred); + out.species_info[3] = mk_gsp_info_entry_helper_( + SpLUT::FeM_dust, OnlyGrainSpLUT::FeM_dust, "FeM_dust", false, + 1500.0, FeM_dust_ingred); const GrainGrowthIngredient Mg2SiO4_dust_ingred[] = { // {coef, species_idx, particle mass} @@ -255,8 +255,8 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { {3, SpLUT::H2O, 18.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; out.species_info[4] = mk_gsp_info_entry_helper_( - SpLUT::Mg2SiO4_dust, OnlyGrainSpLUT::Mg2SiO4_dust, "Mg2SiO4_dust", - false, 1277.0, Mg2SiO4_dust_ingred); + SpLUT::Mg2SiO4_dust, OnlyGrainSpLUT::Mg2SiO4_dust, "Mg2SiO4_dust", false, + 1277.0, Mg2SiO4_dust_ingred); const GrainGrowthIngredient Fe3O4_dust_ingred[] = { // {coef, species_idx, particle mass} @@ -271,27 +271,27 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { // {coef, species_idx, particle mass} {1, SpLUT::SiO2I, 60.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[6] = - mk_gsp_info_entry_helper_(SpLUT::SiO2_dust, OnlyGrainSpLUT::SiO2_dust, - "SiO2_dust", false, 1500.0, SiO2_dust_ingred); + out.species_info[6] = mk_gsp_info_entry_helper_( + SpLUT::SiO2_dust, OnlyGrainSpLUT::SiO2_dust, "SiO2_dust", false, + 1500.0, SiO2_dust_ingred); const GrainGrowthIngredient MgO_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::Mg, 24.}, {1, SpLUT::H2O, 18.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[7] = - mk_gsp_info_entry_helper_(SpLUT::MgO_dust, OnlyGrainSpLUT::MgO_dust, - "MgO_dust", false, 1500.0, MgO_dust_ingred); + out.species_info[7] = mk_gsp_info_entry_helper_( + SpLUT::MgO_dust, OnlyGrainSpLUT::MgO_dust, "MgO_dust", false, + 1500.0, MgO_dust_ingred); const GrainGrowthIngredient FeS_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::Fe, 56.}, {1, SpLUT::S, 32.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[8] = - mk_gsp_info_entry_helper_(SpLUT::FeS_dust, OnlyGrainSpLUT::FeS_dust, - "FeS_dust", false, 680.0, FeS_dust_ingred); + out.species_info[8] = mk_gsp_info_entry_helper_( + SpLUT::FeS_dust, OnlyGrainSpLUT::FeS_dust, "FeS_dust", false, + 680.0, FeS_dust_ingred); const GrainGrowthIngredient Al2O3_dust_ingred[] = { // {coef, species_idx, particle mass} From 19fb8f5ea1fc1a52047d9228001def5c723f96a7 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sun, 7 Dec 2025 12:18:07 -0500 Subject: [PATCH 05/11] make the formatting of GrainSpeciesInfo initializers more legible --- src/clib/dust/grain_species_info.hpp | 104 ++++++++++++++++++++------- 1 file changed, 78 insertions(+), 26 deletions(-) diff --git a/src/clib/dust/grain_species_info.hpp b/src/clib/dust/grain_species_info.hpp index dc55037a2..c1f6a8073 100644 --- a/src/clib/dust/grain_species_info.hpp +++ b/src/clib/dust/grain_species_info.hpp @@ -219,16 +219,24 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { {2, SpLUT::H2O, 18.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; out.species_info[0] = mk_gsp_info_entry_helper_( - SpLUT::MgSiO3_dust, OnlyGrainSpLUT::MgSiO3_dust, "MgSiO3_dust", false, - 1222.0, MgSiO3_dust_ingred); + /* species_idx = */ SpLUT::MgSiO3_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::MgSiO3_dust, + /* name = */ "MgSiO3_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1222.0, + /* growth_ingredients = */ MgSiO3_dust_ingred); const GrainGrowthIngredient AC_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::CI, 12.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; out.species_info[1] = mk_gsp_info_entry_helper_( - SpLUT::AC_dust, OnlyGrainSpLUT::AC_dust, "AC_dust", true, - 1800.0, AC_dust_ingred); + /* species_idx = */ SpLUT::AC_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::AC_dust, + /* name = */ "AC_dust", + /* h2dust_uses_carbonaceous_table = */ true, + /* sublimation_temperature = */ 1800.0, + /* growth_ingredients = */ AC_dust_ingred); } if (dust_species_parameter > 1) { @@ -237,16 +245,24 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { {1, SpLUT::SiI, 28.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; out.species_info[2] = mk_gsp_info_entry_helper_( - SpLUT::SiM_dust, OnlyGrainSpLUT::SiM_dust, "SiM_dust", false, - 1500.0, SiM_dust_ingred); + /* species_idx = */ SpLUT::SiM_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::SiM_dust, + /* name = */ "SiM_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* growth_ingredients = */ SiM_dust_ingred); const GrainGrowthIngredient FeM_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::Fe, 56.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; out.species_info[3] = mk_gsp_info_entry_helper_( - SpLUT::FeM_dust, OnlyGrainSpLUT::FeM_dust, "FeM_dust", false, - 1500.0, FeM_dust_ingred); + /* species_idx = */ SpLUT::FeM_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::FeM_dust, + /* name = */ "FeM_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* growth_ingredients = */ FeM_dust_ingred); const GrainGrowthIngredient Mg2SiO4_dust_ingred[] = { // {coef, species_idx, particle mass} @@ -255,8 +271,12 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { {3, SpLUT::H2O, 18.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; out.species_info[4] = mk_gsp_info_entry_helper_( - SpLUT::Mg2SiO4_dust, OnlyGrainSpLUT::Mg2SiO4_dust, "Mg2SiO4_dust", false, - 1277.0, Mg2SiO4_dust_ingred); + /* species_idx = */ SpLUT::Mg2SiO4_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::Mg2SiO4_dust, + /* name = */ "Mg2SiO4_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1277.0, + /* growth_ingredients = */ Mg2SiO4_dust_ingred); const GrainGrowthIngredient Fe3O4_dust_ingred[] = { // {coef, species_idx, particle mass} @@ -264,16 +284,24 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { {4, SpLUT::H2O, 18.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; out.species_info[5] = mk_gsp_info_entry_helper_( - SpLUT::Fe3O4_dust, OnlyGrainSpLUT::Fe3O4_dust, "Fe3O4_dust", false, - 1500.0, Fe3O4_dust_ingred); + /* species_idx = */ SpLUT::Fe3O4_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::Fe3O4_dust, + /* name = */ "Fe3O4_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* growth_ingredients = */ Fe3O4_dust_ingred); const GrainGrowthIngredient SiO2_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::SiO2I, 60.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; out.species_info[6] = mk_gsp_info_entry_helper_( - SpLUT::SiO2_dust, OnlyGrainSpLUT::SiO2_dust, "SiO2_dust", false, - 1500.0, SiO2_dust_ingred); + /* species_idx = */ SpLUT::SiO2_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::SiO2_dust, + /* name = */ "SiO2_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* growth_ingredients = */ SiO2_dust_ingred); const GrainGrowthIngredient MgO_dust_ingred[] = { // {coef, species_idx, particle mass} @@ -281,8 +309,12 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { {1, SpLUT::H2O, 18.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; out.species_info[7] = mk_gsp_info_entry_helper_( - SpLUT::MgO_dust, OnlyGrainSpLUT::MgO_dust, "MgO_dust", false, - 1500.0, MgO_dust_ingred); + /* species_idx = */ SpLUT::MgO_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::MgO_dust, + /* name = */ "MgO_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* growth_ingredients = */ MgO_dust_ingred); const GrainGrowthIngredient FeS_dust_ingred[] = { // {coef, species_idx, particle mass} @@ -290,8 +322,12 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { {1, SpLUT::S, 32.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; out.species_info[8] = mk_gsp_info_entry_helper_( - SpLUT::FeS_dust, OnlyGrainSpLUT::FeS_dust, "FeS_dust", false, - 680.0, FeS_dust_ingred); + /* species_idx = */ SpLUT::FeS_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::FeS_dust, + /* name = */ "FeS_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 680.0, + /* growth_ingredients = */ FeS_dust_ingred); const GrainGrowthIngredient Al2O3_dust_ingred[] = { // {coef, species_idx, particle mass} @@ -299,8 +335,12 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { {3, SpLUT::H2O, 18.}, GRIMPL_INGREDIENT_LIST_SENTINEL}; out.species_info[9] = mk_gsp_info_entry_helper_( - SpLUT::Al2O3_dust, OnlyGrainSpLUT::Al2O3_dust, "Al2O3_dust", false, - 1500.0, Al2O3_dust_ingred); + /* species_idx = */ SpLUT::Al2O3_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::Al2O3_dust, + /* name = */ "Al2O3_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* growth_ingredients = */ Al2O3_dust_ingred); } if (dust_species_parameter > 2) { @@ -309,14 +349,26 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { // are low (100-600 K). They sublimate before the growth occurs. out.species_info[10] = mk_gsp_info_entry_helper_( - SpLUT::ref_org_dust, OnlyGrainSpLUT::ref_org_dust, "ref_org_dust", - false, 575.0, nullptr); + /* species_idx = */ SpLUT::ref_org_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::ref_org_dust, + /* name = */ "ref_org_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 575.0, + /* growth_ingredients = */ nullptr); out.species_info[11] = mk_gsp_info_entry_helper_( - SpLUT::vol_org_dust, OnlyGrainSpLUT::vol_org_dust, "vol_org_dust", - false, 375.0, nullptr); + /* species_idx = */ SpLUT::vol_org_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::vol_org_dust, + /* name = */ "vol_org_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 375.0, + /* growth_ingredients = */ nullptr); out.species_info[12] = mk_gsp_info_entry_helper_( - SpLUT::H2O_ice_dust, OnlyGrainSpLUT::H2O_ice_dust, "H2O_ice_dust", - false, 153.0, nullptr); + /* species_idx = */ SpLUT::H2O_ice_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::H2O_ice_dust, + /* name = */ "H2O_ice_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 153.0, + /* growth_ingredients = */ nullptr); } return out; From 50273dc839653a77e734c2ccf75e5482bbdccbfc Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sun, 7 Dec 2025 13:26:25 -0500 Subject: [PATCH 06/11] Start tracking bulk density of each grain species within GrainSpeciesInfoEntry This information will only be used after we transcribe `calc_grain_size_increment_1d` (at which point we will remove the corresponding constants from phys_constants.h and can entirely delete the dust_const.def file --- src/clib/dust/grain_species_info.hpp | 26 ++++++++++++++++++++++++-- src/clib/phys_constants.h | 7 ++++++- 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/src/clib/dust/grain_species_info.hpp b/src/clib/dust/grain_species_info.hpp index c1f6a8073..1ae31191a 100644 --- a/src/clib/dust/grain_species_info.hpp +++ b/src/clib/dust/grain_species_info.hpp @@ -71,9 +71,17 @@ struct GrainSpeciesInfoEntry { /// to track the information in this manner. bool h2dust_uses_carbonaceous_table; - /// The sublimation temperature + /// The sublimation temperature in units of Kelvin double sublimation_temperature; + /// specifies the density of a single grain in units of g/cm^3 + /// + /// @note + /// The values are consistent with the values quoted within table 2 of + /// [Chiaki+ 2015](https://ui.adsabs.harvard.edu/abs/2015MNRAS.446.2659C) + /// assuming that all grains are spherical. + double bulk_density_cgs; + /// The number of growth ingredients int n_growth_ingredients; @@ -149,7 +157,7 @@ inline constexpr int max_ingredients_per_grain_species = 3; inline GrainSpeciesInfoEntry mk_gsp_info_entry_helper_( int species_idx, int onlygrainsp_idx, const char* name, bool h2dust_uses_carbonaceous_table, double sublimation_temperature, - const GrainGrowthIngredient* growth_ingredients) { + double bulk_density_cgs, const GrainGrowthIngredient* growth_ingredients) { GrainGrowthIngredient* out_ingredient_ptr = nullptr; int n_ingredients = 0; @@ -170,6 +178,7 @@ inline GrainSpeciesInfoEntry mk_gsp_info_entry_helper_( name, h2dust_uses_carbonaceous_table, sublimation_temperature, + bulk_density_cgs, n_ingredients, out_ingredient_ptr}; } @@ -224,6 +233,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "MgSiO3_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 1222.0, + /* bulk_density_cgs = */ 3.20185, /* growth_ingredients = */ MgSiO3_dust_ingred); const GrainGrowthIngredient AC_dust_ingred[] = { @@ -236,6 +246,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "AC_dust", /* h2dust_uses_carbonaceous_table = */ true, /* sublimation_temperature = */ 1800.0, + /* bulk_density_cgs = */ 2.27949, /* growth_ingredients = */ AC_dust_ingred); } @@ -250,6 +261,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "SiM_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 2.34118, /* growth_ingredients = */ SiM_dust_ingred); const GrainGrowthIngredient FeM_dust_ingred[] = { @@ -262,6 +274,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "FeM_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 7.95995, /* growth_ingredients = */ FeM_dust_ingred); const GrainGrowthIngredient Mg2SiO4_dust_ingred[] = { @@ -276,6 +289,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "Mg2SiO4_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 1277.0, + /* bulk_density_cgs = */ 3.22133, /* growth_ingredients = */ Mg2SiO4_dust_ingred); const GrainGrowthIngredient Fe3O4_dust_ingred[] = { @@ -289,6 +303,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "Fe3O4_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 5.25096, /* growth_ingredients = */ Fe3O4_dust_ingred); const GrainGrowthIngredient SiO2_dust_ingred[] = { @@ -301,6 +316,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "SiO2_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 2.66235, /* growth_ingredients = */ SiO2_dust_ingred); const GrainGrowthIngredient MgO_dust_ingred[] = { @@ -314,6 +330,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "MgO_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 3.58157, /* growth_ingredients = */ MgO_dust_ingred); const GrainGrowthIngredient FeS_dust_ingred[] = { @@ -327,6 +344,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "FeS_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 680.0, + /* bulk_density_cgs = */ 4.87265, /* growth_ingredients = */ FeS_dust_ingred); const GrainGrowthIngredient Al2O3_dust_ingred[] = { @@ -340,6 +358,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "Al2O3_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 4.01610, /* growth_ingredients = */ Al2O3_dust_ingred); } @@ -354,6 +373,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "ref_org_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 575.0, + /* bulk_density_cgs = */ 1.5, /* growth_ingredients = */ nullptr); out.species_info[11] = mk_gsp_info_entry_helper_( /* species_idx = */ SpLUT::vol_org_dust, @@ -361,6 +381,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "vol_org_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 375.0, + /* bulk_density_cgs = */ 1.0, /* growth_ingredients = */ nullptr); out.species_info[12] = mk_gsp_info_entry_helper_( /* species_idx = */ SpLUT::H2O_ice_dust, @@ -368,6 +389,7 @@ inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { /* name = */ "H2O_ice_dust", /* h2dust_uses_carbonaceous_table = */ false, /* sublimation_temperature = */ 153.0, + /* bulk_density_cgs = */ 0.92, /* growth_ingredients = */ nullptr); } diff --git a/src/clib/phys_constants.h b/src/clib/phys_constants.h index bcd916aa8..a2362296a 100644 --- a/src/clib/phys_constants.h +++ b/src/clib/phys_constants.h @@ -115,7 +115,12 @@ /************************************************/ -/* (we may want to relocate these to a different file in the future) */ +/* TODO: After we finish transcribe calc_grain_size_increment_1d to C++, we + * should remove these constants and instead use the values stored by + * `GrainSpeciesInfoEntry::bulk_density_cgs` + * -> to be clear, the values within GrainSpeciesInfoEntry::bulk_density_cgs + * are hardcoded based on the values currently held by these constants. + */ #define sSiM 2.34118e0 #define sFeM 7.95995e0 From b7b13fcd018361e2d3d034b303de67ff998ff2cc Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sun, 7 Dec 2025 14:20:31 -0500 Subject: [PATCH 07/11] split out part of grain_species_info.hpp into a source file to improve readability (and add some comments) --- src/clib/CMakeLists.txt | 2 +- src/clib/Make.config.objects | 1 + src/clib/dust/grain_species_info.cpp | 297 +++++++++++++++++++++++++++ src/clib/dust/grain_species_info.hpp | 254 +---------------------- 4 files changed, 302 insertions(+), 252 deletions(-) create mode 100644 src/clib/dust/grain_species_info.cpp diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index 8c43e456e..c9ae9f141 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -112,7 +112,7 @@ add_library(Grackle_Grackle cool1d_multi_g.cpp cool1d_multi_g.hpp cool_multi_time_g.cpp cool_multi_time_g.h dust_props.hpp - dust/grain_species_info.hpp + dust/grain_species_info.cpp dust/grain_species_info.hpp init_misc_species_cool_rates.cpp init_misc_species_cool_rates.hpp initialize_chemistry_data.cpp initialize_dust_yields.cpp initialize_dust_yields.hpp diff --git a/src/clib/Make.config.objects b/src/clib/Make.config.objects index 2fa8a89ee..7a4bb87af 100644 --- a/src/clib/Make.config.objects +++ b/src/clib/Make.config.objects @@ -25,6 +25,7 @@ OBJS_CONFIG_LIB = \ cool1d_cloudy_old_tables_g.lo \ cool1d_multi_g.lo \ cool_multi_time_g.lo \ + dust/grain_species_info.lo \ dynamic_api.lo \ grackle_units.lo \ index_helper.lo \ diff --git a/src/clib/dust/grain_species_info.cpp b/src/clib/dust/grain_species_info.cpp new file mode 100644 index 000000000..e06a4d84c --- /dev/null +++ b/src/clib/dust/grain_species_info.cpp @@ -0,0 +1,297 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Implements functions associated with GrainSpeciesInfo (namely, the +/// new_GrainSpeciesInfo function) +/// +//===----------------------------------------------------------------------===// + +#include // memcpy + +#include "LUT.hpp" +#include "grain_species_info.hpp" + +// The following logic effectively does 2 (related things): +// 1. it serves as a human-readable registry of all known grain species and +// their associated properties. +// 2. we must be able to use this information to construct an instance of the +// GrainSpeciesInfo data structure (so that the properties of each species +// can be queried programatically) + +// Encoding the gas-phase ingredients (for grain growth) of each grain species +// poses a bit of an obstacle. Essentially, we need a way to specify a +// variable-length list of ingredients and we should be able to infer the +// length. Thus, similar to the way a c-string has a sentinel value signalling +// the end of the list (i.e. the '\0' character), we use a sentinel for the +// hardcoded ingredient lists. This is intended as a stopgap solution! +// +// Note: the version of the lists accessible through GrainSpeciesInfo NEVER +// have the sentinel! + +/// Expands to an initializer list for initializing an instance of the +/// GrainGrowthIngredient struct that is used to signal the termination of +/// a variable length ingredient list (analogous to the role that '\0' plays +/// in a c-string). +#define GRIMPL_INGREDIENT_LIST_SENTINEL {-1, -1, -1.0} + +namespace { // stuff inside an anonymous namespace is local to this file + +/// Constructs a GrainSpeciesInfoEntry instance +/// +/// The primary work that this function must complete is determine the number +/// of ingredients that a grain species has (if any) and then construct a +/// heap-allocated copy of the ingredient list that is stored by the returned +/// GrainSpeciesInfoEntry instance. Importantly: +/// - the @p growth_ingredients argument is terminated by the sentinel +/// - the ingredient list in the returned instance is **NOT** terminated by +/// the sentinel +grackle::impl::GrainSpeciesInfoEntry mk_gsp_info_entry_helper_( + int species_idx, int onlygrainsp_idx, const char* name, + bool h2dust_uses_carbonaceous_table, double sublimation_temperature, + double bulk_density_cgs, + const grackle::impl::GrainGrowthIngredient* growth_ingredients) { + using grackle::impl::GrainGrowthIngredient; + + // the main "work" is to determine the number of specified growth_ingredients + // and copy them (if there are no ingredients, we use a nullptr) + GrainGrowthIngredient* out_ingredient_ptr = nullptr; + int n_ingredients = 0; + + if (growth_ingredients != nullptr) { + n_ingredients = 0; + // increment the counter until we reach the sentinel + while (growth_ingredients[n_ingredients].coef != -1) { + n_ingredients++; + } + // allocate and initialize the pointer + out_ingredient_ptr = + new grackle::impl::GrainGrowthIngredient[n_ingredients]; + std::memcpy(out_ingredient_ptr, growth_ingredients, + (std::size_t)n_ingredients); + } + + return grackle::impl::GrainSpeciesInfoEntry{species_idx, + onlygrainsp_idx, + name, + h2dust_uses_carbonaceous_table, + sublimation_temperature, + bulk_density_cgs, + n_ingredients, + out_ingredient_ptr}; +} + +} // anonymous namespace + +grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( + int dust_species_parameter) { + GrainSpeciesInfo out{-1, nullptr}; // indicates an error + out.n_species = get_n_grain_species(dust_species_parameter); + + if (out.n_species <= 0) { + return out; + } + + out.species_info = new GrainSpeciesInfoEntry[out.n_species]; + + // At the time of writing: + // - we **only** use h2rate_carbonaceous_coef_table for the AC_dust + // species (amorphous carbon grains) + // - we use h2rate_silicate_coef_table for **ALL** other grain species + // (including species that are obviously NOT silicates) + // + // It is not obvious at all why this is the case... + // + // TODO: we should add documentation clearly addressing each of the + // following questions and replace this todo-item with a reference to the + // part of the documentation that discusses this! The documentation + // should make sure to address: + // - if this choice physically motivated or a common convention? (or if + // it was a quick & dirty choice because Gen didn't know what to do) + // - do we have a sense for how "wrong" the choice is? (e.g. will it + // generally overpredict/underpredict?) + // - why don't we use the generic my_rates->h2dusta rate as a generic + // fallback for non-silicate grains instead of using h2dustS? I realize + // that h2dusta has very different units... + + if (dust_species_parameter > 0) { + const GrainGrowthIngredient MgSiO3_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::Mg, 24.}, + {1, SpLUT::SiOI, 44.}, + {2, SpLUT::H2O, 18.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[0] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::MgSiO3_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::MgSiO3_dust, + /* name = */ "MgSiO3_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1222.0, + /* bulk_density_cgs = */ 3.20185, + /* growth_ingredients = */ MgSiO3_dust_ingred); + + const GrainGrowthIngredient AC_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::CI, 12.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[1] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::AC_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::AC_dust, + /* name = */ "AC_dust", + /* h2dust_uses_carbonaceous_table = */ true, + /* sublimation_temperature = */ 1800.0, + /* bulk_density_cgs = */ 2.27949, + /* growth_ingredients = */ AC_dust_ingred); + } + + if (dust_species_parameter > 1) { + const GrainGrowthIngredient SiM_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::SiI, 28.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[2] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::SiM_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::SiM_dust, + /* name = */ "SiM_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 2.34118, + /* growth_ingredients = */ SiM_dust_ingred); + + const GrainGrowthIngredient FeM_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::Fe, 56.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[3] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::FeM_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::FeM_dust, + /* name = */ "FeM_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 7.95995, + /* growth_ingredients = */ FeM_dust_ingred); + + const GrainGrowthIngredient Mg2SiO4_dust_ingred[] = { + // {coef, species_idx, particle mass} + {2, SpLUT::Mg, 24.}, + {1, SpLUT::SiOI, 44.}, + {3, SpLUT::H2O, 18.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[4] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::Mg2SiO4_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::Mg2SiO4_dust, + /* name = */ "Mg2SiO4_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1277.0, + /* bulk_density_cgs = */ 3.22133, + /* growth_ingredients = */ Mg2SiO4_dust_ingred); + + const GrainGrowthIngredient Fe3O4_dust_ingred[] = { + // {coef, species_idx, particle mass} + {3, SpLUT::Fe, 56.}, + {4, SpLUT::H2O, 18.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[5] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::Fe3O4_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::Fe3O4_dust, + /* name = */ "Fe3O4_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 5.25096, + /* growth_ingredients = */ Fe3O4_dust_ingred); + + const GrainGrowthIngredient SiO2_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::SiO2I, 60.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[6] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::SiO2_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::SiO2_dust, + /* name = */ "SiO2_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 2.66235, + /* growth_ingredients = */ SiO2_dust_ingred); + + const GrainGrowthIngredient MgO_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::Mg, 24.}, + {1, SpLUT::H2O, 18.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[7] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::MgO_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::MgO_dust, + /* name = */ "MgO_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 3.58157, + /* growth_ingredients = */ MgO_dust_ingred); + + const GrainGrowthIngredient FeS_dust_ingred[] = { + // {coef, species_idx, particle mass} + {1, SpLUT::Fe, 56.}, + {1, SpLUT::S, 32.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[8] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::FeS_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::FeS_dust, + /* name = */ "FeS_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 680.0, + /* bulk_density_cgs = */ 4.87265, + /* growth_ingredients = */ FeS_dust_ingred); + + const GrainGrowthIngredient Al2O3_dust_ingred[] = { + // {coef, species_idx, particle mass} + {2, SpLUT::Al, 27.}, + {3, SpLUT::H2O, 18.}, + GRIMPL_INGREDIENT_LIST_SENTINEL}; + out.species_info[9] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::Al2O3_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::Al2O3_dust, + /* name = */ "Al2O3_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 1500.0, + /* bulk_density_cgs = */ 4.01610, + /* growth_ingredients = */ Al2O3_dust_ingred); + } + + if (dust_species_parameter > 2) { + // We do not consider the growth of refractory organics, volatile + // organics, and water ice because their sublimation temperatures + // are low (100-600 K). They sublimate before the growth occurs. + + out.species_info[10] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::ref_org_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::ref_org_dust, + /* name = */ "ref_org_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 575.0, + /* bulk_density_cgs = */ 1.5, + /* growth_ingredients = */ nullptr); + out.species_info[11] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::vol_org_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::vol_org_dust, + /* name = */ "vol_org_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 375.0, + /* bulk_density_cgs = */ 1.0, + /* growth_ingredients = */ nullptr); + out.species_info[12] = mk_gsp_info_entry_helper_( + /* species_idx = */ SpLUT::H2O_ice_dust, + /* onlygrainsp_idx = */ OnlyGrainSpLUT::H2O_ice_dust, + /* name = */ "H2O_ice_dust", + /* h2dust_uses_carbonaceous_table = */ false, + /* sublimation_temperature = */ 153.0, + /* bulk_density_cgs = */ 0.92, + /* growth_ingredients = */ nullptr); + } + + return out; +} + +#undef GRIMPL_INGREDIENT_LIST_SENTINEL diff --git a/src/clib/dust/grain_species_info.hpp b/src/clib/dust/grain_species_info.hpp index 1ae31191a..1f20c52ea 100644 --- a/src/clib/dust/grain_species_info.hpp +++ b/src/clib/dust/grain_species_info.hpp @@ -6,17 +6,13 @@ //===----------------------------------------------------------------------===// /// /// @file -/// Defines details about the dust grain species +/// Declares GrainSpeciesInfo as well as other associated types and functions /// //===----------------------------------------------------------------------===// #ifndef GRAIN_SPECIES_INFO_HPP #define GRAIN_SPECIES_INFO_HPP -#include "LUT.hpp" -#include "grackle_macros.h" -#include // memcpy - namespace grackle::impl { /// holds information about a single gas species that is an ingredient for @@ -148,255 +144,11 @@ inline int get_n_grain_species(int dust_species_parameter) { /// The correctness of this constant is explicitly checked in a unit test inline constexpr int max_ingredients_per_grain_species = 3; -// define a few macros to help implement new_GrainSpeciesInfo -// - we undef all of these macros after they're used for cleanliness -// - we may want to consider an alternative that avoids macros in the future -#define GRIMPL_INGREDIENT_LIST_SENTINEL {-1, -1, -1.0} - -// does the heavy-lifting for implementing GR_MK_GRAIN_INFO_ENTRY -inline GrainSpeciesInfoEntry mk_gsp_info_entry_helper_( - int species_idx, int onlygrainsp_idx, const char* name, - bool h2dust_uses_carbonaceous_table, double sublimation_temperature, - double bulk_density_cgs, const GrainGrowthIngredient* growth_ingredients) { - GrainGrowthIngredient* out_ingredient_ptr = nullptr; - int n_ingredients = 0; - - if (growth_ingredients != nullptr) { - n_ingredients = 0; - // increment the counter until we reach the sentinel - while (growth_ingredients[n_ingredients].coef != -1) { - n_ingredients++; - } - // allocate and initialize the pointer - out_ingredient_ptr = new GrainGrowthIngredient[n_ingredients]; - std::memcpy(out_ingredient_ptr, growth_ingredients, - (std::size_t)n_ingredients); - } - - return GrainSpeciesInfoEntry{species_idx, - onlygrainsp_idx, - name, - h2dust_uses_carbonaceous_table, - sublimation_temperature, - bulk_density_cgs, - n_ingredients, - out_ingredient_ptr}; -} - -/// fills the specified out.species_infoay with the number of grain species. -/// -/// The out.species_infoay is expected to have space for the number of entries -/// returned by get_n_grain_species(dust_species_paramter) +/// Constructs an returns a fully initialized GrainSpeciesInfo instance. /// /// @param[in] dust_species_parameter The parameter tracked by #chemistry_data /// @returns A fully initialized GrainSpeciesInfo instance -inline GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter) { - GrainSpeciesInfo out{-1, nullptr}; // indicates an error - out.n_species = get_n_grain_species(dust_species_parameter); - - if (out.n_species <= 0) { - return out; - } - - out.species_info = new GrainSpeciesInfoEntry[out.n_species]; - - // At the time of writing: - // - we **only** use h2rate_carbonaceous_coef_table for the AC_dust - // species (amorphous carbon grains) - // - we use h2rate_silicate_coef_table for **ALL** other grain species - // (including species that are obviously NOT silicates) - // - // It is not obvious at all why this is the case... - // - // TODO: we should add documentation clearly addressing each of the - // following questions and replace this todo-item with a reference to the - // part of the documentation that discusses this! The documentation - // should make sure to address: - // - if this choice physically motivated or a common convention? (or if - // it was a quick & dirty choice because Gen didn't know what to do) - // - do we have a sense for how "wrong" the choice is? (e.g. will it - // generally overpredict/underpredict?) - // - why don't we use the generic my_rates->h2dusta rate as a generic - // fallback for non-silicate grains instead of using h2dustS? I realize - // that h2dusta has very different units... - - if (dust_species_parameter > 0) { - const GrainGrowthIngredient MgSiO3_dust_ingred[] = { - // {coef, species_idx, particle mass} - {1, SpLUT::Mg, 24.}, - {1, SpLUT::SiOI, 44.}, - {2, SpLUT::H2O, 18.}, - GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[0] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::MgSiO3_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::MgSiO3_dust, - /* name = */ "MgSiO3_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 1222.0, - /* bulk_density_cgs = */ 3.20185, - /* growth_ingredients = */ MgSiO3_dust_ingred); - - const GrainGrowthIngredient AC_dust_ingred[] = { - // {coef, species_idx, particle mass} - {1, SpLUT::CI, 12.}, - GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[1] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::AC_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::AC_dust, - /* name = */ "AC_dust", - /* h2dust_uses_carbonaceous_table = */ true, - /* sublimation_temperature = */ 1800.0, - /* bulk_density_cgs = */ 2.27949, - /* growth_ingredients = */ AC_dust_ingred); - } - - if (dust_species_parameter > 1) { - const GrainGrowthIngredient SiM_dust_ingred[] = { - // {coef, species_idx, particle mass} - {1, SpLUT::SiI, 28.}, - GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[2] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::SiM_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::SiM_dust, - /* name = */ "SiM_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 1500.0, - /* bulk_density_cgs = */ 2.34118, - /* growth_ingredients = */ SiM_dust_ingred); - - const GrainGrowthIngredient FeM_dust_ingred[] = { - // {coef, species_idx, particle mass} - {1, SpLUT::Fe, 56.}, - GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[3] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::FeM_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::FeM_dust, - /* name = */ "FeM_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 1500.0, - /* bulk_density_cgs = */ 7.95995, - /* growth_ingredients = */ FeM_dust_ingred); - - const GrainGrowthIngredient Mg2SiO4_dust_ingred[] = { - // {coef, species_idx, particle mass} - {2, SpLUT::Mg, 24.}, - {1, SpLUT::SiOI, 44.}, - {3, SpLUT::H2O, 18.}, - GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[4] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::Mg2SiO4_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::Mg2SiO4_dust, - /* name = */ "Mg2SiO4_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 1277.0, - /* bulk_density_cgs = */ 3.22133, - /* growth_ingredients = */ Mg2SiO4_dust_ingred); - - const GrainGrowthIngredient Fe3O4_dust_ingred[] = { - // {coef, species_idx, particle mass} - {3, SpLUT::Fe, 56.}, - {4, SpLUT::H2O, 18.}, - GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[5] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::Fe3O4_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::Fe3O4_dust, - /* name = */ "Fe3O4_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 1500.0, - /* bulk_density_cgs = */ 5.25096, - /* growth_ingredients = */ Fe3O4_dust_ingred); - - const GrainGrowthIngredient SiO2_dust_ingred[] = { - // {coef, species_idx, particle mass} - {1, SpLUT::SiO2I, 60.}, - GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[6] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::SiO2_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::SiO2_dust, - /* name = */ "SiO2_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 1500.0, - /* bulk_density_cgs = */ 2.66235, - /* growth_ingredients = */ SiO2_dust_ingred); - - const GrainGrowthIngredient MgO_dust_ingred[] = { - // {coef, species_idx, particle mass} - {1, SpLUT::Mg, 24.}, - {1, SpLUT::H2O, 18.}, - GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[7] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::MgO_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::MgO_dust, - /* name = */ "MgO_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 1500.0, - /* bulk_density_cgs = */ 3.58157, - /* growth_ingredients = */ MgO_dust_ingred); - - const GrainGrowthIngredient FeS_dust_ingred[] = { - // {coef, species_idx, particle mass} - {1, SpLUT::Fe, 56.}, - {1, SpLUT::S, 32.}, - GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[8] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::FeS_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::FeS_dust, - /* name = */ "FeS_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 680.0, - /* bulk_density_cgs = */ 4.87265, - /* growth_ingredients = */ FeS_dust_ingred); - - const GrainGrowthIngredient Al2O3_dust_ingred[] = { - // {coef, species_idx, particle mass} - {2, SpLUT::Al, 27.}, - {3, SpLUT::H2O, 18.}, - GRIMPL_INGREDIENT_LIST_SENTINEL}; - out.species_info[9] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::Al2O3_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::Al2O3_dust, - /* name = */ "Al2O3_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 1500.0, - /* bulk_density_cgs = */ 4.01610, - /* growth_ingredients = */ Al2O3_dust_ingred); - } - - if (dust_species_parameter > 2) { - // We do not consider the growth of refractory organics, volatile - // organics, and water ice because their sublimation temperatures - // are low (100-600 K). They sublimate before the growth occurs. - - out.species_info[10] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::ref_org_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::ref_org_dust, - /* name = */ "ref_org_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 575.0, - /* bulk_density_cgs = */ 1.5, - /* growth_ingredients = */ nullptr); - out.species_info[11] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::vol_org_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::vol_org_dust, - /* name = */ "vol_org_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 375.0, - /* bulk_density_cgs = */ 1.0, - /* growth_ingredients = */ nullptr); - out.species_info[12] = mk_gsp_info_entry_helper_( - /* species_idx = */ SpLUT::H2O_ice_dust, - /* onlygrainsp_idx = */ OnlyGrainSpLUT::H2O_ice_dust, - /* name = */ "H2O_ice_dust", - /* h2dust_uses_carbonaceous_table = */ false, - /* sublimation_temperature = */ 153.0, - /* bulk_density_cgs = */ 0.92, - /* growth_ingredients = */ nullptr); - } - - return out; -} - -#undef GRIMPL_INGREDIENT_LIST_SENTINEL +GrainSpeciesInfo new_GrainSpeciesInfo(int dust_species_parameter); /// performs cleanup of the contents of GrainSpeciesInfo /// From 3fdcdc8b14aaa8a44bf94a5efc933b4166102d24 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Sun, 7 Dec 2025 15:01:54 -0500 Subject: [PATCH 08/11] mv grain-growth-rxn comments I moved them from initialize_dust_yields.cpp where the comments were out of place to dust/grain_species_info.cpp --- src/clib/dust/grain_species_info.cpp | 17 +++++++++++++++++ src/clib/initialize_dust_yields.cpp | 19 ++----------------- 2 files changed, 19 insertions(+), 17 deletions(-) diff --git a/src/clib/dust/grain_species_info.cpp b/src/clib/dust/grain_species_info.cpp index e06a4d84c..3a6a76f80 100644 --- a/src/clib/dust/grain_species_info.cpp +++ b/src/clib/dust/grain_species_info.cpp @@ -119,6 +119,7 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( // that h2dusta has very different units... if (dust_species_parameter > 0) { + // growth rxn: "Mg + SiO + 2H2O -> MgSiO3_dust + 2H2I" const GrainGrowthIngredient MgSiO3_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::Mg, 24.}, @@ -134,6 +135,7 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( /* bulk_density_cgs = */ 3.20185, /* growth_ingredients = */ MgSiO3_dust_ingred); + // growth rxn: "C -> AC_dust" const GrainGrowthIngredient AC_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::CI, 12.}, @@ -149,6 +151,7 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( } if (dust_species_parameter > 1) { + // growth rxn: "Si -> SiM_dust" const GrainGrowthIngredient SiM_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::SiI, 28.}, @@ -162,6 +165,7 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( /* bulk_density_cgs = */ 2.34118, /* growth_ingredients = */ SiM_dust_ingred); + // growth rxn: "Fe -> FeM_dust" const GrainGrowthIngredient FeM_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::Fe, 56.}, @@ -175,6 +179,7 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( /* bulk_density_cgs = */ 7.95995, /* growth_ingredients = */ FeM_dust_ingred); + // growth rxn: "2Mg + SiO + 3H2O -> Mg2SiO4_dust + 3H2I" const GrainGrowthIngredient Mg2SiO4_dust_ingred[] = { // {coef, species_idx, particle mass} {2, SpLUT::Mg, 24.}, @@ -190,6 +195,7 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( /* bulk_density_cgs = */ 3.22133, /* growth_ingredients = */ Mg2SiO4_dust_ingred); + // growth rxn: "3Fe + 4H2O -> Fe3O4_dust + 4H2I" const GrainGrowthIngredient Fe3O4_dust_ingred[] = { // {coef, species_idx, particle mass} {3, SpLUT::Fe, 56.}, @@ -204,6 +210,7 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( /* bulk_density_cgs = */ 5.25096, /* growth_ingredients = */ Fe3O4_dust_ingred); + // growth rxn: "SiO2 -> SiO2_dust" const GrainGrowthIngredient SiO2_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::SiO2I, 60.}, @@ -217,6 +224,7 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( /* bulk_density_cgs = */ 2.66235, /* growth_ingredients = */ SiO2_dust_ingred); + // growth rxn: "Mg + H2O -> MgO_dust + H2I" const GrainGrowthIngredient MgO_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::Mg, 24.}, @@ -231,6 +239,7 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( /* bulk_density_cgs = */ 3.58157, /* growth_ingredients = */ MgO_dust_ingred); + // growth rxn: "Fe + S -> FeS_dust" const GrainGrowthIngredient FeS_dust_ingred[] = { // {coef, species_idx, particle mass} {1, SpLUT::Fe, 56.}, @@ -245,6 +254,7 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( /* bulk_density_cgs = */ 4.87265, /* growth_ingredients = */ FeS_dust_ingred); + // growth rxn: "2Al + 3H2O -> Al2O3_dust + 3H2I" const GrainGrowthIngredient Al2O3_dust_ingred[] = { // {coef, species_idx, particle mass} {2, SpLUT::Al, 27.}, @@ -265,6 +275,8 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( // organics, and water ice because their sublimation temperatures // are low (100-600 K). They sublimate before the growth occurs. + // nominal growth rxn: "0.5CO + 0.5CH2 + 1.2N -> ref_org_dust" + // nuclide ratios: C:H:O:N = 1:1:0.5:1.2 out.species_info[10] = mk_gsp_info_entry_helper_( /* species_idx = */ SpLUT::ref_org_dust, /* onlygrainsp_idx = */ OnlyGrainSpLUT::ref_org_dust, @@ -273,6 +285,9 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( /* sublimation_temperature = */ 575.0, /* bulk_density_cgs = */ 1.5, /* growth_ingredients = */ nullptr); + + // nominal growth rxn: "CO + 2H2I -> vol_org_dust" + // effective formula: CH3OH out.species_info[11] = mk_gsp_info_entry_helper_( /* species_idx = */ SpLUT::vol_org_dust, /* onlygrainsp_idx = */ OnlyGrainSpLUT::vol_org_dust, @@ -281,6 +296,8 @@ grackle::impl::GrainSpeciesInfo grackle::impl::new_GrainSpeciesInfo( /* sublimation_temperature = */ 375.0, /* bulk_density_cgs = */ 1.0, /* growth_ingredients = */ nullptr); + + // nominal growth rxn: "H2O -> H2O_ice_dust" out.species_info[12] = mk_gsp_info_entry_helper_( /* species_idx = */ SpLUT::H2O_ice_dust, /* onlygrainsp_idx = */ OnlyGrainSpLUT::H2O_ice_dust, diff --git a/src/clib/initialize_dust_yields.cpp b/src/clib/initialize_dust_yields.cpp index 8c0906d54..41f137ca5 100644 --- a/src/clib/initialize_dust_yields.cpp +++ b/src/clib/initialize_dust_yields.cpp @@ -43,24 +43,9 @@ int grackle::impl::initialize_dust_yields(chemistry_data *my_chemistry, code_units *my_units) { - if (my_chemistry->metal_chemistry == 0) + if (my_chemistry->metal_chemistry == 0) { return SUCCESS; - -//-------kdMgSiO3 : Mg + SiO + 2H2O -> MgSiO3 + 2H2I -//-------kdAC : C -> AC - -//-------kdSiM : Si -> SiM -//-------kdFeM : Fe -> FeM -//-------kdMg2SiO4 : 2Mg + SiO + 3H2O -> Mg2SiO4 + 3H2I -//-------kdFe3O4 : 3Fe + 4H2O -> Fe3O4 + 4H2I -//-------kdSiO2D : SiO2 -> SiO2D -//-------kdMgO : Mg + H2O -> MgO + H2I -//-------kdFeS : Fe + S -> FeS -//-------kdAl2O3 : 2Al + 3H2O -> Al2O3 + 3H2I - -//-------kdreforg : 0.5CO + 0.5CH2 + 1.2N -> reforg (C:H:O:N = 1:1:0.5:1.2) -//-------kdvolorg : CO + 2H2I -> volorg (CH3OH) -//-------kdH2Oice : H2O -> H2Oice + } int NSN, NTd, Nmom; double Td0, dTd; From 9d2b9cf4dd94ae09d8bb7a78f25e464188b03e5b Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 8 Dec 2025 12:17:32 -0500 Subject: [PATCH 09/11] convert interp_table_utils.{h->hpp} --- .clang-format-ignore | 2 +- src/clib/CMakeLists.txt | 2 +- src/clib/init_misc_species_cool_rates.cpp | 2 +- src/clib/initialize_chemistry_data.cpp | 12 +++--- src/clib/interp_table_utils.h | 43 --------------------- src/clib/interp_table_utils.hpp | 46 +++++++++++++++++++++++ 6 files changed, 56 insertions(+), 51 deletions(-) delete mode 100644 src/clib/interp_table_utils.h create mode 100644 src/clib/interp_table_utils.hpp diff --git a/.clang-format-ignore b/.clang-format-ignore index ab4b6b6f1..4cda5afa3 100644 --- a/.clang-format-ignore +++ b/.clang-format-ignore @@ -36,7 +36,7 @@ src/clib/initialize_dust_yields.cpp src/clib/initialize_rates.cpp src/clib/internal_types.hpp src/clib/internal_units.h -src/clib/interp_table_utils.h +src/clib/interp_table_utils.hpp src/clib/interpolate.hpp src/clib/phys_constants.h src/clib/rate_functions.c diff --git a/src/clib/CMakeLists.txt b/src/clib/CMakeLists.txt index c9ae9f141..07734e2d6 100644 --- a/src/clib/CMakeLists.txt +++ b/src/clib/CMakeLists.txt @@ -118,6 +118,7 @@ add_library(Grackle_Grackle initialize_dust_yields.cpp initialize_dust_yields.hpp initialize_rates.cpp initialize_rates.hpp internal_types.cpp internal_types.hpp + interp_table_utils.hpp lookup_cool_rates1d.hpp make_consistent.cpp make_consistent.hpp opaque_storage.hpp @@ -173,7 +174,6 @@ add_library(Grackle_Grackle grackle_chemistry_data_fields.def # <-- acts as a C header grackle_macros.h index_helper.h - interp_table_utils.h phys_constants.h collisional_rxn_rate_members.def # <-- acts as a C header utils.h diff --git a/src/clib/init_misc_species_cool_rates.cpp b/src/clib/init_misc_species_cool_rates.cpp index 83e1dd9c4..90393711a 100644 --- a/src/clib/init_misc_species_cool_rates.cpp +++ b/src/clib/init_misc_species_cool_rates.cpp @@ -17,7 +17,7 @@ #include "grackle.h" #include "grackle_macros.h" -#include "interp_table_utils.h" // free_interp_grid_ +#include "interp_table_utils.hpp" // free_interp_grid_ #include "init_misc_species_cool_rates.hpp" // forward declarations #include "internal_units.h" #include "phys_constants.h" diff --git a/src/clib/initialize_chemistry_data.cpp b/src/clib/initialize_chemistry_data.cpp index b8b39fabf..44d253348 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -19,7 +19,7 @@ #include "grackle.h" #include "grackle_macros.h" #include "auto_general.h" -#include "interp_table_utils.h" // free_interp_grid_ +#include "interp_table_utils.hpp" // free_interp_grid_ #include "init_misc_species_cool_rates.hpp" // free_misc_species_cool_rates #include "initialize_cloudy_data.h" #include "initialize_dust_yields.hpp" // free_dust_yields @@ -574,13 +574,13 @@ extern "C" int local_free_chemistry_data(chemistry_data *my_chemistry, GRACKLE_FREE(my_rates->gas_grain); GRACKLE_FREE(my_rates->gas_grain2); - free_interp_grid_(&my_rates->LH2); - free_interp_grid_(&my_rates->LHD); + grackle::impl::free_interp_grid_(&my_rates->LH2); + grackle::impl::free_interp_grid_(&my_rates->LHD); // we deal with freeing other interp grids inside of // free_misc_species_cool_rates - free_interp_grid_(&my_rates->alphap); + grackle::impl::free_interp_grid_(&my_rates->alphap); GRACKLE_FREE(my_rates->gr_N); @@ -647,7 +647,9 @@ extern "C" int local_free_chemistry_data(chemistry_data *my_chemistry, // delete contents of h2dust_grain_interp_props (automatically handles the // case where we didn't allocate anything) - free_interp_grid_props_(&my_rates->opaque_storage->h2dust_grain_interp_props); + grackle::impl::free_interp_grid_props_( + &my_rates->opaque_storage->h2dust_grain_interp_props, + /* use_delete = */ false); // since h2dust_grain_interp_props isn't a pointer, there is nothing more to // allocate right here diff --git a/src/clib/interp_table_utils.h b/src/clib/interp_table_utils.h deleted file mode 100644 index a9ad59147..000000000 --- a/src/clib/interp_table_utils.h +++ /dev/null @@ -1,43 +0,0 @@ -/*********************************************************************** -/ -/ Declare utility functions related to the interpolation tables -/ -/ -/ 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 INTERP_TABLE_UTILS_H -#define INTERP_TABLE_UTILS_H - -#include "grackle.h" // gr_interp_grid, GRACKLE_CLOUDY_TABLE_MAX_DIMENSION -#include "grackle_macros.h" // GRACKLE_FREE - -#ifdef __cplusplus -extern "C" { -#endif - -/// Free memory associated with a #gr_interp_grid_props instance -static inline void free_interp_grid_props_(gr_interp_grid_props* props) -{ - for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) { - GRACKLE_FREE(props->parameters[i]); - } -} - -/// Free memory associated with a #gr_interp_grid -static inline void free_interp_grid_(gr_interp_grid* grid) -{ - free_interp_grid_props_(&(grid->props)); - GRACKLE_FREE(grid->data); -} - -#ifdef __cplusplus -} // extern "C" -#endif - -#endif /* INTERP_TABLE_UTILS_H */ diff --git a/src/clib/interp_table_utils.hpp b/src/clib/interp_table_utils.hpp new file mode 100644 index 000000000..243d61876 --- /dev/null +++ b/src/clib/interp_table_utils.hpp @@ -0,0 +1,46 @@ +//===----------------------------------------------------------------------===// +// +// See the LICENSE file for license and copyright information +// SPDX-License-Identifier: NCSA AND BSD-3-Clause +// +//===----------------------------------------------------------------------===// +/// +/// @file +/// Declare utility functions related to the interpolation tables +/// +//===----------------------------------------------------------------------===// + +#ifndef INTERP_TABLE_UTILS_HPP +#define INTERP_TABLE_UTILS_HPP + +#include "grackle.h" // gr_interp_grid, GRACKLE_CLOUDY_TABLE_MAX_DIMENSION +#include "grackle_macros.h" // GRACKLE_FREE + +namespace grackle::impl { + +/// Free memory associated with a #gr_interp_grid_props instance +inline void free_interp_grid_props_(gr_interp_grid_props* props, + bool use_delete) +{ + for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) { + if (use_delete) { + if (props->parameters[i] != nullptr) { + delete[] props->parameters[i]; + props->parameters[i] = nullptr; + } + } else { + GRACKLE_FREE(props->parameters[i]); + } + } +} + +/// Free memory associated with a #gr_interp_grid +inline void free_interp_grid_(gr_interp_grid* grid) +{ + free_interp_grid_props_(&(grid->props), /* use_delete = */ false); + GRACKLE_FREE(grid->data); +} + +} // namespace grackle::impl + +#endif /* INTERP_TABLE_UTILS_HPP */ From d7727b33a6a666b3f68f297f915c64c7b15eb835 Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 8 Dec 2025 12:18:56 -0500 Subject: [PATCH 10/11] apply clang-format to src/clib/interp_table_utils.hpp --- .clang-format-ignore | 1 - src/clib/interp_table_utils.hpp | 10 ++++------ 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/.clang-format-ignore b/.clang-format-ignore index 4cda5afa3..d15536e08 100644 --- a/.clang-format-ignore +++ b/.clang-format-ignore @@ -36,7 +36,6 @@ src/clib/initialize_dust_yields.cpp src/clib/initialize_rates.cpp src/clib/internal_types.hpp src/clib/internal_units.h -src/clib/interp_table_utils.hpp src/clib/interpolate.hpp src/clib/phys_constants.h src/clib/rate_functions.c diff --git a/src/clib/interp_table_utils.hpp b/src/clib/interp_table_utils.hpp index 243d61876..c799c851e 100644 --- a/src/clib/interp_table_utils.hpp +++ b/src/clib/interp_table_utils.hpp @@ -13,15 +13,14 @@ #ifndef INTERP_TABLE_UTILS_HPP #define INTERP_TABLE_UTILS_HPP -#include "grackle.h" // gr_interp_grid, GRACKLE_CLOUDY_TABLE_MAX_DIMENSION -#include "grackle_macros.h" // GRACKLE_FREE +#include "grackle.h" // gr_interp_grid, GRACKLE_CLOUDY_TABLE_MAX_DIMENSION +#include "grackle_macros.h" // GRACKLE_FREE namespace grackle::impl { /// Free memory associated with a #gr_interp_grid_props instance inline void free_interp_grid_props_(gr_interp_grid_props* props, - bool use_delete) -{ + bool use_delete) { for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) { if (use_delete) { if (props->parameters[i] != nullptr) { @@ -35,8 +34,7 @@ inline void free_interp_grid_props_(gr_interp_grid_props* props, } /// Free memory associated with a #gr_interp_grid -inline void free_interp_grid_(gr_interp_grid* grid) -{ +inline void free_interp_grid_(gr_interp_grid* grid) { free_interp_grid_props_(&(grid->props), /* use_delete = */ false); GRACKLE_FREE(grid->data); } From 0cd4082fdd3104e9d3b28ef2c78c54af5b5fa52b Mon Sep 17 00:00:00 2001 From: Matthew Abruzzo Date: Mon, 8 Dec 2025 13:01:45 -0500 Subject: [PATCH 11/11] shift 2 functions out of src/clib/initialize_chemistry_data into interp_table_utils.hpp --- src/clib/initialize_chemistry_data.cpp | 39 ++++++++------------------ src/clib/interp_table_utils.hpp | 17 +++++++++++ 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/src/clib/initialize_chemistry_data.cpp b/src/clib/initialize_chemistry_data.cpp index 44d253348..e9089c782 100644 --- a/src/clib/initialize_chemistry_data.cpp +++ b/src/clib/initialize_chemistry_data.cpp @@ -19,7 +19,7 @@ #include "grackle.h" #include "grackle_macros.h" #include "auto_general.h" -#include "interp_table_utils.hpp" // free_interp_grid_ +#include "interp_table_utils.hpp" #include "init_misc_species_cool_rates.hpp" // free_misc_species_cool_rates #include "initialize_cloudy_data.h" #include "initialize_dust_yields.hpp" // free_dust_yields @@ -49,23 +49,6 @@ static void show_version(FILE *fp) fprintf (fp, "\n"); } -/// initialize an empty #gr_interp_grid_props -static void init_empty_interp_grid_props_(gr_interp_grid_props* props) { - props->rank = 0; - for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++){ - props->dimension[i] = 0; - props->parameters[i] = nullptr; - props->parameter_spacing[i] = 0.0; - } - props->data_size = 0; -} - -/// Initialize an empty #gr_interp_grid -static void initialize_empty_interp_grid_(gr_interp_grid* grid) -{ - init_empty_interp_grid_props_(&(grid->props)); - grid->data=NULL; -} /** * Initializes an empty #chemistry_data_storage struct with zeros and NULLs. @@ -152,18 +135,18 @@ static void initialize_empty_chemistry_data_storage_struct(chemistry_data_storag my_rates->cieY06 = NULL; - initialize_empty_interp_grid_(&my_rates->LH2); - initialize_empty_interp_grid_(&my_rates->LHD); + grackle::impl::initialize_empty_interp_grid_(&my_rates->LH2); + grackle::impl::initialize_empty_interp_grid_(&my_rates->LHD); - initialize_empty_interp_grid_(&my_rates->LCI); - initialize_empty_interp_grid_(&my_rates->LCII); - initialize_empty_interp_grid_(&my_rates->LOI); + grackle::impl::initialize_empty_interp_grid_(&my_rates->LCI); + grackle::impl::initialize_empty_interp_grid_(&my_rates->LCII); + grackle::impl::initialize_empty_interp_grid_(&my_rates->LOI); - initialize_empty_interp_grid_(&my_rates->LCO); - initialize_empty_interp_grid_(&my_rates->LOH); - initialize_empty_interp_grid_(&my_rates->LH2O); + grackle::impl::initialize_empty_interp_grid_(&my_rates->LCO); + grackle::impl::initialize_empty_interp_grid_(&my_rates->LOH); + grackle::impl::initialize_empty_interp_grid_(&my_rates->LH2O); - initialize_empty_interp_grid_(&my_rates->alphap); + grackle::impl::initialize_empty_interp_grid_(&my_rates->alphap); my_rates->gr_N = NULL; my_rates->gr_Size = 0; @@ -400,7 +383,7 @@ extern "C" int local_initialize_chemistry_data(chemistry_data *my_chemistry, my_rates->opaque_storage->kcol_rate_tables = nullptr; my_rates->opaque_storage->used_kcol_rate_indices = nullptr; my_rates->opaque_storage->n_kcol_rate_indices = 0; - init_empty_interp_grid_props_( + grackle::impl::init_empty_interp_grid_props_( &my_rates->opaque_storage->h2dust_grain_interp_props); my_rates->opaque_storage->grain_species_info = nullptr; diff --git a/src/clib/interp_table_utils.hpp b/src/clib/interp_table_utils.hpp index c799c851e..121767c09 100644 --- a/src/clib/interp_table_utils.hpp +++ b/src/clib/interp_table_utils.hpp @@ -18,6 +18,23 @@ namespace grackle::impl { +/// initialize an empty #gr_interp_grid_props +inline void init_empty_interp_grid_props_(gr_interp_grid_props* props) { + props->rank = 0; + for (int i = 0; i < GRACKLE_CLOUDY_TABLE_MAX_DIMENSION; i++) { + props->dimension[i] = 0; + props->parameters[i] = nullptr; + props->parameter_spacing[i] = 0.0; + } + props->data_size = 0; +} + +/// Initialize an empty #gr_interp_grid +inline void initialize_empty_interp_grid_(gr_interp_grid* grid) { + init_empty_interp_grid_props_(&(grid->props)); + grid->data = nullptr; +} + /// Free memory associated with a #gr_interp_grid_props instance inline void free_interp_grid_props_(gr_interp_grid_props* props, bool use_delete) {