diff --git a/src/network/private/netgraph.cpp b/src/network/private/netgraph.cpp index a855dfdb..9104da04 100644 --- a/src/network/private/netgraph.cpp +++ b/src/network/private/netgraph.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include @@ -212,6 +213,12 @@ namespace serif::network { // Check if the requested network has already been cached. // PERF: Rebuilding this should be pretty fast but it may be a good point of optimization in the future. const reaclib::REACLIBReactionSet validationReactionSet = build_reaclib_nuclear_network(composition, culling, T9); + // TODO: need some more robust method here to + // A. Build the basic network from the composition's species with non zero mass fractions. + // B. rebuild a new composition from the reaction set's reactants + products (with the mass fractions from the things that are only products set to 0) + // C. Rebuild the reaction set from the new composition + // D. Cull reactions where all reactants have mass fractions below the culling threshold. + // E. Be careful about maintaining caching through all of this // This allows for dynamic network modification while retaining caching for networks which are very similar. @@ -328,6 +335,53 @@ namespace serif::network { } } + void GraphNetwork::exportToDot(const std::string &filename) const { + LOG_INFO(m_logger, "Exporting network graph to DOT file: {}", filename); + + std::ofstream dotFile(filename); + if (!dotFile.is_open()) { + LOG_ERROR(m_logger, "Failed to open file for writing: {}", filename); + throw std::runtime_error("Failed to open file for writing: " + filename); + } + + dotFile << "digraph NuclearReactionNetwork {\n"; + dotFile << " graph [rankdir=LR, splines=true, overlap=false, bgcolor=\"#f0f0f0\"];\n"; + dotFile << " node [shape=circle, style=filled, fillcolor=\"#a7c7e7\", fontname=\"Helvetica\"];\n"; + dotFile << " edge [fontname=\"Helvetica\", fontsize=10];\n\n"; + + // 1. Define all species as nodes + dotFile << " // --- Species Nodes ---\n"; + for (const auto& species : m_networkSpecies) { + dotFile << " \"" << species.name() << "\" [label=\"" << species.name() << "\"];\n"; + } + dotFile << "\n"; + + // 2. Define all reactions as intermediate nodes and connect them + dotFile << " // --- Reaction Edges ---\n"; + for (const auto& reaction : m_reactions) { + // Create a unique ID for the reaction node + std::string reactionNodeId = "reaction_" + std::string(reaction.id()); + + // Define the reaction node (small, black dot) + dotFile << " \"" << reactionNodeId << "\" [shape=point, fillcolor=black, width=0.1, height=0.1, label=\"\"];\n"; + + // Draw edges from reactants to the reaction node + for (const auto& reactant : reaction.reactants()) { + dotFile << " \"" << reactant.name() << "\" -> \"" << reactionNodeId << "\";\n"; + } + + // Draw edges from the reaction node to products + for (const auto& product : reaction.products()) { + dotFile << " \"" << reactionNodeId << "\" -> \"" << product.name() << "\" [label=\"" << reaction.qValue() << " MeV\"];\n"; + } + dotFile << "\n"; + } + + dotFile << "}\n"; + dotFile.close(); + LOG_INFO(m_logger, "Successfully exported network to {}", filename); + } + NetOut GraphNetwork::evaluate(const NetIn &netIn) { namespace ublas = boost::numeric::ublas; namespace odeint = boost::numeric::odeint; diff --git a/src/network/public/netgraph.h b/src/network/public/netgraph.h index 88252cf7..76240f72 100644 --- a/src/network/public/netgraph.h +++ b/src/network/public/netgraph.h @@ -15,81 +15,252 @@ #include "quill/LogMacros.h" - // PERF: The function getNetReactionStoichiometry returns a map of species to their stoichiometric coefficients for a given reaction. // this makes extra copies of the species, which is not ideal and could be optimized further. // Even more relevant is the member m_reactionIDMap which makes copies of a REACLIBReaction for each reaction ID. // REACLIBReactions are quite large data structures, so this could be a performance bottleneck. +/** + * @file netgraph.h + * @brief Defines the GraphNetwork class for representing and evolving nuclear reaction networks as a heavily connected graph. + * + * This file provides the GraphNetwork class, which implements a graph-based nuclear reaction network + * using REACLIB reactions and supports both stiff and non-stiff ODE integration. The network is constructed + * from a composition and can be queried for species, reactions, and stoichiometry. + * + * This is a general nuclear reaction network; however, it does not currently manage reverse reactions, weak reactions, + * or other reactions which become much more relevant for more extreme astrophysical sources. + * + * @see serif::network::GraphNetwork + */ + namespace serif::network { - // Define a concept to check if a type is one of our allowed scalar types + + /** + * @brief Concept to check if a type is either double or CppAD::AD. + * + * Used to enable templated functions for both arithmetic and automatic differentiation types. + */ template concept IsArithmeticOrAD = std::is_same_v || std::is_same_v>; + /// Minimum density threshold (g/cm^3) below which reactions are ignored. static constexpr double MIN_DENSITY_THRESHOLD = 1e-18; + /// Minimum abundance threshold below which reactions are ignored. static constexpr double MIN_ABUNDANCE_THRESHOLD = 1e-18; + /// Minimum Jacobian entry threshold for sparsity. static constexpr double MIN_JACOBIAN_THRESHOLD = 1e-24; + /** + * @brief Graph-based nuclear reaction network using REACLIB reactions. + * + * GraphNetwork constructs a reaction network from a given composition, builds the associated + * stoichiometry and Jacobian matrices, and provides methods for evaluating the network's evolution + * using ODE solvers. It supports both stiff and non-stiff integration and can be queried for + * network species, reactions, and stoichiometry. + * + * @note Currently does not handle reverse reactions, weak reactions, or other complex reaction types. These will be added in future versions. + * + * Example usage: + * @code + * serif::composition::Composition comp = ...; + * serif::network::GraphNetwork net(comp); + * serif::network::NetIn input; + * input.composition = comp; + * input.temperature = 1.5e9; + * input.density = 1e6; + * input.tMax = 1.0; + * input.dt0 = 1e-3; + * serif::network::NetOut result = net.evaluate(input); + * @endcode + */ class GraphNetwork final : public Network { public: + /** + * @brief Construct a GraphNetwork from a composition. + * @param composition The composition specifying the initial isotopic abundances. + */ explicit GraphNetwork(const serif::composition::Composition &composition); + + /** + * @brief Construct a GraphNetwork from a composition with culling and temperature. + * @param composition The composition specifying the initial isotopic abundances. + * @param cullingThreshold specific reaction rate threshold to cull reactions below. + * @param T9 Temperature in units of 10^9 K where culling is evaluated at. + * + * @see serif::network::build_reaclib_nuclear_network + */ explicit GraphNetwork(const serif::composition::Composition &composition, const double cullingThreshold, const double T9); + /** + * @brief Evolve the network for the given input conditions. + * + * This is the primary entry point for users of GridFire. This function integrates the network ODEs using either a stiff or non-stiff solver, + * depending on the detected stiffness of the system. + * + * @param netIn The input structure containing composition, temperature, density, and integration parameters. + * @return NetOut The output structure containing the final composition, number of steps, and energy. + * + * Example: + * @code + * serif::network::NetIn input; + * // ... set up input ... + * serif::network::NetOut result = net.evaluate(input); + * @endcode + */ NetOut evaluate(const NetIn &netIn) override; + /** + * @brief Get the vector of unique species in the network. + * @return Reference to the vector of species. + * + * Example: + * @code + * const auto& species = net.getNetworkSpecies(); + * for (const auto& sp : species) { + * std::cout << sp.name() << std::endl; + * } + * @endcode + */ [[nodiscard]] const std::vector& getNetworkSpecies() const; + + /** + * @brief Get the set of REACLIB reactions in the network. + * @return Reference to the set of reactions. + */ [[nodiscard]] const reaclib::REACLIBReactionSet& getNetworkReactions() const; + + /** + * @brief Get the net stoichiometric coefficients for all species in a reaction. + * + * Returns a map from species to their net stoichiometric coefficient (products minus reactants). + * + * @param reaction The REACLIB reaction to analyze. + * @return Map from species to stoichiometric coefficient. + * + * @throws std::runtime_error If a species in the reaction is not found in the network. + * + * Example: + * @code + * for (const auto& reaction : net.getNetworkReactions()) { + * auto stoichiometry = net.getNetReactionStoichiometry(reaction); + * // ... + * } + * @endcode + */ [[nodiscard]] std::unordered_map getNetReactionStoichiometry( const reaclib::REACLIBReaction &reaction) const; + /** + * @brief Check if a species is present in the network. + * @param species The species to check. + * @return True if the species is present, false otherwise. + * + * Example: + * @code + * if (net.involvesSpecies(mySpecies)) { ... } + * @endcode + */ [[nodiscard]] bool involvesSpecies(const serif::atomic::Species& species) const; - [[deprecated("not implimented")]] std::vector> detectCycles() const; - [[deprecated("not implimented")]] std::vector topologicalSortReactions() const; + /** + * @brief Detect cycles in the reaction network (not implemented). + * @return Vector of cycles, each represented as a vector of reaction IDs. + * + * @note This is not yet implemented; however, it will allow for easy detection of things like the CNO cycle. + * @deprecated Not implemented. + */ + [[deprecated("not implimented")]] std::vector> detectCycles() const = delete; + + /** + * @brief Perform a topological sort of the reactions (not implemented). + * @return Vector of reaction IDs in topologically sorted order. + * @deprecated Not implemented. + */ + [[deprecated("not implimented")]] std::vector topologicalSortReactions() const = delete; + + /** + * @brief Export the network to a DOT file for visualization. + * @param filename The name of the output DOT file. + */ + void exportToDot(const std::string& filename) const; + private: - reaclib::REACLIBReactionSet m_reactions; ///< Set of REACLIB reactions for the network. + reaclib::REACLIBReactionSet m_reactions; ///< Set of REACLIB reactions in the network. + std::unordered_map m_reactionIDMap; ///< Map from reaction ID to REACLIBReaction. //PERF: This makes copies of REACLIBReaction and could be a performance bottleneck. - std::vector m_networkSpecies; ///< The species in the network. - std::unordered_map m_networkSpeciesMap; - std::unordered_map m_reactionIDMap; ///< Map of reaction IDs to REACLIB reactions. + std::vector m_networkSpecies; ///< Vector of unique species in the network. + std::unordered_map m_networkSpeciesMap; ///< Map from species name to Species object. + std::unordered_map m_speciesToIndexMap; ///< Map from species to their index in the stoichiometry matrix. - boost::numeric::ublas::compressed_matrix m_stoichiometryMatrix; ///< Stoichiometry matrix for the network. - boost::numeric::ublas::compressed_matrix m_jacobianMatrix; ///< Jacobian matrix for the network. - std::unordered_map m_speciesToIndexMap; ///< Map of species to their index in the stoichiometry matrix. + boost::numeric::ublas::compressed_matrix m_stoichiometryMatrix; ///< Stoichiometry matrix (species x reactions). + boost::numeric::ublas::compressed_matrix m_jacobianMatrix; ///< Jacobian matrix (species x species). - CppAD::ADFun m_rhsADFun; ///< AD tape + CppAD::ADFun m_rhsADFun; ///< CppAD function for the right-hand side of the ODE. + /** + * @brief Functor for ODE right-hand side evaluation. + * + * Used by ODE solvers to compute dY/dt and energy generation rate. This is the only functor used in the non-NSE case. + */ struct ODETerm { GraphNetwork& m_network; ///< Reference to the network - const double m_T9; - const double m_rho; - const size_t m_numSpecies; + const double m_T9; ///< Temperature in 10^9 K + const double m_rho; ///< Density in g/cm^3 + const size_t m_numSpecies; ///< Number of species + /** + * @brief Construct an ODETerm functor. + * @param network Reference to the GraphNetwork. + * @param T9 Temperature in 10^9 K. + * @param rho Density in g/cm^3. + */ ODETerm(GraphNetwork& network, const double T9, double rho) : m_network(network), m_T9(T9), m_rho(rho), m_numSpecies(m_network.m_networkSpecies.size()) {} - void operator()(const boost::numeric::ublas::vector&Y, boost::numeric::ublas::vector& dYdt, double /* t */) { + /** + * @brief Compute dY/dt and energy rate for the ODE solver. + * @param Y Input vector of abundances. + * @param dYdt Output vector for derivatives (resized). + * + * @note this functor does not need auto differentiation to it called the version of calculateAllDerivatives. + */ + void operator()(const boost::numeric::ublas::vector&Y, boost::numeric::ublas::vector& dYdt, double /* t */) const { const std::vector y(Y.begin(), Y.begin() + m_numSpecies); StepDerivatives derivatives = m_network.calculateAllDerivatives(y, m_T9, m_rho); dYdt.resize(m_numSpecies + 1); - std::copy(derivatives.dydt.begin(), derivatives.dydt.end(), dYdt.begin()); + std::ranges::copy(derivatives.dydt, dYdt.begin()); dYdt(m_numSpecies) = derivatives.specificEnergyRate; // Last element is the specific energy rate } }; + /** + * @brief Functor for Jacobian evaluation for stiff ODE solvers. (used in the NSE case) + */ struct JacobianTerm { - GraphNetwork& m_network; - const double m_T9; - const double m_rho; - const size_t m_numSpecies; + GraphNetwork& m_network; ///< Reference to the network + const double m_T9; ///< Temperature in 10^9 K + const double m_rho; ///< Density in g/cm^3 + const size_t m_numSpecies; ///< Number of species + /** + * @brief Construct a JacobianTerm functor. + * @param network Reference to the GraphNetwork. + * @param T9 Temperature in 10^9 K. + * @param rho Density in g/cm^3. + */ JacobianTerm(GraphNetwork& network, const double T9, const double rho) : m_network(network), m_T9(T9), m_rho(rho), m_numSpecies(m_network.m_networkSpecies.size()) {} - void operator()(const boost::numeric::ublas::vector& Y, boost::numeric::ublas::matrix& J, double /* t */, boost::numeric::ublas::vector& /* dfdt */) { + /** + * @brief Compute the Jacobian matrix for the ODE solver. + * @param Y Input vector of abundances. + * @param J Output matrix for the Jacobian (resized). + */ + void operator()(const boost::numeric::ublas::vector& Y, boost::numeric::ublas::matrix& J, double /* t */, boost::numeric::ublas::vector& /* dfdt */) const { const std::vector y_species(Y.begin(), Y.begin() + m_numSpecies); m_network.generateJacobianMatrix(y_species, m_T9, m_rho); @@ -106,48 +277,179 @@ namespace serif::network { } }; + /** + * @brief Struct holding derivatives for a single ODE step. + * @tparam T Scalar type (double or CppAD::AD). + */ template struct StepDerivatives { - std::vector dydt; - T specificEnergyRate = T(0.0); + std::vector dydt; ///< Derivatives of abundances. + T specificEnergyRate = T(0.0); ///< Specific energy generation rate. }; - private: + /** + * @brief Synchronize all internal maps and matrices after network changes. + * + * Called after the reaction set is updated to rebuild all derived data structures. + * + * @note This method must be called any time the network topology changes. + */ void syncInternalMaps(); + + /** + * @brief Collect all unique species from the reaction set. + * + * Populates m_networkSpecies and m_networkSpeciesMap. + * @throws std::runtime_error If a species is not found in the global atomic species database. + * + * @note Called by syncInternalMaps() to ensure the species list is up-to-date. + */ void collectNetworkSpecies(); + + /** + * @brief Populate the reaction ID map from the reaction set. + * + * @note Called by syncInternalMaps() to ensure the reaction ID map is up-to-date. + */ void populateReactionIDMap(); + + /** + * @brief Populate the species-to-index map for matrix construction. + * + * @note Called by syncInternalMaps() to ensure the species-to-index map is up-to-date. + */ void populateSpeciesToIndexMap(); + + /** + * @brief Reserve and resize the Jacobian matrix. + * + * @note Called by syncInternalMaps() to ensure the Jacobian matrix is ready for use. + */ void reserveJacobianMatrix(); + + /** + * @brief Record the CppAD tape for automatic differentiation. + * @throws std::runtime_error If there are no species in the network. + * + * @note Called by syncInternalMaps() to ensure the AD tape is recorded for the current network state. + */ void recordADTape(); // --- Validation methods --- + + /** + * @brief Validate mass and charge conservation for all reactions. + * @return True if all reactions conserve mass and charge, false otherwise. + * + * @note Logs errors for any violations. + */ bool validateConservation() const; + /** + * @brief Validate and update the network composition if needed. + * + * If the composition or culling/temperature has changed, rebuilds the reaction set and synchronizes maps. + * This is primarily used to enable caching in dynamic network situations where the composition, temperature, and density + * may change but the relevant reaction set remains equivalent. + * + * @param composition The composition to validate. + * @param culling Culling threshold. + * @param T9 Temperature in 10^9 K. + */ void validateComposition(const serif::composition::Composition &composition, double culling, double T9); // --- Simple Derivative Calculations --- + + /** + * @brief Calculate all derivatives (dY/dt and energy rate) for the current state. + * @tparam T Scalar type (double or CppAD::AD). + * @param Y Vector of abundances. + * @param T9 Temperature in 10^9 K. + * @param rho Density in g/cm^3. + * @return StepDerivatives containing dY/dt and energy rate. + */ template StepDerivatives calculateAllDerivatives(const std::vector& Y, T T9, T rho) const; // --- Generate the system matrices --- + + /** + * @brief Generate the stoichiometry matrix for the network. + * + * Populates m_stoichiometryMatrix based on the current reactions and species. + * @throws std::runtime_error If a species is not found in the species-to-index map. + */ void generateStoichiometryMatrix(); + + /** + * @brief Generate the Jacobian matrix for the network. + * + * Populates m_jacobianMatrix using automatic differentiation via the precomputed CppAD tape. + * + * @param Y Vector of abundances. + * @param T9 Temperature in 10^9 K. + * @param rho Density in g/cm^3. + */ void generateJacobianMatrix(const std::vector& Y, const double T9, const double rho); + /** + * @brief Calculate the right-hand side (dY/dt) for the ODE system. + * @tparam GeneralScalarType Scalar type (double or CppAD::AD). + * @param Y Vector of abundances. + * @param T9 Temperature in 10^9 K. + * @param rho Density in g/cm^3. + * @return Vector of dY/dt values. + */ template std::vector calculateRHS(const std::vector &Y, const GeneralScalarType T9, const GeneralScalarType rho) const; + + /** + * @brief Calculate the reaction rate for a given reaction in dimensions of particles per cm^3 per second. + * @tparam GeneralScalarType Scalar type (double or CppAD::AD). + * @param reaction The REACLIB reaction. + * @param Y Vector of abundances. + * @param T9 Temperature in 10^9 K. + * @param rho Density in g/cm^3. + * @return Reaction rate. + * @throws std::runtime_error If a reactant species is not found in the species-to-index map. + * + * @note reaclib uses molar units that vary depending on the number of reactants, It's pretty straightforward + * to convert these into particle based units. Specifically, we just need to divide the final result by + * Avogadro's number raised to the number of reactants - 1; + */ template GeneralScalarType calculateReactionRate(const reaclib::REACLIBReaction &reaction, const std::vector &Y, const GeneralScalarType T9, const GeneralScalarType rho) const; + /** + * @brief Detect if the network is stiff and select the appropriate ODE solver. + * + * Heuristically determines stiffness based on the ratio of timescales. The stiffness heuristic is as follows: + * + * 1. For each species, compute the timescale as |Y_i / (dY_i/dt)|, where Y_i is the abundance and dY_i/dt is its time derivative. + * 2. Find the minimum and maximum timescales across all species. + * 3. Compute the stiffness ratio as (maximum timescale) / (minimum timescale). + * 4. If the stiffness ratio exceeds a threshold (default: 1e6), the system is considered stiff and a stiff ODE solver is used. + * Otherwise, a non-stiff ODE solver is used. + * + * This heuristic is based on the observation that stiff systems have widely varying timescales, which can cause explicit + * solvers to become unstable or inefficient. The threshold can be tuned based on the characteristics of the network. + * + * @param netIn The input structure. + * @param T9 Temperature in 10^9 K. + * @param numSpecies Number of species. + * @param Y Vector of abundances. + */ void detectStiff(const NetIn &netIn, double T9, double numSpecies, const boost::numeric::ublas::vector& Y); + }; template - typename GraphNetwork::template StepDerivatives GraphNetwork::calculateAllDerivatives( + GraphNetwork::StepDerivatives GraphNetwork::calculateAllDerivatives( const std::vector &Y, T T9, T rho) const { StepDerivatives result; result.dydt.resize(m_networkSpecies.size(), static_cast(0.0));