Graph-based nuclear reaction network using REACLIB reactions. More...

#include <netgraph.h>

Inheritance diagram for serif::network::GraphNetwork:

Classes
struct	JacobianTerm
	Functor for Jacobian evaluation for stiff ODE solvers. (used in the NSE case) More...

struct	ODETerm
	Functor for ODE right-hand side evaluation. More...

struct	StepDerivatives
	Struct holding derivatives for a single ODE step. More...

Public Member Functions
	GraphNetwork (const serif::composition::Composition &composition)
	Construct a GraphNetwork from a composition.

	GraphNetwork (const serif::composition::Composition &composition, const double cullingThreshold, const double T9)
	Construct a GraphNetwork from a composition with culling and temperature.

NetOut	evaluate (const NetIn &netIn) override
	Evolve the network for the given input conditions.

const std::vector< serif::atomic::Species > &	getNetworkSpecies () const
	Get the vector of unique species in the network.

const reaclib::REACLIBReactionSet &	getNetworkReactions () const
	Get the set of REACLIB reactions in the network.

std::unordered_map< serif::atomic::Species, int >	getNetReactionStoichiometry (const reaclib::REACLIBReaction &reaction) const
	Get the net stoichiometric coefficients for all species in a reaction.

bool	involvesSpecies (const serif::atomic::Species &species) const
	Check if a species is present in the network.

std::vector< std::vector< std::string > >	detectCycles () const
	Detect cycles in the reaction network (not implemented).

std::vector< std::string >	topologicalSortReactions () const
	Perform a topological sort of the reactions (not implemented).

Public Member Functions inherited from serif::network::Network
	Network (const NetworkFormat format=NetworkFormat::APPROX8)

virtual	~Network ()=default

NetworkFormat	getFormat () const

NetworkFormat	setFormat (const NetworkFormat format)

virtual bool	isStiff () const

virtual void	setStiff (const bool stiff)

Private Member Functions
void	syncInternalMaps ()
	Synchronize all internal maps and matrices after network changes.

void	collectNetworkSpecies ()
	Collect all unique species from the reaction set.

void	populateReactionIDMap ()
	Populate the reaction ID map from the reaction set.

void	populateSpeciesToIndexMap ()
	Populate the species-to-index map for matrix construction.

void	reserveJacobianMatrix ()
	Reserve and resize the Jacobian matrix.

void	recordADTape ()
	Record the CppAD tape for automatic differentiation.

bool	validateConservation () const
	Validate mass and charge conservation for all reactions.

void	validateComposition (const serif::composition::Composition &composition, double culling, double T9)
	Validate and update the network composition if needed.

template<IsArithmeticOrAD T>
StepDerivatives< T >	calculateAllDerivatives (const std::vector< T > &Y, T T9, T rho) const
	Calculate all derivatives (dY/dt and energy rate) for the current state.

void	generateStoichiometryMatrix ()
	Generate the stoichiometry matrix for the network.

void	generateJacobianMatrix (const std::vector< double > &Y, const double T9, const double rho)
	Generate the Jacobian matrix for the network.

template<IsArithmeticOrAD GeneralScalarType>
std::vector< GeneralScalarType >	calculateRHS (const std::vector< GeneralScalarType > &Y, const GeneralScalarType T9, const GeneralScalarType rho) const
	Calculate the right-hand side (dY/dt) for the ODE system.

template<IsArithmeticOrAD GeneralScalarType>
GeneralScalarType	calculateReactionRate (const reaclib::REACLIBReaction &reaction, const std::vector< GeneralScalarType > &Y, const GeneralScalarType T9, const GeneralScalarType rho) const
	Calculate the reaction rate for a given reaction in dimensions of particles per cm^3 per second.

void	detectStiff (const NetIn &netIn, double T9, double numSpecies, const boost::numeric::ublas::vector< double > &Y)
	Detect if the network is stiff and select the appropriate ODE solver.

Private Attributes
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
	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 (species x reactions).

boost::numeric::ublas::compressed_matrix< double >	m_jacobianMatrix
	Jacobian matrix (species x species).

CppAD::ADFun< double >	m_rhsADFun
	CppAD function for the right-hand side of the ODE.

Additional Inherited Members
Protected Attributes inherited from serif::network::Network
serif::config::Config &	m_config
	Configuration instance.

