From 29af4c3bab6473a1b1dee2f659e2eeb699c962ce Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sun, 29 Jun 2025 14:53:39 -0400 Subject: [PATCH] feat(network): added half lifes, spin flip parity, better reaction acritecture --- .../include/gridfire/engine/engine_abstract.h | 166 ++++- .../include/gridfire/engine/engine_adaptive.h | 265 ++++++- .../include/gridfire/engine/engine_graph.h | 336 ++++++++- .../gridfire/engine/engine_view_abstract.h | 102 ++- src/network/include/gridfire/network.h | 50 +- .../include/gridfire/reaction/reaction.h | 687 +++++++++++++----- src/network/include/gridfire/solver/solver.h | 354 ++++++++- src/network/lib/engine/engine_adaptive.cpp | 264 ++++++- src/network/lib/engine/engine_graph.cpp | 83 ++- src/network/lib/network.cpp | 89 +-- src/network/lib/reaction/reaclib.cpp | 149 +++- src/network/lib/reaction/reaction.cpp | 224 +++--- src/network/lib/solver/solver.cpp | 132 ++-- src/network/meson.build | 6 +- 14 files changed, 2270 insertions(+), 637 deletions(-) diff --git a/src/network/include/gridfire/engine/engine_abstract.h b/src/network/include/gridfire/engine/engine_abstract.h index a62e7c88..9aed9628 100644 --- a/src/network/include/gridfire/engine/engine_abstract.h +++ b/src/network/include/gridfire/engine/engine_abstract.h @@ -1,30 +1,99 @@ #pragma once -#include "gridfire/network.h" // For NetIn, NetOut -#include "../reaction/reaction.h" - -#include "fourdst/composition/composition.h" -#include "fourdst/config/config.h" -#include "fourdst/logging/logging.h" +#include "gridfire/reaction/reaction.h" #include #include +/** + * @file engine_abstract.h + * @brief Abstract interfaces for reaction network engines in GridFire. + * + * This header defines the abstract base classes and concepts for implementing + * reaction network solvers in the GridFire framework. It provides the contract + * for calculating right-hand sides, energy generation, Jacobians, stoichiometry, + * and other core operations required for time integration of nuclear reaction networks. + * + * @author + * Emily M. Boudreaux + */ + namespace gridfire { + /** + * @brief Concept for types allowed in engine calculations. + * + * This concept restricts template parameters to either double or CppAD::AD, + * enabling both standard and automatic differentiation types. + */ template concept IsArithmeticOrAD = std::is_same_v || std::is_same_v>; + /** + * @brief Structure holding derivatives and energy generation for a network step. + * + * @tparam T Numeric type (double or CppAD::AD). + * + * This struct is used to return both the time derivatives of all species abundances + * and the specific nuclear energy generation rate for a single network evaluation. + * + * Example usage: + * @code + * StepDerivatives result = engine.calculateRHSAndEnergy(Y, T9, rho); + * for (double dydt_i : result.dydt) { + * // Use derivative + * } + * double energyRate = result.nuclearEnergyGenerationRate; + * @endcode + */ template struct StepDerivatives { - std::vector dydt; ///< Derivatives of abundances. - T nuclearEnergyGenerationRate = T(0.0); ///< Specific energy generation rate. + std::vector dydt; ///< Derivatives of abundances (dY/dt for each species). + T nuclearEnergyGenerationRate = T(0.0); ///< Specific energy generation rate (e.g., erg/g/s). }; + /** + * @brief Abstract base class for a reaction network engine. + * + * This class defines the minimal interface for a reaction network engine, + * which is responsible for evaluating the right-hand side (dY/dt) and + * energy generation for a given set of abundances, temperature, and density. + * + * Intended usage: Derive from this class to implement a concrete engine + * for a specific network or integration method. + * + * Example: + * @code + * class MyEngine : public gridfire::Engine { + * // Implement required methods... + * }; + * @endcode + */ class Engine { public: + /** + * @brief Virtual destructor. + */ virtual ~Engine() = default; + + /** + * @brief Get the list of species in the network. + * @return Vector of Species objects representing all network species. + */ virtual const std::vector& getNetworkSpecies() const = 0; + + /** + * @brief Calculate the right-hand side (dY/dt) and energy generation. + * + * @param Y Vector of current abundances for all species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return StepDerivatives containing dY/dt and energy generation rate. + * + * This function must be implemented by derived classes to compute the + * time derivatives of all species and the specific nuclear energy generation + * rate for the current state. + */ virtual StepDerivatives calculateRHSAndEnergy( const std::vector& Y, double T9, @@ -32,23 +101,85 @@ namespace gridfire { ) const = 0; }; + /** + * @brief Abstract class for engines supporting Jacobian and stoichiometry operations. + * + * Extends Engine with additional methods for: + * - Generating and accessing the Jacobian matrix (for implicit solvers). + * - Generating and accessing the stoichiometry matrix. + * - Calculating molar reaction flows for individual reactions. + * - Accessing the set of logical reactions in the network. + * - Computing timescales for each species. + * + * Intended usage: Derive from this class to implement engines that support + * advanced solver features such as implicit integration, sensitivity analysis, + * QSE (Quasi-Steady-State Equilibrium) handling, and more. + */ class DynamicEngine : public Engine { public: + /** + * @brief Generate the Jacobian matrix for the current state. + * + * @param Y Vector of current abundances. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * + * This method must compute and store the Jacobian matrix (∂(dY/dt)_i/∂Y_j) + * for the current state. The matrix can then be accessed via getJacobianMatrixEntry(). + */ virtual void generateJacobianMatrix( const std::vector& Y, double T9, double rho ) = 0; + + /** + * @brief Get an entry from the previously generated Jacobian matrix. + * + * @param i Row index (species index). + * @param j Column index (species index). + * @return Value of the Jacobian matrix at (i, j). + * + * The Jacobian must have been generated by generateJacobianMatrix() before calling this. + */ virtual double getJacobianMatrixEntry( int i, int j ) const = 0; + /** + * @brief Generate the stoichiometry matrix for the network. + * + * This method must compute and store the stoichiometry matrix, + * which encodes the net change of each species in each reaction. + */ virtual void generateStoichiometryMatrix() = 0; + + /** + * @brief Get an entry from the stoichiometry matrix. + * + * @param speciesIndex Index of the species. + * @param reactionIndex Index of the reaction. + * @return Stoichiometric coefficient for the species in the reaction. + * + * The stoichiometry matrix must have been generated by generateStoichiometryMatrix(). + */ virtual int getStoichiometryMatrixEntry( int speciesIndex, int reactionIndex ) const = 0; + /** + * @brief Calculate the molar reaction flow for a given reaction. + * + * @param reaction The reaction for which to calculate the flow. + * @param Y Vector of current abundances. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return Molar flow rate for the reaction (e.g., mol/g/s). + * + * This method computes the net rate at which the given reaction proceeds + * under the current state. + */ virtual double calculateMolarReactionFlow( const reaction::Reaction& reaction, const std::vector& Y, @@ -56,7 +187,24 @@ namespace gridfire { double rho ) const = 0; - virtual const reaction::REACLIBLogicalReactionSet& getNetworkReactions() const = 0; + /** + * @brief Get the set of logical reactions in the network. + * + * @return Reference to the LogicalReactionSet containing all reactions. + */ + virtual const reaction::LogicalReactionSet& getNetworkReactions() const = 0; + + /** + * @brief Compute timescales for all species in the network. + * + * @param Y Vector of current abundances. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return Map from Species to their characteristic timescales (s). + * + * This method estimates the timescale for abundance change of each species, + * which can be used for timestep control, diagnostics, and reaction network culling. + */ virtual std::unordered_map getSpeciesTimescales( const std::vector& Y, double T9, diff --git a/src/network/include/gridfire/engine/engine_adaptive.h b/src/network/include/gridfire/engine/engine_adaptive.h index 0d6c4bde..974a43a7 100644 --- a/src/network/include/gridfire/engine/engine_adaptive.h +++ b/src/network/include/gridfire/engine/engine_adaptive.h @@ -1,75 +1,316 @@ #pragma once #include "gridfire/engine/engine_abstract.h" #include "gridfire/engine/engine_view_abstract.h" +#include "gridfire/network.h" #include "fourdst/composition/atomicSpecies.h" +#include "fourdst/config/config.h" #include "fourdst/logging/logging.h" #include "quill/Logger.h" namespace gridfire { + /** + * @class AdaptiveEngineView + * @brief An engine view that dynamically adapts the reaction network based on runtime conditions. + * + * This class implements an EngineView that dynamically culls species and reactions from the + * full reaction network based on their reaction flow rates and connectivity. This allows for + * efficient simulation of reaction networks by focusing computational effort on the most + * important species and reactions. + * + * The AdaptiveEngineView maintains a subset of "active" species and reactions, and maps + * between the full network indices and the active subset indices. This allows the base engine + * to operate on the full network data, while the AdaptiveEngineView provides a reduced view + * for external clients. + * + * The adaptation process is driven by the `update()` method, which performs the following steps: + * 1. **Reaction Flow Calculation:** Calculates the molar reaction flow rate for each reaction + * in the full network based on the current temperature, density, and composition. + * 2. **Reaction Culling:** Culls reactions with flow rates below a threshold, determined by + * a relative culling threshold multiplied by the maximum flow rate. + * 3. **Connectivity Analysis:** Performs a connectivity analysis to identify species that are + * reachable from the initial fuel species through the culled reaction network. + * 4. **Species Culling:** Culls species that are not reachable from the initial fuel. + * 5. **Index Map Construction:** Constructs index maps to map between the full network indices + * and the active subset indices for species and reactions. + * + * @implements DynamicEngine + * @implements EngineView + * + * @see engine_abstract.h + * @see engine_view_abstract.h + * @see AdaptiveEngineView::update() + */ class AdaptiveEngineView final : public DynamicEngine, public EngineView { public: + /** + * @brief Constructs an AdaptiveEngineView. + * + * @param baseEngine The underlying DynamicEngine to which this view delegates calculations. + * + * Initializes the active species and reactions to the full network, and constructs the + * initial index maps. + */ explicit AdaptiveEngineView(DynamicEngine& baseEngine); + /** + * @brief Updates the active species and reactions based on the current conditions. + * + * @param netIn The current network input, containing temperature, density, and composition. + * + * This method performs the reaction flow calculation, reaction culling, connectivity analysis, + * and index map construction steps described above. + * + * The culling thresholds are read from the configuration using the following keys: + * - `gridfire:AdaptiveEngineView:RelativeCullingThreshold` (default: 1e-75) + * + * @throws std::runtime_error If there is a mismatch between the active reactions and the base engine. + * @post The active species and reactions are updated, and the index maps are reconstructed. + * @see AdaptiveEngineView + * @see AdaptiveEngineView::constructSpeciesIndexMap() + * @see AdaptiveEngineView::constructReactionIndexMap() + */ void update(const NetIn& netIn); + /** + * @brief Gets the list of active species in the network. + * @return A const reference to the vector of active species. + */ const std::vector& getNetworkSpecies() const override; + + /** + * @brief Calculates the right-hand side (dY/dt) and energy generation for the active species. + * + * @param Y_culled 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. + * @return A StepDerivatives struct containing the derivatives of the active species and the + * nuclear energy generation rate. + * + * This method maps the culled abundances to the full network abundances, calls the base engine + * to calculate the RHS and energy generation, and then maps the full network derivatives back + * to the culled derivatives. + * + * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). + * @see AdaptiveEngineView::update() + */ StepDerivatives calculateRHSAndEnergy( - const std::vector &Y, + const std::vector &Y_culled, const double T9, const double rho ) const override; + /** + * @brief Generates the Jacobian matrix for the active species. + * + * @param Y_culled 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. + * + * This method maps the culled abundances to the full network abundances and calls the base engine + * to generate the Jacobian matrix. + * + * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). + * @see AdaptiveEngineView::update() + */ void generateJacobianMatrix( - const std::vector &Y, + const std::vector &Y_culled, const double T9, const double rho ) override; + + /** + * @brief Gets an entry from the Jacobian matrix for the active species. + * + * @param i_culled The row index (species index) in the culled matrix. + * @param j_culled The column index (species index) in the culled matrix. + * @return The value of the Jacobian matrix at (i_culled, j_culled). + * + * This method maps the culled indices to the full network indices and calls the base engine + * to get the Jacobian matrix entry. + * + * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). + * @throws std::out_of_range If the culled index is out of bounds for the species index map. + * @see AdaptiveEngineView::update() + */ double getJacobianMatrixEntry( - const int i, - const int j + const int i_culled, + const int j_culled ) const override; + /** + * @brief Generates the stoichiometry matrix for the active reactions and species. + * + * This method calls the base engine to generate the stoichiometry matrix. + * + * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). + * @note The stoichiometry matrix generated by the base engine is assumed to be consistent with + * the active species and reactions in this view. + */ void generateStoichiometryMatrix() override; + + /** + * @brief Gets an entry from the stoichiometry matrix for the active species and reactions. + * + * @param speciesIndex_culled The index of the species in the culled species list. + * @param reactionIndex_culled The index of the reaction in the culled reaction list. + * @return The stoichiometric coefficient for the given species and reaction. + * + * This method maps the culled indices to the full network indices and calls the base engine + * to get the stoichiometry matrix entry. + * + * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). + * @throws std::out_of_range If the culled index is out of bounds for the species or reaction index map. + * @see AdaptiveEngineView::update() + */ int getStoichiometryMatrixEntry( - const int speciesIndex, - const int reactionIndex + const int speciesIndex_culled, + const int reactionIndex_culled ) const override; + /** + * @brief Calculates the molar reaction flow for a given reaction in the active network. + * + * @param reaction The reaction for which to calculate the flow. + * @param Y_culled Vector of current abundances for the active species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return Molar flow rate for the reaction (e.g., mol/g/s). + * + * This method maps the culled abundances to the full network abundances and calls the base engine + * to calculate the molar reaction flow. + * + * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). + * @throws std::runtime_error If the reaction is not part of the active reactions in the adaptive engine view. + */ double calculateMolarReactionFlow( const reaction::Reaction &reaction, - const std::vector &Y, + const std::vector &Y_culled, double T9, double rho ) const override; - const reaction::REACLIBLogicalReactionSet& getNetworkReactions() const override; + /** + * @brief Gets the set of active logical reactions in the network. + * + * @return Reference to the LogicalReactionSet containing all active reactions. + */ + const reaction::LogicalReactionSet& getNetworkReactions() const override; + + /** + * @brief Computes timescales for all active species in the network. + * + * @param Y_culled Vector of current abundances for the active species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return Map from Species to their characteristic timescales (s). + * + * This method maps the culled abundances to the full network abundances and calls the base engine + * to compute the species timescales. + * + * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). + */ std::unordered_map getSpeciesTimescales( - const std::vector &Y, + const std::vector &Y_culled, double T9, double rho ) const override; + /** + * @brief Gets the base engine. + * @return A const reference to the base engine. + */ const DynamicEngine& getBaseEngine() const override { return m_baseEngine; } private: using Config = fourdst::config::Config; using LogManager = fourdst::logging::LogManager; + Config& m_config = Config::getInstance(); + quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); DynamicEngine& m_baseEngine; std::vector m_activeSpecies; - reaction::REACLIBLogicalReactionSet m_activeReactions; + reaction::LogicalReactionSet m_activeReactions; std::vector m_speciesIndexMap; + std::vector m_reactionIndexMap; bool m_isStale = true; - Config& m_config = Config::getInstance(); - quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); private: + /** + * @brief Constructs the species index map. + * + * @return A vector mapping culled species indices to full species indices. + * + * This method creates a map from the indices of the active species to the indices of the + * corresponding species in the full network. + * + * @see AdaptiveEngineView::update() + */ std::vector constructSpeciesIndexMap() const; + + /** + * @brief Constructs the reaction index map. + * + * @return A vector mapping culled reaction indices to full reaction indices. + * + * This method creates a map from the indices of the active reactions to the indices of the + * corresponding reactions in the full network. + * + * @see AdaptiveEngineView::update() + */ + std::vector constructReactionIndexMap() const; + + /** + * @brief Maps a vector of culled abundances to a vector of full abundances. + * + * @param culled A vector of abundances for the active species. + * @return A vector of abundances for the full network, with the abundances of the active + * species copied from the culled vector. + */ + std::vector mapCulledToFull(const std::vector& culled) const; + + /** + * @brief Maps a vector of full abundances to a vector of culled abundances. + * + * @param full A vector of abundances for the full network. + * @return A vector of abundances for the active species, with the abundances of the active + * species copied from the full vector. + */ + std::vector mapFullToCulled(const std::vector& full) const; + + /** + * @brief Maps a culled species index to a full species index. + * + * @param culledSpeciesIndex The index of the species in the culled species list. + * @return The index of the corresponding species in the full network. + * + * @throws std::out_of_range If the culled index is out of bounds for the species index map. + */ + size_t mapCulledToFullSpeciesIndex(size_t culledSpeciesIndex) const; + + /** + * @brief Maps a culled reaction index to a full reaction index. + * + * @param culledReactionIndex The index of the reaction in the culled reaction list. + * @return The index of the corresponding reaction in the full network. + * + * @throws std::out_of_range If the culled index is out of bounds for the reaction index map. + */ + size_t mapCulledToFullReactionIndex(size_t culledReactionIndex) const; + + /** + * @brief Validates that the AdaptiveEngineView is not stale. + * + * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). + */ + void validateState() const; private: + /** + * @brief A struct to hold a reaction and its flow rate. + */ struct ReactionFlow { const reaction::Reaction* reactionPtr; double flowRate; diff --git a/src/network/include/gridfire/engine/engine_graph.h b/src/network/include/gridfire/engine/engine_graph.h index ac444186..0c712e41 100644 --- a/src/network/include/gridfire/engine/engine_graph.h +++ b/src/network/include/gridfire/engine/engine_graph.h @@ -23,35 +23,148 @@ // REACLIBReactions are quite large data structures, so this could be a performance bottleneck. namespace gridfire { - typedef CppAD::AD ADDouble; ///< Alias for CppAD AD type for double precision. + /** + * @brief Alias for CppAD AD type for double precision. + * + * This alias simplifies the use of the CppAD automatic differentiation type. + */ + typedef CppAD::AD ADDouble; using fourdst::config::Config; using fourdst::logging::LogManager; using fourdst::constant::Constants; + /** + * @brief Minimum density threshold below which reactions are ignored. + * + * Reactions are not calculated if the density falls below this threshold. + * This helps to improve performance by avoiding unnecessary calculations + * in very low-density regimes. + */ static constexpr double MIN_DENSITY_THRESHOLD = 1e-18; + + /** + * @brief Minimum abundance threshold below which species are ignored. + * + * Species with abundances below this threshold are treated as zero in + * reaction rate calculations. This helps to improve performance by + * avoiding unnecessary calculations for trace species. + */ static constexpr double MIN_ABUNDANCE_THRESHOLD = 1e-18; + + /** + * @brief Minimum value for Jacobian matrix entries. + * + * Jacobian matrix entries with absolute values below this threshold are + * treated as zero to maintain sparsity and improve performance. + */ static constexpr double MIN_JACOBIAN_THRESHOLD = 1e-24; + /** + * @class GraphEngine + * @brief A reaction network engine that uses a graph-based representation. + * + * The GraphEngine class implements the DynamicEngine interface using a + * graph-based representation of the reaction network. It uses sparse + * matrices for efficient storage and computation of the stoichiometry + * and Jacobian matrices. Automatic differentiation (AD) is used to + * calculate the Jacobian matrix. + * + * The engine supports: + * - Calculation of the right-hand side (dY/dt) and energy generation rate. + * - Generation and access to the Jacobian matrix. + * - Generation and access to the stoichiometry matrix. + * - Calculation of molar reaction flows. + * - Access to the set of logical reactions in the network. + * - Computation of timescales for each species. + * - Exporting the network to DOT and CSV formats for visualization and analysis. + * + * @implements DynamicEngine + * + * @see engine_abstract.h + */ class GraphEngine final : public DynamicEngine{ public: + /** + * @brief Constructs a GraphEngine from a composition. + * + * @param composition The composition of the material. + * + * This constructor builds the reaction network from the given composition + * using the `build_reaclib_nuclear_network` function. + * + * @see build_reaclib_nuclear_network + */ explicit GraphEngine(const fourdst::composition::Composition &composition); - explicit GraphEngine(reaction::REACLIBLogicalReactionSet reactions); + /** + * @brief Constructs a GraphEngine from a set of reactions. + * + * @param reactions The set of reactions to use in the network. + * + * This constructor uses the given set of reactions to construct the + * reaction network. + */ + explicit GraphEngine(reaction::LogicalReactionSet reactions); + + /** + * @brief Calculates the right-hand side (dY/dt) and energy generation rate. + * + * @param Y Vector of current abundances for all species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return StepDerivatives containing dY/dt and energy generation rate. + * + * This method calculates the time derivatives of all species and the + * specific nuclear energy generation rate for the current state. + * + * @see StepDerivatives + */ StepDerivatives calculateRHSAndEnergy( const std::vector& Y, const double T9, const double rho ) const override; + /** + * @brief Generates the Jacobian matrix for the current state. + * + * @param Y Vector of current abundances. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * + * This method computes and stores the Jacobian matrix (∂(dY/dt)_i/∂Y_j) + * for the current state using automatic differentiation. The matrix can + * then be accessed via `getJacobianMatrixEntry()`. + * + * @see getJacobianMatrixEntry() + */ void generateJacobianMatrix( const std::vector& Y, const double T9, const double rho ) override; + /** + * @brief Generates the stoichiometry matrix for the network. + * + * This method computes and stores the stoichiometry matrix, + * which encodes the net change of each species in each reaction. + */ void generateStoichiometryMatrix() override; + /** + * @brief Calculates the molar reaction flow for a given reaction. + * + * @param reaction The reaction for which to calculate the flow. + * @param Y Vector of current abundances. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return Molar flow rate for the reaction (e.g., mol/g/s). + * + * This method computes the net rate at which the given reaction proceeds + * under the current state. + */ double calculateMolarReactionFlow( const reaction::Reaction& reaction, const std::vector&Y, @@ -59,39 +172,130 @@ namespace gridfire { const double rho ) const override; + /** + * @brief Gets the list of species in the network. + * @return Vector of Species objects representing all network species. + */ [[nodiscard]] const std::vector& getNetworkSpecies() const override; - [[nodiscard]] const reaction::REACLIBLogicalReactionSet& getNetworkReactions() const override; + + /** + * @brief Gets the set of logical reactions in the network. + * @return Reference to the LogicalReactionSet containing all reactions. + */ + [[nodiscard]] const reaction::LogicalReactionSet& getNetworkReactions() const override; + + /** + * @brief Gets an entry from the previously generated Jacobian matrix. + * + * @param i Row index (species index). + * @param j Column index (species index). + * @return Value of the Jacobian matrix at (i, j). + * + * The Jacobian must have been generated by `generateJacobianMatrix()` before calling this. + * + * @see generateJacobianMatrix() + */ [[nodiscard]] double getJacobianMatrixEntry( const int i, const int j ) const override; - [[nodiscard]] std::unordered_map getNetReactionStoichiometry( + + /** + * @brief Gets the net stoichiometry for a given reaction. + * + * @param reaction The reaction for which to get the stoichiometry. + * @return Map of species to their stoichiometric coefficients. + */ + [[nodiscard]] static std::unordered_map getNetReactionStoichiometry( const reaction::Reaction& reaction - ) const; + ); + + /** + * @brief Gets an entry from the stoichiometry matrix. + * + * @param speciesIndex Index of the species. + * @param reactionIndex Index of the reaction. + * @return Stoichiometric coefficient for the species in the reaction. + * + * The stoichiometry matrix must have been generated by `generateStoichiometryMatrix()`. + * + * @see generateStoichiometryMatrix() + */ [[nodiscard]] int getStoichiometryMatrixEntry( const int speciesIndex, const int reactionIndex ) const override; + + /** + * @brief Computes timescales for all species in the network. + * + * @param Y Vector of current abundances. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return Map from Species to their characteristic timescales (s). + * + * This method estimates the timescale for abundance change of each species, + * which can be used for timestep control or diagnostics. + */ [[nodiscard]] std::unordered_map getSpeciesTimescales( const std::vector& Y, double T9, double rho ) const override; + /** + * @brief Checks if a given species is involved in the network. + * + * @param species The species to check. + * @return True if the species is involved in the network, false otherwise. + */ [[nodiscard]] bool involvesSpecies( const fourdst::atomic::Species& species ) const; + /** + * @brief Exports the network to a DOT file for visualization. + * + * @param filename The name of the DOT file to create. + * + * This method generates a DOT file that can be used to visualize the + * reaction network as a graph. The DOT file can be converted to a + * graphical image using Graphviz. + * + * @throws std::runtime_error If the file cannot be opened for writing. + * + * Example usage: + * @code + * engine.exportToDot("network.dot"); + * @endcode + */ void exportToDot( const std::string& filename ) const; + + /** + * @brief Exports the network to a CSV file for analysis. + * + * @param filename The name of the CSV file to create. + * + * This method generates a CSV file containing information about the + * reactions in the network, including the reactants, products, Q-value, + * and reaction rate coefficients. + * + * @throws std::runtime_error If the file cannot be opened for writing. + * + * Example usage: + * @code + * engine.exportToCSV("network.csv"); + * @endcode + */ void exportToCSV( const std::string& filename ) const; private: - reaction::REACLIBLogicalReactionSet m_reactions; ///< Set of REACLIB reactions in the network. + reaction::LogicalReactionSet m_reactions; ///< Set of REACLIB reactions in the network. std::unordered_map m_reactionIDMap; ///< Map from reaction ID to REACLIBReaction. //PERF: This makes copies of REACLIBReaction and could be a performance bottleneck. std::vector m_networkSpecies; ///< Vector of unique species in the network. @@ -108,20 +312,100 @@ namespace gridfire { quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); private: + /** + * @brief Synchronizes the internal maps. + * + * This method synchronizes the internal maps used by the engine, + * including the species map, reaction ID map, and species-to-index map. + * It also generates the stoichiometry matrix and records the AD tape. + */ void syncInternalMaps(); + + /** + * @brief Collects the unique species in the network. + * + * This method collects the unique species in the network from the + * reactants and products of all reactions. + */ void collectNetworkSpecies(); + + /** + * @brief Populates the reaction ID map. + * + * This method populates the reaction ID map, which maps reaction IDs + * to REACLIBReaction objects. + */ void populateReactionIDMap(); + + /** + * @brief Populates the species-to-index map. + * + * This method populates the species-to-index map, which maps species + * to their index in the stoichiometry matrix. + */ void populateSpeciesToIndexMap(); + + /** + * @brief Reserves space for the Jacobian matrix. + * + * This method reserves space for the Jacobian matrix, which is used + * to store the partial derivatives of the right-hand side of the ODE + * with respect to the species abundances. + */ void reserveJacobianMatrix(); + + /** + * @brief Records the AD tape for the right-hand side of the ODE. + * + * This method records the AD tape for the right-hand side of the ODE, + * which is used to calculate the Jacobian matrix using automatic + * differentiation. + * + * @throws std::runtime_error If there are no species in the network. + */ void recordADTape(); + /** + * @brief Validates mass and charge conservation across all reactions. + * + * @return True if all reactions conserve mass and charge, false otherwise. + * + * This method checks that all reactions in the network conserve mass + * and charge. If any reaction does not conserve mass or charge, an + * error message is logged and false is returned. + */ [[nodiscard]] bool validateConservation() const; + + /** + * @brief Validates the composition against the current reaction set. + * + * @param composition The composition to validate. + * @param culling The culling threshold to use. + * @param T9 The temperature to use. + * + * This method validates the composition against the current reaction set. + * If the composition is not compatible with the reaction set, the + * reaction set is rebuilt from the composition. + */ void validateComposition( const fourdst::composition::Composition &composition, double culling, double T9 ); + /** + * @brief Calculates the molar reaction flow for a given reaction. + * + * @tparam T The numeric type to use for the calculation. + * @param reaction The reaction for which to calculate the flow. + * @param Y Vector of current abundances. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return Molar flow rate for the reaction (e.g., mol/g/s). + * + * This method computes the net rate at which the given reaction proceeds + * under the current state. + */ template T calculateMolarReactionFlow( const reaction::Reaction &reaction, @@ -130,6 +414,18 @@ namespace gridfire { const T rho ) const; + /** + * @brief Calculates all derivatives (dY/dt) and the energy generation rate. + * + * @tparam T The numeric type to use for the calculation. + * @param Y_in Vector of current abundances for all species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return StepDerivatives containing dY/dt and energy generation rate. + * + * This method calculates the time derivatives of all species and the + * specific nuclear energy generation rate for the current state. + */ template StepDerivatives calculateAllDerivatives( const std::vector &Y_in, @@ -137,16 +433,40 @@ namespace gridfire { T rho ) const; + /** + * @brief Calculates all derivatives (dY/dt) and the energy generation rate (double precision). + * + * @param Y_in Vector of current abundances for all species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return StepDerivatives containing dY/dt and energy generation rate. + * + * This method calculates the time derivatives of all species and the + * specific nuclear energy generation rate for the current state using + * double precision arithmetic. + */ StepDerivatives calculateAllDerivatives( const std::vector& Y_in, const double T9, const double rho ) const; + /** + * @brief Calculates all derivatives (dY/dt) and the energy generation rate (automatic differentiation). + * + * @param Y_in Vector of current abundances for all species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return StepDerivatives containing dY/dt and energy generation rate. + * + * This method calculates the time derivatives of all species and the + * specific nuclear energy generation rate for the current state using + * automatic differentiation. + */ StepDerivatives calculateAllDerivatives( const std::vector& Y_in, - const ADDouble T9, - const ADDouble rho + const ADDouble &T9, + const ADDouble &rho ) const; }; diff --git a/src/network/include/gridfire/engine/engine_view_abstract.h b/src/network/include/gridfire/engine/engine_view_abstract.h index c754fe00..6397d369 100644 --- a/src/network/include/gridfire/engine/engine_view_abstract.h +++ b/src/network/include/gridfire/engine/engine_view_abstract.h @@ -1,8 +1,98 @@ -// -// Created by Emily Boudreaux on 6/27/25. -// +#pragma once -#ifndef ENGINE_VIEW_ABSTRACT_H -#define ENGINE_VIEW_ABSTRACT_H +#include "gridfire/engine/engine_abstract.h" -#endif //ENGINE_VIEW_ABSTRACT_H +/** + * @file engine_view_abstract.h + * @brief Abstract interfaces for engine "views" in GridFire. + * + * This header defines the abstract base classes and concepts for "views" of reaction network engines. + * The primary purpose of an EngineView is to enable dynamic or adaptive network topologies + * (such as species/reaction culling, masking, or remapping) without modifying the underlying + * physics engine or its implementation. Engine views act as a flexible interface layer, + * allowing the network structure to be changed at runtime while preserving the core + * physics and solver logic in the base engine. + * + * Typical use cases include: + * - Adaptive or reduced networks for computational efficiency. + * - Dynamic masking or culling of species/reactions based on runtime criteria. + * - Providing a filtered or remapped view of the network for integration or analysis. + * - Supporting dynamic network reconfiguration in multi-zone or multi-physics contexts. + * + * The base engine types referenced here are defined in @ref engine_abstract.h. + * + * @author + * Emily M. Boudreaux + */ + +namespace gridfire { + + /** + * @brief Concept for types allowed as engine bases in EngineView. + * + * This concept restricts template parameters to types derived from either + * gridfire::Engine or gridfire::DynamicEngine, as defined in engine_abstract.h. + * + * Example usage: + * @code + * static_assert(EngineType); + * @endcode + */ + template + concept EngineType = std::is_base_of_v || std::is_base_of_v; + + /** + * @brief Abstract base class for a "view" of a reaction network engine. + * + * @tparam EngineT The engine type being viewed (must satisfy EngineType). + * + * EngineView provides an interface for accessing an underlying engine instance, + * while presenting a potentially modified or reduced network structure to the user. + * This enables dynamic or adaptive network topologies (e.g., culling, masking, or + * remapping of species and reactions) without altering the core physics engine. + * + * Intended usage: Derive from this class to implement a custom view or wrapper + * that manages a dynamic or adaptive network structure, delegating core calculations + * to the base engine. The contract is that getBaseEngine() must return a reference + * to the underlying engine instance, which remains responsible for the full physics. + * + * Example (see also AdaptiveEngineView): + * @code + * class MyAdaptiveView : public gridfire::EngineView { + * public: + * MyAdaptiveView(DynamicEngine& engine) : engine_(engine) {} + * const DynamicEngine& getBaseEngine() const override { return engine_; } + * // Implement dynamic masking/culling logic... + * private: + * DynamicEngine& engine_; + * }; + * @endcode + * + * @see gridfire::AdaptiveEngineView for a concrete example of dynamic culling. + */ + template + class EngineView { + public: + /** + * @brief Virtual destructor. + */ + virtual ~EngineView() = default; + + /** + * @brief Access the underlying engine instance. + * + * @return Const reference to the underlying engine. + * + * This method must be implemented by derived classes to provide access + * to the base engine. The returned reference should remain valid for the + * lifetime of the EngineView. + * + * Example: + * @code + * const DynamicEngine& engine = myView.getBaseEngine(); + * @endcode + */ + virtual const EngineT& getBaseEngine() const = 0; + }; + +} \ No newline at end of file diff --git a/src/network/include/gridfire/network.h b/src/network/include/gridfire/network.h index a9d56283..74c87e65 100644 --- a/src/network/include/gridfire/network.h +++ b/src/network/include/gridfire/network.h @@ -50,23 +50,6 @@ namespace gridfire { {UNKNOWN, "Unknown"} }; - /** - * @struct NetIn - * @brief Input structure for the network evaluation. - * - * This structure holds the input parameters required for the network evaluation. - * - * Example usage: - * @code - * nuclearNetwork::NetIn netIn; - * netIn.composition = {1.0, 0.0, 0.0}; - * netIn.tmax = 1.0e6; - * netIn.dt0 = 1.0e-3; - * netIn.temperature = 1.0e8; - * netIn.density = 1.0e5; - * netIn.energy = 1.0e12; - * @endcode - */ struct NetIn { fourdst::composition::Composition composition; ///< Composition of the network double tMax; ///< Maximum time @@ -75,22 +58,10 @@ namespace gridfire { double density; ///< Density in g/cm^3 double energy; ///< Energy in ergs double culling = 0.0; ///< Culling threshold for reactions (default is 0.0, meaning no culling) + + std::vector MolarAbundance() const; }; - /** - * @struct NetOut - * @brief Output structure for the network evaluation. - * - * This structure holds the output results from the network evaluation. - * - * Example usage: - * @code - * nuclearNetwork::NetOut netOut = network.evaluate(netIn); - * std::vector composition = netOut.composition; - * int steps = netOut.num_steps; - * double energy = netOut.energy; - * @endcode - */ struct NetOut { fourdst::composition::Composition composition; ///< Composition of the network after evaluation int num_steps; ///< Number of steps taken in the evaluation @@ -102,20 +73,6 @@ namespace gridfire { } }; - /** - * @class Network - * @brief Class for network evaluation. - * - * This class provides methods to evaluate the network based on the input parameters. - * - * Example usage: - * @code - * nuclearNetwork::Network network; - * nuclearNetwork::NetIn netIn; - * // Set netIn parameters... - * nuclearNetwork::NetOut netOut = network.evaluate(netIn); - * @endcode - */ class Network { public: explicit Network(const NetworkFormat format = NetworkFormat::APPROX8); @@ -147,8 +104,7 @@ namespace gridfire { }; - reaction::REACLIBLogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, bool reverse); - reaction::REACLIBLogicalReactionSet build_reaclib_nuclear_network_from_file(const std::string& filename, bool reverse); + reaction::LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, bool reverse); } // namespace nuclearNetwork diff --git a/src/network/include/gridfire/reaction/reaction.h b/src/network/include/gridfire/reaction/reaction.h index 250b32f9..6ed7a6db 100644 --- a/src/network/include/gridfire/reaction/reaction.h +++ b/src/network/include/gridfire/reaction/reaction.h @@ -7,154 +7,87 @@ #include "quill/Logger.h" #include #include -#include #include -#include #include "cppad/cppad.hpp" +/** + * @file reaction.h + * @brief Defines classes for representing and managing nuclear reactions. + * + * This file contains the core data structures for handling nuclear reactions, + * including individual reactions from specific sources (`Reaction`), collections + * of reactions (`ReactionSet`), and logical reactions that aggregate rates from +* multiple sources (`LogicalReaction`, `LogicalReactionSet`). + */ namespace gridfire::reaction { /** - * @struct REACLIBRateCoefficientSet - * @brief Holds the seven fitting parameters for a single REACLIB rate set. - * @details The thermonuclear reaction rate for a single set is calculated as: - * rate = exp(a0 + a1/T9 + a2/T9^(-1/3) + a3*T9^(1/3) + a4*T9 + a5*T9^(5/3) + a6*ln(T9)) - * where T9 is the temperature in billions of Kelvin. The total rate for a - * reaction is the sum of the rates from all its sets. + * @struct RateCoefficientSet + * @brief Holds the seven coefficients for the REACLIB rate equation. + * + * This struct stores the parameters (a0-a6) used to calculate reaction rates + * as a function of temperature. */ + struct RateCoefficientSet { + double a0; ///< Coefficient a0 + double a1; ///< Coefficient a1 + double a2; ///< Coefficient a2 + double a3; ///< Coefficient a3 + double a4; ///< Coefficient a4 + double a5; ///< Coefficient a5 + double a6; ///< Coefficient a6 - struct REACLIBRateCoefficientSet { - const double a0; - const double a1; - const double a2; - const double a3; - const double a4; - const double a5; - const double a6; - - friend std::ostream& operator<<(std::ostream& os, const REACLIBRateCoefficientSet& r) { + /** + * @brief Overloads the stream insertion operator for easy printing. + * @param os The output stream. + * @param r The RateCoefficientSet to print. + * @return The output stream. + */ + friend std::ostream& operator<<(std::ostream& os, const RateCoefficientSet& r) { os << "[" << r.a0 << ", " << r.a1 << ", " << r.a2 << ", " << r.a3 << ", " << r.a4 << ", " << r.a5 << ", " << r.a6 << "]"; return os; } }; + /** + * @class Reaction + * @brief Represents a single nuclear reaction from a specific data source. + * + * This class encapsulates all properties of a single nuclear reaction as defined + * in formats like REACLIB, including reactants, products, Q-value, and rate + * coefficients from a particular evaluation (source). + * + * Example: + * @code + * // Assuming species and rate coefficients are defined + * Reaction p_gamma_d( + * "H_1_H_1_to_H_2", "p(p,g)d", 1, {H_1, H_1}, {H_2}, 5.493, "st08", rate_coeffs + * ); + * double rate = p_gamma_d.calculate_rate(0.1); // T9 = 0.1 + * @endcode + */ class Reaction { public: - Reaction( - const std::string_view id, - const double qValue, - const std::vector& reactants, - const std::vector& products, - const bool reverse = false - ); + /** + * @brief Virtual destructor. + */ virtual ~Reaction() = default; - virtual std::unique_ptr clone() const = 0; - - virtual double calculate_rate(double T9) const = 0; - virtual CppAD::AD calculate_rate(const CppAD::AD T9) const = 0; - - virtual std::string_view peName() const { return ""; } - - [[nodiscard]] bool contains(const fourdst::atomic::Species& species) const; - [[nodiscard]] bool contains_reactant(const fourdst::atomic::Species& species) const; - [[nodiscard]] bool contains_product(const fourdst::atomic::Species& species) const; - - std::unordered_set all_species() const; - std::unordered_set reactant_species() const; - std::unordered_set product_species() const; - - size_t num_species() const; - - int stoichiometry(const fourdst::atomic::Species& species) const; - std::unordered_map stoichiometry() const; - - std::string_view id() const { return m_id; } - double qValue() const { return m_qValue; } - const std::vector& reactants() const { return m_reactants; } - const std::vector& products() const { return m_products; } - bool is_reverse() const { return m_reverse; } - - double excess_energy() const; - - bool operator==(const Reaction& other) const { return m_id == other.m_id; } - bool operator!=(const Reaction& other) const { return !(*this == other); } - [[nodiscard]] uint64_t hash(uint64_t seed) const; - - protected: - quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - std::string m_id; - double m_qValue = 0.0; ///< Q-value of the reaction - std::vector m_reactants; ///< Reactants of the reaction - std::vector m_products; ///< Products of the reaction - bool m_reverse = false; - }; - - class ReactionSet { - public: - explicit ReactionSet(std::vector> reactions); - ReactionSet(const ReactionSet& other); - ReactionSet& operator=(const ReactionSet& other); - virtual ~ReactionSet() = default; - - virtual void add_reaction(std::unique_ptr reaction); - virtual void remove_reaction(const std::unique_ptr& reaction); - - bool contains(const std::string_view& id) const; - bool contains(const Reaction& reaction) const; - - size_t size() const { return m_reactions.size(); } - - void sort(double T9=1.0); - - bool contains_species(const fourdst::atomic::Species& species) const; - bool contains_reactant(const fourdst::atomic::Species& species) const; - bool contains_product(const fourdst::atomic::Species& species) const; - - [[nodiscard]] const Reaction& operator[](size_t index) const; - [[nodiscard]] const Reaction& operator[](const std::string_view& id) const; - - bool operator==(const ReactionSet& other) const; - bool operator!=(const ReactionSet& other) const; - - [[nodiscard]] uint64_t hash(uint64_t seed = 0) const; - - auto begin() { return m_reactions.begin(); } - auto begin() const { return m_reactions.cbegin(); } - auto end() { return m_reactions.end(); } - auto end() const { return m_reactions.cend(); } - protected: - quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - std::vector> m_reactions; - std::string m_id; - std::unordered_map m_reactionNameMap; ///< Maps reaction IDs to Reaction objects for quick lookup - - }; - - /** - * @struct REACLIBReaction - * @brief Represents a single nuclear reaction from the JINA REACLIB database. - * @details This struct is designed to be constructed at compile time (constexpr) from - * the data parsed by the Python generation script. It stores all necessary - * information to identify a reaction and calculate its rate. - */ - class REACLIBReaction final : public Reaction { - public: /** - * @brief Constructs a REACLIBReaction object at compile time. - * @param id A unique string identifier generated by the Python script. - * @param peName - * @param chapter The REACLIB chapter number, defining the reaction structure. - * @param reactants A vector of strings with the names of the reactant species. - * @param products A vector of strings with the names of the product species. + * @brief Constructs a Reaction object. + * @param id A unique identifier for the reaction. + * @param peName The name in (projectile, ejectile) notation (e.g., "p(p,g)d"). + * @param chapter The REACLIB chapter number, defining reaction structure. + * @param reactants A vector of reactant species. + * @param products A vector of product species. * @param qValue The Q-value of the reaction in MeV. - * @param label The source label for the rate data (e.g., "wc12w", "st08"). - * @param sets A vector of RateFitSet, containing the fitting coefficients for the rate. - * @param reverse A boolean indicating if the reaction is reversed (default is false). + * @param label The source label for the rate data (e.g., "wc12", "st08"). + * @param sets The set of rate coefficients. + * @param reverse True if this is a reverse reaction rate. */ - REACLIBReaction( + Reaction( const std::string_view id, const std::string_view peName, const int chapter, @@ -162,67 +95,433 @@ namespace gridfire::reaction { const std::vector &products, const double qValue, const std::string_view label, - const REACLIBRateCoefficientSet &sets, + const RateCoefficientSet &sets, const bool reverse = false); - [[nodiscard]] std::unique_ptr clone() const override; + /** + * @brief Calculates the reaction rate for a given temperature. + * @param T9 The temperature in units of 10^9 K. + * @return The calculated reaction rate. + */ + [[nodiscard]] virtual double calculate_rate(const double T9) const; - [[nodiscard]] double calculate_rate(const double T9) const override; - [[nodiscard]] CppAD::AD calculate_rate(const CppAD::AD T9) const override; + /** + * @brief Calculates the reaction rate for a given temperature using CppAD types. + * @param T9 The temperature in units of 10^9 K, as a CppAD::AD. + * @return The calculated reaction rate, as a CppAD::AD. + */ + [[nodiscard]] virtual CppAD::AD calculate_rate(const CppAD::AD T9) const; - template - [[nodiscard]] GeneralScalarType calculate_rate(const GeneralScalarType T9) const { - const GeneralScalarType T913 = CppAD::pow(T9, 1.0/3.0); - const GeneralScalarType rateExponent = m_rateCoefficients.a0 + - m_rateCoefficients.a1 / T9 + - m_rateCoefficients.a2 / T913 + - m_rateCoefficients.a3 * T913 + - m_rateCoefficients.a4 * T9 + - m_rateCoefficients.a5 * CppAD::pow(T9, 5.0/3.0) + - m_rateCoefficients.a6 * CppAD::log(T9); - return CppAD::exp(rateExponent); - } - - [[nodiscard]] std::string_view peName() const override { return m_peName; } + /** + * @brief Gets the reaction name in (projectile, ejectile) notation. + * @return The reaction name (e.g., "p(p,g)d"). + */ + [[nodiscard]] virtual std::string_view peName() const { return m_peName; } + /** + * @brief Gets the REACLIB chapter number. + * @return The chapter number. + */ [[nodiscard]] int chapter() const { return m_chapter; } + /** + * @brief Gets the source label for the rate data. + * @return The source label (e.g., "wc12w", "st08"). + */ [[nodiscard]] std::string_view sourceLabel() const { return m_sourceLabel; } - [[nodiscard]] const REACLIBRateCoefficientSet& rateCoefficients() const { return m_rateCoefficients; } + /** + * @brief Gets the set of rate coefficients. + * @return A const reference to the RateCoefficientSet. + */ + [[nodiscard]] const RateCoefficientSet& rateCoefficients() const { return m_rateCoefficients; } - friend std::ostream& operator<<(std::ostream& os, const REACLIBReaction& reaction); - private: - std::string m_peName; ///< Name of the reaction in (projectile, ejectile) notation (e.g. p(p, g)d) + /** + * @brief Checks if the reaction involves a given species as a reactant or product. + * @param species The species to check for. + * @return True if the species is involved, false otherwise. + */ + [[nodiscard]] bool contains(const fourdst::atomic::Species& species) const; + + /** + * @brief Checks if the reaction involves a given species as a reactant. + * @param species The species to check for. + * @return True if the species is a reactant, false otherwise. + */ + [[nodiscard]] bool contains_reactant(const fourdst::atomic::Species& species) const; + + /** + * @brief Checks if the reaction involves a given species as a product. + * @param species The species to check for. + * @return True if the species is a product, false otherwise. + */ + [[nodiscard]] bool contains_product(const fourdst::atomic::Species& species) const; + + /** + * @brief Gets a set of all unique species involved in the reaction. + * @return An unordered_set of all reactant and product species. + */ + [[nodiscard]] std::unordered_set all_species() const; + + /** + * @brief Gets a set of all unique reactant species. + * @return An unordered_set of reactant species. + */ + [[nodiscard]] std::unordered_set reactant_species() const; + + /** + * @brief Gets a set of all unique product species. + * @return An unordered_set of product species. + */ + [[nodiscard]] std::unordered_set product_species() const; + + /** + * @brief Gets the number of unique species involved in the reaction. + * @return The count of unique species. + */ + [[nodiscard]] size_t num_species() const; + + /** + * @brief Calculates the stoichiometric coefficient for a given species. + * @param species The species for which to find the coefficient. + * @return The stoichiometric coefficient (negative for reactants, positive for products). + */ + [[nodiscard]] int stoichiometry(const fourdst::atomic::Species& species) const; + + /** + * @brief Gets a map of all species to their stoichiometric coefficients. + * @return An unordered_map from species to their integer coefficients. + */ + [[nodiscard]] std::unordered_map stoichiometry() const; + + /** + * @brief Gets the unique identifier of the reaction. + * @return The reaction ID. + */ + [[nodiscard]] std::string_view id() const { return m_id; } + + /** + * @brief Gets the Q-value of the reaction. + * @return The Q-value in whatever units the reaction was defined in (usually MeV). + */ + [[nodiscard]] double qValue() const { return m_qValue; } + + /** + * @brief Gets the vector of reactant species. + * @return A const reference to the vector of reactants. + */ + [[nodiscard]] const std::vector& reactants() const { return m_reactants; } + + /** + * @brief Gets the vector of product species. + * @return A const reference to the vector of products. + */ + [[nodiscard]] const std::vector& products() const { return m_products; } + + /** + * @brief Checks if this is a reverse reaction rate. + * @return True if it is a reverse rate, false otherwise. + */ + [[nodiscard]] bool is_reverse() const { return m_reverse; } + + /** + * @brief Calculates the excess energy from the mass difference of reactants and products. + * @return The excess energy in MeV. + */ + [[nodiscard]] double excess_energy() const; + + /** + * @brief Compares this reaction with another for equality based on their IDs. + * @param other The other Reaction to compare with. + * @return True if the reaction IDs are the same. + */ + bool operator==(const Reaction& other) const { return m_id == other.m_id; } + + /** + * @brief Compares this reaction with another for inequality. + * @param other The other Reaction to compare with. + * @return True if the reactions are not equal. + */ + bool operator!=(const Reaction& other) const { return !(*this == other); } + + /** + * @brief Computes a hash for the reaction based on its ID. + * @param seed The seed for the hash function. + * @return A 64-bit hash value. + * @details Uses the XXHash64 algorithm on the reaction's ID string. + */ + [[nodiscard]] uint64_t hash(uint64_t seed = 0) const; + + protected: + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + std::string m_id; ///< Unique identifier for the reaction (e.g., "h1+h1=>h2+e+nu"). + std::string m_peName; ///< Name of the reaction in (projectile, ejectile) notation (e.g. "p(p,g)d"). int m_chapter; ///< Chapter number from the REACLIB database, defining the reaction structure. - std::string m_sourceLabel; ///< Source label for the rate data, indicating the origin of the rate coefficients (e.g., "wc12w", "st08"). - REACLIBRateCoefficientSet m_rateCoefficients; + double m_qValue = 0.0; ///< Q-value of the reaction in MeV. + std::vector m_reactants; ///< Reactants of the reaction. + std::vector m_products; ///< Products of the reaction. + std::string m_sourceLabel; ///< Source label for the rate data (e.g., "wc12w", "st08"). + RateCoefficientSet m_rateCoefficients; ///< The seven rate coefficients. + bool m_reverse = false; ///< Flag indicating if this is a reverse reaction rate. + private: + /** + * @brief Template implementation for calculating the reaction rate. + * @tparam T The numeric type (double or CppAD::AD). + * @param T9 The temperature in units of 10^9 K. + * @return The calculated reaction rate. + * @details The rate is calculated using the standard REACLIB formula: + * `rate = exp(a0 + a1/T9 + a2/T9^(1/3) + a3*T9^(1/3) + a4*T9 + a5*T9^(5/3) + a6*ln(T9))` + */ + template + [[nodiscard]] T calculate_rate(const T T9) const { + const T T913 = CppAD::pow(T9, 1.0/3.0); + const T T953 = CppAD::pow(T9, 5.0/3.0); + const T logT9 = CppAD::log(T9); + const T exponent = m_rateCoefficients.a0 + + m_rateCoefficients.a1 / T9 + + m_rateCoefficients.a2 / T913 + + m_rateCoefficients.a3 * T913 + + m_rateCoefficients.a4 * T9 + + m_rateCoefficients.a5 * T953 + + m_rateCoefficients.a6 * logT9; + return CppAD::exp(exponent); + } }; - class REACLIBReactionSet final : public ReactionSet { + /** + * @class ReactionSet + * @brief A collection of Reaction objects. + * + * This class manages a set of individual `Reaction` objects, providing + * efficient lookup by ID and functionality to query the entire set. + * + * Example: + * @code + * ReactionSet my_set({reaction1, reaction2}); + * my_set.add_reaction(reaction3); + * if (my_set.contains("h1(p,g)h2")) { + * const Reaction& r = my_set["h1(p,g)h2"]; + * } + * @endcode + */ + class ReactionSet { public: - explicit REACLIBReactionSet(std::vector); - std::unordered_set peNames() const; - friend std::ostream& operator<<(std::ostream& os, const REACLIBReactionSet& set); + /** + * @brief Constructs a ReactionSet from a vector of reactions. + * @param reactions The initial vector of Reaction objects. + */ + explicit ReactionSet(std::vector reactions); + + /** + * @brief Copy constructor. + * @param other The ReactionSet to copy. + */ + ReactionSet(const ReactionSet& other); + + /** + * @brief Copy assignment operator. + * @param other The ReactionSet to assign from. + * @return A reference to this ReactionSet. + */ + ReactionSet& operator=(const ReactionSet& other); + + /** + * @brief Virtual destructor. + */ + virtual ~ReactionSet() = default; + + /** + * @brief Adds a reaction to the set. + * @param reaction The Reaction to add. + */ + virtual void add_reaction(Reaction reaction); + + /** + * @brief Removes a reaction from the set. + * @param reaction The Reaction to remove. + */ + virtual void remove_reaction(const Reaction& reaction); + + /** + * @brief Checks if the set contains a reaction with the given ID. + * @param id The ID of the reaction to find. + * @return True if the reaction is in the set, false otherwise. + */ + [[nodiscard]] bool contains(const std::string_view& id) const; + + /** + * @brief Checks if the set contains the given reaction. + * @param reaction The Reaction to find. + * @return True if the reaction is in the set, false otherwise. + */ + [[nodiscard]] bool contains(const Reaction& reaction) const; + + /** + * @brief Gets the number of reactions in the set. + * @return The size of the set. + */ + [[nodiscard]] virtual size_t size() const { return m_reactions.size(); } + + /** + * @brief Removes all reactions from the set. + */ + void clear(); + + /** + * @brief Checks if any reaction in the set involves the given species. + * @param species The species to check for. + * @return True if the species is involved in any reaction. + */ + [[nodiscard]] bool contains_species(const fourdst::atomic::Species& species) const; + + /** + * @brief Checks if any reaction in the set contains the given species as a reactant. + * @param species The species to check for. + * @return True if the species is a reactant in any reaction. + */ + [[nodiscard]] bool contains_reactant(const fourdst::atomic::Species& species) const; + + /** + * @brief Checks if any reaction in the set contains the given species as a product. + * @param species The species to check for. + * @return True if the species is a product in any reaction. + */ + [[nodiscard]] bool contains_product(const fourdst::atomic::Species& species) const; + + /** + * @brief Accesses a reaction by its index. + * @param index The index of the reaction to access. + * @return A const reference to the Reaction. + * @throws std::out_of_range if the index is out of bounds. + */ + [[nodiscard]] virtual const Reaction& operator[](size_t index) const; + + /** + * @brief Accesses a reaction by its ID. + * @param id The ID of the reaction to access. + * @return A const reference to the Reaction. + * @throws std::out_of_range if no reaction with the given ID exists. + */ + [[nodiscard]] const Reaction& operator[](const std::string_view& id) const; + + /** + * @brief Compares this set with another for equality. + * @param other The other ReactionSet to compare with. + * @return True if the sets are equal (same size and hash). + */ + bool operator==(const ReactionSet& other) const; + + /** + * @brief Compares this set with another for inequality. + * @param other The other ReactionSet to compare with. + * @return True if the sets are not equal. + */ + bool operator!=(const ReactionSet& other) const; + + /** + * @brief Computes a hash for the entire set. + * @param seed The seed for the hash function. + * @return A 64-bit hash value. + * @details The algorithm computes the hash of each individual reaction, + * sorts the hashes, and then computes a final hash over the sorted list + * of hashes. This ensures the hash is order-independent. + */ + [[nodiscard]] uint64_t hash(uint64_t seed = 0) const; + + /** @name Iterators + * Provides iterators to loop over the reactions in the set. + */ + ///@{ + auto begin() { return m_reactions.begin(); } + [[nodiscard]] auto begin() const { return m_reactions.cbegin(); } + auto end() { return m_reactions.end(); } + [[nodiscard]] auto end() const { return m_reactions.cend(); } + ///@} + private: + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + std::vector m_reactions; + std::string m_id; + std::unordered_map m_reactionNameMap; ///< Maps reaction IDs to Reaction objects for quick lookup. + }; - class REACLIBLogicalReaction final : public Reaction { + + /** + * @class LogicalReaction + * @brief Represents a "logical" reaction that aggregates rates from multiple sources. + * + * A LogicalReaction shares the same reactants and products but combines rates + * from different evaluations (e.g., "wc12" and "st08" for the same physical + * reaction). The total rate is the sum of the individual rates. + * It inherits from Reaction, using the properties of the first provided reaction + * as its base properties (reactants, products, Q-value, etc.). + */ + class LogicalReaction final : public Reaction { public: - explicit REACLIBLogicalReaction(const std::vector &reactions); - explicit REACLIBLogicalReaction(const REACLIBReaction &reaction); - void add_reaction(const REACLIBReaction& reaction); + /** + * @brief Constructs a LogicalReaction from a vector of `Reaction` objects. + * @param reactions A vector of reactions that represent the same logical process. + * @throws std::runtime_error if the provided reactions have inconsistent Q-values. + */ + explicit LogicalReaction(const std::vector &reactions); - [[nodiscard]] std::unique_ptr clone() const override; + /** + * @brief Adds another `Reaction` source to this logical reaction. + * @param reaction The reaction to add. + * @throws std::runtime_error if the reaction has a different `peName`, a duplicate + * source label, or an inconsistent Q-value. + */ + void add_reaction(const Reaction& reaction); - [[nodiscard]] std::string_view peName() const override { return m_id; }; + /** + * @brief Gets the number of source rates contributing to this logical reaction. + * @return The number of aggregated rates. + */ [[nodiscard]] size_t size() const { return m_rates.size(); } + /** + * @brief Gets the list of source labels for the aggregated rates. + * @return A vector of source label strings. + */ [[nodiscard]] std::vector sources() const { return m_sources; } + /** + * @brief Calculates the total reaction rate by summing all source rates. + * @param T9 The temperature in units of 10^9 K. + * @return The total calculated reaction rate. + */ [[nodiscard]] double calculate_rate(const double T9) const override; + /** + * @brief Calculates the total reaction rate using CppAD types. + * @param T9 The temperature in units of 10^9 K, as a CppAD::AD. + * @return The total calculated reaction rate, as a CppAD::AD. + */ [[nodiscard]] CppAD::AD calculate_rate(const CppAD::AD T9) const override; + /** @name Iterators + * Provides iterators to loop over the rate coefficient sets. + */ + ///@{ + auto begin() { return m_rates.begin(); } + [[nodiscard]] auto begin() const { return m_rates.cbegin(); } + auto end() { return m_rates.end(); } + [[nodiscard]] auto end() const { return m_rates.cend(); } + ///@} + + private: + std::vector m_sources; ///< List of source labels. + std::vector m_rates; ///< List of rate coefficient sets from each source. + + private: + /** + * @brief Template implementation for calculating the total reaction rate. + * @tparam T The numeric type (double or CppAD::AD). + * @param T9 The temperature in units of 10^9 K. + * @return The total calculated reaction rate. + * @details This method iterates through all stored `RateCoefficientSet`s, + * calculates the rate for each, and returns their sum. + */ template [[nodiscard]] T calculate_rate(const T T9) const { T sum = static_cast(0.0); @@ -242,29 +541,59 @@ namespace gridfire::reaction { } return sum; } - - [[nodiscard]] int chapter() const { return m_chapter; } - - auto begin() { return m_rates.begin(); } - auto begin() const { return m_rates.cbegin(); } - auto end() { return m_rates.end(); } - auto end() const { return m_rates.cend(); } - - private: - int m_chapter; - std::vector m_sources; - std::vector m_rates; - }; - class REACLIBLogicalReactionSet final : public ReactionSet { + /** + * @class LogicalReactionSet + * @brief A collection of LogicalReaction objects. + * + * This class takes a `ReactionSet` and groups individual `Reaction` objects + * into `LogicalReaction` objects based on their `peName`. This provides a + * view of the network where all rates for the same physical process are combined. + */ + class LogicalReactionSet final : public ReactionSet { public: - REACLIBLogicalReactionSet() = delete; - explicit REACLIBLogicalReactionSet(const REACLIBReactionSet& reactionSet); + /** + * @brief Deleted default constructor. + */ + LogicalReactionSet() = delete; - [[nodiscard]] std::unordered_set peNames() const; + /** + * @brief Constructs a LogicalReactionSet from a ReactionSet. + * @param reactionSet The set of individual reactions to group. + * @details This constructor iterates through the provided `ReactionSet`, + * groups reactions by their `peName`, and creates a `LogicalReaction` for each group. + */ + explicit LogicalReactionSet(const ReactionSet& reactionSet); + + /** @name Iterators + * Provides iterators to loop over the logical reactions in the set. + */ + ///@{ + auto begin() { return m_reactions.begin(); } + [[nodiscard]] auto begin() const { return m_reactions.cbegin(); } + auto end() { return m_reactions.end(); } + [[nodiscard]] auto end() const { return m_reactions.cend(); } + ///@} + + /** + * @brief Gets the number of logical reactions in the set. + * @return The size of the set. + */ + [[nodiscard]] size_t size() const { return m_reactions.size(); } + + /** + * @brief Accesses a logical reaction by its index. + * @param index The index of the logical reaction. + * @return A const reference to the LogicalReaction. + */ + [[nodiscard]] const LogicalReaction& operator[](size_t index) const { return m_reactions[index]; } private: - std::unordered_set m_peNames; + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + std::vector m_reactions; + std::string m_id; + std::unordered_map m_reactionNameMap; ///< Maps reaction IDs to LogicalReaction objects for quick lookup. }; -} \ No newline at end of file +} + diff --git a/src/network/include/gridfire/solver/solver.h b/src/network/include/gridfire/solver/solver.h index 50181e20..86ff74bb 100644 --- a/src/network/include/gridfire/solver/solver.h +++ b/src/network/include/gridfire/solver/solver.h @@ -2,6 +2,7 @@ #include "gridfire/engine/engine_graph.h" #include "gridfire/engine/engine_abstract.h" +#include "gridfire/engine/engine_adaptive.h" #include "gridfire/network.h" #include "fourdst/logging/logging.h" @@ -15,53 +16,207 @@ namespace gridfire::solver { + /** + * @struct dynamicQSESpeciesIndices + * @brief Structure to hold indices of dynamic and QSE species. + * + * This structure is used by the QSENetworkSolver to store the indices of species + * that are treated dynamically and those that are assumed to be in Quasi-Steady-State + * Equilibrium (QSE). + */ struct dynamicQSESpeciesIndices { - std::vector dynamicSpeciesIndices; // Slow species that are not in QSE - std::vector QSESpeciesIndices; // Fast species that are in QSE + std::vector dynamicSpeciesIndices; ///< Indices of slow species that are not in QSE. + std::vector QSESpeciesIndices; ///< Indices of fast species that are in QSE. }; + /** + * @class NetworkSolverStrategy + * @brief Abstract base class for network solver strategies. + * + * This class defines the interface for network solver strategies, which are responsible + * for integrating the reaction network ODEs over a given timestep. It is templated on the + * engine type to allow for different engine implementations to be used with the same solver. + * + * @tparam EngineT The type of engine to use with this solver strategy. Must inherit from Engine. + */ template class NetworkSolverStrategy { public: + /** + * @brief Constructor for the NetworkSolverStrategy. + * @param engine The engine to use for evaluating the network. + */ explicit NetworkSolverStrategy(EngineT& engine) : m_engine(engine) {}; + + /** + * @brief Virtual destructor. + */ virtual ~NetworkSolverStrategy() = default; + + /** + * @brief Evaluates the network for a given timestep. + * @param netIn The input conditions for the network. + * @return The output conditions after the timestep. + */ virtual NetOut evaluate(const NetIn& netIn) = 0; protected: - EngineT& m_engine; + EngineT& m_engine; ///< The engine used by this solver strategy. }; + /** + * @brief Type alias for a network solver strategy that uses a DynamicEngine. + */ using DynamicNetworkSolverStrategy = NetworkSolverStrategy; + + /** + * @brief Type alias for a network solver strategy that uses an AdaptiveEngineView. + */ + using AdaptiveNetworkSolverStrategy = NetworkSolverStrategy; + + /** + * @brief Type alias for a network solver strategy that uses a static Engine. + */ using StaticNetworkSolverStrategy = NetworkSolverStrategy; - - class QSENetworkSolver final : public DynamicNetworkSolverStrategy { + /** + * @class QSENetworkSolver + * @brief A network solver that uses a Quasi-Steady-State Equilibrium (QSE) approach. + * + * This solver partitions the network into "fast" species in QSE and "slow" (dynamic) species. + * The abundances of the fast species are determined by solving a system of algebraic + * equations, while the abundances of the slow species are integrated using an ODE solver. + * This hybrid approach is highly effective for stiff networks with disparate timescales. + * + * The QSE solver uses an AdaptiveEngineView to dynamically cull unimportant species and + * reactions, which significantly improves performance for large networks. + * + * @implements AdaptiveNetworkSolverStrategy + * + * @see AdaptiveEngineView + * @see DynamicEngine::getSpeciesTimescales() + */ + class QSENetworkSolver final : public AdaptiveNetworkSolverStrategy { public: - using DynamicNetworkSolverStrategy::DynamicNetworkSolverStrategy; + /** + * @brief Constructor for the QSENetworkSolver. + * @param engine The adaptive engine view to use for evaluating the network. + */ + using AdaptiveNetworkSolverStrategy::AdaptiveNetworkSolverStrategy; + + /** + * @brief Evaluates the network for a given timestep using the QSE approach. + * @param netIn The input conditions for the network. + * @return The output conditions after the timestep. + * + * This method performs the following steps: + * 1. Updates the adaptive engine view (if necessary). + * 2. Partitions the species into dynamic and QSE species based on their timescales. + * 3. Calculates the steady-state abundances of the QSE species. + * 4. Integrates the ODEs for the dynamic species using a Runge-Kutta solver. + * 5. Marshals the output variables into a NetOut struct. + * + * @throws std::runtime_error If the steady-state abundances cannot be calculated. + * + * @see AdaptiveEngineView::update() + * @see packSpeciesTypeIndexVectors() + * @see calculateSteadyStateAbundances() + */ NetOut evaluate(const NetIn& netIn) override; private: // methods + /** + * @brief Packs the species indices into vectors based on their type (dynamic or QSE). + * @param Y Vector of current abundances for all species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @return A dynamicQSESpeciesIndices struct containing the indices of the dynamic and QSE species. + * + * This method determines whether each species should be treated dynamically or as + * being in QSE based on its timescale and abundance. Species with short timescales + * or low abundances are assumed to be in QSE. + * + * @see DynamicEngine::getSpeciesTimescales() + */ dynamicQSESpeciesIndices packSpeciesTypeIndexVectors( const std::vector& Y, const double T9, const double rho ) const; + + /** + * @brief Calculates the steady-state abundances of the QSE species. + * @param Y Vector of current abundances for all species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + * @param indices A dynamicQSESpeciesIndices struct containing the indices of the dynamic and QSE species. + * @return An Eigen::VectorXd containing the steady-state abundances of the QSE species. + * + * This method solves a system of algebraic equations to determine the steady-state + * abundances of the QSE species. + * + * @throws std::runtime_error If the steady-state abundances cannot be calculated. + */ Eigen::VectorXd calculateSteadyStateAbundances( const std::vector& Y, const double T9, const double rho, const dynamicQSESpeciesIndices& indices ) const; + + /** + * @brief Initializes the network with a short ignition phase. + * @param netIn The input conditions for the network. + * @return The output conditions after the ignition phase. + * + * This method performs a short integration of the network at a high temperature and + * density to ignite the network and bring it closer to equilibrium. This can improve + * the convergence of the QSE solver. + * + * @see DirectNetworkSolver::evaluate() + */ NetOut initializeNetworkWithShortIgnition( const NetIn& netIn ) const; - private: // Nested functors for ODE integration - struct RHSFunctor { - DynamicEngine& m_engine; - const std::vector& m_dynamicSpeciesIndices; - const std::vector& m_QSESpeciesIndices; - const Eigen::VectorXd& m_Y_QSE; - const double m_T9; - const double m_rho; + /** + * @brief Determines whether the adaptive engine view should be updated. + * @param conditions The current input conditions. + * @return True if the view should be updated, false otherwise. + * + * This method implements a policy for determining when the adaptive engine view + * should be updated. The view is updated if the temperature or density has changed + * significantly, or if a primary fuel source has been depleted. + * + * @see AdaptiveEngineView::update() + */ + bool shouldUpdateView(const NetIn& conditions) const; + private: // Nested functors for ODE integration + /** + * @struct RHSFunctor + * @brief Functor for calculating the right-hand side of the ODEs for the dynamic species. + * + * This functor is used by the ODE solver to calculate the time derivatives of the + * dynamic species. It takes the current abundances of the dynamic species as input + * and returns the time derivatives of those abundances. + */ + struct RHSFunctor { + DynamicEngine& m_engine; ///< The engine used to evaluate the network. + const std::vector& m_dynamicSpeciesIndices; ///< Indices of the dynamic species. + const std::vector& m_QSESpeciesIndices; ///< Indices of the QSE species. + 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. + + bool m_isViewInitialized = false; + + /** + * @brief Constructor for the RHSFunctor. + * @param engine The engine used to evaluate the network. + * @param dynamicSpeciesIndices Indices of the dynamic species. + * @param QSESpeciesIndices Indices of the QSE species. + * @param Y_QSE Steady-state abundances of the QSE species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + */ RHSFunctor( DynamicEngine& engine, const std::vector& dynamicSpeciesIndices, @@ -77,6 +232,12 @@ namespace gridfire::solver { m_T9(T9), m_rho(rho) {} + /** + * @brief Calculates the time derivatives of the dynamic species. + * @param YDynamic Vector of current abundances for the dynamic species. + * @param dYdtDynamic Vector to store the time derivatives of the dynamic species. + * @param t Current time. + */ void operator()( const boost::numeric::ublas::vector& YDynamic, boost::numeric::ublas::vector& dYdtDynamic, @@ -85,13 +246,31 @@ namespace gridfire::solver { }; + /** + * @struct JacobianFunctor + * @brief Functor for calculating the Jacobian matrix of the ODEs for the dynamic species. + * + * This functor is used by the ODE solver to calculate the Jacobian matrix of the + * ODEs for the dynamic species. It takes the current abundances of the dynamic + * species as input and returns the Jacobian matrix. + * + * @todo Implement the Jacobian functor. + */ struct JacobianFunctor { - DynamicEngine& m_engine; - const std::vector& m_dynamicSpeciesIndices; - const std::vector& m_QSESpeciesIndices; - const double m_T9; - const double m_rho; + DynamicEngine& m_engine; ///< The engine used to evaluate the network. + 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. + /** + * @brief Constructor for the JacobianFunctor. + * @param engine The engine used to evaluate the network. + * @param dynamicSpeciesIndices Indices of the dynamic species. + * @param QSESpeciesIndices Indices of the QSE species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + */ JacobianFunctor( DynamicEngine& engine, const std::vector& dynamicSpeciesIndices, @@ -105,6 +284,13 @@ namespace gridfire::solver { m_T9(T9), m_rho(rho) {} + /** + * @brief Calculates the Jacobian matrix of the ODEs for the dynamic species. + * @param YDynamic Vector of current abundances for the dynamic species. + * @param JDynamic Matrix to store the Jacobian matrix. + * @param t Current time. + * @param dfdt Vector to store the time derivatives of the dynamic species (not used). + */ void operator()( const boost::numeric::ublas::vector& YDynamic, boost::numeric::ublas::matrix& JDynamic, @@ -113,19 +299,34 @@ namespace gridfire::solver { ) const; }; + /** + * @struct EigenFunctor + * @brief Functor for calculating the residual and Jacobian for the QSE species using Eigen. + * + * @tparam T The numeric type to use for the calculation (double or ADDouble). + */ template struct EigenFunctor { using InputType = Eigen::Matrix; using OutputType = Eigen::Matrix; using JacobianType = Eigen::Matrix; - DynamicEngine& m_engine; - const std::vector& m_YDynamic; - const std::vector& m_dynamicSpeciesIndices; - const std::vector& m_QSESpeciesIndices; - const double m_T9; - const double m_rho; + DynamicEngine& m_engine; ///< The engine used to evaluate the network. + const std::vector& m_YDynamic; ///< Abundances of the dynamic species. + 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. + /** + * @brief Constructor for the EigenFunctor. + * @param engine The engine used to evaluate the network. + * @param YDynamic Abundances of the dynamic species. + * @param dynamicSpeciesIndices Indices of the dynamic species. + * @param QSESpeciesIndices Indices of the QSE species. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + */ EigenFunctor( DynamicEngine& engine, const std::vector& YDynamic, @@ -141,25 +342,75 @@ namespace gridfire::solver { m_T9(T9), m_rho(rho) {} + /** + * @brief Calculates the residual vector for the QSE species. + * @param v_QSE Input vector of QSE species abundances (logarithmic). + * @param f_QSE Output vector of residuals. + * @return 0 for success. + */ int operator()(const InputType& v_QSE, OutputType& f_QSE) const; + + /** + * @brief Calculates the Jacobian matrix for the QSE species. + * @param v_QSE Input vector of QSE species abundances (logarithmic). + * @param J_QSE Output Jacobian matrix. + * @return 0 for success. + */ int df(const InputType& v_QSE, JacobianType& J_QSE) const; }; private: - quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - fourdst::config::Config& m_config = fourdst::config::Config::getInstance(); + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance. + fourdst::config::Config& m_config = fourdst::config::Config::getInstance(); ///< Configuration instance. + + bool m_isViewInitialized = false; ///< Flag indicating whether the adaptive engine view has been initialized. + NetIn m_lastSeenConditions; ///< The last seen input conditions. }; + /** + * @class DirectNetworkSolver + * @brief A network solver that directly integrates the reaction network ODEs. + * + * This solver uses a Runge-Kutta method to directly integrate the reaction network + * ODEs. It is simpler than the QSENetworkSolver, but it can be less efficient for + * stiff networks with disparate timescales. + * + * @implements DynamicNetworkSolverStrategy + */ class DirectNetworkSolver final : public DynamicNetworkSolverStrategy { public: + /** + * @brief Constructor for the DirectNetworkSolver. + * @param engine The dynamic engine to use for evaluating the network. + */ using DynamicNetworkSolverStrategy::DynamicNetworkSolverStrategy; + + /** + * @brief Evaluates the network for a given timestep using direct integration. + * @param netIn The input conditions for the network. + * @return The output conditions after the timestep. + */ NetOut evaluate(const NetIn& netIn) override; private: + /** + * @struct RHSFunctor + * @brief Functor for calculating the right-hand side of the ODEs. + * + * This functor is used by the ODE solver to calculate the time derivatives of the + * species abundances. It takes the current abundances as input and returns the + * time derivatives. + */ struct RHSFunctor { - DynamicEngine& m_engine; - const double m_T9; - const double m_rho; - const size_t m_numSpecies; + DynamicEngine& m_engine; ///< The engine used to evaluate the network. + const double m_T9; ///< Temperature in units of 10^9 K. + const double m_rho; ///< Density in g/cm^3. + const size_t m_numSpecies; ///< The number of species in the network. + /** + * @brief Constructor for the RHSFunctor. + * @param engine The engine used to evaluate the network. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + */ RHSFunctor( DynamicEngine& engine, const double T9, @@ -170,18 +421,38 @@ namespace gridfire::solver { m_rho(rho), m_numSpecies(engine.getNetworkSpecies().size()) {} + /** + * @brief Calculates the time derivatives of the species abundances. + * @param Y Vector of current abundances. + * @param dYdt Vector to store the time derivatives. + * @param t Current time. + */ void operator()( const boost::numeric::ublas::vector& Y, boost::numeric::ublas::vector& dYdt, double t ) const; }; - struct JacobianFunctor { - DynamicEngine& m_engine; - const double m_T9; - const double m_rho; - const size_t m_numSpecies; + /** + * @struct JacobianFunctor + * @brief Functor for calculating the Jacobian matrix. + * + * This functor is used by the ODE solver to calculate the Jacobian matrix of the + * ODEs. It takes the current abundances as input and returns the Jacobian matrix. + */ + struct JacobianFunctor { + DynamicEngine& m_engine; ///< The engine used to evaluate the network. + const double m_T9; ///< Temperature in units of 10^9 K. + const double m_rho; ///< Density in g/cm^3. + const size_t m_numSpecies; ///< The number of species in the network. + + /** + * @brief Constructor for the JacobianFunctor. + * @param engine The engine used to evaluate the network. + * @param T9 Temperature in units of 10^9 K. + * @param rho Density in g/cm^3. + */ JacobianFunctor( DynamicEngine& engine, const double T9, @@ -192,6 +463,13 @@ namespace gridfire::solver { m_rho(rho), m_numSpecies(engine.getNetworkSpecies().size()) {} + /** + * @brief Calculates the Jacobian matrix. + * @param Y Vector of current abundances. + * @param J Matrix to store the Jacobian matrix. + * @param t Current time. + * @param dfdt Vector to store the time derivatives (not used). + */ void operator()( const boost::numeric::ublas::vector& Y, boost::numeric::ublas::matrix& J, @@ -202,8 +480,8 @@ namespace gridfire::solver { }; private: - quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - fourdst::config::Config& m_config = fourdst::config::Config::getInstance(); + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance. + fourdst::config::Config& m_config = fourdst::config::Config::getInstance(); ///< Configuration instance. }; template diff --git a/src/network/lib/engine/engine_adaptive.cpp b/src/network/lib/engine/engine_adaptive.cpp index 333ef1ae..cdeb1045 100644 --- a/src/network/lib/engine/engine_adaptive.cpp +++ b/src/network/lib/engine/engine_adaptive.cpp @@ -1,10 +1,12 @@ -#include "gridfire/engine/engine_culled.h" +#include "gridfire/engine/engine_adaptive.h" #include +#include #include "gridfire/network.h" #include "quill/LogMacros.h" +#include "quill/Logger.h" namespace gridfire { using fourdst::atomic::Species; @@ -14,7 +16,8 @@ namespace gridfire { m_baseEngine(baseEngine), m_activeSpecies(baseEngine.getNetworkSpecies()), m_activeReactions(baseEngine.getNetworkReactions()), - m_speciesIndexMap(constructSpeciesIndexMap()) + m_speciesIndexMap(constructSpeciesIndexMap()), + m_reactionIndexMap(constructReactionIndexMap()) { } @@ -38,14 +41,47 @@ namespace gridfire { speciesIndexMap.push_back(it->second); } else { LOG_ERROR(m_logger, "Species '{}' not found in full species map.", active_species.name()); + m_logger -> flush_log(); throw std::runtime_error("Species not found in full species map: " + std::string(active_species.name())); } } - LOG_TRACE_L1(m_logger, "Successfully constructed species index map with {} entries.", speciesIndexMap.size()); + LOG_TRACE_L1(m_logger, "Species index map constructed with {} entries.", speciesIndexMap.size()); return speciesIndexMap; } + std::vector AdaptiveEngineView::constructReactionIndexMap() const { + LOG_TRACE_L1(m_logger, "Constructing reaction index map for adaptive engine view..."); + + // --- Step 1: Create a reverse map using the reaction's unique ID as the key. --- + std::unordered_map fullReactionReverseMap; + const auto& fullReactionSet = m_baseEngine.getNetworkReactions(); + fullReactionReverseMap.reserve(fullReactionSet.size()); + + for (size_t i_full = 0; i_full < fullReactionSet.size(); ++i_full) { + fullReactionReverseMap[fullReactionSet[i_full].id()] = i_full; + } + + // --- Step 2: Build the final index map using the active reaction set. --- + std::vector reactionIndexMap; + reactionIndexMap.reserve(m_activeReactions.size()); + + for (const auto& active_reaction_ptr : m_activeReactions) { + auto it = fullReactionReverseMap.find(active_reaction_ptr.id()); + + if (it != fullReactionReverseMap.end()) { + reactionIndexMap.push_back(it->second); + } else { + LOG_ERROR(m_logger, "Active reaction '{}' not found in base engine during reaction index map construction.", active_reaction_ptr.id()); + m_logger->flush_log(); + throw std::runtime_error("Mismatch between active reactions and base engine."); + } + } + + LOG_TRACE_L1(m_logger, "Reaction index map constructed with {} entries.", reactionIndexMap.size()); + return reactionIndexMap; + } + void AdaptiveEngineView::update(const NetIn& netIn) { LOG_TRACE_L1(m_logger, "Updating adaptive engine view..."); @@ -57,7 +93,7 @@ namespace gridfire { if (netIn.composition.contains(species)) { Y_full.push_back(netIn.composition.getMolarAbundance(std::string(species.name()))); } else { - LOG_DEBUG(m_logger, "Species '{}' not found in composition. Setting abundance to 0.0.", species.name()); + LOG_TRACE_L2(m_logger, "Species '{}' not found in composition. Setting abundance to 0.0.", species.name()); Y_full.push_back(0.0); } } @@ -65,25 +101,119 @@ namespace gridfire { const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) const double rho = netIn.density; // Density in g/cm^3 - m_isStale = false; std::vector reactionFlows; const auto& fullReactionSet = m_baseEngine.getNetworkReactions(); reactionFlows.reserve(fullReactionSet.size()); for (const auto& reactionPtr : fullReactionSet) { - const double flow = m_baseEngine.calculateMolarReactionFlow(*reactionPtr, Y_full, T9, rho); - reactionFlows.push_back({reactionPtr.get(), flow}); + const double flow = m_baseEngine.calculateMolarReactionFlow(reactionPtr, Y_full, T9, rho); + reactionFlows.push_back({&reactionPtr, flow}); } double max_flow = 0.0; + double min_flow = std::numeric_limits::max(); + double flowSum = 0.0; for (const auto&[reactionPtr, flowRate] : reactionFlows) { if (flowRate > max_flow) { max_flow = flowRate; + } else if (flowRate < min_flow) { + min_flow = flowRate; + } + flowSum += flowRate; + LOG_TRACE_L2(m_logger, "Reaction '{}' has flow rate: {:0.3E} [mol/s]", reactionPtr->id(), flowRate); + } + flowSum /= reactionFlows.size(); + + LOG_DEBUG(m_logger, "Maximum reaction flow rate in adaptive engine view: {:0.3E} [mol/s]", max_flow); + LOG_DEBUG(m_logger, "Minimum reaction flow rate in adaptive engine view: {:0.3E} [mol/s]", min_flow); + LOG_DEBUG(m_logger, "Average reaction flow rate in adaptive engine view: {:0.3E} [mol/s]", flowSum); + + const double relative_culling_threshold = m_config.get("gridfire:AdaptiveEngineView:RelativeCullingThreshold", 1e-75); + + double absoluteCullingThreshold = relative_culling_threshold * max_flow; + LOG_DEBUG(m_logger, "Relative culling threshold: {:0.3E} ({})", relative_culling_threshold, absoluteCullingThreshold); + + // --- Reaction Culling --- + LOG_TRACE_L1(m_logger, "Culling reactions based on reaction flow rates..."); + std::vector flowCulledReactions; + for (const auto&[reactionPtr, flowRate] : reactionFlows) { + if (flowRate > absoluteCullingThreshold) { + LOG_TRACE_L2(m_logger, "Maintaining reaction '{}' with relative (abs) flow rate: {:0.3E} ({:0.3E} [mol/s])", reactionPtr->id(), flowRate/max_flow, flowRate); + flowCulledReactions.push_back(reactionPtr); + } + } + LOG_DEBUG(m_logger, "Selected {} (total: {}, culled: {}) reactions based on flow rates.", flowCulledReactions.size(), fullReactionSet.size(), fullReactionSet.size() - flowCulledReactions.size()); + + // --- Connectivity Analysis --- + std::queue species_to_visit; + std::unordered_set reachable_species; + + constexpr double ABUNDANCE_FLOOR = 1e-12; // Abundance floor for a species to be considered part of the initial fuel + for (const auto& species : fullSpeciesList) { + if (netIn.composition.contains(species) && netIn.composition.getMassFraction(std::string(species.name())) > ABUNDANCE_FLOOR) { + species_to_visit.push(species); + reachable_species.insert(species); + LOG_TRACE_L2(m_logger, "Species '{}' is part of the initial fuel.", species.name()); + } + } + std::unordered_map> reactant_to_reactions_map; + for (const auto* reaction_ptr : flowCulledReactions) { + for (const auto& reactant : reaction_ptr->reactants()) { + reactant_to_reactions_map[reactant].push_back(reaction_ptr); } } - LOG_DEBUG(m_logger, "Maximum reaction flow rate in adaptive engine view: {:0.3E} [mol/s]", max_flow); + while (!species_to_visit.empty()) { + Species currentSpecies = species_to_visit.front(); + species_to_visit.pop(); + + auto it = reactant_to_reactions_map.find(currentSpecies); + if (it == reactant_to_reactions_map.end()) { + continue; // The species does not initiate any further reactions + } + + const auto& reactions = it->second; + for (const auto* reaction_ptr : reactions) { + for (const auto& product : reaction_ptr->products()) { + if (!reachable_species.contains(product)) { + reachable_species.insert(product); + species_to_visit.push(product); + LOG_TRACE_L2(m_logger, "Species '{}' is reachable via reaction '{}'.", product.name(), reaction_ptr->id()); + } + } + } + } + LOG_DEBUG(m_logger, "Reachable species count: {}", reachable_species.size()); + + m_activeSpecies.assign(reachable_species.begin(), reachable_species.end()); + std::ranges::sort(m_activeSpecies, + [](const Species &a, const Species &b) { return a.mass() < b.mass(); }); + + m_activeReactions.clear(); + for (const auto* reaction_ptr : flowCulledReactions) { + bool all_reactants_present = true; + for (const auto& reactant : reaction_ptr->reactants()) { + if (!reachable_species.contains(reactant)) { + all_reactants_present = false; + break; + } + } + + if (all_reactants_present) { + m_activeReactions.add_reaction(*reaction_ptr); + LOG_TRACE_L2(m_logger, "Maintaining reaction '{}' with all reactants present.", reaction_ptr->id()); + } else { + LOG_TRACE_L1(m_logger, "Culling reaction '{}' due to missing reactants.", reaction_ptr->id()); + } + } + LOG_DEBUG(m_logger, "Active reactions count: {} (total: {}, culled: {}, culled due to connectivity: {})", m_activeReactions.size(), + fullReactionSet.size(), fullReactionSet.size() - m_activeReactions.size(), flowCulledReactions.size() - m_activeReactions.size()); + + m_speciesIndexMap = constructSpeciesIndexMap(); + m_reactionIndexMap = constructReactionIndexMap(); + + m_isStale = false; } const std::vector & AdaptiveEngineView::getNetworkSpecies() const { @@ -91,65 +221,143 @@ namespace gridfire { } StepDerivatives AdaptiveEngineView::calculateRHSAndEnergy( - const std::vector &Y, + const std::vector &Y_culled, const double T9, const double rho ) const { - return m_baseEngine.calculateRHSAndEnergy(Y, T9, rho); + validateState(); + + const auto Y_full = mapCulledToFull(Y_culled); + + const auto [dydt, nuclearEnergyGenerationRate] = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); + + StepDerivatives culledResults; + culledResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate; + culledResults.dydt = mapFullToCulled(dydt); + + return culledResults; } void AdaptiveEngineView::generateJacobianMatrix( - const std::vector &Y, + const std::vector &Y_culled, const double T9, const double rho ) { - m_baseEngine.generateJacobianMatrix(Y, T9, rho); + validateState(); + const auto Y_full = mapCulledToFull(Y_culled); + + m_baseEngine.generateJacobianMatrix(Y_full, T9, rho); } double AdaptiveEngineView::getJacobianMatrixEntry( - const int i, - const int j + const int i_culled, + const int j_culled ) const { - return m_baseEngine.getJacobianMatrixEntry(i, j); + validateState(); + const size_t i_full = mapCulledToFullSpeciesIndex(i_culled); + const size_t j_full = mapCulledToFullSpeciesIndex(j_culled); + + return m_baseEngine.getJacobianMatrixEntry(i_full, j_full); } void AdaptiveEngineView::generateStoichiometryMatrix() { + validateState(); m_baseEngine.generateStoichiometryMatrix(); } int AdaptiveEngineView::getStoichiometryMatrixEntry( - const int speciesIndex, - const int reactionIndex + const int speciesIndex_culled, + const int reactionIndex_culled ) const { - return m_baseEngine.getStoichiometryMatrixEntry(speciesIndex, reactionIndex); + validateState(); + const size_t speciesIndex_full = mapCulledToFullSpeciesIndex(speciesIndex_culled); + const size_t reactionIndex_full = mapCulledToFullReactionIndex(reactionIndex_culled); + return m_baseEngine.getStoichiometryMatrixEntry(speciesIndex_full, reactionIndex_full); } double AdaptiveEngineView::calculateMolarReactionFlow( const reaction::Reaction &reaction, - const std::vector &Y, + const std::vector &Y_culled, const double T9, const double rho ) const { + validateState(); + if (!m_activeReactions.contains(reaction)) { + LOG_ERROR(m_logger, "Reaction '{}' is not part of the active reactions in the adaptive engine view.", reaction.id()); + m_logger -> flush_log(); + throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id())); + } + const auto Y = mapCulledToFull(Y_culled); + return m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho); } - const reaction::REACLIBLogicalReactionSet & AdaptiveEngineView::getNetworkReactions() const { + const reaction::LogicalReactionSet & AdaptiveEngineView::getNetworkReactions() const { return m_activeReactions; } - std::unordered_map AdaptiveEngineView::getSpeciesTimescales( - const std::vector &Y, + std::unordered_map AdaptiveEngineView::getSpeciesTimescales( + const std::vector &Y_culled, const double T9, const double rho ) const { - auto timescales = m_baseEngine.getSpeciesTimescales(Y, T9, rho); - for (const auto &species: timescales | std::views::keys) { - // remove species that are not in the active species list - if (std::ranges::find(m_activeSpecies, species) == m_activeSpecies.end()) { - timescales.erase(species); + validateState(); + const auto Y_full = mapCulledToFull(Y_culled); + const auto fullTimescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + + std::unordered_map culledTimescales; + culledTimescales.reserve(m_activeSpecies.size()); + for (const auto& active_species : m_activeSpecies) { + if (fullTimescales.contains(active_species)) { + culledTimescales[active_species] = fullTimescales.at(active_species); } } - return timescales; + return culledTimescales; + + } + + 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) { + const size_t i_full = m_speciesIndexMap[i_culled]; + full[i_full] += culled[i_culled]; + } + return full; + } + + std::vector AdaptiveEngineView::mapFullToCulled(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]; + culled[i_culled] = full[i_full]; + } + return culled; + } + + size_t AdaptiveEngineView::mapCulledToFullSpeciesIndex(size_t culledSpeciesIndex) const { + if (culledSpeciesIndex < 0 || culledSpeciesIndex >= static_cast(m_speciesIndexMap.size())) { + LOG_ERROR(m_logger, "Culled index {} is out of bounds for species index map of size {}.", culledSpeciesIndex, m_speciesIndexMap.size()); + m_logger->flush_log(); + throw std::out_of_range("Culled index " + std::to_string(culledSpeciesIndex) + " is out of bounds for species index map of size " + std::to_string(m_speciesIndexMap.size()) + "."); + } + return m_speciesIndexMap[culledSpeciesIndex]; + } + + size_t AdaptiveEngineView::mapCulledToFullReactionIndex(size_t culledReactionIndex) const { + if (culledReactionIndex < 0 || culledReactionIndex >= static_cast(m_reactionIndexMap.size())) { + LOG_ERROR(m_logger, "Culled index {} is out of bounds for reaction index map of size {}.", culledReactionIndex, m_reactionIndexMap.size()); + m_logger->flush_log(); + throw std::out_of_range("Culled index " + std::to_string(culledReactionIndex) + " is out of bounds for reaction index map of size " + std::to_string(m_reactionIndexMap.size()) + "."); + } + return m_reactionIndexMap[culledReactionIndex]; + } + + void AdaptiveEngineView::validateState() const { + if (m_isStale) { + LOG_ERROR(m_logger, "AdaptiveEngineView is stale. Please call update() before calculating RHS and energy."); + m_logger->flush_log(); + throw std::runtime_error("AdaptiveEngineView is stale. Please call update() before calculating RHS and energy."); + } } } diff --git a/src/network/lib/engine/engine_graph.cpp b/src/network/lib/engine/engine_graph.cpp index 5d0bfafd..92c01006 100644 --- a/src/network/lib/engine/engine_graph.cpp +++ b/src/network/lib/engine/engine_graph.cpp @@ -18,7 +18,6 @@ #include #include -#include #include @@ -30,7 +29,7 @@ namespace gridfire { syncInternalMaps(); } - GraphEngine::GraphEngine(reaction::REACLIBLogicalReactionSet reactions) : + GraphEngine::GraphEngine(reaction::LogicalReactionSet reactions) : m_reactions(std::move(reactions)) { syncInternalMaps(); } @@ -61,21 +60,22 @@ namespace gridfire { std::set uniqueSpeciesNames; for (const auto& reaction: m_reactions) { - for (const auto& reactant: reaction->reactants()) { + for (const auto& reactant: reaction.reactants()) { uniqueSpeciesNames.insert(reactant.name()); } - for (const auto& product: reaction->products()) { + for (const auto& product: reaction.products()) { uniqueSpeciesNames.insert(product.name()); } } for (const auto& name: uniqueSpeciesNames) { - auto it = fourdst::atomic::species.find(name); + auto it = fourdst::atomic::species.find(std::string(name)); if (it != fourdst::atomic::species.end()) { m_networkSpecies.push_back(it->second); m_networkSpeciesMap.insert({name, it->second}); } else { LOG_ERROR(m_logger, "Species '{}' not found in global atomic species database.", name); + m_logger->flush_log(); throw std::runtime_error("Species not found in global atomic species database: " + std::string(name)); } } @@ -85,8 +85,8 @@ namespace gridfire { void GraphEngine::populateReactionIDMap() { LOG_TRACE_L1(m_logger, "Populating reaction ID map for REACLIB graph network (serif::network::GraphNetwork)..."); m_reactionIDMap.clear(); - for (const auto& reaction: m_reactions) { - m_reactionIDMap.emplace(reaction->id(), reaction.get()); + for (auto& reaction: m_reactions) { + m_reactionIDMap.emplace(reaction.id(), &reaction); } LOG_TRACE_L1(m_logger, "Populated {} reactions in the reaction ID map.", m_reactionIDMap.size()); } @@ -111,13 +111,13 @@ namespace gridfire { // --- Basic Accessors and Queries --- const std::vector& GraphEngine::getNetworkSpecies() const { // Returns a constant reference to the vector of unique species in the network. - LOG_DEBUG(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size()); + LOG_TRACE_L3(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size()); return m_networkSpecies; } - const reaction::REACLIBLogicalReactionSet& GraphEngine::getNetworkReactions() const { + const reaction::LogicalReactionSet& GraphEngine::getNetworkReactions() const { // Returns a constant reference to the set of reactions in the network. - LOG_DEBUG(m_logger, "Providing access to network reactions set. Size: {}.", m_reactions.size()); + LOG_TRACE_L3(m_logger, "Providing access to network reactions set. Size: {}.", m_reactions.size()); return m_reactions; } @@ -139,7 +139,7 @@ namespace gridfire { uint64_t totalProductZ = 0; // Calculate total A and Z for reactants - for (const auto& reactant : reaction->reactants()) { + for (const auto& reactant : reaction.reactants()) { auto it = m_networkSpeciesMap.find(reactant.name()); if (it != m_networkSpeciesMap.end()) { totalReactantA += it->second.a(); @@ -148,13 +148,13 @@ namespace gridfire { // This scenario indicates a severe data integrity issue: // a reactant is part of a reaction but not in the network's species map. LOG_ERROR(m_logger, "CRITICAL ERROR: Reactant species '{}' in reaction '{}' not found in network species map during conservation validation.", - reactant.name(), reaction->id()); + reactant.name(), reaction.id()); return false; } } // Calculate total A and Z for products - for (const auto& product : reaction->products()) { + for (const auto& product : reaction.products()) { auto it = m_networkSpeciesMap.find(product.name()); if (it != m_networkSpeciesMap.end()) { totalProductA += it->second.a(); @@ -162,7 +162,7 @@ namespace gridfire { } else { // Similar critical error for product species LOG_ERROR(m_logger, "CRITICAL ERROR: Product species '{}' in reaction '{}' not found in network species map during conservation validation.", - product.name(), reaction->id()); + product.name(), reaction.id()); return false; } } @@ -170,12 +170,12 @@ namespace gridfire { // Compare totals for conservation if (totalReactantA != totalProductA) { LOG_ERROR(m_logger, "Mass number (A) not conserved for reaction '{}': Reactants A={} vs Products A={}.", - reaction->id(), totalReactantA, totalProductA); + reaction.id(), totalReactantA, totalProductA); return false; } if (totalReactantZ != totalProductZ) { LOG_ERROR(m_logger, "Atomic number (Z) not conserved for reaction '{}': Reactants Z={} vs Products Z={}.", - reaction->id(), totalReactantZ, totalProductZ); + reaction.id(), totalReactantZ, totalProductZ); return false; } } @@ -187,7 +187,7 @@ namespace gridfire { void GraphEngine::validateComposition(const fourdst::composition::Composition &composition, double culling, double T9) { // Check if the requested network has already been cached. // PERF: Rebuilding this should be pretty fast but it may be a good point of optimization in the future. - const reaction::REACLIBLogicalReactionSet validationReactionSet = build_reaclib_nuclear_network(composition, false); + const reaction::LogicalReactionSet validationReactionSet = build_reaclib_nuclear_network(composition, false); // TODO: need some more robust method here to // A. Build the basic network from the composition's species with non zero mass fractions. // B. rebuild a new composition from the reaction set's reactants + products (with the mass fractions from the things that are only products set to 0) @@ -199,7 +199,7 @@ namespace gridfire { // This allows for dynamic network modification while retaining caching for networks which are very similar. if (validationReactionSet != m_reactions) { LOG_DEBUG(m_logger, "Reaction set not cached. Rebuilding the reaction set for T9={} and culling={}.", T9, culling); - m_reactions = validationReactionSet; + m_reactions = std::move(validationReactionSet); syncInternalMaps(); // Re-sync internal maps after updating reactions. Note this will also retrace the AD tape. } } @@ -221,7 +221,7 @@ namespace gridfire { size_t reactionColumnIndex = 0; for (const auto& reaction : m_reactions) { // Get the net stoichiometry for the current reaction - std::unordered_map netStoichiometry = reaction->stoichiometry(); + std::unordered_map netStoichiometry = reaction.stoichiometry(); // Iterate through the species and their coefficients in the stoichiometry map for (const auto& [species, coefficient] : netStoichiometry) { @@ -234,7 +234,8 @@ namespace gridfire { } else { // This scenario should ideally not happen if m_networkSpeciesMap and m_speciesToIndexMap are correctly synced LOG_ERROR(m_logger, "CRITICAL ERROR: Species '{}' from reaction '{}' stoichiometry not found in species to index map.", - species.name(), reaction->id()); + species.name(), reaction.id()); + m_logger -> flush_log(); throw std::runtime_error("Species not found in species to index map: " + std::string(species.name())); } } @@ -255,8 +256,8 @@ namespace gridfire { StepDerivatives GraphEngine::calculateAllDerivatives( const std::vector &Y_in, - const ADDouble T9, - const ADDouble rho + const ADDouble &T9, + const ADDouble &rho ) const { return calculateAllDerivatives(Y_in, T9, rho); } @@ -300,7 +301,7 @@ namespace gridfire { } } } - LOG_DEBUG(m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); + 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 { @@ -309,7 +310,7 @@ namespace gridfire { std::unordered_map GraphEngine::getNetReactionStoichiometry( const reaction::Reaction &reaction - ) const { + ) { return reaction.stoichiometry(); } @@ -326,6 +327,7 @@ namespace gridfire { std::ofstream dotFile(filename); if (!dotFile.is_open()) { LOG_ERROR(m_logger, "Failed to open file for writing: {}", filename); + m_logger->flush_log(); throw std::runtime_error("Failed to open file for writing: " + filename); } @@ -345,19 +347,19 @@ namespace gridfire { dotFile << " // --- Reaction Edges ---\n"; for (const auto& reaction : m_reactions) { // Create a unique ID for the reaction node - std::string reactionNodeId = "reaction_" + std::string(reaction->id()); + std::string reactionNodeId = "reaction_" + std::string(reaction.id()); // Define the reaction node (small, black dot) dotFile << " \"" << reactionNodeId << "\" [shape=point, fillcolor=black, width=0.1, height=0.1, label=\"\"];\n"; // Draw edges from reactants to the reaction node - for (const auto& reactant : reaction->reactants()) { + for (const auto& reactant : reaction.reactants()) { dotFile << " \"" << reactant.name() << "\" -> \"" << reactionNodeId << "\";\n"; } // Draw edges from the reaction node to products - for (const auto& product : reaction->products()) { - dotFile << " \"" << reactionNodeId << "\" -> \"" << product.name() << "\" [label=\"" << reaction->qValue() << " MeV\"];\n"; + for (const auto& product : reaction.products()) { + dotFile << " \"" << reactionNodeId << "\" -> \"" << product.name() << "\" [label=\"" << reaction.qValue() << " MeV\"];\n"; } dotFile << "\n"; } @@ -373,36 +375,32 @@ namespace gridfire { std::ofstream csvFile(filename, std::ios::out | std::ios::trunc); if (!csvFile.is_open()) { LOG_ERROR(m_logger, "Failed to open file for writing: {}", filename); + m_logger->flush_log(); throw std::runtime_error("Failed to open file for writing: " + filename); } csvFile << "Reaction;Reactants;Products;Q-value;sources;rates\n"; for (const auto& reaction : m_reactions) { // Dynamic cast to REACLIBReaction to access specific properties - csvFile << reaction->id() << ";"; + csvFile << reaction.id() << ";"; // Reactants int count = 0; - for (const auto& reactant : reaction->reactants()) { + for (const auto& reactant : reaction.reactants()) { csvFile << reactant.name(); - if (++count < reaction->reactants().size()) { + if (++count < reaction.reactants().size()) { csvFile << ","; } } csvFile << ";"; count = 0; - for (const auto& product : reaction->products()) { + for (const auto& product : reaction.products()) { csvFile << product.name(); - if (++count < reaction->products().size()) { + if (++count < reaction.products().size()) { csvFile << ","; } } - csvFile << ";" << reaction->qValue() << ";"; + csvFile << ";" << reaction.qValue() << ";"; // Reaction coefficients - auto* reaclibReaction = dynamic_cast(reaction.get()); - if (!reaclibReaction) { - LOG_ERROR(m_logger, "Failed to cast Reaction to REACLIBLogicalReaction in GraphNetwork::exportToCSV()."); - throw std::runtime_error("Failed to cast Reaction to REACLIBLogicalReaction in GraphNetwork::exportToCSV(). This should not happen, please check your reaction setup."); - } - auto sources = reaclibReaction->sources(); + auto sources = reaction.sources(); count = 0; for (const auto& source : sources) { csvFile << source; @@ -413,9 +411,9 @@ namespace gridfire { csvFile << ";"; // Reaction coefficients count = 0; - for (const auto& rates : *reaclibReaction) { + for (const auto& rates : reaction) { csvFile << rates; - if (++count < reaclibReaction->size()) { + if (++count < reaction.size()) { csvFile << ","; } } @@ -448,6 +446,7 @@ namespace gridfire { const size_t numSpecies = m_networkSpecies.size(); if (numSpecies == 0) { LOG_ERROR(m_logger, "Cannot record AD tape: No species in the network."); + m_logger->flush_log(); throw std::runtime_error("Cannot record AD tape: No species in the network."); } const size_t numADInputs = numSpecies + 2; // Note here that by not letting T9 and rho be independent variables, we are constraining the network to a constant temperature and density during each evaluation. diff --git a/src/network/lib/network.cpp b/src/network/lib/network.cpp index 6793cf90..c6dbaabd 100644 --- a/src/network/lib/network.cpp +++ b/src/network/lib/network.cpp @@ -19,15 +19,25 @@ // // *********************************************************************** */ #include "gridfire/network.h" -#include "gridfire/reactions.h" -#include "../include/gridfire/reaction/reaction.h" +#include "gridfire/reaction/reaclib.h" +#include "gridfire/reaction/reaction.h" #include -#include #include "quill/LogMacros.h" namespace gridfire { + std::vector NetIn::MolarAbundance() const { + std::vector y; + y.reserve(composition.getRegisteredSymbols().size()); + const auto [fst, snd] = composition.getComposition(); + for (const auto &name: fst | std::views::keys) { + y.push_back(composition.getMolarAbundance(name)); + } + return y; + } + + Network::Network(const NetworkFormat format) : m_config(fourdst::config::Config::getInstance()), m_logManager(fourdst::logging::LogManager::getInstance()), @@ -36,6 +46,7 @@ namespace gridfire { m_constants(fourdst::constant::Constants::getInstance()){ if (format == NetworkFormat::UNKNOWN) { LOG_ERROR(m_logger, "nuclearNetwork::Network::Network() called with UNKNOWN format"); + m_logger->flush_log(); throw std::runtime_error("nuclearNetwork::Network::Network() called with UNKNOWN format"); } } @@ -50,17 +61,12 @@ namespace gridfire { return oldFormat; } - reaction::REACLIBLogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, bool reverse) { + reaction::LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, bool reverse) { using namespace reaction; - std::vector reaclibReactions; + std::vector reaclibReactions; auto logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - if (!reaclib::s_initialized) { - LOG_DEBUG(logger, "REACLIB reactions not initialized. Calling initializeAllReaclibReactions()..."); - reaclib::initializeAllReaclibReactions(); - } - - for (const auto &reaction: reaclib::s_all_reaclib_reactions | std::views::values) { + for (const auto &reaction: reaclib::get_all_reactions()) { if (reaction.is_reverse() != reverse) { continue; // Skip reactions that do not match the requested direction } @@ -77,8 +83,8 @@ namespace gridfire { reaclibReactions.push_back(reaction); } } - const REACLIBReactionSet reactionSet(reaclibReactions); - return REACLIBLogicalReactionSet(reactionSet); + const ReactionSet reactionSet(reaclibReactions); + return LogicalReactionSet(reactionSet); } // Trim whitespace from both ends of a string @@ -97,61 +103,4 @@ namespace gridfire { return std::string(startIt, ritr.base()); } - reaction::REACLIBLogicalReactionSet build_reaclib_nuclear_network_from_file(const std::string& filename, bool reverse) { - const auto logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - std::vector reactionPENames; - std::ifstream infile(filename, std::ios::in); - if (!infile.is_open()) { - LOG_ERROR(logger, "Failed to open network definition file: {}", filename); - throw std::runtime_error("Failed to open network definition file: " + filename); - } - std::string line; - while (std::getline(infile, line)) { - std::string trimmedLine = trim_whitespace(line); - if (trimmedLine.empty() || trimmedLine.starts_with('#')) { - continue; // Skip empty lines and comments - } - reactionPENames.push_back(trimmedLine); - } - infile.close(); - std::vector reaclibReactions; - - if (reactionPENames.empty()) { - LOG_ERROR(logger, "No reactions found in the network definition file: {}", filename); - throw std::runtime_error("No reactions found in the network definition file: " + filename); - } - - if (!reaclib::s_initialized) { - LOG_DEBUG(logger, "REACLIB reactions not initialized. Calling initializeAllReaclibReactions()..."); - reaclib::initializeAllReaclibReactions(); - } else { - LOG_DEBUG(logger, "REACLIB reactions already initialized."); - } - - for (const auto& peName : reactionPENames) { - bool found = false; - for (const auto& reaction : reaclib::s_all_reaclib_reactions | std::views::values) { - if (reaction.peName() == peName && reaction.is_reverse() == reverse) { - reaclibReactions.push_back(reaction); - found = true; - LOG_TRACE_L3(logger, "Found reaction {} (version {}) in REACLIB database.", peName, reaction.sourceLabel()); - } - } - if (!found) { - LOG_ERROR(logger, "Reaction {} not found in REACLIB database. Skipping...", peName); - throw std::runtime_error("Reaction not found in REACLIB database."); - } - } - - - // Log any reactions that were not found in the REACLIB database - for (const auto& peName : reactionPENames) { - if (std::ranges::find(reaclibReactions, peName, &reaction::REACLIBReaction::peName) == reaclibReactions.end()) { - LOG_WARNING(logger, "Reaction {} not found in REACLIB database.", peName); - } - } - const reaction::REACLIBReactionSet reactionSet(reaclibReactions); - return reaction::REACLIBLogicalReactionSet(reactionSet); - } - } diff --git a/src/network/lib/reaction/reaclib.cpp b/src/network/lib/reaction/reaclib.cpp index b7694f86..f29e11b8 100644 --- a/src/network/lib/reaction/reaclib.cpp +++ b/src/network/lib/reaction/reaclib.cpp @@ -1,3 +1,146 @@ -// -// Created by Emily Boudreaux on 6/28/25. -// +#include "fourdst/composition/atomicSpecies.h" +#include "fourdst/composition/species.h" + +#include "gridfire/reaction/reaclib.h" +#include "gridfire/reaction/reactions_data.h" +#include "gridfire/network.h" + +#include +#include +#include +#include + +std::string trim_whitespace(const std::string& str) { + auto startIt = str.begin(); + auto endIt = str.end(); + + while (startIt != endIt && std::isspace(static_cast(*startIt))) { + ++startIt; + } + if (startIt == endIt) { + return ""; + } + auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt), + [](unsigned char ch){ return !std::isspace(ch); }); + return std::string(startIt, ritr.base()); +} + +namespace gridfire::reaclib { + static reaction::LogicalReactionSet* s_all_reaclib_reactions_ptr = nullptr; + + #pragma pack(push, 1) + struct ReactionRecord { + int32_t chapter; + double qValue; + double coeffs[7]; + bool reverse; + char label[8]; + char rpName[64]; + char reactants_str[128]; + char products_str[128]; + }; + #pragma pack(pop) + + std::ostream& operator<<(std::ostream& os, const ReactionRecord& r) { + os << "Chapter: " << r.chapter + << ", Q-value: " << r.qValue + << ", Coefficients: [" << r.coeffs[0] << ", " << r.coeffs[1] << ", " + << r.coeffs[2] << ", " << r.coeffs[3] << ", " << r.coeffs[4] << ", " + << r.coeffs[5] << ", " << r.coeffs[6] << "]" + << ", Reverse: " << (r.reverse ? "true" : "false") + << ", Label: '" << std::string(r.label, strnlen(r.label, sizeof(r.label))) << "'" + << ", RP Name: '" << std::string(r.rpName, strnlen(r.rpName, sizeof(r.rpName))) << "'" + << ", Reactants: '" << std::string(r.reactants_str, strnlen(r.reactants_str, sizeof(r.reactants_str))) << "'" + << ", Products: '" << std::string(r.products_str, strnlen(r.products_str, sizeof(r.products_str))) << "'"; + return os; + } + + static std::vector parseSpeciesString(const std::string_view str) { + std::vector result; + std::stringstream ss{std::string(str)}; + std::string name; + + while (ss >> name) { + // Trim whitespace that might be left over from the fixed-width char arrays + const auto trimmed_name = trim_whitespace(name); + if (trimmed_name.empty()) continue; + + auto it = fourdst::atomic::species.find(trimmed_name); + if (it != fourdst::atomic::species.end()) { + result.push_back(it->second); + } else { + // If a species is not found, it's a critical data error. + throw std::runtime_error("Unknown species in reaction data: " + std::string(trimmed_name)); + } + } + return result; + } + + static void initializeAllReaclibReactions() { + if (s_initialized) { + return; + } + + // Cast the raw byte data to our structured record format. + const auto* records = reinterpret_cast(raw_reactions_data); + const size_t num_reactions = raw_reactions_data_len / sizeof(ReactionRecord); + + std::vector reaction_list; + reaction_list.reserve(num_reactions); + + for (size_t i = 0; i < num_reactions; ++i) { + const auto& record = records[i]; + + // The char arrays from the binary are not guaranteed to be null-terminated + // if the string fills the entire buffer. We create null-terminated string_views. + const std::string_view label_sv(record.label, strnlen(record.label, sizeof(record.label))); + const std::string_view rpName_sv(record.rpName, strnlen(record.rpName, sizeof(record.rpName))); + const std::string_view reactants_sv(record.reactants_str, strnlen(record.reactants_str, sizeof(record.reactants_str))); + const std::string_view products_sv(record.products_str, strnlen(record.products_str, sizeof(record.products_str))); + + auto reactants = parseSpeciesString(reactants_sv); + auto products = parseSpeciesString(products_sv); + + const reaction::RateCoefficientSet rate_coeffs = { + record.coeffs[0], record.coeffs[1], record.coeffs[2], + record.coeffs[3], record.coeffs[4], record.coeffs[5], + record.coeffs[6] + }; + + // Construct the Reaction object. We use rpName for both the unique ID and the human-readable name. + reaction_list.emplace_back( + rpName_sv, + rpName_sv, + record.chapter, + reactants, + products, + record.qValue, + label_sv, + rate_coeffs, + record.reverse + ); + } + + // The ReactionSet takes the vector of all individual reactions. + reaction::ReactionSet reaction_set(std::move(reaction_list)); + + // The LogicalReactionSet groups reactions by their peName, which is what we want. + s_all_reaclib_reactions_ptr = new reaction::LogicalReactionSet(reaction_set); + + s_initialized = true; + } + + + // --- Public Interface Implementation --- + + const reaction::LogicalReactionSet& get_all_reactions() { + // This ensures that the initialization happens only on the first call. + if (!s_initialized) { + initializeAllReaclibReactions(); + } + if (s_all_reaclib_reactions_ptr == nullptr) { + throw std::runtime_error("Reaclib reactions have not been initialized."); + } + return *s_all_reaclib_reactions_ptr; + } +} // namespace gridfire::reaclib \ No newline at end of file diff --git a/src/network/lib/reaction/reaction.cpp b/src/network/lib/reaction/reaction.cpp index ce68e9bf..b72c4d62 100644 --- a/src/network/lib/reaction/reaction.cpp +++ b/src/network/lib/reaction/reaction.cpp @@ -3,9 +3,9 @@ #include #include #include -#include #include #include +#include #include "quill/LogMacros.h" @@ -18,15 +18,31 @@ namespace gridfire::reaction { Reaction::Reaction( const std::string_view id, - const double qValue, + const std::string_view peName, + const int chapter, const std::vector& reactants, const std::vector& products, + const double qValue, + const std::string_view label, + const RateCoefficientSet& sets, const bool reverse) : - m_id(id), - m_qValue(qValue), - m_reactants(std::move(reactants)), - m_products(std::move(products)), - m_reverse(reverse) {} + m_id(id), + m_peName(peName), + m_chapter(chapter), + m_qValue(qValue), + m_reactants(reactants), + m_products(products), + m_sourceLabel(label), + m_rateCoefficients(sets), + m_reverse(reverse) {} + + double Reaction::calculate_rate(const double T9) const { + return calculate_rate(T9); + } + + CppAD::AD Reaction::calculate_rate(const CppAD::AD T9) const { + return calculate_rate>(T9); + } bool Reaction::contains(const Species &species) const { return contains_reactant(species) || contains_product(species); @@ -122,27 +138,28 @@ namespace gridfire::reaction { } ReactionSet::ReactionSet( - std::vector> reactions) : - m_reactions(std::move(reactions)) { + std::vector reactions + ) : + m_reactions(std::move(reactions)) { if (m_reactions.empty()) { return; // Case where the reactions will be added later. } m_reactionNameMap.reserve(reactions.size()); for (const auto& reaction : m_reactions) { - m_id += reaction->id(); - m_reactionNameMap.emplace(reaction->id(), reaction.get()); + m_id += reaction.id(); + m_reactionNameMap.emplace(reaction.id(), reaction); } } ReactionSet::ReactionSet(const ReactionSet &other) { m_reactions.reserve(other.m_reactions.size()); for (const auto& reaction_ptr: other.m_reactions) { - m_reactions.push_back(reaction_ptr->clone()); + m_reactions.push_back(reaction_ptr); } m_reactionNameMap.reserve(other.m_reactionNameMap.size()); for (const auto& reaction_ptr : m_reactions) { - m_reactionNameMap.emplace(reaction_ptr->id(), reaction_ptr.get()); + m_reactionNameMap.emplace(reaction_ptr.id(), reaction_ptr); } } @@ -155,28 +172,27 @@ namespace gridfire::reaction { return *this; } - void ReactionSet::add_reaction(std::unique_ptr reaction) { - m_reactions.emplace_back(std::move(reaction)); - m_id += m_reactions.back()->id(); - m_reactionNameMap.emplace(m_reactions.back()->id(), m_reactions.back().get()); + void ReactionSet::add_reaction(Reaction reaction) { + m_reactions.emplace_back(reaction); + m_id += m_reactions.back().id(); + m_reactionNameMap.emplace(m_reactions.back().id(), m_reactions.back()); } - void ReactionSet::remove_reaction(const std::unique_ptr& reaction) { - if (!m_reactionNameMap.contains(std::string(reaction->id()))) { - // LOG_INFO(m_logger, "Attempted to remove reaction {} that does not exist in ReactionSet. Skipping...", reaction->id()); + void ReactionSet::remove_reaction(const Reaction& reaction) { + if (!m_reactionNameMap.contains(std::string(reaction.id()))) { return; } - m_reactionNameMap.erase(std::string(reaction->id())); + m_reactionNameMap.erase(std::string(reaction.id())); - std::erase_if(m_reactions, [&reaction](const std::unique_ptr& r) { - return *r == *reaction; + std::erase_if(m_reactions, [&reaction](const Reaction& r) { + return r == reaction; }); } bool ReactionSet::contains(const std::string_view& id) const { for (const auto& reaction : m_reactions) { - if (reaction->id() == id) { + if (reaction.id() == id) { return true; } } @@ -185,23 +201,21 @@ namespace gridfire::reaction { bool ReactionSet::contains(const Reaction& reaction) const { for (const auto& r : m_reactions) { - if (*r == reaction) { + if (r == reaction) { return true; } } return false; } - void ReactionSet::sort(double T9) { - std::ranges::sort(m_reactions, - [&T9](const std::unique_ptr& r1, const std::unique_ptr& r2) { - return r1->calculate_rate(T9) < r2->calculate_rate(T9); - }); + void ReactionSet::clear() { + m_reactions.clear(); + m_reactionNameMap.clear(); } bool ReactionSet::contains_species(const Species& species) const { for (const auto& reaction : m_reactions) { - if (reaction->contains(species)) { + if (reaction.contains(species)) { return true; } } @@ -210,7 +224,7 @@ namespace gridfire::reaction { bool ReactionSet::contains_reactant(const Species& species) const { for (const auto& r : m_reactions) { - if (r->contains_reactant(species)) { + if (r.contains_reactant(species)) { return true; } } @@ -219,7 +233,7 @@ namespace gridfire::reaction { bool ReactionSet::contains_product(const Species& species) const { for (const auto& r : m_reactions) { - if (r->contains_product(species)) { + if (r.contains_product(species)) { return true; } } @@ -228,15 +242,17 @@ namespace gridfire::reaction { const Reaction& ReactionSet::operator[](const size_t index) const { if (index >= m_reactions.size()) { + m_logger -> flush_log(); throw std::out_of_range("Index" + std::to_string(index) + " out of range for ReactionSet of size " + std::to_string(m_reactions.size()) + "."); } - return *m_reactions[index]; + return m_reactions[index]; } const Reaction& ReactionSet::operator[](const std::string_view& id) const { if (auto it = m_reactionNameMap.find(std::string(id)); it != m_reactionNameMap.end()) { - return *it->second; + return it->second; } + m_logger -> flush_log(); throw std::out_of_range("Species " + std::string(id) + " does not exist in ReactionSet."); } @@ -258,7 +274,7 @@ namespace gridfire::reaction { std::vector individualReactionHashes; individualReactionHashes.reserve(m_reactions.size()); for (const auto& reaction : m_reactions) { - individualReactionHashes.push_back(reaction->hash(seed)); + individualReactionHashes.push_back(reaction.hash(seed)); } std::ranges::sort(individualReactionHashes); @@ -268,146 +284,82 @@ namespace gridfire::reaction { return XXHash64::hash(data, sizeInBytes, seed); } - REACLIBReaction::REACLIBReaction( - const std::string_view id, - const std::string_view peName, - const int chapter, - const std::vector &reactants, - const std::vector &products, - const double qValue, - const std::string_view label, - const REACLIBRateCoefficientSet &sets, - const bool reverse) : - Reaction(id, qValue, reactants, products, reverse), - m_peName(peName), - m_chapter(chapter), - m_sourceLabel(label), - m_rateCoefficients(sets) {} - std::unique_ptr REACLIBReaction::clone() const { - return std::make_unique(*this); - } - - - double REACLIBReaction::calculate_rate(const double T9) const { - return calculate_rate(T9); - } - - CppAD::AD REACLIBReaction::calculate_rate(const CppAD::AD T9) const { - return calculate_rate>(T9); - } - - REACLIBReactionSet::REACLIBReactionSet(std::vector reactions) : - ReactionSet(std::vector>()) { - // Convert REACLIBReaction to unique_ptr and store in m_reactions - m_reactions.reserve(reactions.size()); - m_reactionNameMap.reserve(reactions.size()); - for (auto& reaction : reactions) { - m_reactions.emplace_back(std::make_unique(std::move(reaction))); - m_reactionNameMap.emplace(std::string(reaction.id()), m_reactions.back().get()); - } - } - - std::unordered_set REACLIBReactionSet::peNames() const { - std::unordered_set peNames; - for (const auto& reactionPtr: m_reactions) { - if (const auto* reaction = dynamic_cast(reactionPtr.get())) { - peNames.insert(std::string(reaction->peName())); - } else { - // LOG_ERROR(m_logger, "Failed to cast Reaction to REACLIBReaction in REACLIBReactionSet::peNames()."); - throw std::runtime_error("Failed to cast Reaction to REACLIBReaction in REACLIBReactionSet::peNames(). This should not happen, please check your reaction setup."); - } - } - return peNames; - } - - REACLIBLogicalReaction::REACLIBLogicalReaction(const std::vector& reactants) : + LogicalReaction::LogicalReaction(const std::vector& reactants) : Reaction(reactants.front().peName(), - reactants.front().qValue(), + reactants.front().peName(), + reactants.front().chapter(), reactants.front().reactants(), reactants.front().products(), - reactants.front().is_reverse()), - m_chapter(reactants.front().chapter()) { + reactants.front().qValue(), + reactants.front().sourceLabel(), + reactants.front().rateCoefficients(), + reactants.front().is_reverse()) { m_sources.reserve(reactants.size()); m_rates.reserve(reactants.size()); for (const auto& reaction : reactants) { - if (std::abs(reaction.qValue() - m_qValue) > 1e-6) { - LOG_ERROR(m_logger, "REACLIBLogicalReaction constructed with reactions having different Q-values. Expected {} got {}.", m_qValue, reaction.qValue()); - throw std::runtime_error("REACLIBLogicalReaction constructed with reactions having different Q-values. Expected " + std::to_string(m_qValue) + " got " + std::to_string(reaction.qValue()) + "."); + if (std::abs(std::abs(reaction.qValue()) - std::abs(m_qValue)) > 1e-6) { + LOG_ERROR( + m_logger, + "LogicalReaction constructed with reactions having different Q-values. Expected {} got {}.", + m_qValue, + reaction.qValue() + ); + m_logger -> flush_log(); + throw std::runtime_error("LogicalReaction constructed with reactions having different Q-values. Expected " + std::to_string(m_qValue) + " got " + std::to_string(reaction.qValue()) + " (difference : " + std::to_string(std::abs(reaction.qValue() - m_qValue)) + ")."); } m_sources.push_back(std::string(reaction.sourceLabel())); m_rates.push_back(reaction.rateCoefficients()); } } - REACLIBLogicalReaction::REACLIBLogicalReaction(const REACLIBReaction& reaction) : - Reaction(reaction.peName(), - reaction.qValue(), - reaction.reactants(), - reaction.products(), - reaction.is_reverse()), - m_chapter(reaction.chapter()) { - m_sources.reserve(1); - m_rates.reserve(1); - m_sources.push_back(std::string(reaction.sourceLabel())); - m_rates.push_back(reaction.rateCoefficients()); - } - - void REACLIBLogicalReaction::add_reaction(const REACLIBReaction& reaction) { + void LogicalReaction::add_reaction(const Reaction& reaction) { if (reaction.peName() != m_id) { - LOG_ERROR(m_logger, "Cannot add reaction with different peName to REACLIBLogicalReaction. Expected {} got {}.", m_id, reaction.peName()); - throw std::runtime_error("Cannot add reaction with different peName to REACLIBLogicalReaction. Expected " + std::string(m_id) + " got " + std::string(reaction.peName()) + "."); + LOG_ERROR(m_logger, "Cannot add reaction with different peName to LogicalReaction. Expected {} got {}.", m_id, reaction.peName()); + m_logger -> flush_log(); + throw std::runtime_error("Cannot add reaction with different peName to LogicalReaction. Expected " + std::string(m_id) + " got " + std::string(reaction.peName()) + "."); } for (const auto& source : m_sources) { if (source == reaction.sourceLabel()) { - LOG_ERROR(m_logger, "Cannot add reaction with duplicate source label {} to REACLIBLogicalReaction.", reaction.sourceLabel()); - throw std::runtime_error("Cannot add reaction with duplicate source label " + std::string(reaction.sourceLabel()) + " to REACLIBLogicalReaction."); + LOG_ERROR(m_logger, "Cannot add reaction with duplicate source label {} to LogicalReaction.", reaction.sourceLabel()); + m_logger -> flush_log(); + throw std::runtime_error("Cannot add reaction with duplicate source label " + std::string(reaction.sourceLabel()) + " to LogicalReaction."); } } if (std::abs(reaction.qValue() - m_qValue) > 1e-6) { - LOG_ERROR(m_logger, "REACLIBLogicalReaction constructed with reactions having different Q-values. Expected {} got {}.", m_qValue, reaction.qValue()); - throw std::runtime_error("REACLIBLogicalReaction constructed with reactions having different Q-values. Expected " + std::to_string(m_qValue) + " got " + std::to_string(reaction.qValue()) + "."); + LOG_ERROR(m_logger, "LogicalReaction constructed with reactions having different Q-values. Expected {} got {}.", m_qValue, reaction.qValue()); + m_logger -> flush_log(); + throw std::runtime_error("LogicalReaction constructed with reactions having different Q-values. Expected " + std::to_string(m_qValue) + " got " + std::to_string(reaction.qValue()) + "."); } m_sources.push_back(std::string(reaction.sourceLabel())); m_rates.push_back(reaction.rateCoefficients()); } - std::unique_ptr REACLIBLogicalReaction::clone() const { - return std::make_unique(*this); - } - - double REACLIBLogicalReaction::calculate_rate(const double T9) const { + double LogicalReaction::calculate_rate(const double T9) const { return calculate_rate(T9); } - CppAD::AD REACLIBLogicalReaction::calculate_rate(const CppAD::AD T9) const { + CppAD::AD LogicalReaction::calculate_rate(const CppAD::AD T9) const { return calculate_rate>(T9); } - REACLIBLogicalReactionSet::REACLIBLogicalReactionSet(const REACLIBReactionSet &reactionSet) : - ReactionSet(std::vector>()) { + LogicalReactionSet::LogicalReactionSet(const ReactionSet &reactionSet) : + ReactionSet(std::vector()) { - std::unordered_map> grouped_reactions; + std::unordered_map> grouped_reactions; - for (const auto& reaction_ptr : reactionSet) { - if (const auto* reaclib_reaction = dynamic_cast(reaction_ptr.get())) { - grouped_reactions[reaclib_reaction->peName()].push_back(*reaclib_reaction); - } + for (const auto& reaction : reactionSet) { + grouped_reactions[reaction.peName()].push_back(reaction); } m_reactions.reserve(grouped_reactions.size()); m_reactionNameMap.reserve(grouped_reactions.size()); - for (const auto& [peName, reactions_for_peName] : grouped_reactions) { - m_peNames.insert(std::string(peName)); - auto logical_reaction = std::make_unique(reactions_for_peName); - m_reactionNameMap.emplace(logical_reaction->id(), logical_reaction.get()); + for (const auto &reactions_for_peName: grouped_reactions | std::views::values) { + LogicalReaction logical_reaction(reactions_for_peName); + m_reactionNameMap.emplace(logical_reaction.id(), logical_reaction); m_reactions.push_back(std::move(logical_reaction)); } } - - std::unordered_set REACLIBLogicalReactionSet::peNames() const { - return m_peNames; - } } namespace std { diff --git a/src/network/lib/solver/solver.cpp b/src/network/lib/solver/solver.cpp index 620787de..970e5188 100644 --- a/src/network/lib/solver/solver.cpp +++ b/src/network/lib/solver/solver.cpp @@ -21,6 +21,16 @@ namespace gridfire::solver { NetOut QSENetworkSolver::evaluate(const NetIn &netIn) { + // --- Use the policy to decide whether to update the view --- + if (shouldUpdateView(netIn)) { + LOG_DEBUG(m_logger, "Solver update policy triggered, network view updating..."); + m_engine.update(netIn); + LOG_DEBUG(m_logger, "Network view updated!"); + + m_lastSeenConditions = netIn; + m_isViewInitialized = true; + } + m_engine.generateJacobianMatrix(netIn.MolarAbundance(), netIn.temperature / 1e9, netIn.density); using state_type = boost::numeric::ublas::vector; using namespace boost::numeric::odeint; @@ -124,10 +134,19 @@ namespace gridfire::solver { for (size_t i = 0; i < m_engine.getNetworkSpecies().size(); ++i) { const auto& species = m_engine.getNetworkSpecies()[i]; - const double timescale = speciesTimescale[species]; + const double network_timescale = speciesTimescale.at(species); const double abundance = Y[i]; - if (std::isinf(timescale) || abundance < abundanceCutoff || timescale <= timescaleCutoff) { + double decay_timescale = std::numeric_limits::infinity(); + const double half_life = species.halfLife(); + if (half_life > 0 && !std::isinf(half_life)) { + constexpr double LN2 = 0.69314718056; + decay_timescale = half_life / LN2; + } + + const double final_timescale = std::min(network_timescale, decay_timescale); + + if (std::isinf(final_timescale) || abundance < abundanceCutoff || final_timescale <= timescaleCutoff) { QSESpeciesIndices.push_back(i); } else { dynamicSpeciesIndices.push_back(i); @@ -171,59 +190,13 @@ namespace gridfire::solver { const double rho, const dynamicQSESpeciesIndices &indices ) const { - std::vector Y_dyn; - Eigen::VectorXd Y_qse_initial(indices.QSESpeciesIndices.size()); - for (const auto& i : indices.dynamicSpeciesIndices) { Y_dyn.push_back(Y[i]); } - for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) { - Y_qse_initial(i) = Y[indices.QSESpeciesIndices[i]]; - if (Y_qse_initial(i) < 1.0e-99) { Y_qse_initial(i) = 1.0e-99; } - } - - Eigen::VectorXd v_qse = Y_qse_initial.array().log(); - - EigenFunctor qse_problem(m_engine, Y_dyn, indices.dynamicSpeciesIndices, indices.QSESpeciesIndices, T9, rho); - LOG_INFO(m_logger, "--- QSE Pre-Solve Diagnostics ---"); - - Eigen::VectorXd f_initial(indices.QSESpeciesIndices.size()); - qse_problem(v_qse, f_initial); - LOG_INFO(m_logger, "Initial Guess ||f||: {:0.4e}", f_initial.norm()); - - Eigen::MatrixXd J_initial(indices.QSESpeciesIndices.size(), indices.QSESpeciesIndices.size()); - qse_problem.df(v_qse, J_initial); - const Eigen::JacobiSVD svd(J_initial); - double cond = svd.singularValues().maxCoeff() / svd.singularValues().minCoeff(); - LOG_INFO(m_logger, "Initial Jacobian Condition Number: {:0.4e}", cond); - LOG_INFO(m_logger, "Starting QSE solve..."); - - Eigen::HybridNonLinearSolver> solver(qse_problem); - solver.parameters.xtol = 1.0e-8; // Set tolerance - - // 5. Run the solver. It will modify v_qse in place. - const int eigenStatus = solver.solve(v_qse); - - // 6. Check for convergence and return the result - if(eigenStatus != Eigen::Success) { - LOG_WARNING(m_logger, "--- QSE SOLVER FAILED ---"); - LOG_WARNING(m_logger, "Eigen status code: {}", eigenStatus); - LOG_WARNING(m_logger, "Iterations performed: {}", solver.iter); - - // Log the final state that caused the failure - Eigen::VectorXd Y_qse_final_fail = v_qse.array().exp(); - for(long i=0; i Y_qse[{}]: {:0.4e}", i, v_qse(i), i, Y_qse_final_fail(i)); - } - - // Log the residual at the final state - Eigen::VectorXd f_final(indices.QSESpeciesIndices.size()); - qse_problem(v_qse, f_final); - LOG_WARNING(m_logger, "Final ||f||: {:0.4e}", f_final.norm()); - - throw std::runtime_error("Eigen QSE solver did not converge."); - } - LOG_INFO(m_logger, "Eigen QSE solver converged in {} iterations.", solver.iter); - - return v_qse.array().exp(); + LOG_TRACE_L1(m_logger, "Calculating steady state abundances for QSE species..."); + LOG_WARNING(m_logger, "QSE solver logic not yet implemented, assuming all QSE species have 0 abundance."); + // --- Prepare the QSE species vector --- + Eigen::VectorXd v_qse(indices.QSESpeciesIndices.size()); + v_qse.setZero(); + return v_qse.array(); } NetOut QSENetworkSolver::initializeNetworkWithShortIgnition(const NetIn &netIn) const { @@ -265,6 +238,48 @@ namespace gridfire::solver { return postIgnition; } + bool QSENetworkSolver::shouldUpdateView(const NetIn &conditions) const { + // Policy 1: If the view has never been initialized, we must update. + if (!m_isViewInitialized) { + return true; + } + + // Policy 2: Check for significant relative change in temperature. + // Reaction rates are exponentially sensitive to temperature, so we use a tight threshold. + const double temp_threshold = m_config.get("gridfire:solver:policy:temp_threshold", 0.05); // 5% + const double temp_relative_change = std::abs(conditions.temperature - m_lastSeenConditions.temperature) / m_lastSeenConditions.temperature; + if (temp_relative_change > temp_threshold) { + LOG_DEBUG(m_logger, "Temperature changed by {:.1f}%, triggering view update.", temp_relative_change * 100); + return true; + } + + // Policy 3: Check for significant relative change in density. + const double rho_threshold = m_config.get("gridfire:solver:policy:rho_threshold", 0.10); // 10% + const double rho_relative_change = std::abs(conditions.density - m_lastSeenConditions.density) / m_lastSeenConditions.density; + if (rho_relative_change > rho_threshold) { + LOG_DEBUG(m_logger, "Density changed by {:.1f}%, triggering view update.", rho_relative_change * 100); + return true; + } + + // Policy 4: Check for fuel depletion. + // If a primary fuel source changes significantly, the network structure may change. + const double fuel_threshold = m_config.get("gridfire:solver:policy:fuel_threshold", 0.15); // 15% + + // Example: Check hydrogen abundance + const double h1_old = m_lastSeenConditions.composition.getMassFraction("H-1"); + const double h1_new = conditions.composition.getMassFraction("H-1"); + if (h1_old > 1e-12) { // Avoid division by zero + const double h1_relative_change = std::abs(h1_new - h1_old) / h1_old; + if (h1_relative_change > fuel_threshold) { + LOG_DEBUG(m_logger, "H-1 mass fraction changed by {:.1f}%, triggering view update.", h1_relative_change * 100); + return true; + } + } + + // If none of the above conditions are met, the current view is still good enough. + return false; + } + void QSENetworkSolver::RHSFunctor::operator()( const boost::numeric::ublas::vector &YDynamic, boost::numeric::ublas::vector &dYdtDynamic, @@ -295,6 +310,7 @@ namespace gridfire::solver { namespace odeint = boost::numeric::odeint; using fourdst::composition::Composition; + const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) const unsigned long numSpecies = m_engine.getNetworkSpecies().size(); @@ -363,7 +379,7 @@ namespace gridfire::solver { const std::vector y(Y.begin(), m_numSpecies + Y.begin()); auto [dydt, eps] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho); dYdt.resize(m_numSpecies + 1); - std::ranges::copy(dydt, dydt.begin()); + std::ranges::copy(dydt, dYdt.begin()); dYdt(m_numSpecies) = eps; } @@ -373,10 +389,10 @@ namespace gridfire::solver { double t, boost::numeric::ublas::vector &dfdt ) const { - J.resize(m_numSpecies + 1, m_numSpecies + 1); + J.resize(m_numSpecies+1, m_numSpecies+1); J.clear(); - for (int i = 0; i < m_numSpecies + 1; ++i) { - for (int j = 0; j < m_numSpecies + 1; ++j) { + for (int i = 0; i < m_numSpecies; ++i) { + for (int j = 0; j < m_numSpecies; ++j) { J(i, j) = m_engine.getJacobianMatrixEntry(i, j); } } diff --git a/src/network/meson.build b/src/network/meson.build index 5b68487f..91a3e9e3 100644 --- a/src/network/meson.build +++ b/src/network/meson.build @@ -3,7 +3,9 @@ network_sources = files( 'lib/network.cpp', 'lib/engine/engine_approx8.cpp', 'lib/engine/engine_graph.cpp', + 'lib/engine/engine_adaptive.cpp', 'lib/reaction/reaction.cpp', + 'lib/reaction/reaclib.cpp', 'lib/solver/solver.cpp', ) @@ -13,7 +15,6 @@ dependencies = [ const_dep, config_dep, composition_dep, - reaclib_reactions_dep, cppad_dep, log_dep, xxhash_dep, @@ -38,9 +39,12 @@ network_dep = declare_dependency( network_headers = files( 'include/gridfire/network.h', 'include/gridfire/engine/engine_abstract.h', + 'include/gridfire/engine/engine_view_abstract.h', 'include/gridfire/engine/engine_approx8.h', 'include/gridfire/engine/engine_graph.h', + 'include/gridfire/engine/engine_adaptive.h', 'include/gridfire/reaction/reaction.h', + 'include/gridfire/reaction/reaclib.h', 'include/gridfire/solver/solver.h', ) install_headers(network_headers, subdir : 'gridfire')