From 8cfa067ad071515df703d869c8f002444bac7e1e Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sun, 7 Dec 2025 12:34:12 -0500 Subject: [PATCH] perf(GridFire) More preformance improvmnets 1. Switch to mimalloc which gave a roughly 10% improvment accross the board 2. Use much faster compososition caching 3. Reusing work vector --- build-config/meson.build | 4 + build-config/mimalloc/meson.build | 2 + meson_options.txt | 3 +- src/include/gridfire/engine/engine_abstract.h | 4 +- src/include/gridfire/engine/engine_graph.h | 8 +- .../gridfire/engine/views/engine_adaptive.h | 6 +- .../gridfire/engine/views/engine_defined.h | 4 +- .../gridfire/engine/views/engine_multiscale.h | 7 +- src/include/gridfire/reaction/reaction.h | 2 + .../solver/strategies/CVODE_solver_strategy.h | 7 +- src/include/gridfire/utils/hashing.h | 28 +++++-- src/lib/engine/engine_graph.cpp | 40 ++++++---- src/lib/engine/views/engine_adaptive.cpp | 14 +++- src/lib/engine/views/engine_defined.cpp | 2 +- src/lib/engine/views/engine_multiscale.cpp | 38 +++++---- src/lib/reaction/reaction.cpp | 17 +++- .../strategies/CVODE_solver_strategy.cpp | 79 ++++++++++++++----- src/meson.build | 4 + subprojects/fourdst.wrap | 2 +- subprojects/mimalloc.wrap | 13 +++ tests/graphnet_sandbox/main.cpp | 25 ++++-- validation/vv/GridFireValidationSuite.py | 37 ++++++--- validation/vv/testsuite.py | 57 +++++++++++-- 23 files changed, 306 insertions(+), 97 deletions(-) create mode 100644 build-config/mimalloc/meson.build create mode 100644 subprojects/mimalloc.wrap diff --git a/build-config/meson.build b/build-config/meson.build index a433637d..f4cc39f1 100644 --- a/build-config/meson.build +++ b/build-config/meson.build @@ -15,3 +15,7 @@ subdir('json') subdir('CLI11') + +if get_option('use_mimalloc') + subdir('mimalloc') +endif diff --git a/build-config/mimalloc/meson.build b/build-config/mimalloc/meson.build new file mode 100644 index 00000000..09819469 --- /dev/null +++ b/build-config/mimalloc/meson.build @@ -0,0 +1,2 @@ +mimalloc_proj = subproject('mimalloc') +mimalloc_dep = mimalloc_proj.get_variable('mi_dep') diff --git a/meson_options.txt b/meson_options.txt index cbee156c..88da0631 100644 --- a/meson_options.txt +++ b/meson_options.txt @@ -9,4 +9,5 @@ option('plugin_support', type: 'boolean', value: false, description: 'Enable sup option('python_target_version', type: 'string', value: '3.13', description: 'Target version for python compilation, only used for cross compilation') option('build_c_api', type: 'boolean', value: true, description: 'compile the C API') option('build_tools', type: 'boolean', value: true, description: 'build the GridFire command line tools') -option('openmp_support', type: 'boolean', value: false, description: 'Enable OpenMP support for parallelization') \ No newline at end of file +option('openmp_support', type: 'boolean', value: false, description: 'Enable OpenMP support for parallelization') +option('use_mimalloc', type: 'boolean', value: true, description: 'Use mimalloc as the memory allocator for GridFire. Generally this is ~10% faster than the system allocator.') \ No newline at end of file diff --git a/src/include/gridfire/engine/engine_abstract.h b/src/include/gridfire/engine/engine_abstract.h index e372afd9..4ee96930 100644 --- a/src/include/gridfire/engine/engine_abstract.h +++ b/src/include/gridfire/engine/engine_abstract.h @@ -144,6 +144,7 @@ namespace gridfire::engine { * @param comp Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. + * @param trust If true, indicates that the engine should trust the passed composition has already been collected. * @return expected> containing either dY/dt and energy generation rate or a stale engine * error indicating that the engine must be updated * @@ -154,7 +155,8 @@ namespace gridfire::engine { [[nodiscard]] virtual std::expected, EngineStatus> calculateRHSAndEnergy( const fourdst::composition::CompositionAbstract &comp, double T9, - double rho + double rho, + bool trust ) const = 0; }; diff --git a/src/include/gridfire/engine/engine_graph.h b/src/include/gridfire/engine/engine_graph.h index 4a1f2264..0f57ca5c 100644 --- a/src/include/gridfire/engine/engine_graph.h +++ b/src/include/gridfire/engine/engine_graph.h @@ -143,6 +143,7 @@ namespace gridfire::engine { * @param comp Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. + * @param trust * @return StepDerivatives containing dY/dt and energy generation rate. * * This method calculates the time derivatives of all species and the @@ -153,7 +154,8 @@ namespace gridfire::engine { [[nodiscard]] std::expected, engine::EngineStatus> calculateRHSAndEnergy( const fourdst::composition::CompositionAbstract &comp, double T9, - double rho + double rho, + bool trust ) const override; /** @@ -883,6 +885,8 @@ namespace gridfire::engine { mutable CppAD::sparse_jac_work m_jac_work; ///< Work object for sparse Jacobian calculations. mutable std::vector m_local_abundance_cache; mutable std::unordered_map> m_stepDerivativesCache; + mutable std::unordered_map, std::vector>> m_jacobianSubsetCache; + mutable std::unordered_map m_jacWorkCache; bool m_has_been_primed = false; ///< Flag indicating if the engine has been primed. @@ -895,7 +899,7 @@ namespace gridfire::engine { std::unique_ptr m_screeningModel = screening::selectScreeningModel(m_screeningType); bool m_usePrecomputation = true; ///< Flag to enable or disable using precomputed reactions for efficiency. Mathematically, this should not change the results. Generally end users should not need to change this. - bool m_useReverseReactions = true; ///< Flag to enable or disable reverse reactions. If false, only forward reactions are considered. + bool m_useReverseReactions = false; ///< Flag to enable or disable reverse reactions. If false, only forward reactions are considered. bool m_store_intermediate_reaction_contributions = false; ///< Flag to enable or disable storing intermediate reaction contributions for debugging. BuildDepthType m_depth; diff --git a/src/include/gridfire/engine/views/engine_adaptive.h b/src/include/gridfire/engine/views/engine_adaptive.h index 59d461cc..f0e8aac2 100644 --- a/src/include/gridfire/engine/views/engine_adaptive.h +++ b/src/include/gridfire/engine/views/engine_adaptive.h @@ -92,6 +92,7 @@ namespace gridfire::engine { * @param comp The current composition of the system. * @param T9 The temperature in units of 10^9 K. * @param rho The density in g/cm^3. + * @param trust * @return A StepDerivatives struct containing the derivatives of the active species and the * nuclear energy generation rate. * @@ -105,7 +106,8 @@ namespace gridfire::engine { [[nodiscard]] std::expected, engine::EngineStatus> calculateRHSAndEnergy( const fourdst::composition::CompositionAbstract &comp, double T9, - double rho + double rho, + bool trust ) const override; @@ -406,6 +408,8 @@ namespace gridfire::engine { /** @brief A flag indicating whether the view is stale and needs to be updated. */ bool m_isStale = true; + mutable std::unordered_map m_collected_composition_cache; + private: /** * @brief A struct to hold a reaction and its flow rate. diff --git a/src/include/gridfire/engine/views/engine_defined.h b/src/include/gridfire/engine/views/engine_defined.h index 16da8efd..9abadbbd 100644 --- a/src/include/gridfire/engine/views/engine_defined.h +++ b/src/include/gridfire/engine/views/engine_defined.h @@ -39,6 +39,7 @@ namespace gridfire::engine { * @param comp A Composition object containing the current composition of the system * @param T9 The temperature in units of 10^9 K. * @param rho The density in g/cm^3. + * @param trust * @return A StepDerivatives struct containing the derivatives of the active species and the * nuclear energy generation rate. * @@ -47,7 +48,8 @@ namespace gridfire::engine { [[nodiscard]] std::expected, engine::EngineStatus> calculateRHSAndEnergy( const fourdst::composition::CompositionAbstract &comp, double T9, - double rho + double rho, + bool trust ) const override; [[nodiscard]] EnergyDerivatives calculateEpsDerivatives( diff --git a/src/include/gridfire/engine/views/engine_multiscale.h b/src/include/gridfire/engine/views/engine_multiscale.h index 8565763e..8c6fed74 100644 --- a/src/include/gridfire/engine/views/engine_multiscale.h +++ b/src/include/gridfire/engine/views/engine_multiscale.h @@ -97,6 +97,7 @@ namespace gridfire::engine { * @param comp The current composition. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. + * @param trust * @return A `std::expected` containing `StepDerivatives` on success, or a * `StaleEngineError` if the engine's QSE cache does not contain a solution * for the given state. @@ -121,7 +122,8 @@ namespace gridfire::engine { [[nodiscard]] std::expected, engine::EngineStatus> calculateRHSAndEnergy( const fourdst::composition::CompositionAbstract &comp, double T9, - double rho + double rho, + bool trust ) const override; /** @@ -586,6 +588,7 @@ namespace gridfire::engine { * @param comp The input composition. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. + * @param trust * @return A new `Composition` object with algebraic species set to their equilibrium values. * * @par Purpose @@ -598,7 +601,7 @@ namespace gridfire::engine { * @pre The engine must have a valid QSE partition for the given state. * @throws StaleEngineError If the QSE cache misses. */ - fourdst::composition::Composition getNormalizedEquilibratedComposition(const fourdst::composition::CompositionAbstract& comp, double T9, double rho) const; + fourdst::composition::Composition getNormalizedEquilibratedComposition(const fourdst::composition::CompositionAbstract& comp, double T9, double rho, bool trust) const; /** * @brief Collect the composition from this and sub engines. diff --git a/src/include/gridfire/reaction/reaction.h b/src/include/gridfire/reaction/reaction.h index b6ce7345..c01f1621 100644 --- a/src/include/gridfire/reaction/reaction.h +++ b/src/include/gridfire/reaction/reaction.h @@ -637,6 +637,7 @@ namespace gridfire::reaction { mutable std::optional> m_reactantsVec; mutable std::optional> m_productsVec; + mutable std::optional m_hashCache = std::nullopt; ///< Cached hash value for the reaction. std::string m_sourceLabel; ///< Source label for the rate data (e.g., "wc12w", "st08"). RateCoefficientSet m_rateCoefficients; ///< The seven rate coefficients. @@ -1006,6 +1007,7 @@ namespace gridfire::reaction { std::string m_id; std::unordered_map m_reactionNameMap; ///< Maps reaction IDs to Reaction objects for quick lookup. std::unordered_set m_reactionHashes; + mutable std::optional m_hashCache = std::nullopt; }; diff --git a/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h b/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h index 42c41580..babcc5b7 100644 --- a/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h +++ b/src/include/gridfire/solver/strategies/CVODE_solver_strategy.h @@ -126,6 +126,7 @@ namespace gridfire::solver { * @brief Call to evaluate which will let the user control if the trigger reasoning is displayed * @param netIn Inputs: temperature [K], density [g cm^-3], tMax [s], composition. * @param displayTrigger Boolean flag to control if trigger reasoning is displayed + * @param forceReinitialize Boolean flag to force reinitialization of CVODE resources at the start * @return NetOut containing final Composition, accumulated energy [erg/g], step count, * and dEps/dT, dEps/dRho. * @throws std::runtime_error If any CVODE or SUNDIALS call fails (negative return codes), @@ -133,7 +134,7 @@ namespace gridfire::solver { * @throws exceptions::StaleEngineTrigger Propagated if the engine signals a stale state * during RHS evaluation (captured in the wrapper then rethrown here). */ - NetOut evaluate(const NetIn& netIn, bool displayTrigger); + NetOut evaluate(const NetIn& netIn, bool displayTrigger, bool forceReinitialize = false); /** * @brief Install a timestep callback. @@ -324,5 +325,9 @@ namespace gridfire::solver { std::optional m_relTol; ///< User-specified relative tolerance. bool m_detailed_step_logging = false; ///< If true, log detailed step diagnostics (error ratios, Jacobian, species balance). + + mutable size_t m_last_size = 0; + mutable size_t m_last_composition_hash = 0ULL; + mutable sunrealtype m_last_good_time_step = 0ULL; }; } \ No newline at end of file diff --git a/src/include/gridfire/utils/hashing.h b/src/include/gridfire/utils/hashing.h index 58922cf4..7017ab55 100644 --- a/src/include/gridfire/utils/hashing.h +++ b/src/include/gridfire/utils/hashing.h @@ -71,20 +71,32 @@ namespace gridfire::utils { return seed; } + inline std::size_t fast_mix(std::size_t h) noexcept { + h ^= h >> 33; + h *= 0xff51afd7ed558ccdULL; + h ^= h >> 33; + h *= 0xc4ceb9fe1a85ec53ULL; + h ^= h >> 33; + return h; + } + inline std::size_t hash_state( const fourdst::composition::CompositionAbstract& comp, const double T9, const double rho, const reaction::ReactionSet& reactions ) noexcept { - constexpr std::size_t seed = 0; - std::size_t comp_hash = fourdst::composition::utils::CompositionHash::hash_exact(comp); - for (const auto& reaction : reactions) { - comp_hash = hash_combine(comp_hash, hash_reaction(*reaction)); - } - std::size_t hash = hash_combine(seed, comp_hash); - hash = hash_combine(hash, std::bit_cast(T9)); - hash = hash_combine(hash, std::bit_cast(rho)); + std::size_t hash = comp.hash(); + const std::size_t topology_hash = reactions.hash(0); + + hash ^= topology_hash + 0x517cc1b727220a95 + (hash << 6) + (hash >> 2); + + const std::uint64_t t9_bits = std::bit_cast(T9); + const std::uint64_t rho_bits = std::bit_cast(rho); + + hash ^= fast_mix(t9_bits) + 0x9e3779b9 + (hash << 6) + (hash >> 2); + hash ^= fast_mix(rho_bits) + 0x9e3779b9 + (hash << 6) + (hash >> 2); + return hash; } } diff --git a/src/lib/engine/engine_graph.cpp b/src/lib/engine/engine_graph.cpp index 34cf7152..d29474bf 100644 --- a/src/lib/engine/engine_graph.cpp +++ b/src/lib/engine/engine_graph.cpp @@ -133,7 +133,8 @@ namespace gridfire::engine { std::expected, EngineStatus> GraphEngine::calculateRHSAndEnergy( const fourdst::composition::CompositionAbstract &comp, const double T9, - const double rho + const double rho, + bool trust ) const { return calculateRHSAndEnergy(comp, T9, rho, m_reactions); } @@ -744,6 +745,7 @@ namespace gridfire::engine { void GraphEngine::setUseReverseReactions(const bool useReverse) { m_useReverseReactions = useReverse; + syncInternalMaps(); } size_t GraphEngine::getSpeciesIndex(const fourdst::atomic::Species &species) const { @@ -1034,6 +1036,7 @@ namespace gridfire::engine { const double rho, const SparsityPattern &sparsityPattern ) const { + // --- Compute the intersection of the requested sparsity pattern with the full sparsity pattern --- SparsityPattern intersectionSparsityPattern; for (const auto& entry : sparsityPattern) { if (m_full_sparsity_set.contains(entry)) { @@ -1044,10 +1047,6 @@ namespace gridfire::engine { // --- Pack the input variables into a vector for CppAD --- const size_t numSpecies = m_networkSpecies.size(); std::vector x(numSpecies + 2, 0.0); - // const std::vector& Y_dynamic = comp.getMolarAbundanceVector(); - // for (size_t i = 0; i < numSpecies; ++i) { - // x[i] = Y_dynamic[i]; - // } size_t i = 0; for (const auto& species: m_networkSpecies) { double Yi = 0.0; // Small floor to avoid issues with zero abundances @@ -1075,18 +1074,25 @@ namespace gridfire::engine { const size_t num_cols_jac = numSpecies + 2; // +2 for T9 and rho CppAD::sparse_rc> CppAD_sparsity_pattern(num_rows_jac, num_cols_jac, nnz); + std::size_t sparsity_hash = 0; for (size_t k = 0; k < nnz; ++k) { + size_t local_intersection_hash = utils::hash_combine(intersectionSparsityPattern[k].first, intersectionSparsityPattern[k].second); + sparsity_hash = utils::hash_combine(sparsity_hash, local_intersection_hash); + CppAD_sparsity_pattern.set(k, intersectionSparsityPattern[k].first, intersectionSparsityPattern[k].second); } - CppAD::sparse_rcv, std::vector> jac_subset(CppAD_sparsity_pattern); - - // PERF: one of *the* most pressing things that needs to be done is remove the need for this call every - // time the jacobian is needed since coloring is expensive and we are throwing away the caching - // power of CppAD by clearing the work vector each time. We do this since we make a new subset every - // time. However, a better solution would be to make the subset stateful so it only changes if the requested - // sparsity pattern changes. This way we could reuse the work vector. - m_jac_work.clear(); + // --- Check cache for existing subset --- + if (!m_jacobianSubsetCache.contains(sparsity_hash)) { + m_jacobianSubsetCache.emplace(sparsity_hash, CppAD_sparsity_pattern); + m_jac_work.clear(); + } else { + if (m_jacWorkCache.contains(sparsity_hash)) { + m_jac_work.clear(); + m_jac_work = m_jacWorkCache.at(sparsity_hash); + } + } + auto& jac_subset = m_jacobianSubsetCache.at(sparsity_hash); m_rhsADFun.sparse_jac_rev( x, jac_subset, // Sparse Jacobian output @@ -1095,6 +1101,11 @@ namespace gridfire::engine { m_jac_work // Work vector for CppAD ); + // --- Stash the now populated work vector in the cache if not already present --- + if (!m_jacWorkCache.contains(sparsity_hash)) { + m_jacWorkCache.emplace(sparsity_hash, m_jac_work); + } + Eigen::SparseMatrix jacobianMatrix(numSpecies, numSpecies); std::vector > triplets; for (size_t k = 0; k < nnz; ++k) { @@ -1391,6 +1402,7 @@ namespace gridfire::engine { dependentVector.push_back(result.nuclearEnergyGenerationRate); m_rhsADFun.Dependent(adInput, dependentVector); + m_rhsADFun.optimize(); LOG_TRACE_L1(m_logger, "AD tape recorded successfully for the RHS and Eps calculation. Number of independent variables: {}.", adInput.size()); } @@ -1400,7 +1412,7 @@ namespace gridfire::engine { m_atomicReverseRates.reserve(m_reactions.size()); for (const auto& reaction: m_reactions) { - if (reaction->qValue() != 0.0) { + if (reaction->qValue() != 0.0 and m_useReverseReactions) { m_atomicReverseRates.push_back(std::make_unique(*reaction, *this)); } else { m_atomicReverseRates.push_back(nullptr); diff --git a/src/lib/engine/views/engine_adaptive.cpp b/src/lib/engine/views/engine_adaptive.cpp index a1956d3e..d8eb27c5 100644 --- a/src/lib/engine/views/engine_adaptive.cpp +++ b/src/lib/engine/views/engine_adaptive.cpp @@ -7,6 +7,7 @@ #include "gridfire/types/types.h" #include "gridfire/exceptions/error_engine.h" +#include "gridfire/utils/hashing.h" #include "quill/LogMacros.h" #include "quill/Logger.h" @@ -80,7 +81,7 @@ namespace gridfire::engine { std::expected, EngineStatus> AdaptiveEngineView::calculateRHSAndEnergy( const fourdst::composition::CompositionAbstract &comp, const double T9, - const double rho + const double rho, bool trust ) const { LOG_TRACE_L2(m_logger, "Calculating RHS and Energy in AdaptiveEngineView at T9 = {}, rho = {}.", T9, rho); validateState(); @@ -99,7 +100,14 @@ namespace gridfire::engine { } return ss.str(); }()); - fourdst::composition::Composition collectedComp = collectComposition(comp, T9, rho); + fourdst::composition::Composition collectedComp; + std::size_t state_hash = utils::hash_state(comp, T9, rho, m_activeReactions); + if (m_collected_composition_cache.contains(state_hash)) { + collectedComp = m_collected_composition_cache.at(state_hash); + } else { + collectedComp = collectComposition(comp, T9, rho); + m_collected_composition_cache[state_hash] = collectedComp; + } LOG_TRACE_L2( m_logger, "Composition Collected prior to passing to base engine. Collected Composition: {}", @@ -118,7 +126,7 @@ namespace gridfire::engine { } return ss.str(); }()); - auto result = m_baseEngine.calculateRHSAndEnergy(collectedComp, T9, rho); + auto result = m_baseEngine.calculateRHSAndEnergy(collectedComp, T9, rho, true); LOG_TRACE_L2(m_logger, "Base engine calculation of RHS and Energy complete."); if (!result) { diff --git a/src/lib/engine/views/engine_defined.cpp b/src/lib/engine/views/engine_defined.cpp index fd318b3f..69b1624c 100644 --- a/src/lib/engine/views/engine_defined.cpp +++ b/src/lib/engine/views/engine_defined.cpp @@ -43,7 +43,7 @@ namespace gridfire::engine { std::expected, EngineStatus> DefinedEngineView::calculateRHSAndEnergy( const fourdst::composition::CompositionAbstract &comp, const double T9, - const double rho + const double rho, bool trust ) const { validateNetworkState(); diff --git a/src/lib/engine/views/engine_multiscale.cpp b/src/lib/engine/views/engine_multiscale.cpp index 3bab8b9c..e3a55f63 100644 --- a/src/lib/engine/views/engine_multiscale.cpp +++ b/src/lib/engine/views/engine_multiscale.cpp @@ -216,7 +216,8 @@ namespace gridfire::engine { std::expected, EngineStatus> MultiscalePartitioningEngineView::calculateRHSAndEnergy( const fourdst::composition::CompositionAbstract &comp, const double T9, - const double rho + const double rho, + bool trust ) const { LOG_TRACE_L2(m_logger, "Calculating RHS and Energy in MultiscalePartitioningEngineView at T9 = {}, rho = {}.", T9, rho); LOG_TRACE_L2(m_logger, "Input composition is {}", [&comp]() -> std::string { @@ -231,7 +232,8 @@ namespace gridfire::engine { } return ss.str(); }()); - const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho); + // TODO: Figure out why setting trust -> trust causes issues. The only place I think I am setting that to true is in AdaptiveEngineView which has just called getNormalizedEquilibratedComposition... + const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho, false); LOG_TRACE_L2(m_logger, "Equilibrated composition prior to calling base engine is {}", [&qseComposition, &comp]() -> std::string { std::stringstream ss; size_t i = 0; @@ -248,7 +250,7 @@ namespace gridfire::engine { return ss.str(); }()); - const auto result = m_baseEngine.calculateRHSAndEnergy(qseComposition, T9, rho); + const auto result = m_baseEngine.calculateRHSAndEnergy(qseComposition, T9, rho, false); LOG_TRACE_L2(m_logger, "Base engine calculation of RHS and Energy complete."); if (!result) { @@ -271,7 +273,7 @@ namespace gridfire::engine { const double T9, const double rho ) const { - const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho); + const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho, false); return m_baseEngine.calculateEpsDerivatives(qseComposition, T9, rho); } @@ -280,7 +282,7 @@ namespace gridfire::engine { const double T9, const double rho ) const { - const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho); + const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho, false); return m_baseEngine.generateJacobianMatrix(qseComposition, T9, rho, m_dynamic_species); } @@ -318,7 +320,7 @@ namespace gridfire::engine { } } - const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho); + const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho, false); return m_baseEngine.generateJacobianMatrix(qseComposition, T9, rho, dynamicActiveSpeciesIntersection); } @@ -329,7 +331,7 @@ namespace gridfire::engine { const double rho, const SparsityPattern &sparsityPattern ) const { - const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho); + const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho, false); return m_baseEngine.generateJacobianMatrix(qseComposition, T9, rho, sparsityPattern); } @@ -350,7 +352,7 @@ namespace gridfire::engine { const double T9, const double rho ) const { - const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho); + const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho, false); return m_baseEngine.calculateMolarReactionFlow(reaction, qseComposition, T9, rho); } @@ -369,7 +371,7 @@ namespace gridfire::engine { const double T9, const double rho ) const { - const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho); + const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho, false); const auto result = m_baseEngine.getSpeciesTimescales(qseComposition, T9, rho); if (!result) { return std::unexpected{result.error()}; @@ -386,7 +388,7 @@ namespace gridfire::engine { const double T9, const double rho ) const { - const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho); + const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho, false); const auto result = m_baseEngine.getSpeciesDestructionTimescales(qseComposition, T9, rho); if (!result) { return std::unexpected{result.error()}; @@ -803,7 +805,7 @@ namespace gridfire::engine { LOG_TRACE_L1(m_logger, "{} QSE solvers created.", m_qse_solvers.size()); LOG_TRACE_L1(m_logger, "Calculating final equilibrated composition..."); - fourdst::composition::Composition result = getNormalizedEquilibratedComposition(comp, T9, rho); + fourdst::composition::Composition result = getNormalizedEquilibratedComposition(comp, T9, rho, false); LOG_TRACE_L1(m_logger, "Final equilibrated composition calculated..."); return result; @@ -836,7 +838,7 @@ namespace gridfire::engine { } } - const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho); + const fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(comp, T9, rho, false); // Calculate reaction flows and find min/max for logarithmic scaling of transparency std::vector reaction_flows; reaction_flows.reserve(all_reactions.size()); @@ -1072,8 +1074,12 @@ namespace gridfire::engine { fourdst::composition::Composition MultiscalePartitioningEngineView::getNormalizedEquilibratedComposition( const fourdst::composition::CompositionAbstract& comp, const double T9, - const double rho + const double rho, + const bool trust ) const { + if (trust) { + return fourdst::composition::Composition(comp); + } // Caching mechanism to avoid redundant QSE solves const std::array hashes = { fourdst::composition::utils::CompositionHash::hash_exact(comp), @@ -1108,7 +1114,7 @@ namespace gridfire::engine { ) const { const fourdst::composition::Composition result = m_baseEngine.collectComposition(comp, T9, rho); - fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(result, T9, rho); + fourdst::composition::Composition qseComposition = getNormalizedEquilibratedComposition(result, T9, rho, false); return qseComposition; } @@ -1893,7 +1899,7 @@ namespace gridfire::engine { scale_data[i] = 1.0 / Y; } - auto initial_rhs = m_engine.calculateRHSAndEnergy(result, T9, rho); + auto initial_rhs = m_engine.calculateRHSAndEnergy(result, T9, rho, false); if (!initial_rhs) { throw std::runtime_error("In QSE solver failed to calculate initial RHS"); } @@ -2068,7 +2074,7 @@ namespace gridfire::engine { data->comp.setMolarAbundance(species, y_data[index]); } - const auto result = data->engine.calculateRHSAndEnergy(data->comp, data->T9, data->rho); + const auto result = data->engine.calculateRHSAndEnergy(data->comp, data->T9, data->rho, false); if (!result) { return 1; // Potentially recoverable error diff --git a/src/lib/reaction/reaction.cpp b/src/lib/reaction/reaction.cpp index 2cb9ef2e..f1ab34ff 100644 --- a/src/lib/reaction/reaction.cpp +++ b/src/lib/reaction/reaction.cpp @@ -185,7 +185,11 @@ namespace gridfire::reaction { } uint64_t ReaclibReaction::hash(const uint64_t seed) const { - return XXHash64::hash(m_id.data(), m_id.size(), seed); + if (m_hashCache.has_value()) { + return m_hashCache.value(); + } + m_hashCache = XXHash64::hash(m_id.data(), m_id.size(), seed); + return m_hashCache.value(); } std::unique_ptr ReaclibReaction::clone() const { @@ -416,6 +420,7 @@ namespace gridfire::reaction { std::swap(m_reactions, temp.m_reactions); std::swap(m_reactionNameMap, temp.m_reactionNameMap); } + m_hashCache = std::nullopt; return *this; } @@ -430,6 +435,7 @@ namespace gridfire::reaction { m_reactionNameMap.emplace(std::move(reaction_id), new_index); m_reactionHashes.insert(reaction.hash(0)); + m_hashCache = std::nullopt; } void ReactionSet::add_reaction(std::unique_ptr&& reaction) { @@ -445,6 +451,7 @@ namespace gridfire::reaction { m_reactionNameMap.emplace(std::move(reaction_id), new_index); m_reactionHashes.insert(reaction_hash); + m_hashCache = std::nullopt; } void ReactionSet::extend(const ReactionSet &other) { @@ -482,6 +489,7 @@ namespace gridfire::reaction { } m_reactionHashes.erase(rh); + m_hashCache = std::nullopt; } bool ReactionSet::contains(const std::string_view& id) const { @@ -496,6 +504,7 @@ namespace gridfire::reaction { void ReactionSet::clear() { m_reactions.clear(); m_reactionNameMap.clear(); + m_hashCache = std::nullopt; } bool ReactionSet::contains_species(const Species& species) const { @@ -554,6 +563,9 @@ namespace gridfire::reaction { } uint64_t ReactionSet::hash(const uint64_t seed) const { + if (m_hashCache.has_value()) { + return m_hashCache.value(); + } if (m_reactions.empty()) { return XXHash64::hash(nullptr, 0, seed); } @@ -567,7 +579,8 @@ namespace gridfire::reaction { const auto data = static_cast(individualReactionHashes.data()); const size_t sizeInBytes = individualReactionHashes.size() * sizeof(uint64_t); - return XXHash64::hash(data, sizeInBytes, seed); + m_hashCache = XXHash64::hash(data, sizeInBytes, seed); + return m_hashCache.value(); } std::unordered_set ReactionSet::getReactionSetSpecies() const { diff --git a/src/lib/solver/strategies/CVODE_solver_strategy.cpp b/src/lib/solver/strategies/CVODE_solver_strategy.cpp index 01e2727e..b9e3c412 100644 --- a/src/lib/solver/strategies/CVODE_solver_strategy.cpp +++ b/src/lib/solver/strategies/CVODE_solver_strategy.cpp @@ -97,7 +97,8 @@ namespace gridfire::solver { NetOut CVODESolverStrategy::evaluate( const NetIn &netIn, - bool displayTrigger + bool displayTrigger, + bool forceReinitialize ) { LOG_TRACE_L1(m_logger, "Starting solver evaluation with T9: {} and rho: {}", netIn.temperature/1e9, netIn.density); LOG_TRACE_L1(m_logger, "Building engine update trigger...."); @@ -122,20 +123,54 @@ namespace gridfire::solver { relTol = *m_relTol; } - LOG_TRACE_L1(m_logger, "Starting engine update chain..."); - fourdst::composition::Composition equilibratedComposition = m_engine.update(netIn); - LOG_TRACE_L1(m_logger, "Engine updated and equilibrated composition found!"); + bool resourcesExist = (m_cvode_mem != nullptr) && (m_Y != nullptr); + + bool inconsistentComposition = netIn.composition.hash() != m_last_composition_hash; + fourdst::composition::Composition equilibratedComposition; + + if (forceReinitialize || !resourcesExist || inconsistentComposition) { + cleanup_cvode_resources(true); + LOG_INFO( + m_logger, + "Preforming full CVODE initialization (Reason: {})", + forceReinitialize ? "Forced reinitialization" : + (!resourcesExist ? "CVODE resources do not exist" : + "Input composition inconsistent with previous state")); + LOG_TRACE_L1(m_logger, "Starting engine update chain..."); + equilibratedComposition = m_engine.update(netIn); + LOG_TRACE_L1(m_logger, "Engine updated and equilibrated composition found!"); + + size_t numSpecies = m_engine.getNetworkSpecies().size(); + uint64_t N = numSpecies + 1; + + LOG_TRACE_L1(m_logger, "Number of species: {} ({} independent variables)", numSpecies, N); + LOG_TRACE_L1(m_logger, "Initializing CVODE resources"); + m_cvode_mem = CVodeCreate(CV_BDF, m_sun_ctx); + utils::check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate"); + + initialize_cvode_integration_resources(N, numSpecies, 0.0, equilibratedComposition, absTol, relTol, 0.0); + m_last_size = N; + } else { + LOG_INFO(m_logger, "Reusing existing CVODE resources (size: {})", m_last_size); + + const size_t numSpecies = m_engine.getNetworkSpecies().size(); + sunrealtype *y_data = N_VGetArrayPointer(m_Y); + for (size_t i = 0; i < numSpecies; i++) { + const auto& species = m_engine.getNetworkSpecies()[i]; + if (netIn.composition.contains(species)) { + y_data[i] = netIn.composition.getMolarAbundance(species); + } else { + y_data[i] = std::numeric_limits::min(); + } + } + y_data[numSpecies] = 0.0; // Reset energy accumulator + utils::check_cvode_flag(CVodeSStolerances(m_cvode_mem, relTol, absTol), "CVodeSStolerances"); + utils::check_cvode_flag(CVodeReInit(m_cvode_mem, 0.0, m_Y), "CVodeReInit"); + + equilibratedComposition = netIn.composition; // Use the provided composition as-is if we already have validated CVODE resources and that the composition is consistent with the previous state + } size_t numSpecies = m_engine.getNetworkSpecies().size(); - uint64_t N = numSpecies + 1; - - LOG_TRACE_L1(m_logger, "Number of species: {} ({} independent variables)", numSpecies, N); - LOG_TRACE_L1(m_logger, "Initializing CVODE resources"); - m_cvode_mem = CVodeCreate(CV_BDF, m_sun_ctx); - utils::check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate"); - - initialize_cvode_integration_resources(N, numSpecies, 0.0, equilibratedComposition, absTol, relTol, 0.0); - CVODEUserData user_data; user_data.solver_instance = this; user_data.engine = &m_engine; @@ -217,16 +252,16 @@ namespace gridfire::solver { postStep.setMolarAbundance(species, y_data[i]); } } - fourdst::composition::Composition collectedComposition = m_engine.collectComposition(postStep, netIn.temperature/1e9, netIn.density); - for (size_t i = 0; i < numSpecies; ++i) { - y_data[i] = collectedComposition.getMolarAbundance(m_engine.getNetworkSpecies()[i]); - } + // fourdst::composition::Composition collectedComposition = m_engine.collectComposition(postStep, netIn.temperature/1e9, netIn.density); + // for (size_t i = 0; i < numSpecies; ++i) { + // y_data[i] = collectedComposition.getMolarAbundance(m_engine.getNetworkSpecies()[i]); + // } LOG_INFO(m_logger, "Completed {:5} steps to time {:10.4E} [s] (dt = {:15.6E} [s]). Current specific energy: {:15.6E} [erg/g]", total_steps + n_steps, current_time, last_step_size, current_energy); LOG_DEBUG(m_logger, "Current composition (molar abundance): {}", [&]() -> std::string { std::stringstream ss; for (size_t i = 0; i < numSpecies; ++i) { const auto& species = m_engine.getNetworkSpecies()[i]; - ss << species.name() << ": (y_data = " << y_data[i] << ", collected = " << collectedComposition.getMolarAbundance(species) << ")"; + ss << species.name() << ": (y_data = " << y_data[i] << ", collected = " << postStep.getMolarAbundance(species) << ")"; if (i < numSpecies - 1) { ss << ", "; } @@ -428,7 +463,7 @@ namespace gridfire::solver { ); numSpecies = m_engine.getNetworkSpecies().size(); - N = numSpecies + 1; + size_t N = numSpecies + 1; LOG_INFO(m_logger, "Starting CVODE reinitialization after engine update..."); cleanup_cvode_resources(true); @@ -516,7 +551,9 @@ namespace gridfire::solver { LOG_TRACE_L2(m_logger, "Output data built!"); LOG_TRACE_L2(m_logger, "Solver evaluation complete!."); - + m_last_composition_hash = netOut.composition.hash(); + m_last_size = netOut.composition.size() + 1; + CVodeGetLastStep(m_cvode_mem, &m_last_good_time_step); return netOut; } @@ -730,7 +767,7 @@ namespace gridfire::solver { fourdst::composition::Composition composition(m_engine.getNetworkSpecies(), y_vec); LOG_TRACE_L2(m_logger, "Calculating RHS at time {} with {} species in composition", t, composition.size()); - const auto result = m_engine.calculateRHSAndEnergy(composition, data->T9, data->rho); + const auto result = m_engine.calculateRHSAndEnergy(composition, data->T9, data->rho, false); if (!result) { LOG_CRITICAL(m_logger, "Failed to calculate RHS at time {}: {}", t, EngineStatus_to_string(result.error())); throw exceptions::BadRHSEngineError(std::format("Failed to calculate RHS at time {}: {}", t, EngineStatus_to_string(result.error()))); diff --git a/src/meson.build b/src/meson.build index 4fa4fcea..17ea2b21 100644 --- a/src/meson.build +++ b/src/meson.build @@ -42,6 +42,10 @@ gridfire_build_dependencies = [ json_dep, ] +if get_option('use_mimalloc') + gridfire_build_dependencies += [mimalloc_dep] +endif + if get_option('plugin_support') gridfire_build_dependencies += [plugin_dep] endif diff --git a/subprojects/fourdst.wrap b/subprojects/fourdst.wrap index a805e371..bce058af 100644 --- a/subprojects/fourdst.wrap +++ b/subprojects/fourdst.wrap @@ -1,4 +1,4 @@ [wrap-git] url = https://github.com/4D-STAR/fourdst -revision = v0.9.14 +revision = v0.9.16 depth = 1 diff --git a/subprojects/mimalloc.wrap b/subprojects/mimalloc.wrap new file mode 100644 index 00000000..555cdc7e --- /dev/null +++ b/subprojects/mimalloc.wrap @@ -0,0 +1,13 @@ +[wrap-file] +directory = mimalloc-3.1.5 +source_url = https://github.com/microsoft/mimalloc/archive/refs/tags/v3.1.5.tar.gz +source_filename = mimalloc-3.1.5.tar.gz +source_hash = 1c6949032069d5ebea438ec5cedd602d06f40a92ddf0f0d9dcff0993e5f6635c +patch_filename = mimalloc_3.1.5-1_patch.zip +patch_url = https://wrapdb.mesonbuild.com/v2/mimalloc_3.1.5-1/get_patch +patch_hash = 321b4507c1adda5b7aa9954a5f1748e17bf30a11142b4f6c3d52929523565e80 +source_fallback_url = https://github.com/mesonbuild/wrapdb/releases/download/mimalloc_3.1.5-1/mimalloc-3.1.5.tar.gz +wrapdb_version = 3.1.5-1 + +[provide] +mimalloc = mi_dep diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp index f42ffa67..f44f5cb0 100644 --- a/tests/graphnet_sandbox/main.cpp +++ b/tests/graphnet_sandbox/main.cpp @@ -1,5 +1,7 @@ #include #include +#include +#include #include "gridfire/gridfire.h" @@ -228,9 +230,10 @@ int main(int argc, char** argv) { CLI::App app{"GridFire Sandbox Application."}; + constexpr size_t breaks = 100; double temp = 1.5e7; double rho = 1.5e2; - double tMax = 3.1536e+17; + double tMax = 3.1536e+17/breaks; app.add_option("-t,--temp", temp, "Temperature in K (Default 1.5e7K)"); app.add_option("-r,--rho", rho, "Density in g/cm^3 (Default 1.5e2g/cm^3)"); @@ -238,17 +241,29 @@ int main(int argc, char** argv) { CLI11_PARSE(app, argc, argv); - const NetIn netIn = init(temp, rho, tMax); + NetIn netIn = init(temp, rho, tMax); policy::MainSequencePolicy stellarPolicy(netIn.composition); stellarPolicy.construct(); engine::DynamicEngine& engine = stellarPolicy.construct(); solver::CVODESolverStrategy solver(engine); - solver.set_callback(solver::CVODESolverStrategy::TimestepCallback(callback_main)); + solver.set_stdout_logging_enabled(false); + // solver.set_callback(solver::CVODESolverStrategy::TimestepCallback(callback_main)); + + fourdst::composition::Composition reinputComp = netIn.composition; + NetOut netOut; + const auto timer = std::chrono::high_resolution_clock::now(); + for (int i = 0; i < breaks; ++i) { + NetIn in({.composition = reinputComp, .temperature = temp, .density = rho, .tMax = tMax, .dt0 = 1e-12}); + netOut = solver.evaluate(in, false, false); + reinputComp = netOut.composition; + } + const auto duration = std::chrono::duration_cast(std::chrono::high_resolution_clock::now() - timer).count(); + std::cout << "Average execution time over run: " << duration/breaks << " ms" << std::endl; + std::cout << "Total execution time over " << breaks << " runs: " << duration << " ms" << std::endl; - const NetOut netOut = solver.evaluate(netIn, false); log_results(netOut, netIn); - log_callback_data(temp); + // log_callback_data(temp); } diff --git a/validation/vv/GridFireValidationSuite.py b/validation/vv/GridFireValidationSuite.py index 7c31d986..2a81672a 100644 --- a/validation/vv/GridFireValidationSuite.py +++ b/validation/vv/GridFireValidationSuite.py @@ -1,7 +1,8 @@ from gridfire.policy import MainSequencePolicy, NetworkPolicy -from gridfire.engine import DynamicEngine, GraphEngine +from gridfire.engine import DynamicEngine, GraphEngine, EngineTypes from gridfire.type import NetIn +from typing import Dict from fourdst.composition import Composition from testsuite import TestSuite @@ -9,6 +10,12 @@ from utils import init_netIn, init_composition, years_to_seconds from enum import Enum +EngineNameToType: Dict[str, EngineTypes] = { + "graphengine": EngineTypes.GRAPH_ENGINE, + "multiscalepartitioningengineview": EngineTypes.MULTISCALE_PARTITIONING_ENGINE_VIEW, + "adaptiveengineview": EngineTypes.ADAPTIVE_ENGINE_VIEW +} + class SolarLikeStar_QSE_Suite(TestSuite): def __init__(self): initialComposition : Composition = init_composition() @@ -22,11 +29,11 @@ class SolarLikeStar_QSE_Suite(TestSuite): notes="Thermodynamically Static, MultiscalePartitioning Engine View" ) - def __call__(self): + def __call__(self, pynucastro_compare: bool = False, pync_engine: str = "AdaptiveEngineView"): policy : MainSequencePolicy = MainSequencePolicy(self.composition) engine : DynamicEngine = policy.construct() netIn : NetIn = init_netIn(self.temperature, self.density, self.tMax, self.composition) - self.evolve(engine, netIn) + self.evolve(engine, netIn, pynucastro_compare = pynucastro_compare, engine_type=EngineNameToType[pync_engine.lower()]) class MetalEnhancedSolarLikeStar_QSE_Suite(TestSuite): def __init__(self): @@ -41,7 +48,7 @@ class MetalEnhancedSolarLikeStar_QSE_Suite(TestSuite): notes="Thermodynamically Static, MultiscalePartitioning Engine View, Z enhanced by 1 dex, temperature reduced to 80% of solar core" ) - def __call__(self): + def __call__(self, pynucastro_compare: bool = False, pync_engine: str = "AdaptiveEngineView"): policy : MainSequencePolicy = MainSequencePolicy(self.composition) engine : GraphEngine = policy.construct() netIn : NetIn = init_netIn(self.temperature, self.density, self.tMax, self.composition) @@ -59,7 +66,7 @@ class MetalDepletedSolarLikeStar_QSE_Suite(TestSuite): notes="Thermodynamically Static, MultiscalePartitioning Engine View, Z depleted by 1 dex, temperature increased to 120% of solar core" ) - def __call__(self): + def __call__(self, pynucastro_compare: bool = False, pync_engine: str = "AdaptiveEngineView"): policy : MainSequencePolicy = MainSequencePolicy(self.composition) engine : GraphEngine = policy.construct() netIn : NetIn = init_netIn(self.temperature, self.density, self.tMax, self.composition) @@ -78,7 +85,7 @@ class SolarLikeStar_No_QSE_Suite(TestSuite): notes="Thermodynamically Static, No MultiscalePartitioning Engine View" ) - def __call__(self): + def __call__(self, pynucastro_compare: bool = False, pync_engine: str = "AdaptiveEngineView"): engine : GraphEngine = GraphEngine(self.composition, 3) netIn : NetIn = init_netIn(self.temperature, self.density, self.tMax, self.composition) self.evolve(engine, netIn) @@ -94,9 +101,19 @@ if __name__ == "__main__": import argparse parser = argparse.ArgumentParser(description="Run some subset of the GridFire validation suite.") parser.add_argument('--suite', type=str, choices=[suite.name for suite in ValidationSuites], nargs="+", help="The validation suite to run.") + parser.add_argument("--all", action="store_true", help="Run all validation suites.") + parser.add_argument("--pynucastro-compare", action="store_true", help="Generate pynucastro comparison data.") + parser.add_argument("--pync-engine", type=str, choices=["GraphEngine", "MultiscalePartitioningEngineView", "AdaptiveEngineView"], default="AdaptiveEngineView", help="The GridFire engine to use to select the reactions for pyuncastro comparison.") args = parser.parse_args() - for suite_name in args.suite: - suite = ValidationSuites[suite_name] - instance : TestSuite = suite.value() - instance() + if args.all: + for suite in ValidationSuites: + instance : TestSuite = suite.value() + instance(args.pynucastro_compare, args.pync_engine) + + else: + for suite_name in args.suite: + suite = ValidationSuites[suite_name] + instance : TestSuite = suite.value() + instance(args.pynucastro_compare, args.pync_engine) + diff --git a/validation/vv/testsuite.py b/validation/vv/testsuite.py index a4f0ed79..397831f8 100644 --- a/validation/vv/testsuite.py +++ b/validation/vv/testsuite.py @@ -2,8 +2,11 @@ from abc import ABC, abstractmethod import fourdst.atomic import scipy.integrate +import gridfire from fourdst.composition import Composition -from gridfire.engine import DynamicEngine, GraphEngine +from gridfire.engine import DynamicEngine, GraphEngine, AdaptiveEngineView, MultiscalePartitioningEngineView +from gridfire.engine import EngineTypes +from gridfire.policy import MainSequencePolicy from gridfire.type import NetIn, NetOut from gridfire.exceptions import GridFireError from gridfire.solver import CVODESolverStrategy @@ -21,6 +24,12 @@ import numpy as np import json import time +EngineTypeLookup : Dict[EngineTypes, Any] = { + EngineTypes.ADAPTIVE_ENGINE_VIEW: AdaptiveEngineView, + EngineTypes.MULTISCALE_PARTITIONING_ENGINE_VIEW: MultiscalePartitioningEngineView, + EngineTypes.GRAPH_ENGINE: GraphEngine +} + def load_network_module(filepath): module_name = os.path.basename(filepath).replace(".py", "") if module_name in sys.modules: # clear any existing module with the same name @@ -103,12 +112,16 @@ class TestSuite(ABC): self.composition : Composition = composition self.notes : str = notes - def evolve_pynucastro(self, engine: GraphEngine): + def evolve_pynucastro(self, engine: DynamicEngine): print("Evolution complete. Now building equivalent pynucastro network...") # Build equivalent pynucastro network for comparison reaclib_library : pyna.ReacLibLibrary = pyna.ReacLibLibrary() rate_names = [r.id().replace("e+","").replace("e-","").replace(", ", ",") for r in engine.getNetworkReactions()] + with open(f"{self.name}_rate_names_pynuc.txt", "w") as f: + for r_name in rate_names: + f.write(f"{r_name}\n") + goodRates : List[pyna.rates.reaclib_rate.ReacLibRate] = [] missingRates = [] @@ -156,7 +169,23 @@ class TestSuite(ABC): atol=1e-8 ) endTime = time.time() + initial_duration = endTime - startTime print("Pynucastro integration complete. Writing results to JSON...") + print("Running pynucastro a second time to account for any JIT compilation overhead...") + startTime = time.time() + sol = scipy.integrate.solve_ivp( + net.rhs, + [0, self.tMax], + Y0, + args=(self.density, self.temperature), + method="BDF", + jac=net.jacobian, + rtol=1e-5, + atol=1e-8 + ) + endTime = time.time() + final_duration = endTime - startTime + print(f"Pynucastro second integration complete. Initial run time: {initial_duration: .4f} s, Second run time: {final_duration: .4f} s") data: List[Dict[str, Union[float, Dict[str, float]]]] = [] @@ -182,7 +211,8 @@ class TestSuite(ABC): "Temperature": self.temperature, "Density": self.density, "tMax": self.tMax, - "ElapsedTime": endTime - startTime, + "RunTime0": initial_duration, + "RunTime1": final_duration, "DateCreated": datetime.now().isoformat() }, "Steps": data @@ -191,7 +221,7 @@ class TestSuite(ABC): with open(f"GridFireValidationSuite_{self.name}_pynucastro.json", "w") as f: json.dump(pynucastro_json, f, indent=4) - def evolve(self, engine: GraphEngine, netIn: NetIn, pynucastro_compare: bool = True): + def evolve(self, engine: DynamicEngine, netIn: NetIn, pynucastro_compare: bool = True, engine_type: EngineTypes | None = None): solver : CVODESolverStrategy = CVODESolverStrategy(engine) stepLogger : StepLogger = StepLogger() @@ -232,10 +262,23 @@ class TestSuite(ABC): ) if pynucastro_compare: - self.evolve_pynucastro(engine) - + if engine_type is not None: + if engine_type == EngineTypes.ADAPTIVE_ENGINE_VIEW: + print("Pynucastro comparison using AdaptiveEngineView...") + self.evolve_pynucastro(engine) + elif engine_type == EngineTypes.MULTISCALE_PARTITIONING_ENGINE_VIEW: + print("Pynucastro comparison using MultiscalePartitioningEngineView...") + graphEngine : GraphEngine = GraphEngine(self.composition, depth=3) + multiScaleEngine : MultiscalePartitioningEngineView = MultiscalePartitioningEngineView(graphEngine) + self.evolve_pynucastro(multiScaleEngine) + elif engine_type == EngineTypes.GRAPH_ENGINE: + print("Pynucastro comparison using GraphEngine...") + graphEngine : GraphEngine = GraphEngine(self.composition, depth=3) + self.evolve_pynucastro(graphEngine) + else: + print(f"Pynucastro comparison not implemented for engine type: {engine_type}") @abstractmethod - def __call__(self): + def __call__(self, pynucastro_compare: bool = False, pync_engine: str = "AdaptiveEngineView"): pass