From 7012eb819ab34d3367b62e1502be43deb9c4b8b4 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Thu, 10 Jul 2025 09:36:05 -0400 Subject: [PATCH] feat(MultiscalePartitioningEngineView): added *much* more robust qse group identifiction and solving --- .../include/gridfire/engine/engine_abstract.h | 8 +- .../include/gridfire/engine/engine_graph.h | 258 ++++++-- .../gridfire/engine/procedures/priming.h | 31 + .../gridfire/engine/views/engine_adaptive.h | 8 +- .../gridfire/engine/views/engine_defined.h | 137 +--- .../gridfire/engine/views/engine_multiscale.h | 188 ++++++ .../gridfire/engine/views/engine_priming.h | 38 ++ .../include/gridfire/io/network_file.h | 18 +- .../include/gridfire/reaction/reaction.h | 5 + .../gridfire/screening/screening_weak.h | 2 +- src/network/include/gridfire/solver/solver.h | 65 +- src/network/lib/engine/engine_graph.cpp | 224 ++++++- src/network/lib/engine/procedures/priming.cpp | 84 +++ .../lib/engine/views/engine_adaptive.cpp | 23 +- .../lib/engine/views/engine_defined.cpp | 238 ++++--- .../lib/engine/views/engine_multiscale.cpp | 623 ++++++++++++++++++ .../lib/engine/views/engine_priming.cpp | 72 ++ src/network/lib/io/network_file.cpp | 4 +- .../composite/partition_composite.cpp | 8 +- src/network/lib/solver/solver.cpp | 150 ++++- src/network/meson.build | 7 + subprojects/libcomposition.wrap | 2 +- tests/graphnet_sandbox/main.cpp | 48 +- 23 files changed, 1863 insertions(+), 378 deletions(-) create mode 100644 src/network/include/gridfire/engine/procedures/priming.h create mode 100644 src/network/include/gridfire/engine/views/engine_multiscale.h create mode 100644 src/network/include/gridfire/engine/views/engine_priming.h create mode 100644 src/network/lib/engine/procedures/priming.cpp create mode 100644 src/network/lib/engine/views/engine_multiscale.cpp create mode 100644 src/network/lib/engine/views/engine_priming.cpp diff --git a/src/network/include/gridfire/engine/engine_abstract.h b/src/network/include/gridfire/engine/engine_abstract.h index 69abdafe..b7aa6a9d 100644 --- a/src/network/include/gridfire/engine/engine_abstract.h +++ b/src/network/include/gridfire/engine/engine_abstract.h @@ -123,7 +123,7 @@ namespace gridfire { /** * @brief Generate the Jacobian matrix for the current state. * - * @param Y Vector of current abundances. + * @param Y_dynamic Vector of current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @@ -131,7 +131,7 @@ namespace gridfire { * for the current state. The matrix can then be accessed via getJacobianMatrixEntry(). */ virtual void generateJacobianMatrix( - const std::vector& Y, + const std::vector& Y_dynamic, double T9, double rho ) = 0; @@ -265,5 +265,9 @@ namespace gridfire { * @endcode */ [[nodiscard]] virtual screening::ScreeningType getScreeningModel() const = 0; + + [[nodiscard]] virtual int getSpeciesIndex(const fourdst::atomic::Species &species) const = 0; + + [[nodiscard]] virtual std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const = 0; }; } \ 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 8eef8cc1..100c49c9 100644 --- a/src/network/include/gridfire/engine/engine_graph.h +++ b/src/network/include/gridfire/engine/engine_graph.h @@ -17,10 +17,12 @@ #include #include #include +#include #include #include "cppad/cppad.hpp" +#include "quill/LogMacros.h" // PERF: The function getNetReactionStoichiometry returns a map of species to their stoichiometric coefficients for a given reaction. // this makes extra copies of the species, which is not ideal and could be optimized further. @@ -139,7 +141,7 @@ namespace gridfire { /** * @brief Generates the Jacobian matrix for the current state. * - * @param Y Vector of current abundances. + * @param Y_dynamic Vector of current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. * @@ -150,7 +152,7 @@ namespace gridfire { * @see getJacobianMatrixEntry() */ void generateJacobianMatrix( - const std::vector& Y, + const std::vector& Y_dynamic, const double T9, const double rho ) override; @@ -317,18 +319,48 @@ namespace gridfire { [[nodiscard]] double calculateReverseRate( const reaction::Reaction &reaction, - double T9, - double expFactor + double T9 ) const; + double calculateReverseRateTwoBody( + const reaction::Reaction &reaction, + const double T9, + const double forwardRate, + const double expFactor + ) const; + + double calculateReverseRateTwoBodyDerivative( + const reaction::Reaction &reaction, + const double T9, + const double reverseRate + ) const; + + bool isUsingReverseReactions() const; + + void setUseReverseReactions(bool useReverse); + + int getSpeciesIndex( + const fourdst::atomic::Species& species + ) const override; + + std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const override; + + + private: struct PrecomputedReaction { + // Forward cacheing size_t reaction_index; std::vector unique_reactant_indices; std::vector reactant_powers; double symmetry_factor; std::vector affected_species_indices; std::vector stoichiometric_coefficients; + + // Reverse cacheing + std::vector unique_product_indices; ///< Unique product indices for reverse reactions. + std::vector product_powers; ///< Powers of each unique product in the reverse reaction. + double reverse_symmetry_factor; ///< Symmetry factor for reverse reactions. }; struct constants { @@ -337,7 +369,47 @@ namespace gridfire { const double c = Constants::getInstance().get("c").value; ///< Speed of light in cm/s. const double kB = Constants::getInstance().get("kB").value; ///< Boltzmann constant in erg/K. }; + private: + class AtomicReverseRate final : public CppAD::atomic_base { + public: + AtomicReverseRate( + const reaction::Reaction& reaction, + const GraphEngine& engine + ): + atomic_base("AtomicReverseRate"), + m_reaction(reaction), + m_engine(engine) {} + bool forward( + size_t p, + size_t q, + const CppAD::vector& vx, + CppAD::vector& vy, + const CppAD::vector& tx, + CppAD::vector& ty + ) override; + bool reverse( + size_t q, + const CppAD::vector& tx, + const CppAD::vector& ty, + CppAD::vector& px, + const CppAD::vector& py + ) override; + bool for_sparse_jac( + size_t q, + const CppAD::vector>&r, + CppAD::vector>& s + ) override; + bool rev_sparse_jac( + size_t q, + const CppAD::vector>&rt, + CppAD::vector>& st + ) override; + + private: + const reaction::Reaction& m_reaction; + const GraphEngine& m_engine; + }; private: Config& m_config = Config::getInstance(); quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); @@ -355,12 +427,15 @@ 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. + std::vector> m_atomicReverseRates; screening::ScreeningType m_screeningType = screening::ScreeningType::BARE; ///< Screening type for the reaction network. Default to no screening. std::unique_ptr m_screeningModel = screening::selectScreeningModel(m_screeningType); bool m_usePrecomputation = true; ///< Flag to enable or disable using precomputed reactions for efficiency. Mathematically, this should not change the results. Generally end users should not need to change this. + bool m_useReverseReactions = true; ///< Flag to enable or disable reverse reactions. If false, only forward reactions are considered. + std::vector m_precomputedReactions; ///< Precomputed reactions for efficiency. std::unique_ptr m_partitionFunction; ///< Partition function for the network. @@ -418,8 +493,11 @@ namespace gridfire { */ void recordADTape(); + void collectAtomicReverseRateAtomicBases(); + void precomputeNetwork(); + /** * @brief Validates mass and charge conservation across all reactions. * @@ -449,24 +527,12 @@ namespace gridfire { ); - double calculateReverseRateTwoBody( - const reaction::Reaction &reaction, - const double T9, - const double forwardRate, - const double expFactor - ) const; - - double GraphEngine::calculateReverseRateTwoBodyDerivative( - const reaction::Reaction &reaction, - const double T9, - const double reverseRate - ) const; [[nodiscard]] StepDerivatives calculateAllDerivativesUsingPrecomputation( const std::vector &Y_in, const std::vector& bare_rates, - double T9, - double rho + const std::vector &bare_reverse_rates, + double T9, double rho ) const; /** @@ -490,6 +556,16 @@ namespace gridfire { const T rho ) const; + template + T calculateReverseMolarReactionFlow( + T T9, + T rho, + std::vector screeningFactors, + std::vector Y, + size_t reactionIndex, + const reaction::LogicalReaction &reaction + ) const; + /** * @brief Calculates all derivatives (dY/dt) and the energy generation rate. * @@ -547,6 +623,74 @@ namespace gridfire { }; + + template + T GraphEngine::calculateReverseMolarReactionFlow( + T T9, + T rho, + std::vector screeningFactors, + std::vector Y, + size_t reactionIndex, + const reaction::LogicalReaction &reaction + ) const { + if (!m_useReverseReactions) { + return static_cast(0.0); // If reverse reactions are not used, return zero + } + T reverseMolarFlow = static_cast(0.0); + + if (reaction.qValue() != 0.0) { + T reverseRateConstant = static_cast(0.0); + if constexpr (std::is_same_v) { // Check if T is an AD type at compile time + const auto& atomic_func_ptr = m_atomicReverseRates[reactionIndex]; + if (atomic_func_ptr != nullptr) { + // A. Instantiate the atomic operator for the specific reaction + // B. Marshal the input vector + std::vector ax = { T9 }; + + std::vector ay(1); + (*atomic_func_ptr)(ax, ay); + reverseRateConstant = static_cast(ay[0]); + } else { + return reverseMolarFlow; // If no atomic function is available, return zero + } + } else { + // A,B If not calling with an AD type, calculate the reverse rate directly + reverseRateConstant = calculateReverseRate(reaction, T9); + } + + // C. Get product multiplicities + std::unordered_map productCounts; + for (const auto& product : reaction.products()) { + productCounts[product]++; + } + + // D. Calculate the symmetry factor + T reverseSymmetryFactor = static_cast(1.0); + for (const auto &count: productCounts | std::views::values) { + reverseSymmetryFactor /= static_cast(std::tgamma(static_cast(count + 1))); // Gamma function for factorial + } + + // E. Calculate the abundance term + T productAbundanceTerm = static_cast(1.0); + for (const auto& [species, count] : productCounts) { + const int speciesIndex = m_speciesToIndexMap.at(species); + productAbundanceTerm *= CppAD::pow(Y[speciesIndex], count); + } + + // F. Determine the power for the density term + const size_t num_products = reaction.products().size(); + const T rho_power = CppAD::pow(rho, static_cast(num_products > 1 ? num_products - 1 : 0)); // Density raised to the power of (N-1) for N products + + // G. Assemble the reverse molar flow rate + reverseMolarFlow = screeningFactors[reactionIndex] * + reverseRateConstant * + productAbundanceTerm * + reverseSymmetryFactor * + rho_power; + } + return reverseMolarFlow; + } + template StepDerivatives GraphEngine::calculateAllDerivatives( const std::vector &Y_in, T T9, T rho) const { @@ -594,13 +738,39 @@ namespace gridfire { for (size_t reactionIndex = 0; reactionIndex < m_reactions.size(); ++reactionIndex) { const auto& reaction = m_reactions[reactionIndex]; - // 1. Calculate reaction rate - const T molarReactionFlow = screeningFactors[reactionIndex] * calculateMolarReactionFlow(reaction, Y, T9, rho); + // 1. Calculate forward reaction rate + const T forwardMolarReactionFlow = screeningFactors[reactionIndex] * + calculateMolarReactionFlow(reaction, Y, T9, rho); - // 2. Use the rate to update all relevant species derivatives (dY/dt) + // 2. Calculate reverse reaction rate + T reverseMolarFlow = calculateReverseMolarReactionFlow( + T9, + rho, + screeningFactors, + Y, + reactionIndex, + reaction + ); + + const T molarReactionFlow = forwardMolarReactionFlow - reverseMolarFlow; // Net molar reaction flow + std::stringstream ss; + ss << "Forward: " << forwardMolarReactionFlow + << ", Reverse: " << reverseMolarFlow + << ", Net: " << molarReactionFlow; + LOG_DEBUG( + m_logger, + "Reaction: {}, {}", + reaction.peName(), + ss.str() + ); + // std::cout << "Forward molar flow for reaction " << reaction.peName() << ": " << forwardMolarReactionFlow << std::endl; + // std::cout << "Reverse molar flow for reaction " << reaction.peName() << ": " << reverseMolarFlow << std::endl; + // std::cout << "Net molar flow for reaction " << reaction.peName() << ": " << molarReactionFlow << std::endl; + + // 3. Use the rate to update all relevant species derivatives (dY/dt) for (size_t speciesIndex = 0; speciesIndex < m_networkSpecies.size(); ++speciesIndex) { const T nu_ij = static_cast(m_stoichiometryMatrix(speciesIndex, reactionIndex)); - result.dydt[speciesIndex] += threshold_flag * nu_ij * molarReactionFlow / rho; + result.dydt[speciesIndex] += threshold_flag * nu_ij * molarReactionFlow; } } @@ -644,6 +814,7 @@ namespace gridfire { for (const auto& reactant : reaction.reactants()) { reactant_counts[std::string(reactant.name())]++; } + const int totalReactants = static_cast(reaction.reactants().size()); // --- Accumulator for the molar concentration --- auto molar_concentration_product = static_cast(1.0); @@ -657,56 +828,23 @@ namespace gridfire { const T Yi = Y[species_index]; // --- Check if the species abundance is below the threshold where we ignore reactions --- - threshold_flag *= CppAD::CondExpLt(Yi, Y_threshold, zero, one); - - // --- Convert from molar abundance to molar concentration --- - T molar_concentration = Yi * rho; + // threshold_flag *= CppAD::CondExpLt(Yi, Y_threshold, zero, one); // --- If count is > 1 , we need to raise the molar concentration to the power of count since there are really count bodies in that reaction --- - molar_concentration_product *= CppAD::pow(molar_concentration, static_cast(count)); // ni^count + molar_concentration_product *= CppAD::pow(Yi, static_cast(count)); // ni^count // --- Apply factorial correction for identical reactions --- if (count > 1) { molar_concentration_product /= static_cast(std::tgamma(static_cast(count + 1))); // Gamma function for factorial } } - // --- Final reaction flow calculation [mol][s^-1][cm^-3] --- + // --- Final reaction flow calculation [mol][s^-1][g^-1] --- // Note: If the threshold flag ever gets set to zero this will return zero. // This will result basically in multiple branches being written to the AD tape, which will make // the tape more expensive to record, but it will also mean that we only need to record it once for // the entire network. - return molar_concentration_product * k_reaction * threshold_flag; + const T densityTerm = CppAD::pow(rho, totalReactants > 1 ? static_cast(totalReactants - 1) : zero); // Density raised to the power of (N-1) for N reactants + return molar_concentration_product * k_reaction * threshold_flag * densityTerm; } - class AtomicReverseRate final : public CppAD::atomic_base { - public: - AtomicReverseRate( - const reaction::Reaction& reaction, - const GraphEngine& engine - ): - atomic_base("AtomicReverseRate"), - m_reaction(reaction), - m_engine(engine) {} - - bool forward( - size_t p, - size_t q, - const CppAD::vector& vx, - CppAD::vector& vy, - const CppAD::vector& tx, - CppAD::vector& ty - ) override; - bool reverse( - size_t id, - size_t an, - const CppAD::vector& tx, - const CppAD::vector& ty, - CppAD::vector& px, - const CppAD::vector& py - ); - private: - const double m_kB = Constants::getInstance().get("k_b").value; ///< Boltzmann constant in erg/K. - const reaction::Reaction& m_reaction; - const GraphEngine& m_engine; - }; }; \ 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 new file mode 100644 index 00000000..d3cc91e1 --- /dev/null +++ b/src/network/include/gridfire/engine/procedures/priming.h @@ -0,0 +1,31 @@ +#pragma once + +#include "gridfire/engine/engine_abstract.h" +#include "gridfire/network.h" + +#include "fourdst/composition/composition.h" +#include "fourdst/composition/atomicSpecies.h" + + +namespace gridfire { + fourdst::composition::Composition primeNetwork( + const NetIn&, + DynamicEngine& engine + ); + + double calculateDestructionRateConstant( + const DynamicEngine& engine, + const fourdst::atomic::Species& species, + const std::vector& Y, + double T9, + double rho + ); + + double calculateCreationRate( + const DynamicEngine& engine, + const fourdst::atomic::Species& species, + const std::vector& Y, + double T9, + double rho + ); +} \ 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 da0d76ed..510bc6c3 100644 --- a/src/network/include/gridfire/engine/views/engine_adaptive.h +++ b/src/network/include/gridfire/engine/views/engine_adaptive.h @@ -106,7 +106,7 @@ namespace gridfire { /** * @brief Generates the Jacobian matrix for the active species. * - * @param Y_culled A vector of abundances for the active species. + * @param Y_dynamic A vector of abundances for the active species. * @param T9 The temperature in units of 10^9 K. * @param rho The density in g/cm^3. * @@ -117,7 +117,7 @@ namespace gridfire { * @see AdaptiveEngineView::update() */ void generateJacobianMatrix( - const std::vector &Y_culled, + const std::vector &Y_dynamic, const double T9, const double rho ) override; @@ -256,6 +256,10 @@ namespace gridfire { * @endcode */ [[nodiscard]] screening::ScreeningType getScreeningModel() const override; + + [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override; + + [[nodiscard]] std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const 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 5ec910f4..454199b5 100644 --- a/src/network/include/gridfire/engine/views/engine_defined.h +++ b/src/network/include/gridfire/engine/views/engine_defined.h @@ -13,56 +13,9 @@ #include namespace gridfire{ - /** - * @class FileDefinedEngineView - * @brief An engine view that uses a user-defined reaction network from a file. - * - * This class implements an EngineView that restricts the reaction network to a specific set of - * reactions defined in an external file. It acts as a filter or a view on a larger, more - * comprehensive base engine. The file provides a list of reaction identifiers, and this view - * will only consider those reactions and the species involved in them. - * - * This is useful for focusing on a specific sub-network for analysis, debugging, or performance - * reasons, without modifying the underlying full network. - * - * The view maintains mappings between the indices of its active (defined) species and reactions - * and the corresponding indices in the full network of the base engine. All calculations are - * delegated to the base engine after mapping the inputs from the view's context to the full - * network context, and the results are mapped back. - * - * @implements DynamicEngine - * @implements EngineView - */ - class FileDefinedEngineView final: public DynamicEngine, public EngineView { + class DefinedEngineView : public DynamicEngine, public EngineView { public: - /** - * @brief Constructs a FileDefinedEngineView. - * - * @param baseEngine The underlying DynamicEngine to which this view delegates calculations. - * @param fileName The path to the file that defines the reaction network for this view. - * @param parser A reference to a parser object capable of parsing the network file. - * - * @par Usage Example: - * @code - * MyParser parser; - * DynamicEngine baseEngine(...); - * FileDefinedEngineView view(baseEngine, "my_network.net", parser); - * @endcode - * - * @post The view is initialized with the reactions and species from the specified file. - * @throws std::runtime_error If a reaction from the file is not found in the base engine. - */ - explicit FileDefinedEngineView( - DynamicEngine& baseEngine, - const std::string& fileName, - const io::NetworkFileParser& parser - ); - - // --- EngineView Interface --- - /** - * @brief Gets the base engine. - * @return A const reference to the base engine. - */ + DefinedEngineView(const std::vector& peNames, DynamicEngine& baseEngine); const DynamicEngine& getBaseEngine() const override; // --- Engine Interface --- @@ -92,14 +45,14 @@ namespace gridfire{ /** * @brief Generates the Jacobian matrix for the active species. * - * @param Y_defined A vector of abundances for the active species. + * @param Y_dynamic A vector of abundances for the active species. * @param T9 The temperature in units of 10^9 K. * @param rho The density in g/cm^3. * * @throws std::runtime_error If the view is stale. */ void generateJacobianMatrix( - const std::vector& Y_defined, + const std::vector& Y_dynamic, const double T9, const double rho ) override; @@ -191,21 +144,6 @@ namespace gridfire{ */ void update(const NetIn &netIn) override; - /** - * @brief Sets a new network file to define the active reactions. - * - * @param fileName The path to the new network definition file. - * - * @par Usage Example: - * @code - * view.setNetworkFile("another_network.net"); - * view.update(netIn); // Must be called before using the view again - * @endcode - * - * @post The view is marked as stale. `update()` must be called before further use. - */ - void setNetworkFile(const std::string& fileName); - /** * @brief Sets the screening model for the base engine. * @@ -219,21 +157,15 @@ namespace gridfire{ * @return The current screening model type. */ [[nodiscard]] screening::ScreeningType getScreeningModel() const override; - private: - using Config = fourdst::config::Config; - using LogManager = fourdst::logging::LogManager; - /** @brief A reference to the singleton Config instance. */ - Config& m_config = Config::getInstance(); - /** @brief A pointer to the logger instance. */ - quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); - /** @brief The underlying engine to which this view delegates calculations. */ + [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override; + + [[nodiscard]] std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const override; + protected: + bool m_isStale = true; DynamicEngine& m_baseEngine; - ///< Name of the file defining the reaction set considered by the engine view. - std::string m_fileName; - ///< Parser for the network file. - const io::NetworkFileParser& m_parser; - + private: + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance for trace and debug information. ///< Active species in the defined engine. std::vector m_activeSpecies; ///< Active reactions in the defined engine. @@ -243,30 +175,7 @@ namespace gridfire{ std::vector m_speciesIndexMap; ///< Maps indices of active reactions to indices in the full network. std::vector m_reactionIndexMap; - - /** @brief A flag indicating whether the view is stale and needs to be updated. */ - bool m_isStale = true; - private: - /** - * @brief Builds the active species and reaction sets from a file. - * - * This method uses the provided parser to read reaction names from the given file. - * It then finds these reactions in the base engine's full network and populates - * the `m_activeReactions` and `m_activeSpecies` members. Finally, it constructs - * the index maps for the active sets. - * - * @param fileName The path to the network definition file. - * - * @post - * - `m_activeReactions` and `m_activeSpecies` are populated. - * - `m_speciesIndexMap` and `m_reactionIndexMap` are constructed. - * - `m_isStale` is set to false. - * - * @throws std::runtime_error If a reaction from the file is not found in the base engine. - */ - void buildFromFile(const std::string& fileName); - /** * @brief Constructs the species index map. * @@ -329,11 +238,25 @@ namespace gridfire{ */ size_t mapViewToFullReactionIndex(size_t definedReactionIndex) const; - /** - * @brief Validates that the FileDefinedEngineView is not stale. - * - * @throws std::runtime_error If the view is stale (i.e., `update()` has not been called after the view was made stale). - */ void validateNetworkState() const; }; + + class FileDefinedEngineView final: public DefinedEngineView { + public: + explicit FileDefinedEngineView( + DynamicEngine& baseEngine, + const std::string& fileName, + const io::NetworkFileParser& parser + ); + std::string getNetworkFile() const { return m_fileName; } + const io::NetworkFileParser& getParser() const { return m_parser; } + private: + using Config = fourdst::config::Config; + using LogManager = fourdst::logging::LogManager; + Config& m_config = Config::getInstance(); + quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); + std::string m_fileName; + ///< Parser for the network file. + const io::NetworkFileParser& m_parser; + }; } \ No newline at end of file diff --git a/src/network/include/gridfire/engine/views/engine_multiscale.h b/src/network/include/gridfire/engine/views/engine_multiscale.h new file mode 100644 index 00000000..bf66980c --- /dev/null +++ b/src/network/include/gridfire/engine/views/engine_multiscale.h @@ -0,0 +1,188 @@ +#pragma once + +#include "gridfire/engine/engine_abstract.h" +#include "gridfire/engine/views/engine_view_abstract.h" +#include "gridfire/engine/engine_graph.h" + +#include "Eigen/Dense" +#include "unsupported/Eigen/NonLinearOptimization" + +namespace gridfire { + class MultiscalePartitioningEngineView final: public DynamicEngine, public EngineView { + public: + explicit MultiscalePartitioningEngineView(GraphEngine& baseEngine); + + [[nodiscard]] const std::vector & getNetworkSpecies() const override; + + [[nodiscard]] StepDerivatives calculateRHSAndEnergy( + const std::vector &Y, + double T9, + double rho + ) const override; + + void generateJacobianMatrix( + const std::vector &Y_dynamic, + double T9, + double rho + ) override; + + [[nodiscard]] double getJacobianMatrixEntry( + int i, + int j + ) const override; + + void generateStoichiometryMatrix() override; + + [[nodiscard]] int getStoichiometryMatrixEntry( + int speciesIndex, + int reactionIndex + ) const override; + + [[nodiscard]] double calculateMolarReactionFlow( + const reaction::Reaction &reaction, + const std::vector &Y, + double T9, + double rho + ) const override; + + [[nodiscard]] const reaction::LogicalReactionSet & getNetworkReactions() const override; + + [[nodiscard]] std::unordered_map getSpeciesTimescales( + const std::vector &Y, + double T9, + double rho + ) const override; + + void update( + const NetIn &netIn + ) override; + + void setScreeningModel( + screening::ScreeningType model + ) override; + + [[nodiscard]] screening::ScreeningType getScreeningModel() const override; + + const DynamicEngine & getBaseEngine() const override; + + void partitionNetwork( + const std::vector& Y, + double T9, + double rho, + double dt_control + ); + + void partitionNetwork( + const NetIn& netIn, + double dt_control + ); + + void exportToDot(const std::string& filename) const; + + [[nodiscard]] int getSpeciesIndex(const fourdst::atomic::Species &species) const override; + + [[nodiscard]] std::vector mapNetInToMolarAbundanceVector(const NetIn &netIn) const override; + + std::vector getFastSpecies() const; + const std::vector& getDynamicSpecies() const; + + void equilibrateNetwork( + const std::vector& Y, + double T9, + double rho, + double dt_control + ); + + void equilibrateNetwork( + const NetIn& netIn, + const double dt_control + ); + + + private: + struct QSEGroup { + std::vector species_indices; ///< Indices of all species in this group. + size_t seed_nucleus_index; ///< Index of the one species that will be evolved dynamically. + bool is_in_equilibrium = false; ///< Flag set by flux analysis. + }; + + struct EigenFunctor { + using InputType = Eigen::Matrix; + using OutputType = Eigen::Matrix; + using JacobianType = Eigen::Matrix; + enum { + InputsAtCompileTime = Eigen::Dynamic, + ValuesAtCompileTime = Eigen::Dynamic + }; + + MultiscalePartitioningEngineView* m_view; + const std::vector& m_qse_solve_indices; + const std::vector& m_Y_full_initial; + const double m_T9; + const double m_rho; + const Eigen::VectorXd& m_Y_scale; + + EigenFunctor( + MultiscalePartitioningEngineView& view, + const std::vector& qse_solve_indices, + const std::vector& Y_full_initial, + const double T9, + const double rho, + const Eigen::VectorXd& Y_scale + ) : + m_view(&view), + m_qse_solve_indices(qse_solve_indices), + m_Y_full_initial(Y_full_initial), + m_T9(T9), + 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(); } + + int operator()(const InputType& v_qse, OutputType& f_qse) const; + int df(const InputType& v_qse, JacobianType& J_qse) const; + }; + + private: + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + GraphEngine& m_baseEngine; ///< The base engine to which this view delegates calculations. + std::vector m_qse_groups; ///< The list of identified equilibrium groups. + std::vector m_dynamic_species; ///< The simplified set of species presented to the solver. + std::vector m_dynamic_species_indices; ///< Indices mapping the dynamic species back to the master engine's list. + std::unordered_map> m_connectivity_graph; + + private: + std::unordered_set identifyFastReactions( + const std::vector& Y_full, + double T9, + double rho, + double dt_control + ) const; + + void buildConnectivityGraph( + const std::unordered_set& fast_reaction_indices + ); + + void findConnectedComponents(); + + void validateGroupsWithFluxAnalysis( + const std::vector& Y, + double T9, + double rho + ); + + std::pair, std::vector> identifyDynamicSpecies( + const std::vector& Y, + const std::vector& qse_groups, + double T9, + double rho + ) const; + + std::vector solveQSEAbundances( + const std::vector &Y_full, + double T9, + double rho + ); + }; +} diff --git a/src/network/include/gridfire/engine/views/engine_priming.h b/src/network/include/gridfire/engine/views/engine_priming.h new file mode 100644 index 00000000..d38471df --- /dev/null +++ b/src/network/include/gridfire/engine/views/engine_priming.h @@ -0,0 +1,38 @@ +#pragma once + +#include "gridfire/engine/engine_abstract.h" +#include "gridfire/engine/views/engine_defined.h" + +#include "gridfire/network.h" + +#include "fourdst/logging/logging.h" +#include "fourdst/composition/atomicSpecies.h" +#include "fourdst/composition/composition.h" + +#include "quill/Logger.h" + +#include +#include + +namespace gridfire { + + class NetworkPrimingEngineView final : public DefinedEngineView { + 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( + const fourdst::atomic::Species& primingSpecies, + const DynamicEngine& baseEngine + ) const; + }; + +} \ No newline at end of file diff --git a/src/network/include/gridfire/io/network_file.h b/src/network/include/gridfire/io/network_file.h index d742ec7b..c9dd276d 100644 --- a/src/network/include/gridfire/io/network_file.h +++ b/src/network/include/gridfire/io/network_file.h @@ -9,23 +9,7 @@ #include namespace gridfire::io { - - /** - * @struct ParsedNetworkData - * @brief Holds the data parsed from a network file. - * - * This struct is used to return the results of parsing a reaction network - * file. It contains the list of reaction names that define the network. - */ - struct ParsedNetworkData { - /** - * @brief A vector of reaction names in their PEN-style format. - * - * Projectile, Ejectile style names p(p,e+)d is a common format for representing - * nuclear reactions as strings. - */ - std::vector reactionPENames; - }; + typedef std::vector ParsedNetworkData; /** * @class NetworkFileParser diff --git a/src/network/include/gridfire/reaction/reaction.h b/src/network/include/gridfire/reaction/reaction.h index 1adc193b..f964e9e4 100644 --- a/src/network/include/gridfire/reaction/reaction.h +++ b/src/network/include/gridfire/reaction/reaction.h @@ -411,6 +411,8 @@ namespace gridfire::reaction { */ explicit TemplatedReactionSet(std::vector reactions); + TemplatedReactionSet(); + /** * @brief Copy constructor. * @param other The ReactionSet to copy. @@ -577,6 +579,9 @@ namespace gridfire::reaction { } } + template + TemplatedReactionSet::TemplatedReactionSet() {} + template TemplatedReactionSet::TemplatedReactionSet(const TemplatedReactionSet &other) { m_reactions.reserve(other.m_reactions.size()); diff --git a/src/network/include/gridfire/screening/screening_weak.h b/src/network/include/gridfire/screening/screening_weak.h index 94c82e4f..f5e25dea 100644 --- a/src/network/include/gridfire/screening/screening_weak.h +++ b/src/network/include/gridfire/screening/screening_weak.h @@ -145,7 +145,7 @@ namespace gridfire::screening { const T T9, const T rho ) const { - LOG_TRACE_L1( + LOG_TRACE_L3( m_logger, "Calculating weak screening factors for {} reactions...", reactions.size() diff --git a/src/network/include/gridfire/solver/solver.h b/src/network/include/gridfire/solver/solver.h index 2f65ae60..93c8c496 100644 --- a/src/network/include/gridfire/solver/solver.h +++ b/src/network/include/gridfire/solver/solver.h @@ -205,6 +205,7 @@ namespace gridfire::solver { const Eigen::VectorXd& m_Y_QSE; ///< Steady-state abundances of the QSE species. const double m_T9; ///< Temperature in units of 10^9 K. const double m_rho; ///< Density in g/cm^3. + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance for trace and debug information. bool m_isViewInitialized = false; @@ -316,11 +317,13 @@ namespace gridfire::solver { }; DynamicEngine& m_engine; ///< The engine used to evaluate the network. + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance for trace and debug information. const std::vector& m_YFull; ///< The full, initial abundance vector const std::vector& m_dynamicSpeciesIndices; ///< Indices of the dynamic species. const std::vector& m_QSESpeciesIndices; ///< Indices of the QSE species. const double m_T9; ///< Temperature in units of 10^9 K. const double m_rho; ///< Density in g/cm^3. + const Eigen::VectorXd& m_YScale; ///< Scaling vector for the QSE species for asinh scaling. /** * @brief Constructor for the EigenFunctor. @@ -337,14 +340,16 @@ namespace gridfire::solver { const std::vector& dynamicSpeciesIndices, const std::vector& QSESpeciesIndices, const double T9, - const double rho + const double rho, + const Eigen::VectorXd& YScale ) : m_engine(engine), m_YFull(YFull), m_dynamicSpeciesIndices(dynamicSpeciesIndices), m_QSESpeciesIndices(QSESpeciesIndices), m_T9(T9), - m_rho(rho) {} + m_rho(rho), + m_YScale(YScale) {} int values() const { return m_QSESpeciesIndices.size(); } int inputs() const { return m_QSESpeciesIndices.size(); } @@ -493,9 +498,20 @@ namespace gridfire::solver { }; template - int QSENetworkSolver::EigenFunctor::operator()(const InputType &v_QSE_log, OutputType &f_QSE) const { + int QSENetworkSolver::EigenFunctor::operator()(const InputType &v_QSE, OutputType &f_QSE) const { std::vector y = m_YFull; // Full vector of species abundances - Eigen::VectorXd Y_QSE = v_QSE_log.array().exp(); + Eigen::VectorXd Y_QSE = m_YScale.array() * v_QSE.array().sinh(); // Convert to physical abundances using asinh scaling + LOG_DEBUG( + m_logger, + "Calling QSENetworkSolver::EigenFunctor::operator() with QSE abundances {}", + [&]() -> std::string { + std::stringstream ss; + ss << std::scientific << std::setprecision(5); + for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { + ss << std::string(m_engine.getNetworkSpecies()[m_QSESpeciesIndices[i]].name()) << ": " << Y_QSE(i) << ", "; + } + return ss.str(); + }()); for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { y[m_QSESpeciesIndices[i]] = Y_QSE(i); @@ -503,6 +519,17 @@ namespace gridfire::solver { const auto [dydt, specificEnergyRate] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho); + LOG_DEBUG( + m_logger, + "Calling QSENetworkSolver::EigenFunctor::operator() found QSE dydt {}", + [&]() -> std::string { + std::stringstream ss; + ss << std::scientific << std::setprecision(5); + for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { + ss << std::string(m_engine.getNetworkSpecies()[m_QSESpeciesIndices[i]].name()) << ": " << dydt[m_QSESpeciesIndices[i]] << ", "; + } + return ss.str(); + }()); f_QSE.resize(m_QSESpeciesIndices.size()); for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { f_QSE(i) = dydt[m_QSESpeciesIndices[i]]; @@ -511,9 +538,9 @@ namespace gridfire::solver { } template - int QSENetworkSolver::EigenFunctor::df(const InputType& v_QSE_log, JacobianType& J_QSE) const { + int QSENetworkSolver::EigenFunctor::df(const InputType& v_QSE, JacobianType& J_QSE) const { std::vector y = m_YFull; - Eigen::VectorXd Y_QSE = v_QSE_log.array().exp(); + Eigen::VectorXd Y_QSE = m_YScale.array() * v_QSE.array().sinh(); for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { y[m_QSESpeciesIndices[i]] = Y_QSE(i); @@ -528,8 +555,32 @@ namespace gridfire::solver { } } + std::string format_jacobian = [&]() -> std::string { + std::stringstream ss; + ss << std::scientific << std::setprecision(5); + for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { + for (size_t j = 0; j < m_QSESpeciesIndices.size(); ++j) { + ss << J_QSE(i, j) << " "; + } + ss << "\n"; + } + return ss.str(); + }(); + + for (long j = 0; j < J_QSE.cols(); ++j) { - J_QSE(j, j) *= Y_QSE(j); // chain rule for log space transformation + LOG_DEBUG( + m_logger, + "Jacobian ({},{}) before chain rule: {}", + j, j, J_QSE(j, j) + ); + const double dY_dv = m_YScale(j) * CppAD::cosh(v_QSE(j)); + J_QSE.col(j) *= dY_dv; // chain rule for log space transformation + LOG_DEBUG( + m_logger, + "Jacobian ({},{}) after chain rule: {} (dY/dv = {})", + j, j, J_QSE(j, j), dY_dv + ); } return 0; } diff --git a/src/network/lib/engine/engine_graph.cpp b/src/network/lib/engine/engine_graph.cpp index 8845372c..2de28e86 100644 --- a/src/network/lib/engine/engine_graph.cpp +++ b/src/network/lib/engine/engine_graph.cpp @@ -26,12 +26,7 @@ namespace gridfire { GraphEngine::GraphEngine( const fourdst::composition::Composition &composition - ): - m_reactions(build_reaclib_nuclear_network(composition, false)), - m_partitionFunction(std::make_unique()){ - syncInternalMaps(); - precomputeNetwork(); - } + ): GraphEngine(composition, partition::GroundStatePartitionFunction()) {} GraphEngine::GraphEngine( const fourdst::composition::Composition &composition, @@ -40,7 +35,6 @@ namespace gridfire { m_partitionFunction(partitionFunction.clone()) // Clone the partition function to ensure ownership { syncInternalMaps(); - precomputeNetwork(); } GraphEngine::GraphEngine( @@ -48,7 +42,6 @@ namespace gridfire { ) : m_reactions(reactions) { syncInternalMaps(); - precomputeNetwork(); } StepDerivatives GraphEngine::calculateRHSAndEnergy( @@ -58,13 +51,17 @@ namespace gridfire { ) const { 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()); + for (const auto& reaction: m_reactions) { bare_rates.push_back(reaction.calculate_rate(T9)); + bare_reverse_rates.push_back(calculateReverseRate(reaction, T9)); } // --- The public facing interface can always use the precomputed version since taping is done internally --- - return calculateAllDerivativesUsingPrecomputation(Y, bare_rates, T9, rho); + return calculateAllDerivativesUsingPrecomputation(Y, bare_rates, bare_reverse_rates, T9, rho); } else { return calculateAllDerivatives(Y, T9, rho); } @@ -72,12 +69,17 @@ namespace gridfire { void GraphEngine::syncInternalMaps() { + LOG_INFO(m_logger, "Synchronizing internal maps for REACLIB graph network (serif::network::GraphNetwork)..."); collectNetworkSpecies(); populateReactionIDMap(); populateSpeciesToIndexMap(); + collectAtomicReverseRateAtomicBases(); generateStoichiometryMatrix(); reserveJacobianMatrix(); recordADTape(); + precomputeNetwork(); + LOG_INFO(m_logger, "Internal maps synchronized. Network contains {} species and {} reactions.", + m_networkSpecies.size(), m_reactions.size()); } // --- Network Graph Construction Methods --- @@ -234,17 +236,22 @@ namespace gridfire { double GraphEngine::calculateReverseRate( const reaction::Reaction &reaction, - const double T9, - const double expFactor + const double T9 ) const { + if (!m_useReverseReactions) { + LOG_TRACE_L2(m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id()); + return 0.0; // If reverse reactions are not used, return 0.0 + } + const double expFactor = std::exp(-reaction.qValue() / (m_constants.kB * T9)); double reverseRate = 0.0; const double forwardRate = reaction.calculate_rate(T9); if (reaction.reactants().size() == 2 && reaction.products().size() == 2) { reverseRate = calculateReverseRateTwoBody(reaction, T9, forwardRate, expFactor); } else { - LOG_WARNING(m_logger, "Reverse rate calculation for reactions with more than two reactants or products is not implemented."); + LOG_WARNING_LIMIT_EVERY_N(1000000, m_logger, "Reverse rate calculation for reactions with more than two reactants or products is not implemented (reaction id {}).", reaction.peName()); } + LOG_TRACE_L2(m_logger, "Calculated reverse rate for reaction '{}': {:.3E} at T9={:.3E}.", reaction.id(), reverseRate, T9); return reverseRate; } @@ -298,7 +305,9 @@ namespace gridfire { ); // Accumulate partition functions - auto pf_op = [&](double acc, const auto& species) { return acc * m_partitionFunction->evaluate(species.z(), species.a(), T9); }; + auto pf_op = [&](double acc, const auto& species) { + return acc * m_partitionFunction->evaluate(species.z(), species.a(), T9); + }; const double partitionFunctionNumerator = std::accumulate( reaction.reactants().begin(), reaction.reactants().end(), @@ -312,9 +321,14 @@ namespace gridfire { pf_op ); - const double CT = std::pow(massNumerator/massDenominator, 1.5) * (partitionFunctionNumerator/partitionFunctionDenominator); + const double CT = std::pow(massNumerator/massDenominator, 1.5) * + (partitionFunctionNumerator/partitionFunctionDenominator); - return forwardRate * symmetryFactor * CT * expFactor; + double reverseRate = forwardRate * symmetryFactor * CT * expFactor; + if (!std::isfinite(reverseRate)) { + return 0.0; // If the reverse rate is not finite, return 0.0 + } + return reverseRate; // Return the calculated reverse rate } @@ -323,6 +337,10 @@ namespace gridfire { const double T9, const double reverseRate ) const { + if (!m_useReverseReactions) { + LOG_TRACE_L2(m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate derivative of reaction '{}'.", reaction.id()); + return 0.0; // If reverse reactions are not used, return 0.0 + } const double d_log_kFwd = reaction.calculate_forward_rate_log_derivative(T9); auto log_deriv_pf_op = [&](double acc, const auto& species) { @@ -355,9 +373,30 @@ namespace gridfire { } + bool GraphEngine::isUsingReverseReactions() const { + return m_useReverseReactions; + } + + void GraphEngine::setUseReverseReactions(const bool useReverse) { + m_useReverseReactions = useReverse; + } + + int GraphEngine::getSpeciesIndex(const fourdst::atomic::Species &species) const { + return m_speciesToIndexMap.at(species); // Returns the index of the species in the stoichiometry matrix + } + + std::vector GraphEngine::mapNetInToMolarAbundanceVector(const NetIn &netIn) const { + std::vector Y(m_networkSpecies.size(), 0.0); // Initialize with zeros + for (const auto& [symbol, entry] : netIn.composition) { + Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance + } + return Y; // Return the vector of molar abundances + } + StepDerivatives GraphEngine::calculateAllDerivativesUsingPrecomputation( const std::vector &Y_in, const std::vector &bare_rates, + const std::vector &bare_reverse_rates, const double T9, const double rho ) const { @@ -396,15 +435,31 @@ namespace gridfire { 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 molarReactionFlow = + const double forwardMolarReactionFlow = screeningFactor * bare_rate * precomp.symmetry_factor * abundanceProduct * - std::pow(rho, numReactants); + std::pow(rho, numReactants > 1 ? numReactants - 1 : 0); + + 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]; + double reverseAbundanceProduct = 1.0; + for (size_t i = 0; i < precomp.unique_product_indices.size(); ++i) { + reverseAbundanceProduct *= std::pow(Y_in[precomp.unique_product_indices[i]], precomp.product_powers[i]); + } + reverseMolarReactionFlow = screeningFactor * + bare_reverse_rate * + precomp.reverse_symmetry_factor * + reverseAbundanceProduct * + std::pow(rho, numProducts > 1 ? numProducts - 1 : 0); + } + + molarReactionFlows.push_back(forwardMolarReactionFlow - reverseMolarReactionFlow); - molarReactionFlows.push_back(molarReactionFlow); } // --- Assemble molar abundance derivatives --- @@ -523,7 +578,7 @@ namespace gridfire { } void GraphEngine::generateJacobianMatrix( - const std::vector &Y, + const std::vector &Y_dynamic, const double T9, const double rho ) { @@ -534,10 +589,24 @@ namespace gridfire { // 1. Pack the input variables into a vector for CppAD std::vector adInput(numSpecies + 2, 0.0); // +2 for T9 and rho for (size_t i = 0; i < numSpecies; ++i) { - adInput[i] = Y[i]; + adInput[i] = Y_dynamic[i]; } adInput[numSpecies] = T9; // T9 adInput[numSpecies + 1] = rho; // rho + LOG_DEBUG( + m_logger, + "AD Input to jacobian {}", + [&]() -> std::string { + std::stringstream ss; + ss << std::scientific << std::setprecision(5); + for (size_t i = 0; i < adInput.size(); ++i) { + ss << adInput[i]; + if (i < adInput.size() - 1) { + ss << ", "; + } + } + return ss.str(); + }()); // 2. Calculate the full jacobian const std::vector dotY = m_rhsADFun.Jacobian(adInput); @@ -552,10 +621,28 @@ namespace gridfire { } } } + LOG_DEBUG( + m_logger, + "Final Jacobian is:\n{}", + [&]() -> std::string { + std::stringstream ss; + ss << std::scientific << std::setprecision(5); + for (size_t i = 0; i < m_jacobianMatrix.size1(); ++i) { + for (size_t j = 0; j < m_jacobianMatrix.size2(); ++j) { + ss << m_jacobianMatrix(i, j); + if (j < m_jacobianMatrix.size2() - 1) { + ss << ", "; + } + } + ss << "\n"; + } + return ss.str(); + }()); LOG_TRACE_L1(m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); } 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); } @@ -674,8 +761,11 @@ namespace gridfire { LOG_TRACE_L1(m_logger, "Successfully exported network graph to {}", filename); } - std::unordered_map GraphEngine::getSpeciesTimescales(const std::vector &Y, const double T9, - const double rho) const { + std::unordered_map GraphEngine::getSpeciesTimescales( + const std::vector &Y, + const double T9, + const double rho + ) const { auto [dydt, _] = calculateAllDerivatives(Y, T9, rho); std::unordered_map speciesTimescales; speciesTimescales.reserve(m_networkSpecies.size()); @@ -739,6 +829,19 @@ namespace gridfire { adInput.size()); } + void GraphEngine::collectAtomicReverseRateAtomicBases() { + m_atomicReverseRates.clear(); + m_atomicReverseRates.reserve(m_reactions.size()); + + for (const auto& reaction: m_reactions) { + if (reaction.qValue() != 0.0) { + m_atomicReverseRates.push_back(std::make_unique(reaction, *this)); + } else { + m_atomicReverseRates.push_back(nullptr); + } + } + } + void GraphEngine::precomputeNetwork() { LOG_TRACE_L1(m_logger, "Pre-computing constant components of GraphNetwork state..."); @@ -756,7 +859,7 @@ namespace gridfire { PrecomputedReaction precomp; precomp.reaction_index = i; - // --- Precompute reactant information --- + // --- Precompute forward reaction information --- // Count occurrences for each reactant to determine powers and symmetry std::unordered_map reactantCounts; for (const auto& reactant: reaction.reactants()) { @@ -769,10 +872,30 @@ namespace gridfire { precomp.unique_reactant_indices.push_back(index); precomp.reactant_powers.push_back(count); - symmetryDenominator *= 1.0/std::tgamma(count + 1); + symmetryDenominator *= std::tgamma(count + 1); } - precomp.symmetry_factor = symmetryDenominator; + precomp.symmetry_factor = 1.0/symmetryDenominator; + + // --- Precompute reverse reaction information --- + if (reaction.qValue() != 0.0) { + std::unordered_map productCounts; + for (const auto& product : reaction.products()) { + productCounts[speciesIndexMap.at(product)]++; + } + double reverseSymmetryDenominator = 1.0; + for (const auto& [index, count] : productCounts) { + precomp.unique_product_indices.push_back(index); + precomp.product_powers.push_back(count); + reverseSymmetryDenominator *= std::tgamma(count + 1); + } + + precomp.reverse_symmetry_factor = 1.0/reverseSymmetryDenominator; + } else { + precomp.unique_product_indices.clear(); + precomp.product_powers.clear(); + precomp.reverse_symmetry_factor = 0.0; // No reverse reaction for Q = 0 reactions + } // --- Precompute stoichiometry information --- const auto stoichiometryMap = reaction.stoichiometry(); @@ -788,7 +911,7 @@ namespace gridfire { } } - bool AtomicReverseRate::forward( + bool GraphEngine::AtomicReverseRate::forward( const size_t p, const size_t q, const CppAD::vector &vx, @@ -796,12 +919,53 @@ namespace gridfire { const CppAD::vector &tx, CppAD::vector &ty ) { + + if ( p != 0) { return false; } const double T9 = tx[0]; - const double expFactor = std::exp(-m_reaction.qValue() / (m_kB * T9 * 1e9)); // Convert MeV to erg - if (p == 0) { - // --- Zeroth order forward sweep --- - const auto k_rev = m_engine.calculateReverseRate(m_reaction, T9, expFactor); + + const double reverseRate = m_engine.calculateReverseRate(m_reaction, T9); + // std::cout << m_reaction.peName() << " reverseRate: " << reverseRate << " at T9: " << T9 << "\n"; + ty[0] = reverseRate; // Store the reverse rate in the output vector + + if (vx.size() > 0) { + vy[0] = vx[0]; } + return true; } + bool GraphEngine::AtomicReverseRate::reverse( + size_t q, + const CppAD::vector &tx, + const CppAD::vector &ty, + CppAD::vector &px, + const CppAD::vector &py + ) { + const double T9 = tx[0]; + const double reverseRate = ty[0]; + + const double derivative = m_engine.calculateReverseRateTwoBodyDerivative(m_reaction, T9, reverseRate); + // std::cout << m_reaction.peName() << " reverseRate Derivative: " << derivative << "\n"; + + px[0] = py[0] * derivative; // Return the derivative of the reverse rate with respect to T9 + + return true; + } + + bool GraphEngine::AtomicReverseRate::for_sparse_jac( + size_t q, + const CppAD::vector> &r, + CppAD::vector> &s + ) { + s[0] = r[0]; + return true; + } + + bool GraphEngine::AtomicReverseRate::rev_sparse_jac( + size_t q, + const CppAD::vector> &rt, + CppAD::vector> &st + ) { + st[0] = rt[0]; + return true; + } } diff --git a/src/network/lib/engine/procedures/priming.cpp b/src/network/lib/engine/procedures/priming.cpp new file mode 100644 index 00000000..b1850a30 --- /dev/null +++ b/src/network/lib/engine/procedures/priming.cpp @@ -0,0 +1,84 @@ +#include "gridfire/engine/procedures/priming.h" +#include "gridfire/engine/views/engine_priming.h" +#include "gridfire/solver/solver.h" + +#include "gridfire/engine/engine_abstract.h" +#include "gridfire/network.h" + +namespace gridfire { + fourdst::composition::Composition primeNetwork(const NetIn& netIn, DynamicEngine& engine) { + using fourdst::composition::Composition; + using fourdst::atomic::Species; + + + NetIn currentConditions = netIn; + std::vector speciesToPrime; + for (const auto &entry: netIn.composition | std::views::values) { + if (entry.mass_fraction() == 0.0) { + speciesToPrime.push_back(entry.isotope()); + } + } + + + if (speciesToPrime.empty()) { + return netIn.composition; // No priming needed + } + + const double T9 = netIn.temperature / 1e9; // Convert temperature to T9 + const double rho = netIn.density; // Density in g/cm^3 + + 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(); + + currentConditions.composition.setMassFraction( + std::string(primingSpecies.name()), + equilibriumMassFraction + ); + } + + return currentConditions.composition; + } + + double calculateDestructionRateConstant( + const DynamicEngine& engine, + const fourdst::atomic::Species& species, + const std::vector& Y, + const double T9, + const double rho + ) { + const int speciesIndex = engine.getSpeciesIndex(species); + std::vector Y_scaled(Y.begin(), Y.end()); + Y_scaled[speciesIndex] = 1.0; // Set the abundance of the species to 1.0 for rate constant calculation + double destructionRateConstant = 0.0; + for (const auto& reaction: engine.getNetworkReactions()) { + if (reaction.contains_reactant(species)) { + const int stoichiometry = reaction.stoichiometry(species); + destructionRateConstant += std::abs(stoichiometry) * engine.calculateMolarReactionFlow(reaction, Y_scaled, T9, rho); + } + } + return destructionRateConstant; + } + + double calculateCreationRate( + const DynamicEngine& engine, + const fourdst::atomic::Species& species, + const std::vector& Y, + const double T9, + const double rho + ) { + double creationRate = 0.0; + for (const auto& reaction: engine.getNetworkReactions()) { + const int stoichiometry = reaction.stoichiometry(species); + if (stoichiometry > 0) { + creationRate += stoichiometry * engine.calculateMolarReactionFlow(reaction, Y, T9, rho); + } + } + return creationRate; + } +} \ No newline at end of file diff --git a/src/network/lib/engine/views/engine_adaptive.cpp b/src/network/lib/engine/views/engine_adaptive.cpp index 5a0763c5..9cd09ebf 100644 --- a/src/network/lib/engine/views/engine_adaptive.cpp +++ b/src/network/lib/engine/views/engine_adaptive.cpp @@ -136,12 +136,12 @@ namespace gridfire { } void AdaptiveEngineView::generateJacobianMatrix( - const std::vector &Y_culled, + const std::vector &Y_dynamic, const double T9, const double rho ) { validateState(); - const auto Y_full = mapCulledToFull(Y_culled); + const auto Y_full = mapCulledToFull(Y_dynamic); m_baseEngine.generateJacobianMatrix(Y_full, T9, rho); } @@ -221,6 +221,25 @@ namespace gridfire { return m_baseEngine.getScreeningModel(); } + std::vector AdaptiveEngineView::mapNetInToMolarAbundanceVector(const NetIn &netIn) const { + std::vector Y(m_activeSpecies.size(), 0.0); // Initialize with zeros + for (const auto& [symbol, entry] : netIn.composition) { + Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance + } + return Y; // Return the vector of molar abundances + } + + 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()) { + return static_cast(std::distance(m_activeSpecies.begin(), it)); + } else { + LOG_ERROR(m_logger, "Species '{}' not found in active species list.", species.name()); + m_logger->flush_log(); + throw std::runtime_error("Species not found in active species list: " + std::string(species.name())); + } + } + std::vector AdaptiveEngineView::mapCulledToFull(const std::vector& culled) const { std::vector full(m_baseEngine.getNetworkSpecies().size(), 0.0); for (size_t i_culled = 0; i_culled < culled.size(); ++i_culled) { diff --git a/src/network/lib/engine/views/engine_defined.cpp b/src/network/lib/engine/views/engine_defined.cpp index c8fc781f..d5a922de 100644 --- a/src/network/lib/engine/views/engine_defined.cpp +++ b/src/network/lib/engine/views/engine_defined.cpp @@ -1,32 +1,88 @@ #include "gridfire/engine/views/engine_defined.h" +#include + #include "quill/LogMacros.h" +#include +#include +#include +#include +#include +#include + namespace gridfire { using fourdst::atomic::Species; - FileDefinedEngineView::FileDefinedEngineView( - DynamicEngine &baseEngine, - const std::string &fileName, - const io::NetworkFileParser &parser - ): - m_baseEngine(baseEngine), - m_fileName(fileName), - m_parser(parser), - m_activeSpecies(baseEngine.getNetworkSpecies()), - m_activeReactions(baseEngine.getNetworkReactions()) { - buildFromFile(fileName); + DefinedEngineView::DefinedEngineView(const std::vector& peNames, DynamicEngine& baseEngine) : + m_baseEngine(baseEngine) { + m_activeReactions.clear(); + m_activeSpecies.clear(); + + std::unordered_set seenSpecies; + + const auto& fullNetworkReactionSet = m_baseEngine.getNetworkReactions(); + for (const auto& peName : peNames) { + if (!fullNetworkReactionSet.contains(peName)) { + LOG_ERROR(m_logger, "Reaction with name '{}' not found in the base engine's network reactions. Aborting...", peName); + m_logger->flush_log(); + throw std::runtime_error("Reaction with name '" + std::string(peName) + "' not found in the base engine's network reactions."); + } + auto reaction = fullNetworkReactionSet[peName]; + for (const auto& reactant : reaction.reactants()) { + if (!seenSpecies.contains(reactant)) { + seenSpecies.insert(reactant); + m_activeSpecies.push_back(reactant); + } + } + for (const auto& product : reaction.products()) { + if (!seenSpecies.contains(product)) { + seenSpecies.insert(product); + m_activeSpecies.push_back(product); + } + } + m_activeReactions.add_reaction(reaction); + } + LOG_TRACE_L1(m_logger, "DefinedEngineView built with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size()); + LOG_DEBUG(m_logger, "Active species: {}", [this]() -> std::string { + std::string result; + for (const auto& species : m_activeSpecies) { + result += std::string(species.name()) + ", "; + } + if (!result.empty()) { + result.pop_back(); // Remove last space + result.pop_back(); // Remove last comma + } + return result; + }()); + LOG_DEBUG(m_logger, "Active reactions: {}", [this]() -> std::string { + std::string result; + for (const auto& reaction : m_activeReactions) { + result += std::string(reaction.id()) + ", "; + } + if (!result.empty()) { + result.pop_back(); // Remove last space + result.pop_back(); // Remove last comma + } + return result; + }()); + m_speciesIndexMap = constructSpeciesIndexMap(); + m_reactionIndexMap = constructReactionIndexMap(); + m_isStale = false; + } - const DynamicEngine & FileDefinedEngineView::getBaseEngine() const { + const DynamicEngine & DefinedEngineView::getBaseEngine() const { return m_baseEngine; } - const std::vector & FileDefinedEngineView::getNetworkSpecies() const { + + + const std::vector & DefinedEngineView::getNetworkSpecies() const { return m_activeSpecies; } - StepDerivatives FileDefinedEngineView::calculateRHSAndEnergy( + StepDerivatives DefinedEngineView::calculateRHSAndEnergy( const std::vector &Y_defined, const double T9, const double rho @@ -42,18 +98,18 @@ namespace gridfire { return definedResults; } - void FileDefinedEngineView::generateJacobianMatrix( - const std::vector &Y_defined, + void DefinedEngineView::generateJacobianMatrix( + const std::vector &Y_dynamic, const double T9, const double rho ) { validateNetworkState(); - const auto Y_full = mapViewToFull(Y_defined); + const auto Y_full = mapViewToFull(Y_dynamic); m_baseEngine.generateJacobianMatrix(Y_full, T9, rho); } - double FileDefinedEngineView::getJacobianMatrixEntry( + double DefinedEngineView::getJacobianMatrixEntry( const int i_defined, const int j_defined ) const { @@ -65,13 +121,13 @@ namespace gridfire { return m_baseEngine.getJacobianMatrixEntry(i_full, j_full); } - void FileDefinedEngineView::generateStoichiometryMatrix() { + void DefinedEngineView::generateStoichiometryMatrix() { validateNetworkState(); m_baseEngine.generateStoichiometryMatrix(); } - int FileDefinedEngineView::getStoichiometryMatrixEntry( + int DefinedEngineView::getStoichiometryMatrixEntry( const int speciesIndex_defined, const int reactionIndex_defined ) const { @@ -82,7 +138,7 @@ namespace gridfire { return m_baseEngine.getStoichiometryMatrixEntry(i_full, j_full); } - double FileDefinedEngineView::calculateMolarReactionFlow( + double DefinedEngineView::calculateMolarReactionFlow( const reaction::Reaction &reaction, const std::vector &Y_defined, const double T9, @@ -91,7 +147,7 @@ namespace gridfire { validateNetworkState(); if (!m_activeReactions.contains(reaction)) { - LOG_ERROR(m_logger, "Reaction '{}' is not part of the active reactions in the file defined engine view.", reaction.id()); + LOG_ERROR(m_logger, "Reaction '{}' is not part of the active reactions in the DefinedEngineView.", reaction.id()); m_logger -> flush_log(); throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id())); } @@ -99,13 +155,13 @@ namespace gridfire { return m_baseEngine.calculateMolarReactionFlow(reaction, Y_full, T9, rho); } - const reaction::LogicalReactionSet & FileDefinedEngineView::getNetworkReactions() const { + const reaction::LogicalReactionSet & DefinedEngineView::getNetworkReactions() const { validateNetworkState(); return m_activeReactions; } - std::unordered_map FileDefinedEngineView::getSpeciesTimescales( + std::unordered_map DefinedEngineView::getSpeciesTimescales( const std::vector &Y_defined, const double T9, const double rho @@ -124,28 +180,45 @@ namespace gridfire { return definedTimescales; } - void FileDefinedEngineView::update(const NetIn &netIn) { - if (m_isStale) { - buildFromFile(m_fileName); - } + void DefinedEngineView::update(const NetIn &netIn) { + return; } - void FileDefinedEngineView::setNetworkFile(const std::string &fileName) { - m_isStale = true; - m_fileName = fileName; - LOG_DEBUG(m_logger, "File '{}' set to '{}'. FileDefinedNetworkView is now stale! You MUST call update() before you use it!", m_fileName, fileName); - } - void FileDefinedEngineView::setScreeningModel(const screening::ScreeningType model) { + void DefinedEngineView::setScreeningModel(const screening::ScreeningType model) { m_baseEngine.setScreeningModel(model); } - screening::ScreeningType FileDefinedEngineView::getScreeningModel() const { + screening::ScreeningType DefinedEngineView::getScreeningModel() const { return m_baseEngine.getScreeningModel(); } - std::vector FileDefinedEngineView::constructSpeciesIndexMap() const { - LOG_TRACE_L1(m_logger, "Constructing species index map for file defined engine view..."); + int DefinedEngineView::getSpeciesIndex(const Species &species) const { + validateNetworkState(); + + auto it = std::find(m_activeSpecies.begin(), m_activeSpecies.end(), species); + if (it != m_activeSpecies.end()) { + return static_cast(std::distance(m_activeSpecies.begin(), it)); + } else { + LOG_ERROR(m_logger, "Species '{}' not found in active species list.", species.name()); + m_logger->flush_log(); + throw std::runtime_error("Species not found in active species list: " + std::string(species.name())); + } + } + + std::vector DefinedEngineView::mapNetInToMolarAbundanceVector(const NetIn &netIn) const { + std::vector Y(m_activeSpecies.size(), 0.0); // Initialize with zeros + for (const auto& [symbol, entry] : netIn.composition) { + auto it = std::ranges::find(m_activeSpecies, entry.isotope()); + if (it != m_activeSpecies.end()) { + Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance + } + } + return Y; // Return the vector of molar abundances + } + + std::vector DefinedEngineView::constructSpeciesIndexMap() const { + LOG_TRACE_L1(m_logger, "Constructing species index map for DefinedEngineView..."); std::unordered_map fullSpeciesReverseMap; const auto& fullSpeciesList = m_baseEngine.getNetworkSpecies(); @@ -173,8 +246,8 @@ namespace gridfire { } - std::vector FileDefinedEngineView::constructReactionIndexMap() const { - LOG_TRACE_L1(m_logger, "Constructing reaction index map for file defined engine view..."); + std::vector DefinedEngineView::constructReactionIndexMap() const { + LOG_TRACE_L1(m_logger, "Constructing reaction index map for DefinedEngineView..."); // --- Step 1: Create a reverse map using the reaction's unique ID as the key. --- std::unordered_map fullReactionReverseMap; @@ -205,66 +278,7 @@ namespace gridfire { return reactionIndexMap; } - void FileDefinedEngineView::buildFromFile(const std::string &fileName) { - LOG_TRACE_L1(m_logger, "Building file defined engine view from {}...", fileName); - auto [reactionPENames] = m_parser.parse(fileName); - - m_activeReactions.clear(); - m_activeSpecies.clear(); - - std::unordered_set seenSpecies; - - const auto& fullNetworkReactionSet = m_baseEngine.getNetworkReactions(); - for (const auto& peName : reactionPENames) { - if (!fullNetworkReactionSet.contains(peName)) { - LOG_ERROR(m_logger, "Reaction with name '{}' not found in the base engine's network reactions. Aborting...", peName); - m_logger->flush_log(); - throw std::runtime_error("Reaction with name '" + std::string(peName) + "' not found in the base engine's network reactions."); - } - auto reaction = fullNetworkReactionSet[peName]; - for (const auto& reactant : reaction.reactants()) { - if (!seenSpecies.contains(reactant)) { - seenSpecies.insert(reactant); - m_activeSpecies.push_back(reactant); - } - } - for (const auto& product : reaction.products()) { - if (!seenSpecies.contains(product)) { - seenSpecies.insert(product); - m_activeSpecies.push_back(product); - } - } - m_activeReactions.add_reaction(reaction); - } - LOG_TRACE_L1(m_logger, "File defined engine view built with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size()); - LOG_DEBUG(m_logger, "Active species: {}", [this]() -> std::string { - std::string result; - for (const auto& species : m_activeSpecies) { - result += std::string(species.name()) + ", "; - } - if (!result.empty()) { - result.pop_back(); // Remove last space - result.pop_back(); // Remove last comma - } - return result; - }()); - LOG_DEBUG(m_logger, "Active reactions: {}", [this]() -> std::string { - std::string result; - for (const auto& reaction : m_activeReactions) { - result += std::string(reaction.id()) + ", "; - } - if (!result.empty()) { - result.pop_back(); // Remove last space - result.pop_back(); // Remove last comma - } - return result; - }()); - m_speciesIndexMap = constructSpeciesIndexMap(); - m_reactionIndexMap = constructReactionIndexMap(); - m_isStale = false; - } - - std::vector FileDefinedEngineView::mapViewToFull(const std::vector& culled) const { + std::vector DefinedEngineView::mapViewToFull(const std::vector& culled) const { std::vector full(m_baseEngine.getNetworkSpecies().size(), 0.0); for (size_t i_culled = 0; i_culled < culled.size(); ++i_culled) { const size_t i_full = m_speciesIndexMap[i_culled]; @@ -273,7 +287,7 @@ namespace gridfire { return full; } - std::vector FileDefinedEngineView::mapFullToView(const std::vector& full) const { + std::vector DefinedEngineView::mapFullToView(const std::vector& full) const { std::vector culled(m_activeSpecies.size(), 0.0); for (size_t i_culled = 0; i_culled < m_activeSpecies.size(); ++i_culled) { const size_t i_full = m_speciesIndexMap[i_culled]; @@ -282,7 +296,7 @@ namespace gridfire { return culled; } - size_t FileDefinedEngineView::mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const { + size_t DefinedEngineView::mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const { if (culledSpeciesIndex < 0 || culledSpeciesIndex >= static_cast(m_speciesIndexMap.size())) { LOG_ERROR(m_logger, "Defined index {} is out of bounds for species index map of size {}.", culledSpeciesIndex, m_speciesIndexMap.size()); m_logger->flush_log(); @@ -291,7 +305,7 @@ namespace gridfire { return m_speciesIndexMap[culledSpeciesIndex]; } - size_t FileDefinedEngineView::mapViewToFullReactionIndex(size_t culledReactionIndex) const { + size_t DefinedEngineView::mapViewToFullReactionIndex(size_t culledReactionIndex) const { if (culledReactionIndex < 0 || culledReactionIndex >= static_cast(m_reactionIndexMap.size())) { LOG_ERROR(m_logger, "Defined index {} is out of bounds for reaction index map of size {}.", culledReactionIndex, m_reactionIndexMap.size()); m_logger->flush_log(); @@ -300,11 +314,25 @@ namespace gridfire { return m_reactionIndexMap[culledReactionIndex]; } - void FileDefinedEngineView::validateNetworkState() const { + void DefinedEngineView::validateNetworkState() const { if (m_isStale) { - LOG_ERROR(m_logger, "File defined engine view is stale. Please call update() with a valid NetIn object."); + LOG_ERROR(m_logger, "DefinedEngineView is stale. Please call update() with a valid NetIn object."); m_logger->flush_log(); - throw std::runtime_error("File defined engine view is stale. Please call update() with a valid NetIn object."); + throw std::runtime_error("DefinedEngineView is stale. Please call update() with a valid NetIn object."); } } + + + //////////////////////////////////////////// + /// FileDefinedEngineView Implementation /// + ///////////////////////////////////////////// + + FileDefinedEngineView::FileDefinedEngineView( + DynamicEngine &baseEngine, + const std::string &fileName, + const io::NetworkFileParser &parser + ): + DefinedEngineView(parser.parse(fileName), baseEngine), + m_fileName(fileName), + m_parser(parser) {} } diff --git a/src/network/lib/engine/views/engine_multiscale.cpp b/src/network/lib/engine/views/engine_multiscale.cpp new file mode 100644 index 00000000..4d8f5509 --- /dev/null +++ b/src/network/lib/engine/views/engine_multiscale.cpp @@ -0,0 +1,623 @@ +#include "gridfire/engine/views/engine_multiscale.h" + +#include +#include +#include +#include + +namespace gridfire { + using fourdst::atomic::Species; + + MultiscalePartitioningEngineView::MultiscalePartitioningEngineView( + GraphEngine& baseEngine + ) : m_baseEngine(baseEngine) {} + + const std::vector & MultiscalePartitioningEngineView::getNetworkSpecies() const { + return m_baseEngine.getNetworkSpecies(); + } + + StepDerivatives MultiscalePartitioningEngineView::calculateRHSAndEnergy( + const std::vector &Y, + const double T9, + const double rho + ) const { + return m_baseEngine.calculateRHSAndEnergy(Y, T9, rho); + } + + void MultiscalePartitioningEngineView::generateJacobianMatrix( + const std::vector &Y_dynamic, + const double T9, + const double rho + ) { + std::vector Y_full(m_baseEngine.getNetworkSpecies().size(), 0.0); + for (size_t i = 0; i < m_dynamic_species_indices.size(); ++i) { + Y_full[m_dynamic_species_indices[i]] = Y_dynamic[i]; + } + // solveQSEAbundances(Y_full, T9, rho); + m_baseEngine.generateJacobianMatrix(Y_dynamic, T9, rho); + } + + double MultiscalePartitioningEngineView::getJacobianMatrixEntry( + const int i, + const int j + ) const { + return m_baseEngine.getJacobianMatrixEntry(i, j); + } + + void MultiscalePartitioningEngineView::generateStoichiometryMatrix() { + m_baseEngine.generateStoichiometryMatrix(); + } + + int MultiscalePartitioningEngineView::getStoichiometryMatrixEntry( + const int speciesIndex, + const int reactionIndex + ) const { + return m_baseEngine.getStoichiometryMatrixEntry(speciesIndex, reactionIndex); + } + + double MultiscalePartitioningEngineView::calculateMolarReactionFlow( + const reaction::Reaction &reaction, + const std::vector &Y, + double T9, + double rho + ) const { + return m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho); + } + + const reaction::LogicalReactionSet & MultiscalePartitioningEngineView::getNetworkReactions() const { + return m_baseEngine.getNetworkReactions(); + } + + std::unordered_map MultiscalePartitioningEngineView::getSpeciesTimescales( + const std::vector &Y, + const double T9, + const double rho + ) const { + return m_baseEngine.getSpeciesTimescales(Y, T9, rho); + } + + void MultiscalePartitioningEngineView::update(const NetIn &netIn) { + return; // call partition eventually + } + + void MultiscalePartitioningEngineView::setScreeningModel( + const screening::ScreeningType model + ) { + m_baseEngine.setScreeningModel(model); + } + + screening::ScreeningType MultiscalePartitioningEngineView::getScreeningModel() const { + return m_baseEngine.getScreeningModel(); + } + + const DynamicEngine & MultiscalePartitioningEngineView::getBaseEngine() const { + return m_baseEngine; + } + + void MultiscalePartitioningEngineView::partitionNetwork( + const std::vector &Y, + const double T9, + const double rho, + const double dt_control + ) { + // --- Step 0. Clear previous state --- + m_qse_groups.clear(); + m_dynamic_species.clear(); + m_dynamic_species_indices.clear(); + m_connectivity_graph.clear(); + + // --- Step 1. Identify fast reactions --- + const std::unordered_set fast_reaction_indices = identifyFastReactions(Y, T9, rho, dt_control); + + // --- Step 2. Build connectivity graph based on fast reactions --- + buildConnectivityGraph(fast_reaction_indices); + + // --- Step 3. Find connected components (candidate QSE clusters) --- + findConnectedComponents(); + + // --- Step 4. Validate clusters via flux analysis --- + validateGroupsWithFluxAnalysis(Y, T9, rho); + + // --- Step 5. Identify dynamic species --- + const auto [fst, snd] = identifyDynamicSpecies(Y, m_qse_groups, T9, rho); + m_dynamic_species = fst; + m_dynamic_species_indices = snd; + + } + + void MultiscalePartitioningEngineView::partitionNetwork( + const NetIn &netIn, + const double dt_control + ) { + std::vector Y(m_baseEngine.getNetworkSpecies().size(), 0.0); + for (const auto& [symbol, entry]: netIn.composition ) { + if (!m_baseEngine.involvesSpecies(entry.isotope())) { + LOG_ERROR(m_logger, "Species {} is not part of the network. Exiting...", entry.isotope().name()); + throw std::runtime_error("Species " + std::string(entry.isotope().name()) + " is not part of the network."); + } + const size_t species_index = m_baseEngine.getSpeciesIndex(entry.isotope()); + Y[species_index] = netIn.composition.getMolarAbundance(symbol); + } + + 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); + } + + void MultiscalePartitioningEngineView::exportToDot(const std::string &filename) const { + std::ofstream dotFile(filename); + if (!dotFile.is_open()) { + LOG_ERROR(m_logger, "Failed to open file for writing: {}", filename); + throw std::runtime_error("Failed to open file for writing: " + filename); + } + + dotFile << "digraph PartitionedNetwork {\n"; + dotFile << " graph [rankdir=LR, splines=true, overlap=false, bgcolor=\"#fdfdfd\", label=\"Multiscale Partitioned Network View\"];\n"; + dotFile << " node [shape=circle, style=filled, fillcolor=\"#eafaf1\", fontname=\"Helvetica\"];\n"; + dotFile << " edge [fontname=\"Helvetica\", fontsize=10, color=\"#5d6d7e\"];\n\n"; + + // 1. Define all species from the base network as nodes + const auto& all_species = m_baseEngine.getNetworkSpecies(); + dotFile << " // --- Species Nodes ---\n"; + for (const auto& species : all_species) { + dotFile << " \"" << species.name() << "\" [label=\"" << species.name() << "\"];\n"; + } + dotFile << "\n"; + + // 2. Define all reactions and their connections from the base network + dotFile << " // --- Reaction Edges ---\n"; + for (const auto& reaction : m_baseEngine.getNetworkReactions()) { + std::string reactionNodeId = "reaction_" + std::string(reaction.id()); + dotFile << " \"" << reactionNodeId << "\" [shape=point, fillcolor=black, width=0.08, height=0.08];\n"; + for (const auto& reactant : reaction.reactants()) { + dotFile << " \"" << reactant.name() << "\" -> \"" << reactionNodeId << "\";\n"; + } + for (const auto& product : reaction.products()) { + dotFile << " \"" << reactionNodeId << "\" -> \"" << product.name() << "\";\n"; + } + dotFile << "\n"; + } + + // 3. Draw clusters around the identified QSE groups + dotFile << " // --- QSE Group Clusters ---\n"; + int group_counter = 0; + for (const auto& group : m_qse_groups) { + dotFile << " subgraph cluster_" << group_counter << " {\n"; + dotFile << " style = \"filled,rounded\";\n"; + if (group.is_in_equilibrium) { + dotFile << " label = \"Valid QSE Group" << group_counter + 1 << "\";\n"; + dotFile << " color = \"#fdebd0\";\n"; + } else { + dotFile << " label = \"Non valid QSE Group" << group_counter + 1 << "\";\n"; + dotFile << " color = \"#d5f5e3\";\n"; // Different color for non-equilibrium groups + } + + for (const size_t species_idx : group.species_indices) { + dotFile << " \"" << all_species[species_idx].name() << "\";\n"; + } + + dotFile << " }\n"; + group_counter++; + } + + dotFile << "}\n"; + dotFile.close(); + } + + std::vector MultiscalePartitioningEngineView::mapNetInToMolarAbundanceVector(const NetIn &netIn) const { + std::vector Y(m_dynamic_species.size(), 0.0); // Initialize with zeros + for (const auto& [symbol, entry] : netIn.composition) { + Y[getSpeciesIndex(entry.isotope())] = netIn.composition.getMolarAbundance(symbol); // Map species to their molar abundance + } + return Y; // Return the vector of molar abundances + } + + std::vector MultiscalePartitioningEngineView::getFastSpecies() const { + const auto& all_species = m_baseEngine.getNetworkSpecies(); + std::vector fast_species; + fast_species.reserve(all_species.size() - m_dynamic_species.size()); + for (const auto& species : all_species) { + auto it = std::ranges::find(m_dynamic_species, species); + if (it == m_dynamic_species.end()) { + fast_species.push_back(species); + } + } + return fast_species; + } + + const std::vector & MultiscalePartitioningEngineView::getDynamicSpecies() const { + return m_dynamic_species; + } + + void 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); + } + + void MultiscalePartitioningEngineView::equilibrateNetwork( + const NetIn& netIn, + const double dt_control + ) { + std::vector Y(m_baseEngine.getNetworkSpecies().size(), 0.0); + for (const auto& [symbol, entry]: netIn.composition ) { + if (!m_baseEngine.involvesSpecies(entry.isotope())) { + LOG_ERROR(m_logger, "Species {} is not part of the network. Exiting...", entry.isotope().name()); + throw std::runtime_error("Species " + std::string(entry.isotope().name()) + " is not part of the network."); + } + const size_t species_index = m_baseEngine.getSpeciesIndex(entry.isotope()); + Y[species_index] = netIn.composition.getMolarAbundance(symbol); + } + + 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); + + + } + + int MultiscalePartitioningEngineView::getSpeciesIndex(const fourdst::atomic::Species &species) const { + return m_baseEngine.getSpeciesIndex(species); + } + + std::unordered_set MultiscalePartitioningEngineView::identifyFastReactions( + const std::vector &Y_full, + const double T9, + const double rho, + const double dt_control + ) const { + std::unordered_set fast_reaction_indices; + std::unordered_set fast_species_indices; + + const double timescale_threshold = 0.01 * dt_control; + constexpr int FAST_BOOTSTRAPPING_ITERATIONS = 2; + + const auto& all_reactions = m_baseEngine.getNetworkReactions(); + + // Loop a few times to refine the fast species and reactions. + for (int iteration = 0; iteration < FAST_BOOTSTRAPPING_ITERATIONS; ++iteration) { + // --- Step A: Identify fast species based on their individual reaction timescales --- + for (size_t i = 0; i < all_reactions.size(); ++i) { + const auto& reaction = all_reactions[i]; + + // Calculate the forward molar flow, which represents the rate of destruction. + // Note: We use the base engine's direct calculation method. + const double forward_flow = m_baseEngine.calculateMolarReactionFlow(reaction, Y_full, T9, rho); + + if (forward_flow == 0.0) continue; + + // Determine the timescale for each reactant against this specific reaction. + for (const auto& reactant : reaction.reactants()) { + // ~ Timescale = Abundance / Destruction Rate + + // TODO: Should this call base engine getSpeciesIndex or this view's getSpeciesIndex? + const size_t reactant_idx = m_baseEngine.getSpeciesIndex(reactant); + const double abundance = Y_full[reactant_idx]; + + if (abundance > 0) { + const double timescale = abundance / forward_flow; + if (timescale < timescale_threshold) { + fast_species_indices.insert(reactant_idx); + } + } + } + } + + // --- Step B: Identify fast reactions based on reactants --- + for (size_t i = 0; i < all_reactions.size(); ++i) { + const auto& reaction = all_reactions[i]; + bool all_reactants_are_fast = true; + if (reaction.reactants().empty()) { + all_reactants_are_fast = false; + } + for (const auto& reactant : reaction.reactants()) { + if (!fast_species_indices.contains(m_baseEngine.getSpeciesIndex(reactant))) { + all_reactants_are_fast = false; + break; // No need to check further if one reactant is not fast + } + } + if (all_reactants_are_fast) { + fast_reaction_indices.insert(i); + } + } + + // --- Step C: Propagate "fast" status to products of the fast reactions (this handles things like n-1) --- + for (size_t reaction_idx : fast_reaction_indices) { + for (const auto& product : all_reactions[reaction_idx].products()) { + fast_species_indices.insert(m_baseEngine.getSpeciesIndex(product)); + } + } + } + return fast_reaction_indices; + } + + void MultiscalePartitioningEngineView::buildConnectivityGraph( + const std::unordered_set &fast_reaction_indices + ) { + const auto& all_reactions = m_baseEngine.getNetworkReactions(); + for (const size_t reaction_idx : fast_reaction_indices) { + const auto& reaction = all_reactions[reaction_idx]; + const auto& reactants = reaction.reactants(); + const auto& products = reaction.products(); + + // For each fast reaction, create edges between all reactants and all products. + // This represents that nucleons can flow quickly between these species. + for (const auto& reactant : reactants) { + const size_t reactant_idx = m_baseEngine.getSpeciesIndex(reactant); + for (const auto& product : products) { + const size_t product_idx = m_baseEngine.getSpeciesIndex(product); + + // Add a two-way edge to the adjacency list. + m_connectivity_graph[reactant_idx].push_back(product_idx); + m_connectivity_graph[product_idx].push_back(reactant_idx); + } + } + } + + + } + + void MultiscalePartitioningEngineView::findConnectedComponents() { + const size_t num_total_species = m_baseEngine.getNetworkSpecies().size(); + std::vector visited(num_total_species, false); + + for (size_t i = 0; i < num_total_species; ++i) { + if (!visited[i]) { + // Start a new BFS traversal from this unvisited node. + // This will find one complete connected component. + QSEGroup new_group; + std::queue q; + + q.push(i); + visited[i] = true; + + while (!q.empty()) { + const size_t current_species_idx = q.front(); + q.pop(); + new_group.species_indices.push_back(current_species_idx); + + // Check if the node exists in the graph (it might be an isolated species) + if (m_connectivity_graph.contains(current_species_idx)) { + for (const size_t neighbor_idx : m_connectivity_graph.at(current_species_idx)) { + if (!visited[neighbor_idx]) { + visited[neighbor_idx] = true; + q.push(neighbor_idx); + } + } + } + } + + // A "group" must contain more than one species to be a cluster. + // Isolated species are treated as dynamic by default. + if (new_group.species_indices.size() > 1) { + m_qse_groups.push_back(new_group); + } + } + } + + } + + void MultiscalePartitioningEngineView::validateGroupsWithFluxAnalysis( + const std::vector &Y, + const double T9, + const double rho + ) { + constexpr double FLUX_RATIO_THRESHOLD = 100; + for (auto& group : m_qse_groups) { + double internal_flux = 0.0; + double external_flux = 0.0; + + const std::unordered_set group_members( + group.species_indices.begin(), + group.species_indices.end() + ); + + for (const auto& reaction: m_baseEngine.getNetworkReactions()) { + const double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho)); + if (flow == 0.0) { + continue; // Skip reactions with zero flow + } + bool has_internal_reactant = false; + bool has_external_reactant = false; + + for (const auto& reactant : reaction.reactants()) { + if (group_members.contains(m_baseEngine.getSpeciesIndex(reactant))) { + has_internal_reactant = true; + } else { + has_external_reactant = true; + } + } + + bool has_internal_product = false; + bool has_external_product = false; + + for (const auto& product : reaction.products()) { + if (group_members.contains(m_baseEngine.getSpeciesIndex(product))) { + has_internal_product = true; + } else { + has_external_product = true; + } + } + + // Classify the reaction based on its participants + if ((has_internal_reactant && has_internal_product) && !(has_external_reactant || has_external_product)) { + internal_flux += flow; // Internal flux within the group + } else if ((has_internal_reactant || has_internal_product) && (has_external_reactant || has_external_product)) { + external_flux += flow; // External flux to/from the group + } + // otherwise the reaction is fully decoupled from the QSE group and can be ignored. + } + if (internal_flux > FLUX_RATIO_THRESHOLD * external_flux) { + group.is_in_equilibrium = true; // This group is in equilibrium if internal flux is significantly larger than external flux. + } else { + group.is_in_equilibrium = false; + } + } + } + + std::pair, std::vector> MultiscalePartitioningEngineView::identifyDynamicSpecies( + const std::vector& Y, + const std::vector& qse_groups, + const double T9, + const double rho + ) const { + // This set will contain the indices of species that are truly in qse (as opposed to those which might be in a qse group but are not in equilibrium). + std::unordered_set algebraic_qse_species_indices; + auto species_timescales = m_baseEngine.getSpeciesTimescales(Y, T9, rho); + + for (auto& group : qse_groups) { + if (!group.is_in_equilibrium || group.species_indices.empty()) continue; + + for (size_t species_idx : group.species_indices) { + const double timescale = species_timescales.at(m_baseEngine.getNetworkSpecies()[species_idx]); + + if (timescale < 1) { // If the timescale is less than 1 second, we consider it algebraic. + algebraic_qse_species_indices.insert(species_idx); + } + } + } + + std::vector dynamic_species; + std::unordered_set qse_species_indices_set; + std::vector dynamic_species_indices; + const auto& all_base_species = m_baseEngine.getNetworkSpecies(); + + for (size_t i = 0; i < all_base_species.size(); ++i) { + if (!algebraic_qse_species_indices.contains(i)) { + dynamic_species.push_back(all_base_species[i]); + dynamic_species_indices.push_back(i); + } + } + return {dynamic_species, dynamic_species_indices}; + } + + std::vector MultiscalePartitioningEngineView::solveQSEAbundances( + const std::vector &Y_full, + const double T9, + const double rho + ) { + auto Y_out = Y_full; + auto species_timescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + + for (const auto& qse_group : m_qse_groups) { + if (!qse_group.is_in_equilibrium || qse_group.species_indices.empty()) continue; + std::vector qse_solve_indices; + std::vector seed_indices; + + for (size_t species_idx : qse_group.species_indices) { + const double timescale = species_timescales.at(m_baseEngine.getNetworkSpecies()[species_idx]); + if (timescale < 1.0) { // If the timescale is less than 1 second, we consider it algebraic. + qse_solve_indices.push_back(species_idx); + } else { + seed_indices.push_back(species_idx); + } + } + + if (qse_solve_indices.empty()) continue; + + Eigen::VectorXd Y_scale(qse_solve_indices.size()); + Eigen::VectorXd v_initial(qse_solve_indices.size()); + for (size_t i = 0; i < qse_solve_indices.size(); ++i) { + constexpr double abundance_floor = 1.0e-15; + const double initial_abundance = Y_full[qse_solve_indices[i]]; + Y_scale(i) = std::max(initial_abundance, abundance_floor); + v_initial(i) = std::asinh(initial_abundance / Y_scale(i)); // Scale the initial abundances using asinh + } + + EigenFunctor functor(*this, qse_solve_indices, Y_full, T9, rho, Y_scale); + Eigen::LevenbergMarquardt lm(functor); + lm.parameters.ftol = 1.0e-10; + lm.parameters.xtol = 1.0e-10; + + Eigen::LevenbergMarquardtSpace::Status status = lm.minimize(v_initial); + + if (status <= 0 || status >= 4) { + std::stringstream msg; + msg << "QSE solver failed with status: " << status << " for group with seed "; + if (seed_indices.size() == 1) { + msg << "nucleus " << m_baseEngine.getNetworkSpecies()[seed_indices[0]].name(); + } else { + msg << "nuclei "; + // TODO: Refactor nice list printing into utility function somewhere + int counter = 0; + for (const auto& seed_idx : seed_indices) { + msg << m_baseEngine.getNetworkSpecies()[seed_idx].name(); + if (counter < seed_indices.size() - 2) { + msg << ", "; + } else if (counter == seed_indices.size() - 2) { + if (seed_indices.size() < 2) { + msg << " and "; + } else { + msg << ", and "; + } + } + ++counter; + } + } + throw std::runtime_error(msg.str()); + } + 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 { + 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 + + for (size_t i = 0; i < m_qse_solve_indices.size(); ++i) { + y_trial[m_qse_solve_indices[i]] = y_qse(i); + } + + const auto [dydt, energy] = m_view->calculateRHSAndEnergy(y_trial, m_T9, m_rho); + f_qse.resize(m_qse_solve_indices.size()); + for (size_t i = 0; i < m_qse_solve_indices.size(); ++i) { + f_qse(i) = dydt[m_qse_solve_indices[i]]; + } + return 0; // Success + } + + int MultiscalePartitioningEngineView::EigenFunctor::df(const InputType &v_qse, JacobianType &J_qse) const { + 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 + + for (size_t i = 0; i < m_qse_solve_indices.size(); ++i) { + y_trial[m_qse_solve_indices[i]] = y_qse(i); + } + + // TODO: Think about if the jacobian matrix should be mutable so that generateJacobianMatrix can be const + m_view->generateJacobianMatrix(y_trial, m_T9, m_rho); + + // TODO: Think very carefully about the indices here. + J_qse.resize(m_qse_solve_indices.size(), m_qse_solve_indices.size()); + for (size_t i = 0; i < m_qse_solve_indices.size(); ++i) { + for (size_t j = 0; j < m_qse_solve_indices.size(); ++j) { + J_qse(i, j) = m_view->getJacobianMatrixEntry( + m_qse_solve_indices[i], + m_qse_solve_indices[j] + ); + } + } + + // Chain rule for asinh scaling: + for (long j = 0; j < J_qse.cols(); ++j) { + const double dY_dv = m_Y_scale(j) * std::cosh(v_qse(j)); + J_qse.col(j) *= dY_dv; // Scale the column by the derivative of the asinh scaling + } + return 0; // Success + } + + +} diff --git a/src/network/lib/engine/views/engine_priming.cpp b/src/network/lib/engine/views/engine_priming.cpp new file mode 100644 index 00000000..1abc3a65 --- /dev/null +++ b/src/network/lib/engine/views/engine_priming.cpp @@ -0,0 +1,72 @@ +#include "gridfire/engine/views/engine_priming.h" +#include "gridfire/solver/solver.h" + +#include "fourdst/composition/species.h" +#include "fourdst/logging/logging.h" + +#include "quill/LogMacros.h" +#include "quill/Logger.h" + +#include +#include +#include +#include +#include +#include +#include +#include + + +namespace gridfire { + using fourdst::atomic::species; + + NetworkPrimingEngineView::NetworkPrimingEngineView( + const std::string &primingSymbol, + DynamicEngine &baseEngine + ) : + DefinedEngineView( + constructPrimingReactionSet( + species.at(primingSymbol), + baseEngine + ), + baseEngine + ), + m_primingSpecies(species.at(primingSymbol)) {} + + NetworkPrimingEngineView::NetworkPrimingEngineView( + const fourdst::atomic::Species &primingSpecies, + DynamicEngine &baseEngine + ) : + DefinedEngineView( + constructPrimingReactionSet( + primingSpecies, + baseEngine + ), + baseEngine + ), + m_primingSpecies(primingSpecies) {} + + + std::vector NetworkPrimingEngineView::constructPrimingReactionSet( + const fourdst::atomic::Species &primingSpecies, + const DynamicEngine &baseEngine + ) const { + std::unordered_set primeReactions; + for (const auto &reaction : baseEngine.getNetworkReactions()) { + if (reaction.contains(primingSpecies)) { + primeReactions.insert(std::string(reaction.peName())); + } + } + if (primeReactions.empty()) { + LOG_ERROR(m_logger, "No priming reactions found for species '{}'.", primingSpecies.name()); + m_logger->flush_log(); + throw std::runtime_error("No priming reactions found for species '" + std::string(primingSpecies.name()) + "'."); + } + std::vector primingReactionSet(primeReactions.begin(), primeReactions.end()); + // LOG_INFO(m_logger, "Constructed priming reaction set with {} reactions for species '{}'.", primingReactionSet.size(), primingSpecies.name()); + return primingReactionSet; + } + + + +} diff --git a/src/network/lib/io/network_file.cpp b/src/network/lib/io/network_file.cpp index 273a6b79..7e3e460f 100644 --- a/src/network/lib/io/network_file.cpp +++ b/src/network/lib/io/network_file.cpp @@ -68,9 +68,9 @@ namespace gridfire::io { if (line.empty()) { continue; // Skip empty lines } - parsed.reactionPENames.push_back(line); + parsed.push_back(line); } - LOG_TRACE_L1(m_logger, "Parsed {} reactions from file: {}", parsed.reactionPENames.size(), filename); + LOG_TRACE_L1(m_logger, "Parsed {} reactions from file: {}", parsed.size(), filename); return parsed; } diff --git a/src/network/lib/partition/composite/partition_composite.cpp b/src/network/lib/partition/composite/partition_composite.cpp index 7f00c7bd..9f02b267 100644 --- a/src/network/lib/partition/composite/partition_composite.cpp +++ b/src/network/lib/partition/composite/partition_composite.cpp @@ -25,13 +25,13 @@ namespace gridfire::partition { } double CompositePartitionFunction::evaluate(int z, int a, double T9) const { - LOG_TRACE_L1(m_logger, "Evaluating partition function for Z={} A={} T9={}", z, a, T9); + LOG_TRACE_L3(m_logger, "Evaluating partition function for Z={} A={} T9={}", z, a, T9); for (const auto& partitionFunction : m_partitionFunctions) { if (partitionFunction->supports(z, a)) { - LOG_TRACE_L2(m_logger, "Partition function of type {} supports Z={} A={}", partitionFunction->type(), z, a); + LOG_TRACE_L3(m_logger, "Partition function of type {} supports Z={} A={}", partitionFunction->type(), z, a); return partitionFunction->evaluate(z, a, T9); } else { - LOG_TRACE_L2(m_logger, "Partition function of type {} does not support Z={} A={}", partitionFunction->type(), z, a); + LOG_TRACE_L3(m_logger, "Partition function of type {} does not support Z={} A={}", partitionFunction->type(), z, a); } } LOG_ERROR( @@ -48,7 +48,7 @@ namespace gridfire::partition { double CompositePartitionFunction::evaluateDerivative(int z, int a, double T9) const { for (const auto& partitionFunction : m_partitionFunctions) { if (partitionFunction->supports(z, a)) { - LOG_TRACE_L2(m_logger, "Evaluating derivative of partition function for Z={} A={} T9={}", z, a, T9); + LOG_TRACE_L3(m_logger, "Evaluating derivative of partition function for Z={} A={} T9={}", z, a, T9); return partitionFunction->evaluateDerivative(z, a, T9); } } diff --git a/src/network/lib/solver/solver.cpp b/src/network/lib/solver/solver.cpp index 34021816..65701d7c 100644 --- a/src/network/lib/solver/solver.cpp +++ b/src/network/lib/solver/solver.cpp @@ -2,7 +2,9 @@ #include "gridfire/engine/engine_graph.h" #include "gridfire/network.h" + #include "gridfire/utils/logging.h" +#include "gridfire/utils/qse_rules.h" #include "fourdst/composition/atomicSpecies.h" #include "fourdst/composition/composition.h" @@ -22,6 +24,7 @@ #include "quill/LogMacros.h" namespace gridfire::solver { + double s_prevTimestep = 0.0; NetOut QSENetworkSolver::evaluate(const NetIn &netIn) { // --- Use the policy to decide whether to update the view --- @@ -149,6 +152,15 @@ namespace gridfire::solver { std::vectorQSESpeciesIndices; // Fast species that are in QSE std::unordered_map speciesTimescale = m_engine.getSpeciesTimescales(Y, T9, rho); + LOG_DEBUG( + m_logger, + "{}", + gridfire::utils::formatNuclearTimescaleLogString( + m_engine, + Y, + T9, + rho + )); for (size_t i = 0; i < m_engine.getNetworkSpecies().size(); ++i) { const auto& species = m_engine.getNetworkSpecies()[i]; @@ -164,9 +176,17 @@ namespace gridfire::solver { const double final_timescale = std::min(network_timescale, decay_timescale); - if (std::isinf(final_timescale) || abundance < abundanceCutoff || final_timescale <= timescaleCutoff) { + const bool isQSE = is_species_in_qse( + network_timescale, + decay_timescale, + abundance + ); + + if (isQSE) { + LOG_TRACE_L2(m_logger, "{} is in QSE based on rules in qse_rules.h", species.name()); QSESpeciesIndices.push_back(i); } else { + LOG_TRACE_L2(m_logger, "{} is dynamic based on rules in qse_rules.h", species.name()); dynamicSpeciesIndices.push_back(i); } } @@ -209,11 +229,37 @@ namespace gridfire::solver { const dynamicQSESpeciesIndices &indices ) const { LOG_TRACE_L1(m_logger, "Calculating steady state abundances for QSE species..."); + LOG_TRACE_L1( + m_logger, + "Initial QSE species abundances: {}", + [&]() -> std::string { + std::stringstream ss; + ss << std::scientific << std::setprecision(5); + int count = 0; + for (const auto& i : indices.QSESpeciesIndices) { + ss << std::string(m_engine.getNetworkSpecies()[i].name()) << ": " << Y[i] << ", "; + if (count < indices.QSESpeciesIndices.size() - 2) { + ss << ", "; + } else if (count == indices.QSESpeciesIndices.size() - 2) { + ss << " and "; + } + count++; + } + return ss.str(); + }()); if (indices.QSESpeciesIndices.empty()) { LOG_DEBUG(m_logger, "No QSE species to solve for."); return Eigen::VectorXd(0); } + + Eigen::VectorXd YScale(indices.QSESpeciesIndices.size()); + const double abundanceFloor = 1e-15; + for (size_t i = 0; i < indices.QSESpeciesIndices.size(); i++) { + const double initial_abundance = Y[indices.QSESpeciesIndices[i]]; + YScale(i) = std::max(initial_abundance, abundanceFloor); + } + // Use the EigenFunctor with Eigen's nonlinear solver EigenFunctor functor( m_engine, @@ -221,13 +267,31 @@ namespace gridfire::solver { indices.dynamicSpeciesIndices, indices.QSESpeciesIndices, T9, - rho + rho, + YScale ); - Eigen::VectorXd v_qse_log_initial(indices.QSESpeciesIndices.size()); + Eigen::VectorXd v_qse_initial(indices.QSESpeciesIndices.size()); for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) { - v_qse_log_initial(i) = std::log(std::max(Y[indices.QSESpeciesIndices[i]], 1e-99)); + const double initial_abundance = Y[indices.QSESpeciesIndices[i]]; + v_qse_initial(i) = std::asinh(initial_abundance/ YScale(i)); // Use asinh for better numerical stability compared to log } + LOG_TRACE_L1( + m_logger, + "v_qse_log_initial: {}", + [&]() -> std::string { + std::stringstream ss; + ss << std::scientific << std::setprecision(5); + for (size_t i = 0; i < v_qse_initial.size(); ++i) { + ss << "log(X_" << std::string(m_engine.getNetworkSpecies()[indices.QSESpeciesIndices[i]].name()) << ") = " << v_qse_initial(i); + if (i < v_qse_initial.size() - 2) { + ss << ", "; + } else if (i == v_qse_initial.size() - 2) { + ss << " and "; + } + } + return ss.str(); + }()); const static std::unordered_map statusMessages = { {Eigen::LevenbergMarquardtSpace::NotStarted, "Not started"}, @@ -245,7 +309,9 @@ namespace gridfire::solver { }; Eigen::LevenbergMarquardt lm(functor); - const Eigen::LevenbergMarquardtSpace::Status info = lm.minimize(v_qse_log_initial); + lm.parameters.xtol = 1e-24; + lm.parameters.ftol = 1e-24; + const Eigen::LevenbergMarquardtSpace::Status info = lm.minimize(v_qse_initial); if (info <= 0 || info >= 4) { LOG_ERROR(m_logger, "QSE species minimization failed with status: {} ({})", @@ -257,7 +323,7 @@ namespace gridfire::solver { } LOG_DEBUG(m_logger, "QSE species minimization completed successfully with status: {} ({})", static_cast(info), statusMessages.at(info)); - return v_qse_log_initial.array().exp(); + return YScale.array() * v_qse_initial.array().sinh(); // Convert back to molar abundances } @@ -287,6 +353,24 @@ namespace gridfire::solver { ignitionTime, ignitionStepSize ); + LOG_INFO( + m_logger, + "Pre-ignition composition: {}", + [&]() -> std::string { + std::stringstream ss; + ss << std::scientific << std::setprecision(3); + int count = 0; + for (const auto& entry : netIn.composition | std::views::values) { + ss << "X_" << std::string(entry.isotope().name()) << ": " << entry.mass_fraction() << " "; + if (count < netIn.composition.getRegisteredSymbols().size() - 2) { + ss << ", "; + } else if (count == netIn.composition.getRegisteredSymbols().size() - 2) { + ss << " and "; + } + count++; + } + return ss.str(); + }()); NetIn preIgnition = netIn; preIgnition.temperature = ignitionTemperature; @@ -300,6 +384,42 @@ namespace gridfire::solver { DirectNetworkSolver ignitionSolver(m_engine); NetOut postIgnition = ignitionSolver.evaluate(preIgnition); LOG_INFO(m_logger, "Network ignition completed in {} steps.", postIgnition.num_steps); + LOG_INFO( + m_logger, + "Post-ignition composition: {}", + [&]() -> std::string { + std::stringstream ss; + ss << std::scientific << std::setprecision(3); + int count = 0; + for (const auto& entry : postIgnition.composition | std::views::values) { + ss << "X_" << std::string(entry.isotope().name()) << ": " << entry.mass_fraction() << " "; + if (count < postIgnition.composition.getRegisteredSymbols().size() - 2) { + ss << ", "; + } else if (count == postIgnition.composition.getRegisteredSymbols().size() - 2) { + ss << " and "; + } + count++; + } + return ss.str(); + }()); + LOG_DEBUG( + m_logger, + "Average change in composition during ignition: {}", + [&]() -> std::string { + std::stringstream ss; + ss << std::scientific << std::setprecision(3); + for (const auto& entry : postIgnition.composition | std::views::values) { + if (!postIgnition.composition.contains(entry.isotope()) || !netIn.composition.contains(entry.isotope())) { + ss << std::string(entry.isotope().name()) << ": inf %, "; + continue; + } + const double initialAbundance = netIn.composition.getMolarAbundance(std::string(entry.isotope().name())); + const double finalAbundance = postIgnition.composition.getMolarAbundance(std::string(entry.isotope().name())); + const double change = (finalAbundance - initialAbundance) / initialAbundance * 100.0; // Percentage change + ss << std::string(entry.isotope().name()) << ": " << change << " %, "; + } + return ss.str(); + }()); m_engine.setScreeningModel(prevScreeningModel); LOG_DEBUG(m_logger, "Restoring previous screening model: {}", static_cast(prevScreeningModel)); return postIgnition; @@ -352,6 +472,15 @@ namespace gridfire::solver { boost::numeric::ublas::vector &dYdtDynamic, double t ) const { + LOG_TRACE_L1( + m_logger, + "Timestepping at t={:0.3E} (dt={:0.3E}) with T9={:0.3E}, ρ={:0.3E}", + t, + t - s_prevTimestep, + m_T9, + m_rho + ); + s_prevTimestep = t; // --- Populate the slow / dynamic species vector --- std::vector YFull(m_engine.getNetworkSpecies().size()); for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) { @@ -445,15 +574,6 @@ namespace gridfire::solver { double t ) const { const std::vector y(Y.begin(), m_numSpecies + Y.begin()); - - // std::string timescales = utils::formatNuclearTimescaleLogString( - // m_engine, - // y, - // m_T9, - // m_rho - // ); - // LOG_TRACE_L2(m_logger, "{}", timescales); - auto [dydt, eps] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho); dYdt.resize(m_numSpecies + 1); std::ranges::copy(dydt, dYdt.begin()); diff --git a/src/network/meson.build b/src/network/meson.build index 43648b6a..f31e8a21 100644 --- a/src/network/meson.build +++ b/src/network/meson.build @@ -5,6 +5,9 @@ network_sources = files( 'lib/engine/engine_graph.cpp', 'lib/engine/views/engine_adaptive.cpp', 'lib/engine/views/engine_defined.cpp', + 'lib/engine/views/engine_multiscale.cpp', + 'lib/engine/views/engine_priming.cpp', + 'lib/engine/procedures/priming.cpp', 'lib/reaction/reaction.cpp', 'lib/reaction/reaclib.cpp', 'lib/io/network_file.cpp', @@ -53,6 +56,9 @@ network_headers = files( 'include/gridfire/engine/engine_graph.h', 'include/gridfire/engine/views/engine_adaptive.h', 'include/gridfire/engine/views/engine_defined.h', + 'include/gridfire/engine/views/engine_multiscale.h', + 'include/gridfire/engine/views/engine_priming.h', + 'include/gridfire/engine/procedures/priming.h', 'include/gridfire/reaction/reaction.h', 'include/gridfire/reaction/reaclib.h', 'include/gridfire/io/network_file.h', @@ -66,5 +72,6 @@ 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') diff --git a/subprojects/libcomposition.wrap b/subprojects/libcomposition.wrap index 04bd9563..2fe7e9cf 100644 --- a/subprojects/libcomposition.wrap +++ b/subprojects/libcomposition.wrap @@ -1,4 +1,4 @@ [wrap-git] url = https://github.com/4D-STAR/libcomposition.git -revision = v1.2.0 +revision = v1.3.0 depth = 1 \ No newline at end of file diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp index 7116dfbe..175a947f 100644 --- a/tests/graphnet_sandbox/main.cpp +++ b/tests/graphnet_sandbox/main.cpp @@ -6,6 +6,8 @@ #include "gridfire/engine/views/engine_adaptive.h" #include "gridfire/partition/partition_types.h" #include "gridfire/engine/views/engine_defined.h" +#include "gridfire/engine/views/engine_multiscale.h" +#include "gridfire/engine/procedures/priming.h" #include "gridfire/io/network_file.h" #include "gridfire/solver/solver.h" @@ -58,17 +60,25 @@ void quill_terminate_handler() int main() { g_previousHandler = std::set_terminate(quill_terminate_handler); quill::Logger* logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - logger->set_log_level(quill::LogLevel::TraceL2); + logger->set_log_level(quill::LogLevel::TraceL1); LOG_DEBUG(logger, "Starting Adaptive Engine View Example..."); using namespace gridfire; const std::vector comp = {0.708, 0.0, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4}; const std::vector symbols = {"H-1", "H-2", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"}; + fourdst::composition::Composition composition; composition.registerSymbol(symbols, true); composition.setMassFraction(symbols, comp); composition.finalize(true); + using partition::BasePartitionType; + const auto partitionFunction = partition::CompositePartitionFunction({ + BasePartitionType::RauscherThielemann, + BasePartitionType::GroundState + }); + GraphEngine ReaclibEngine(composition, partitionFunction); + NetIn netIn; @@ -77,38 +87,30 @@ int main() { netIn.density = 1e2; netIn.energy = 0.0; - netIn.tMax = 3.1536e17; + netIn.tMax = 1; // 1 year in seconds netIn.dt0 = 1e12; NetOut netOut; - - // netIn.dt0 = 1e12; - // approx8::Approx8Network approx8Network; - // measure_execution_time([&]() { - // netOut = approx8Network.evaluate(netIn); - // }, "Approx8 Network Initialization"); - // std::cout << "Approx8 Network H-1: " << netOut.composition.getMassFraction("H-1") << " in " << netOut.num_steps << " steps." << std::endl; - - using partition::BasePartitionType; - const auto partitionFunction = partition::CompositePartitionFunction({ - BasePartitionType::RauscherThielemann, - BasePartitionType::GroundState - }); - std::cout << "Partition Function for Mg-24: " << partitionFunction.evaluate(12, 24, 8) << std::endl; - std::cout << "Partition Function for F-23: " << partitionFunction.evaluate(9, 23, 8) << std::endl; - std::cout << "Partition Function for O-13: " << partitionFunction.evaluate(8, 13, 8) << std::endl; - netIn.dt0 = 1e-15; - GraphEngine ReaclibEngine(composition, partitionFunction); - std::cout << ReaclibEngine.getPartitionFunction().type() << std::endl; + const auto primedComposition = primeNetwork(netIn, ReaclibEngine); + netIn.composition = primedComposition; + netIn.composition.finalize(true); + MultiscalePartitioningEngineView multiscaleEngine(ReaclibEngine); + multiscaleEngine.equilibrateNetwork(netIn, 1e12); + + + // // std::cout << ReaclibEngine.getPartitionFunction().type() << std::endl; // ReaclibEngine.setPrecomputation(true); - // // AdaptiveEngineView adaptiveEngine(ReaclibEngine); + // ReaclibEngine.setUseReverseReactions(false); + // // // AdaptiveEngineView adaptiveEngine(ReaclibEngine); // io::SimpleReactionListFileParser parser{}; // FileDefinedEngineView approx8EngineView(ReaclibEngine, "approx8.net", parser); // approx8EngineView.setScreeningModel(screening::ScreeningType::WEAK); - // solver::QSENetworkSolver solver(approx8EngineView); + // solver::DirectNetworkSolver solver(approx8EngineView); // netOut = solver.evaluate(netIn); + // std::cout << netOut.composition << std::endl; + // measure_execution_time([&]() { // netOut = solver.evaluate(netIn);