tboudreaux/GridFire
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

docs(GraphNetwork): added and rebuilt docs

Browse Source
This commit is contained in:
Emily Boudreaux
2025-06-21 05:04:14 -04:00
parent fe73a021bf
commit 1cd0176b01
2 changed files with 382 additions and 26 deletions
Download Patch File Download Diff File Expand all files Collapse all files

54
src/network/private/netgraph.cpp
View File

@@ -15,6 +15,7 @@
#include <string_view>
#include <unordered_map>
#include <vector>
#include <fstream>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/odeint.hpp>
@@ -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;

354
src/network/public/netgraph.h
View File

@@ -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<double>.
*
* Used to enable templated functions for both arithmetic and automatic differentiation types.
*/
template<typename T>
concept IsArithmeticOrAD = std::is_same_v<T, double> || std::is_same_v<T, CppAD::AD<double>>;
/// 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<serif::atomic::Species>& 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<serif::atomic::Species, int> 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<std::vector<std::string>> detectCycles() const;
[[deprecated("not implimented")]] std::vector<std::string> 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<std::vector<std::string>> 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<std::string> 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<std::string_view, const reaclib::REACLIBReaction> m_reactionIDMap; ///< Map from reaction ID to REACLIBReaction. //PERF: This makes copies of REACLIBReaction and could be a performance bottleneck.
std::vector<serif::atomic::Species> m_networkSpecies; ///< The species in the network.
std::unordered_map<std::string_view, serif::atomic::Species> m_networkSpeciesMap;
std::unordered_map<std::string_view, const reaclib::REACLIBReaction> m_reactionIDMap; ///< Map of reaction IDs to REACLIB reactions.
std::vector<serif::atomic::Species> m_networkSpecies; ///< Vector of unique species in the network.
std::unordered_map<std::string_view, serif::atomic::Species> m_networkSpeciesMap; ///< Map from species name to Species object.
std::unordered_map<serif::atomic::Species, size_t> m_speciesToIndexMap; ///< Map from species to their index in the stoichiometry matrix.
boost::numeric::ublas::compressed_matrix<int> m_stoichiometryMatrix; ///< Stoichiometry matrix for the network.
boost::numeric::ublas::compressed_matrix<double> m_jacobianMatrix; ///< Jacobian matrix for the network.
std::unordered_map<serif::atomic::Species, size_t> m_speciesToIndexMap; ///< Map of species to their index in the stoichiometry matrix.
boost::numeric::ublas::compressed_matrix<int> m_stoichiometryMatrix; ///< Stoichiometry matrix (species x reactions).
boost::numeric::ublas::compressed_matrix<double> m_jacobianMatrix; ///< Jacobian matrix (species x species).
CppAD::ADFun<double> m_rhsADFun; ///< AD tape
CppAD::ADFun<double> 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<double>&Y, boost::numeric::ublas::vector<double>& 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 <double> version of calculateAllDerivatives.
*/
void operator()(const boost::numeric::ublas::vector<double>&Y, boost::numeric::ublas::vector<double>& dYdt, double /* t */) const {
const std::vector<double> y(Y.begin(), Y.begin() + m_numSpecies);
StepDerivatives<double> derivatives = m_network.calculateAllDerivatives<double>(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<double>& Y, boost::numeric::ublas::matrix<double>& J, double /* t */, boost::numeric::ublas::vector<double>& /* 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<double>& Y, boost::numeric::ublas::matrix<double>& J, double /* t */, boost::numeric::ublas::vector<double>& /* dfdt */) const {
const std::vector<double> 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<double>).
*/
template <IsArithmeticOrAD T>
struct StepDerivatives {
std::vector<T> dydt;
T specificEnergyRate = T(0.0);
std::vector<T> 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<double>).
* @param Y Vector of abundances.
* @param T9 Temperature in 10^9 K.
* @param rho Density in g/cm^3.
* @return StepDerivatives<T> containing dY/dt and energy rate.
*/
template <IsArithmeticOrAD T>
StepDerivatives<T> calculateAllDerivatives(const std::vector<T>& 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<double>& 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<double>).
* @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 <IsArithmeticOrAD GeneralScalarType>
std::vector<GeneralScalarType> calculateRHS(const std::vector<GeneralScalarType> &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<double>).
* @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 <IsArithmeticOrAD GeneralScalarType>
GeneralScalarType calculateReactionRate(const reaclib::REACLIBReaction &reaction,
const std::vector<GeneralScalarType> &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<double>& Y);
};
template<IsArithmeticOrAD T>
typename GraphNetwork::template StepDerivatives<T> GraphNetwork::calculateAllDerivatives(
GraphNetwork::StepDerivatives<T> GraphNetwork::calculateAllDerivatives(
const std::vector<T> &Y, T T9, T rho) const {
StepDerivatives<T> result;
result.dydt.resize(m_networkSpecies.size(), static_cast<T>(0.0));