serif::probe::LogManager &	m_logManager
	Log manager instance.

quill::Logger *	m_logger
	Logger instance.

NetworkFormat	m_format
	Format of the network.

serif::constant::Constants &	m_constants

bool	m_stiff = false
	Flag indicating if the network is stiff.

Detailed Description

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:

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);

Definition at line 77 of file netgraph.h.

Constructor & Destructor Documentation

◆ GraphNetwork() [1/2]

serif::network::GraphNetwork::GraphNetwork ( const serif::composition::Composition & composition )

explicit

Construct a GraphNetwork from a composition.

Parameters

composition The composition specifying the initial isotopic abundances.

Definition at line 24 of file netgraph.cpp.

◆ GraphNetwork() [2/2]

serif::network::GraphNetwork::GraphNetwork	(	const serif::composition::Composition &	composition,
		const double	cullingThreshold,
		const double	T9 )

explicit

Construct a GraphNetwork from a composition with culling and temperature.

Parameters

composition	The composition specifying the initial isotopic abundances.
cullingThreshold	specific reaction rate threshold to cull reactions below.
T9	Temperature in units of 10^9 K where culling is evaluated at.

See also: serif::network::build_reaclib_nuclear_network

Definition at line 32 of file netgraph.cpp.

Member Function Documentation

◆ calculateAllDerivatives()

template<IsArithmeticOrAD T>

GraphNetwork::StepDerivatives< T > serif::network::GraphNetwork::calculateAllDerivatives	(	const std::vector< T > &	Y,
		T	T9,
		T	rho ) const

private

Calculate all derivatives (dY/dt and energy rate) for the current state.

Template Parameters

T	Scalar type (double or CppAD::AD<double>).

Parameters

Y	Vector of abundances.
T9	Temperature in 10^9 K.
rho	Density in g/cm^3.

Returns: StepDerivatives<T> containing dY/dt and energy rate.

Definition at line 445 of file netgraph.h.

◆ calculateReactionRate()

template<IsArithmeticOrAD GeneralScalarType>

GeneralScalarType serif::network::GraphNetwork::calculateReactionRate	(	const reaclib::REACLIBReaction &	reaction,
		const std::vector< GeneralScalarType > &	Y,
		const GeneralScalarType	T9,
		const GeneralScalarType	rho ) const

private

Calculate the reaction rate for a given reaction in dimensions of particles per cm^3 per second.

Template Parameters

GeneralScalarType Scalar type (double or CppAD::AD<double>).

Parameters

reaction	The REACLIB reaction.
Y	Vector of abundances.
T9	Temperature in 10^9 K.
rho	Density in g/cm^3.

Returns: Reaction rate.

Exceptions

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 Avagadro's number raised to the number of reactants - 1;

Definition at line 488 of file netgraph.h.

◆ calculateRHS()

template<IsArithmeticOrAD GeneralScalarType>

std::vector< GeneralScalarType > serif::network::GraphNetwork::calculateRHS	(	const std::vector< GeneralScalarType > &	Y,
		const GeneralScalarType	T9,
		const GeneralScalarType	rho ) const

private

Calculate the right-hand side (dY/dt) for the ODE system.

Template Parameters

GeneralScalarType Scalar type (double or CppAD::AD<double>).

Parameters

Y	Vector of abundances.
T9	Temperature in 10^9 K.
rho	Density in g/cm^3.

Returns: Vector of dY/dt values.

Definition at line 555 of file netgraph.h.

◆ collectNetworkSpecies()

void serif::network::GraphNetwork::collectNetworkSpecies ( )

private

Collect all unique species from the reaction set.

Populates m_networkSpecies and m_networkSpeciesMap.

Exceptions

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.

Definition at line 51 of file netgraph.cpp.

◆ detectCycles()

std::vector< std::vector< std::string > > serif::network::GraphNetwork::detectCycles ( ) const

Detect cycles in the reaction network (not implemented).

Returns: 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.

◆ detectStiff()

void serif::network::GraphNetwork::detectStiff	(	const NetIn &	netIn,
		double	T9,
		double	numSpecies,
		const boost::numeric::ublas::vector< double > &	Y )

private

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:

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.
Find the minimum and maximum timescales across all species.
Compute the stiffness ratio as (maximum timescale) / (minimum timescale).
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.

Parameters

