diff --git a/src/network/include/gridfire/engine/engine_abstract.h b/src/network/include/gridfire/engine/engine_abstract.h index 118638c7..4eb0da03 100644 --- a/src/network/include/gridfire/engine/engine_abstract.h +++ b/src/network/include/gridfire/engine/engine_abstract.h @@ -106,7 +106,7 @@ namespace gridfire { * time derivatives of all species and the specific nuclear energy generation * rate for the current state. */ - [[nodiscard]] virtual StepDerivatives calculateRHSAndEnergy( + [[nodiscard]] virtual std::expected, expectations::StaleEngineError> calculateRHSAndEnergy( const std::vector& Y, double T9, double rho @@ -217,6 +217,8 @@ namespace gridfire { */ [[nodiscard]] virtual const reaction::LogicalReactionSet& getNetworkReactions() const = 0; + virtual void setNetworkReactions(const reaction::LogicalReactionSet& reactions) = 0; + /** * @brief Compute timescales for all species in the network. * @@ -228,7 +230,13 @@ namespace gridfire { * This method estimates the timescale for abundance change of each species, * which can be used for timestep control, diagnostics, and reaction network culling. */ - [[nodiscard]] virtual std::unordered_map getSpeciesTimescales( + [[nodiscard]] virtual std::expected, expectations::StaleEngineError> getSpeciesTimescales( + const std::vector& Y, + double T9, + double rho + ) const = 0; + + [[nodiscard]] virtual std::expected, expectations::StaleEngineError> getSpeciesDestructionTimescales( const std::vector& Y, double T9, double rho diff --git a/src/network/include/gridfire/engine/engine_graph.h b/src/network/include/gridfire/engine/engine_graph.h index 6df2ce24..f48bc310 100644 --- a/src/network/include/gridfire/engine/engine_graph.h +++ b/src/network/include/gridfire/engine/engine_graph.h @@ -144,7 +144,7 @@ namespace gridfire { * * @see StepDerivatives */ - [[nodiscard]] StepDerivatives calculateRHSAndEnergy( + [[nodiscard]] std::expected, expectations::StaleEngineError> calculateRHSAndEnergy( const std::vector& Y, const double T9, const double rho @@ -215,6 +215,8 @@ namespace gridfire { */ [[nodiscard]] const reaction::LogicalReactionSet& getNetworkReactions() const override; + void setNetworkReactions(const reaction::LogicalReactionSet& reactions) override; + /** * @brief Gets an entry from the previously generated Jacobian matrix. * @@ -268,7 +270,13 @@ namespace gridfire { * 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( + [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesTimescales( + const std::vector& Y, + double T9, + double rho + ) const override; + + [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesDestructionTimescales( const std::vector& Y, double T9, double rho @@ -511,7 +519,7 @@ namespace gridfire { * to store the partial derivatives of the right-hand side of the ODE * with respect to the species abundances. */ - void reserveJacobianMatrix(); + void reserveJacobianMatrix() const; /** * @brief Records the AD tape for the right-hand side of the ODE. @@ -795,43 +803,12 @@ namespace gridfire { ); const T molarReactionFlow = forwardMolarReactionFlow - reverseMolarFlow; // Net molar reaction flow - // std::stringstream ss; - // ss << "Forward: " << forwardMolarReactionFlow - // << ", Reverse: " << reverseMolarFlow - // << ", Net: " << molarReactionFlow; - // LOG_TRACE_L3( - // m_logger, - // "Reaction: {}, {}", - // reaction.peName(), - // ss.str() - // ); - // std::cout << "Forward molar flow for reaction " << reaction.peName() << ": " << forwardMolarReactionFlow << std::endl; - // std::cout << "Reverse molar flow for reaction " << reaction.peName() << ": " << reverseMolarFlow << std::endl; - // std::cout << "Net molar flow for reaction " << reaction.peName() << ": " << molarReactionFlow << std::endl; // 3. Use the rate to update all relevant species derivatives (dY/dt) for (size_t speciesIndex = 0; speciesIndex < m_networkSpecies.size(); ++speciesIndex) { - const auto& species = m_networkSpecies[speciesIndex]; const T nu_ij = static_cast(m_stoichiometryMatrix(speciesIndex, reactionIndex)); - // bool taping = false; - // if constexpr (std::is_same_v>) { - // taping = true; - // } - // if (species.name() == "F-17" && nu_ij != static_cast(0.0)) { - // T f17dydt = result.dydt[speciesIndex] + threshold_flag * nu_ij * molarReactionFlow; - // std::string tstring = taping ? "During Taping" : "During Evaluation"; - // std::cout << "(" << tstring << ") F-17 Debugging. For Reaction " << reaction.id() << " (" << reaction.peName() << "): " - // << "nu_ij: " << nu_ij << ", molarReactionFlow: " << molarReactionFlow - // << ", threshold_flag: " << threshold_flag << ", dY/dt: " << f17dydt << ", Y: " << Y[speciesIndex] - // << std::endl; - // isF17 = true; - // } result.dydt[speciesIndex] += threshold_flag * nu_ij * molarReactionFlow; } - // if (isF17) { - // std::cout << "=========================================================" << std::endl; - // isF17 = false; - // } } T massProductionRate = static_cast(0.0); // [mol][s^-1] diff --git a/src/network/include/gridfire/engine/types/reporting.h b/src/network/include/gridfire/engine/types/reporting.h index 52a55c3a..8d74cf9b 100644 --- a/src/network/include/gridfire/engine/types/reporting.h +++ b/src/network/include/gridfire/engine/types/reporting.h @@ -33,7 +33,7 @@ namespace gridfire { friend std::ostream& operator<<(std::ostream& os, const PrimingReport& report) { std::stringstream ss; - std::string successStr = report.success ? "true" : "false"; + const std::string successStr = report.success ? "true" : "false"; ss << "PrimingReport(success=" << successStr << ", status=" << PrimingReportStatusStrings[report.status] << ")"; return os << ss.str(); diff --git a/src/network/include/gridfire/engine/views/engine_adaptive.h b/src/network/include/gridfire/engine/views/engine_adaptive.h index b10898ad..997c99e1 100644 --- a/src/network/include/gridfire/engine/views/engine_adaptive.h +++ b/src/network/include/gridfire/engine/views/engine_adaptive.h @@ -102,7 +102,7 @@ namespace gridfire { * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). * @see AdaptiveEngineView::update() */ - [[nodiscard]] StepDerivatives calculateRHSAndEnergy( + [[nodiscard]] std::expected, expectations::StaleEngineError> calculateRHSAndEnergy( const std::vector &Y_culled, const double T9, const double rho @@ -205,6 +205,8 @@ namespace gridfire { */ [[nodiscard]] const reaction::LogicalReactionSet& getNetworkReactions() const override; + void setNetworkReactions(const reaction::LogicalReactionSet& reactions) override; + /** * @brief Computes timescales for all active species in the network. * @@ -218,12 +220,18 @@ namespace gridfire { * * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). */ - [[nodiscard]] std::unordered_map getSpeciesTimescales( + [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesTimescales( const std::vector &Y_culled, double T9, double rho ) const override; + [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesDestructionTimescales( + const std::vector &Y, + double T9, + double rho + ) const override; + /** * @brief Gets the base engine. * @return A const reference to the base engine. @@ -443,7 +451,6 @@ namespace gridfire { typedef std::pair, std::unordered_set> RescueSet; [[nodiscard]] RescueSet rescueEdgeSpeciesDestructionChannel( - const std::vector& allFlows, const std::vector& Y_full, const double T9, const double rho, diff --git a/src/network/include/gridfire/engine/views/engine_defined.h b/src/network/include/gridfire/engine/views/engine_defined.h index 426d9b3a..b3fc8f44 100644 --- a/src/network/include/gridfire/engine/views/engine_defined.h +++ b/src/network/include/gridfire/engine/views/engine_defined.h @@ -37,7 +37,7 @@ namespace gridfire{ * * @throws std::runtime_error If the view is stale (i.e., `update()` has not been called after `setNetworkFile()`). */ - StepDerivatives calculateRHSAndEnergy( + std::expected, expectations::StaleEngineError> calculateRHSAndEnergy( const std::vector& Y_defined, const double T9, const double rho @@ -115,6 +115,8 @@ namespace gridfire{ * @throws std::runtime_error If the view is stale. */ const reaction::LogicalReactionSet& getNetworkReactions() const override; + + void setNetworkReactions(const reaction::LogicalReactionSet& reactions) override; /** * @brief Computes timescales for all active species in the network. * @@ -125,7 +127,13 @@ namespace gridfire{ * * @throws std::runtime_error If the view is stale. */ - std::unordered_map getSpeciesTimescales( + [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesTimescales( + const std::vector& Y_defined, + const double T9, + const double rho + ) const override; + + [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesDestructionTimescales( const std::vector& Y_defined, const double T9, const double rho @@ -243,6 +251,9 @@ namespace gridfire{ size_t mapViewToFullReactionIndex(size_t definedReactionIndex) const; void validateNetworkState() const; + + void collect(const std::vector& peNames); + }; class FileDefinedEngineView final: public DefinedEngineView { diff --git a/src/network/include/gridfire/engine/views/engine_multiscale.h b/src/network/include/gridfire/engine/views/engine_multiscale.h index 7528e87b..631029ea 100644 --- a/src/network/include/gridfire/engine/views/engine_multiscale.h +++ b/src/network/include/gridfire/engine/views/engine_multiscale.h @@ -66,7 +66,7 @@ namespace gridfire { [[nodiscard]] const std::vector & getNetworkSpecies() const override; - [[nodiscard]] StepDerivatives calculateRHSAndEnergy( + [[nodiscard]] std::expected, expectations::StaleEngineError> calculateRHSAndEnergy( const std::vector &Y_full, double T9, double rho @@ -99,7 +99,17 @@ namespace gridfire { [[nodiscard]] const reaction::LogicalReactionSet & getNetworkReactions() const override; - [[nodiscard]] std::unordered_map getSpeciesTimescales( + void setNetworkReactions( + const reaction::LogicalReactionSet &reactions + ) override; + + [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesTimescales( + const std::vector &Y, + double T9, + double rho + ) const override; + + [[nodiscard]] std::expected, expectations::StaleEngineError> getSpeciesDestructionTimescales( const std::vector &Y, double T9, double rho @@ -119,6 +129,13 @@ namespace gridfire { const DynamicEngine & getBaseEngine() const override; + std::vector> analyzeTimescalePoolConnectivity( + const std::vector> ×cale_pools, + const std::vector &Y, + double T9, + double rho + ) const; + void partitionNetwork( const std::vector& Y, double T9, @@ -162,6 +179,12 @@ namespace gridfire { bool is_in_equilibrium = false; ///< Flag set by flux analysis. std::set algebraic_indices; ///< Indices of algebraic species in this group. std::set seed_indices; ///< Indices of dynamic species in this group. + double mean_timescale; ///< Mean timescale of the group. + + bool operator<(const QSEGroup& other) const; + bool operator>(const QSEGroup& other) const; + bool operator==(const QSEGroup& other) const; + bool operator!=(const QSEGroup& other) const; friend std::ostream& operator<<(std::ostream& os, const QSEGroup& group) { os << "QSEGroup(species_indices: ["; @@ -175,7 +198,8 @@ namespace gridfire { } count = 0; os << "], is_in_equilibrium: " << group.is_in_equilibrium - << ", algebraic_indices: ["; + << ", mean_timescale: " << group.mean_timescale + << " s, algebraic_indices: ["; for (const auto& idx : group.algebraic_indices) { os << idx; if (count < group.algebraic_indices.size() - 1) { @@ -242,6 +266,7 @@ namespace gridfire { GenerateJacobianMatrix, CalculateMolarReactionFlow, GetSpeciesTimescales, + GetSpeciesDestructionTimescales, Other, All }; @@ -251,6 +276,7 @@ namespace gridfire { {operators::GenerateJacobianMatrix, "generateJacobianMatrix"}, {operators::CalculateMolarReactionFlow, "calculateMolarReactionFlow"}, {operators::GetSpeciesTimescales, "getSpeciesTimescales"}, + {operators::GetSpeciesDestructionTimescales, "getSpeciesDestructionTimescales"}, {operators::Other, "other"} }; @@ -262,6 +288,7 @@ namespace gridfire { {operators::GenerateJacobianMatrix, 0}, {operators::CalculateMolarReactionFlow, 0}, {operators::GetSpeciesTimescales, 0}, + {operators::GetSpeciesDestructionTimescales, 0}, {operators::Other, 0} }; @@ -351,11 +378,14 @@ namespace gridfire { double rho ) const; + std::unordered_map> buildConnectivityGraph( + const std::vector& species_pool + ) const; + std::vector constructCandidateGroups( - const std::vector>& timescale_pools, + const std::vector>& candidate_pools, const std::vector& Y, - double T9, - double rho + double T9, double rho ) const; }; } diff --git a/src/network/include/gridfire/exceptions/error_engine.h b/src/network/include/gridfire/exceptions/error_engine.h index 59633d01..46beddb4 100644 --- a/src/network/include/gridfire/exceptions/error_engine.h +++ b/src/network/include/gridfire/exceptions/error_engine.h @@ -7,14 +7,67 @@ namespace gridfire::exceptions { class EngineError : public std::exception {}; + class StaleEngineTrigger final : public EngineError { + public: + struct state { + double m_T9; + double m_rho; + std::vector m_Y; + double m_t; + int m_total_steps; + double m_eps_nuc; + }; + explicit StaleEngineTrigger(const state &s) + : m_state(s) {} + + const char* what() const noexcept override{ + return "Engine reports stale state. This means that the caller should trigger a update of the engine state before continuing with the integration. If you as an end user are seeing this error, it is likely a bug in the code that should be reported. Please provide the input parameters and the context in which this error occurred. Thank you for your help!"; + } + + state getState() const { + return m_state; + } + + size_t numSpecies() const { + return m_state.m_Y.size(); + } + + size_t totalSteps() const { + return m_state.m_total_steps; + } + + double energy() const { + return m_state.m_eps_nuc; + } + + double getMolarAbundance(const size_t index) const { + if (index > m_state.m_Y.size() - 1) { + throw std::out_of_range("Index out of bounds for molar abundance vector."); + } + return m_state.m_Y[index]; + } + + double temperature() const { + return m_state.m_T9 * 1e9; // Convert T9 back to Kelvin + } + + double density() const { + return m_state.m_rho; + } + private: + state m_state; + + }; + class StaleEngineError final : public EngineError { public: explicit StaleEngineError(const std::string& message) : m_message(message) {} - const char* what() const noexcept override{ + const char* what() const noexcept override { return m_message.c_str(); } + private: std::string m_message; }; @@ -31,4 +84,29 @@ namespace gridfire::exceptions { std::string m_message; }; + class NetworkResizedError final : public EngineError { + public: + explicit NetworkResizedError(const std::string& message) + : m_message(message) {} + + const char* what() const noexcept override { + return m_message.c_str(); + } + private: + std::string m_message; + }; + + class UnableToSetNetworkReactionsError final : public EngineError { + public: + explicit UnableToSetNetworkReactionsError(const std::string& message) + : m_message(message) {} + + const char* what() const noexcept override { + return m_message.c_str(); + } + + private: + std::string m_message; + }; + } \ No newline at end of file diff --git a/src/network/include/gridfire/expectations/expected_engine.h b/src/network/include/gridfire/expectations/expected_engine.h index d95daaa8..426b001b 100644 --- a/src/network/include/gridfire/expectations/expected_engine.h +++ b/src/network/include/gridfire/expectations/expected_engine.h @@ -10,6 +10,10 @@ namespace gridfire::expectations { STALE }; + enum class StaleEngineErrorTypes { + SYSTEM_RESIZED + }; + struct EngineError { std::string m_message; EngineErrorTypes type = EngineErrorTypes::FAILURE; @@ -30,5 +34,17 @@ namespace gridfire::expectations { struct StaleEngineError : EngineError { EngineErrorTypes type = EngineErrorTypes::STALE; + StaleEngineErrorTypes staleType; + + explicit StaleEngineError(StaleEngineErrorTypes staleType) : staleType(staleType) {} + + explicit operator std::string() const { + switch (staleType) { + case (StaleEngineErrorTypes::SYSTEM_RESIZED): + return "StaleEngineError: System resized, please update the engine."; + default: + return "StaleEngineError: Unknown stale error type."; + } + } }; } \ 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 93c8c496..d32e5c8f 100644 --- a/src/network/include/gridfire/solver/solver.h +++ b/src/network/include/gridfire/solver/solver.h @@ -10,25 +10,9 @@ #include "quill/Logger.h" -#include "unsupported/Eigen/NonLinearOptimization" // Required for LevenbergMarquardt - #include 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; ///< 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. @@ -68,316 +52,6 @@ namespace gridfire::solver { */ 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 - * @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 DynamicNetworkSolverStrategy { - public: - /** - * @brief Constructor for the QSENetworkSolver. - * @param engine The adaptive engine view to use for evaluating the network. - */ - using DynamicNetworkSolverStrategy::DynamicNetworkSolverStrategy; - - /** - * @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; - - /** - * @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. - quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance for trace and debug information. - - 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, - const std::vector& QSESpeciesIndices, - const Eigen::VectorXd& Y_QSE, - const double T9, - const double rho - ) : - m_engine(engine), - m_dynamicSpeciesIndices(dynamicSpeciesIndices), - m_QSESpeciesIndices(QSESpeciesIndices), - m_Y_QSE(Y_QSE), - 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, - double t - ) const; - - }; - - /** - * @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; ///< 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, - const std::vector& QSESpeciesIndices, - const double T9, - const double rho - ) : - m_engine(engine), - m_dynamicSpeciesIndices(dynamicSpeciesIndices), - m_QSESpeciesIndices(QSESpeciesIndices), - 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, - double t, - boost::numeric::ublas::vector& dfdt - ) 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; - enum { - InputsAtCompileTime = Eigen::Dynamic, - ValuesAtCompileTime = Eigen::Dynamic - }; - - DynamicEngine& m_engine; ///< The engine used to evaluate the network. - quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance for trace and debug information. - const std::vector& m_YFull; ///< The full, initial abundance vector - const std::vector& m_dynamicSpeciesIndices; ///< Indices of the dynamic species. - const std::vector& m_QSESpeciesIndices; ///< Indices of the QSE species. - const double m_T9; ///< Temperature in units of 10^9 K. - const double m_rho; ///< Density in g/cm^3. - const Eigen::VectorXd& m_YScale; ///< Scaling vector for the QSE species for asinh scaling. - - /** - * @brief Constructor for the EigenFunctor. - * @param engine The engine used to evaluate the network. - * @param YFull 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& YFull, - const std::vector& dynamicSpeciesIndices, - const std::vector& QSESpeciesIndices, - const double T9, - const double rho, - const Eigen::VectorXd& YScale - ) : - m_engine(engine), - m_YFull(YFull), - m_dynamicSpeciesIndices(dynamicSpeciesIndices), - m_QSESpeciesIndices(QSESpeciesIndices), - m_T9(T9), - m_rho(rho), - m_YScale(YScale) {} - - int values() const { return m_QSESpeciesIndices.size(); } - int inputs() const { return m_QSESpeciesIndices.size(); } - - /** - * @brief Calculates the residual vector for the QSE species. - * @param v_QSE_log Input vector of QSE species abundances (logarithmic). - * @param f_QSE Output vector of residuals. - * @return 0 for success. - */ - int operator()(const InputType& v_QSE_log, OutputType& f_QSE) const; - - /** - * @brief Calculates the Jacobian matrix for the QSE species. - * @param v_QSE_log Input vector of QSE species abundances (logarithmic). - * @param J_QSE Output Jacobian matrix. - * @return 0 for success. - */ - int df(const InputType& v_QSE_log, JacobianType& J_QSE) const; - }; - private: - 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. @@ -411,12 +85,20 @@ namespace gridfire::solver { * species abundances. It takes the current abundances as input and returns the * time derivatives. */ - struct RHSFunctor { + struct RHSManager { 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. - quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); ///< Logger instance. + + mutable double m_cached_time; + mutable std::optional> m_cached_result; + + mutable double m_last_observed_time = 0.0; ///< Last time the state was observed. + + + quill::Logger* m_logger = LogManager::getInstance().newFileLogger("integration.log", "GridFire"); ///< Logger instance. + mutable int m_num_steps = 0; + mutable double m_last_step_time = 1e-20; /** * @brief Constructor for the RHSFunctor. @@ -424,7 +106,7 @@ namespace gridfire::solver { * @param T9 Temperature in units of 10^9 K. * @param rho Density in g/cm^3. */ - RHSFunctor( + RHSManager( DynamicEngine& engine, const double T9, const double rho @@ -432,7 +114,7 @@ namespace gridfire::solver { m_engine(engine), m_T9(T9), m_rho(rho), - m_numSpecies(engine.getNetworkSpecies().size()) {} + m_cached_time(0) {} /** * @brief Calculates the time derivatives of the species abundances. @@ -445,6 +127,10 @@ namespace gridfire::solver { boost::numeric::ublas::vector& dYdt, double t ) const; + + void observe(const boost::numeric::ublas::vector& state, double t) const; + void compute_and_cache(const boost::numeric::ublas::vector& state, double t) const; + }; /** @@ -458,7 +144,6 @@ namespace gridfire::solver { 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. @@ -473,8 +158,7 @@ namespace gridfire::solver { ) : m_engine(engine), m_T9(T9), - m_rho(rho), - m_numSpecies(engine.getNetworkSpecies().size()) {} + m_rho(rho) {} /** * @brief Calculates the Jacobian matrix. @@ -493,97 +177,7 @@ namespace gridfire::solver { }; private: - quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance. - fourdst::config::Config& m_config = fourdst::config::Config::getInstance(); ///< Configuration instance. + quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); ///< Logger instance. + Config& m_config = Config::getInstance(); ///< Configuration instance. }; - - template - int QSENetworkSolver::EigenFunctor::operator()(const InputType &v_QSE, OutputType &f_QSE) const { - std::vector y = m_YFull; // Full vector of species abundances - Eigen::VectorXd Y_QSE = m_YScale.array() * v_QSE.array().sinh(); // Convert to physical abundances using asinh scaling - LOG_DEBUG( - m_logger, - "Calling QSENetworkSolver::EigenFunctor::operator() with QSE abundances {}", - [&]() -> std::string { - std::stringstream ss; - ss << std::scientific << std::setprecision(5); - for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { - ss << std::string(m_engine.getNetworkSpecies()[m_QSESpeciesIndices[i]].name()) << ": " << Y_QSE(i) << ", "; - } - return ss.str(); - }()); - - for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { - y[m_QSESpeciesIndices[i]] = Y_QSE(i); - } - - - const auto [dydt, specificEnergyRate] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho); - LOG_DEBUG( - m_logger, - "Calling QSENetworkSolver::EigenFunctor::operator() found QSE dydt {}", - [&]() -> std::string { - std::stringstream ss; - ss << std::scientific << std::setprecision(5); - for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { - ss << std::string(m_engine.getNetworkSpecies()[m_QSESpeciesIndices[i]].name()) << ": " << dydt[m_QSESpeciesIndices[i]] << ", "; - } - return ss.str(); - }()); - f_QSE.resize(m_QSESpeciesIndices.size()); - for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { - f_QSE(i) = dydt[m_QSESpeciesIndices[i]]; - } - return 0; // Success - } - - template - int QSENetworkSolver::EigenFunctor::df(const InputType& v_QSE, JacobianType& J_QSE) const { - std::vector y = m_YFull; - Eigen::VectorXd Y_QSE = m_YScale.array() * v_QSE.array().sinh(); - - for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { - y[m_QSESpeciesIndices[i]] = Y_QSE(i); - } - - m_engine.generateJacobianMatrix(y, m_T9, m_rho); - - J_QSE.resize(m_QSESpeciesIndices.size(), m_QSESpeciesIndices.size()); - for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { - for (size_t j = 0; j < m_QSESpeciesIndices.size(); ++j) { - J_QSE(i, j) = m_engine.getJacobianMatrixEntry(m_QSESpeciesIndices[i], m_QSESpeciesIndices[j]); - } - } - - std::string format_jacobian = [&]() -> std::string { - std::stringstream ss; - ss << std::scientific << std::setprecision(5); - for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { - for (size_t j = 0; j < m_QSESpeciesIndices.size(); ++j) { - ss << J_QSE(i, j) << " "; - } - ss << "\n"; - } - return ss.str(); - }(); - - - for (long j = 0; j < J_QSE.cols(); ++j) { - LOG_DEBUG( - m_logger, - "Jacobian ({},{}) before chain rule: {}", - j, j, J_QSE(j, j) - ); - const double dY_dv = m_YScale(j) * CppAD::cosh(v_QSE(j)); - J_QSE.col(j) *= dY_dv; // chain rule for log space transformation - LOG_DEBUG( - m_logger, - "Jacobian ({},{}) after chain rule: {} (dY/dv = {})", - j, j, J_QSE(j, j), dY_dv - ); - } - return 0; - } - - } \ No newline at end of file diff --git a/src/network/lib/engine/engine_graph.cpp b/src/network/lib/engine/engine_graph.cpp index 77535f6d..6f9ad5da 100644 --- a/src/network/lib/engine/engine_graph.cpp +++ b/src/network/lib/engine/engine_graph.cpp @@ -54,7 +54,7 @@ namespace gridfire { syncInternalMaps(); } - StepDerivatives GraphEngine::calculateRHSAndEnergy( + std::expected, expectations::StaleEngineError> GraphEngine::calculateRHSAndEnergy( const std::vector &Y, const double T9, const double rho @@ -91,8 +91,8 @@ namespace gridfire { const size_t n = m_rhsADFun.Domain(); const size_t m = m_rhsADFun.Range(); - std::vector select_domain(n, true); - std::vector select_range(m, true); + const std::vector select_domain(n, true); + const std::vector select_range(m, true); m_rhsADFun.subgraph_sparsity(select_domain, select_range, false, m_full_jacobian_sparsity_pattern); m_jac_work.clear(); @@ -148,10 +148,10 @@ namespace gridfire { } } - void GraphEngine::reserveJacobianMatrix() { + void GraphEngine::reserveJacobianMatrix() const { // The implementation of this function (and others) constrains this nuclear network to a constant temperature and density during // each evaluation. - size_t numSpecies = m_networkSpecies.size(); + const size_t numSpecies = m_networkSpecies.size(); m_jacobianMatrix.clear(); m_jacobianMatrix.resize(numSpecies, numSpecies, false); // Sparse matrix, no initial values LOG_TRACE_L2(m_logger, "Jacobian matrix resized to {} rows and {} columns.", @@ -171,6 +171,11 @@ namespace gridfire { return m_reactions; } + void GraphEngine::setNetworkReactions(const reaction::LogicalReactionSet &reactions) { + m_reactions = reactions; + syncInternalMaps(); + } + bool GraphEngine::involvesSpecies(const fourdst::atomic::Species& species) const { // Checks if a given species is present in the network's species map for efficient lookup. const bool found = m_networkSpeciesMap.contains(species.name()); @@ -259,7 +264,7 @@ namespace gridfire { const double T9 ) const { if (!m_useReverseReactions) { - LOG_TRACE_L2(m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id()); + LOG_TRACE_L3_LIMIT_EVERY_N(std::numeric_limits::max(), m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id()); return 0.0; // If reverse reactions are not used, return 0.0 } const double expFactor = std::exp(-reaction.qValue() / (m_constants.kB * T9)); @@ -279,7 +284,7 @@ namespace gridfire { } else { LOG_WARNING_LIMIT_EVERY_N(1000000, m_logger, "Reverse rate calculation for reactions with more than two reactants or products is not implemented (reaction id {}).", reaction.peName()); } - LOG_TRACE_L2(m_logger, "Calculated reverse rate for reaction '{}': {:.3E} at T9={:.3E}.", reaction.id(), reverseRate, T9); + LOG_TRACE_L2_LIMIT_EVERY_N(1000, m_logger, "Calculated reverse rate for reaction '{}': {:.3E} at T9={:.3E}.", reaction.id(), reverseRate, T9); return reverseRate; } @@ -352,7 +357,7 @@ namespace gridfire { const double CT = std::pow(massNumerator/massDenominator, 1.5) * (partitionFunctionNumerator/partitionFunctionDenominator); - double reverseRate = forwardRate * symmetryFactor * CT * expFactor; + const double reverseRate = forwardRate * symmetryFactor * CT * expFactor; if (!std::isfinite(reverseRate)) { return 0.0; // If the reverse rate is not finite, return 0.0 } @@ -366,7 +371,7 @@ namespace gridfire { const double reverseRate ) const { if (!m_useReverseReactions) { - LOG_TRACE_L2(m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate derivative of reaction '{}'.", reaction.id()); + LOG_TRACE_L3_LIMIT_EVERY_N(std::numeric_limits::max(), m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id()); return 0.0; // If reverse reactions are not used, return 0.0 } const double d_log_kFwd = reaction.calculate_forward_rate_log_derivative(T9); @@ -652,7 +657,7 @@ namespace gridfire { const double rho ) const { - LOG_TRACE_L1(m_logger, "Generating jacobian matrix for T9={}, rho={}..", T9, rho); + LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Generating jacobian matrix for T9={}, rho={}..", T9, rho); const size_t numSpecies = m_networkSpecies.size(); // 1. Pack the input variables into a vector for CppAD @@ -662,20 +667,20 @@ namespace gridfire { } adInput[numSpecies] = T9; // T9 adInput[numSpecies + 1] = rho; // rho - LOG_DEBUG( - m_logger, - "AD Input to jacobian {}", - [&]() -> std::string { - std::stringstream ss; - ss << std::scientific << std::setprecision(5); - for (size_t i = 0; i < adInput.size(); ++i) { - ss << adInput[i]; - if (i < adInput.size() - 1) { - ss << ", "; - } - } - return ss.str(); - }()); + // LOG_DEBUG( + // m_logger, + // "AD Input to jacobian {}", + // [&]() -> std::string { + // std::stringstream ss; + // ss << std::scientific << std::setprecision(5); + // for (size_t i = 0; i < adInput.size(); ++i) { + // ss << adInput[i]; + // if (i < adInput.size() - 1) { + // ss << ", "; + // } + // } + // return ss.str(); + // }()); // 2. Calculate the full jacobian const std::vector dotY = m_rhsADFun.Jacobian(adInput); @@ -690,32 +695,32 @@ namespace gridfire { } } } - LOG_DEBUG( - m_logger, - "Final Jacobian is:\n{}", - [&]() -> std::string { - std::stringstream ss; - ss << std::scientific << std::setprecision(5); - for (size_t i = 0; i < m_jacobianMatrix.size1(); ++i) { - ss << getNetworkSpecies()[i].name(); - if (i < m_jacobianMatrix.size1() - 1) { - ss << ", "; - } - } - ss << "\n"; - for (size_t i = 0; i < m_jacobianMatrix.size1(); ++i) { - ss << getNetworkSpecies()[i].name() << ": "; - for (size_t j = 0; j < m_jacobianMatrix.size2(); ++j) { - ss << m_jacobianMatrix(i, j); - if (j < m_jacobianMatrix.size2() - 1) { - ss << ", "; - } - } - ss << "\n"; - } - return ss.str(); - }()); - LOG_TRACE_L1(m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); + // LOG_DEBUG( + // m_logger, + // "Final Jacobian is:\n{}", + // [&]() -> std::string { + // std::stringstream ss; + // ss << std::scientific << std::setprecision(5); + // for (size_t i = 0; i < m_jacobianMatrix.size1(); ++i) { + // ss << getNetworkSpecies()[i].name(); + // if (i < m_jacobianMatrix.size1() - 1) { + // ss << ", "; + // } + // } + // ss << "\n"; + // for (size_t i = 0; i < m_jacobianMatrix.size1(); ++i) { + // ss << getNetworkSpecies()[i].name() << ": "; + // for (size_t j = 0; j < m_jacobianMatrix.size2(); ++j) { + // ss << m_jacobianMatrix(i, j); + // if (j < m_jacobianMatrix.size2() - 1) { + // ss << ", "; + // } + // } + // ss << "\n"; + // } + // return ss.str(); + // }()); + LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); } void GraphEngine::generateJacobianMatrix( @@ -895,7 +900,7 @@ namespace gridfire { LOG_TRACE_L1(m_logger, "Successfully exported network graph to {}", filename); } - std::unordered_map GraphEngine::getSpeciesTimescales( + std::expected, expectations::StaleEngineError> GraphEngine::getSpeciesTimescales( const std::vector &Y, const double T9, const double rho @@ -914,8 +919,41 @@ namespace gridfire { return speciesTimescales; } + std::expected, expectations::StaleEngineError> GraphEngine::getSpeciesDestructionTimescales( + const std::vector &Y, + const double T9, + const double rho + ) const { + auto [dydt, _] = calculateAllDerivatives(Y, T9, rho); + std::unordered_map speciesDestructionTimescales; + speciesDestructionTimescales.reserve(m_networkSpecies.size()); + for (const auto& species : m_networkSpecies) { + double netDestructionFlow = 0.0; + for (const auto& reaction : m_reactions) { + if (reaction.stoichiometry(species) < 0) { + const double flow = calculateMolarReactionFlow(reaction, Y, T9, rho); + netDestructionFlow += flow; + } + } + double timescale = std::numeric_limits::infinity(); + if (netDestructionFlow != 0.0) { + timescale = Y[getSpeciesIndex(species)] / netDestructionFlow; + } + speciesDestructionTimescales.emplace(species, timescale); + } + return speciesDestructionTimescales; + } + fourdst::composition::Composition GraphEngine::update(const NetIn &netIn) { - return netIn.composition; + fourdst::composition::Composition baseUpdatedComposition = netIn.composition; + for (const auto& species : m_networkSpecies) { + if (!netIn.composition.contains(species)) { + baseUpdatedComposition.registerSpecies(species); + baseUpdatedComposition.setMassFraction(species, 0.0); + } + } + baseUpdatedComposition.finalize(false); + return baseUpdatedComposition; } bool GraphEngine::isStale(const NetIn &netIn) { diff --git a/src/network/lib/engine/procedures/priming.cpp b/src/network/lib/engine/procedures/priming.cpp index e899b1ae..295d4f4a 100644 --- a/src/network/lib/engine/procedures/priming.cpp +++ b/src/network/lib/engine/procedures/priming.cpp @@ -45,7 +45,7 @@ namespace gridfire { speciesToPrime.push_back(entry.isotope()); } } - LOG_INFO(logger, "Priming {} species in the network.", speciesToPrime.size()); + LOG_DEBUG(logger, "Priming {} species in the network.", speciesToPrime.size()); PrimingReport report; if (speciesToPrime.empty()) { @@ -57,7 +57,7 @@ namespace gridfire { const double T9 = netIn.temperature / 1e9; const double rho = netIn.density; - const auto initialDepth = engine.getDepth(); + const auto initialReactionSet = engine.getNetworkReactions(); report.status = PrimingReportStatus::FULL_SUCCESS; report.success = true; @@ -76,7 +76,7 @@ namespace gridfire { engine.rebuild(netIn.composition, NetworkBuildDepth::Full); for (const auto& primingSpecies : speciesToPrime) { - LOG_DEBUG(logger, "Priming species: {}", primingSpecies.name()); + LOG_TRACE_L3(logger, "Priming species: {}", primingSpecies.name()); // Create a temporary composition from the current internal state for the primer Composition tempComp; @@ -114,12 +114,12 @@ namespace gridfire { LOG_WARNING(logger, "Equilibrium mass fraction for {} is NaN. Setting to 0.0. This is likely not an issue. It probably originates from all reactions leading to creation and destruction being frozen out. In that case 0.0 should be a good approximation. Hint: This happens often when the network temperature is very the low. ", primingSpecies.name()); equilibriumMassFraction = 0.0; } - LOG_INFO(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction); + LOG_TRACE_L3(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction); const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho); if (dominantChannel) { - LOG_DEBUG(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->peName()); + LOG_TRACE_L3(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->peName()); double totalReactantMass = 0.0; for (const auto& reactant : dominantChannel->reactants()) { @@ -186,7 +186,7 @@ namespace gridfire { report.massFractionChanges.emplace_back(species, change); } - engine.rebuild(netIn.composition, initialDepth); + engine.setNetworkReactions(initialReactionSet); return report; } diff --git a/src/network/lib/engine/views/engine_adaptive.cpp b/src/network/lib/engine/views/engine_adaptive.cpp index f89578a4..c3afb0e1 100644 --- a/src/network/lib/engine/views/engine_adaptive.cpp +++ b/src/network/lib/engine/views/engine_adaptive.cpp @@ -6,6 +6,7 @@ #include "gridfire/network.h" +#include "gridfire/exceptions/error_engine.h" #include "quill/LogMacros.h" #include "quill/Logger.h" @@ -86,11 +87,21 @@ namespace gridfire { fourdst::composition::Composition AdaptiveEngineView::update(const NetIn &netIn) { fourdst::composition::Composition baseUpdatedComposition = m_baseEngine.update(netIn); + NetIn updatedNetIn = netIn; + + // for (const auto &entry: netIn.composition | std::views::values) { + // if (baseUpdatedComposition.contains(entry.isotope())) { + // updatedNetIn.composition.setMassFraction(entry.isotope(), baseUpdatedComposition.getMassFraction(entry.isotope())); + // } + // } + updatedNetIn.composition = baseUpdatedComposition; + + updatedNetIn.composition.finalize(false); LOG_TRACE_L1(m_logger, "Updating AdaptiveEngineView with new network input..."); std::vector Y_Full; - std::vector allFlows = calculateAllReactionFlows(netIn, Y_Full); + std::vector allFlows = calculateAllReactionFlows(updatedNetIn, Y_Full); double maxFlow = 0.0; @@ -101,24 +112,24 @@ namespace gridfire { } LOG_DEBUG(m_logger, "Maximum reaction flow rate in adaptive engine view: {:0.3E} [mol/s]", maxFlow); - const std::unordered_set reachableSpecies = findReachableSpecies(netIn); + const std::unordered_set reachableSpecies = findReachableSpecies(updatedNetIn); LOG_DEBUG(m_logger, "Found {} reachable species in adaptive engine view.", reachableSpecies.size()); const std::vector finalReactions = cullReactionsByFlow(allFlows, reachableSpecies, Y_Full, maxFlow); finalizeActiveSet(finalReactions); - // auto [rescuedReactions, rescuedSpecies] = rescueEdgeSpeciesDestructionChannel(allFlows, Y_Full, netIn.temperature/1e9, netIn.density, m_activeSpecies, m_activeReactions); - // - // for (const auto& reactionPtr : rescuedReactions) { - // m_activeReactions.add_reaction(*reactionPtr); - // } - // - // for (const auto& species : rescuedSpecies) { - // if (!std::ranges::contains(m_activeSpecies, species)) { - // m_activeSpecies.push_back(species); - // } - // } + auto [rescuedReactions, rescuedSpecies] = rescueEdgeSpeciesDestructionChannel(Y_Full, netIn.temperature/1e9, netIn.density, m_activeSpecies, m_activeReactions); + + for (const auto& reactionPtr : rescuedReactions) { + m_activeReactions.add_reaction(*reactionPtr); + } + + for (const auto& species : rescuedSpecies) { + if (!std::ranges::contains(m_activeSpecies, species)) { + m_activeSpecies.push_back(species); + } + } m_speciesIndexMap = constructSpeciesIndexMap(); m_reactionIndexMap = constructReactionIndexMap(); @@ -127,18 +138,18 @@ namespace gridfire { LOG_INFO(m_logger, "AdaptiveEngineView updated successfully with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size()); - return baseUpdatedComposition; + return updatedNetIn.composition; } bool AdaptiveEngineView::isStale(const NetIn &netIn) { - return m_isStale; + return m_isStale || m_baseEngine.isStale(netIn); } const std::vector & AdaptiveEngineView::getNetworkSpecies() const { return m_activeSpecies; } - StepDerivatives AdaptiveEngineView::calculateRHSAndEnergy( + std::expected, expectations::StaleEngineError> AdaptiveEngineView::calculateRHSAndEnergy( const std::vector &Y_culled, const double T9, const double rho @@ -147,8 +158,13 @@ namespace gridfire { const auto Y_full = mapCulledToFull(Y_culled); - const auto [dydt, nuclearEnergyGenerationRate] = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); + auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); + if (!result) { + return std::unexpected{result.error()}; + } + + const auto [dydt, nuclearEnergyGenerationRate] = result.value(); StepDerivatives culledResults; culledResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate; culledResults.dydt = mapFullToCulled(dydt); @@ -214,14 +230,26 @@ namespace gridfire { return m_activeReactions; } - std::unordered_map AdaptiveEngineView::getSpeciesTimescales( + void AdaptiveEngineView::setNetworkReactions(const reaction::LogicalReactionSet &reactions) { + LOG_CRITICAL(m_logger, "AdaptiveEngineView does not support setting network reactions directly. Use update() with NetIn instead. Perhaps you meant to call this on the base engine?"); + throw exceptions::UnableToSetNetworkReactionsError("AdaptiveEngineView does not support setting network reactions directly. Use update() with NetIn instead. Perhaps you meant to call this on the base engine?"); + } + + std::expected, expectations::StaleEngineError> AdaptiveEngineView::getSpeciesTimescales( const std::vector &Y_culled, const double T9, const double rho ) const { validateState(); const auto Y_full = mapCulledToFull(Y_culled); - const auto fullTimescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + + if (!result) { + return std::unexpected{result.error()}; + } + + const std::unordered_map fullTimescales = result.value(); + std::unordered_map culledTimescales; culledTimescales.reserve(m_activeSpecies.size()); @@ -234,6 +262,33 @@ namespace gridfire { } + std::expected, expectations::StaleEngineError> + AdaptiveEngineView::getSpeciesDestructionTimescales( + const std::vector &Y, + double T9, + double rho + ) const { + validateState(); + + const auto Y_full = mapCulledToFull(Y); + const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho); + + if (!result) { + return std::unexpected{result.error()}; + } + + const std::unordered_map destructionTimescales = result.value(); + + std::unordered_map culledTimescales; + culledTimescales.reserve(m_activeSpecies.size()); + for (const auto& active_species : m_activeSpecies) { + if (destructionTimescales.contains(active_species)) { + culledTimescales[active_species] = destructionTimescales.at(active_species); + } + } + return culledTimescales; + } + void AdaptiveEngineView::setScreeningModel(const screening::ScreeningType model) { m_baseEngine.setScreeningModel(model); } @@ -255,7 +310,7 @@ namespace gridfire { } int AdaptiveEngineView::getSpeciesIndex(const fourdst::atomic::Species &species) const { - auto it = std::find(m_activeSpecies.begin(), m_activeSpecies.end(), species); + const auto it = std::ranges::find(m_activeSpecies, species); if (it != m_activeSpecies.end()) { return static_cast(std::distance(m_activeSpecies.begin(), it)); } else { @@ -428,14 +483,18 @@ namespace gridfire { } AdaptiveEngineView::RescueSet AdaptiveEngineView::rescueEdgeSpeciesDestructionChannel( - const std::vector &allFlows, const std::vector &Y_full, const double T9, const double rho, const std::vector &activeSpecies, const reaction::LogicalReactionSet &activeReactions ) const { - const auto timescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + if (!result) { + LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state."); + throw exceptions::StaleEngineError("Failed to get species timescales"); + } + std::unordered_map timescales = result.value(); std::set onlyProducedSpecies; for (const auto& reaction : activeReactions) { const std::vector products = reaction.products(); @@ -521,7 +580,7 @@ namespace gridfire { } int count = 0; ss << "("; - for (const auto& [species, reaction] : reactionsToRescue) { + for (const auto &reaction : reactionsToRescue | std::views::values) { ss << reaction->id(); if (count < reactionsToRescue.size() - 1) { ss << ", "; diff --git a/src/network/lib/engine/views/engine_defined.cpp b/src/network/lib/engine/views/engine_defined.cpp index 1c80d70a..c00e7fa7 100644 --- a/src/network/lib/engine/views/engine_defined.cpp +++ b/src/network/lib/engine/views/engine_defined.cpp @@ -16,60 +16,7 @@ namespace gridfire { DefinedEngineView::DefinedEngineView(const std::vector& peNames, DynamicEngine& baseEngine) : m_baseEngine(baseEngine) { - m_activeReactions.clear(); - m_activeSpecies.clear(); - - std::unordered_set seenSpecies; - - const auto& fullNetworkReactionSet = m_baseEngine.getNetworkReactions(); - for (const auto& peName : peNames) { - if (!fullNetworkReactionSet.contains(peName)) { - LOG_ERROR(m_logger, "Reaction with name '{}' not found in the base engine's network reactions. Aborting...", peName); - m_logger->flush_log(); - throw std::runtime_error("Reaction with name '" + std::string(peName) + "' not found in the base engine's network reactions."); - } - auto reaction = fullNetworkReactionSet[peName]; - for (const auto& reactant : reaction.reactants()) { - if (!seenSpecies.contains(reactant)) { - seenSpecies.insert(reactant); - m_activeSpecies.push_back(reactant); - } - } - for (const auto& product : reaction.products()) { - if (!seenSpecies.contains(product)) { - seenSpecies.insert(product); - m_activeSpecies.push_back(product); - } - } - m_activeReactions.add_reaction(reaction); - } - LOG_TRACE_L1(m_logger, "DefinedEngineView built with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size()); - LOG_DEBUG(m_logger, "Active species: {}", [this]() -> std::string { - std::string result; - for (const auto& species : m_activeSpecies) { - result += std::string(species.name()) + ", "; - } - if (!result.empty()) { - result.pop_back(); // Remove last space - result.pop_back(); // Remove last comma - } - return result; - }()); - LOG_DEBUG(m_logger, "Active reactions: {}", [this]() -> std::string { - std::string result; - for (const auto& reaction : m_activeReactions) { - result += std::string(reaction.id()) + ", "; - } - if (!result.empty()) { - result.pop_back(); // Remove last space - result.pop_back(); // Remove last comma - } - return result; - }()); - m_speciesIndexMap = constructSpeciesIndexMap(); - m_reactionIndexMap = constructReactionIndexMap(); - m_isStale = false; - + collect(peNames); } const DynamicEngine & DefinedEngineView::getBaseEngine() const { @@ -80,7 +27,7 @@ namespace gridfire { return m_activeSpecies; } - StepDerivatives DefinedEngineView::calculateRHSAndEnergy( + std::expected, expectations::StaleEngineError> DefinedEngineView::calculateRHSAndEnergy( const std::vector &Y_defined, const double T9, const double rho @@ -88,8 +35,13 @@ namespace gridfire { validateNetworkState(); const auto Y_full = mapViewToFull(Y_defined); - const auto [dydt, nuclearEnergyGenerationRate] = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); + const auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); + if (!result) { + return std::unexpected{result.error()}; + } + + const auto [dydt, nuclearEnergyGenerationRate] = result.value(); StepDerivatives definedResults; definedResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate; definedResults.dydt = mapFullToView(dydt); @@ -159,7 +111,15 @@ namespace gridfire { return m_activeReactions; } - std::unordered_map DefinedEngineView::getSpeciesTimescales( + void DefinedEngineView::setNetworkReactions(const reaction::LogicalReactionSet &reactions) { + std::vector peNames; + for (const auto& reaction : reactions) { + peNames.push_back(std::string(reaction.id())); + } + collect(peNames); + } + + std::expected, expectations::StaleEngineError> DefinedEngineView::getSpeciesTimescales( const std::vector &Y_defined, const double T9, const double rho @@ -167,7 +127,11 @@ namespace gridfire { validateNetworkState(); const auto Y_full = mapViewToFull(Y_defined); - const auto fullTimescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + if (!result) { + return std::unexpected{result.error()}; + } + const auto& fullTimescales = result.value(); std::unordered_map definedTimescales; for (const auto& active_species : m_activeSpecies) { @@ -178,12 +142,38 @@ namespace gridfire { return definedTimescales; } + std::expected, expectations::StaleEngineError> + DefinedEngineView::getSpeciesDestructionTimescales( + const std::vector &Y_defined, + const double T9, + const double rho + ) const { + validateNetworkState(); + + const auto Y_full = mapViewToFull(Y_defined); + const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho); + + if (!result) { + return std::unexpected{result.error()}; + } + + const auto& destructionTimescales = result.value(); + + std::unordered_map definedTimescales; + for (const auto& active_species : m_activeSpecies) { + if (destructionTimescales.contains(active_species)) { + definedTimescales[active_species] = destructionTimescales.at(active_species); + } + } + return definedTimescales; + } + fourdst::composition::Composition DefinedEngineView::update(const NetIn &netIn) { return m_baseEngine.update(netIn); } bool DefinedEngineView::isStale(const NetIn &netIn) { - return false; + return m_baseEngine.isStale(netIn); } @@ -198,7 +188,7 @@ namespace gridfire { int DefinedEngineView::getSpeciesIndex(const Species &species) const { validateNetworkState(); - auto it = std::find(m_activeSpecies.begin(), m_activeSpecies.end(), species); + const auto it = std::ranges::find(m_activeSpecies, species); if (it != m_activeSpecies.end()) { return static_cast(std::distance(m_activeSpecies.begin(), it)); } else { @@ -224,7 +214,7 @@ namespace gridfire { } std::vector DefinedEngineView::constructSpeciesIndexMap() const { - LOG_TRACE_L1(m_logger, "Constructing species index map for DefinedEngineView..."); + LOG_TRACE_L3(m_logger, "Constructing species index map for DefinedEngineView..."); std::unordered_map fullSpeciesReverseMap; const auto& fullSpeciesList = m_baseEngine.getNetworkSpecies(); @@ -247,13 +237,13 @@ namespace gridfire { throw std::runtime_error("Species not found in full species map: " + std::string(active_species.name())); } } - LOG_TRACE_L1(m_logger, "Species index map constructed with {} entries.", speciesIndexMap.size()); + LOG_TRACE_L3(m_logger, "Species index map constructed with {} entries.", speciesIndexMap.size()); return speciesIndexMap; } std::vector DefinedEngineView::constructReactionIndexMap() const { - LOG_TRACE_L1(m_logger, "Constructing reaction index map for DefinedEngineView..."); + LOG_TRACE_L3(m_logger, "Constructing reaction index map for DefinedEngineView..."); // --- Step 1: Create a reverse map using the reaction's unique ID as the key. --- std::unordered_map fullReactionReverseMap; @@ -280,7 +270,7 @@ namespace gridfire { } } - LOG_TRACE_L1(m_logger, "Reaction index map constructed with {} entries.", reactionIndexMap.size()); + LOG_TRACE_L3(m_logger, "Reaction index map constructed with {} entries.", reactionIndexMap.size()); return reactionIndexMap; } @@ -328,6 +318,59 @@ namespace gridfire { } } + void DefinedEngineView::collect(const std::vector &peNames) { + std::unordered_set seenSpecies; + + const auto& fullNetworkReactionSet = m_baseEngine.getNetworkReactions(); + for (const auto& peName : peNames) { + if (!fullNetworkReactionSet.contains(peName)) { + LOG_ERROR(m_logger, "Reaction with name '{}' not found in the base engine's network reactions. Aborting...", peName); + m_logger->flush_log(); + throw std::runtime_error("Reaction with name '" + std::string(peName) + "' not found in the base engine's network reactions."); + } + auto reaction = fullNetworkReactionSet[peName]; + for (const auto& reactant : reaction.reactants()) { + if (!seenSpecies.contains(reactant)) { + seenSpecies.insert(reactant); + m_activeSpecies.push_back(reactant); + } + } + for (const auto& product : reaction.products()) { + if (!seenSpecies.contains(product)) { + seenSpecies.insert(product); + m_activeSpecies.push_back(product); + } + } + m_activeReactions.add_reaction(reaction); + } + LOG_TRACE_L3(m_logger, "DefinedEngineView built with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size()); + LOG_TRACE_L3(m_logger, "Active species: {}", [this]() -> std::string { + std::string result; + for (const auto& species : m_activeSpecies) { + result += std::string(species.name()) + ", "; + } + if (!result.empty()) { + result.pop_back(); // Remove last space + result.pop_back(); // Remove last comma + } + return result; + }()); + LOG_TRACE_L3(m_logger, "Active reactions: {}", [this]() -> std::string { + std::string result; + for (const auto& reaction : m_activeReactions) { + result += std::string(reaction.id()) + ", "; + } + if (!result.empty()) { + result.pop_back(); // Remove last space + result.pop_back(); // Remove last comma + } + return result; + }()); + m_speciesIndexMap = constructSpeciesIndexMap(); + m_reactionIndexMap = constructReactionIndexMap(); + m_isStale = false; + } + //////////////////////////////////////////// /// FileDefinedEngineView Implementation /// diff --git a/src/network/lib/engine/views/engine_multiscale.cpp b/src/network/lib/engine/views/engine_multiscale.cpp index f7d3ec63..7612c50e 100644 --- a/src/network/lib/engine/views/engine_multiscale.cpp +++ b/src/network/lib/engine/views/engine_multiscale.cpp @@ -10,11 +10,13 @@ #include #include +#include #include "quill/LogMacros.h" #include "quill/Logger.h" namespace { + using namespace fourdst::atomic; std::vector packCompositionToVector(const fourdst::composition::Composition& composition, const gridfire::GraphEngine& engine) { std::vector Y(engine.getNetworkSpecies().size(), 0.0); const auto& allSpecies = engine.getNetworkSpecies(); @@ -34,6 +36,115 @@ namespace { std::hash hashed; seed ^= hashed(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2); } + + std::vector> findConnectedComponentsBFS( + const std::unordered_map>& graph, + const std::vector& nodes + ) { + std::vector> components; + std::unordered_set visited; + + for (const size_t& start_node : nodes) { + if (visited.find(start_node) == visited.end()) { + std::vector current_component; + std::queue q; + + q.push(start_node); + visited.insert(start_node); + + while (!q.empty()) { + size_t u = q.front(); + q.pop(); + current_component.push_back(u); + + if (graph.count(u)) { + for (const auto& v : graph.at(u)) { + if (visited.find(v) == visited.end()) { + visited.insert(v); + q.push(v); + } + } + } + } + components.push_back(current_component); + } + } + return components; + } + + struct SpeciesSetIntersection { + const Species species; + std::size_t count; + }; + + std::expected get_intersection_info ( + const std::unordered_set& setA, + const std::unordered_set& setB + ) { + // Iterate over the smaller of the two + auto* outerSet = &setA; + auto* innerSet = &setB; + if (setA.size() > setB.size()) { + outerSet = &setB; + innerSet = &setA; + } + + std::size_t matchCount = 0; + const Species* firstMatch = nullptr; + + for (const Species& sp : *outerSet) { + if (innerSet->contains(sp)) { + if (matchCount == 0) { + firstMatch = &sp; + } + ++matchCount; + if (matchCount > 1) { + break; + } + } + } + if (!firstMatch) { + // No matches found + return std::unexpected{"Intersection is empty"}; + } + if (matchCount == 0) { + // No matches found + return std::unexpected{"No intersection found"}; + } + + // Return the first match and the count of matches + return SpeciesSetIntersection{*firstMatch, matchCount}; + + } + + bool has_distinct_reactant_and_product_species ( + const std::unordered_set& poolSpecies, + const std::unordered_set& reactants, + const std::unordered_set& products + ) { + const auto reactant_result = get_intersection_info(poolSpecies, reactants); + if (!reactant_result) { + return false; // No reactants found + } + const auto [reactantSample, reactantCount] = reactant_result.value(); + + const auto product_result = get_intersection_info(poolSpecies, products); + if (!product_result) { + return false; // No products found + } + + const auto [productSample, productCount] = product_result.value(); + + // If either side has ≥2 distinct matches, we can always pick + // one from each that differ. + if (reactantCount > 1 || productCount > 1) { + return true; + } + + // Exactly one match on each side → they must differ + return reactantSample != productSample; + } + } namespace gridfire { @@ -48,7 +159,7 @@ namespace gridfire { return m_baseEngine.getNetworkSpecies(); } - StepDerivatives MultiscalePartitioningEngineView::calculateRHSAndEnergy( + std::expected, expectations::StaleEngineError> MultiscalePartitioningEngineView::calculateRHSAndEnergy( const std::vector &Y_full, const double T9, const double rho @@ -75,10 +186,14 @@ namespace gridfire { m_cacheStats.misses(), m_cacheStats.hits() ); - throw exceptions::StaleEngineError("QSE Cache Miss while lacking context for partitioning. This should be caught by the caller and trigger a re-partition stage."); + return std::unexpected{expectations::StaleEngineError(expectations::StaleEngineErrorTypes::SYSTEM_RESIZED)}; } m_cacheStats.hit(CacheStats::operators::CalculateRHSAndEnergy); - StepDerivatives deriv = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); + const auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); + if (!result) { + return std::unexpected{result.error()}; + } + auto deriv = result.value(); for (size_t i = 0; i < m_algebraic_species_indices.size(); ++i) { const size_t species_index = m_algebraic_species_indices[i]; @@ -173,7 +288,12 @@ namespace gridfire { return m_baseEngine.getNetworkReactions(); } - std::unordered_map MultiscalePartitioningEngineView::getSpeciesTimescales( + void MultiscalePartitioningEngineView::setNetworkReactions(const reaction::LogicalReactionSet &reactions) { + LOG_CRITICAL(m_logger, "setNetworkReactions is not supported in MultiscalePartitioningEngineView. Did you mean to call this on the base engine?"); + throw exceptions::UnableToSetNetworkReactionsError("setNetworkReactions is not supported in MultiscalePartitioningEngineView. Did you mean to call this on the base engine?"); + } + + std::expected, expectations::StaleEngineError> MultiscalePartitioningEngineView::getSpeciesTimescales( const std::vector &Y, const double T9, const double rho @@ -192,22 +312,58 @@ namespace gridfire { throw exceptions::StaleEngineError("QSE Cache Miss while lacking context for partitioning. This should be caught by the caller and trigger a re-partition stage."); } m_cacheStats.hit(CacheStats::operators::GetSpeciesTimescales); - std::unordered_map speciesTimescales = m_baseEngine.getSpeciesTimescales(Y, T9, rho); + const auto result = m_baseEngine.getSpeciesTimescales(Y, T9, rho); + if (!result) { + return std::unexpected{result.error()}; + } + std::unordered_map speciesTimescales = result.value(); for (const auto& algebraicSpecies : m_algebraic_species) { speciesTimescales[algebraicSpecies] = std::numeric_limits::infinity(); // Algebraic species have infinite timescales. } return speciesTimescales; } + std::expected, expectations::StaleEngineError> + MultiscalePartitioningEngineView::getSpeciesDestructionTimescales( + const std::vector &Y, + double T9, + double rho + ) const { + const auto key = QSECacheKey(T9, rho, Y); + if (!m_qse_abundance_cache.contains(key)) { + m_cacheStats.miss(CacheStats::operators::GetSpeciesDestructionTimescales); + LOG_ERROR( + m_logger, + "QSE abundance cache miss for T9 = {}, rho = {} (misses: {}, hits: {}). getSpeciesDestructionTimescales does not receive sufficient context to partition and stabilize the network. Throwing an error which should be caught by the caller and trigger a re-partition stage.", + T9, + rho, + m_cacheStats.misses(), + m_cacheStats.hits() + ); + throw exceptions::StaleEngineError("QSE Cache Miss while lacking context for partitioning. This should be caught by the caller and trigger a re-partition stage."); + } + m_cacheStats.hit(CacheStats::operators::GetSpeciesDestructionTimescales); + const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y, T9, rho); + if (!result) { + return std::unexpected{result.error()}; + } + std::unordered_map speciesDestructionTimescales = result.value(); + for (const auto& algebraicSpecies : m_algebraic_species) { + speciesDestructionTimescales[algebraicSpecies] = std::numeric_limits::infinity(); // Algebraic species have infinite destruction timescales. + } + return speciesDestructionTimescales; + } + fourdst::composition::Composition MultiscalePartitioningEngineView::update(const NetIn &netIn) { const fourdst::composition::Composition baseUpdatedComposition = m_baseEngine.update(netIn); + double T9 = netIn.temperature / 1.0e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) - const auto key = QSECacheKey( - netIn.temperature / 1.0e9, // Convert temperature from Kelvin to T9 (T9 = T / 1e9) + const auto preKey = QSECacheKey( + T9, netIn.density, packCompositionToVector(baseUpdatedComposition, m_baseEngine) ); - if (m_qse_abundance_cache.contains(key)) { + if (m_qse_abundance_cache.contains(preKey)) { return baseUpdatedComposition; } NetIn baseUpdatedNetIn = netIn; @@ -218,7 +374,17 @@ namespace gridfire { const size_t species_index = m_algebraic_species_indices[i]; Y_algebraic[i] = equilibratedComposition.getMolarAbundance(m_baseEngine.getNetworkSpecies()[species_index]); } - m_qse_abundance_cache[key] = std::move(Y_algebraic); + + // We store the algebraic abundances in the cache for both pre- and post-conditions to avoid recalculating them. + m_qse_abundance_cache[preKey] = Y_algebraic; + + const auto postKey = QSECacheKey( + T9, + netIn.density, + packCompositionToVector(equilibratedComposition, m_baseEngine) + ); + m_qse_abundance_cache[postKey] = Y_algebraic; + return equilibratedComposition; } @@ -229,7 +395,7 @@ namespace gridfire { packCompositionToVector(netIn.composition, m_baseEngine) ); if (m_qse_abundance_cache.contains(key)) { - return false; // The cache hit indicates the engine is not stale for the given conditions. + return m_baseEngine.isStale(netIn); // The cache hit indicates the engine is not stale for the given conditions. } return true; } @@ -248,6 +414,28 @@ namespace gridfire { return m_baseEngine; } + std::vector> MultiscalePartitioningEngineView::analyzeTimescalePoolConnectivity( + const std::vector> ×cale_pools, + const std::vector &Y, + double T9, + double rho + ) const { + std::vector> final_connected_pools; + + for (const auto& pool : timescale_pools) { + if (pool.empty()) { + continue; // Skip empty pools + } + + // For each timescale pool, we need to analyze connectivity. + auto connectivity_graph = buildConnectivityGraph(pool); + auto components = findConnectedComponentsBFS(connectivity_graph, pool); + final_connected_pools.insert(final_connected_pools.end(), components.begin(), components.end()); + } + return final_connected_pools; + + } + void MultiscalePartitioningEngineView::partitionNetwork( const std::vector &Y, const double T9, @@ -286,40 +474,22 @@ namespace gridfire { candidate_pools.push_back(timescale_pools[i]); } + LOG_TRACE_L1(m_logger, "Preforming connectivity analysis on timescale pools..."); + const std::vector> connected_pools = analyzeTimescalePoolConnectivity(candidate_pools, Y, T9, rho); + LOG_TRACE_L1(m_logger, "Found {} connected pools (compared to {} timescale pools) for QSE analysis.", connected_pools.size(), timescale_pools.size()); + + // --- Step 5. Identify potential seed species for each candidate pool --- LOG_TRACE_L1(m_logger, "Identifying potential seed species for candidate pools..."); - const std::vector candidate_groups = constructCandidateGroups(candidate_pools, Y, T9, rho); + const std::vector candidate_groups = constructCandidateGroups(connected_pools, Y, T9, rho); LOG_TRACE_L1( m_logger, - "Found {} candidate QSE groups for further analysis. {}", + "Found {} candidate QSE groups for further analysis: {}", candidate_groups.size(), [&]() -> std::string { std::stringstream ss; - int count = 0; for (const auto& group : candidate_groups) { - ss << "CandidateGroup(seeds: ["; - int icount = 0; - for (const auto& seed_index : group.seed_indices) { - ss << m_baseEngine.getNetworkSpecies()[seed_index].name(); - if (icount < group.seed_indices.size() - 1) { - ss << ", "; - } - icount++; - } - ss << "], qse species: ["; - icount = 0; - for (const auto& qse_index : group.algebraic_indices) { - ss << m_baseEngine.getNetworkSpecies()[qse_index].name(); - if (icount < group.algebraic_indices.size() - 1) { - ss << ", "; - } - icount++; - } - ss << "])"; - if (count < candidate_groups.size() - 1) { - ss << ", "; - } - count++; + ss << group << " "; } return ss.str(); }() @@ -364,6 +534,27 @@ namespace gridfire { } } + LOG_INFO( + m_logger, + "Partitioning complete. Found {} dynamic species, {} algebraic (QSE) species ({}) spread over {} QSE group{}.", + m_dynamic_species.size(), + m_algebraic_species.size(), + [&]() -> std::string { + std::stringstream ss; + int count = 0; + for (const auto& species : m_algebraic_species) { + ss << species.name(); + if (m_algebraic_species.size() > 1 && count < m_algebraic_species.size() - 1) { + ss << ", "; + } + count++; + } + return ss.str(); + }(), + m_qse_groups.size(), + m_qse_groups.size() == 1 ? "" : "s" + ); + } void MultiscalePartitioningEngineView::partitionNetwork( @@ -674,7 +865,13 @@ namespace gridfire { const double rho ) const { LOG_TRACE_L1(m_logger, "Partitioning by timescale..."); - const auto all_timescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + const auto result= m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho); + if (!result) { + LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state"); + m_logger->flush_log(); + throw exceptions::StaleEngineError("Failed to get species timescales due to stale engine state"); + } + std::unordered_map all_timescales = result.value(); const auto& all_species = m_baseEngine.getNetworkSpecies(); std::vector> sorted_timescales; @@ -699,7 +896,7 @@ namespace gridfire { } constexpr double ABSOLUTE_QSE_TIMESCALE_THRESHOLD = 3.156e7; // Absolute threshold for QSE timescale (1 yr) - constexpr double MIN_GAP_THRESHOLD = 3.0; // Require a 3 order of magnitude gap + constexpr double MIN_GAP_THRESHOLD = 2.0; // Require a 2 order of magnitude gap LOG_TRACE_L1(m_logger, "Found {} species with finite timescales.", sorted_timescales.size()); LOG_TRACE_L1(m_logger, "Absolute QSE timescale threshold: {} seconds ({} years).", @@ -712,11 +909,11 @@ namespace gridfire { // 1. First Pass: Absolute Timescale Cutoff for (const auto& ts_pair : sorted_timescales) { if (ts_pair.first > ABSOLUTE_QSE_TIMESCALE_THRESHOLD) { - LOG_TRACE_L1(m_logger, "Species {} with timescale {} is considered dynamic (slower than qse timescale threshold).", + LOG_TRACE_L3(m_logger, "Species {} with timescale {} is considered dynamic (slower than qse timescale threshold).", all_species[ts_pair.second].name(), ts_pair.first); dynamic_pool_indices.push_back(ts_pair.second); } else { - LOG_TRACE_L1(m_logger, "Species {} with timescale {} is a candidate fast species (faster than qse timescale threshold).", + LOG_TRACE_L3(m_logger, "Species {} with timescale {} is a candidate fast species (faster than qse timescale threshold).", all_species[ts_pair.second].name(), ts_pair.first); fast_candidates.push_back(ts_pair); } @@ -738,7 +935,7 @@ namespace gridfire { const double t1 = fast_candidates[i].first; const double t2 = fast_candidates[i+1].first; if (std::log10(t1) - std::log10(t2) > MIN_GAP_THRESHOLD) { - LOG_TRACE_L1(m_logger, "Detected gap between species {} (timescale {:0.2E}) and {} (timescale {:0.2E}).", + LOG_TRACE_L3(m_logger, "Detected gap between species {} (timescale {:0.2E}) and {} (timescale {:0.2E}).", all_species[fast_candidates[i].second].name(), t1, all_species[fast_candidates[i+1].second].name(), t2); split_points.push_back(i + 1); @@ -787,140 +984,32 @@ namespace gridfire { } - std::unordered_map> MultiscalePartitioningEngineView::buildConnectivityGraph( - const std::unordered_set &fast_reaction_indices - ) const { - const auto& all_reactions = m_baseEngine.getNetworkReactions(); - std::unordered_map> connectivity; - for (const size_t reaction_idx : fast_reaction_indices) { - const auto& reaction = all_reactions[reaction_idx]; - const auto& reactants = reaction.reactants(); - const auto& products = reaction.products(); - - // For each fast reaction, create edges between all reactants and all products. - // This represents that nucleons can flow quickly between these species. - for (const auto& reactant : reactants) { - const size_t reactant_idx = m_baseEngine.getSpeciesIndex(reactant); - for (const auto& product : products) { - const size_t product_idx = m_baseEngine.getSpeciesIndex(product); - - // Add a two-way edge to the adjacency list. - connectivity[reactant_idx].push_back(product_idx); - connectivity[product_idx].push_back(reactant_idx); - } - } - } - return connectivity; - } - - // std::vector MultiscalePartitioningEngineView::refineCandidateGroup( - // const std::vector &candidate_group, - // const std::vector &Y, - // const double T9, - // const double rho + // std::unordered_map> MultiscalePartitioningEngineView::buildConnectivityGraph( + // const std::unordered_set &fast_reaction_indices // ) const { - // LOG_TRACE_L1(m_logger, "Refining candidate group with {} species.", candidate_group.size()); + // const auto& all_reactions = m_baseEngine.getNetworkReactions(); + // std::unordered_map> connectivity; + // for (const size_t reaction_idx : fast_reaction_indices) { + // const auto& reaction = all_reactions[reaction_idx]; + // const auto& reactants = reaction.reactants(); + // const auto& products = reaction.products(); // - // // --- Step 1. Graph Construction --- - // std::unordered_map>> weighted_graph; - // std::vector> edge_list; - // std::unordered_map node_strength; + // // For each fast reaction, create edges between all reactants and all products. + // // This represents that nucleons can flow quickly between these species. + // for (const auto& reactant : reactants) { + // const size_t reactant_idx = m_baseEngine.getSpeciesIndex(reactant); + // for (const auto& product : products) { + // const size_t product_idx = m_baseEngine.getSpeciesIndex(product); // - // for (const auto& node_idx : candidate_group) { - // node_strength[node_idx] = 0.0; - // } - // - // std::unordered_set candidate_set(candidate_group.begin(), candidate_group.end()); - // - // for (const auto& reaction : m_baseEngine.getNetworkReactions()) { - // bool is_internal = true; - // std::vector involved_indices; - // - // // Check if all reactants and products are in the candidate group. - // std::vector participating_species = reaction.reactants(); - // participating_species.insert(participating_species.end(), reaction.products().begin(), reaction.products().end()); - // - // for (const auto& sp : participating_species) { - // size_t idx = m_baseEngine.getSpeciesIndex(sp); - // if (! candidate_set.contains(idx)) { - // is_internal = false; - // break; - // } - // } - // if (!is_internal) { - // continue; // Skip reactions that are not fully internal to the candidate group. - // } - // - // const double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho)); - // if (flow < 1e-99) continue; - // - // for (const auto& reactant : reaction.reactants()) { - // for (const auto& product : reaction.products()) { - // size_t u = m_baseEngine.getSpeciesIndex(reactant); - // size_t v = m_baseEngine.getSpeciesIndex(product); - // - // if (u == v) { - // continue; // Skip self-loops - // } - // - // weighted_graph[u].push_back(std::make_tuple(v, flow)); - // weighted_graph[v].push_back(std::make_tuple(u, flow)); - // - // node_strength[u] += flow; - // node_strength[v] += flow; - // - // if (u < v) { - // edge_list.emplace_back(u, v, flow); - // } + // // Add a two-way edge to the adjacency list. + // connectivity[reactant_idx].push_back(product_idx); + // connectivity[product_idx].push_back(reactant_idx); // } // } // } - // - // // --- Phase 2. Prune weak bridges --- - // auto pruned_graph = weighted_graph; - // constexpr double beta = 0.01; // Weak link threshold: flux < 1% of the node strength - // - // for (const auto& edge : edge_list) { - // size_t u = std::get<0>(edge); - // size_t v = std::get<1>(edge); - // double flow = std::get<2>(edge); - // - // if (flow < (beta * node_strength[u]) && flow < (beta * node_strength[v])) { - // LOG_TRACE_L1(m_logger, "Pruning weak link ({}, {}) with flow {} below threshold {}", - // u, v, flow, beta * std::min(node_strength[u], node_strength[v])); - // - // auto& u_adj = pruned_graph.at(u); - // std::erase_if( - // u_adj, - // [v](const auto& p) { - // return p.first == v; - // } - // ); - // auto& v_adj = pruned_graph.at(v); - // std::erase_if( - // v_adj, - // [u](const auto& p) { - // return p.first == u; - // } - // ); - // } - // } - // - // // --- Phase 3. Find final connected components --- - // auto final_components = findConnectedComponentsBFS(pruned_graph, candidate_group); - // - // std::vector refined_groups; - // for (const auto& component : final_components) { - // if (!component.empty()) { - // std::set tempSet(component.begin(), component.end()); - // refined_groups.push_back({tempSet, false}); // Initialize with is_in_equilibrium = false - // } - // } - // - // return refined_groups; + // return connectivity; // } - std::vector MultiscalePartitioningEngineView::validateGroupsWithFluxAnalysis( const std::vector &candidate_groups, @@ -967,15 +1056,89 @@ namespace gridfire { // Classify the reaction based on its participants if ((has_internal_reactant && has_internal_product) && !(has_external_reactant || has_external_product)) { + LOG_TRACE_L3( + m_logger, + "Reaction {} is internal to the group containing {} and contributes to internal flux by {}", + reaction.id(), + [&]() -> std::string { + std::stringstream ss; + int count = 0; + for (const auto& idx : group.algebraic_indices) { + ss << m_baseEngine.getNetworkSpecies()[idx].name(); + if (count < group.species_indices.size() - 1) { + ss << ", "; + } + count++; + } + return ss.str(); + }(), + flow + ); internal_flux += flow; // Internal flux within the group } else if ((has_internal_reactant || has_internal_product) && (has_external_reactant || has_external_product)) { + LOG_TRACE_L3( + m_logger, + "Reaction {} is external to the group containing {} and contributes to external flux by {}", + reaction.id(), + [&]() -> std::string { + std::stringstream ss; + int count = 0; + for (const auto& idx : group.algebraic_indices) { + ss << m_baseEngine.getNetworkSpecies()[idx].name(); + if (count < group.species_indices.size() - 1) { + ss << ", "; + } + count++; + } + return ss.str(); + }(), + flow + ); external_flux += flow; // External flux to/from the group } // otherwise the reaction is fully decoupled from the QSE group and can be ignored. } if (internal_flux > FLUX_RATIO_THRESHOLD * external_flux) { + LOG_TRACE_L1( + m_logger, + "Group containing {} is in equilibrium: internal flux = {}, external flux = {}, ratio = {}", + [&]() -> std::string { + std::stringstream ss; + int count = 0; + for (const auto& idx : group.algebraic_indices) { + ss << m_baseEngine.getNetworkSpecies()[idx].name(); + if (count < group.species_indices.size() - 1) { + ss << ", "; + } + count++; + } + return ss.str(); + }(), + internal_flux, + external_flux, + internal_flux / external_flux + ); group.is_in_equilibrium = true; // This group is in equilibrium if internal flux is significantly larger than external flux. } else { + LOG_TRACE_L1( + m_logger, + "Group containing {} is NOT in equilibrium: internal flux = {}, external flux = {}, ratio = {}", + [&]() -> std::string { + std::stringstream ss; + int count = 0; + for (const auto& idx : group.algebraic_indices) { + ss << m_baseEngine.getNetworkSpecies()[idx].name(); + if (count < group.species_indices.size() - 1) { + ss << ", "; + } + count++; + } + return ss.str(); + }(), + internal_flux, + external_flux, + internal_flux / external_flux + ); group.is_in_equilibrium = false; } } @@ -989,16 +1152,45 @@ namespace gridfire { ) { LOG_TRACE_L1(m_logger, "Solving for QSE abundances..."); auto Y_out = Y_full; - auto species_timescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); - for (const auto& qse_group : m_qse_groups) { - if (!qse_group.is_in_equilibrium || qse_group.species_indices.empty()) { + + // Sort by timescale to solve fastest QSE groups first (these can seed slower groups) + std::ranges::sort(m_qse_groups, [](const QSEGroup& a, const QSEGroup& b) { + return a.mean_timescale < b.mean_timescale; + }); + + for (const auto&[species_indices, is_in_equilibrium, algebraic_indices, seed_indices, mean_timescale] : m_qse_groups) { + if (!is_in_equilibrium || species_indices.empty()) { LOG_TRACE_L1( m_logger, - "Skipping QSE group with {} species ({} algebraic, {} seeds) as it is not in equilibrium.", - qse_group.species_indices.size(), - qse_group.algebraic_indices.size(), - qse_group.seed_indices.size() + "Skipping QSE group with {} species ({} algebraic ({}), {} seeds ({})) as it is not in equilibrium.", + species_indices.size(), + algebraic_indices.size(), + [&]() -> std::string { + std::ostringstream os; + int count = 0; + for (const auto& idx : algebraic_indices) { + os << m_baseEngine.getNetworkSpecies()[idx].name(); + if (count < algebraic_indices.size() - 1) { + os << ", "; + } + count++; + } + return os.str(); + }(), + seed_indices.size(), + [&]() -> std::string { + std::ostringstream os; + int count = 0; + for (const auto& idx : seed_indices) { + os << m_baseEngine.getNetworkSpecies()[idx].name(); + if (count < seed_indices.size() - 1) { + os << ", "; + } + count++; + } + return os.str(); + }() ); continue; } @@ -1006,13 +1198,13 @@ namespace gridfire { LOG_TRACE_L1( m_logger, "Solving for QSE abundances for group with {} species ([{}] algebraic, [{}] seeds).", - qse_group.species_indices.size(), + species_indices.size(), [&]() -> std::string { std::stringstream ss; int count = 0; - for (const auto& idx : qse_group.algebraic_indices) { + for (const auto& idx : algebraic_indices) { ss << m_baseEngine.getNetworkSpecies()[idx].name(); - if (count < qse_group.algebraic_indices.size() - 1) { + if (count < algebraic_indices.size() - 1) { ss << ", "; } count++; @@ -1022,9 +1214,9 @@ namespace gridfire { [&]() -> std::string { std::stringstream ss; int count = 0; - for (const auto& idx : qse_group.seed_indices) { + for (const auto& idx : seed_indices) { ss << m_baseEngine.getNetworkSpecies()[idx].name(); - if (count < qse_group.seed_indices.size() - 1) { + if (count < seed_indices.size() - 1) { ss << ", "; } count++; @@ -1034,16 +1226,16 @@ namespace gridfire { ); std::vector qse_solve_indices; - std::vector seed_indices; + std::vector seed_indices_vec; - seed_indices.reserve(qse_group.species_indices.size()); - qse_solve_indices.reserve(qse_group.species_indices.size()); + seed_indices_vec.reserve(species_indices.size()); + qse_solve_indices.reserve(species_indices.size()); - for (size_t seed_idx : qse_group.seed_indices) { - seed_indices.emplace_back(seed_idx); + for (size_t seed_idx : seed_indices) { + seed_indices_vec.emplace_back(seed_idx); } - for (size_t algebraic_idx : qse_group.algebraic_indices) { + for (size_t algebraic_idx : algebraic_indices) { qse_solve_indices.emplace_back(algebraic_idx); } @@ -1070,7 +1262,7 @@ namespace gridfire { std::stringstream msg; msg << "QSE solver failed with status: " << status << " for group with seed "; if (seed_indices.size() == 1) { - msg << "nucleus " << m_baseEngine.getNetworkSpecies()[seed_indices[0]].name(); + msg << "nucleus " << m_baseEngine.getNetworkSpecies()[seed_indices_vec[0]].name(); } else { msg << "nuclei "; // TODO: Refactor nice list printing into utility function somewhere @@ -1114,17 +1306,45 @@ namespace gridfire { const double T9, const double rho ) const { - const auto& all_timescales = m_baseEngine.getSpeciesTimescales(Y, T9, rho); + const auto& result = m_baseEngine.getSpeciesDestructionTimescales(Y, T9, rho); + if (!result) { + LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state"); + m_logger->flush_log(); + throw exceptions::StaleEngineError("Failed to get species timescales due to stale engine state"); + } + const std::unordered_map all_timescales = result.value(); + + const auto& all_species = m_baseEngine.getNetworkSpecies(); + size_t slowest_pool_index = 0; // Default to the first pool if no valid pool is found double slowest_mean_timescale = std::numeric_limits::min(); + int count = 0; for (const auto& pool : pools) { double mean_timescale = 0.0; for (const auto& species_idx : pool) { - double timescale = all_timescales.at(all_species[species_idx]); + const double timescale = all_timescales.at(all_species[species_idx]); mean_timescale += timescale; } mean_timescale = mean_timescale / pool.size(); + if (std::isinf(mean_timescale)) { + LOG_CRITICAL(m_logger, "Encountered infinite mean timescale for pool {} with species: {}", + count, [&]() -> std::string { + std::stringstream ss; + int iCount = 0; + for (const auto& idx : pool) { + ss << all_species[idx].name() << ": " << all_timescales.at(all_species[idx]); + if (iCount < pool.size() - 1) { + ss << ", "; + } + iCount++; + } + return ss.str(); + }() + ); + m_logger->flush_log(); + throw std::logic_error("Encountered infinite mean destruction timescale for a pool while identifying the mean slowest pool set, indicating a potential issue with species timescales. Check log file for more details on specific pool composition..."); + } if (mean_timescale > slowest_mean_timescale) { slowest_mean_timescale = mean_timescale; slowest_pool_index = &pool - &pools[0]; // Get the index of the pool @@ -1133,18 +1353,88 @@ namespace gridfire { return slowest_pool_index; } + std::unordered_map> MultiscalePartitioningEngineView::buildConnectivityGraph( + const std::vector &species_pool + ) const { + std::unordered_map> connectivity_graph; + const std::set pool_set(species_pool.begin(), species_pool.end()); + const std::unordered_set pool_species = [&]() -> std::unordered_set { + std::unordered_set result; + for (const auto& species_idx : species_pool) { + Species species = m_baseEngine.getNetworkSpecies()[species_idx]; + result.insert(species); + } + return result; + }(); + + std::map> speciesReactionMap; + std::vector candidate_reactions; + + auto getSpeciesIdx = [&](const std::vector &species) -> std::vector { + std::vector result; + result.reserve(species.size()); + for (const auto& s : species) { + size_t idx = m_baseEngine.getSpeciesIndex(s); + result.push_back(idx); + } + return result; + }; + + for (const auto& reaction : m_baseEngine.getNetworkReactions()) { + const std::vector &reactants = reaction.reactants(); + const std::vector &products = reaction.products(); + + std::unordered_set reactant_set(reactants.begin(), reactants.end()); + std::unordered_set product_set(products.begin(), products.end()); + + // Only consider reactions where at least one distinct reactant and product are in the pool + if (has_distinct_reactant_and_product_species(pool_species, reactant_set, product_set)) { + std::vector involvedIDs = getSpeciesIdx(reactants); + std::vector involvedProducts = getSpeciesIdx(products); + involvedIDs.insert(involvedIDs.end(), involvedProducts.begin(), involvedProducts.end()); + std::set involvedSet(involvedIDs.begin(), involvedIDs.end()); + + std::vector intersection; + intersection.reserve(involvedSet.size()); + + std::ranges::set_intersection(pool_set, involvedSet, std::back_inserter(intersection)); + + // Add clique + for (const size_t& u : intersection) { + const auto& uSpecies = m_baseEngine.getNetworkSpecies()[u]; + for (const size_t& v : intersection) { + const auto & vSpecies = m_baseEngine.getNetworkSpecies()[v]; + if (u != v) { // Avoid self-loops + connectivity_graph[u].push_back(v); + } + } + } + } + } + + + return connectivity_graph; + } + std::vector MultiscalePartitioningEngineView::constructCandidateGroups( - const std::vector> ×cale_pools, + const std::vector> &candidate_pools, const std::vector &Y, - const double T9, - const double rho + const double T9, const double rho ) const { const auto& all_species = m_baseEngine.getNetworkSpecies(); const auto& all_reactions = m_baseEngine.getNetworkReactions(); + const auto& result = m_baseEngine.getSpeciesDestructionTimescales(Y, T9, rho); + if (!result) { + LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state"); + m_logger->flush_log(); + throw exceptions::StaleEngineError("Failed to get species timescales due to stale engine state"); + } + const std::unordered_map destruction_timescales = result.value(); std::vector candidate_groups; - for (const auto& pool : timescale_pools) { + for (const auto& pool : candidate_pools) { if (pool.empty()) continue; // Skip empty pools + // For each pool first identify all topological bridge connections std::vector> bridge_reactions; for (const auto& species_idx : pool) { @@ -1157,15 +1447,13 @@ namespace gridfire { for (const auto& reactant : reaction.reactants()) { if (std::ranges::find(pool, m_baseEngine.getSpeciesIndex(reactant)) == pool.end()) { has_external_reactant = true; - LOG_TRACE_L2(m_logger, "Found external reactant {} in reaction {} for species {}.", - reactant.name(), reaction.id(), ash.name()); + LOG_TRACE_L3(m_logger, "Found external reactant {} in reaction {} for species {}.", reactant.name(), reaction.id(), ash.name()); break; // Found an external reactant, no need to check further } } if (has_external_reactant) { double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho)); - LOG_TRACE_L2(m_logger, "Found bridge reaction {} with flow {} for species {}.", - reaction.id(), flow, ash.name()); + LOG_TRACE_L3(m_logger, "Found bridge reaction {} with flow {} for species {}.", reaction.id(), flow, ash.name()); bridge_reactions.push_back({reaction, flow}); } } @@ -1176,23 +1464,6 @@ namespace gridfire { [](const auto& a, const auto& b) { return a.second > b.second; // Sort by flow in descending order }); - LOG_TRACE_L1( - m_logger, - "Found a total of {} bridge reactions for the pool with qse species: {}.", - bridge_reactions.size(), - [&]() -> std::string { - std::stringstream ss; - int count = 0; - for (const auto& idx : pool) { - ss << all_species[idx].name(); - if (count < pool.size() - 1) { - ss << ", "; - } - count++; - } - return ss.str(); - }() - ); constexpr double MIN_GAP_THRESHOLD = 1; // Minimum logarithmic molar flow gap threshold for bridge reactions std::vector split_points; @@ -1200,38 +1471,10 @@ namespace gridfire { const double f1 = bridge_reactions[i].second; const double f2 = bridge_reactions[i + 1].second; if (std::log10(f1) - std::log10(f2) > MIN_GAP_THRESHOLD) { - LOG_TRACE_L2(m_logger, "Detected gap between bridge reactions with flows {} and {}.", - f1, f2); + LOG_TRACE_L3(m_logger, "Detected gap between bridge reactions with flows {} and {}.", f1, f2); split_points.push_back(i + 1); } } - LOG_TRACE_L1(m_logger, "Found {} candidate molar flow split points in bridge reactions.", split_points.size() == 0 ? 1 : split_points.size()); - LOG_TRACE_L1(m_logger, "Split point summary: {}", - [&]() -> std::string { - std::stringstream ss; - int count = 0; - // Print out which reactions are in each molar reaction flow pool - int lastSplitPoint = 0; - for (const auto& split_point : split_points) { - ss << "Reactions in molar flow pool " << count + 1 << ": ["; - int icount = 0; - for (size_t i = lastSplitPoint; i < split_point; ++i) { - ss << bridge_reactions[i].first.id(); - if (icount < (split_point - lastSplitPoint) - 1) { - ss << ", "; - } - icount++; - } - ss << "]"; - if (count < split_points.size() - 1) { - ss << ", "; - } - count++; - lastSplitPoint = split_point; - } - return ss.str(); - }() - ); if (split_points.empty()) { // If no split points were found, we consider the whole set of bridge reactions as one group. split_points.push_back(bridge_reactions.size() - 1); @@ -1253,7 +1496,18 @@ namespace gridfire { } const std::set poolSet(pool.begin(), pool.end()); const std::set seedSet(seed_indices.begin(), seed_indices.end()); - QSEGroup qse_group(all_indices, false, poolSet, seedSet); + + double mean_timescale = 0.0; + for (const auto& pool_idx : poolSet) { + const auto& species = all_species[pool_idx]; + if (destruction_timescales.contains(species)) { + mean_timescale += std::min(destruction_timescales.at(species), species.halfLife()); // Use the minimum of destruction timescale and half-life + } else { + mean_timescale += species.halfLife(); + } + } + mean_timescale /= poolSet.size(); + QSEGroup qse_group(all_indices, false, poolSet, seedSet, mean_timescale); candidate_groups.push_back(qse_group); } return candidate_groups; @@ -1269,7 +1523,11 @@ namespace gridfire { y_trial[m_qse_solve_indices[i]] = y_qse(i); } - const auto [dydt, energy] = m_view->getBaseEngine().calculateRHSAndEnergy(y_trial, m_T9, m_rho); + const auto result = m_view->getBaseEngine().calculateRHSAndEnergy(y_trial, m_T9, m_rho); + if (!result) { + throw exceptions::StaleEngineError("Failed to calculate RHS and energy due to stale engine state"); + } + const auto&[dydt, nuclearEnergyGenerationRate] = result.value(); f_qse.resize(m_qse_solve_indices.size()); for (size_t i = 0; i < m_qse_solve_indices.size(); ++i) { f_qse(i) = dydt[m_qse_solve_indices[i]]; @@ -1348,4 +1606,21 @@ namespace gridfire { bool QSECacheKey::operator==(const QSECacheKey &other) const { return m_hash == other.m_hash; } + + bool MultiscalePartitioningEngineView::QSEGroup::operator==(const QSEGroup &other) const { + return mean_timescale == other.mean_timescale; + } + + bool MultiscalePartitioningEngineView::QSEGroup::operator<(const QSEGroup &other) const { + return mean_timescale < other.mean_timescale; + } + + bool MultiscalePartitioningEngineView::QSEGroup::operator>(const QSEGroup &other) const { + return mean_timescale > other.mean_timescale; + } + + bool MultiscalePartitioningEngineView::QSEGroup::operator!=(const QSEGroup &other) const { + return !(*this == other); + } + } diff --git a/src/network/lib/network.cpp b/src/network/lib/network.cpp index 850a3043..145e5a62 100644 --- a/src/network/lib/network.cpp +++ b/src/network/lib/network.cpp @@ -64,7 +64,7 @@ namespace gridfire { // Trim whitespace from both ends of a string std::string trim_whitespace(const std::string& str) { auto startIt = str.begin(); - auto endIt = str.end(); + const auto endIt = str.end(); while (startIt != endIt && std::isspace(static_cast(*startIt))) { ++startIt; @@ -72,8 +72,8 @@ namespace gridfire { if (startIt == endIt) { return ""; } - auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt), - [](unsigned char ch){ return !std::isspace(ch); }); + const auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt), + [](const unsigned char ch){ return !std::isspace(ch); }); return std::string(startIt, ritr.base()); } diff --git a/src/network/lib/partition/partition_rauscher_thielemann.cpp b/src/network/lib/partition/partition_rauscher_thielemann.cpp index 6a825241..e3a79a51 100644 --- a/src/network/lib/partition/partition_rauscher_thielemann.cpp +++ b/src/network/lib/partition/partition_rauscher_thielemann.cpp @@ -22,17 +22,17 @@ namespace gridfire::partition { m_partitionData.reserve(numRecords); const auto* records = reinterpret_cast(rauscher_thielemann_partition_data); for (size_t i = 0; i < numRecords; ++i) { - const auto& record = records[i]; + const auto&[z, a, ground_state_spin, normalized_g_values] = records[i]; IsotopeData data; - data.ground_state_spin = record.ground_state_spin; - std::ranges::copy(record.normalized_g_values, data.normalized_g_values.begin()); - const int key = make_key(record.z, record.a); + data.ground_state_spin = ground_state_spin; + std::ranges::copy(normalized_g_values, data.normalized_g_values.begin()); + const int key = make_key(z, a); LOG_TRACE_L3_LIMIT_EVERY_N( 100, m_logger, "(EVERY 100) Adding Rauscher-Thielemann partition data for Z={} A={} (key={})", - record.z, - record.a, + z, + a, key ); diff --git a/src/network/lib/reaction/reaclib.cpp b/src/network/lib/reaction/reaclib.cpp index 87e54368..0f84eeac 100644 --- a/src/network/lib/reaction/reaclib.cpp +++ b/src/network/lib/reaction/reaclib.cpp @@ -12,7 +12,7 @@ std::string trim_whitespace(const std::string& str) { auto startIt = str.begin(); - auto endIt = str.end(); + const auto endIt = str.end(); while (startIt != endIt && std::isspace(static_cast(*startIt))) { ++startIt; @@ -20,8 +20,8 @@ std::string trim_whitespace(const std::string& str) { if (startIt == endIt) { return ""; } - auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt), - [](unsigned char ch){ return !std::isspace(ch); }); + const auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt), + [](const unsigned char ch){ return !std::isspace(ch); }); return std::string(startIt, ritr.base()); } @@ -83,41 +83,41 @@ namespace gridfire::reaclib { // 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); + constexpr 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]; + const auto&[chapter, qValue, coeffs, reverse, label, rpName, reactants_str, products_str] = 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))); + const std::string_view label_sv(label, strnlen(label, sizeof(label))); + const std::string_view rpName_sv(rpName, strnlen(rpName, sizeof(rpName))); + const std::string_view reactants_sv(reactants_str, strnlen(reactants_str, sizeof(reactants_str))); + const std::string_view products_sv(products_str, strnlen(products_str, sizeof(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] + coeffs[0], coeffs[1], coeffs[2], + coeffs[3], coeffs[4], coeffs[5], + 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, + chapter, reactants, products, - record.qValue, + qValue, label_sv, rate_coeffs, - record.reverse + reverse ); } diff --git a/src/network/lib/reaction/reaction.cpp b/src/network/lib/reaction/reaction.cpp index 1387eff6..8917a293 100644 --- a/src/network/lib/reaction/reaction.cpp +++ b/src/network/lib/reaction/reaction.cpp @@ -154,7 +154,7 @@ namespace gridfire::reaction { return (reactantMass - productMass) * AMU2MeV; } - uint64_t Reaction::hash(uint64_t seed) const { + uint64_t Reaction::hash(const uint64_t seed) const { return XXHash64::hash(m_id.data(), m_id.size(), seed); } diff --git a/src/network/lib/solver/solver.cpp b/src/network/lib/solver/solver.cpp index a41077cf..57df8e22 100644 --- a/src/network/lib/solver/solver.cpp +++ b/src/network/lib/solver/solver.cpp @@ -3,518 +3,22 @@ #include "gridfire/network.h" #include "gridfire/exceptions/error_engine.h" - -#include "gridfire/utils/logging.h" - #include "fourdst/composition/atomicSpecies.h" #include "fourdst/composition/composition.h" #include "fourdst/config/config.h" -#include "Eigen/Dense" #include "unsupported/Eigen/NonLinearOptimization" #include #include -#include #include #include #include -#include "gridfire/engine/views/engine_multiscale.h" - #include "quill/LogMacros.h" -struct adaptive_integrator_observer { - double m_last_t = 0.0; - int m_total_steps = 0; - - void operator()(boost::numeric::ublas::vector &state, double t) { - // std::cout << "(Step: " << m_total_steps << ") Current Time: " << t << " (dt: " << t - m_last_t << ")\n"; - m_last_t = t; - m_total_steps++; - - } -}; - namespace gridfire::solver { - double s_prevTimestep = 0.0; - - 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; - - NetOut postIgnition = initializeNetworkWithShortIgnition(netIn); - - constexpr double abundance_floor = 1.0e-30; - std::vectorY_sanitized_initial; - Y_sanitized_initial.reserve(m_engine.getNetworkSpecies().size()); - - LOG_DEBUG(m_logger, "Sanitizing initial abundances with a floor of {:0.3E}...", abundance_floor); - for (const auto& species : m_engine.getNetworkSpecies()) { - double molar_abundance = 0.0; - if (postIgnition.composition.contains(species)) { - molar_abundance = postIgnition.composition.getMolarAbundance(std::string(species.name())); - } - - if (molar_abundance < abundance_floor) { - molar_abundance = abundance_floor; - } - Y_sanitized_initial.push_back(molar_abundance); - } - 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 - - const auto indices = packSpeciesTypeIndexVectors(Y_sanitized_initial, T9, rho); - Eigen::VectorXd Y_QSE; - try { - Y_QSE = calculateSteadyStateAbundances(Y_sanitized_initial, T9, rho, indices); - LOG_DEBUG(m_logger, "QSE Abundances: {}", [*this](const dynamicQSESpeciesIndices& indices, const Eigen::VectorXd& Y_QSE) -> std::string { - std::stringstream ss; - ss << std::scientific << std::setprecision(5); - for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) { - ss << std::string(m_engine.getNetworkSpecies()[indices.QSESpeciesIndices[i]].name()) + ": "; - ss << Y_QSE(i); - if (i < indices.QSESpeciesIndices.size() - 2) { - ss << ", "; - } else if (i == indices.QSESpeciesIndices.size() - 2) { - ss << ", and "; - } - - } - return ss.str(); - }(indices, Y_QSE)); - } catch (const std::runtime_error& e) { - LOG_ERROR(m_logger, "Failed to calculate steady state abundances. Aborting QSE evaluation."); - m_logger->flush_log(); - throw std::runtime_error("Failed to calculate steady state abundances: " + std::string(e.what())); - } - - state_type YDynamic_ublas(indices.dynamicSpeciesIndices.size() + 1); - for (size_t i = 0; i < indices.dynamicSpeciesIndices.size(); ++i) { - YDynamic_ublas(i) = Y_sanitized_initial[indices.dynamicSpeciesIndices[i]]; - } - YDynamic_ublas(indices.dynamicSpeciesIndices.size()) = 0.0; // Placeholder for specific energy rate - - const RHSFunctor rhs_functor(m_engine, indices.dynamicSpeciesIndices, indices.QSESpeciesIndices, Y_QSE, T9, rho); - const auto stepper = make_controlled>(1.0e-8, 1.0e-8); - - size_t stepCount = integrate_adaptive( - stepper, - rhs_functor, - YDynamic_ublas, - 0.0, // Start time - netIn.tMax, - netIn.dt0 - ); - - std::vector YFinal = Y_sanitized_initial; - for (size_t i = 0; i speciesNames(m_engine.getNetworkSpecies().size()); - std::vector finalMassFractions(m_engine.getNetworkSpecies().size()); - double massFractionSum = 0.0; - - for (size_t i = 0; i < speciesNames.size(); ++i) { - const auto& species = m_engine.getNetworkSpecies()[i]; - speciesNames[i] = species.name(); - finalMassFractions[i] = YFinal[i] * species.mass(); // Convert from molar abundance to mass fraction - massFractionSum += finalMassFractions[i]; - } - for (auto& mf : finalMassFractions) { - mf /= massFractionSum; // Normalize to get mass fractions - } - - fourdst::composition::Composition outputComposition(speciesNames, finalMassFractions); - NetOut netOut; - netOut.composition = outputComposition; - netOut.energy = finalSpecificEnergyRate; // Specific energy rate - netOut.num_steps = stepCount; - return netOut; - } - - dynamicQSESpeciesIndices QSENetworkSolver::packSpeciesTypeIndexVectors( - const std::vector& Y, - const double T9, - const double rho - ) const { - constexpr double timescaleCutoff = 1.0e-5; - constexpr double abundanceCutoff = 1.0e-15; - - LOG_INFO(m_logger, "Partitioning species using T9={:0.2f} and ρ={:0.2e}", T9, rho); - LOG_INFO(m_logger, "Timescale Cutoff: {:.1e} s, Abundance Cutoff: {:.1e}", timescaleCutoff, abundanceCutoff); - - std::vectordynamicSpeciesIndices; // Slow species that are not in QSE - std::vectorQSESpeciesIndices; // Fast species that are in QSE - - std::unordered_map speciesTimescale = m_engine.getSpeciesTimescales(Y, T9, rho); - LOG_DEBUG( - m_logger, - "{}", - gridfire::utils::formatNuclearTimescaleLogString( - m_engine, - Y, - T9, - rho - )); - - for (size_t i = 0; i < m_engine.getNetworkSpecies().size(); ++i) { - const auto& species = m_engine.getNetworkSpecies()[i]; - const double network_timescale = speciesTimescale.at(species); - const double abundance = Y[i]; - - 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); - - const bool isQSE = false; - // network_timescale, - // decay_timescale, - // abundance - // ); - - if (isQSE) { - LOG_TRACE_L2(m_logger, "{} is in QSE based on rules in qse_rules.h", species.name()); - QSESpeciesIndices.push_back(i); - } else { - LOG_TRACE_L2(m_logger, "{} is dynamic based on rules in qse_rules.h", species.name()); - dynamicSpeciesIndices.push_back(i); - } - } - LOG_INFO(m_logger, "Partitioning complete. Dynamical species: {}, QSE species: {}.", dynamicSpeciesIndices.size(), QSESpeciesIndices.size()); - LOG_INFO(m_logger, "Dynamic species: {}", [dynamicSpeciesIndices](const DynamicEngine& engine_wrapper) -> std::string { - std::string result; - int count = 0; - for (const auto& i : dynamicSpeciesIndices) { - result += std::string(engine_wrapper.getNetworkSpecies()[i].name()); - if (count < dynamicSpeciesIndices.size() - 2) { - result += ", "; - } else if (count == dynamicSpeciesIndices.size() - 2) { - result += " and "; - } - count++; - } - return result; - }(m_engine)); - LOG_INFO(m_logger, "QSE species: {}", [QSESpeciesIndices](const DynamicEngine& engine_wrapper) -> std::string { - std::string result; - int count = 0; - for (const auto& i : QSESpeciesIndices) { - result += std::string(engine_wrapper.getNetworkSpecies()[i].name()); - if (count < QSESpeciesIndices.size() - 2) { - result += ", "; - } else if (count == QSESpeciesIndices.size() - 2) { - result += " and "; - } - count++; - } - return result; - }(m_engine)); - return {dynamicSpeciesIndices, QSESpeciesIndices}; - } - - Eigen::VectorXd QSENetworkSolver::calculateSteadyStateAbundances( - const std::vector &Y, - const double T9, - const double rho, - const dynamicQSESpeciesIndices &indices - ) const { - LOG_TRACE_L1(m_logger, "Calculating steady state abundances for QSE species..."); - LOG_TRACE_L1( - m_logger, - "Initial QSE species abundances: {}", - [&]() -> std::string { - std::stringstream ss; - ss << std::scientific << std::setprecision(5); - int count = 0; - for (const auto& i : indices.QSESpeciesIndices) { - ss << std::string(m_engine.getNetworkSpecies()[i].name()) << ": " << Y[i] << ", "; - if (count < indices.QSESpeciesIndices.size() - 2) { - ss << ", "; - } else if (count == indices.QSESpeciesIndices.size() - 2) { - ss << " and "; - } - count++; - } - return ss.str(); - }()); - - if (indices.QSESpeciesIndices.empty()) { - LOG_DEBUG(m_logger, "No QSE species to solve for."); - return Eigen::VectorXd(0); - } - - Eigen::VectorXd YScale(indices.QSESpeciesIndices.size()); - const double abundanceFloor = 1e-15; - for (size_t i = 0; i < indices.QSESpeciesIndices.size(); i++) { - const double initial_abundance = Y[indices.QSESpeciesIndices[i]]; - YScale(i) = std::max(initial_abundance, abundanceFloor); - } - - // Use the EigenFunctor with Eigen's nonlinear solver - EigenFunctor functor( - m_engine, - Y, - indices.dynamicSpeciesIndices, - indices.QSESpeciesIndices, - T9, - rho, - YScale - ); - - Eigen::VectorXd v_qse_initial(indices.QSESpeciesIndices.size()); - for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) { - const double initial_abundance = Y[indices.QSESpeciesIndices[i]]; - v_qse_initial(i) = std::asinh(initial_abundance/ YScale(i)); // Use asinh for better numerical stability compared to log - } - LOG_TRACE_L1( - m_logger, - "v_qse_log_initial: {}", - [&]() -> std::string { - std::stringstream ss; - ss << std::scientific << std::setprecision(5); - for (size_t i = 0; i < v_qse_initial.size(); ++i) { - ss << "log(X_" << std::string(m_engine.getNetworkSpecies()[indices.QSESpeciesIndices[i]].name()) << ") = " << v_qse_initial(i); - if (i < v_qse_initial.size() - 2) { - ss << ", "; - } else if (i == v_qse_initial.size() - 2) { - ss << " and "; - } - } - return ss.str(); - }()); - - const static std::unordered_map statusMessages = { - {Eigen::LevenbergMarquardtSpace::NotStarted, "Not started"}, - {Eigen::LevenbergMarquardtSpace::Running, "Running"}, - {Eigen::LevenbergMarquardtSpace::ImproperInputParameters, "Improper input parameters"}, - {Eigen::LevenbergMarquardtSpace::RelativeReductionTooSmall, "Relative reduction too small"}, - {Eigen::LevenbergMarquardtSpace::RelativeErrorTooSmall, "Relative error too small"}, - {Eigen::LevenbergMarquardtSpace::RelativeErrorAndReductionTooSmall, "Relative error and reduction too small"}, - {Eigen::LevenbergMarquardtSpace::CosinusTooSmall, "Cosine too small"}, - {Eigen::LevenbergMarquardtSpace::TooManyFunctionEvaluation, "Too many function evaluations"}, - {Eigen::LevenbergMarquardtSpace::FtolTooSmall, "Function tolerance too small"}, - {Eigen::LevenbergMarquardtSpace::XtolTooSmall, "X tolerance too small"}, - {Eigen::LevenbergMarquardtSpace::GtolTooSmall, "Gradient tolerance too small"}, - {Eigen::LevenbergMarquardtSpace::UserAsked, "User asked to stop"} - }; - - Eigen::LevenbergMarquardt lm(functor); - lm.parameters.xtol = 1e-24; - lm.parameters.ftol = 1e-24; - const Eigen::LevenbergMarquardtSpace::Status info = lm.minimize(v_qse_initial); - - if (info <= 0 || info >= 4) { - LOG_ERROR(m_logger, "QSE species minimization failed with status: {} ({})", - static_cast(info), statusMessages.at(info)); - throw std::runtime_error( - "QSE species minimization failed with status: " + std::to_string(static_cast(info)) + - " (" + std::string(statusMessages.at(info)) + ")" - ); - } - LOG_DEBUG(m_logger, "QSE species minimization completed successfully with status: {} ({})", - static_cast(info), statusMessages.at(info)); - return YScale.array() * v_qse_initial.array().sinh(); // Convert back to molar abundances - - } - - NetOut QSENetworkSolver::initializeNetworkWithShortIgnition(const NetIn &netIn) const { - const auto ignitionTemperature = m_config.get( - "gridfire:solver:QSE:ignition:temperature", - 2e8 - ); // 0.2 GK - const auto ignitionDensity = m_config.get( - "gridfire:solver:QSE:ignition:density", - 1e6 - ); // 1e6 g/cm^3 - const auto ignitionTime = m_config.get( - "gridfire:solver:QSE:ignition:tMax", - 1e-7 - ); // 0.1 μs - const auto ignitionStepSize = m_config.get( - "gridfire:solver:QSE:ignition:dt0", - 1e-15 - ); // 1e-15 seconds - - LOG_INFO( - m_logger, - "Igniting network with T={:<5.3E}, ρ={:<5.3E}, tMax={:<5.3E}, dt0={:<5.3E}...", - ignitionTemperature, - ignitionDensity, - ignitionTime, - ignitionStepSize - ); - LOG_INFO( - m_logger, - "Pre-ignition composition: {}", - [&]() -> std::string { - std::stringstream ss; - ss << std::scientific << std::setprecision(3); - int count = 0; - for (const auto& entry : netIn.composition | std::views::values) { - ss << "X_" << std::string(entry.isotope().name()) << ": " << entry.mass_fraction() << " "; - if (count < netIn.composition.getRegisteredSymbols().size() - 2) { - ss << ", "; - } else if (count == netIn.composition.getRegisteredSymbols().size() - 2) { - ss << " and "; - } - count++; - } - return ss.str(); - }()); - - NetIn preIgnition = netIn; - preIgnition.temperature = ignitionTemperature; - preIgnition.density = ignitionDensity; - preIgnition.tMax = ignitionTime; - preIgnition.dt0 = ignitionStepSize; - - const auto prevScreeningModel = m_engine.getScreeningModel(); - LOG_DEBUG(m_logger, "Setting screening model to BARE for high temperature and density ignition."); - m_engine.setScreeningModel(screening::ScreeningType::BARE); - DirectNetworkSolver ignitionSolver(m_engine); - NetOut postIgnition = ignitionSolver.evaluate(preIgnition); - LOG_INFO(m_logger, "Network ignition completed in {} steps.", postIgnition.num_steps); - LOG_INFO( - m_logger, - "Post-ignition composition: {}", - [&]() -> std::string { - std::stringstream ss; - ss << std::scientific << std::setprecision(3); - int count = 0; - for (const auto& entry : postIgnition.composition | std::views::values) { - ss << "X_" << std::string(entry.isotope().name()) << ": " << entry.mass_fraction() << " "; - if (count < postIgnition.composition.getRegisteredSymbols().size() - 2) { - ss << ", "; - } else if (count == postIgnition.composition.getRegisteredSymbols().size() - 2) { - ss << " and "; - } - count++; - } - return ss.str(); - }()); - LOG_DEBUG( - m_logger, - "Average change in composition during ignition: {}", - [&]() -> std::string { - std::stringstream ss; - ss << std::scientific << std::setprecision(3); - for (const auto& entry : postIgnition.composition | std::views::values) { - if (!postIgnition.composition.contains(entry.isotope()) || !netIn.composition.contains(entry.isotope())) { - ss << std::string(entry.isotope().name()) << ": inf %, "; - continue; - } - const double initialAbundance = netIn.composition.getMolarAbundance(std::string(entry.isotope().name())); - const double finalAbundance = postIgnition.composition.getMolarAbundance(std::string(entry.isotope().name())); - const double change = (finalAbundance - initialAbundance) / initialAbundance * 100.0; // Percentage change - ss << std::string(entry.isotope().name()) << ": " << change << " %, "; - } - return ss.str(); - }()); - m_engine.setScreeningModel(prevScreeningModel); - LOG_DEBUG(m_logger, "Restoring previous screening model: {}", static_cast(prevScreeningModel)); - return postIgnition; - } - - 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, - double t - ) const { - LOG_TRACE_L1( - m_logger, - "Timestepping at t={:0.3E} (dt={:0.3E}) with T9={:0.3E}, ρ={:0.3E}", - t, - t - s_prevTimestep, - m_T9, - m_rho - ); - s_prevTimestep = t; - // --- Populate the slow / dynamic species vector --- - std::vector YFull(m_engine.getNetworkSpecies().size()); - for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) { - YFull[m_dynamicSpeciesIndices[i]] = YDynamic(i); - } - - // --- Populate the QSE species vector --- - for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { - YFull[m_QSESpeciesIndices[i]] = m_Y_QSE(i); - } - - auto [full_dYdt, specificEnergyRate] = m_engine.calculateRHSAndEnergy(YFull, m_T9, m_rho); - - dYdtDynamic.resize(m_dynamicSpeciesIndices.size() + 1); - for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) { - dYdtDynamic(i) = full_dYdt[m_dynamicSpeciesIndices[i]]; - } - dYdtDynamic[m_dynamicSpeciesIndices.size()] = specificEnergyRate; - } - NetOut DirectNetworkSolver::evaluate(const NetIn &netIn) { namespace ublas = boost::numeric::ublas; namespace odeint = boost::numeric::odeint; @@ -526,44 +30,97 @@ namespace gridfire::solver { const auto absTol = m_config.get("gridfire:solver:DirectNetworkSolver:absTol", 1.0e-8); const auto relTol = m_config.get("gridfire:solver:DirectNetworkSolver:relTol", 1.0e-8); - size_t stepCount = 0; - Composition equilibratedComposition = m_engine.update(netIn); - const unsigned long numSpecies = m_engine.getNetworkSpecies().size(); + size_t numSpecies = m_engine.getNetworkSpecies().size(); ublas::vector Y(numSpecies + 1); - RHSFunctor rhsFunctor(m_engine, T9, netIn.density); + RHSManager manager(m_engine, T9, netIn.density); JacobianFunctor jacobianFunctor(m_engine, T9, netIn.density); + auto populateY = [&](const Composition& comp) { + const size_t numSpeciesInternal = m_engine.getNetworkSpecies().size(); + Y.resize(numSpeciesInternal + 1); + for (size_t i = 0; i < numSpeciesInternal; i++) { + const auto& species = m_engine.getNetworkSpecies()[i]; + if (!comp.contains(species)) { + double lim = std::numeric_limits::min(); + LOG_DEBUG(m_logger, "Species '{}' not found in composition. Setting abundance to {:0.3E}.", species.name(), lim); + Y(i) = lim; // Species not in the composition, set to zero + } else { + Y(i) = comp.getMolarAbundance(species); + } + } + + // TODO: a good starting point to make the temperature, density, and energy self consistent would be to turn this into an accumulator + Y(numSpeciesInternal) = 0.0; // Specific energy rate, initialized to zero + }; // This is a quick debug that can be turned on. For solar code input parameters (T~1.5e7K, ρ~1.5e3 g/cm^3) this should be near 8e-17 // std::cout << "D/H: " << equilibratedComposition.getMolarAbundance("H-2") / equilibratedComposition.getMolarAbundance("H-1") << std::endl; - for (size_t i = 0; i < numSpecies; ++i) { - const auto& species = m_engine.getNetworkSpecies()[i]; - if (!equilibratedComposition.contains(species)) { - LOG_DEBUG(m_logger, "Species '{}' not found in composition. Setting abundance to 1e-99.", species.name()); - Y(i) = 1e-99; - continue; - } - double abun = equilibratedComposition.getMolarAbundance(species); - // if (abun == 0.0) { - // abun = 1e-99; // Avoid zero abundances - // } - Y(i) = abun; - } - Y(numSpecies) = 0.0; + populateY(equilibratedComposition); const auto stepper = odeint::make_controlled>(absTol, relTol); - adaptive_integrator_observer observer; - stepCount = odeint::integrate_adaptive( - stepper, - std::make_pair(rhsFunctor, jacobianFunctor), - Y, - 0.0, - netIn.tMax, - netIn.dt0, - observer - ); + + double current_time = 0.0; + double current_initial_timestep = netIn.dt0; + double accumulated_energy = 0.0; + // size_t total_update_stages_triggered = 0; + + while (current_time < netIn.tMax) { + try { + odeint::integrate_adaptive( + stepper, + std::make_pair(manager, jacobianFunctor), + Y, + current_time, + netIn.tMax, + current_initial_timestep, + [&](const auto& state, double t) { + current_time = t; + manager.observe(state, t); + } + ); + current_time = netIn.tMax; + } catch (const exceptions::StaleEngineTrigger &e) { + LOG_INFO(m_logger, "Catching StaleEngineTrigger at t = {:0.3E} with T9 = {:0.3E}, rho = {:0.3E}. Triggering update stage (last stage took {} steps).", current_time, T9, netIn.density, e.totalSteps()); + exceptions::StaleEngineTrigger::state staleState = e.getState(); + accumulated_energy += e.energy(); // Add the specific energy rate to the accumulated energy + // total_update_stages_triggered++; + + Composition temp_comp; + std::vector mass_fractions; + size_t num_species_at_stop = e.numSpecies(); + + if (num_species_at_stop != m_engine.getNetworkSpecies().size()) { + throw std::runtime_error( + "StaleEngineError state has a different number of species than the engine. This should not happen." + ); + } + mass_fractions.reserve(num_species_at_stop); + + for (size_t i = 0; i < num_species_at_stop; ++i) { + const auto& species = m_engine.getNetworkSpecies()[i]; + temp_comp.registerSpecies(species); + mass_fractions.push_back(e.getMolarAbundance(i) * species.mass()); // Convert from molar abundance to mass fraction + } + temp_comp.setMassFraction(m_engine.getNetworkSpecies(), mass_fractions); + temp_comp.finalize(true); + + NetIn netInTemp = netIn; + netInTemp.temperature = e.temperature(); + netInTemp.density = e.density(); + netInTemp.composition = std::move(temp_comp); + + Composition currentComposition = m_engine.update(netInTemp); + populateY(currentComposition); + Y(Y.size() - 1) = e.energy(); // Set the specific energy rate from the stale state + numSpecies = m_engine.getNetworkSpecies().size(); + + // current_initial_timestep = 0.001 * manager.m_last_step_time; // set the new timestep to the last successful timestep before repartitioning + } + } + + accumulated_energy += Y(Y.size() - 1); // Add the specific energy rate to the accumulated energy std::vector finalMassFractions(numSpecies); for (size_t i = 0; i < numSpecies; ++i) { @@ -586,72 +143,73 @@ namespace gridfire::solver { NetOut netOut; netOut.composition = std::move(outputComposition); - netOut.energy = Y(numSpecies); // Specific energy rate - netOut.num_steps = stepCount; + netOut.energy = accumulated_energy; // Specific energy rate + netOut.num_steps = manager.m_num_steps; return netOut; } - void DirectNetworkSolver::RHSFunctor::operator()( + void DirectNetworkSolver::RHSManager::operator()( const boost::numeric::ublas::vector &Y, boost::numeric::ublas::vector &dYdt, const double t ) const { - const std::vector y(Y.begin(), m_numSpecies + Y.begin()); - // const auto timescales = m_engine.getSpeciesTimescales(y, m_T9, m_rho); - StepDerivatives deriv; - - // TODO: Refactor this to use std::expected - try { - deriv = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho); - } catch (const exceptions::StaleEngineError& e) { // If the engine reports that it is stale, we need to update it - LOG_INFO(m_logger, "Engine reports stale state. Calling engine update method..."); - - NetIn netInTemp; - fourdst::composition::Composition compositionTemp; - - for (const auto& species : m_engine.getNetworkSpecies()) { - compositionTemp.registerSpecies(species); - } - - std::vector X(y.size(), 0.0); - for (size_t i = 0; i < y.size(); ++i) { - double Xi = y[i] * m_engine.getNetworkSpecies()[i].mass(); // Convert from molar abundance to mass fraction - if (Xi < 0 && std::abs(Xi) < 1e-20) { - Xi = 0.0; // Avoid negative mass fractions - } - X[i] = Xi; - } - - compositionTemp.setMassFraction(m_engine.getNetworkSpecies(), X); - compositionTemp.finalize(true); - - netInTemp.composition = std::move(compositionTemp); - netInTemp.temperature = m_T9 * 1e9; // Convert T9 back to Kelvin - netInTemp.density = m_rho; // Density in g/cm^3 - - m_engine.update(netInTemp); - - LOG_INFO(m_logger, "Engine update complete. Recalculating RHS and energy..."); - deriv = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho); - LOG_INFO(m_logger, "RHS and energy recalculated successfully!"); + const size_t numSpecies = m_engine.getNetworkSpecies().size(); + if (t != m_cached_time || !m_cached_result.has_value() || m_cached_result.value().dydt.size() != numSpecies + 1) { + compute_and_cache(Y, t); } + const auto&[dydt, nuclearEnergyGenerationRate] = m_cached_result.value(); + dYdt.resize(numSpecies + 1); + std::ranges::copy(dydt, dYdt.begin()); + dYdt(numSpecies) = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate + } - dYdt.resize(m_numSpecies + 1); - // for (size_t i = 0; i < m_numSpecies; ++i) { - // std::cout << "Species: " << m_engine.getNetworkSpecies()[i].name() << ", dYdt: " << deriv.dydt[i] << '\n'; - // } + void DirectNetworkSolver::RHSManager::observe( + const boost::numeric::ublas::vector &state, + const double t + ) const { + double dt = t - m_last_observed_time; + compute_and_cache(state, t); + LOG_INFO( + m_logger, + "(Step {}) Observed state at t = {:0.3E} (dt = {:0.3E})", + m_num_steps, + t, + dt + ); + std::ostringstream oss; + oss << std::scientific << std::setprecision(3); + oss << "(Step: " << std::setw(10) << m_num_steps << ") t = " << t << " (dt = " << dt << ", eps_nuc: " << state(state.size() - 1) << " [erg])\n"; + std::cout << oss.str(); + m_last_observed_time = t; + m_last_step_time = dt; - // std::vector S; - // S.reserve(m_numSpecies); - // for (size_t i = 0; i < m_numSpecies; ++i) { - // S.push_back(std::abs(deriv.dydt[i]) / (1e-8 + 1e-8 * y[i])); - // } - // for (size_t i = 0; i < m_numSpecies; ++i) { - // std::cout << "Species: " << m_engine.getNetworkSpecies()[i].name() << ", S: " << S[i] << '\n'; - // } - std::ranges::copy(deriv.dydt, dYdt.begin()); - dYdt(m_numSpecies) = deriv.nuclearEnergyGenerationRate; + } + + void DirectNetworkSolver::RHSManager::compute_and_cache( + const boost::numeric::ublas::vector &state, + double t + ) const { + std::vector y_vec(state.begin(), state.end() - 1); + std::ranges::replace_if( + y_vec, + [](const double yi){ + return yi < 0.0; + }, + 0.0 // Avoid negative abundances + ); + + const auto result = m_engine.calculateRHSAndEnergy(y_vec, m_T9, m_rho); + if (!result) { + LOG_INFO(m_logger, + "Triggering update stage due to stale engine detected at t = {:0.3E} with T9 = {:0.3E}, rho = {:0.3E}, y_vec (size: {})", + t, m_T9, m_rho, y_vec.size()); + throw exceptions::StaleEngineTrigger({m_T9, m_rho, y_vec, t, m_num_steps, state(state.size() - 1)}); + } + m_cached_result = result.value(); + m_cached_time = t; + + m_num_steps++; } void DirectNetworkSolver::JacobianFunctor::operator()( @@ -660,13 +218,14 @@ namespace gridfire::solver { double t, boost::numeric::ublas::vector &dfdt ) const { - J.resize(m_numSpecies+1, m_numSpecies+1); + size_t numSpecies = m_engine.getNetworkSpecies().size(); + J.resize(numSpecies+1, numSpecies+1); J.clear(); - for (int i = 0; i < m_numSpecies; ++i) { - for (int j = 0; j < m_numSpecies; ++j) { + for (int i = 0; i < numSpecies; ++i) { + for (int j = 0; j < numSpecies; ++j) { J(i, j) = m_engine.getJacobianMatrixEntry(i, j); - // std::cout << "J(" << m_engine.getNetworkSpecies()[i].name() << ", " << m_engine.getNetworkSpecies()[j].name() << ") = " << J(i, j) << '\n'; } } } + } \ No newline at end of file diff --git a/src/network/lib/utils/logging.cpp b/src/network/lib/utils/logging.cpp index 19b618da..af3f1cce 100644 --- a/src/network/lib/utils/logging.cpp +++ b/src/network/lib/utils/logging.cpp @@ -16,7 +16,13 @@ std::string gridfire::utils::formatNuclearTimescaleLogString( const double T9, const double rho ) { - auto const& timescales = engine.getSpeciesTimescales(Y, T9, rho); + auto const& result = engine.getSpeciesTimescales(Y, T9, rho); + if (!result) { + std::ostringstream ss; + ss << "Failed to get species timescales: " << result.error(); + return ss.str(); + } + const std::unordered_map& timescales = result.value(); // Figure out how wide the "Species" column needs to be: std::size_t maxNameLen = std::string_view("Species").size(); diff --git a/subprojects/libcomposition.wrap b/subprojects/libcomposition.wrap index 268c9605..cd3facd8 100644 --- a/subprojects/libcomposition.wrap +++ b/subprojects/libcomposition.wrap @@ -1,4 +1,4 @@ [wrap-git] url = https://github.com/4D-STAR/libcomposition.git -revision = v1.4.0 +revision = v1.4.1 depth = 1 \ No newline at end of file diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp index 27a71498..1f93c2c6 100644 --- a/tests/graphnet_sandbox/main.cpp +++ b/tests/graphnet_sandbox/main.cpp @@ -72,24 +72,30 @@ int main() { NetIn netIn; netIn.composition = composition; netIn.temperature = 1.5e7; - netIn.density = 1.5e3; + netIn.density = 1.6e2; netIn.energy = 0; - netIn.tMax = 3.546e17; + netIn.tMax = 3.1536e17; // ~ 10Gyr netIn.dt0 = 1e-12; - - GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder); - ReaclibEngine.setPrecomputation(true); - ReaclibEngine.setUseReverseReactions(false); - MultiscalePartitioningEngineView partitioningView(ReaclibEngine); - AdaptiveEngineView adaptiveView(partitioningView); - - solver::DirectNetworkSolver solver(adaptiveView); - NetOut netOut; - measure_execution_time([&](){netOut = solver.evaluate(netIn);}, "DirectNetworkSolver Evaluation"); - std::cout << "DirectNetworkSolver completed in " << netOut.num_steps << " steps.\n"; - std::cout << "Final composition:\n"; - for (const auto& [symbol, entry] : netOut.composition) { + for (const auto& [symbol, entry] : netIn.composition) { std::cout << symbol << ": " << entry.mass_fraction() << "\n"; } + // GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder); + // + // ReaclibEngine.setPrecomputation(true); + // ReaclibEngine.setUseReverseReactions(false); + // ReaclibEngine.setScreeningModel(screening::ScreeningType::WEAK); + // + // MultiscalePartitioningEngineView partitioningView(ReaclibEngine); + // AdaptiveEngineView adaptiveView(partitioningView); + // + // solver::DirectNetworkSolver solver(adaptiveView); + // NetOut netOut; + // measure_execution_time([&](){netOut = solver.evaluate(netIn);}, "DirectNetworkSolver Evaluation"); + // std::cout << "DirectNetworkSolver completed in " << netOut.num_steps << " steps.\n"; + // std::cout << "Final composition:\n"; + // for (const auto& [symbol, entry] : netOut.composition) { + // std::cout << symbol << ": " << entry.mass_fraction() << "\n"; + // } + } \ No newline at end of file