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

fix(GraphNetwork): working on loads of small bugs

Browse Source 
Fized stoichiometry matrix initialization, added penames to reablib reactions, began work on LogicalReaction to sum the contributions of different fitting functions provided by reaclib
This commit is contained in:
Emily Boudreaux
2025-06-23 15:18:56 -04:00
parent 9f6e360b0f
commit dd03873bc9
31 changed files with 101449 additions and 19120 deletions
Download Patch File Download Diff File Expand all files Collapse all files

331
src/network/include/gridfire/approx8.h Normal file
View File

@@ -0,0 +1,331 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <array>
#include <boost/numeric/odeint.hpp>
#include "gridfire/network.h"
/**
* @file approx8.h
* @brief Header file for the Approx8 nuclear reaction network.
*
* This file contains the definitions and declarations for the Approx8 nuclear reaction network.
* The network is based on Frank Timmes' "approx8" and includes 8 isotopes and various nuclear reactions.
* The rates are evaluated using a fitting function with coefficients from reaclib.jinaweb.org.
*/
namespace gridfire::approx8{
/**
* @typedef vector_type
* @brief Alias for a vector of doubles using Boost uBLAS.
*/
typedef boost::numeric::ublas::vector< double > vector_type;
/**
* @typedef matrix_type
* @brief Alias for a matrix of doubles using Boost uBLAS.
*/
typedef boost::numeric::ublas::matrix< double > matrix_type;
/**
* @typedef vec7
* @brief Alias for a std::array of 7 doubles.
*/
typedef std::array<double,7> vec7;
/**
* @struct Approx8Net
* @brief Contains constants and arrays related to the nuclear network.
*/
struct Approx8Net{
static constexpr int ih1=0;
static constexpr int ihe3=1;
static constexpr int ihe4=2;
static constexpr int ic12=3;
static constexpr int in14=4;
static constexpr int io16=5;
static constexpr int ine20=6;
static constexpr int img24=7;
static constexpr int iTemp=img24+1;
static constexpr int iDensity =iTemp+1;
static constexpr int iEnergy=iDensity+1;
static constexpr int nIso=img24+1; // number of isotopes
static constexpr int nVar=iEnergy+1; // number of variables
static constexpr std::array<int,nIso> aIon = {
1,
3,
4,
12,
14,
16,
20,
24
};
static constexpr std::array<double,nIso> mIon = {
1.67262164e-24,
5.00641157e-24,
6.64465545e-24,
1.99209977e-23,
2.32462686e-23,
2.65528858e-23,
3.31891077e-23,
3.98171594e-23
};
};
/**
* @brief Multiplies two arrays and sums the resulting elements.
* @param a First array.
* @param b Second array.
* @return Sum of the product of the arrays.
* @example
* @code
* vec7 a = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0};
* vec7 b = {0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5};
* double result = sum_product(a, b);
* @endcode
*/
double sum_product( const vec7 &a, const vec7 &b);
/**
* @brief Returns an array of T9 terms for the nuclear reaction rate fit.
* @param T Temperature in GigaKelvin.
* @return Array of T9 terms.
* @example
* @code
* double T = 1.5;
* vec7 T9_array = get_T9_array(T);
* @endcode
*/
vec7 get_T9_array(const double &T);
/**
* @brief Evaluates the nuclear reaction rate given the T9 array and coefficients.
* @param T9 Array of T9 terms.
* @param coef Array of coefficients.
* @return Evaluated rate.
* @example
* @code
* vec7 T9 = get_T9_array(1.5);
* vec7 coef = {1.0, 0.1, 0.01, 0.001, 0.0001, 0.00001, 0.000001};
* double rate = rate_fit(T9, coef);
* @endcode
*/
double rate_fit(const vec7 &T9, const vec7 &coef);
/**
* @brief Calculates the rate for the reaction p + p -> d.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double pp_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction p + d -> he3.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double dp_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction he3 + he3 -> he4 + 2p.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double he3he3_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction he3(he3,2p)he4.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double he3he4_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction he4 + he4 + he4 -> c12.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double triple_alpha_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12 + p -> n13.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12p_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12 + he4 -> o16.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n14(p,g)o15 - o15 + p -> c12 + he4.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n14p_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n14(a,g)f18 assumed to go on to ne20.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n14a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n15(p,a)c12 (CNO I).
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n15pa_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction n15(p,g)o16 (CNO II).
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double n15pg_rate(const vec7 &T9);
/**
* @brief Calculates the fraction for the reaction n15(p,g)o16.
* @param T9 Array of T9 terms.
* @return Fraction of the reaction.
*/
double n15pg_frac(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction o16(p,g)f17 then f17 -> o17(p,a)n14.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double o16p_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction o16(a,g)ne20.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double o16a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction ne20(a,g)mg24.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double ne20a_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12(c12,a)ne20.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12c12_rate(const vec7 &T9);
/**
* @brief Calculates the rate for the reaction c12(o16,a)mg24.
* @param T9 Array of T9 terms.
* @return Rate of the reaction.
*/
double c12o16_rate(const vec7 &T9);
/**
* @struct Jacobian
* @brief Functor to calculate the Jacobian matrix for implicit solvers.
*/
struct Jacobian {
/**
* @brief Calculates the Jacobian matrix.
* @param y State vector.
* @param J Jacobian matrix.
* @param dfdt Derivative of the state vector.
*/
void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const;
};
/**
* @struct ODE
* @brief Functor to calculate the derivatives for the ODE solver.
*/
struct ODE {
/**
* @brief Calculates the derivatives of the state vector.
* @param y State vector.
* @param dydt Derivative of the state vector.
*/
void operator() ( const vector_type &y, vector_type &dydt, double /* t */) const;
};
/**
* @class Approx8Network
* @brief Class for the Approx8 nuclear reaction network.
*/
class Approx8Network final : public Network {
public:
Approx8Network();
/**
* @brief Evaluates the nuclear network.
* @param netIn Input parameters for the network.
* @return Output results from the network.
*/
NetOut evaluate(const NetIn &netIn) override;
/**
* @brief Sets whether the solver should use a stiff method.
* @param stiff Boolean indicating if a stiff method should be used.
*/
void setStiff(bool stiff) override;
/**
* @brief Checks if the solver is using a stiff method.
* @return Boolean indicating if a stiff method is being used.
*/
bool isStiff() const override { return m_stiff; }
private:
vector_type m_y;
double m_tMax = 0;
double m_dt0 = 0;
bool m_stiff = false;
/**
* @brief Converts the input parameters to the internal state vector.
* @param netIn Input parameters for the network.
* @return Internal state vector.
*/
static vector_type convert_netIn(const NetIn &netIn);
};
} // namespace nnApprox8

571
src/network/include/gridfire/netgraph.h Normal file
View File

@@ -0,0 +1,571 @@
#pragma once
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/composition/composition.h"
#include "gridfire/network.h"
#include "gridfire/reaclib.h"
#include <string>
#include <unordered_map>
#include <vector>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include "cppad/cppad.hpp"
#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 gridfire::GraphNetwork
*/
namespace gridfire {
/**
* @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 fourdst::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 fourdst::composition::Composition &composition,
const double cullingThreshold, const double T9);
explicit GraphNetwork(const reaclib::REACLIBReactionSet& reactions);
/**
* @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<fourdst::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<fourdst::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 fourdst::atomic::Species& species) 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")]] [[nodiscard]] 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")]] [[nodiscard]] 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 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<fourdst::atomic::Species> m_networkSpecies; ///< Vector of unique species in the network.
std::unordered_map<std::string_view, fourdst::atomic::Species> m_networkSpeciesMap; ///< Map from species name to Species object.
std::unordered_map<fourdst::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.
/**
* @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; ///< 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()) {}
/**
* @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(), m_numSpecies + Y.begin());
StepDerivatives<double> derivatives = m_network.calculateAllDerivatives<double>(y, m_T9, m_rho);
dYdt.resize(m_numSpecies + 1);
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; ///< 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()) {}
/**
* @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);
J.resize(m_numSpecies + 1, m_numSpecies + 1);
J.clear(); // Zero out the matrix
// PERF: Could probably be optimized
for (auto it1 = m_network.m_jacobianMatrix.begin1(); it1 != m_network.m_jacobianMatrix.end1(); ++it1) {
for (auto it2 = it1.begin(); it2 != it1.end(); ++it2) {
J(it2.index1(), it2.index2()) = *it2;
}
}
}
};
/**
* @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; ///< 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.
*/
[[nodiscard]] 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 fourdst::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, 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,
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, unsigned long numSpecies, const boost::numeric::ublas::vector<double>& Y);
};
template<IsArithmeticOrAD T>
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));
if (rho < MIN_DENSITY_THRESHOLD) {
return result; // Return zero for everything if density is too low
}
const T u = static_cast<T>(m_constants.get("u").value); // Atomic mass unit in grams
const T MeV_to_erg = static_cast<T>(m_constants.get("MeV_to_erg").value);
T volumetricEnergyRate = static_cast<T>(0.0); // Accumulator for erg / (cm^3 * s)
// --- SINGLE LOOP OVER ALL REACTIONS ---
for (size_t reactionIndex = 0; reactionIndex < m_reactions.size(); ++reactionIndex) {
const auto& reaction = m_reactions[reactionIndex];
// 1. Calculate reaction rate
const T reactionRate = calculateReactionRate(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) {
const T nu_ij = static_cast<T>(m_stoichiometryMatrix(speciesIndex, reactionIndex));
if (nu_ij != 0) {
const T speciesAtomicMassAMU = static_cast<T>(m_networkSpecies[speciesIndex].mass());
const T speciesAtomicMassGrams = speciesAtomicMassAMU * u;
result.dydt[speciesIndex] += (nu_ij * reactionRate * speciesAtomicMassGrams) / rho;
}
}
// 3. Use the same rate to update the energy generation rate
const T q_value_ergs = static_cast<T>(reaction.qValue()) * MeV_to_erg;
volumetricEnergyRate += reactionRate * q_value_ergs;
}
result.specificEnergyRate = volumetricEnergyRate / rho;
return result;
}
template <IsArithmeticOrAD GeneralScalarType>
GeneralScalarType GraphNetwork::calculateReactionRate(const reaclib::REACLIBReaction &reaction, const std::vector<GeneralScalarType> &Y,
const GeneralScalarType T9, const GeneralScalarType rho) const {
const auto &constants = fourdst::constant::Constants::getInstance();
const auto u = constants.get("u"); // Atomic mass unit in g/mol
const auto uValue = static_cast<GeneralScalarType>(u.value); // Convert to double for calculations
const GeneralScalarType k_reaction = reaction.calculate_rate(T9);
auto reactant_product_or_particle_densities = static_cast<GeneralScalarType>(1.0);
std::unordered_map<std::string, int> reactant_counts;
reactant_counts.reserve(reaction.reactants().size());
for (const auto& reactant : reaction.reactants()) {
reactant_counts[std::string(reactant.name())]++;
}
const auto minAbundanceThreshold = static_cast<GeneralScalarType>(MIN_ABUNDANCE_THRESHOLD);
const auto minDensityThreshold = static_cast<GeneralScalarType>(MIN_DENSITY_THRESHOLD);
if (rho < minDensityThreshold) {
// LOG_INFO(m_logger, "Density is below the minimum threshold ({} g/cm^3), returning zero reaction rate for reaction '{}'.",
// CppAD::Value(rho), reaction.id()); // CppAD::Value causes a compilation error here, not clear why...
return static_cast<GeneralScalarType>(0.0); // If density is below a threshold, return zero rate.
}
for (const auto& [species_name, count] : reactant_counts) {
auto species_it = m_speciesToIndexMap.find(m_networkSpeciesMap.at(species_name));
if (species_it == m_speciesToIndexMap.end()) {
LOG_ERROR(m_logger, "Reactant species '{}' not found in species to index map for reaction '{}'.",
species_name, reaction.id());
throw std::runtime_error("Reactant species not found in species to index map: " + species_name);
}
const size_t species_index = species_it->second;
const GeneralScalarType Yi = Y[species_index];
if (Yi < minAbundanceThreshold) {
return static_cast<GeneralScalarType>(0.0); // If any reactant is below a threshold, return zero rate.
}
auto atomicMassAMU = static_cast<GeneralScalarType>(m_networkSpecies[species_index].mass());
// Convert to number density
GeneralScalarType ni;
const GeneralScalarType denominator = atomicMassAMU * uValue;
assert (denominator > 0.0);
ni = (Yi * rho) / (denominator);
reactant_product_or_particle_densities *= ni;
// Apply factorial correction for identical reactions
if (count > 1) {
reactant_product_or_particle_densities /= static_cast<GeneralScalarType>(std::tgamma(static_cast<double>(count + 1))); // Gamma function for factorial
}
}
const GeneralScalarType Na = static_cast<GeneralScalarType>(constants.get("N_a").value); // Avogadro's number in mol^-1
const int numReactants = static_cast<int>(reaction.reactants().size());
auto molarCorrectionFactor = static_cast<GeneralScalarType>(1.0); // No correction needed for single reactant reactions
if (numReactants > 1) {
molarCorrectionFactor = CppAD::pow(Na, static_cast<GeneralScalarType>(reaction.reactants().size() - 1));
}
return (reactant_product_or_particle_densities * k_reaction) / molarCorrectionFactor; // reaction rate in per volume per time (particles/cm^3/s)
}
template <IsArithmeticOrAD GeneralScalarType>
std::vector<GeneralScalarType> GraphNetwork::calculateRHS(
const std::vector<GeneralScalarType>& Y,
const GeneralScalarType T9,
const GeneralScalarType rho) const {
auto derivatives = calculateAllDerivatives<GeneralScalarType>(Y, T9, rho);
return derivatives.dydt;
}
};

