diff --git a/src/include/gridfire/engine/engine_abstract.h b/src/include/gridfire/engine/engine_abstract.h index d74d2368..3af6aad9 100644 --- a/src/include/gridfire/engine/engine_abstract.h +++ b/src/include/gridfire/engine/engine_abstract.h @@ -157,14 +157,19 @@ namespace gridfire { double rho ) const = 0; + virtual void generateJacobianMatrix( + const fourdst::composition::Composition& comp, + double T9, + double rho, + const std::vector& activeSpecies + ) const = 0; + virtual void generateJacobianMatrix( const fourdst::composition::Composition& comp, double T9, double rho, const SparsityPattern& sparsityPattern - ) const { - throw std::logic_error("Sparsity pattern not supported by this engine."); - } + ) const = 0; /** * @brief Get an entry from the previously generated Jacobian matrix. diff --git a/src/include/gridfire/engine/engine_graph.h b/src/include/gridfire/engine/engine_graph.h index ba2347c8..7b0b8480 100644 --- a/src/include/gridfire/engine/engine_graph.h +++ b/src/include/gridfire/engine/engine_graph.h @@ -12,7 +12,6 @@ #include "gridfire/screening/screening_types.h" #include "gridfire/partition/partition_abstract.h" #include "gridfire/engine/procedures/construction.h" -#include "gridfire/utils/general_composition.h" #include #include @@ -235,6 +234,25 @@ namespace gridfire { double rho ) const override; + /** + * @brief Generates the Jacobian matrix for the current state with a specified set of active species. + * generally this will be much faster than the full matrix generation. Here we use forward mode + * to generate the Jacobian only for the active species. + * @param comp The Composition object containing current abundances. + * @param T9 The temperature in units of 10^9 K. + * @param rho The density in g/cm^3. + * @param activeSpecies A vector of Species objects representing the active species. + * + * @see getJacobianMatrixEntry() + * @see generateJacobianMatrix() + */ + void generateJacobianMatrix( + const fourdst::composition::Composition& comp, + double T9, + double rho, + const std::vector& activeSpecies + ) const override; + /** * @brief Generates the Jacobian matrix for the current state with a specified sparsity pattern. * @@ -717,6 +735,7 @@ namespace gridfire { // Forward cacheing size_t reaction_index; reaction::ReactionType reaction_type; + uint64_t reaction_hash; std::vector unique_reactant_indices; std::vector reactant_powers; double symmetry_factor; @@ -800,6 +819,7 @@ namespace gridfire { mutable CppAD::ADFun m_epsADFun; ///< CppAD function for the energy generation rate. mutable CppAD::sparse_jac_work m_jac_work; ///< Work object for sparse Jacobian calculations. CppAD::sparse_rc> m_full_jacobian_sparsity_pattern; ///< Full sparsity pattern for the Jacobian matrix. + std::set> m_full_sparsity_set; ///< For quick lookups of the base sparsity pattern std::vector> m_atomicReverseRates; @@ -813,6 +833,7 @@ namespace gridfire { BuildDepthType m_depth; std::vector m_precomputedReactions; ///< Precomputed reactions for efficiency. + std::unordered_map m_precomputedReactionIndexMap; ///< Set of hashed precomputed reactions for quick lookup. std::unique_ptr m_partitionFunction; ///< Partition function for the network. private: @@ -890,10 +911,10 @@ namespace gridfire { [[nodiscard]] StepDerivatives calculateAllDerivativesUsingPrecomputation( const fourdst::composition::Composition &comp, - const std::vector& bare_rates, + const std::vector &bare_rates, const std::vector &bare_reverse_rates, double T9, - double rho + double rho, const reaction::ReactionSet &activeReactions ) const; /** diff --git a/src/include/gridfire/engine/views/engine_adaptive.h b/src/include/gridfire/engine/views/engine_adaptive.h index a65e9808..b1645d69 100644 --- a/src/include/gridfire/engine/views/engine_adaptive.h +++ b/src/include/gridfire/engine/views/engine_adaptive.h @@ -140,6 +140,20 @@ namespace gridfire { double rho ) const override; + void generateJacobianMatrix( + const fourdst::composition::Composition &comp, + double T9, + double rho, + const std::vector &activeSpecies + ) const override; + + void generateJacobianMatrix( + const fourdst::composition::Composition &comp, + double T9, + double rho, + const SparsityPattern &sparsityPattern + ) const override; + /** * @brief Gets an entry from the Jacobian matrix for the active species. * @@ -409,7 +423,7 @@ namespace gridfire { * 5. For each reaction, it calls the base engine's `calculateMolarReactionFlow` to get the flow rate. * 6. Stores the reaction pointer and its flow rate in a `ReactionFlow` struct and adds it to the returned vector. */ - std::pair, fourdst::composition::Composition> calculateAllReactionFlows( + [[nodiscard]] std::pair, fourdst::composition::Composition> calculateAllReactionFlows( const NetIn& netIn ) const; /** diff --git a/src/include/gridfire/engine/views/engine_defined.h b/src/include/gridfire/engine/views/engine_defined.h index b743edaf..68210683 100644 --- a/src/include/gridfire/engine/views/engine_defined.h +++ b/src/include/gridfire/engine/views/engine_defined.h @@ -65,9 +65,44 @@ namespace gridfire{ */ void generateJacobianMatrix( const fourdst::composition::Composition& comp, - const double T9, - const double rho + double T9, + double rho ) const override; + + /** + * @brief Generates the Jacobian matrix for the active species. + * + * @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 activeSpecies The vector of active species to include in the Jacobian. + * + * @throws std::runtime_error If the view is stale. + */ + void generateJacobianMatrix( + const fourdst::composition::Composition &comp, + double T9, + double rho, + const std::vector &activeSpecies + ) const override; + + /** + * @brief Generates the Jacobian matrix for a given sparsity pattern + * + * @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 sparsityPattern The sparsity pattern to use for the Jacobian matrix. + * + * @throws std::runtime_error If the view is stale. + */ + void generateJacobianMatrix( + const fourdst::composition::Composition &comp, + double T9, + double rho, + const SparsityPattern &sparsityPattern + ) const override; + /** * @brief Gets an entry from the Jacobian matrix for the active species. * diff --git a/src/include/gridfire/engine/views/engine_multiscale.h b/src/include/gridfire/engine/views/engine_multiscale.h index 33b4842d..4ffb608c 100644 --- a/src/include/gridfire/engine/views/engine_multiscale.h +++ b/src/include/gridfire/engine/views/engine_multiscale.h @@ -271,6 +271,64 @@ namespace gridfire { double rho ) const override; + /** + * @brief Generates the Jacobian matrix for a subset of active species. + * + * @param comp The current composition. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @param activeSpecies The subset of species to include in the Jacobian. + * + * @par Purpose + * To compute a reduced Jacobian matrix for implicit solvers that only + * consider a subset of species. + * + * @par How + * Similar to the full Jacobian generation, it first checks the QSE cache. + * On a hit, it calls the base engine's `generateJacobianMatrix` with + * the specified active species. The returned Jacobian still reflects + * the full network, but only for the active species subset. + * + * @pre The engine must have a valid QSE cache entry for the given state. + * @post The base engine's internal Jacobian is updated for the active species. + * + * @throws exceptions::StaleEngineError If the QSE cache misses. + */ + void generateJacobianMatrix( + const fourdst::composition::Composition &comp, + double T9, + double rho, + const std::vector &activeSpecies + ) const override; + + /** + * @brief Generates the Jacobian matrix using a sparsity pattern. + * + * @param comp The current composition. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @param sparsityPattern The sparsity pattern to use for the Jacobian. + * + * @par Purpose + * To compute the Jacobian matrix while leveraging a known sparsity pattern + * for efficiency. This is effectively a lower level version of the active species method. + * + * @par How + * It first checks the QSE cache. On a hit, it delegates to the base engine's + * `generateJacobianMatrix` method with the provided sparsity pattern. + * + * @pre The engine must have a valid QSE cache entry for the given state. + * @post The base engine's internal Jacobian is updated according to the sparsity pattern. + * + * @throws exceptions::StaleEngineError If the QSE cache misses. + */ + void generateJacobianMatrix( + const fourdst::composition::Composition &comp, + double T9, + double rho, + const SparsityPattern &sparsityPattern + ) const override; + /** * @brief Gets an entry from the previously generated Jacobian matrix. * @@ -720,6 +778,12 @@ namespace gridfire { const NetIn &netIn ); + bool involvesSpecies(const fourdst::atomic::Species &species) const; + + bool involvesSpeciesInQSE(const fourdst::atomic::Species &species) const; + + bool involvesSpeciesInDynamic(const fourdst::atomic::Species &species) const; + private: /** @@ -806,7 +870,7 @@ namespace gridfire { */ MultiscalePartitioningEngineView* m_view; /** - * @brief Indices of the species to solve for in the QSE group. + * @brief The set of species to solve for in the QSE group. */ const std::set& m_qse_solve_species; /** diff --git a/src/include/gridfire/engine/views/engine_priming.h b/src/include/gridfire/engine/views/engine_priming.h index d5604e2b..fed5131b 100644 --- a/src/include/gridfire/engine/views/engine_priming.h +++ b/src/include/gridfire/engine/views/engine_priming.h @@ -1,9 +1,7 @@ #pragma once -#include "gridfire/engine/engine_abstract.h" #include "gridfire/engine/views/engine_defined.h" - #include "fourdst/logging/logging.h" #include "fourdst/composition/atomicSpecies.h" @@ -51,7 +49,7 @@ namespace gridfire { private: - quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); fourdst::atomic::Species m_primingSpecies; ///< The priming species, if specified. private: /** diff --git a/src/include/gridfire/exceptions/error_utils.h b/src/include/gridfire/exceptions/error_utils.h new file mode 100644 index 00000000..4af77b9d --- /dev/null +++ b/src/include/gridfire/exceptions/error_utils.h @@ -0,0 +1,23 @@ +#pragma once + +#include +#include +#include + +namespace gridfire::exceptions { + class UtilityError : public std::exception {}; + + class HashingError final : public UtilityError { + public: + explicit HashingError() = default; + + explicit HashingError(std::string message) + : m_message(std::move(message)) {} + + [[nodiscard]] const char* what() const noexcept override { + return m_message.c_str(); + } + private: + std::string m_message; + }; +} diff --git a/src/include/gridfire/exceptions/exceptions.h b/src/include/gridfire/exceptions/exceptions.h index 8cbedcd2..8aad7c52 100644 --- a/src/include/gridfire/exceptions/exceptions.h +++ b/src/include/gridfire/exceptions/exceptions.h @@ -1,3 +1,4 @@ #pragma once -#include "gridfire/exceptions/error_engine.h" \ No newline at end of file +#include "gridfire/exceptions/error_engine.h" +#include "gridfire/exceptions/error_utils.h" \ No newline at end of file diff --git a/src/include/gridfire/utils/hashing.h b/src/include/gridfire/utils/hashing.h index dd4100b2..66eee8f3 100644 --- a/src/include/gridfire/utils/hashing.h +++ b/src/include/gridfire/utils/hashing.h @@ -2,6 +2,9 @@ #include +#include "gridfire/exceptions/exceptions.h" +#include "gridfire/reaction/reaction.h" + namespace gridfire::utils { /** * @brief Generate a unique hash for an isotope given its mass number (A) and atomic number (Z). @@ -16,4 +19,47 @@ namespace gridfire::utils { return (static_cast(a) << 8) | static_cast(z); } -} \ No newline at end of file + namespace hashing::reaction { + static std::uint64_t splitmix64(std::uint64_t x) noexcept { + x += 0x9E3779B97F4A7C15ULL; + x = (x ^ (x >> 30)) * 0xBF58476D1CE4E5B9ULL; + x = (x ^ (x >> 27)) * 0x94D049BB133111EBULL; + x ^= (x >> 31); + return x; + } + + static std::uint64_t mix_species(const unsigned a, const unsigned z) noexcept { + const std::uint64_t code = (static_cast(a) << 7) | static_cast(z); + return splitmix64(code); + } + + static std::uint64_t multiset_combine(std::uint64_t acc, const std::uint64_t x) noexcept { + acc += x; + acc ^= (x << 23) | (x >> (64 - 23)); + acc = splitmix64(acc); + return acc; + } + } + + inline std::uint64_t hash_reaction(const reaction::Reaction& reaction) noexcept { + using namespace hashing::reaction; + + std::uint64_t hR = 0; + + for (const auto& s : reaction.reactants()) { + hR = multiset_combine(hR, mix_species(static_cast(s.a()), + static_cast(s.z()))); + } + + std::uint64_t hP = 0; + for (const auto& s : reaction.products()) { + hP = multiset_combine(hP, mix_species(static_cast(s.a()), + static_cast(s.z()))); + } + + std::uint64_t h = splitmix64(hR ^ 0xC3A5C85C97CB3127ULL); + h ^= splitmix64((hP << 1) | (hP >> 63)); + return splitmix64(h); + } + +} diff --git a/src/lib/engine/engine_graph.cpp b/src/lib/engine/engine_graph.cpp index 31f84dc7..cabf2fa5 100644 --- a/src/lib/engine/engine_graph.cpp +++ b/src/lib/engine/engine_graph.cpp @@ -5,6 +5,7 @@ #include "gridfire/engine/procedures/priming.h" #include "gridfire/partition/partition_ground.h" #include "gridfire/engine/procedures/construction.h" +#include "gridfire/utils/hashing.h" #include "fourdst/composition/species.h" #include "fourdst/composition/atomicSpecies.h" @@ -75,12 +76,11 @@ namespace gridfire { if (m_usePrecomputation) { std::vector bare_rates; std::vector bare_reverse_rates; - bare_rates.reserve(m_reactions.size()); - bare_reverse_rates.reserve(m_reactions.size()); + bare_rates.reserve(activeReactions.size()); + bare_reverse_rates.reserve(activeReactions.size()); - // TODO: Add cache to this - - for (const auto& reaction: m_reactions) { + for (const auto& reaction: activeReactions) { + assert(m_reactions.contains(*reaction)); // A bug which results in this failing indicates a serious internal inconsistency and should only be present during development. bare_rates.push_back(reaction->calculate_rate(T9, rho, Ye, mue, comp.getMolarAbundanceVector(), m_indexToSpeciesMap)); if (reaction->type() != reaction::ReactionType::WEAK) { bare_reverse_rates.push_back(calculateReverseRate(*reaction, T9, rho, comp)); @@ -88,7 +88,7 @@ namespace gridfire { } // --- The public facing interface can always use the precomputed version since taping is done internally --- - return calculateAllDerivativesUsingPrecomputation(comp, bare_rates, bare_reverse_rates, T9, rho); + return calculateAllDerivativesUsingPrecomputation(comp, bare_rates, bare_reverse_rates, T9, rho, activeReactions); } else { return calculateAllDerivatives( comp.getMolarAbundanceVector(), @@ -171,6 +171,8 @@ namespace gridfire { collectAtomicReverseRateAtomicBases(); generateStoichiometryMatrix(); reserveJacobianMatrix(); + + // PERF: These do *a lot* of the same work* can this be optimized to only do the common work once? recordADTape(); // Record the AD tape for the RHS function recordEpsADTape(); // Record the AD tape for the energy generation rate function @@ -183,6 +185,17 @@ namespace gridfire { m_rhsADFun.subgraph_sparsity(select_domain, select_range, false, m_full_jacobian_sparsity_pattern); m_jac_work.clear(); + m_full_sparsity_set.clear(); + const auto& rows = m_full_jacobian_sparsity_pattern.row(); + const auto& cols = m_full_jacobian_sparsity_pattern.col(); + const size_t nnz = m_full_jacobian_sparsity_pattern.nnz(); + + for (size_t k = 0; k < nnz; ++k) { + if (cols[k] < m_networkSpecies.size()) { + m_full_sparsity_set.insert(std::make_pair(rows[k], cols[k])); + } + } + precomputeNetwork(); LOG_INFO(m_logger, "Internal maps synchronized. Network contains {} species and {} reactions.", m_networkSpecies.size(), m_reactions.size()); @@ -190,7 +203,6 @@ namespace gridfire { // --- Network Graph Construction Methods --- void GraphEngine::collectNetworkSpecies() { - //TODO: Ensure consistent ordering in the m_networkSpecies vector so that it is ordered by species mass. m_networkSpecies.clear(); m_networkSpeciesMap.clear(); @@ -217,7 +229,7 @@ namespace gridfire { } } // TODO: Currently this works. We sort the vector based on mass so that for the same set of species we always get the same ordering and we get the same ordering as a composition with the same set of species - // However, we need some checks so that when we get a composition we confirm that it is the same ordering / contains teh same species. This is important for the ODE integrator to work properly. + // However, we need some checks so that when we get a composition we confirm that it is the same ordering / contains the same species. This is important for the ODE integrator to work properly. std::ranges::sort(m_networkSpecies, [](const fourdst::atomic::Species& a, const fourdst::atomic::Species& b) -> bool { return a.mass() < b.mass(); // Otherwise, sort by mass }); @@ -255,14 +267,10 @@ namespace gridfire { // --- Basic Accessors and Queries --- const std::vector& GraphEngine::getNetworkSpecies() const { - // Returns a constant reference to the vector of unique species in the network. - // LOG_TRACE_L3(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size()); return m_networkSpecies; } const reaction::ReactionSet& GraphEngine::getNetworkReactions() const { - // Returns a constant reference to the set of reactions in the network. - // LOG_TRACE_L3(m_logger, "Providing access to network reactions set. Size: {}.", m_reactions.size()); return m_reactions; } @@ -573,64 +581,75 @@ namespace gridfire { const std::vector &bare_rates, const std::vector &bare_reverse_rates, const double T9, - const double rho + const double rho, + const reaction::ReactionSet &activeReactions ) const { // --- Calculate screening factors --- const std::vector screeningFactors = m_screeningModel->calculateScreeningFactors( - m_reactions, + activeReactions, m_networkSpecies, comp.getMolarAbundanceVector(), T9, rho ); - // TODO: Fix up the precomputation to use the new comp in interface as opposed to a raw vector of molar abundances - // This will require carefully checking the way the precomputation is stashed. - // --- Optimized loop --- std::vector molarReactionFlows; molarReactionFlows.reserve(m_precomputedReactions.size()); - for (const auto& precomp : m_precomputedReactions) { + size_t reactionCounter = 0; + for (const auto& reaction : activeReactions) { + // --- Efficient lookup of only the active reactions --- + uint64_t reactionHash = utils::hash_reaction(*reaction); + const size_t reactionIndex = m_precomputedReactionIndexMap.at(reactionHash); + PrecomputedReaction precomputedReaction = m_precomputedReactions[reactionIndex]; + + // --- Forward abundance product --- double forwardAbundanceProduct = 1.0; - for (size_t i = 0; i < precomp.unique_reactant_indices.size(); ++i) { - const size_t reactantIndex = precomp.unique_reactant_indices[i]; + for (size_t i = 0; i < precomputedReaction.unique_reactant_indices.size(); ++i) { + const size_t reactantIndex = precomputedReaction.unique_reactant_indices[i]; const fourdst::atomic::Species& reactant = m_networkSpecies[reactantIndex]; - const int power = precomp.reactant_powers[i]; + const int power = precomputedReaction.reactant_powers[i]; forwardAbundanceProduct *= std::pow(comp.getMolarAbundance(reactant), power); } - const double bare_rate = bare_rates[precomp.reaction_index]; - const double screeningFactor = screeningFactors[precomp.reaction_index]; - const size_t numReactants = m_reactions[precomp.reaction_index].reactants().size(); - const size_t numProducts = m_reactions[precomp.reaction_index].products().size(); + const double bare_rate = bare_rates.at(reactionCounter); + const double screeningFactor = screeningFactors[reactionCounter]; + const size_t numReactants = m_reactions[reactionIndex].reactants().size(); + const size_t numProducts = m_reactions[reactionIndex].products().size(); + + // --- Forward reaction flow --- const double forwardMolarReactionFlow = screeningFactor * bare_rate * - precomp.symmetry_factor * + precomputedReaction.symmetry_factor * forwardAbundanceProduct * std::pow(rho, numReactants > 1 ? static_cast(numReactants) - 1 : 0.0); + // --- Reverse reaction flow --- + // Only do this is the reaction has a non-zero reverse symmetry factor (i.e. is reversible) double reverseMolarReactionFlow = 0.0; - if (precomp.reverse_symmetry_factor != 0.0 and m_useReverseReactions) { - const double bare_reverse_rate = bare_reverse_rates[precomp.reaction_index]; + if (precomputedReaction.reverse_symmetry_factor != 0.0 and m_useReverseReactions) { + const double bare_reverse_rate = bare_reverse_rates.at(reactionCounter); + double reverseAbundanceProduct = 1.0; - for (size_t i = 0; i < precomp.unique_product_indices.size(); ++i) { - const size_t productIndex = precomp.unique_product_indices[i]; + for (size_t i = 0; i < precomputedReaction.unique_product_indices.size(); ++i) { + const size_t productIndex = precomputedReaction.unique_product_indices[i]; const fourdst::atomic::Species& product = m_networkSpecies[productIndex]; - reverseAbundanceProduct *= std::pow(comp.getMolarAbundance(product), precomp.product_powers[i]); + reverseAbundanceProduct *= std::pow(comp.getMolarAbundance(product), precomputedReaction.product_powers[i]); } + reverseMolarReactionFlow = screeningFactor * bare_reverse_rate * - precomp.reverse_symmetry_factor * + precomputedReaction.reverse_symmetry_factor * reverseAbundanceProduct * std::pow(rho, numProducts > 1 ? static_cast(numProducts) - 1 : 0.0); } molarReactionFlows.push_back(forwardMolarReactionFlow - reverseMolarReactionFlow); - + reactionCounter++; } // --- Assemble molar abundance derivatives --- @@ -638,9 +657,12 @@ namespace gridfire { for (const auto& species: m_networkSpecies) { result.dydt[species] = 0.0; // Initialize the change in abundance for each network species to 0 } - for (size_t j = 0; j < m_precomputedReactions.size(); ++j) { + + reactionCounter = 0; + for (const auto& reaction: activeReactions) { + size_t j = m_precomputedReactionIndexMap.at(utils::hash_reaction(*reaction)); const auto& precomp = m_precomputedReactions[j]; - const double R_j = molarReactionFlows[j]; + const double R_j = molarReactionFlows[reactionCounter]; for (size_t i = 0; i < precomp.affected_species_indices.size(); ++i) { const size_t speciesIndex = precomp.affected_species_indices[i]; @@ -651,6 +673,7 @@ namespace gridfire { // Update the derivative for this species result.dydt.at(species) += static_cast(stoichiometricCoefficient) * R_j; } + reactionCounter++; } // --- Calculate the nuclear energy generation rate --- @@ -802,12 +825,51 @@ namespace gridfire { LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); } + void GraphEngine::generateJacobianMatrix( + const fourdst::composition::Composition &comp, + const double T9, + const double rho, + const std::vector &activeSpecies + ) const { + // PERF: For small k it may make sense to implement a purley forward mode AD computation, some heuristic could be used to switch between the two methods based on k and total network species + const size_t k_active = activeSpecies.size(); + + // --- 1. Get the list of global indices --- + std::vector active_indices; + active_indices.reserve(k_active); + + for (const auto& species : activeSpecies) { + assert(involvesSpecies(species)); + active_indices.push_back(getSpeciesIndex(species)); + } + + // --- 2. Build the k x k sparsity pattern --- + SparsityPattern sparsityPattern; + sparsityPattern.reserve(k_active * k_active); + + for (const size_t i_global : active_indices) { // k rows + for (const size_t j_global : active_indices) { // k columns + sparsityPattern.emplace_back(i_global, j_global); + } + } + + // --- 3. Call the sparse reverse-mode implementation --- + generateJacobianMatrix(comp, T9, rho, sparsityPattern); + } + void GraphEngine::generateJacobianMatrix( const fourdst::composition::Composition &comp, const double T9, const double rho, const SparsityPattern &sparsityPattern ) const { + SparsityPattern intersectionSparsityPattern; + for (const auto& entry : sparsityPattern) { + if (m_full_sparsity_set.contains(entry)) { + intersectionSparsityPattern.push_back(entry); + } + } + // --- Pack the input variables into a vector for CppAD --- const size_t numSpecies = m_networkSpecies.size(); std::vector x(numSpecies + 2, 0.0); @@ -819,13 +881,13 @@ namespace gridfire { x[numSpecies + 1] = rho; // --- Convert into CppAD Sparsity pattern --- - const size_t nnz = sparsityPattern.size(); // Number of non-zero entries in the sparsity pattern + const size_t nnz = intersectionSparsityPattern.size(); // Number of non-zero entries in the sparsity pattern std::vector row_indices(nnz); std::vector col_indices(nnz); for (size_t k = 0; k < nnz; ++k) { - row_indices[k] = sparsityPattern[k].first; - col_indices[k] = sparsityPattern[k].second; + row_indices[k] = intersectionSparsityPattern[k].first; + col_indices[k] = intersectionSparsityPattern[k].second; } std::vector values(nnz); @@ -834,7 +896,7 @@ namespace gridfire { CppAD::sparse_rc> CppAD_sparsity_pattern(num_rows_jac, num_cols_jac, nnz); for (size_t k = 0; k < nnz; ++k) { - CppAD_sparsity_pattern.set(k, sparsityPattern[k].first, sparsityPattern[k].second); + CppAD_sparsity_pattern.set(k, intersectionSparsityPattern[k].first, intersectionSparsityPattern[k].second); } CppAD::sparse_rcv, std::vector> jac_subset(CppAD_sparsity_pattern); @@ -854,7 +916,7 @@ namespace gridfire { const size_t col = jac_subset.col()[k]; const double value = jac_subset.val()[k]; - if (std::abs(value) > MIN_JACOBIAN_THRESHOLD) { + if (std::abs(value) > MIN_JACOBIAN_THRESHOLD || row == col) { // Always keep diagonal elements to avoid pathological stiffness m_jacobianMatrix(row, col) = value; // Insert into the sparse matrix } } @@ -1241,12 +1303,18 @@ namespace gridfire { m_precomputedReactions.clear(); m_precomputedReactions.reserve(m_reactions.size()); + m_precomputedReactionIndexMap.clear(); + m_precomputedReactionIndexMap.reserve(m_reactions.size()); for (size_t i = 0; i < m_reactions.size(); ++i) { const auto& reaction = m_reactions[i]; PrecomputedReaction precomp; precomp.reaction_index = i; precomp.reaction_type = reaction.type(); + uint64_t reactionHash = utils::hash_reaction(reaction); + + precomp.reaction_hash = reactionHash; + m_precomputedReactionIndexMap[reactionHash] = i; // --- Precompute forward reaction information --- // Count occurrences for each reactant to determine powers and symmetry @@ -1298,6 +1366,7 @@ namespace gridfire { m_precomputedReactions.push_back(std::move(precomp)); } + LOG_TRACE_L1(m_logger, "Pre-computation complete. Precomputed data for {} reactions.", m_precomputedReactions.size()); } bool GraphEngine::AtomicReverseRate::forward( diff --git a/src/lib/engine/procedures/priming.cpp b/src/lib/engine/procedures/priming.cpp index 83a5ce0e..93d37655 100644 --- a/src/lib/engine/procedures/priming.cpp +++ b/src/lib/engine/procedures/priming.cpp @@ -231,11 +231,11 @@ namespace gridfire { // Store the request instead of applying it immediately. requests.push_back({primingSpecies, equilibriumMassFraction, dominantChannel->reactants()}); } else { - LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name()); + LOG_TRACE_L1(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name()); report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL; } } else { - LOG_WARNING(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name()); + LOG_TRACE_L2(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name()); // For species with no destruction, we can't calculate an equilibrium. // We add a request with a tiny fallback abundance to ensure it exists in the network. requests.push_back({primingSpecies, 1e-40, {}}); diff --git a/src/lib/engine/views/engine_adaptive.cpp b/src/lib/engine/views/engine_adaptive.cpp index d66c8a3c..5987e2cd 100644 --- a/src/lib/engine/views/engine_adaptive.cpp +++ b/src/lib/engine/views/engine_adaptive.cpp @@ -176,9 +176,29 @@ namespace gridfire { const fourdst::composition::Composition &comp, const double T9, const double rho + ) const { + generateJacobianMatrix(comp, T9, rho, m_activeSpecies); + } + + void AdaptiveEngineView::generateJacobianMatrix( + const fourdst::composition::Composition &comp, + const double T9, + const double rho, + const std::vector &activeSpecies ) const { validateState(); - m_baseEngine.generateJacobianMatrix(comp, T9, rho); + m_baseEngine.generateJacobianMatrix(comp, T9, rho, activeSpecies); + + } + + void AdaptiveEngineView::generateJacobianMatrix( + const fourdst::composition::Composition &comp, + const double T9, + const double rho, + const SparsityPattern &sparsityPattern + ) const { + validateState(); + m_baseEngine.generateJacobianMatrix(comp, T9, rho, sparsityPattern); } double AdaptiveEngineView::getJacobianMatrixEntry( diff --git a/src/lib/engine/views/engine_defined.cpp b/src/lib/engine/views/engine_defined.cpp index b371a371..51369f80 100644 --- a/src/lib/engine/views/engine_defined.cpp +++ b/src/lib/engine/views/engine_defined.cpp @@ -86,11 +86,39 @@ namespace gridfire { const double rho ) const { validateNetworkState(); - + if (!m_activeSpeciesVectorCache.has_value()) { + m_activeSpeciesVectorCache = std::vector(m_activeSpecies.begin(), m_activeSpecies.end()); + } const MaskedComposition masked(comp, m_activeSpecies); + m_baseEngine.generateJacobianMatrix(masked, T9, rho, m_activeSpeciesVectorCache.value()); + } - // TODO: We likely want to be able to think more carefully about this so that the jacobian matches the active species/reactions - m_baseEngine.generateJacobianMatrix(masked, T9, rho); + void DefinedEngineView::generateJacobianMatrix( + const fourdst::composition::Composition &comp, + const double T9, + const double rho, + const std::vector &activeSpecies + ) const { + validateNetworkState(); + + const std::set activeSpeciesSet( + activeSpecies.begin(), + activeSpecies.end() + ); + + const MaskedComposition masked(comp, activeSpeciesSet); + m_baseEngine.generateJacobianMatrix(masked, T9, rho, activeSpecies); + } + + void DefinedEngineView::generateJacobianMatrix( + const fourdst::composition::Composition &comp, + const double T9, + const double rho, + const SparsityPattern &sparsityPattern + ) const { + validateNetworkState(); + const MaskedComposition masked(comp, m_activeSpecies); + m_baseEngine.generateJacobianMatrix(masked, T9, rho, sparsityPattern); } double DefinedEngineView::getJacobianMatrixEntry( diff --git a/src/lib/engine/views/engine_multiscale.cpp b/src/lib/engine/views/engine_multiscale.cpp index 3a8abebd..b3bb0ecd 100644 --- a/src/lib/engine/views/engine_multiscale.cpp +++ b/src/lib/engine/views/engine_multiscale.cpp @@ -17,8 +17,6 @@ #include "quill/LogMacros.h" #include "quill/Logger.h" -static std::ofstream debug_multiscale_log("debug_multiscale_log.txt"); - namespace { using namespace fourdst::atomic; //TODO: Replace all calls to this function with composition.getMolarAbundanceVector() so that @@ -214,8 +212,41 @@ namespace gridfire { const double T9, const double rho ) const { - // TODO: Add sparsity pattern to this to prevent base engine from doing unnecessary work. - m_baseEngine.generateJacobianMatrix(comp, T9, rho); + // We do not need to generate the jacobian for QSE species since those entries are by definition 0 + m_baseEngine.generateJacobianMatrix(comp, T9, rho, m_dynamic_species); + } + + void MultiscalePartitioningEngineView::generateJacobianMatrix( + const fourdst::composition::Composition &comp, + const double T9, + const double rho, + const std::vector &activeSpecies + ) const { + const bool activeSpeciesIsSubset = std::ranges::any_of(activeSpecies, [&](const auto& species) -> bool { + return !involvesSpecies(species); + }); + if (activeSpeciesIsSubset) { + LOG_CRITICAL(m_logger, "Active species contains species not in the network partition. Cannot generate Jacobian matrix for active species."); + throw std::runtime_error("Active species contains species not in the network partition. Cannot generate Jacobian matrix for active species."); + } + + std::vector dynamicActiveSpeciesIntersection; + for (const auto& species : activeSpecies) { + if (involvesSpeciesInDynamic(species)) { + dynamicActiveSpeciesIntersection.push_back(species); + } + } + + m_baseEngine.generateJacobianMatrix(comp, T9, rho, dynamicActiveSpeciesIntersection); + } + + void MultiscalePartitioningEngineView::generateJacobianMatrix( + const fourdst::composition::Composition &comp, + const double T9, + const double rho, + const SparsityPattern &sparsityPattern + ) const { + return m_baseEngine.generateJacobianMatrix(comp, T9, rho, sparsityPattern); } double MultiscalePartitioningEngineView::getJacobianMatrixEntry( @@ -811,6 +842,26 @@ namespace gridfire { return equilibrateNetwork(primingReport.primedComposition, T9, rho); } + bool MultiscalePartitioningEngineView::involvesSpecies( + const Species &species + ) const { + if (involvesSpeciesInQSE(species)) return true; // check this first since the vector is likely to be smaller so short circuit cost is less on average + if (involvesSpeciesInDynamic(species)) return true; + return false; + } + + bool MultiscalePartitioningEngineView::involvesSpeciesInQSE( + const Species &species + ) const { + return std::ranges::find(m_algebraic_species, species) != m_algebraic_species.end(); + } + + bool MultiscalePartitioningEngineView::involvesSpeciesInDynamic( + const Species &species + ) const { + return std::ranges::find(m_dynamic_species, species) != m_dynamic_species.end(); + } + size_t MultiscalePartitioningEngineView::getSpeciesIndex(const Species &species) const { return m_baseEngine.getSpeciesIndex(species); } @@ -1083,7 +1134,7 @@ namespace gridfire { } - if (bool group_is_coupled = (coupling_flux / leakage_flux) > FLUX_RATIO_THRESHOLD) { + if (coupling_flux / leakage_flux > FLUX_RATIO_THRESHOLD) { LOG_TRACE_L1( m_logger, "Group containing {} is in equilibrium due to high coupling flux and balanced creation and destruction rate: , ",