diff --git a/src/network/include/gridfire/engine/engine_abstract.h b/src/network/include/gridfire/engine/engine_abstract.h index 9aed9628..a99ad8b2 100644 --- a/src/network/include/gridfire/engine/engine_abstract.h +++ b/src/network/include/gridfire/engine/engine_abstract.h @@ -1,6 +1,9 @@ #pragma once #include "gridfire/reaction/reaction.h" +#include "gridfire/network.h" +#include "gridfire/screening/screening_abstract.h" +#include "gridfire/screening/screening_types.h" #include #include @@ -80,7 +83,7 @@ namespace gridfire { * @brief Get the list of species in the network. * @return Vector of Species objects representing all network species. */ - virtual const std::vector& getNetworkSpecies() const = 0; + [[nodiscard]] virtual const std::vector& getNetworkSpecies() const = 0; /** * @brief Calculate the right-hand side (dY/dt) and energy generation. @@ -94,7 +97,7 @@ namespace gridfire { * time derivatives of all species and the specific nuclear energy generation * rate for the current state. */ - virtual StepDerivatives calculateRHSAndEnergy( + [[nodiscard]] virtual StepDerivatives calculateRHSAndEnergy( const std::vector& Y, double T9, double rho @@ -141,7 +144,7 @@ namespace gridfire { * * The Jacobian must have been generated by generateJacobianMatrix() before calling this. */ - virtual double getJacobianMatrixEntry( + [[nodiscard]] virtual double getJacobianMatrixEntry( int i, int j ) const = 0; @@ -163,7 +166,7 @@ namespace gridfire { * * The stoichiometry matrix must have been generated by generateStoichiometryMatrix(). */ - virtual int getStoichiometryMatrixEntry( + [[nodiscard]] virtual int getStoichiometryMatrixEntry( int speciesIndex, int reactionIndex ) const = 0; @@ -180,7 +183,7 @@ namespace gridfire { * This method computes the net rate at which the given reaction proceeds * under the current state. */ - virtual double calculateMolarReactionFlow( + [[nodiscard]] virtual double calculateMolarReactionFlow( const reaction::Reaction& reaction, const std::vector& Y, double T9, @@ -192,7 +195,7 @@ namespace gridfire { * * @return Reference to the LogicalReactionSet containing all reactions. */ - virtual const reaction::LogicalReactionSet& getNetworkReactions() const = 0; + [[nodiscard]] virtual const reaction::LogicalReactionSet& getNetworkReactions() const = 0; /** * @brief Compute timescales for all species in the network. @@ -205,10 +208,16 @@ 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. */ - virtual std::unordered_map getSpeciesTimescales( + [[nodiscard]] virtual std::unordered_map getSpeciesTimescales( const std::vector& Y, double T9, double rho ) const = 0; + + virtual void update(const NetIn& netIn) = 0; + + virtual void setScreeningModel(screening::ScreeningType model) = 0; + + [[nodiscard]] virtual screening::ScreeningType getScreeningModel() const = 0; }; } \ No newline at end of file diff --git a/src/network/include/gridfire/engine/engine_graph.h b/src/network/include/gridfire/engine/engine_graph.h index 0c712e41..2d0c6e76 100644 --- a/src/network/include/gridfire/engine/engine_graph.h +++ b/src/network/include/gridfire/engine/engine_graph.h @@ -8,10 +8,13 @@ #include "gridfire/network.h" #include "gridfire/reaction/reaction.h" #include "gridfire/engine/engine_abstract.h" +#include "gridfire/screening/screening_abstract.h" +#include "gridfire/screening/screening_types.h" #include #include #include +#include #include @@ -120,7 +123,7 @@ namespace gridfire { * * @see StepDerivatives */ - StepDerivatives calculateRHSAndEnergy( + [[nodiscard]] StepDerivatives calculateRHSAndEnergy( const std::vector& Y, const double T9, const double rho @@ -165,7 +168,7 @@ namespace gridfire { * This method computes the net rate at which the given reaction proceeds * under the current state. */ - double calculateMolarReactionFlow( + [[nodiscard]] double calculateMolarReactionFlow( const reaction::Reaction& reaction, const std::vector&Y, const double T9, @@ -243,6 +246,8 @@ namespace gridfire { double rho ) const override; + void update(const NetIn& netIn) override; + /** * @brief Checks if a given species is involved in the network. * @@ -293,6 +298,10 @@ namespace gridfire { const std::string& filename ) const; + void setScreeningModel(screening::ScreeningType) override; + + [[nodiscard]] screening::ScreeningType getScreeningModel() const override; + private: reaction::LogicalReactionSet m_reactions; ///< Set of REACLIB reactions in the network. @@ -307,6 +316,9 @@ namespace gridfire { CppAD::ADFun m_rhsADFun; ///< CppAD function for the right-hand side of the ODE. + screening::ScreeningType m_screeningType = screening::ScreeningType::BARE; ///< Screening type for the reaction network. Default to no screening. + std::unique_ptr m_screeningModel = screening::selectScreeningModel(m_screeningType); + Config& m_config = Config::getInstance(); Constants& m_constants = Constants::getInstance(); ///< Access to physical constants. quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); @@ -427,7 +439,7 @@ namespace gridfire { * specific nuclear energy generation rate for the current state. */ template - StepDerivatives calculateAllDerivatives( + [[nodiscard]] StepDerivatives calculateAllDerivatives( const std::vector &Y_in, T T9, T rho @@ -445,7 +457,7 @@ namespace gridfire { * specific nuclear energy generation rate for the current state using * double precision arithmetic. */ - StepDerivatives calculateAllDerivatives( + [[nodiscard]] StepDerivatives calculateAllDerivatives( const std::vector& Y_in, const double T9, const double rho @@ -463,7 +475,7 @@ namespace gridfire { * specific nuclear energy generation rate for the current state using * automatic differentiation. */ - StepDerivatives calculateAllDerivatives( + [[nodiscard]] StepDerivatives calculateAllDerivatives( const std::vector& Y_in, const ADDouble &T9, const ADDouble &rho @@ -474,6 +486,13 @@ namespace gridfire { template StepDerivatives GraphEngine::calculateAllDerivatives( const std::vector &Y_in, T T9, T rho) const { + std::vector screeningFactors = m_screeningModel->calculateScreeningFactors( + m_reactions, + m_networkSpecies, + Y_in, + T9, + rho + ); // --- Setup output derivatives structure --- StepDerivatives result; @@ -512,7 +531,7 @@ namespace gridfire { const auto& reaction = m_reactions[reactionIndex]; // 1. Calculate reaction rate - const T molarReactionFlow = calculateMolarReactionFlow(reaction, Y, T9, rho); + const T molarReactionFlow = screeningFactors[reactionIndex] * calculateMolarReactionFlow(reaction, Y, T9, rho); // 2. Use the rate to update all relevant species derivatives (dY/dt) for (size_t speciesIndex = 0; speciesIndex < m_networkSpecies.size(); ++speciesIndex) { diff --git a/src/network/include/gridfire/engine/views/engine_adaptive.h b/src/network/include/gridfire/engine/views/engine_adaptive.h index 34d97f2a..1cc8cb2f 100644 --- a/src/network/include/gridfire/engine/views/engine_adaptive.h +++ b/src/network/include/gridfire/engine/views/engine_adaptive.h @@ -1,6 +1,8 @@ #pragma once #include "gridfire/engine/engine_abstract.h" -#include "engine_view_abstract.h" +#include "gridfire/engine/views/engine_view_abstract.h" +#include "gridfire/screening/screening_abstract.h" +#include "gridfire/screening/screening_types.h" #include "gridfire/network.h" #include "fourdst/composition/atomicSpecies.h" @@ -71,13 +73,13 @@ namespace gridfire { * @see AdaptiveEngineView::constructSpeciesIndexMap() * @see AdaptiveEngineView::constructReactionIndexMap() */ - void update(const NetIn& netIn); + void update(const NetIn& netIn) override; /** * @brief Gets the list of active species in the network. * @return A const reference to the vector of active species. */ - const std::vector& getNetworkSpecies() const override; + [[nodiscard]] const std::vector& getNetworkSpecies() const override; /** * @brief Calculates the right-hand side (dY/dt) and energy generation for the active species. @@ -95,7 +97,7 @@ namespace gridfire { * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). * @see AdaptiveEngineView::update() */ - StepDerivatives calculateRHSAndEnergy( + [[nodiscard]] StepDerivatives calculateRHSAndEnergy( const std::vector &Y_culled, const double T9, const double rho @@ -134,7 +136,7 @@ namespace gridfire { * @throws std::out_of_range If the culled index is out of bounds for the species index map. * @see AdaptiveEngineView::update() */ - double getJacobianMatrixEntry( + [[nodiscard]] double getJacobianMatrixEntry( const int i_culled, const int j_culled ) const override; @@ -164,7 +166,7 @@ namespace gridfire { * @throws std::out_of_range If the culled index is out of bounds for the species or reaction index map. * @see AdaptiveEngineView::update() */ - int getStoichiometryMatrixEntry( + [[nodiscard]] int getStoichiometryMatrixEntry( const int speciesIndex_culled, const int reactionIndex_culled ) const override; @@ -184,7 +186,7 @@ namespace gridfire { * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). * @throws std::runtime_error If the reaction is not part of the active reactions in the adaptive engine view. */ - double calculateMolarReactionFlow( + [[nodiscard]] double calculateMolarReactionFlow( const reaction::Reaction &reaction, const std::vector &Y_culled, double T9, @@ -196,7 +198,7 @@ namespace gridfire { * * @return Reference to the LogicalReactionSet containing all active reactions. */ - const reaction::LogicalReactionSet& getNetworkReactions() const override; + [[nodiscard]] const reaction::LogicalReactionSet& getNetworkReactions() const override; /** * @brief Computes timescales for all active species in the network. @@ -211,7 +213,7 @@ namespace gridfire { * * @throws std::runtime_error If the AdaptiveEngineView is stale (i.e., `update()` has not been called). */ - std::unordered_map getSpeciesTimescales( + [[nodiscard]] std::unordered_map getSpeciesTimescales( const std::vector &Y_culled, double T9, double rho @@ -221,7 +223,11 @@ namespace gridfire { * @brief Gets the base engine. * @return A const reference to the base engine. */ - const DynamicEngine& getBaseEngine() const override { return m_baseEngine; } + [[nodiscard]] const DynamicEngine& getBaseEngine() const override { return m_baseEngine; } + + void setScreeningModel(screening::ScreeningType model) override; + + [[nodiscard]] screening::ScreeningType getScreeningModel() const override; private: using Config = fourdst::config::Config; using LogManager = fourdst::logging::LogManager; @@ -257,7 +263,7 @@ namespace gridfire { * * @see AdaptiveEngineView::update() */ - std::vector constructSpeciesIndexMap() const; + [[nodiscard]] std::vector constructSpeciesIndexMap() const; /** * @brief Constructs the reaction index map. @@ -269,7 +275,7 @@ namespace gridfire { * * @see AdaptiveEngineView::update() */ - std::vector constructReactionIndexMap() const; + [[nodiscard]] std::vector constructReactionIndexMap() const; /** * @brief Maps a vector of culled abundances to a vector of full abundances. @@ -278,7 +284,7 @@ namespace gridfire { * @return A vector of abundances for the full network, with the abundances of the active * species copied from the culled vector. */ - std::vector mapCulledToFull(const std::vector& culled) const; + [[nodiscard]] std::vector mapCulledToFull(const std::vector& culled) const; /** * @brief Maps a vector of full abundances to a vector of culled abundances. @@ -287,7 +293,7 @@ namespace gridfire { * @return A vector of abundances for the active species, with the abundances of the active * species copied from the full vector. */ - std::vector mapFullToCulled(const std::vector& full) const; + [[nodiscard]] std::vector mapFullToCulled(const std::vector& full) const; /** * @brief Maps a culled species index to a full species index. @@ -297,7 +303,7 @@ namespace gridfire { * * @throws std::out_of_range If the culled index is out of bounds for the species index map. */ - size_t mapCulledToFullSpeciesIndex(size_t culledSpeciesIndex) const; + [[nodiscard]] size_t mapCulledToFullSpeciesIndex(size_t culledSpeciesIndex) const; /** * @brief Maps a culled reaction index to a full reaction index. @@ -307,7 +313,7 @@ namespace gridfire { * * @throws std::out_of_range If the culled index is out of bounds for the reaction index map. */ - size_t mapCulledToFullReactionIndex(size_t culledReactionIndex) const; + [[nodiscard]] size_t mapCulledToFullReactionIndex(size_t culledReactionIndex) const; /** * @brief Validates that the AdaptiveEngineView is not stale. @@ -320,10 +326,10 @@ namespace gridfire { const NetIn& netIn, std::vector& out_Y_Full ) const; - std::unordered_set findReachableSpecies( + [[nodiscard]] std::unordered_set findReachableSpecies( const NetIn& netIn ) const; - std::vector cullReactionsByFlow( + [[nodiscard]] std::vector cullReactionsByFlow( const std::vector& allFlows, const std::unordered_set& reachableSpecies, const std::vector& Y_full, diff --git a/src/network/include/gridfire/engine/views/engine_defined.h b/src/network/include/gridfire/engine/views/engine_defined.h index b5aae428..90d61ffc 100644 --- a/src/network/include/gridfire/engine/views/engine_defined.h +++ b/src/network/include/gridfire/engine/views/engine_defined.h @@ -1,7 +1,9 @@ #pragma once -#include "engine_view_abstract.h" -#include "../engine_abstract.h" +#include "gridfire/engine/views/engine_view_abstract.h" +#include "gridfire/engine/engine_abstract.h" +#include "gridfire/io/network_file.h" +#include "gridfire/network.h" #include "fourdst/config/config.h" #include "fourdst/logging/logging.h" @@ -13,7 +15,11 @@ namespace gridfire{ class FileDefinedEngineView final: public DynamicEngine, public EngineView { public: - explicit FileDefinedEngineView(DynamicEngine& baseEngine, const std::string& fileName); + explicit FileDefinedEngineView( + DynamicEngine& baseEngine, + const std::string& fileName, + const io::NetworkFileParser& parser + ); // --- EngineView Interface --- const DynamicEngine& getBaseEngine() const override; @@ -53,6 +59,14 @@ namespace gridfire{ const double T9, const double rho ) const override; + + void update(const NetIn &netIn) override; + + void setNetworkFile(const std::string& fileName); + + void setScreeningModel(screening::ScreeningType model) override; + + [[nodiscard]] screening::ScreeningType getScreeningModel() const override; private: using Config = fourdst::config::Config; using LogManager = fourdst::logging::LogManager; @@ -60,12 +74,17 @@ namespace gridfire{ quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); DynamicEngine& m_baseEngine; + std::string m_fileName; ///< Name of the file defining the reaction set considered by the engine view. + const io::NetworkFileParser& m_parser; ///< Parser for the network file. std::vector m_activeSpecies; ///< Active species in the defined engine. reaction::LogicalReactionSet m_activeReactions; ///< Active reactions in the defined engine. std::vector m_speciesIndexMap; ///< Maps indices of active species to indices in the full network. std::vector m_reactionIndexMap; ///< Maps indices of active reactions to indices in the full network. + + bool m_isStale = true; + private: void buildFromFile(const std::string& fileName); @@ -100,7 +119,7 @@ namespace gridfire{ * @return A vector of abundances for the full network, with the abundances of the active * species copied from the culled vector. */ - std::vector mapCulledToFull(const std::vector& culled) const; + std::vector mapViewToFull(const std::vector& culled) const; /** * @brief Maps a vector of full abundances to a vector of culled abundances. @@ -109,7 +128,7 @@ namespace gridfire{ * @return A vector of abundances for the active species, with the abundances of the active * species copied from the full vector. */ - std::vector mapFullToCulled(const std::vector& full) const; + std::vector mapFullToView(const std::vector& full) const; /** * @brief Maps a culled species index to a full species index. @@ -119,7 +138,7 @@ namespace gridfire{ * * @throws std::out_of_range If the culled index is out of bounds for the species index map. */ - size_t mapCulledToFullSpeciesIndex(size_t culledSpeciesIndex) const; + size_t mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const; /** * @brief Maps a culled reaction index to a full reaction index. @@ -129,7 +148,8 @@ namespace gridfire{ * * @throws std::out_of_range If the culled index is out of bounds for the reaction index map. */ - size_t mapCulledToFullReactionIndex(size_t culledReactionIndex) const; + size_t mapViewToFullReactionIndex(size_t culledReactionIndex) const; + void validateNetworkState() const; }; } \ No newline at end of file diff --git a/src/network/include/gridfire/io/network_file.h b/src/network/include/gridfire/io/network_file.h index e69de29b..4ed81f00 100644 --- a/src/network/include/gridfire/io/network_file.h +++ b/src/network/include/gridfire/io/network_file.h @@ -0,0 +1,48 @@ +#pragma once + +#include "fourdst/config/config.h" +#include "fourdst/logging/logging.h" + +#include "quill/Logger.h" + +#include +#include + +namespace gridfire::io { + + struct ParsedNetworkData { + std::vector reactionPENames; + }; + + class NetworkFileParser { + public: + virtual ~NetworkFileParser() = default; + + [[nodiscard]] virtual ParsedNetworkData parse(const std::string& filename) const = 0; + + }; + + class SimpleReactionListFileParser final : public NetworkFileParser { + public: + explicit SimpleReactionListFileParser(); + ParsedNetworkData parse(const std::string& filename) const override; + private: + using Config = fourdst::config::Config; + using LogManager = fourdst::logging::LogManager; + Config& m_config = Config::getInstance(); + quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); + }; + + class MESANetworkFileParser final : public NetworkFileParser { + public: + explicit MESANetworkFileParser(const std::string& filename); + ParsedNetworkData parse(const std::string& filename) const override; + private: + using Config = fourdst::config::Config; + using LogManager = fourdst::logging::LogManager; + Config& m_config = Config::getInstance(); + quill::Logger* m_logger = LogManager::getInstance().getLogger("log"); + + std::string m_filename; + }; +} \ No newline at end of file diff --git a/src/network/include/gridfire/reaction/reaction.h b/src/network/include/gridfire/reaction/reaction.h index 6ed7a6db..b51840e7 100644 --- a/src/network/include/gridfire/reaction/reaction.h +++ b/src/network/include/gridfire/reaction/reaction.h @@ -11,6 +11,7 @@ #include "cppad/cppad.hpp" +#include "xxhash64.h" /** * @file reaction.h @@ -252,6 +253,10 @@ namespace gridfire::reaction { */ [[nodiscard]] uint64_t hash(uint64_t seed = 0) const; + friend std::ostream& operator<<(std::ostream& os, const Reaction& r) { + return os << "(Reaction:" << r.m_id << ")"; + } + protected: quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); std::string m_id; ///< Unique identifier for the reaction (e.g., "h1+h1=>h2+e+nu"). @@ -288,162 +293,6 @@ namespace gridfire::reaction { } }; - /** - * @class ReactionSet - * @brief A collection of Reaction objects. - * - * This class manages a set of individual `Reaction` objects, providing - * efficient lookup by ID and functionality to query the entire set. - * - * Example: - * @code - * ReactionSet my_set({reaction1, reaction2}); - * my_set.add_reaction(reaction3); - * if (my_set.contains("h1(p,g)h2")) { - * const Reaction& r = my_set["h1(p,g)h2"]; - * } - * @endcode - */ - class ReactionSet { - public: - /** - * @brief Constructs a ReactionSet from a vector of reactions. - * @param reactions The initial vector of Reaction objects. - */ - explicit ReactionSet(std::vector reactions); - - /** - * @brief Copy constructor. - * @param other The ReactionSet to copy. - */ - ReactionSet(const ReactionSet& other); - - /** - * @brief Copy assignment operator. - * @param other The ReactionSet to assign from. - * @return A reference to this ReactionSet. - */ - ReactionSet& operator=(const ReactionSet& other); - - /** - * @brief Virtual destructor. - */ - virtual ~ReactionSet() = default; - - /** - * @brief Adds a reaction to the set. - * @param reaction The Reaction to add. - */ - virtual void add_reaction(Reaction reaction); - - /** - * @brief Removes a reaction from the set. - * @param reaction The Reaction to remove. - */ - virtual void remove_reaction(const Reaction& reaction); - - /** - * @brief Checks if the set contains a reaction with the given ID. - * @param id The ID of the reaction to find. - * @return True if the reaction is in the set, false otherwise. - */ - [[nodiscard]] bool contains(const std::string_view& id) const; - - /** - * @brief Checks if the set contains the given reaction. - * @param reaction The Reaction to find. - * @return True if the reaction is in the set, false otherwise. - */ - [[nodiscard]] bool contains(const Reaction& reaction) const; - - /** - * @brief Gets the number of reactions in the set. - * @return The size of the set. - */ - [[nodiscard]] virtual size_t size() const { return m_reactions.size(); } - - /** - * @brief Removes all reactions from the set. - */ - void clear(); - - /** - * @brief Checks if any reaction in the set involves the given species. - * @param species The species to check for. - * @return True if the species is involved in any reaction. - */ - [[nodiscard]] bool contains_species(const fourdst::atomic::Species& species) const; - - /** - * @brief Checks if any reaction in the set contains the given species as a reactant. - * @param species The species to check for. - * @return True if the species is a reactant in any reaction. - */ - [[nodiscard]] bool contains_reactant(const fourdst::atomic::Species& species) const; - - /** - * @brief Checks if any reaction in the set contains the given species as a product. - * @param species The species to check for. - * @return True if the species is a product in any reaction. - */ - [[nodiscard]] bool contains_product(const fourdst::atomic::Species& species) const; - - /** - * @brief Accesses a reaction by its index. - * @param index The index of the reaction to access. - * @return A const reference to the Reaction. - * @throws std::out_of_range if the index is out of bounds. - */ - [[nodiscard]] virtual const Reaction& operator[](size_t index) const; - - /** - * @brief Accesses a reaction by its ID. - * @param id The ID of the reaction to access. - * @return A const reference to the Reaction. - * @throws std::out_of_range if no reaction with the given ID exists. - */ - [[nodiscard]] const Reaction& operator[](const std::string_view& id) const; - - /** - * @brief Compares this set with another for equality. - * @param other The other ReactionSet to compare with. - * @return True if the sets are equal (same size and hash). - */ - bool operator==(const ReactionSet& other) const; - - /** - * @brief Compares this set with another for inequality. - * @param other The other ReactionSet to compare with. - * @return True if the sets are not equal. - */ - bool operator!=(const ReactionSet& other) const; - - /** - * @brief Computes a hash for the entire set. - * @param seed The seed for the hash function. - * @return A 64-bit hash value. - * @details The algorithm computes the hash of each individual reaction, - * sorts the hashes, and then computes a final hash over the sorted list - * of hashes. This ensures the hash is order-independent. - */ - [[nodiscard]] uint64_t hash(uint64_t seed = 0) const; - - /** @name Iterators - * Provides iterators to loop over the reactions in the set. - */ - ///@{ - auto begin() { return m_reactions.begin(); } - [[nodiscard]] auto begin() const { return m_reactions.cbegin(); } - auto end() { return m_reactions.end(); } - [[nodiscard]] auto end() const { return m_reactions.cend(); } - ///@} - private: - quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - std::vector m_reactions; - std::string m_id; - std::unordered_map m_reactionNameMap; ///< Maps reaction IDs to Reaction objects for quick lookup. - - }; /** @@ -508,6 +357,12 @@ namespace gridfire::reaction { auto end() { return m_rates.end(); } [[nodiscard]] auto end() const { return m_rates.cend(); } ///@} + /// + + friend std::ostream& operator<<(std::ostream& os, const LogicalReaction& r) { + os << "(LogicalReaction: " << r.id() << ", reverse: " << r.is_reverse() << ")"; + return os; + } private: std::vector m_sources; ///< List of source labels. @@ -543,31 +398,128 @@ namespace gridfire::reaction { } }; - /** - * @class LogicalReactionSet - * @brief A collection of LogicalReaction objects. - * - * This class takes a `ReactionSet` and groups individual `Reaction` objects - * into `LogicalReaction` objects based on their `peName`. This provides a - * view of the network where all rates for the same physical process are combined. - */ - class LogicalReactionSet final : public ReactionSet { + template + class TemplatedReactionSet final { public: /** - * @brief Deleted default constructor. + * @brief Constructs a ReactionSet from a vector of reactions. + * @param reactions The initial vector of Reaction objects. */ - LogicalReactionSet() = delete; + explicit TemplatedReactionSet(std::vector reactions); /** - * @brief Constructs a LogicalReactionSet from a ReactionSet. - * @param reactionSet The set of individual reactions to group. - * @details This constructor iterates through the provided `ReactionSet`, - * groups reactions by their `peName`, and creates a `LogicalReaction` for each group. + * @brief Copy constructor. + * @param other The ReactionSet to copy. */ - explicit LogicalReactionSet(const ReactionSet& reactionSet); + TemplatedReactionSet(const TemplatedReactionSet& other); + + /** + * @brief Copy assignment operator. + * @param other The ReactionSet to assign from. + * @return A reference to this ReactionSet. + */ + TemplatedReactionSet& operator=(const TemplatedReactionSet& other); + + /** + * @brief Adds a reaction to the set. + * @param reaction The Reaction to add. + */ + void add_reaction(ReactionT reaction); + + /** + * @brief Removes a reaction from the set. + * @param reaction The Reaction to remove. + */ + void remove_reaction(const ReactionT& reaction); + + /** + * @brief Checks if the set contains a reaction with the given ID. + * @param id The ID of the reaction to find. + * @return True if the reaction is in the set, false otherwise. + */ + [[nodiscard]] bool contains(const std::string_view& id) const; + + /** + * @brief Checks if the set contains the given reaction. + * @param reaction The Reaction to find. + * @return True if the reaction is in the set, false otherwise. + */ + [[nodiscard]] bool contains(const Reaction& reaction) const; + + /** + * @brief Gets the number of reactions in the set. + * @return The size of the set. + */ + [[nodiscard]] size_t size() const { return m_reactions.size(); } + + /** + * @brief Removes all reactions from the set. + */ + void clear(); + + /** + * @brief Checks if any reaction in the set involves the given species. + * @param species The species to check for. + * @return True if the species is involved in any reaction. + */ + [[nodiscard]] bool contains_species(const fourdst::atomic::Species& species) const; + + /** + * @brief Checks if any reaction in the set contains the given species as a reactant. + * @param species The species to check for. + * @return True if the species is a reactant in any reaction. + */ + [[nodiscard]] bool contains_reactant(const fourdst::atomic::Species& species) const; + + /** + * @brief Checks if any reaction in the set contains the given species as a product. + * @param species The species to check for. + * @return True if the species is a product in any reaction. + */ + [[nodiscard]] bool contains_product(const fourdst::atomic::Species& species) const; + + /** + * @brief Accesses a reaction by its index. + * @param index The index of the reaction to access. + * @return A const reference to the Reaction. + * @throws std::out_of_range if the index is out of bounds. + */ + [[nodiscard]] const ReactionT& operator[](size_t index) const; + + /** + * @brief Accesses a reaction by its ID. + * @param id The ID of the reaction to access. + * @return A const reference to the Reaction. + * @throws std::out_of_range if no reaction with the given ID exists. + */ + [[nodiscard]] const ReactionT& operator[](const std::string_view& id) const; + + /** + * @brief Compares this set with another for equality. + * @param other The other ReactionSet to compare with. + * @return True if the sets are equal (same size and hash). + */ + bool operator==(const TemplatedReactionSet& other) const; + + /** + * @brief Compares this set with another for inequality. + * @param other The other ReactionSet to compare with. + * @return True if the sets are not equal. + */ + bool operator!=(const TemplatedReactionSet& other) const; + + /** + * @brief Computes a hash for the entire set. + * @param seed The seed for the hash function. + * @return A 64-bit hash value. + * @details The algorithm computes the hash of each individual reaction, + * sorts the hashes, and then computes a final hash over the sorted list + * of hashes. This ensures the hash is order-independent. + */ + [[nodiscard]] uint64_t hash(uint64_t seed = 0) const; /** @name Iterators - * Provides iterators to loop over the logical reactions in the set. + * Provides iterators to loop over the reactions in the set. */ ///@{ auto begin() { return m_reactions.begin(); } @@ -575,25 +527,208 @@ namespace gridfire::reaction { auto end() { return m_reactions.end(); } [[nodiscard]] auto end() const { return m_reactions.cend(); } ///@} + /// + friend std::ostream& operator<<(std::ostream& os, const TemplatedReactionSet& r) { + os << "(ReactionSet: ["; + int counter = 0; + for (const auto& reaction : r.m_reactions) { + os << reaction; + if (counter < r.m_reactions.size() - 2) { + os << ", "; + } else if (counter == r.m_reactions.size() - 2) { + os << " and "; + } + ++counter; + } + os << "])"; + return os; + } - /** - * @brief Gets the number of logical reactions in the set. - * @return The size of the set. - */ - [[nodiscard]] size_t size() const { return m_reactions.size(); } - - /** - * @brief Accesses a logical reaction by its index. - * @param index The index of the logical reaction. - * @return A const reference to the LogicalReaction. - */ - [[nodiscard]] const LogicalReaction& operator[](size_t index) const { return m_reactions[index]; } + [[nodiscard]] std::unordered_set getReactionSetSpecies() const; private: quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - std::vector m_reactions; + std::vector m_reactions; std::string m_id; - std::unordered_map m_reactionNameMap; ///< Maps reaction IDs to LogicalReaction objects for quick lookup. + std::unordered_map m_reactionNameMap; ///< Maps reaction IDs to Reaction objects for quick lookup. + }; + using ReactionSet = TemplatedReactionSet; ///< A set of reactions, typically from a single source like REACLIB. + using LogicalReactionSet = TemplatedReactionSet; ///< A set of logical reactions. + + LogicalReactionSet packReactionSetToLogicalReactionSet(const ReactionSet& reactionSet); + + template + TemplatedReactionSet::TemplatedReactionSet( + std::vector reactions + ) : + m_reactions(std::move(reactions)) { + if (m_reactions.empty()) { + return; // Case where the reactions will be added later. + } + m_reactionNameMap.reserve(reactions.size()); + for (const auto& reaction : m_reactions) { + m_id += reaction.id(); + m_reactionNameMap.emplace(reaction.id(), reaction); + } + } + + template + TemplatedReactionSet::TemplatedReactionSet(const TemplatedReactionSet &other) { + m_reactions.reserve(other.m_reactions.size()); + for (const auto& reaction_ptr: other.m_reactions) { + m_reactions.push_back(reaction_ptr); + } + + m_reactionNameMap.reserve(other.m_reactionNameMap.size()); + for (const auto& reaction_ptr : m_reactions) { + m_reactionNameMap.emplace(reaction_ptr.id(), reaction_ptr); + } + } + + template + TemplatedReactionSet& TemplatedReactionSet::operator=(const TemplatedReactionSet &other) { + if (this != &other) { + TemplatedReactionSet temp(other); + std::swap(m_reactions, temp.m_reactions); + std::swap(m_reactionNameMap, temp.m_reactionNameMap); + } + return *this; + } + + template + void TemplatedReactionSet::add_reaction(ReactionT reaction) { + m_reactions.emplace_back(reaction); + m_id += m_reactions.back().id(); + m_reactionNameMap.emplace(m_reactions.back().id(), m_reactions.back()); + } + + template + void TemplatedReactionSet::remove_reaction(const ReactionT& reaction) { + if (!m_reactionNameMap.contains(std::string(reaction.id()))) { + return; + } + + m_reactionNameMap.erase(std::string(reaction.id())); + + std::erase_if(m_reactions, [&reaction](const Reaction& r) { + return r == reaction; + }); + } + + template + bool TemplatedReactionSet::contains(const std::string_view& id) const { + for (const auto& reaction : m_reactions) { + if (reaction.id() == id) { + return true; + } + } + return false; + } + + template + bool TemplatedReactionSet::contains(const Reaction& reaction) const { + for (const auto& r : m_reactions) { + if (r == reaction) { + return true; + } + } + return false; + } + + template + void TemplatedReactionSet::clear() { + m_reactions.clear(); + m_reactionNameMap.clear(); + } + + template + bool TemplatedReactionSet::contains_species(const fourdst::atomic::Species& species) const { + for (const auto& reaction : m_reactions) { + if (reaction.contains(species)) { + return true; + } + } + return false; + } + + template + bool TemplatedReactionSet::contains_reactant(const fourdst::atomic::Species& species) const { + for (const auto& r : m_reactions) { + if (r.contains_reactant(species)) { + return true; + } + } + return false; + } + + template + bool TemplatedReactionSet::contains_product(const fourdst::atomic::Species& species) const { + for (const auto& r : m_reactions) { + if (r.contains_product(species)) { + return true; + } + } + return false; + } + + template + const ReactionT& TemplatedReactionSet::operator[](const size_t index) const { + if (index >= m_reactions.size()) { + m_logger -> flush_log(); + throw std::out_of_range("Index" + std::to_string(index) + " out of range for ReactionSet of size " + std::to_string(m_reactions.size()) + "."); + } + return m_reactions[index]; + } + + template + const ReactionT& TemplatedReactionSet::operator[](const std::string_view& id) const { + if (auto it = m_reactionNameMap.find(std::string(id)); it != m_reactionNameMap.end()) { + return it->second; + } + m_logger -> flush_log(); + throw std::out_of_range("Species " + std::string(id) + " does not exist in ReactionSet."); + } + + template + bool TemplatedReactionSet::operator==(const TemplatedReactionSet& other) const { + if (size() != other.size()) { + return false; + } + return hash() == other.hash(); + } + + template + bool TemplatedReactionSet::operator!=(const TemplatedReactionSet& other) const { + return !(*this == other); + } + + template + uint64_t TemplatedReactionSet::hash(uint64_t seed) const { + if (m_reactions.empty()) { + return XXHash64::hash(nullptr, 0, seed); + } + std::vector individualReactionHashes; + individualReactionHashes.reserve(m_reactions.size()); + for (const auto& reaction : m_reactions) { + individualReactionHashes.push_back(reaction.hash(seed)); + } + + std::ranges::sort(individualReactionHashes); + + const auto data = static_cast(individualReactionHashes.data()); + const size_t sizeInBytes = individualReactionHashes.size() * sizeof(uint64_t); + return XXHash64::hash(data, sizeInBytes, seed); + } + + template + std::unordered_set TemplatedReactionSet::getReactionSetSpecies() const { + std::unordered_set species; + for (const auto& reaction : m_reactions) { + const auto reactionSpecies = reaction.all_species(); + species.insert(reactionSpecies.begin(), reactionSpecies.end()); + } + return species; + } } diff --git a/src/network/include/gridfire/screening/screening_abstract.h b/src/network/include/gridfire/screening/screening_abstract.h new file mode 100644 index 00000000..0529fcda --- /dev/null +++ b/src/network/include/gridfire/screening/screening_abstract.h @@ -0,0 +1,33 @@ +#pragma once + +#include "gridfire/reaction/reaction.h" + +#include "fourdst/composition/atomicSpecies.h" + +#include "cppad/cppad.hpp" + +#include + +namespace gridfire::screening { + class ScreeningModel { + public: + using ADDouble = CppAD::AD; + virtual ~ScreeningModel() = default; + + virtual std::vector calculateScreeningFactors( + const reaction::LogicalReactionSet& reactions, + const std::vector& species, + const std::vector& Y, + const double T9, + const double rho + ) const = 0; + + virtual std::vector calculateScreeningFactors( + const reaction::LogicalReactionSet& reactions, + const std::vector& species, + const std::vector& Y, + const ADDouble T9, + const ADDouble rho + ) const = 0; + }; +} \ No newline at end of file diff --git a/src/network/include/gridfire/screening/screening_bare.h b/src/network/include/gridfire/screening/screening_bare.h new file mode 100644 index 00000000..52bff705 --- /dev/null +++ b/src/network/include/gridfire/screening/screening_bare.h @@ -0,0 +1,48 @@ +#pragma once + +#include "gridfire/screening/screening_abstract.h" +#include "gridfire/reaction/reaction.h" + +#include "cppad/cppad.hpp" + +namespace gridfire::screening { + class BareScreeningModel final : public ScreeningModel { + using ADDouble = CppAD::AD; + public: + [[nodiscard]] std::vector calculateScreeningFactors( + const reaction::LogicalReactionSet& reactions, + const std::vector& species, + const std::vector& Y, + const double T9, + const double rho + ) const override; + + [[nodiscard]] std::vector calculateScreeningFactors( + const reaction::LogicalReactionSet& reactions, + const std::vector& species, + const std::vector& Y, + const ADDouble T9, + const ADDouble rho + ) const override; + private: + template + [[nodiscard]] std::vector calculateFactors_impl( + const reaction::LogicalReactionSet& reactions, + const std::vector& species, + const std::vector& Y, + const T T9, + const T rho + ) const; + }; + + template + std::vector BareScreeningModel::calculateFactors_impl( + const reaction::LogicalReactionSet &reactions, + const std::vector &species, + const std::vector &Y, + const T T9, + const T rho + ) const { + return std::vector(reactions.size(), T(1.0)); // Bare screening returns 1.0 for all reactions + } +} diff --git a/src/network/include/gridfire/screening/screening_types.h b/src/network/include/gridfire/screening/screening_types.h new file mode 100644 index 00000000..77c2168e --- /dev/null +++ b/src/network/include/gridfire/screening/screening_types.h @@ -0,0 +1,14 @@ +#pragma once + +#include "gridfire/screening/screening_abstract.h" + +#include + +namespace gridfire::screening { + enum class ScreeningType { + BARE, ///< No screening applied + WEAK, ///< Weak screening model + }; + + std::unique_ptr selectScreeningModel(ScreeningType type); +} \ No newline at end of file diff --git a/src/network/include/gridfire/screening/screening_weak.h b/src/network/include/gridfire/screening/screening_weak.h new file mode 100644 index 00000000..a1c1286a --- /dev/null +++ b/src/network/include/gridfire/screening/screening_weak.h @@ -0,0 +1,117 @@ +#pragma once + +#include "gridfire/screening/screening_abstract.h" +#include "gridfire/reaction/reaction.h" + +#include "fourdst/logging/logging.h" +#include "quill/Logger.h" +#include "quill/LogMacros.h" + +#include "cppad/cppad.hpp" + +namespace gridfire::screening { + class WeakScreeningModel final : public ScreeningModel { + public: + [[nodiscard]] std::vector calculateScreeningFactors( + const reaction::LogicalReactionSet& reactions, + const std::vector& species, + const std::vector& Y, + const double T9, + const double rho + ) const override; + + [[nodiscard]] std::vector> calculateScreeningFactors( + const reaction::LogicalReactionSet& reactions, + const std::vector& species, + const std::vector>& Y, + const CppAD::AD T9, + const CppAD::AD rho + ) const override; + private: + quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); + + private: + template + [[nodiscard]] std::vector calculateFactors_impl( + const reaction::LogicalReactionSet& reactions, + const std::vector& species, + const std::vector& Y, + const T T9, + const T rho + ) const; + }; + + template + std::vector WeakScreeningModel::calculateFactors_impl( + const reaction::LogicalReactionSet& reactions, + const std::vector& species, + const std::vector& Y, + const T T9, + const T rho + ) const { + LOG_TRACE_L1( + m_logger, + "Calculating weak screening factors for {} reactions...", + reactions.size() + ); + // --- CppAD Safe low temp checking --- + const T zero(0.0); + const T one(1.0); + const T low_temp_threshold(1e-9); + + const T low_T_flag = CppAD::CondExpLt(T9, low_temp_threshold, zero, one); + + // --- Calculate composition-dependent terms --- + // ζ = ∑(Z_i^2 + Z_i) * X_i / A_i + // This simplifies somewhat to ζ = ∑ (Z_i^2 + Z_i) * Y_i + // Where Y_i is the molar abundance (mol/g) + T zeta(0.0); + for (size_t i = 0; i < species.size(); ++i) { + const T Z = species[i].m_z; + zeta += (Z * Z + Z) * Y[i]; + } + + // --- Constant prefactors --- + const T T7 = T9 * static_cast(100.00); + const T T7_safe = CppAD::CondExpLe(T7, low_temp_threshold, low_temp_threshold, T7); + const T prefactor = static_cast(0.188) * CppAD::sqrt(rho / (T7_safe * T7_safe * T7_safe)) * CppAD::sqrt(zeta); + + // --- Loop through reactions and calculate screening factors for each --- + std::vector factors; + factors.reserve(reactions.size()); + for (const auto& reaction : reactions) { + T H_12(0.0); // screening abundance term + const auto& reactants = reaction.reactants(); + const bool isTripleAlpha = ( + reactants.size() == 3 && + reactants[0].m_z == 2 && + reactants[1].m_z == 2 && + reactants[2].m_z == 2 && + reactants[0] == reactants[1] && + reactants[1] == reactants[2] + ); + if (reactants.size() == 2) { + LOG_TRACE_L3(m_logger, "Calculating screening factor for reaction: {}", reaction.peName()); + const T Z1 = static_cast(reactants[0].m_z); + const T Z2 = static_cast(reactants[1].m_z); + H_12 = prefactor * Z1 * Z2; + } + else if (isTripleAlpha) { + LOG_TRACE_L3(m_logger, "Special case for triple alpha process in reaction: {}", reaction.peName()); + // Special case for triple alpha process + const T Z_alpha = static_cast(2.0); + const T H_alpha_alpha = prefactor * Z_alpha * Z_alpha; + H_12 = static_cast(3.0) * H_alpha_alpha; // Triple alpha process + } + // For 1 body reactions H_12 remains 0 so e^H_12 will be 1.0 (screening does not apply) + // Aside from triple alpha, all other astrophysically relevant reactions are 2-body in the weak screening regime + + H_12 *= low_T_flag; // Apply low temperature flag to screening factor + H_12 = CppAD::CondExpGe(H_12, static_cast(2.0), static_cast(2.0), H_12); // Caps the screening factor at 10 to avoid numerical issues + factors.push_back(CppAD::exp(H_12)); + // std::cout << "Screening factor: " << reaction.peName() << " : " << factors.back() << "(" << H_12 << ")" << std::endl; + } + return factors; + } + +} \ 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 86ff74bb..2f65ae60 100644 --- a/src/network/include/gridfire/solver/solver.h +++ b/src/network/include/gridfire/solver/solver.h @@ -2,7 +2,7 @@ #include "gridfire/engine/engine_graph.h" #include "gridfire/engine/engine_abstract.h" -#include "gridfire/engine/engine_adaptive.h" +#include "../engine/views/engine_adaptive.h" #include "gridfire/network.h" #include "fourdst/logging/logging.h" @@ -10,7 +10,7 @@ #include "quill/Logger.h" -#include "Eigen/Dense" +#include "unsupported/Eigen/NonLinearOptimization" // Required for LevenbergMarquardt #include @@ -95,13 +95,13 @@ namespace gridfire::solver { * @see AdaptiveEngineView * @see DynamicEngine::getSpeciesTimescales() */ - class QSENetworkSolver final : public AdaptiveNetworkSolverStrategy { + class QSENetworkSolver final : public DynamicNetworkSolverStrategy { public: /** * @brief Constructor for the QSENetworkSolver. * @param engine The adaptive engine view to use for evaluating the network. */ - using AdaptiveNetworkSolverStrategy::AdaptiveNetworkSolverStrategy; + using DynamicNetworkSolverStrategy::DynamicNetworkSolverStrategy; /** * @brief Evaluates the network for a given timestep using the QSE approach. @@ -310,9 +310,13 @@ namespace gridfire::solver { 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. - const std::vector& m_YDynamic; ///< Abundances of the dynamic species. + 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. @@ -321,7 +325,7 @@ namespace gridfire::solver { /** * @brief Constructor for the EigenFunctor. * @param engine The engine used to evaluate the network. - * @param YDynamic Abundances of the dynamic species. + * @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. @@ -329,34 +333,37 @@ namespace gridfire::solver { */ EigenFunctor( DynamicEngine& engine, - const std::vector& YDynamic, + const std::vector& YFull, const std::vector& dynamicSpeciesIndices, const std::vector& QSESpeciesIndices, const double T9, const double rho ) : m_engine(engine), - m_YDynamic(YDynamic), + m_YFull(YFull), m_dynamicSpeciesIndices(dynamicSpeciesIndices), m_QSESpeciesIndices(QSESpeciesIndices), m_T9(T9), m_rho(rho) {} + 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 Input vector of QSE species abundances (logarithmic). + * @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, OutputType& f_QSE) const; + int operator()(const InputType& v_QSE_log, OutputType& f_QSE) const; /** * @brief Calculates the Jacobian matrix for the QSE species. - * @param v_QSE Input vector of QSE species abundances (logarithmic). + * @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, JacobianType& J_QSE) const; + int df(const InputType& v_QSE_log, JacobianType& J_QSE) const; }; private: quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log"); ///< Logger instance. @@ -404,6 +411,7 @@ namespace gridfire::solver { 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. /** * @brief Constructor for the RHSFunctor. @@ -485,49 +493,45 @@ namespace gridfire::solver { }; template - int QSENetworkSolver::EigenFunctor::operator()(const InputType &v_QSE, OutputType &f_QSE) const { - std::vector YFull(m_engine.getNetworkSpecies().size(), 0.0); - Eigen::VectorXd Y_QSE(v_QSE.array().exp()); - for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) { - YFull[m_dynamicSpeciesIndices[i]] = m_YDynamic[i]; - } + int QSENetworkSolver::EigenFunctor::operator()(const InputType &v_QSE_log, OutputType &f_QSE) const { + std::vector y = m_YFull; // Full vector of species abundances + Eigen::VectorXd Y_QSE = v_QSE_log.array().exp(); + for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { - YFull[m_QSESpeciesIndices[i]] = Y_QSE(i); - } - const auto [full_dYdt, specificEnergyGenerationRate] = m_engine.calculateRHSAndEnergy(YFull, m_T9, m_rho); - f_QSE.resize(m_QSESpeciesIndices.size()); - for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { - f_QSE(i) = full_dYdt[m_QSESpeciesIndices[i]]; + y[m_QSESpeciesIndices[i]] = Y_QSE(i); } + + const auto [dydt, specificEnergyRate] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho); + 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 YFull(m_engine.getNetworkSpecies().size(), 0.0); - Eigen::VectorXd Y_QSE(v_QSE.array().exp()); - for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) { - YFull[m_dynamicSpeciesIndices[i]] = m_YDynamic[i]; - } + int QSENetworkSolver::EigenFunctor::df(const InputType& v_QSE_log, JacobianType& J_QSE) const { + std::vector y = m_YFull; + Eigen::VectorXd Y_QSE = v_QSE_log.array().exp(); + for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) { - YFull[m_QSESpeciesIndices[i]] = Y_QSE(i); + y[m_QSESpeciesIndices[i]] = Y_QSE(i); } - m_engine.generateJacobianMatrix(YFull, m_T9, m_rho); + m_engine.generateJacobianMatrix(y, m_T9, m_rho); - Eigen::MatrixXd J_orig(m_QSESpeciesIndices.size(), m_QSESpeciesIndices.size()); + 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_orig(i, j) = m_engine.getJacobianMatrixEntry(m_QSESpeciesIndices[i], m_QSESpeciesIndices[j]); + J_QSE(i, j) = m_engine.getJacobianMatrixEntry(m_QSESpeciesIndices[i], m_QSESpeciesIndices[j]); } } - J_QSE = J_orig; for (long j = 0; j < J_QSE.cols(); ++j) { - J_QSE.col(j) *= Y_QSE(j); // Chain rule for log space + J_QSE(j, j) *= Y_QSE(j); // chain rule for log space transformation } - return 0; // Success + return 0; } diff --git a/src/network/include/gridfire/utils/logging.h b/src/network/include/gridfire/utils/logging.h new file mode 100644 index 00000000..e3e6f96f --- /dev/null +++ b/src/network/include/gridfire/utils/logging.h @@ -0,0 +1,15 @@ +#pragma once + +#include "gridfire/engine/engine_abstract.h" + +#include +#include + +namespace gridfire::utils { + std::string formatNuclearTimescaleLogString( + const DynamicEngine& engine, + const std::vector& Y, + const double T9, + const double rho + ); +} \ No newline at end of file diff --git a/src/network/lib/engine/engine_approx8.cpp b/src/network/lib/engine/engine_approx8.cpp index a4b58169..56cbc50d 100644 --- a/src/network/lib/engine/engine_approx8.cpp +++ b/src/network/lib/engine/engine_approx8.cpp @@ -509,7 +509,6 @@ namespace gridfire::approx8{ vector_type Approx8Network::convert_netIn(const NetIn &netIn) { vector_type y(Approx8Net::nVar, 0.0); y[Approx8Net::ih1] = netIn.composition.getNumberFraction("H-1"); - std::cout << "Approx8::convert_netIn -> H-1 fraction: " << y[Approx8Net::ih1] << std::endl; y[Approx8Net::ihe3] = netIn.composition.getNumberFraction("He-3"); y[Approx8Net::ihe4] = netIn.composition.getNumberFraction("He-4"); y[Approx8Net::ic12] = netIn.composition.getNumberFraction("C-12"); diff --git a/src/network/lib/engine/engine_graph.cpp b/src/network/lib/engine/engine_graph.cpp index 92c01006..bbadaf34 100644 --- a/src/network/lib/engine/engine_graph.cpp +++ b/src/network/lib/engine/engine_graph.cpp @@ -1,6 +1,7 @@ #include "gridfire/engine/engine_graph.h" #include "gridfire/reaction/reaction.h" #include "gridfire/network.h" +#include "gridfire/screening/screening_types.h" #include "fourdst/composition/species.h" #include "fourdst/composition/atomicSpecies.h" @@ -262,6 +263,15 @@ namespace gridfire { return calculateAllDerivatives(Y_in, T9, rho); } + void GraphEngine::setScreeningModel(const screening::ScreeningType model) { + m_screeningModel = screening::selectScreeningModel(model); + m_screeningType = model; + } + + screening::ScreeningType GraphEngine::getScreeningModel() const { + return m_screeningType; + } + double GraphEngine::calculateMolarReactionFlow( const reaction::Reaction &reaction, const std::vector &Y, @@ -439,6 +449,10 @@ namespace gridfire { return speciesTimescales; } + void GraphEngine::update(const NetIn &netIn) { + return; // No-op for GraphEngine, as it does not support manually triggering updates. + } + void GraphEngine::recordADTape() { LOG_TRACE_L1(m_logger, "Recording AD tape for the RHS calculation..."); diff --git a/src/network/lib/engine/views/engine_adaptive.cpp b/src/network/lib/engine/views/engine_adaptive.cpp index 3fa54399..5a0763c5 100644 --- a/src/network/lib/engine/views/engine_adaptive.cpp +++ b/src/network/lib/engine/views/engine_adaptive.cpp @@ -213,6 +213,14 @@ namespace gridfire { } + void AdaptiveEngineView::setScreeningModel(const screening::ScreeningType model) { + m_baseEngine.setScreeningModel(model); + } + + screening::ScreeningType AdaptiveEngineView::getScreeningModel() const { + return m_baseEngine.getScreeningModel(); + } + std::vector AdaptiveEngineView::mapCulledToFull(const std::vector& culled) const { std::vector full(m_baseEngine.getNetworkSpecies().size(), 0.0); for (size_t i_culled = 0; i_culled < culled.size(); ++i_culled) { @@ -338,7 +346,7 @@ namespace gridfire { const double maxFlow ) const { LOG_TRACE_L1(m_logger, "Culling reactions based on flow rates..."); - const double relative_culling_threshold = m_config.get("gridfire:AdaptiveEngineView:RelativeCullingThreshold", 1e-75); + const auto relative_culling_threshold = m_config.get("gridfire:AdaptiveEngineView:RelativeCullingThreshold", 1e-75); double absoluteCullingThreshold = relative_culling_threshold * maxFlow; LOG_DEBUG(m_logger, "Relative culling threshold: {:0.3E} ({})", relative_culling_threshold, absoluteCullingThreshold); std::vector culledReactions; diff --git a/src/network/lib/engine/views/engine_defined.cpp b/src/network/lib/engine/views/engine_defined.cpp index 9520af78..c8fc781f 100644 --- a/src/network/lib/engine/views/engine_defined.cpp +++ b/src/network/lib/engine/views/engine_defined.cpp @@ -1,3 +1,310 @@ -// -// Created by Emily Boudreaux on 6/30/25. -// +#include "gridfire/engine/views/engine_defined.h" + +#include "quill/LogMacros.h" + +namespace gridfire { + using fourdst::atomic::Species; + + FileDefinedEngineView::FileDefinedEngineView( + DynamicEngine &baseEngine, + const std::string &fileName, + const io::NetworkFileParser &parser + ): + m_baseEngine(baseEngine), + m_fileName(fileName), + m_parser(parser), + m_activeSpecies(baseEngine.getNetworkSpecies()), + m_activeReactions(baseEngine.getNetworkReactions()) { + buildFromFile(fileName); + } + + const DynamicEngine & FileDefinedEngineView::getBaseEngine() const { + return m_baseEngine; + } + + const std::vector & FileDefinedEngineView::getNetworkSpecies() const { + return m_activeSpecies; + } + + StepDerivatives FileDefinedEngineView::calculateRHSAndEnergy( + const std::vector &Y_defined, + const double T9, + const double rho + ) const { + validateNetworkState(); + + const auto Y_full = mapViewToFull(Y_defined); + const auto [dydt, nuclearEnergyGenerationRate] = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho); + + StepDerivatives definedResults; + definedResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate; + definedResults.dydt = mapFullToView(dydt); + return definedResults; + } + + void FileDefinedEngineView::generateJacobianMatrix( + const std::vector &Y_defined, + const double T9, + const double rho + ) { + validateNetworkState(); + + const auto Y_full = mapViewToFull(Y_defined); + m_baseEngine.generateJacobianMatrix(Y_full, T9, rho); + } + + double FileDefinedEngineView::getJacobianMatrixEntry( + const int i_defined, + const int j_defined + ) const { + validateNetworkState(); + + const size_t i_full = mapViewToFullSpeciesIndex(i_defined); + const size_t j_full = mapViewToFullSpeciesIndex(j_defined); + + return m_baseEngine.getJacobianMatrixEntry(i_full, j_full); + } + + void FileDefinedEngineView::generateStoichiometryMatrix() { + validateNetworkState(); + + m_baseEngine.generateStoichiometryMatrix(); + } + + int FileDefinedEngineView::getStoichiometryMatrixEntry( + const int speciesIndex_defined, + const int reactionIndex_defined + ) const { + validateNetworkState(); + + const size_t i_full = mapViewToFullSpeciesIndex(speciesIndex_defined); + const size_t j_full = mapViewToFullReactionIndex(reactionIndex_defined); + return m_baseEngine.getStoichiometryMatrixEntry(i_full, j_full); + } + + double FileDefinedEngineView::calculateMolarReactionFlow( + const reaction::Reaction &reaction, + const std::vector &Y_defined, + const double T9, + const double rho + ) const { + validateNetworkState(); + + if (!m_activeReactions.contains(reaction)) { + LOG_ERROR(m_logger, "Reaction '{}' is not part of the active reactions in the file defined engine view.", reaction.id()); + m_logger -> flush_log(); + throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id())); + } + const auto Y_full = mapViewToFull(Y_defined); + return m_baseEngine.calculateMolarReactionFlow(reaction, Y_full, T9, rho); + } + + const reaction::LogicalReactionSet & FileDefinedEngineView::getNetworkReactions() const { + validateNetworkState(); + + return m_activeReactions; + } + + std::unordered_map FileDefinedEngineView::getSpeciesTimescales( + const std::vector &Y_defined, + const double T9, + const double rho + ) const { + validateNetworkState(); + + const auto Y_full = mapViewToFull(Y_defined); + const auto fullTimescales = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho); + + std::unordered_map definedTimescales; + for (const auto& active_species : m_activeSpecies) { + if (fullTimescales.contains(active_species)) { + definedTimescales[active_species] = fullTimescales.at(active_species); + } + } + return definedTimescales; + } + + void FileDefinedEngineView::update(const NetIn &netIn) { + if (m_isStale) { + buildFromFile(m_fileName); + } + } + + void FileDefinedEngineView::setNetworkFile(const std::string &fileName) { + m_isStale = true; + m_fileName = fileName; + LOG_DEBUG(m_logger, "File '{}' set to '{}'. FileDefinedNetworkView is now stale! You MUST call update() before you use it!", m_fileName, fileName); + } + + void FileDefinedEngineView::setScreeningModel(const screening::ScreeningType model) { + m_baseEngine.setScreeningModel(model); + } + + screening::ScreeningType FileDefinedEngineView::getScreeningModel() const { + return m_baseEngine.getScreeningModel(); + } + + std::vector FileDefinedEngineView::constructSpeciesIndexMap() const { + LOG_TRACE_L1(m_logger, "Constructing species index map for file defined engine view..."); + std::unordered_map fullSpeciesReverseMap; + const auto& fullSpeciesList = m_baseEngine.getNetworkSpecies(); + + fullSpeciesReverseMap.reserve(fullSpeciesList.size()); + + for (size_t i = 0; i < fullSpeciesList.size(); ++i) { + fullSpeciesReverseMap[fullSpeciesList[i]] = i; + } + + std::vector speciesIndexMap; + speciesIndexMap.reserve(m_activeSpecies.size()); + + for (const auto& active_species : m_activeSpecies) { + auto it = fullSpeciesReverseMap.find(active_species); + if (it != fullSpeciesReverseMap.end()) { + speciesIndexMap.push_back(it->second); + } else { + LOG_ERROR(m_logger, "Species '{}' not found in full species map.", active_species.name()); + m_logger -> flush_log(); + throw std::runtime_error("Species not found in full species map: " + std::string(active_species.name())); + } + } + LOG_TRACE_L1(m_logger, "Species index map constructed with {} entries.", speciesIndexMap.size()); + return speciesIndexMap; + + } + + std::vector FileDefinedEngineView::constructReactionIndexMap() const { + LOG_TRACE_L1(m_logger, "Constructing reaction index map for file defined engine view..."); + + // --- Step 1: Create a reverse map using the reaction's unique ID as the key. --- + std::unordered_map fullReactionReverseMap; + const auto& fullReactionSet = m_baseEngine.getNetworkReactions(); + fullReactionReverseMap.reserve(fullReactionSet.size()); + + for (size_t i_full = 0; i_full < fullReactionSet.size(); ++i_full) { + fullReactionReverseMap[fullReactionSet[i_full].id()] = i_full; + } + + // --- Step 2: Build the final index map using the active reaction set. --- + std::vector reactionIndexMap; + reactionIndexMap.reserve(m_activeReactions.size()); + + for (const auto& active_reaction_ptr : m_activeReactions) { + auto it = fullReactionReverseMap.find(active_reaction_ptr.id()); + + if (it != fullReactionReverseMap.end()) { + reactionIndexMap.push_back(it->second); + } else { + LOG_ERROR(m_logger, "Active reaction '{}' not found in base engine during reaction index map construction.", active_reaction_ptr.id()); + m_logger->flush_log(); + throw std::runtime_error("Mismatch between active reactions and base engine."); + } + } + + LOG_TRACE_L1(m_logger, "Reaction index map constructed with {} entries.", reactionIndexMap.size()); + return reactionIndexMap; + } + + void FileDefinedEngineView::buildFromFile(const std::string &fileName) { + LOG_TRACE_L1(m_logger, "Building file defined engine view from {}...", fileName); + auto [reactionPENames] = m_parser.parse(fileName); + + m_activeReactions.clear(); + m_activeSpecies.clear(); + + std::unordered_set seenSpecies; + + const auto& fullNetworkReactionSet = m_baseEngine.getNetworkReactions(); + for (const auto& peName : reactionPENames) { + if (!fullNetworkReactionSet.contains(peName)) { + LOG_ERROR(m_logger, "Reaction with name '{}' not found in the base engine's network reactions. Aborting...", peName); + m_logger->flush_log(); + throw std::runtime_error("Reaction with name '" + std::string(peName) + "' not found in the base engine's network reactions."); + } + auto reaction = fullNetworkReactionSet[peName]; + for (const auto& reactant : reaction.reactants()) { + if (!seenSpecies.contains(reactant)) { + seenSpecies.insert(reactant); + m_activeSpecies.push_back(reactant); + } + } + for (const auto& product : reaction.products()) { + if (!seenSpecies.contains(product)) { + seenSpecies.insert(product); + m_activeSpecies.push_back(product); + } + } + m_activeReactions.add_reaction(reaction); + } + LOG_TRACE_L1(m_logger, "File defined engine view built with {} active species and {} active reactions.", m_activeSpecies.size(), m_activeReactions.size()); + LOG_DEBUG(m_logger, "Active species: {}", [this]() -> std::string { + std::string result; + for (const auto& species : m_activeSpecies) { + result += std::string(species.name()) + ", "; + } + if (!result.empty()) { + result.pop_back(); // Remove last space + result.pop_back(); // Remove last comma + } + return result; + }()); + LOG_DEBUG(m_logger, "Active reactions: {}", [this]() -> std::string { + std::string result; + for (const auto& reaction : m_activeReactions) { + result += std::string(reaction.id()) + ", "; + } + if (!result.empty()) { + result.pop_back(); // Remove last space + result.pop_back(); // Remove last comma + } + return result; + }()); + m_speciesIndexMap = constructSpeciesIndexMap(); + m_reactionIndexMap = constructReactionIndexMap(); + m_isStale = false; + } + + std::vector FileDefinedEngineView::mapViewToFull(const std::vector& culled) const { + std::vector full(m_baseEngine.getNetworkSpecies().size(), 0.0); + for (size_t i_culled = 0; i_culled < culled.size(); ++i_culled) { + const size_t i_full = m_speciesIndexMap[i_culled]; + full[i_full] += culled[i_culled]; + } + return full; + } + + std::vector FileDefinedEngineView::mapFullToView(const std::vector& full) const { + std::vector culled(m_activeSpecies.size(), 0.0); + for (size_t i_culled = 0; i_culled < m_activeSpecies.size(); ++i_culled) { + const size_t i_full = m_speciesIndexMap[i_culled]; + culled[i_culled] = full[i_full]; + } + return culled; + } + + size_t FileDefinedEngineView::mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const { + if (culledSpeciesIndex < 0 || culledSpeciesIndex >= static_cast(m_speciesIndexMap.size())) { + LOG_ERROR(m_logger, "Defined index {} is out of bounds for species index map of size {}.", culledSpeciesIndex, m_speciesIndexMap.size()); + m_logger->flush_log(); + throw std::out_of_range("Defined index " + std::to_string(culledSpeciesIndex) + " is out of bounds for species index map of size " + std::to_string(m_speciesIndexMap.size()) + "."); + } + return m_speciesIndexMap[culledSpeciesIndex]; + } + + size_t FileDefinedEngineView::mapViewToFullReactionIndex(size_t culledReactionIndex) const { + if (culledReactionIndex < 0 || culledReactionIndex >= static_cast(m_reactionIndexMap.size())) { + LOG_ERROR(m_logger, "Defined index {} is out of bounds for reaction index map of size {}.", culledReactionIndex, m_reactionIndexMap.size()); + m_logger->flush_log(); + throw std::out_of_range("Defined index " + std::to_string(culledReactionIndex) + " is out of bounds for reaction index map of size " + std::to_string(m_reactionIndexMap.size()) + "."); + } + return m_reactionIndexMap[culledReactionIndex]; + } + + void FileDefinedEngineView::validateNetworkState() const { + if (m_isStale) { + LOG_ERROR(m_logger, "File defined engine view is stale. Please call update() with a valid NetIn object."); + m_logger->flush_log(); + throw std::runtime_error("File defined engine view is stale. Please call update() with a valid NetIn object."); + } + } +} diff --git a/src/network/lib/io/network_file.cpp b/src/network/lib/io/network_file.cpp index 9520af78..273a6b79 100644 --- a/src/network/lib/io/network_file.cpp +++ b/src/network/lib/io/network_file.cpp @@ -1,3 +1,77 @@ -// -// Created by Emily Boudreaux on 6/30/25. -// +#include "gridfire/io/network_file.h" + +#include +#include +#include +#include +#include + +#include "quill/LogMacros.h" + +namespace gridfire::io { + namespace { + inline void ltrim(std::string &s) { + s.erase( + s.begin(), + std::ranges::find_if(s, + [](const unsigned char ch) { + return !std::isspace(ch); + }) + ); + } + + inline void rtrim(std::string &s) { + s.erase( + std::find_if( + s.rbegin(), + s.rend(), + [](const unsigned char ch) { + return !std::isspace(ch); + }).base(), + s.end() + ); + } + + inline void trim(std::string &s) { + ltrim(s); + rtrim(s); + } + + + } + SimpleReactionListFileParser::SimpleReactionListFileParser() {} + + ParsedNetworkData SimpleReactionListFileParser::parse(const std::string& filename) const { + LOG_TRACE_L1(m_logger, "Parsing simple reaction list file: {}", filename); + + std::ifstream file(filename); + if (!file.is_open()) { + LOG_ERROR(m_logger, "Failed to open file: {}", filename); + m_logger -> flush_log(); + throw std::runtime_error("Could not open file: " + filename); + } + + ParsedNetworkData parsed; + std::string line; + int line_number = 0; + while (std::getline(file, line)) { + line_number++; + LOG_TRACE_L3(m_logger, "Parsing reaction list file {}, line {}: {}", filename, line_number, line); + + const size_t comment_pos = line.find('#'); + if (comment_pos != std::string::npos) { + line = line.substr(0, comment_pos); + } + + trim(line); + + if (line.empty()) { + continue; // Skip empty lines + } + parsed.reactionPENames.push_back(line); + } + LOG_TRACE_L1(m_logger, "Parsed {} reactions from file: {}", parsed.reactionPENames.size(), filename); + return parsed; + } + +} \ No newline at end of file diff --git a/src/network/lib/network.cpp b/src/network/lib/network.cpp index c6dbaabd..1c161069 100644 --- a/src/network/lib/network.cpp +++ b/src/network/lib/network.cpp @@ -84,7 +84,7 @@ namespace gridfire { } } const ReactionSet reactionSet(reaclibReactions); - return LogicalReactionSet(reactionSet); + return packReactionSetToLogicalReactionSet(reactionSet); } // Trim whitespace from both ends of a string diff --git a/src/network/lib/reaction/reaclib.cpp b/src/network/lib/reaction/reaclib.cpp index f29e11b8..87e54368 100644 --- a/src/network/lib/reaction/reaclib.cpp +++ b/src/network/lib/reaction/reaclib.cpp @@ -122,10 +122,12 @@ namespace gridfire::reaclib { } // The ReactionSet takes the vector of all individual reactions. - reaction::ReactionSet reaction_set(std::move(reaction_list)); + const reaction::ReactionSet reaction_set(std::move(reaction_list)); // The LogicalReactionSet groups reactions by their peName, which is what we want. - s_all_reaclib_reactions_ptr = new reaction::LogicalReactionSet(reaction_set); + s_all_reaclib_reactions_ptr = new reaction::LogicalReactionSet( + reaction::packReactionSetToLogicalReactionSet(reaction_set) + ); s_initialized = true; } diff --git a/src/network/lib/reaction/reaction.cpp b/src/network/lib/reaction/reaction.cpp index b72c4d62..9cb03ac8 100644 --- a/src/network/lib/reaction/reaction.cpp +++ b/src/network/lib/reaction/reaction.cpp @@ -137,152 +137,6 @@ namespace gridfire::reaction { return XXHash64::hash(m_id.data(), m_id.size(), seed); } - ReactionSet::ReactionSet( - std::vector reactions - ) : - m_reactions(std::move(reactions)) { - if (m_reactions.empty()) { - return; // Case where the reactions will be added later. - } - m_reactionNameMap.reserve(reactions.size()); - for (const auto& reaction : m_reactions) { - m_id += reaction.id(); - m_reactionNameMap.emplace(reaction.id(), reaction); - } - } - - ReactionSet::ReactionSet(const ReactionSet &other) { - m_reactions.reserve(other.m_reactions.size()); - for (const auto& reaction_ptr: other.m_reactions) { - m_reactions.push_back(reaction_ptr); - } - - m_reactionNameMap.reserve(other.m_reactionNameMap.size()); - for (const auto& reaction_ptr : m_reactions) { - m_reactionNameMap.emplace(reaction_ptr.id(), reaction_ptr); - } - } - - ReactionSet & ReactionSet::operator=(const ReactionSet &other) { - if (this != &other) { - ReactionSet temp(other); - std::swap(m_reactions, temp.m_reactions); - std::swap(m_reactionNameMap, temp.m_reactionNameMap); - } - return *this; - } - - void ReactionSet::add_reaction(Reaction reaction) { - m_reactions.emplace_back(reaction); - m_id += m_reactions.back().id(); - m_reactionNameMap.emplace(m_reactions.back().id(), m_reactions.back()); - } - - void ReactionSet::remove_reaction(const Reaction& reaction) { - if (!m_reactionNameMap.contains(std::string(reaction.id()))) { - return; - } - - m_reactionNameMap.erase(std::string(reaction.id())); - - std::erase_if(m_reactions, [&reaction](const Reaction& r) { - return r == reaction; - }); - } - - bool ReactionSet::contains(const std::string_view& id) const { - for (const auto& reaction : m_reactions) { - if (reaction.id() == id) { - return true; - } - } - return false; - } - - bool ReactionSet::contains(const Reaction& reaction) const { - for (const auto& r : m_reactions) { - if (r == reaction) { - return true; - } - } - return false; - } - - void ReactionSet::clear() { - m_reactions.clear(); - m_reactionNameMap.clear(); - } - - bool ReactionSet::contains_species(const Species& species) const { - for (const auto& reaction : m_reactions) { - if (reaction.contains(species)) { - return true; - } - } - return false; - } - - bool ReactionSet::contains_reactant(const Species& species) const { - for (const auto& r : m_reactions) { - if (r.contains_reactant(species)) { - return true; - } - } - return false; - } - - bool ReactionSet::contains_product(const Species& species) const { - for (const auto& r : m_reactions) { - if (r.contains_product(species)) { - return true; - } - } - return false; - } - - const Reaction& ReactionSet::operator[](const size_t index) const { - if (index >= m_reactions.size()) { - m_logger -> flush_log(); - throw std::out_of_range("Index" + std::to_string(index) + " out of range for ReactionSet of size " + std::to_string(m_reactions.size()) + "."); - } - return m_reactions[index]; - } - - const Reaction& ReactionSet::operator[](const std::string_view& id) const { - if (auto it = m_reactionNameMap.find(std::string(id)); it != m_reactionNameMap.end()) { - return it->second; - } - m_logger -> flush_log(); - throw std::out_of_range("Species " + std::string(id) + " does not exist in ReactionSet."); - } - - bool ReactionSet::operator==(const ReactionSet& other) const { - if (size() != other.size()) { - return false; - } - return hash() == other.hash(); - } - - bool ReactionSet::operator!=(const ReactionSet& other) const { - return !(*this == other); - } - - uint64_t ReactionSet::hash(uint64_t seed) const { - if (m_reactions.empty()) { - return XXHash64::hash(nullptr, 0, seed); - } - std::vector individualReactionHashes; - individualReactionHashes.reserve(m_reactions.size()); - for (const auto& reaction : m_reactions) { - individualReactionHashes.push_back(reaction.hash(seed)); - } - - std::ranges::sort(individualReactionHashes); - - const void* data = static_cast(individualReactionHashes.data()); - size_t sizeInBytes = individualReactionHashes.size() * sizeof(uint64_t); - return XXHash64::hash(data, sizeInBytes, seed); - } LogicalReaction::LogicalReaction(const std::vector& reactants) : @@ -344,21 +198,21 @@ namespace gridfire::reaction { return calculate_rate>(T9); } - LogicalReactionSet::LogicalReactionSet(const ReactionSet &reactionSet) : - ReactionSet(std::vector()) { + LogicalReactionSet packReactionSetToLogicalReactionSet(const ReactionSet& reactionSet) { + std::unordered_map> groupedReactions; - std::unordered_map> grouped_reactions; + for (const auto& reaction: reactionSet) { + groupedReactions[reaction.peName()].push_back(reaction); + } - for (const auto& reaction : reactionSet) { - grouped_reactions[reaction.peName()].push_back(reaction); - } - m_reactions.reserve(grouped_reactions.size()); - m_reactionNameMap.reserve(grouped_reactions.size()); - for (const auto &reactions_for_peName: grouped_reactions | std::views::values) { - LogicalReaction logical_reaction(reactions_for_peName); - m_reactionNameMap.emplace(logical_reaction.id(), logical_reaction); - m_reactions.push_back(std::move(logical_reaction)); + std::vector reactions; + reactions.reserve(groupedReactions.size()); + + for (const auto &reactionsGroup: groupedReactions | std::views::values) { + LogicalReaction logicalReaction(reactionsGroup); + reactions.push_back(logicalReaction); } + return LogicalReactionSet(std::move(reactions)); } } @@ -376,4 +230,11 @@ namespace std { return s.hash(0); } }; + + template<> + struct hash { + size_t operator()(const gridfire::reaction::LogicalReactionSet& s) const noexcept { + return s.hash(0); + } + }; } // namespace std diff --git a/src/network/lib/screening/screening_bare.cpp b/src/network/lib/screening/screening_bare.cpp new file mode 100644 index 00000000..e817ea40 --- /dev/null +++ b/src/network/lib/screening/screening_bare.cpp @@ -0,0 +1,31 @@ +#include "gridfire/screening/screening_bare.h" + +#include "fourdst/composition/atomicSpecies.h" + +#include "cppad/cppad.hpp" + +#include + + +namespace gridfire::screening { + using ADDouble = CppAD::AD; + std::vector BareScreeningModel::calculateScreeningFactors( + const reaction::LogicalReactionSet &reactions, + const std::vector& species, + const std::vector &Y, + const ADDouble T9, + const ADDouble rho + ) const { + return calculateFactors_impl(reactions, species, Y, T9, rho); + } + + std::vector BareScreeningModel::calculateScreeningFactors( + const reaction::LogicalReactionSet &reactions, + const std::vector& species, + const std::vector &Y, + const double T9, + const double rho + ) const { + return calculateFactors_impl(reactions, species, Y, T9, rho); + } +} diff --git a/src/network/lib/screening/screening_types.cpp b/src/network/lib/screening/screening_types.cpp new file mode 100644 index 00000000..4d366b81 --- /dev/null +++ b/src/network/lib/screening/screening_types.cpp @@ -0,0 +1,19 @@ +#include "gridfire/screening/screening_abstract.h" +#include "gridfire/screening/screening_types.h" +#include "gridfire/screening/screening_weak.h" +#include "gridfire/screening/screening_bare.h" + +#include + +namespace gridfire::screening { + std::unique_ptr selectScreeningModel(const ScreeningType type) { + switch (type) { + case ScreeningType::WEAK: + return std::make_unique(); + case ScreeningType::BARE: + return std::make_unique(); + default: + return std::make_unique(); + } + } +} diff --git a/src/network/lib/screening/screening_weak.cpp b/src/network/lib/screening/screening_weak.cpp new file mode 100644 index 00000000..e617cdbf --- /dev/null +++ b/src/network/lib/screening/screening_weak.cpp @@ -0,0 +1,31 @@ +#include "gridfire/screening/screening_weak.h" + +#include "fourdst/composition/atomicSpecies.h" + +#include "cppad/cppad.hpp" + +#include + + +namespace gridfire::screening { + using ADDouble = CppAD::AD; + std::vector WeakScreeningModel::calculateScreeningFactors( + const reaction::LogicalReactionSet &reactions, + const std::vector& species, + const std::vector &Y, + const ADDouble T9, + const ADDouble rho + ) const { + return calculateFactors_impl(reactions, species, Y, T9, rho); + } + + std::vector WeakScreeningModel::calculateScreeningFactors( + const reaction::LogicalReactionSet &reactions, + const std::vector& species, + const std::vector &Y, + const double T9, + const double rho + ) const { + return calculateFactors_impl(reactions, species, Y, T9, rho); + } +} diff --git a/src/network/lib/solver/solver.cpp b/src/network/lib/solver/solver.cpp index 970e5188..34021816 100644 --- a/src/network/lib/solver/solver.cpp +++ b/src/network/lib/solver/solver.cpp @@ -2,6 +2,8 @@ #include "gridfire/engine/engine_graph.h" #include "gridfire/network.h" +#include "gridfire/utils/logging.h" + #include "fourdst/composition/atomicSpecies.h" #include "fourdst/composition/composition.h" #include "fourdst/config/config.h" @@ -15,6 +17,7 @@ #include #include #include +#include #include "quill/LogMacros.h" @@ -59,6 +62,21 @@ namespace gridfire::solver { 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(); @@ -185,18 +203,62 @@ namespace gridfire::solver { } Eigen::VectorXd QSENetworkSolver::calculateSteadyStateAbundances( - const std::vector &Y, - const double T9, - const double rho, - const dynamicQSESpeciesIndices &indices + 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_WARNING(m_logger, "QSE solver logic not yet implemented, assuming all QSE species have 0 abundance."); - // --- Prepare the QSE species vector --- - Eigen::VectorXd v_qse(indices.QSESpeciesIndices.size()); - v_qse.setZero(); - return v_qse.array(); + if (indices.QSESpeciesIndices.empty()) { + LOG_DEBUG(m_logger, "No QSE species to solve for."); + return Eigen::VectorXd(0); + } + // Use the EigenFunctor with Eigen's nonlinear solver + EigenFunctor functor( + m_engine, + Y, + indices.dynamicSpeciesIndices, + indices.QSESpeciesIndices, + T9, + rho + ); + + Eigen::VectorXd v_qse_log_initial(indices.QSESpeciesIndices.size()); + for (size_t i = 0; i < indices.QSESpeciesIndices.size(); ++i) { + v_qse_log_initial(i) = std::log(std::max(Y[indices.QSESpeciesIndices[i]], 1e-99)); + } + + const 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); + const Eigen::LevenbergMarquardtSpace::Status info = lm.minimize(v_qse_log_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 v_qse_log_initial.array().exp(); + } NetOut QSENetworkSolver::initializeNetworkWithShortIgnition(const NetIn &netIn) const { @@ -232,9 +294,14 @@ namespace gridfire::solver { 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); + m_engine.setScreeningModel(prevScreeningModel); + LOG_DEBUG(m_logger, "Restoring previous screening model: {}", static_cast(prevScreeningModel)); return postIgnition; } @@ -360,7 +427,8 @@ namespace gridfire::solver { speciesNames.push_back(std::string(species.name())); } - Composition outputComposition(speciesNames, finalMassFractions); + Composition outputComposition(speciesNames); + outputComposition.setMassFraction(speciesNames, finalMassFractions); outputComposition.finalize(true); NetOut netOut; @@ -377,6 +445,15 @@ namespace gridfire::solver { double t ) const { const std::vector y(Y.begin(), m_numSpecies + Y.begin()); + + // std::string timescales = utils::formatNuclearTimescaleLogString( + // m_engine, + // y, + // m_T9, + // m_rho + // ); + // LOG_TRACE_L2(m_logger, "{}", timescales); + auto [dydt, eps] = m_engine.calculateRHSAndEnergy(y, m_T9, m_rho); dYdt.resize(m_numSpecies + 1); std::ranges::copy(dydt, dYdt.begin()); diff --git a/src/network/lib/utils/logging.cpp b/src/network/lib/utils/logging.cpp new file mode 100644 index 00000000..19b618da --- /dev/null +++ b/src/network/lib/utils/logging.cpp @@ -0,0 +1,60 @@ +#include "gridfire/utils/logging.h" +#include "gridfire/engine/engine_abstract.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +std::string gridfire::utils::formatNuclearTimescaleLogString( + const DynamicEngine& engine, + std::vector const& Y, + const double T9, + const double rho +) { + auto const& timescales = engine.getSpeciesTimescales(Y, T9, rho); + + // Figure out how wide the "Species" column needs to be: + std::size_t maxNameLen = std::string_view("Species").size(); + for (const auto &key: timescales | std::views::keys) { + std::string_view name = key.name(); + maxNameLen = std::max(maxNameLen, name.size()); + } + + // Pick a fixed width for the timescale column: + constexpr int timescaleWidth = 12; + + std::ostringstream ss; + ss << "== Timescales (s) ==\n"; + + // Header row + ss << std::left << std::setw(static_cast(maxNameLen) + 2) << "Species" + << std::right << std::setw(timescaleWidth) << "Timescale (s)" << "\n"; + + // Underline + ss << std::string(static_cast(maxNameLen) + 2 + timescaleWidth, '=') << "\n"; + + ss << std::scientific; + + // Data rows + for (auto const& [species, timescale] : timescales) { + const std::string_view name = species.name(); + ss << std::left << std::setw(static_cast(maxNameLen) + 2) << name; + + if (std::isinf(timescale)) { + ss << std::right << std::setw(timescaleWidth) << "inf" << "\n"; + } else { + ss << std::right << std::setw(timescaleWidth) + << std::scientific << std::setprecision(3) << timescale << "\n"; + } + } + + // Footer underline + ss << std::string(static_cast(maxNameLen) + 2 + timescaleWidth, '=') << "\n"; + + return ss.str(); +} diff --git a/src/network/meson.build b/src/network/meson.build index 91a3e9e3..3c54b543 100644 --- a/src/network/meson.build +++ b/src/network/meson.build @@ -3,10 +3,16 @@ network_sources = files( 'lib/network.cpp', 'lib/engine/engine_approx8.cpp', 'lib/engine/engine_graph.cpp', - 'lib/engine/engine_adaptive.cpp', + 'lib/engine/views/engine_adaptive.cpp', + 'lib/engine/views/engine_defined.cpp', 'lib/reaction/reaction.cpp', 'lib/reaction/reaclib.cpp', + 'lib/io/network_file.cpp', 'lib/solver/solver.cpp', + 'lib/screening/screening_types.cpp', + 'lib/screening/screening_weak.cpp', + 'lib/screening/screening_bare.cpp', + 'lib/utils/logging.cpp', ) @@ -39,12 +45,19 @@ network_dep = declare_dependency( network_headers = files( 'include/gridfire/network.h', 'include/gridfire/engine/engine_abstract.h', - 'include/gridfire/engine/engine_view_abstract.h', + 'include/gridfire/engine/views/engine_view_abstract.h', 'include/gridfire/engine/engine_approx8.h', 'include/gridfire/engine/engine_graph.h', - 'include/gridfire/engine/engine_adaptive.h', + 'include/gridfire/engine/views/engine_adaptive.h', + 'include/gridfire/engine/views/engine_defined.h', 'include/gridfire/reaction/reaction.h', 'include/gridfire/reaction/reaclib.h', + 'include/gridfire/io/network_file.h', 'include/gridfire/solver/solver.h', + 'include/gridfire/screening/screening_abstract.h', + 'include/gridfire/screening/screening_bare.h', + 'include/gridfire/screening/screening_weak.h', + 'include/gridfire/screening/screening_types.h', + 'include/gridfire/utils/logging.h', ) install_headers(network_headers, subdir : 'gridfire')