netIn	The input structure.
T9	Temperature in 10^9 K.
numSpecies	Number of species.
Y	Vector of abundances.

Definition at line 305 of file netgraph.cpp.

◆ evaluate()

NetOut serif::network::GraphNetwork::evaluate ( const NetIn & netIn )

overridevirtual

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.

Parameters

netIn The input structure containing composition, temperature, density, and integration parameters.

Returns: NetOut The output structure containing the final composition, number of steps, and energy.

Example:

serif::network::NetIn input;
// ... set up input ...
serif::network::NetOut result = net.evaluate(input);

Reimplemented from serif::network::Network.

Definition at line 337 of file netgraph.cpp.

◆ generateJacobianMatrix()

void serif::network::GraphNetwork::generateJacobianMatrix	(	const std::vector< double > &	Y,
		const double	T9,
		const double	rho )

private

Generate the Jacobian matrix for the network.

Populates m_jacobianMatrix using automatic differentiation via the precomputed CppAD tape.

Parameters

Y	Vector of abundances.
T9	Temperature in 10^9 K.
rho	Density in g/cm^3.

Definition at line 275 of file netgraph.cpp.

◆ generateStoichiometryMatrix()

void serif::network::GraphNetwork::generateStoichiometryMatrix ( )

private

Generate the stoichiometry matrix for the network.

Populates m_stoichiometryMatrix based on the current reactions and species.

Exceptions

std::runtime_error If a species is not found in the species-to-index map.

Definition at line 232 of file netgraph.cpp.

◆ getNetReactionStoichiometry()

std::unordered_map< serif::atomic::Species, int > serif::network::GraphNetwork::getNetReactionStoichiometry ( const reaclib::REACLIBReaction & reaction ) const

nodiscard

Get the net stoichiometric coefficients for all species in a reaction.

Returns a map from species to their net stoichiometric coefficient (products minus reactants).

Parameters

reaction The REACLIB reaction to analyze.

Returns: Map from species to stoichiometric coefficient.

Exceptions

std::runtime_error If a species in the reaction is not found in the network.

Example:

for (const auto& reaction : net.getNetworkReactions()) {
    auto stoichiometry = net.getNetReactionStoichiometry(reaction);
    // ...
}

Definition at line 125 of file netgraph.cpp.

◆ getNetworkReactions()

const reaclib::REACLIBReactionSet & serif::network::GraphNetwork::getNetworkReactions ( ) const

nodiscard

Get the set of REACLIB reactions in the network.

Returns: Reference to the set of reactions.

Definition at line 112 of file netgraph.cpp.

◆ getNetworkSpecies()

const std::vector< serif::atomic::Species > & serif::network::GraphNetwork::getNetworkSpecies ( ) const

nodiscard

Get the vector of unique species in the network.

Returns: Reference to the vector of species.

Example:

const auto& species = net.getNetworkSpecies();
for (const auto& sp : species) {
    std::cout << sp.name() << std::endl;
}

Definition at line 106 of file netgraph.cpp.

◆ involvesSpecies()

bool serif::network::GraphNetwork::involvesSpecies ( const serif::atomic::Species & species ) const

nodiscard

Check if a species is present in the network.

Parameters

species The species to check.

Returns: True if the species is present, false otherwise.

Example:

if (net.involvesSpecies(mySpecies)) { ... }

Definition at line 118 of file netgraph.cpp.

◆ populateReactionIDMap()

void serif::network::GraphNetwork::populateReactionIDMap ( )

private

Populate the reaction ID map from the reaction set.

Note: Called by syncInternalMaps() to ensure the reaction ID map is up-to-date.

Definition at line 79 of file netgraph.cpp.

◆ populateSpeciesToIndexMap()

void serif::network::GraphNetwork::populateSpeciesToIndexMap ( )

private

Populate the species-to-index map for matrix construction.

Note: Called by syncInternalMaps() to ensure the species-to-index map is up-to-date.

Definition at line 88 of file netgraph.cpp.

◆ recordADTape()

void serif::network::GraphNetwork::recordADTape ( )

private

Record the CppAD tape for automatic differentiation.

Exceptions

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.

Definition at line 419 of file netgraph.cpp.

◆ reserveJacobianMatrix()

void serif::network::GraphNetwork::reserveJacobianMatrix ( )

private

Reserve and resize the Jacobian matrix.

