From 69bd2cd4667f2ea4ea4642e3d2f077f58ff8924f Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 14 Jul 2025 14:50:49 -0400 Subject: [PATCH] feat(GridFire): Much more robust starting network GraphEngine now can initialize with a much more robust set of reactions (including the entire reaction set). The jacobian can still be efficiently evaluated using CppAD's sparse jacobian feature. Further, the primeing network has been signifiganty enhanced to handle much hotter termperatures --- .../include/gridfire/engine/engine_abstract.h | 29 ++- .../include/gridfire/engine/engine_graph.h | 64 +++++-- .../gridfire/engine/procedures/construction.h | 21 +- .../gridfire/engine/procedures/priming.h | 8 +- .../include/gridfire/engine/types/building.h | 16 ++ .../include/gridfire/engine/types/reporting.h | 46 ++++- .../gridfire/engine/views/engine_adaptive.h | 5 + .../gridfire/engine/views/engine_defined.h | 2 + .../gridfire/engine/views/engine_multiscale.h | 18 +- .../gridfire/engine/views/engine_priming.h | 3 - src/network/include/gridfire/network.h | 2 - src/network/lib/engine/engine_graph.cpp | 168 ++++++++++++++-- .../lib/engine/procedures/construction.cpp | 116 ++++++++++- src/network/lib/engine/procedures/priming.cpp | 181 ++++++++++++++++-- .../lib/engine/views/engine_adaptive.cpp | 4 + .../lib/engine/views/engine_defined.cpp | 6 +- .../lib/engine/views/engine_multiscale.cpp | 49 +++-- .../lib/engine/views/engine_priming.cpp | 3 +- src/network/lib/network.cpp | 26 --- src/network/lib/solver/solver.cpp | 11 +- src/network/meson.build | 3 +- 21 files changed, 650 insertions(+), 131 deletions(-) create mode 100644 src/network/include/gridfire/engine/types/building.h diff --git a/src/network/include/gridfire/engine/engine_abstract.h b/src/network/include/gridfire/engine/engine_abstract.h index b7aa6a9d..584062dd 100644 --- a/src/network/include/gridfire/engine/engine_abstract.h +++ b/src/network/include/gridfire/engine/engine_abstract.h @@ -5,8 +5,12 @@ #include "gridfire/screening/screening_abstract.h" #include "gridfire/screening/screening_types.h" +#include "gridfire/engine/types/reporting.h" +#include "gridfire/engine/types/building.h" + #include #include +#include /** * @file engine_abstract.h @@ -55,6 +59,8 @@ namespace gridfire { T nuclearEnergyGenerationRate = T(0.0); ///< Specific energy generation rate (e.g., erg/g/s). }; + using SparsityPattern = std::vector>; + /** * @brief Abstract base class for a reaction network engine. * @@ -132,9 +138,19 @@ namespace gridfire { */ virtual void generateJacobianMatrix( const std::vector& Y_dynamic, - double T9, double rho + double T9, + double rho ) = 0; + virtual void generateJacobianMatrix( + const std::vector& Y_dynamic, + double T9, + double rho, + const SparsityPattern& sparsityPattern + ) { + throw std::logic_error("Sparsity pattern not supported by this engine."); + } + /** * @brief Get an entry from the previously generated Jacobian matrix. * @@ -149,6 +165,7 @@ namespace gridfire { int j ) const = 0; + /** * @brief Generate the stoichiometry matrix for the network. * @@ -269,5 +286,15 @@ namespace gridfire { [[nodiscard]] virtual int getSpeciesIndex(const fourdst::atomic::Species &species) const = 0; [[nodiscard]] virtual std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const = 0; + + [[nodiscard]] virtual PrimingReport primeEngine(const NetIn &netIn) = 0; + + [[nodiscard]] virtual BuildDepthType getDepth() const { + throw std::logic_error("Network depth not supported by this engine."); + } + + virtual void rebuild(const fourdst::composition::Composition& comp, BuildDepthType depth) { + throw std::logic_error("Setting network depth not supported by this engine."); + } }; } \ No newline at end of file diff --git a/src/network/include/gridfire/engine/engine_graph.h b/src/network/include/gridfire/engine/engine_graph.h index 100c49c9..ac162d42 100644 --- a/src/network/include/gridfire/engine/engine_graph.h +++ b/src/network/include/gridfire/engine/engine_graph.h @@ -11,7 +11,7 @@ #include "gridfire/screening/screening_abstract.h" #include "gridfire/screening/screening_types.h" #include "gridfire/partition/partition_abstract.h" -#include "gridfire/partition/partition_ground.h" +#include "gridfire/engine/procedures/construction.h" #include #include @@ -22,6 +22,12 @@ #include #include "cppad/cppad.hpp" +#include "cppad/utility/sparse_rc.hpp" +#include "cppad/speed/sparse_jac_fun.hpp" + +#include "procedures/priming.h" + + #include "quill/LogMacros.h" // PERF: The function getNetReactionStoichiometry returns a map of species to their stoichiometric coefficients for a given reaction. @@ -30,6 +36,7 @@ // REACLIBReactions are quite large data structures, so this could be a performance bottleneck. namespace gridfire { + static bool s_debug = false; // Global debug flag for the GraphEngine /** * @brief Alias for CppAD AD type for double precision. * @@ -103,11 +110,16 @@ namespace gridfire { * * @see build_reaclib_nuclear_network */ - explicit GraphEngine(const fourdst::composition::Composition &composition); + explicit GraphEngine( + const fourdst::composition::Composition &composition, + const BuildDepthType = NetworkBuildDepth::Full + ); - explicit GraphEngine(const fourdst::composition::Composition &composition, - const partition::PartitionFunction& partitionFunction - ); + explicit GraphEngine( + const fourdst::composition::Composition &composition, + const partition::PartitionFunction& partitionFunction, + const BuildDepthType buildDepth = NetworkBuildDepth::Full + ); /** * @brief Constructs a GraphEngine from a set of reactions. @@ -157,6 +169,13 @@ namespace gridfire { const double rho ) override; + void generateJacobianMatrix( + const std::vector &Y_dynamic, + double T9, + double rho, + const SparsityPattern &sparsityPattern + ) override; + /** * @brief Generates the stoichiometry matrix for the network. * @@ -322,28 +341,34 @@ namespace gridfire { double T9 ) const; - double calculateReverseRateTwoBody( + [[nodiscard]] double calculateReverseRateTwoBody( const reaction::Reaction &reaction, const double T9, const double forwardRate, const double expFactor ) const; - double calculateReverseRateTwoBodyDerivative( + [[nodiscard]] double calculateReverseRateTwoBodyDerivative( const reaction::Reaction &reaction, const double T9, const double reverseRate ) const; - bool isUsingReverseReactions() const; + [[nodiscard]] bool isUsingReverseReactions() const; void setUseReverseReactions(bool useReverse); - int getSpeciesIndex( + [[nodiscard]] int getSpeciesIndex( const fourdst::atomic::Species& species ) const override; - std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const override; + [[nodiscard]] std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const override; + + [[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override; + + [[nodiscard]] BuildDepthType getDepth() const override; + + void rebuild(const fourdst::composition::Composition& comp, const BuildDepthType depth) override; @@ -427,6 +452,8 @@ namespace gridfire { boost::numeric::ublas::compressed_matrix m_jacobianMatrix; ///< Jacobian matrix (species x species). CppAD::ADFun m_rhsADFun; ///< CppAD function for the right-hand side of the ODE. + 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::vector> m_atomicReverseRates; screening::ScreeningType m_screeningType = screening::ScreeningType::BARE; ///< Screening type for the reaction network. Default to no screening. @@ -436,6 +463,8 @@ namespace gridfire { bool m_useReverseReactions = true; ///< Flag to enable or disable reverse reactions. If false, only forward reactions are considered. + BuildDepthType m_depth; + std::vector m_precomputedReactions; ///< Precomputed reactions for efficiency. std::unique_ptr m_partitionFunction; ///< Partition function for the network. @@ -636,6 +665,14 @@ namespace gridfire { if (!m_useReverseReactions) { return static_cast(0.0); // If reverse reactions are not used, return zero } + s_debug = false; + if (reaction.peName() == "p(p,e+)d" || reaction.peName() =="d(d,n)he3" || reaction.peName() == "c12(p,g)n13") { + if constexpr (!std::is_same_v) { + s_debug = true; + std::cout << "Calculating reverse molar flow for reaction: " << reaction.peName() << std::endl; + std::cout << "\tT9: " << T9 << ", rho: " << rho << std::endl; + } + } T reverseMolarFlow = static_cast(0.0); if (reaction.qValue() != 0.0) { @@ -655,6 +692,9 @@ namespace gridfire { } } else { // A,B If not calling with an AD type, calculate the reverse rate directly + if (s_debug) { + std::cout << "\tUsing double overload\n"; + } reverseRateConstant = calculateReverseRate(reaction, T9); } @@ -673,7 +713,7 @@ namespace gridfire { // E. Calculate the abundance term T productAbundanceTerm = static_cast(1.0); for (const auto& [species, count] : productCounts) { - const int speciesIndex = m_speciesToIndexMap.at(species); + const unsigned long speciesIndex = m_speciesToIndexMap.at(species); productAbundanceTerm *= CppAD::pow(Y[speciesIndex], count); } @@ -757,7 +797,7 @@ namespace gridfire { ss << "Forward: " << forwardMolarReactionFlow << ", Reverse: " << reverseMolarFlow << ", Net: " << molarReactionFlow; - LOG_DEBUG( + LOG_TRACE_L2( m_logger, "Reaction: {}, {}", reaction.peName(), diff --git a/src/network/include/gridfire/engine/procedures/construction.h b/src/network/include/gridfire/engine/procedures/construction.h index 39859633..525eea53 100644 --- a/src/network/include/gridfire/engine/procedures/construction.h +++ b/src/network/include/gridfire/engine/procedures/construction.h @@ -1,8 +1,17 @@ -// -// Created by Emily Boudreaux on 7/14/25. -// +#pragma once -#ifndef CONSTRUCTION_H -#define CONSTRUCTION_H +#include "gridfire/reaction/reaction.h" +#include "gridfire/engine/types/building.h" -#endif //CONSTRUCTION_H +#include "fourdst/composition/composition.h" + +#include + +namespace gridfire { + + reaction::LogicalReactionSet build_reaclib_nuclear_network( + const fourdst::composition::Composition &composition, + BuildDepthType maxLayers = NetworkBuildDepth::Full, + bool reverse = false + ); +} \ No newline at end of file diff --git a/src/network/include/gridfire/engine/procedures/priming.h b/src/network/include/gridfire/engine/procedures/priming.h index d3cc91e1..22f365a0 100644 --- a/src/network/include/gridfire/engine/procedures/priming.h +++ b/src/network/include/gridfire/engine/procedures/priming.h @@ -6,9 +6,15 @@ #include "fourdst/composition/composition.h" #include "fourdst/composition/atomicSpecies.h" +#include +#include +#include + namespace gridfire { - fourdst::composition::Composition primeNetwork( + + + PrimingReport primeNetwork( const NetIn&, DynamicEngine& engine ); diff --git a/src/network/include/gridfire/engine/types/building.h b/src/network/include/gridfire/engine/types/building.h new file mode 100644 index 00000000..1fd1a7c7 --- /dev/null +++ b/src/network/include/gridfire/engine/types/building.h @@ -0,0 +1,16 @@ +#pragma once + +#include + +namespace gridfire { + enum class NetworkBuildDepth { + Full = -1, + Shallow = 1, + SecondOrder = 2, + ThirdOrder = 3, + FourthOrder = 4, + FifthOrder = 5 + }; + + using BuildDepthType = std::variant; +} diff --git a/src/network/include/gridfire/engine/types/reporting.h b/src/network/include/gridfire/engine/types/reporting.h index 238582df..52a55c3a 100644 --- a/src/network/include/gridfire/engine/types/reporting.h +++ b/src/network/include/gridfire/engine/types/reporting.h @@ -1,8 +1,42 @@ -// -// Created by Emily Boudreaux on 7/14/25. -// +#pragma once -#ifndef REPORTING_H -#define REPORTING_H +#include +#include +#include -#endif //REPORTING_H +namespace gridfire { + 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 + }; + + 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"} + }; + + struct PrimingReport { + fourdst::composition::Composition primedComposition; + std::vector> massFractionChanges; ///< Species and their destruction/creation rates + bool success; + PrimingReportStatus status; + + friend std::ostream& operator<<(std::ostream& os, const PrimingReport& report) { + std::stringstream ss; + std::string successStr = report.success ? "true" : "false"; + ss << "PrimingReport(success=" << successStr + << ", status=" << PrimingReportStatusStrings[report.status] << ")"; + return os << ss.str(); + } + }; +} \ No newline at end of file diff --git a/src/network/include/gridfire/engine/views/engine_adaptive.h b/src/network/include/gridfire/engine/views/engine_adaptive.h index 510bc6c3..71f2e8b9 100644 --- a/src/network/include/gridfire/engine/views/engine_adaptive.h +++ b/src/network/include/gridfire/engine/views/engine_adaptive.h @@ -9,6 +9,9 @@ #include "fourdst/config/config.h" #include "fourdst/logging/logging.h" +#include "gridfire/engine/procedures/priming.h" +#include "gridfire/engine/procedures/construction.h" + #include "quill/Logger.h" namespace gridfire { @@ -260,6 +263,8 @@ namespace gridfire { [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override; [[nodiscard]] std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const override; + + [[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override; private: using Config = fourdst::config::Config; using LogManager = fourdst::logging::LogManager; diff --git a/src/network/include/gridfire/engine/views/engine_defined.h b/src/network/include/gridfire/engine/views/engine_defined.h index 454199b5..76303196 100644 --- a/src/network/include/gridfire/engine/views/engine_defined.h +++ b/src/network/include/gridfire/engine/views/engine_defined.h @@ -161,6 +161,8 @@ namespace gridfire{ [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override; [[nodiscard]] std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const override; + + [[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override; protected: bool m_isStale = true; DynamicEngine& m_baseEngine; diff --git a/src/network/include/gridfire/engine/views/engine_multiscale.h b/src/network/include/gridfire/engine/views/engine_multiscale.h index bf66980c..7eaeeb28 100644 --- a/src/network/include/gridfire/engine/views/engine_multiscale.h +++ b/src/network/include/gridfire/engine/views/engine_multiscale.h @@ -83,18 +83,20 @@ namespace gridfire { [[nodiscard]] std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const override; - std::vector getFastSpecies() const; - const std::vector& getDynamicSpecies() const; + [[nodiscard]] PrimingReport primeEngine(const NetIn &netIn) override; - void equilibrateNetwork( - const std::vector& Y, + [[nodiscard]] std::vector getFastSpecies() const; + [[nodiscard]] const std::vector& getDynamicSpecies() const; + + fourdst::composition::Composition equilibrateNetwork( + const std::vector &Y, double T9, double rho, double dt_control ); - void equilibrateNetwork( - const NetIn& netIn, + fourdst::composition::Composition equilibrateNetwork( + const NetIn &netIn, const double dt_control ); @@ -137,8 +139,8 @@ namespace gridfire { m_rho(rho), m_Y_scale(Y_scale) {} - int values() const { return m_qse_solve_indices.size(); } - int inputs() const { return m_qse_solve_indices.size(); } + [[nodiscard]] int values() const { return m_qse_solve_indices.size(); } + [[nodiscard]] int inputs() const { return m_qse_solve_indices.size(); } int operator()(const InputType& v_qse, OutputType& f_qse) const; int df(const InputType& v_qse, JacobianType& J_qse) const; diff --git a/src/network/include/gridfire/engine/views/engine_priming.h b/src/network/include/gridfire/engine/views/engine_priming.h index d38471df..93c88706 100644 --- a/src/network/include/gridfire/engine/views/engine_priming.h +++ b/src/network/include/gridfire/engine/views/engine_priming.h @@ -20,13 +20,10 @@ namespace gridfire { public: NetworkPrimingEngineView(const std::string& primingSymbol, DynamicEngine& baseEngine); NetworkPrimingEngineView(const fourdst::atomic::Species& primingSpecies, DynamicEngine& baseEngine); - const std::vector& getPrimingReactionNames() const { return m_peNames; } - const fourdst::atomic::Species& getPrimingSpecies() const { return m_primingSpecies; } private: quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - std::vector m_peNames; ///< Names of the priming reactions. fourdst::atomic::Species m_primingSpecies; ///< The priming species, if specified. private: std::vector constructPrimingReactionSet( diff --git a/src/network/include/gridfire/network.h b/src/network/include/gridfire/network.h index 74c87e65..c5ee5f33 100644 --- a/src/network/include/gridfire/network.h +++ b/src/network/include/gridfire/network.h @@ -104,7 +104,5 @@ namespace gridfire { }; - reaction::LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, bool reverse); - } // namespace nuclearNetwork diff --git a/src/network/lib/engine/engine_graph.cpp b/src/network/lib/engine/engine_graph.cpp index 2de28e86..5e762c4b 100644 --- a/src/network/lib/engine/engine_graph.cpp +++ b/src/network/lib/engine/engine_graph.cpp @@ -2,6 +2,9 @@ #include "gridfire/reaction/reaction.h" #include "gridfire/network.h" #include "gridfire/screening/screening_types.h" +#include "gridfire/engine/procedures/priming.h" +#include "gridfire/partition/partition_ground.h" +#include "gridfire/engine/procedures/construction.h" #include "fourdst/composition/species.h" #include "fourdst/composition/atomicSpecies.h" @@ -22,17 +25,24 @@ #include +#include "cppad/cppad.hpp" +#include "cppad/utility/sparse_rc.hpp" +#include "cppad/utility/sparse_rcv.hpp" + namespace gridfire { GraphEngine::GraphEngine( - const fourdst::composition::Composition &composition - ): GraphEngine(composition, partition::GroundStatePartitionFunction()) {} + const fourdst::composition::Composition &composition, + const BuildDepthType buildDepth + ): GraphEngine(composition, partition::GroundStatePartitionFunction(), buildDepth) {} GraphEngine::GraphEngine( const fourdst::composition::Composition &composition, - const partition::PartitionFunction& partitionFunction) : - m_reactions(build_reaclib_nuclear_network(composition, false)), - m_partitionFunction(partitionFunction.clone()) // Clone the partition function to ensure ownership + const partition::PartitionFunction& partitionFunction, + const BuildDepthType buildDepth) : + m_reactions(build_reaclib_nuclear_network(composition, buildDepth, false)), + m_partitionFunction(partitionFunction.clone()), + m_depth(buildDepth) { syncInternalMaps(); } @@ -77,6 +87,16 @@ namespace gridfire { generateStoichiometryMatrix(); reserveJacobianMatrix(); recordADTape(); + + const size_t n = m_rhsADFun.Domain(); + const size_t m = m_rhsADFun.Range(); + + std::vector select_domain(n, true); + std::vector select_range(m, true); + + m_rhsADFun.subgraph_sparsity(select_domain, select_range, false, m_full_jacobian_sparsity_pattern); + m_jac_work.clear(); + precomputeNetwork(); LOG_INFO(m_logger, "Internal maps synchronized. Network contains {} species and {} reactions.", m_networkSpecies.size(), m_reactions.size()); @@ -243,6 +263,14 @@ namespace gridfire { return 0.0; // If reverse reactions are not used, return 0.0 } const double expFactor = std::exp(-reaction.qValue() / (m_constants.kB * T9)); + if (s_debug) { + std::cout << "\texpFactor = exp(-qValue/(kB * T9))\n"; + std::cout << "\texpFactor: " << expFactor << " for reaction: " << reaction.peName() << std::endl; + std::cout << "\tQ-value: " << reaction.qValue() << " at T9: " << T9 << std::endl; + std::cout << "\tT9: " << T9 << std::endl; + std::cout << "\tkB * T9: " << m_constants.kB * T9 << std::endl; + std::cout << "\tqValue/(kB * T9): " << reaction.qValue() / (m_constants.kB * T9) << std::endl; + } double reverseRate = 0.0; const double forwardRate = reaction.calculate_rate(T9); @@ -393,6 +421,47 @@ namespace gridfire { return Y; // Return the vector of molar abundances } + PrimingReport GraphEngine::primeEngine(const NetIn &netIn) { + NetIn fullNetIn; + fourdst::composition::Composition composition; + + std::vector symbols; + symbols.reserve(m_networkSpecies.size()); + for (const auto &symbol: m_networkSpecies) { + symbols.emplace_back(symbol.name()); + } + composition.registerSymbol(symbols); + for (const auto& [symbol, entry] : netIn.composition) { + if (m_networkSpeciesMap.contains(symbol)) { + composition.setMassFraction(symbol, entry.mass_fraction()); + } else { + composition.setMassFraction(symbol, 0.0); + } + } + composition.finalize(true); + fullNetIn.composition = composition; + fullNetIn.temperature = netIn.temperature; + fullNetIn.density = netIn.density; + + auto primingReport = primeNetwork(fullNetIn, *this); + + return primingReport; + } + + BuildDepthType GraphEngine::getDepth() const { + return m_depth; + } + + void GraphEngine::rebuild(const fourdst::composition::Composition& comp, const BuildDepthType depth) { + if (depth != m_depth) { + m_depth = depth; + m_reactions = build_reaclib_nuclear_network(comp, m_depth, false); + 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."); + } + } + StepDerivatives GraphEngine::calculateAllDerivativesUsingPrecomputation( const std::vector &Y_in, const std::vector &bare_rates, @@ -414,23 +483,23 @@ namespace gridfire { molarReactionFlows.reserve(m_precomputedReactions.size()); for (const auto& precomp : m_precomputedReactions) { - double abundanceProduct = 1.0; - bool below_threshold = false; + double forwardAbundanceProduct = 1.0; + // bool below_threshold = false; for (size_t i = 0; i < precomp.unique_reactant_indices.size(); ++i) { const size_t reactantIndex = precomp.unique_reactant_indices[i]; const int power = precomp.reactant_powers[i]; - const double abundance = Y_in[reactantIndex]; - if (abundance < MIN_ABUNDANCE_THRESHOLD) { - below_threshold = true; - break; - } + // const double abundance = Y_in[reactantIndex]; + // if (abundance < MIN_ABUNDANCE_THRESHOLD) { + // below_threshold = true; + // break; + // } - abundanceProduct *= std::pow(Y_in[reactantIndex], power); - } - if (below_threshold) { - molarReactionFlows.push_back(0.0); - continue; // Skip this reaction if any reactant is below the abundance threshold + forwardAbundanceProduct *= std::pow(Y_in[reactantIndex], power); } + // if (below_threshold) { + // molarReactionFlows.push_back(0.0); + // continue; // Skip this reaction if any reactant is below the abundance threshold + // } const double bare_rate = bare_rates[precomp.reaction_index]; const double screeningFactor = screeningFactors[precomp.reaction_index]; @@ -441,8 +510,8 @@ namespace gridfire { screeningFactor * bare_rate * precomp.symmetry_factor * - abundanceProduct * - std::pow(rho, numReactants > 1 ? numReactants - 1 : 0); + forwardAbundanceProduct * + std::pow(rho, numReactants > 1 ? numReactants - 1 : 0.0); double reverseMolarReactionFlow = 0.0; if (precomp.reverse_symmetry_factor != 0.0 and m_useReverseReactions) { @@ -455,7 +524,7 @@ namespace gridfire { bare_reverse_rate * precomp.reverse_symmetry_factor * reverseAbundanceProduct * - std::pow(rho, numProducts > 1 ? numProducts - 1 : 0); + std::pow(rho, numProducts > 1 ? numProducts - 1 : 0.0); } molarReactionFlows.push_back(forwardMolarReactionFlow - reverseMolarReactionFlow); @@ -474,7 +543,7 @@ namespace gridfire { const int stoichiometricCoefficient = precomp.stoichiometric_coefficients[i]; // Update the derivative for this species - result.dydt[speciesIndex] += static_cast(stoichiometricCoefficient) * R_j / rho; + result.dydt[speciesIndex] += static_cast(stoichiometricCoefficient) * R_j; } } @@ -641,6 +710,63 @@ namespace gridfire { LOG_TRACE_L1(m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); } + void GraphEngine::generateJacobianMatrix( + const std::vector &Y_dynamic, + const double T9, + const double rho, + const SparsityPattern &sparsityPattern + ) { + // --- Pack the input variables into a vector for CppAD --- + const size_t numSpecies = m_networkSpecies.size(); + std::vector x(numSpecies + 2, 0.0); + for (size_t i = 0; i < numSpecies; ++i) { + x[i] = Y_dynamic[i]; + } + x[numSpecies] = T9; + x[numSpecies + 1] = rho; + + // --- Convert into CppAD Sparsity pattern --- + const size_t nnz = sparsityPattern.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; + } + + std::vector values(nnz); + const size_t num_rows_jac = numSpecies; + 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); + for (size_t k = 0; k < nnz; ++k) { + CppAD_sparsity_pattern.set(k, sparsityPattern[k].first, sparsityPattern[k].second); + } + + CppAD::sparse_rcv, std::vector> jac_subset(CppAD_sparsity_pattern); + + m_rhsADFun.sparse_jac_rev( + x, + jac_subset, // Sparse Jacobian output + m_full_jacobian_sparsity_pattern, + "cppad", + m_jac_work // Work vector for CppAD + ); + + // --- Convert the sparse Jacobian back to the Boost uBLAS format --- + m_jacobianMatrix.clear(); + 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) { + m_jacobianMatrix(row, col) = value; // Insert into the sparse matrix + } + } + } + double GraphEngine::getJacobianMatrixEntry(const int i, const int j) const { LOG_TRACE_L3(m_logger, "Getting jacobian matrix entry for {},{} = {}", i, j, m_jacobianMatrix(i, j)); return m_jacobianMatrix(i, j); diff --git a/src/network/lib/engine/procedures/construction.cpp b/src/network/lib/engine/procedures/construction.cpp index 93dfc79b..d891c67e 100644 --- a/src/network/lib/engine/procedures/construction.cpp +++ b/src/network/lib/engine/procedures/construction.cpp @@ -1,3 +1,113 @@ -// -// Created by Emily Boudreaux on 7/14/25. -// +#include "gridfire/engine/procedures/construction.h" + +#include +#include + +#include "gridfire/reaction/reaction.h" +#include "gridfire/reaction/reaclib.h" + +#include "fourdst/composition/composition.h" + +#include "fourdst/logging/logging.h" + +#include "quill/Logger.h" +#include "quill/LogMacros.h" + +namespace gridfire { + using reaction::LogicalReactionSet; + using reaction::ReactionSet; + using reaction::Reaction; + using fourdst::composition::Composition; + using fourdst::atomic::Species; + + + LogicalReactionSet build_reaclib_nuclear_network( + const Composition &composition, + BuildDepthType maxLayers, + bool reverse + ) { + int depth; + if (std::holds_alternative(maxLayers)) { + depth = static_cast(std::get(maxLayers)); + } else { + depth = std::get(maxLayers); + } + auto logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + if (depth == 0) { + LOG_ERROR(logger, "Network build depth is set to 0. No reactions will be collected."); + throw std::logic_error("Network build depth is set to 0. No reactions will be collected."); + } + + const auto allReactions = reaclib::get_all_reactions(); + std::vector remainingReactions; + for (const auto& reaction : allReactions) { + if (reaction.is_reverse() == reverse) { + remainingReactions.push_back(reaction); + } + } + + if (depth == static_cast(NetworkBuildDepth::Full)) { + LOG_INFO(logger, "Building full nuclear network with a total of {} reactions.", allReactions.size()); + const ReactionSet reactionSet(remainingReactions); + return reaction::packReactionSetToLogicalReactionSet(reactionSet); + } + + std::unordered_set availableSpecies; + for (const auto &entry: composition | std::views::values) { + if (entry.mass_fraction() > 0.0) { + availableSpecies.insert(entry.isotope()); + } + } + + + std::vector collectedReactions; + + LOG_INFO(logger, "Starting network construction with {} available species.", availableSpecies.size()); + for (int layer = 0; layer < depth && !remainingReactions.empty(); ++layer) { + LOG_TRACE_L1(logger, "Collecting reactions for layer {} with {} remaining reactions. Currently there are {} available species", layer, remainingReactions.size(), availableSpecies.size()); + std::vector reactionsForNextPass; + std::unordered_set newProductsThisLayer; + bool newReactionsAdded = false; + + reactionsForNextPass.reserve(remainingReactions.size()); + + for (const auto &reaction : remainingReactions) { + bool allReactantsAvailable = true; + for (const auto& reactant : reaction.reactants()) { + if (!availableSpecies.contains(reactant)) { + allReactantsAvailable = false; + break; + } + } + + if (allReactantsAvailable) { + collectedReactions.push_back(reaction); + newReactionsAdded = true; + + for (const auto& product : reaction.products()) { + newProductsThisLayer.insert(product); + } + } else { + reactionsForNextPass.push_back(reaction); + } + } + + if (!newReactionsAdded) { + LOG_INFO(logger, "No new reactions added in layer {}. Stopping network construction with {} reactions collected.", layer, collectedReactions.size()); + break; + } + + LOG_TRACE_L1(logger, "Layer {}: Collected {} reactions. New products this layer: {}", layer, collectedReactions.size(), newProductsThisLayer.size()); + availableSpecies.insert(newProductsThisLayer.begin(), newProductsThisLayer.end()); + + remainingReactions = std::move(reactionsForNextPass); + } + + LOG_INFO(logger, "Network construction completed with {} reactions collected.", collectedReactions.size()); + const ReactionSet reactionSet(collectedReactions); + return reaction::packReactionSetToLogicalReactionSet(reactionSet); + } + + + +} \ No newline at end of file diff --git a/src/network/lib/engine/procedures/priming.cpp b/src/network/lib/engine/procedures/priming.cpp index b1850a30..c57b4f73 100644 --- a/src/network/lib/engine/procedures/priming.cpp +++ b/src/network/lib/engine/procedures/priming.cpp @@ -1,51 +1,192 @@ #include "gridfire/engine/procedures/priming.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 "fourdst/logging/logging.h" +#include "quill/Logger.h" +#include "quill/LogMacros.h" + namespace gridfire { - fourdst::composition::Composition primeNetwork(const NetIn& netIn, DynamicEngine& engine) { - using fourdst::composition::Composition; - using fourdst::atomic::Species; + using fourdst::composition::Composition; + using fourdst::atomic::Species; + + const reaction::Reaction* findDominantCreationChannel ( + const DynamicEngine& engine, + const Species& species, + const std::vector& Y, + const double T9, + const double rho + ) { + const reaction::Reaction* dominateReaction = nullptr; + double maxFlow = -1.0; + for (const auto& reaction : engine.getNetworkReactions()) { + if (reaction.contains(species) && reaction.stoichiometry(species) > 0) { + const double flow = engine.calculateMolarReactionFlow(reaction, Y, T9, rho); + if (flow > maxFlow) { + maxFlow = flow; + dominateReaction = &reaction; + } + } + } + return dominateReaction; + } - NetIn currentConditions = netIn; + PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) { + auto logger = LogManager::getInstance().getLogger("log"); + std::vector speciesToPrime; for (const auto &entry: netIn.composition | std::views::values) { if (entry.mass_fraction() == 0.0) { speciesToPrime.push_back(entry.isotope()); } } + LOG_INFO(logger, "Priming {} species in the network.", speciesToPrime.size()); - + PrimingReport report; if (speciesToPrime.empty()) { - return netIn.composition; // No priming needed + report.primedComposition = netIn.composition; + report.success = true; + report.status = PrimingReportStatus::NO_SPECIES_TO_PRIME; + return report; } - const double T9 = netIn.temperature / 1e9; // Convert temperature to T9 - const double rho = netIn.density; // Density in g/cm^3 + const double T9 = netIn.temperature / 1e9; + const double rho = netIn.density; + const auto initialDepth = engine.getDepth(); + + report.status = PrimingReportStatus::FULL_SUCCESS; + report.success = true; + + // --- 1: pack composition into internal map --- + std::unordered_map currentMassFractions; + for (const auto& entry : netIn.composition | std::views::values) { + currentMassFractions[entry.isotope()] = entry.mass_fraction(); + } + for (const auto& entry : speciesToPrime) { + currentMassFractions[entry] = 0.0; // Initialize priming species with 0 mass fraction + } + + std::unordered_map totalMassFractionChanges; + + engine.rebuild(netIn.composition, NetworkBuildDepth::Full); + + for (const auto& primingSpecies : speciesToPrime) { + LOG_DEBUG(logger, "Priming species: {}", primingSpecies.name()); + + // Create a temporary composition from the current internal state for the primer + Composition tempComp; + for(const auto& [sp, mf] : currentMassFractions) { + tempComp.registerSymbol(std::string(sp.name())); + if (mf < 0.0 && std::abs(mf) < 1e-16) { + tempComp.setMassFraction(sp, 0.0); // Avoid negative mass fractions + } else { + tempComp.setMassFraction(sp, mf); + } + } + tempComp.finalize(true); // Finalize with normalization + + NetIn tempNetIn = netIn; + tempNetIn.composition = tempComp; - for (const auto& primingSpecies :speciesToPrime) { NetworkPrimingEngineView primer(primingSpecies, engine); - const auto Y = primer.mapNetInToMolarAbundanceVector(currentConditions); - const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho); - const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho); - const double equilibriumMolarAbundance = creationRate / destructionRateConstant; - const double equilibriumMassFraction = equilibriumMolarAbundance * primingSpecies.mass(); + 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; + } - currentConditions.composition.setMassFraction( - std::string(primingSpecies.name()), - equilibriumMassFraction - ); + const auto Y = primer.mapNetInToMolarAbundanceVector(tempNetIn); + const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho); + + double equilibriumMassFraction = 0.0; + + if (destructionRateConstant > 1e-99) { + const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho); + equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass(); + LOG_INFO(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction); + + const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho); + + if (dominantChannel) { + LOG_DEBUG(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->peName()); + + double totalReactantMass = 0.0; + for (const auto& reactant : dominantChannel->reactants()) { + totalReactantMass += reactant.mass(); + } + + double scalingFactor = 1.0; + for (const auto& reactant : dominantChannel->reactants()) { + const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass); + double availableMass = 0.0; + if (currentMassFractions.contains(reactant)) { + availableMass = currentMassFractions.at(reactant); + } + if (massToSubtract > availableMass && availableMass > 0) { + scalingFactor = std::min(scalingFactor, availableMass / massToSubtract); + } + } + + if (scalingFactor < 1.0) { + LOG_WARNING(logger, "Priming for {} was limited by reactant availability. Scaling transfer by {:.4e}", primingSpecies.name(), scalingFactor); + equilibriumMassFraction *= scalingFactor; + } + + // Update the internal mass fraction map and accumulate total changes + totalMassFractionChanges[primingSpecies] += equilibriumMassFraction; + currentMassFractions[primingSpecies] += equilibriumMassFraction; + + for (const auto& reactant : dominantChannel->reactants()) { + const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass); + totalMassFractionChanges[reactant] -= massToSubtract; + currentMassFractions[reactant] -= massToSubtract; + } + } else { + LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name()); + report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL; + totalMassFractionChanges[primingSpecies] += 1e-40; + currentMassFractions[primingSpecies] += 1e-40; + } + } else { + LOG_WARNING(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name()); + totalMassFractionChanges[primingSpecies] += 1e-40; + currentMassFractions[primingSpecies] += 1e-40; + report.status = PrimingReportStatus::BASE_NETWORK_TOO_SHALLOW; + } } - return currentConditions.composition; + // --- Final Composition Construction --- + std::vector final_symbols; + std::vector final_mass_fractions; + for(const auto& [species, mass_fraction] : currentMassFractions) { + final_symbols.push_back(std::string(species.name())); + if (mass_fraction < 0.0 && std::abs(mass_fraction) < 1e-16) { + final_mass_fractions.push_back(0.0); // Avoid negative mass fractions + } else { + final_mass_fractions.push_back(mass_fraction); + } + } + + // Create the final composition object from the pre-normalized mass fractions + Composition primedComposition(final_symbols, final_mass_fractions, true); + + report.primedComposition = primedComposition; + for (const auto& [species, change] : totalMassFractionChanges) { + report.massFractionChanges.emplace_back(species, change); + } + + engine.rebuild(netIn.composition, initialDepth); + return report; } - double calculateDestructionRateConstant( + double calculateDestructionRateConstant( const DynamicEngine& engine, const fourdst::atomic::Species& species, const std::vector& Y, diff --git a/src/network/lib/engine/views/engine_adaptive.cpp b/src/network/lib/engine/views/engine_adaptive.cpp index 9cd09ebf..d40f36cb 100644 --- a/src/network/lib/engine/views/engine_adaptive.cpp +++ b/src/network/lib/engine/views/engine_adaptive.cpp @@ -229,6 +229,10 @@ namespace gridfire { return Y; // Return the vector of molar abundances } + PrimingReport AdaptiveEngineView::primeEngine(const NetIn &netIn) { + return m_baseEngine.primeEngine(netIn); + } + int AdaptiveEngineView::getSpeciesIndex(const fourdst::atomic::Species &species) const { auto it = std::find(m_activeSpecies.begin(), m_activeSpecies.end(), species); if (it != m_activeSpecies.end()) { diff --git a/src/network/lib/engine/views/engine_defined.cpp b/src/network/lib/engine/views/engine_defined.cpp index d5a922de..f2818bac 100644 --- a/src/network/lib/engine/views/engine_defined.cpp +++ b/src/network/lib/engine/views/engine_defined.cpp @@ -76,8 +76,6 @@ namespace gridfire { return m_baseEngine; } - - const std::vector & DefinedEngineView::getNetworkSpecies() const { return m_activeSpecies; } @@ -217,6 +215,10 @@ namespace gridfire { return Y; // Return the vector of molar abundances } + PrimingReport DefinedEngineView::primeEngine(const NetIn &netIn) { + return m_baseEngine.primeEngine(netIn); + } + std::vector DefinedEngineView::constructSpeciesIndexMap() const { LOG_TRACE_L1(m_logger, "Constructing species index map for DefinedEngineView..."); std::unordered_map fullSpeciesReverseMap; diff --git a/src/network/lib/engine/views/engine_multiscale.cpp b/src/network/lib/engine/views/engine_multiscale.cpp index 4d8f5509..003dad8a 100644 --- a/src/network/lib/engine/views/engine_multiscale.cpp +++ b/src/network/lib/engine/views/engine_multiscale.cpp @@ -6,6 +6,7 @@ #include namespace gridfire { + static int s_operator_parens_called = 0; using fourdst::atomic::Species; MultiscalePartitioningEngineView::MultiscalePartitioningEngineView( @@ -58,8 +59,8 @@ namespace gridfire { double MultiscalePartitioningEngineView::calculateMolarReactionFlow( const reaction::Reaction &reaction, const std::vector &Y, - double T9, - double rho + const double T9, + const double rho ) const { return m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho); } @@ -230,18 +231,45 @@ namespace gridfire { return m_dynamic_species; } - void MultiscalePartitioningEngineView::equilibrateNetwork( + PrimingReport MultiscalePartitioningEngineView::primeEngine(const NetIn &netIn) { + return m_baseEngine.primeEngine(netIn); + } + + fourdst::composition::Composition MultiscalePartitioningEngineView::equilibrateNetwork( const std::vector &Y, const double T9, const double rho, const double dt_control ) { partitionNetwork(Y, T9, rho, dt_control); - std::vector Y_equilibrated = solveQSEAbundances(Y, T9, rho); + const std::vector Y_equilibrated = solveQSEAbundances(Y, T9, rho); + fourdst::composition::Composition composition; + + std::vector symbols; + symbols.reserve(m_baseEngine.getNetworkSpecies().size()); + for (const auto& species : m_baseEngine.getNetworkSpecies()) { + symbols.emplace_back(species.name()); + } + composition.registerSymbol(symbols); + + std::vector X; + X.reserve(Y_equilibrated.size()); + for (size_t i = 0; i < Y_equilibrated.size(); ++i) { + const double molarMass = m_baseEngine.getNetworkSpecies()[i].mass(); + X.push_back(Y_equilibrated[i] * molarMass); // Convert from molar abundance to mass fraction + } + + for (size_t i = 0; i < Y_equilibrated.size(); ++i) { + const auto& species = m_baseEngine.getNetworkSpecies()[i]; + composition.setMassFraction(std::string(species.name()), X[i]); + } + + composition.finalize(true); + return composition; } - void MultiscalePartitioningEngineView::equilibrateNetwork( - const NetIn& netIn, + fourdst::composition::Composition MultiscalePartitioningEngineView::equilibrateNetwork( + const NetIn &netIn, const double dt_control ) { std::vector Y(m_baseEngine.getNetworkSpecies().size(), 0.0); @@ -257,10 +285,7 @@ namespace gridfire { const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) const double rho = netIn.density; // Density in g/cm^3 - partitionNetwork(Y, T9, rho, dt_control); - std::vector Y_equilibrated = solveQSEAbundances(Y, T9, rho); - - + return equilibrateNetwork(Y, T9, rho, dt_control); } int MultiscalePartitioningEngineView::getSpeciesIndex(const fourdst::atomic::Species &species) const { @@ -566,14 +591,13 @@ namespace gridfire { Eigen::VectorXd Y_final_qse = Y_scale.array() * v_initial.array().sinh(); // Convert back to physical abundances using asinh scaling for (size_t i = 0; i < qse_solve_indices.size(); ++i) { Y_out[qse_solve_indices[i]] = Y_final_qse(i); - std::cout << "Updating species " << m_baseEngine.getNetworkSpecies()[qse_solve_indices[i]].name() - << " to QSE abundance: " << Y_final_qse(i) << std::endl; } } return Y_out; } int MultiscalePartitioningEngineView::EigenFunctor::operator()(const InputType &v_qse, OutputType &f_qse) const { + s_operator_parens_called++; std::vector y_trial = m_Y_full_initial; Eigen::VectorXd y_qse = m_Y_scale.array() * v_qse.array().sinh(); // Convert to physical abundances using asinh scaling @@ -586,6 +610,7 @@ namespace gridfire { for (size_t i = 0; i < m_qse_solve_indices.size(); ++i) { f_qse(i) = dydt[m_qse_solve_indices[i]]; } + return 0; // Success } diff --git a/src/network/lib/engine/views/engine_priming.cpp b/src/network/lib/engine/views/engine_priming.cpp index 1abc3a65..7132eff8 100644 --- a/src/network/lib/engine/views/engine_priming.cpp +++ b/src/network/lib/engine/views/engine_priming.cpp @@ -44,7 +44,8 @@ namespace gridfire { ), baseEngine ), - m_primingSpecies(primingSpecies) {} + m_primingSpecies(primingSpecies) { + } std::vector NetworkPrimingEngineView::constructPrimingReactionSet( diff --git a/src/network/lib/network.cpp b/src/network/lib/network.cpp index 1c161069..850a3043 100644 --- a/src/network/lib/network.cpp +++ b/src/network/lib/network.cpp @@ -61,32 +61,6 @@ namespace gridfire { return oldFormat; } - reaction::LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, bool reverse) { - using namespace reaction; - std::vector reaclibReactions; - auto logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - - for (const auto &reaction: reaclib::get_all_reactions()) { - if (reaction.is_reverse() != reverse) { - continue; // Skip reactions that do not match the requested direction - } - bool gotReaction = true; - const auto& reactants = reaction.reactants(); - for (const auto& reactant : reactants) { - if (!composition.contains(reactant)) { - gotReaction = false; - break; // If any reactant is not in the composition, skip this reaction - } - } - if (gotReaction) { - LOG_TRACE_L3(logger, "Adding reaction {} to REACLIB reaction set.", reaction.peName()); - reaclibReactions.push_back(reaction); - } - } - const ReactionSet reactionSet(reaclibReactions); - return packReactionSetToLogicalReactionSet(reactionSet); - } - // Trim whitespace from both ends of a string std::string trim_whitespace(const std::string& str) { auto startIt = str.begin(); diff --git a/src/network/lib/solver/solver.cpp b/src/network/lib/solver/solver.cpp index 65701d7c..52e1ca46 100644 --- a/src/network/lib/solver/solver.cpp +++ b/src/network/lib/solver/solver.cpp @@ -4,7 +4,6 @@ #include "gridfire/utils/logging.h" -#include "gridfire/utils/qse_rules.h" #include "fourdst/composition/atomicSpecies.h" #include "fourdst/composition/composition.h" @@ -176,11 +175,11 @@ namespace gridfire::solver { const double final_timescale = std::min(network_timescale, decay_timescale); - const bool isQSE = is_species_in_qse( - network_timescale, - decay_timescale, - abundance - ); + const bool isQSE = false; + // network_timescale, + // decay_timescale, + // abundance + // ); if (isQSE) { LOG_TRACE_L2(m_logger, "{} is in QSE based on rules in qse_rules.h", species.name()); diff --git a/src/network/meson.build b/src/network/meson.build index f31e8a21..7016d081 100644 --- a/src/network/meson.build +++ b/src/network/meson.build @@ -8,6 +8,7 @@ network_sources = files( 'lib/engine/views/engine_multiscale.cpp', 'lib/engine/views/engine_priming.cpp', 'lib/engine/procedures/priming.cpp', + 'lib/engine/procedures/construction.cpp', 'lib/reaction/reaction.cpp', 'lib/reaction/reaclib.cpp', 'lib/io/network_file.cpp', @@ -59,6 +60,7 @@ network_headers = files( 'include/gridfire/engine/views/engine_multiscale.h', 'include/gridfire/engine/views/engine_priming.h', 'include/gridfire/engine/procedures/priming.h', + 'include/gridfire/engine/procedures/construction.h', 'include/gridfire/reaction/reaction.h', 'include/gridfire/reaction/reaclib.h', 'include/gridfire/io/network_file.h', @@ -72,6 +74,5 @@ network_headers = files( 'include/gridfire/partition/partition_ground.h', 'include/gridfire/partition/composite/partition_composite.h', 'include/gridfire/utils/logging.h', - 'include/gridfire/utils/qse_rules.h', ) install_headers(network_headers, subdir : 'gridfire')