247
src/network/include/gridfire/network.h Normal file
View File

@@ -0,0 +1,247 @@
/* ***********************************************************************
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Authors: Emily Boudreaux, Aaron Dotter
// Last Modified: March 21, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
// License version 3 (GPLv3) as published by the Free Software Foundation.
//
// 4DSSE is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU Library General Public License for more details.
//
// You should have received a copy of the GNU Library General Public License
// along with this software; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// *********************************************************************** */
#pragma once
#include <vector>
#include "fourdst/logging/logging.h"
#include "fourdst/config/config.h"
#include "fourdst/composition/species.h"
#include "quill/Logger.h"
#include "fourdst/composition/composition.h"
#include "gridfire/reaclib.h"
#include <unordered_map>
#include "fourdst/constants/const.h"
namespace gridfire {
enum NetworkFormat {
APPROX8, ///< Approx8 nuclear reaction network format.
REACLIB, ///< General REACLIB nuclear reaction network format.
UNKNOWN,
};
static inline std::unordered_map<NetworkFormat, std::string> FormatStringLookup = {
{APPROX8, "Approx8"},
{REACLIB, "REACLIB"},
{UNKNOWN, "Unknown"}
};
/**
* @struct NetIn
* @brief Input structure for the network evaluation.
*
* This structure holds the input parameters required for the network evaluation.
*
* Example usage:
* @code
* nuclearNetwork::NetIn netIn;
* netIn.composition = {1.0, 0.0, 0.0};
* netIn.tmax = 1.0e6;
* netIn.dt0 = 1.0e-3;
* netIn.temperature = 1.0e8;
* netIn.density = 1.0e5;
* netIn.energy = 1.0e12;
* @endcode
*/
struct NetIn {
fourdst::composition::Composition composition; ///< Composition of the network
double tMax; ///< Maximum time
double dt0; ///< Initial time step
double temperature; ///< Temperature in Kelvin
double density; ///< Density in g/cm^3
double energy; ///< Energy in ergs
double culling = 0.0; ///< Culling threshold for reactions (default is 0.0, meaning no culling)
};
/**
* @struct NetOut
* @brief Output structure for the network evaluation.
*
* This structure holds the output results from the network evaluation.
*
* Example usage:
* @code
* nuclearNetwork::NetOut netOut = network.evaluate(netIn);
* std::vector<double> composition = netOut.composition;
* int steps = netOut.num_steps;
* double energy = netOut.energy;
* @endcode
*/
struct NetOut {
fourdst::composition::Composition composition; ///< Composition of the network after evaluation
int num_steps; ///< Number of steps taken in the evaluation
double energy; ///< Energy in ergs after evaluation
friend std::ostream& operator<<(std::ostream& os, const NetOut& netOut) {
os << "NetOut(composition=" << netOut.composition << ", num_steps=" << netOut.num_steps << ", energy=" << netOut.energy << ")";
return os;
}
};
/**
* @class Network
* @brief Class for network evaluation.
*
* This class provides methods to evaluate the network based on the input parameters.
*
* Example usage:
* @code
* nuclearNetwork::Network network;
* nuclearNetwork::NetIn netIn;
* // Set netIn parameters...
* nuclearNetwork::NetOut netOut = network.evaluate(netIn);
* @endcode
*/
class Network {
public:
explicit Network(const NetworkFormat format = NetworkFormat::APPROX8);
virtual ~Network() = default;
[[nodiscard]] NetworkFormat getFormat() const;
NetworkFormat setFormat(const NetworkFormat format);
/**
* @brief Evaluate the network based on the input parameters.
*
* @param netIn Input parameters for the network evaluation.
* @return NetOut Output results from the network evaluation.
*/
virtual NetOut evaluate(const NetIn &netIn);
virtual bool isStiff() const { return m_stiff; }
virtual void setStiff(const bool stiff) { m_stiff = stiff; }
protected:
fourdst::config::Config& m_config; ///< Configuration instance
fourdst::logging::LogManager& m_logManager; ///< Log manager instance
quill::Logger* m_logger; ///< Logger instance
NetworkFormat m_format; ///< Format of the network
fourdst::constant::Constants& m_constants;
bool m_stiff = false; ///< Flag indicating if the network is stiff
};
class LogicalReaction {
public:
explicit LogicalReaction(const std::vector<reaclib::REACLIBReaction> &reactions);
explicit LogicalReaction(const reaclib::REACLIBReaction &reaction);
void add_reaction(const reaclib::REACLIBReaction& reaction);
template <typename T>
[[nodiscard]] T calculate_rate(const T T9) const {
T sum = static_cast<T>(0.0);
const T T913 = CppAD::pow(T9, 1.0/3.0);
const T T953 = CppAD::pow(T9, 5.0/3.0);
const T logT9 = CppAD::log(T9);
for (const auto& rate : m_rates) {
const T exponent = rate.a0 +
rate.a1 / T9 +
rate.a2 / T913 +
rate.a3 * T913 +
rate.a4 * T9 +
rate.a5 * T953 +
rate.a6 * logT9;
sum += CppAD::exp(exponent);
}
return sum;
}
[[nodiscard]] const std::string& id() const {return std::string(m_peID); }
[[nodiscard]] std::string_view peName() const { return m_peID; }
[[nodiscard]] int chapter() const { return m_chapter; }
[[nodiscard]] const std::vector<fourdst::atomic::Species>& reactants() const { return m_reactants; }
[[nodiscard]] const std::vector<fourdst::atomic::Species>& products() const { return m_products; }
[[nodiscard]] double qValue() const { return m_qValue; }
[[nodiscard]] bool is_reverse() const { return m_reverse; }
auto begin();
auto begin() const;
auto end();
auto end() const;
private:
const quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
std::string_view m_peID;
std::vector<fourdst::atomic::Species> m_reactants; ///< Reactants of the reaction
std::vector<fourdst::atomic::Species> m_products; ///< Products of the reaction
double m_qValue = 0.0; ///< Q-value of the reaction
bool m_reverse = false; ///< True if the reaction is reverse
int m_chapter;
std::vector<reaclib::RateFitSet> m_rates;
};
class LogicalReactionSet {
public:
LogicalReactionSet() = default;
explicit LogicalReactionSet(const std::vector<LogicalReaction>& reactions);
explicit LogicalReactionSet(const std::vector<reaclib::REACLIBReaction>& reactions);
explicit LogicalReactionSet(const reaclib::REACLIBReactionSet& reactionSet);
void add_reaction(const LogicalReaction& reaction);
void add_reaction(const reaclib::REACLIBReaction& reaction);
void remove_reaction(const LogicalReaction& reaction);
[[nodiscard]] bool contains(const std::string_view& id) const;
[[nodiscard]] bool contains(const LogicalReaction& reactions) const;
[[nodiscard]] bool contains(const reaclib::REACLIBReaction& reaction) const;
[[nodiscard]] size_t size() const;
void sort(double T9=1.0);
bool contains_species(const fourdst::atomic::Species &species) const;
bool contains_reactant(const fourdst::atomic::Species &species) const;
bool contains_product(const fourdst::atomic::Species &species) const;
[[nodiscard]] const LogicalReaction& operator[](size_t index) const;
[[nodiscard]] const LogicalReaction& operator[](const std::string_view& id) const;
auto begin();
auto begin() const;
auto end();
auto end() const;
private:
const quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
std::vector<LogicalReaction> m_reactions;
std::unordered_map<std::string_view, LogicalReaction&> m_reactionNameMap;
};
LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition);
LogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, double culling, double T9 = 1.0);
} // namespace nuclearNetwork