diff --git a/src/include/gridfire/engine/types/reporting.h b/src/include/gridfire/engine/types/reporting.h index ace5b125..a4a4a8cd 100644 --- a/src/include/gridfire/engine/types/reporting.h +++ b/src/include/gridfire/engine/types/reporting.h @@ -27,13 +27,8 @@ namespace gridfire { * @see PrimingReport for data associated with each status. */ enum class PrimingReportStatus { - FULL_SUCCESS = 0, - NO_SPECIES_TO_PRIME = 1, - MAX_ITERATIONS_REACHED = 2, - FAILED_TO_FINALIZE_COMPOSITION = 3, - FAILED_TO_FIND_CREATION_CHANNEL = 4, - FAILED_TO_FIND_PRIMING_REACTIONS = 5, - BASE_NETWORK_TOO_SHALLOW = 6 + SUCCESS, + SOLVER_FAILURE, }; /** @@ -43,13 +38,8 @@ namespace gridfire { * The map contains entries for all PrimingReportStatus values. */ inline std::map PrimingReportStatusStrings = { - {PrimingReportStatus::FULL_SUCCESS, "Full Success"}, - {PrimingReportStatus::NO_SPECIES_TO_PRIME, "No Species to Prime"}, - {PrimingReportStatus::MAX_ITERATIONS_REACHED, "Max Iterations Reached"}, - {PrimingReportStatus::FAILED_TO_FINALIZE_COMPOSITION, "Failed to Finalize Composition"}, - {PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL, "Failed to Find Creation Channel"}, - {PrimingReportStatus::FAILED_TO_FIND_PRIMING_REACTIONS, "Failed to Find Priming Reactions"}, - {PrimingReportStatus::BASE_NETWORK_TOO_SHALLOW, "Base Network Too Shallow"} + {PrimingReportStatus::SUCCESS, "SUCCESS"}, + {PrimingReportStatus::SOLVER_FAILURE, "SOLVER_FAILURE"}, }; /** @@ -65,11 +55,6 @@ namespace gridfire { struct PrimingReport { /** Finalized composition after priming. */ fourdst::composition::Composition primedComposition; - /** - * List of pairs (species, rate change) representing destruction (<0) - * or creation (>0) rates of species during priming. - */ - std::vector> massFractionChanges; /** True if priming completed without error. */ bool success; /** Detailed status code indicating the result. */ @@ -93,4 +78,11 @@ namespace gridfire { } }; + enum class SpeciesStatus { + ACTIVE, + EQUILIBRIUM, + INACTIVE_FLOW, + NOT_PRESENT + }; + } \ No newline at end of file diff --git a/src/lib/engine/engine_graph.cpp b/src/lib/engine/engine_graph.cpp index 64ae284d..151a65ac 100644 --- a/src/lib/engine/engine_graph.cpp +++ b/src/lib/engine/engine_graph.cpp @@ -25,6 +25,7 @@ #include #include +#include #include "cppad/cppad.hpp" #include "cppad/utility/sparse_rc.hpp" @@ -180,7 +181,6 @@ namespace gridfire { populateSpeciesToIndexMap(); collectAtomicReverseRateAtomicBases(); generateStoichiometryMatrix(); - reserveJacobianMatrix(); recordADTape(); // Record the AD tape for the RHS of the ODE (dY/di and dEps/di) for all independent variables i @@ -264,18 +264,6 @@ namespace gridfire { } } - void GraphEngine::reserveJacobianMatrix() const { - // The implementation of this function (and others) constrains this nuclear network to a constant temperature and density during - // each evaluation. - const size_t numSpecies = m_networkSpecies.size(); - m_jacobianMatrix.clear(); - m_jacobianMatrix.resize(numSpecies, numSpecies, false); // Sparse matrix, no initial values - LOG_TRACE_L2(m_logger, "Jacobian matrix resized to {} rows and {} columns.", - m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); - - m_jacobianMatrixState = JacobianMatrixState::UNINITIALIZED; - } - // --- Basic Accessors and Queries --- const std::vector& GraphEngine::getNetworkSpecies() const { return m_networkSpecies; @@ -287,7 +275,6 @@ namespace gridfire { void GraphEngine::setNetworkReactions(const reaction::ReactionSet &reactions) { m_reactions = reactions; - m_jacobianMatrixState = JacobianMatrixState::STALE; syncInternalMaps(); } @@ -521,9 +508,6 @@ namespace gridfire { } void GraphEngine::setUseReverseReactions(const bool useReverse) { - if (useReverse != m_useReverseReactions) { - m_jacobianMatrixState = JacobianMatrixState::STALE; - } m_useReverseReactions = useReverse; } @@ -572,7 +556,6 @@ namespace gridfire { if (depth != m_depth) { m_depth = depth; m_reactions = build_nuclear_network(comp, m_weakRateInterpolator, m_depth); - m_jacobianMatrixState = JacobianMatrixState::STALE; syncInternalMaps(); // Resync internal maps after changing the depth } else { LOG_DEBUG(m_logger, "Rebuild requested with the same depth. No changes made to the network."); @@ -599,6 +582,14 @@ namespace gridfire { return result; } + SpeciesStatus GraphEngine::getSpeciesStatus(const fourdst::atomic::Species &species) const { + if (m_networkSpeciesMap.contains(species.name())) { + return SpeciesStatus::ACTIVE; + } + return SpeciesStatus::NOT_PRESENT; + + } + StepDerivatives GraphEngine::calculateAllDerivativesUsingPrecomputation( const fourdst::composition::CompositionAbstract &comp, const std::vector &bare_rates, @@ -784,7 +775,6 @@ namespace gridfire { void GraphEngine::setScreeningModel(const screening::ScreeningType model) { m_screeningModel = screening::selectScreeningModel(model); m_screeningType = model; - m_jacobianMatrixState = JacobianMatrixState::STALE; // The screening model affects the jacobian so if its changed the jacobian must be made stale } screening::ScreeningType GraphEngine::getScreeningModel() const { @@ -828,7 +818,7 @@ namespace gridfire { ); } - void GraphEngine::generateJacobianMatrix( + NetworkJacobian GraphEngine::generateJacobianMatrix( const fourdst::composition::CompositionAbstract &comp, const double T9, const double rho @@ -847,7 +837,7 @@ namespace gridfire { std::vector adInput(numSpecies + 2, 0.0); // +2 for T9 and rho const std::vector& Y_dynamic = mutableComp.getMolarAbundanceVector(); for (size_t i = 0; i < numSpecies; ++i) { - adInput[i] = std::max(Y_dynamic[i], 1e-99); // regularize the jacobian... + adInput[i] = Y_dynamic[i]; } adInput[numSpecies] = T9; // T9 adInput[numSpecies + 1] = rho; // rho @@ -856,26 +846,34 @@ namespace gridfire { const std::vector dotY = m_rhsADFun.Jacobian(adInput); // 3. Pack jacobian vector into sparse matrix - m_jacobianMatrix.clear(); - // std::vector> columns; + Eigen::SparseMatrix jacobianMatrix(numSpecies, numSpecies); + std::vector > triplets; for (size_t i = 0; i < numSpecies; ++i) { - // std::vector colData; for (size_t j = 0; j < numSpecies; ++j) { - const double value = dotY[i * (numSpecies + 2) + j]; + double value = dotY[i * (numSpecies + 2) + j]; if (std::abs(value) > MIN_JACOBIAN_THRESHOLD || i == j) { // Always keep diagonal elements to avoid pathological stiffness - m_jacobianMatrix(i, j) = value; + if (i == j && value == 0) { + LOG_WARNING(m_logger, "While generating the Jacobian matrix, a zero diagonal element was encountered at index ({}, {}) (species: {}, abundance: {}). This may lead to numerical instability. Setting to -1 to avoid singularity", i, j, m_networkSpecies[i].name(), adInput[i]); + // value = -1.0; + } + triplets.emplace_back(i, j, value); } - // colData.push_back(value); } - // columns.push_back(std::make_unique>(std::to_string(i), colData)); } - // std::cout << utils::format_table("Jacobian after dense calculation", columns) << std::endl; - // exit(0); - LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); - m_jacobianMatrixState = JacobianMatrixState::READY_DENSE; + LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", jacobianMatrix.rows(), jacobianMatrix.cols()); + + auto index_to_species = [this](const size_t index) -> fourdst::atomic::Species { + if (index < m_networkSpecies.size()) { + return m_networkSpecies[index]; + } + throw std::out_of_range("Index out of range in index_to_species mapping."); + }; + jacobianMatrix.setFromTriplets(triplets.begin(), triplets.end()); + NetworkJacobian jac(jacobianMatrix, index_to_species); + return jac; } - void GraphEngine::generateJacobianMatrix( + NetworkJacobian GraphEngine::generateJacobianMatrix( const fourdst::composition::CompositionAbstract &comp, const double T9, const double rho, @@ -904,10 +902,10 @@ namespace gridfire { } // --- 3. Call the sparse reverse-mode implementation --- - generateJacobianMatrix(comp, T9, rho, sparsityPattern); + return generateJacobianMatrix(comp, T9, rho, sparsityPattern); } - void GraphEngine::generateJacobianMatrix( + NetworkJacobian GraphEngine::generateJacobianMatrix( const fourdst::composition::CompositionAbstract &comp, const double T9, const double rho, @@ -929,7 +927,7 @@ namespace gridfire { // } size_t i = 0; for (const auto& species: m_networkSpecies) { - double Yi = 0.0; + double Yi = 0.0; // Small floor to avoid issues with zero abundances if (comp.contains(species)) { Yi = comp.getMolarAbundance(species); } @@ -974,46 +972,26 @@ namespace gridfire { m_jac_work // Work vector for CppAD ); - // --- Convert the sparse Jacobian back to the Boost uBLAS format --- - m_jacobianMatrix.clear(); + Eigen::SparseMatrix jacobianMatrix(numSpecies, numSpecies); + std::vector > triplets; for (size_t k = 0; k < nnz; ++k) { const size_t row = jac_subset.row()[k]; const size_t col = jac_subset.col()[k]; const double value = jac_subset.val()[k]; 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 + triplets.emplace_back(row, col, value); } } - m_jacobianMatrixState = JacobianMatrixState::READY_SPARSE; - } - - double GraphEngine::getJacobianMatrixEntry( - const fourdst::atomic::Species& rowSpecies, - const fourdst::atomic::Species& colSpecies - ) const { - switch (m_jacobianMatrixState) { - case JacobianMatrixState::STALE: { - const std::string staleMsg = std::format("Cannot retrieve jacobian entry for row {}, column {} as jacobian matrix is stale (has not been regenerated since last network topology change)", rowSpecies.name(), colSpecies.name()); - throw exceptions::StaleJacobianError(staleMsg); + jacobianMatrix.setFromTriplets(triplets.begin(), triplets.end()); + auto index_to_species = [this](const size_t index) -> fourdst::atomic::Species { + if (index < m_networkSpecies.size()) { + return m_networkSpecies[index]; } - case JacobianMatrixState::UNINITIALIZED: { - const std::string unInitMsg = std::format("Cannot retrieve jacobian entry for row {}, column {} as jacobian matrix is uninitialized (will return all 0s)", rowSpecies.name(), colSpecies.name()); - throw exceptions::UninitializedJacobianError(unInitMsg); - } - case JacobianMatrixState::READY_DENSE: - [[fallthrough]]; - case JacobianMatrixState::READY_SPARSE: { - const size_t i = getSpeciesIndex(rowSpecies); - const size_t j = getSpeciesIndex(colSpecies); - return m_jacobianMatrix(i, j); - } - default: { - // Code should not be able to get into this state - const std::string msg = std::format("An unknown error has occurred while attempting to retrieve the jacobian element at row {}, column {}. This should be taken as a catastrophic failure and reported to GridFire developers.", rowSpecies.name(), colSpecies.name()); - throw exceptions::UnknownJacobianError(msg); - } - } + throw std::out_of_range("Index out of range in index_to_species mapping."); + }; + NetworkJacobian jac(jacobianMatrix, index_to_species); + return jac; } std::unordered_map GraphEngine::getNetReactionStoichiometry( @@ -1287,7 +1265,6 @@ namespace gridfire { } ); - // Extract the raw vector from the associative map std::vector> dependentVector; dependentVector.reserve(dydt.size() + 1); @@ -1302,8 +1279,7 @@ namespace gridfire { m_rhsADFun.Dependent(adInput, dependentVector); - LOG_TRACE_L1(m_logger, "AD tape recorded successfully for the RHS and Eps calculation. Number of independent variables: {}.", - adInput.size()); + LOG_TRACE_L1(m_logger, "AD tape recorded successfully for the RHS and Eps calculation. Number of independent variables: {}.", adInput.size()); } void GraphEngine::collectAtomicReverseRateAtomicBases() { diff --git a/src/lib/engine/procedures/priming.cpp b/src/lib/engine/procedures/priming.cpp index e0841663..42954d05 100644 --- a/src/lib/engine/procedures/priming.cpp +++ b/src/lib/engine/procedures/priming.cpp @@ -3,317 +3,79 @@ #include "fourdst/atomic/species.h" #include "fourdst/composition/utils.h" #include "gridfire/engine/views/engine_priming.h" -#include "gridfire/engine/procedures/construction.h" #include "gridfire/solver/solver.h" #include "gridfire/engine/engine_abstract.h" #include "gridfire/network.h" +#include "gridfire/exceptions/error_solver.h" #include "fourdst/logging/logging.h" +#include "gridfire/solver/strategies/CVODE_solver_strategy.h" #include "quill/Logger.h" #include "quill/LogMacros.h" -namespace { - bool isReactionIgnorable( - const gridfire::reaction::Reaction& reaction, - const std::optional>& reactionsTypesToIgnore - ) { - if (reactionsTypesToIgnore.has_value()) { - for (const auto& type : reactionsTypesToIgnore.value()) { - if (reaction.type() == type) { - return true; - } - } - } - return false; - } -} namespace gridfire { using fourdst::composition::Composition; using fourdst::atomic::Species; - const reaction::Reaction* findDominantCreationChannel ( - const DynamicEngine& engine, - const Species& species, - const Composition &comp, - const double T9, - const double rho, - const std::optional> &reactionsTypesToIgnore - ) { - const reaction::Reaction* dominateReaction = nullptr; - double maxFlow = -1.0; - for (const auto& reaction : engine.getNetworkReactions()) { - if (isReactionIgnorable(*reaction, reactionsTypesToIgnore)) continue; - - if (reaction->contains(species) && reaction->stoichiometry(species) > 0) { - const double flow = engine.calculateMolarReactionFlow(*reaction, comp, T9, rho); - if (flow > maxFlow) { - maxFlow = flow; - dominateReaction = reaction.get(); - } - } - } - return dominateReaction; - } - - - - /** - * @brief Primes absent species in the network to their equilibrium abundances using a robust, two-stage approach. - * - * @details This function implements a robust network priming algorithm that avoids the pitfalls of - * sequential, one-by-one priming. The previous, brittle method could allow an early priming - * reaction to consume all of a shared reactant, starving later reactions. This new, two-stage - * method ensures that all priming reactions are considered collectively, competing for the - * same limited pool of initial reactants in a physically consistent manner. - * - * The algorithm proceeds in three main stages: - * 1. **Calculation Stage:** It first loops through all species that need priming. For each one, - * it calculates its theoretical equilibrium mass fraction and identifies the dominant - * creation channel. Crucially, it *does not* modify any abundances at this stage. Instead, - * it stores these calculations as a list of "mass transfer requests". - * - * 2. **Collective Scaling Stage:** It then processes the full list of requests to determine the - * total "debit" required from each reactant. By comparing these total debits to the - * initially available mass of each reactant, it calculates a single, global `scalingFactor`. - * If any reactant is overdrawn, this factor will be less than 1.0, ensuring that no - * reactant's abundance can go negative. - * - * 3. **Application Stage:** Finally, it loops through the requests again, applying the mass - * transfers. Each calculated equilibrium mass fraction and corresponding reactant debit is - * multiplied by the global `scalingFactor` before being applied to the final composition. - * This ensures that if resources are limited, all primed species are scaled down proportionally. - * - * @param netIn Input network data containing initial composition, temperature, and density. - * @param engine DynamicEngine used to build and evaluate the reaction network. - * @param ignoredReactionTypes Types of reactions to ignore during priming (e.g., weak reactions). - * @return PrimingReport encapsulating the results of the priming operation, including the new - * robustly primed composition. - */ PrimingReport primeNetwork( - const NetIn& netIn, - GraphEngine& engine, - const std::optional>& ignoredReactionTypes + const NetIn& netIn, + GraphEngine& engine, + const std::optional>& ignoredReactionTypes ) { - auto logger = LogManager::getInstance().getLogger("log"); + const auto logger = LogManager::getInstance().getLogger("log"); + solver::CVODESolverStrategy integrator(engine); + integrator.set_stdout_logging_enabled(false); + NetIn solverInput(netIn); - // --- Initial Setup --- - // Identify all species with zero initial abundance that need to be primed. - std::vector speciesToPrime; - for (const auto &[sp, y]: netIn.composition) { - if (y == 0.0) { - speciesToPrime.push_back(sp); - } - } + solverInput.tMax = 1e-15; + solverInput.temperature = 1e7; - // Sort priming species by mass number, lightest to heaviest. - std::ranges::sort(speciesToPrime, [](const Species& a, const Species& b) { - return a.mass() < b.mass(); - }); - - LOG_DEBUG(logger, "Priming {} species in the network.", speciesToPrime.size()); - - // If no species need priming, return immediately. + LOG_INFO(logger, "Short timescale ({}) network ignition started.", solverInput.tMax); PrimingReport report; - if (speciesToPrime.empty()) { - report.primedComposition = netIn.composition; - report.success = true; - report.status = PrimingReportStatus::NO_SPECIES_TO_PRIME; - return report; - } - - const double T9 = netIn.temperature / 1e9; - const double rho = netIn.density; - const auto initialReactionSet = engine.getNetworkReactions(); - - report.status = PrimingReportStatus::FULL_SUCCESS; - report.success = true; - - // Build full set of species - std::set allSpecies; - for (const auto &sp: netIn.composition | std::views::keys) { - allSpecies.insert(sp); - } - for (const auto& sp : speciesToPrime) { - allSpecies.insert(sp); - } - - // Rebuild the engine with the full network to ensure all possible creation channels are available. - engine.rebuild(netIn.composition, NetworkBuildDepth::Full); - - // Initialize mutable molar abundances for all species - std::map molarAbundances; - for (const auto& sp : allSpecies) { - molarAbundances[sp] = netIn.composition.contains(sp) ? netIn.composition.getMolarAbundance(sp) : 0.0; - } - - // --- Prime Each Species --- - // Since molar abundances are independent, we can directly calculate and apply changes - std::unordered_map totalMolarAbundanceChanges; - - for (const auto& primingSpecies : speciesToPrime) { - // Create a temporary composition reflecting the current state for rate calculations. - Composition tempComp(allSpecies); - for (const auto& [sp, abundance] : molarAbundances) { - tempComp.setMolarAbundance(sp, abundance); - } - - NetworkPrimingEngineView primer(primingSpecies, engine); - if (primer.getNetworkReactions().size() == 0) { - LOG_ERROR(logger, "No priming reactions found for species {}.", primingSpecies.name()); - report.success = false; - report.status = PrimingReportStatus::FAILED_TO_FIND_PRIMING_REACTIONS; - continue; - } - - const double destructionRateConstant = calculateDestructionRateConstant( - primer, - primingSpecies, - tempComp, - T9, - rho, - ignoredReactionTypes + try { + const NetOut netOut = integrator.evaluate(solverInput, false); + LOG_INFO(logger, "Network ignition completed."); + LOG_TRACE_L2( + logger, + "After ignition composition is {}", + [netOut, netIn]() -> std::string { + std::stringstream ss; + size_t i = 0; + for (const auto& [species, abundance] : netOut.composition) { + ss << species.name() << ": " << abundance << " (prior: " << netIn.composition.getMolarAbundance(species); + ss << ", fractional change: " << (abundance - netIn.composition.getMolarAbundance(species)) / netIn.composition.getMolarAbundance(species) * 100.0 << "%)"; + if (i < netOut.composition.size() - 1) { + ss << ", "; + } + ++i; + } + return ss.str(); + }() ); - - double equilibriumMolarAbundance = 0.0; - std::vector reactants; - - if (destructionRateConstant > 1e-99) { - const double creationRate = calculateCreationRate( - primer, - primingSpecies, - tempComp, - T9, - rho, - ignoredReactionTypes - ); - - // Equilibrium: creation rate = destruction rate constant * molar abundance - equilibriumMolarAbundance = creationRate / destructionRateConstant; - - if (std::isnan(equilibriumMolarAbundance)) { - equilibriumMolarAbundance = 0.0; + report.primedComposition = netOut.composition; + std::unordered_set unprimedSpecies; + double minAbundance = std::numeric_limits::max(); + for (const auto& [sp, y] : report.primedComposition) { + if (y == 0) { + unprimedSpecies.insert(sp); } - - if (const reaction::Reaction* dominantChannel = findDominantCreationChannel( - primer, - primingSpecies, - tempComp, - T9, - rho, - ignoredReactionTypes) - ) { - reactants = dominantChannel->reactants(); - } else { - LOG_TRACE_L1(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name()); - report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL; - reactants.clear(); // Use fallback - } - } else { - LOG_TRACE_L2(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name()); - // For species with no destruction, use a tiny fallback abundance - equilibriumMolarAbundance = 1e-40; - } - - // Add the equilibrium molar abundance to the primed species - molarAbundances.at(primingSpecies) += equilibriumMolarAbundance; - totalMolarAbundanceChanges[primingSpecies] += equilibriumMolarAbundance; - - // Subtract from reactants proportionally to their stoichiometry - if (!reactants.empty()) { - const double totalStoichiometry = static_cast(reactants.size()); - const double abundancePerReactant = equilibriumMolarAbundance / totalStoichiometry; - - for (const auto& reactant : reactants) { - LOG_TRACE_L1(logger, "Transferring {:.4e} molar abundance from {} to {} during priming.", - abundancePerReactant, reactant.name(), primingSpecies.name()); - - if (!molarAbundances.contains(reactant)) { - continue; - } - molarAbundances.at(reactant) -= abundancePerReactant; - totalMolarAbundanceChanges[reactant] -= abundancePerReactant; - - // Ensure non-negative abundances - if (molarAbundances.at(reactant) < 0.0) { - LOG_WARNING(logger, "Species {} went negative during priming. Clamping to zero.", reactant.name()); - totalMolarAbundanceChanges[reactant] += molarAbundances.at(reactant); // Adjust change to reflect clamp - molarAbundances.at(reactant) = 0.0; - } + if (y != 0 && y < minAbundance) { + minAbundance = y; } } + double abundanceForUnprimedSpecies = minAbundance / 1e10; + for (const auto& sp : unprimedSpecies) { + LOG_TRACE_L1(logger, "Clamping Species {}: initial abundance {}, primed abundance {} to {}", sp.name(), netIn.composition.getMolarAbundance(sp), report.primedComposition.getMolarAbundance(sp), abundanceForUnprimedSpecies); + report.primedComposition.setMolarAbundance(sp, abundanceForUnprimedSpecies); + } + report.success = true; + report.status = PrimingReportStatus::SUCCESS; + } catch (const exceptions::SolverError& e) { + LOG_ERROR(logger, "Failed to prime network: solver failure encountered: {}", e.what()); + std::rethrow_exception(std::current_exception()); } - - // --- Final Composition Construction --- - std::vector final_species; - std::vector final_molar_abundances; - - for (const auto& [species, abundance] : molarAbundances) { - final_species.emplace_back(species); - final_molar_abundances.push_back(std::max(0.0, abundance)); - - LOG_TRACE_L1(logger, "After priming, species {} has molar abundance {:.4e} (had {:0.4e} prior to priming).", - species.name(), - std::max(0.0, abundance), - netIn.composition.getMolarAbundance(species)); - } - - Composition primedComposition(final_species, final_molar_abundances); - - report.primedComposition = primedComposition; - - // Convert molar abundance changes to mass fraction changes for reporting - for (const auto& [species, molarChange] : totalMolarAbundanceChanges) { - double massFractionChange = molarChange * species.mass(); - report.massFractionChanges.emplace_back(species, massFractionChange); - } - - // Restore the engine to its original, smaller network state. - engine.setNetworkReactions(initialReactionSet); - return report; } - - double calculateDestructionRateConstant( - const DynamicEngine& engine, - const Species& species, - const Composition& composition, - const double T9, - const double rho, - const std::optional> &reactionTypesToIgnore - ) { - Composition unrestrictedComp(composition); - unrestrictedComp.setMolarAbundance(species, 1.0); - double destructionRateConstant = 0.0; - for (const auto& reaction: engine.getNetworkReactions()) { - if (isReactionIgnorable(*reaction, reactionTypesToIgnore)) continue; - - const int stoichiometry = reaction->stoichiometry(species); - if (stoichiometry < 0) { - destructionRateConstant += std::abs(stoichiometry) * engine.calculateMolarReactionFlow(*reaction, unrestrictedComp, T9, rho); - } - } - return destructionRateConstant; - } - - double calculateCreationRate( - const DynamicEngine& engine, - const Species& species, - const Composition& composition, - const double T9, - const double rho, - const std::optional> &reactionTypesToIgnore - ) { - double creationRate = 0.0; - for (const auto& reaction: engine.getNetworkReactions()) { - if (isReactionIgnorable(*reaction, reactionTypesToIgnore)) continue; - - const int stoichiometry = reaction->stoichiometry(species); - if (stoichiometry > 0) { - creationRate += stoichiometry * engine.calculateMolarReactionFlow(*reaction, composition, T9, rho); - } - } - return creationRate; - } } \ No newline at end of file