Note: Called by syncInternalMaps() to ensure the Jacobian matrix is ready for use.

Definition at line 95 of file netgraph.cpp.

◆ syncInternalMaps()

void serif::network::GraphNetwork::syncInternalMaps ( )

private

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.

Definition at line 42 of file netgraph.cpp.

◆ topologicalSortReactions()

std::vector< std::string > serif::network::GraphNetwork::topologicalSortReactions ( ) const

Perform a topological sort of the reactions (not implemented).

Returns: Vector of reaction IDs in topologically sorted order.

Deprecated: Not implemented.

◆ validateComposition()

void serif::network::GraphNetwork::validateComposition	(	const serif::composition::Composition &	composition,
		double	culling,
		double	T9 )

private

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.

Parameters

composition	The composition to validate.
culling	Culling threshold.
T9	Temperature in 10^9 K.

Definition at line 210 of file netgraph.cpp.

◆ validateConservation()

bool serif::network::GraphNetwork::validateConservation ( ) const

private

Validate mass and charge conservation for all reactions.

Returns: True if all reactions conserve mass and charge, false otherwise.

Note: Logs errors for any violations.

Definition at line 155 of file netgraph.cpp.

Member Data Documentation

◆ m_jacobianMatrix

boost::numeric::ublas::compressed_matrix<double> serif::network::GraphNetwork::m_jacobianMatrix

private

Jacobian matrix (species x species).

Definition at line 192 of file netgraph.h.

◆ m_networkSpecies

std::vector<serif::atomic::Species> serif::network::GraphNetwork::m_networkSpecies

private

Vector of unique species in the network.

Definition at line 187 of file netgraph.h.

◆ m_networkSpeciesMap

std::unordered_map<std::string_view, serif::atomic::Species> serif::network::GraphNetwork::m_networkSpeciesMap

private

Map from species name to Species object.

Definition at line 188 of file netgraph.h.

◆ m_reactionIDMap

std::unordered_map<std::string_view, const reaclib::REACLIBReaction> serif::network::GraphNetwork::m_reactionIDMap

private

Map from reaction ID to REACLIBReaction. //PERF: This makes copies of REACLIBReaction and could be a performance bottleneck.

Definition at line 185 of file netgraph.h.

◆ m_reactions

reaclib::REACLIBReactionSet serif::network::GraphNetwork::m_reactions

private

Set of REACLIB reactions in the network.

Definition at line 184 of file netgraph.h.

◆ m_rhsADFun

CppAD::ADFun<double> serif::network::GraphNetwork::m_rhsADFun

private

CppAD function for the right-hand side of the ODE.

Definition at line 194 of file netgraph.h.

◆ m_speciesToIndexMap

std::unordered_map<serif::atomic::Species, size_t> serif::network::GraphNetwork::m_speciesToIndexMap

private

Map from species to their index in the stoichiometry matrix.

Definition at line 189 of file netgraph.h.

◆ m_stoichiometryMatrix

boost::numeric::ublas::compressed_matrix<int> serif::network::GraphNetwork::m_stoichiometryMatrix

private

Stoichiometry matrix (species x reactions).

Definition at line 191 of file netgraph.h.

The documentation for this class was generated from the following files:

src/network/public/netgraph.h
src/network/private/netgraph.cpp

Classes

Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ GraphNetwork() [1/2]

◆ GraphNetwork() [2/2]

Member Function Documentation

◆ calculateAllDerivatives()

◆ calculateReactionRate()

◆ calculateRHS()

◆ collectNetworkSpecies()

◆ detectCycles()

◆ detectStiff()

◆ evaluate()

◆ generateJacobianMatrix()

◆ generateStoichiometryMatrix()

◆ getNetReactionStoichiometry()

◆ getNetworkReactions()

◆ getNetworkSpecies()

◆ involvesSpecies()

◆ populateReactionIDMap()

◆ populateSpeciesToIndexMap()

◆ recordADTape()

◆ reserveJacobianMatrix()

◆ syncInternalMaps()

◆ topologicalSortReactions()

◆ validateComposition()

◆ validateConservation()

Member Data Documentation

◆ m_jacobianMatrix

◆ m_networkSpecies

◆ m_networkSpeciesMap

◆ m_reactionIDMap

◆ m_reactions

◆ m_rhsADFun

◆ m_speciesToIndexMap

◆ m_stoichiometryMatrix