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

feat(GridFire): major design changes

Browse Source 
Switching to an Engine + solver design. Also brought xxHash and Eigen in. Working on QSE and Culling.
This commit is contained in:
Emily Boudreaux
2025-06-26 15:13:46 -04:00
parent dd03873bc9
commit cd191cff23
32 changed files with 2737 additions and 1441 deletions
Download Patch File Download Diff File Expand all files Collapse all files

66
src/network/include/gridfire/engine/engine_abstract.h Normal file
View File

@@ -0,0 +1,66 @@
#pragma once
#include "gridfire/network.h" // For NetIn, NetOut
#include "../reaction/reaction.h"
#include "fourdst/composition/composition.h"
#include "fourdst/config/config.h"
#include "fourdst/logging/logging.h"
#include <vector>
#include <unordered_map>
namespace gridfire {
template<typename T>
concept IsArithmeticOrAD = std::is_same_v<T, double> || std::is_same_v<T, CppAD::AD<double>>;
template <IsArithmeticOrAD T>
struct StepDerivatives {
std::vector<T> dydt; ///< Derivatives of abundances.
T nuclearEnergyGenerationRate = T(0.0); ///< Specific energy generation rate.
};
class Engine {
public:
virtual ~Engine() = default;
virtual const std::vector<fourdst::atomic::Species>& getNetworkSpecies() const = 0;
virtual StepDerivatives<double> calculateRHSAndEnergy(
const std::vector<double>& Y,
double T9,
double rho
) const = 0;
};
class DynamicEngine : public Engine {
public:
virtual void generateJacobianMatrix(
const std::vector<double>& Y,
double T9, double rho
) = 0;
virtual double getJacobianMatrixEntry(
int i,
int j
) const = 0;
virtual void generateStoichiometryMatrix() = 0;
virtual int getStoichiometryMatrixEntry(
int speciesIndex,
int reactionIndex
) const = 0;
virtual double calculateMolarReactionFlow(
const reaction::Reaction& reaction,
const std::vector<double>& Y,
double T9,
double rho
) const = 0;
virtual const reaction::REACLIBLogicalReactionSet& getNetworkReactions() const = 0;
virtual std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
const std::vector<double>& Y,
double T9,
double rho
) const = 0;
};
}

0
src/network/include/gridfire/approx8.h → src/network/include/gridfire/engine/engine_approx8.h
View File

55
src/network/include/gridfire/engine/engine_culled.h Normal file
View File

@@ -0,0 +1,55 @@
#pragma once
#include "gridfire/engine/engine_abstract.h"
#include "fourdst/composition/atomicSpecies.h"
namespace gridfire {
class CulledEngine final : public DynamicEngine {
public:
explicit CulledEngine(DynamicEngine& baseEngine);
const std::vector<fourdst::atomic::Species>& getNetworkSpecies() const override;
StepDerivatives<double> calculateRHSAndEnergy(
const std::vector<double> &Y,
double T9,
double rho
) const override;
void generateJacobianMatrix(
const std::vector<double> &Y,
double T9,
double rho
) override;
double getJacobianMatrixEntry(
int i,
int j
) const override;
void generateStoichiometryMatrix() override;
int getStoichiometryMatrixEntry(
int speciesIndex,
int reactionIndex
) const override;
double calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const std::vector<double> &Y,
double T9,
double rho
) const override;
const reaction::REACLIBLogicalReactionSet& getNetworkReactions() const override;
std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
const std::vector<double> &Y,
double T9,
double rho
) const override;
private:
DynamicEngine& m_baseEngine;
std::unordered_map<fourdst::atomic::Species, size_t> m_fullSpeciesToIndexMap;
std::vector<fourdst::atomic::Species> m_culledSpecies;
private:
std::unordered_map<fourdst::atomic::Species, size_t> populatedSpeciesToIndexMap;
std::vector<fourdst::atomic::Species> determineCullableSpecies;
};
}

277
src/network/include/gridfire/engine/engine_graph.h Normal file
View File

@@ -0,0 +1,277 @@
#pragma once
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/composition/composition.h"
#include "fourdst/logging/logging.h"
#include "fourdst/config/config.h"
#include "gridfire/network.h"
#include "gridfire/reaction/reaction.h"
#include "gridfire/engine/engine_abstract.h"
#include <string>
#include <unordered_map>
#include <vector>
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include "cppad/cppad.hpp"
// 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.
namespace gridfire {
typedef CppAD::AD<double> ADDouble; ///< Alias for CppAD AD type for double precision.
using fourdst::config::Config;
using fourdst::logging::LogManager;
using fourdst::constant::Constants;
static constexpr double MIN_DENSITY_THRESHOLD = 1e-18;
static constexpr double MIN_ABUNDANCE_THRESHOLD = 1e-18;
static constexpr double MIN_JACOBIAN_THRESHOLD = 1e-24;
class GraphEngine final : public DynamicEngine{
public:
explicit GraphEngine(const fourdst::composition::Composition &composition);
explicit GraphEngine(reaction::REACLIBLogicalReactionSet reactions);
StepDerivatives<double> calculateRHSAndEnergy(
const std::vector<double>& Y,
const double T9,
const double rho
) const override;
void generateJacobianMatrix(
const std::vector<double>& Y,
const double T9,
const double rho
) override;
void generateStoichiometryMatrix() override;
double calculateMolarReactionFlow(
const reaction::Reaction& reaction,
const std::vector<double>&Y,
const double T9,
const double rho
) const override;
[[nodiscard]] const std::vector<fourdst::atomic::Species>& getNetworkSpecies() const override;
[[nodiscard]] const reaction::REACLIBLogicalReactionSet& getNetworkReactions() const override;
[[nodiscard]] double getJacobianMatrixEntry(
const int i,
const int j
) const override;
[[nodiscard]] std::unordered_map<fourdst::atomic::Species, int> getNetReactionStoichiometry(
const reaction::Reaction& reaction
) const;
[[nodiscard]] int getStoichiometryMatrixEntry(
const int speciesIndex,
const int reactionIndex
) const override;
[[nodiscard]] std::unordered_map<fourdst::atomic::Species, double> getSpeciesTimescales(
const std::vector<double>& Y,
double T9,
double rho
) const override;
[[nodiscard]] bool involvesSpecies(
const fourdst::atomic::Species& species
) const;
void exportToDot(
const std::string& filename
) const;
void exportToCSV(
const std::string& filename
) const;
private:
reaction::REACLIBLogicalReactionSet m_reactions; ///< Set of REACLIB reactions in the network.
std::unordered_map<std::string_view, reaction::Reaction*> 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.
Config& m_config = Config::getInstance();
Constants& m_constants = Constants::getInstance(); ///< Access to physical constants.
quill::Logger* m_logger = LogManager::getInstance().getLogger("log");
private:
void syncInternalMaps();
void collectNetworkSpecies();
void populateReactionIDMap();
void populateSpeciesToIndexMap();
void reserveJacobianMatrix();
void recordADTape();
[[nodiscard]] bool validateConservation() const;
void validateComposition(
const fourdst::composition::Composition &composition,
double culling,
double T9
);
template <IsArithmeticOrAD T>
T calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const std::vector<T> &Y,
const T T9,
const T rho
) const;
template<IsArithmeticOrAD T>
StepDerivatives<T> calculateAllDerivatives(
const std::vector<T> &Y_in,
T T9,
T rho
) const;
StepDerivatives<double> calculateAllDerivatives(
const std::vector<double>& Y_in,
const double T9,
const double rho
) const;
StepDerivatives<ADDouble> calculateAllDerivatives(
const std::vector<ADDouble>& Y_in,
const ADDouble T9,
const ADDouble rho
) const;
};
template<IsArithmeticOrAD T>
StepDerivatives<T> GraphEngine::calculateAllDerivatives(
const std::vector<T> &Y_in, T T9, T rho) const {
// --- Setup output derivatives structure ---
StepDerivatives<T> result;
result.dydt.resize(m_networkSpecies.size(), static_cast<T>(0.0));
// --- AD Pre-setup (flags to control conditionals in an AD safe / branch aware manner) ---
// ----- Constants for AD safe calculations ---
const T zero = static_cast<T>(0.0);
const T one = static_cast<T>(1.0);
// ----- Initialize variables for molar concentration product and thresholds ---
// Note: the logic here is that we use CppAD::CondExprLt to test thresholds and if they are less we set the flag
// to zero so that the final returned reaction flow is 0. This is as opposed to standard if statements
// which create branches that break the AD tape.
const T rho_threshold = static_cast<T>(MIN_DENSITY_THRESHOLD);
// --- Check if the density is below the threshold where we ignore reactions ---
T threshold_flag = CppAD::CondExpLt(rho, rho_threshold, zero, one); // If rho < threshold, set flag to 0
std::vector<T> Y = Y_in;
for (size_t i = 0; i < m_networkSpecies.size(); ++i) {
// We use CppAD::CondExpLt to handle AD taping and prevent branching
// Note that while this is syntactically more complex this is equivalent to
// if (Y[i] < 0) {Y[i] = 0;}
// The issue is that this would introduce a branch which would require the auto diff tape to be re-recorded
// each timestep, which is very inefficient.
Y[i] = CppAD::CondExpLt(Y[i], zero, zero, Y[i]); // Ensure no negative abundances
}
const T u = static_cast<T>(m_constants.get("u").value); // Atomic mass unit in grams
const T N_A = static_cast<T>(m_constants.get("N_a").value); // Avogadro's number in mol^-1
const T c = static_cast<T>(m_constants.get("c").value); // Speed of light in cm/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 molarReactionFlow = calculateMolarReactionFlow<T>(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));
result.dydt[speciesIndex] += threshold_flag * nu_ij * molarReactionFlow / rho;
}
}
T massProductionRate = static_cast<T>(0.0); // [mol][s^-1]
for (const auto& [species, index] : m_speciesToIndexMap) {
massProductionRate += result.dydt[index] * species.mass() * u;
}
result.nuclearEnergyGenerationRate = -massProductionRate * N_A * c * c; // [cm^2][s^-3] = [erg][s^-1][g^-1]
return result;
}
template <IsArithmeticOrAD T>
T GraphEngine::calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const std::vector<T> &Y,
const T T9,
const T rho
) const {
// --- Pre-setup (flags to control conditionals in an AD safe / branch aware manner) ---
// ----- Constants for AD safe calculations ---
const T zero = static_cast<T>(0.0);
const T one = static_cast<T>(1.0);
// ----- Initialize variables for molar concentration product and thresholds ---
// Note: the logic here is that we use CppAD::CondExprLt to test thresholds and if they are less we set the flag
// to zero so that the final returned reaction flow is 0. This is as opposed to standard if statements
// which create branches that break the AD tape.
const T Y_threshold = static_cast<T>(MIN_ABUNDANCE_THRESHOLD);
T threshold_flag = one;
// --- Calculate the molar reaction rate (in units of [s^-1][cm^3(N-1)][mol^(1-N)] for N reactants) ---
const T k_reaction = reaction.calculate_rate(T9);
// --- Cound the number of each reactant species to account for species multiplicity ---
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())]++;
}
// --- Accumulator for the molar concentration ---
auto molar_concentration_product = static_cast<T>(1.0);
// --- Loop through each unique reactant species and calculate the molar concentration for that species then multiply that into the accumulator ---
for (const auto& [species_name, count] : reactant_counts) {
// --- Resolve species to molar abundance ---
// PERF: Could probably optimize out this lookup
const auto species_it = m_speciesToIndexMap.find(m_networkSpeciesMap.at(species_name));
const size_t species_index = species_it->second;
const T Yi = Y[species_index];
// --- Check if the species abundance is below the threshold where we ignore reactions ---
threshold_flag *= CppAD::CondExpLt(Yi, Y_threshold, zero, one);
// --- Convert from molar abundance to molar concentration ---
T molar_concentration = Yi * rho;
// --- If count is > 1 , we need to raise the molar concentration to the power of count since there are really count bodies in that reaction ---
molar_concentration_product *= CppAD::pow(molar_concentration, static_cast<T>(count)); // ni^count
// --- Apply factorial correction for identical reactions ---
if (count > 1) {
molar_concentration_product /= static_cast<T>(std::tgamma(static_cast<double>(count + 1))); // Gamma function for factorial
}
}
// --- Final reaction flow calculation [mol][s^-1][cm^-3] ---
// Note: If the threshold flag ever gets set to zero this will return zero.
// This will result basically in multiple branches being written to the AD tape, which will make
// the tape more expensive to record, but it will also mean that we only need to record it once for
// the entire network.
return molar_concentration_product * k_reaction * threshold_flag;
}
};

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

@@ -1,571 +0,0 @@
#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;
}
};

111
src/network/include/gridfire/network.h
View File

@@ -25,12 +25,15 @@
#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 "fourdst/constants/const.h"
#include "gridfire/reaction/reaction.h"
#include "quill/Logger.h"
#include <unordered_map>
#include "fourdst/constants/const.h"
namespace gridfire {
@@ -127,7 +130,7 @@ namespace gridfire {
* @param netIn Input parameters for the network evaluation.
* @return NetOut Output results from the network evaluation.
*/
virtual NetOut evaluate(const NetIn &netIn);
virtual NetOut evaluate(const NetIn &netIn) = 0;
virtual bool isStiff() const { return m_stiff; }
virtual void setStiff(const bool stiff) { m_stiff = stiff; }
@@ -143,105 +146,9 @@ namespace gridfire {
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);
reaction::REACLIBLogicalReactionSet build_reaclib_nuclear_network(const fourdst::composition::Composition &composition, bool reverse);
reaction::REACLIBLogicalReactionSet build_reaclib_nuclear_network_from_file(const std::string& filename, bool reverse);
} // namespace nuclearNetwork

270
src/network/include/gridfire/reaction/reaction.h Normal file
View File

@@ -0,0 +1,270 @@
#pragma once
#include <string_view>
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/logging/logging.h"
#include "quill/Logger.h"
#include <unordered_map>
#include <vector>
#include <memory>
#include <unordered_set>
#include <cstdint>
#include "cppad/cppad.hpp"
namespace gridfire::reaction {
/**
* @struct REACLIBRateCoefficientSet
* @brief Holds the seven fitting parameters for a single REACLIB rate set.
* @details The thermonuclear reaction rate for a single set is calculated as:
* rate = exp(a0 + a1/T9 + a2/T9^(-1/3) + a3*T9^(1/3) + a4*T9 + a5*T9^(5/3) + a6*ln(T9))
* where T9 is the temperature in billions of Kelvin. The total rate for a
* reaction is the sum of the rates from all its sets.
*/
struct REACLIBRateCoefficientSet {
const double a0;
const double a1;
const double a2;
const double a3;
const double a4;
const double a5;
const double a6;
friend std::ostream& operator<<(std::ostream& os, const REACLIBRateCoefficientSet& r) {
os << "[" << r.a0 << ", " << r.a1 << ", " << r.a2 << ", "
<< r.a3 << ", " << r.a4 << ", " << r.a5 << ", " << r.a6 << "]";
return os;
}
};
class Reaction {
public:
Reaction(
const std::string_view id,
const double qValue,
const std::vector<fourdst::atomic::Species>& reactants,
const std::vector<fourdst::atomic::Species>& products,
const bool reverse = false
);
virtual ~Reaction() = default;
virtual std::unique_ptr<Reaction> clone() const = 0;
virtual double calculate_rate(double T9) const = 0;
virtual CppAD::AD<double> calculate_rate(const CppAD::AD<double> T9) const = 0;
virtual std::string_view peName() const { return ""; }
[[nodiscard]] bool contains(const fourdst::atomic::Species& species) const;
[[nodiscard]] bool contains_reactant(const fourdst::atomic::Species& species) const;
[[nodiscard]] bool contains_product(const fourdst::atomic::Species& species) const;
std::unordered_set<fourdst::atomic::Species> all_species() const;
std::unordered_set<fourdst::atomic::Species> reactant_species() const;
std::unordered_set<fourdst::atomic::Species> product_species() const;
size_t num_species() const;
int stoichiometry(const fourdst::atomic::Species& species) const;
std::unordered_map<fourdst::atomic::Species, int> stoichiometry() const;
std::string_view id() const { return m_id; }
double qValue() const { return m_qValue; }
const std::vector<fourdst::atomic::Species>& reactants() const { return m_reactants; }
const std::vector<fourdst::atomic::Species>& products() const { return m_products; }
bool is_reverse() const { return m_reverse; }
double excess_energy() const;
bool operator==(const Reaction& other) const { return m_id == other.m_id; }
bool operator!=(const Reaction& other) const { return !(*this == other); }
[[nodiscard]] uint64_t hash(uint64_t seed) const;
protected:
quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
std::string m_id;
double m_qValue = 0.0; ///< Q-value of the reaction
std::vector<fourdst::atomic::Species> m_reactants; ///< Reactants of the reaction
std::vector<fourdst::atomic::Species> m_products; ///< Products of the reaction
bool m_reverse = false;
};
class ReactionSet {
public:
explicit ReactionSet(std::vector<std::unique_ptr<Reaction>> reactions);
ReactionSet(const ReactionSet& other);
ReactionSet& operator=(const ReactionSet& other);
virtual ~ReactionSet() = default;
virtual void add_reaction(std::unique_ptr<Reaction> reaction);
virtual void remove_reaction(const std::unique_ptr<Reaction>& reaction);
bool contains(const std::string_view& id) const;
bool contains(const Reaction& reaction) const;
size_t size() const { return m_reactions.size(); }
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 Reaction& operator[](size_t index) const;
[[nodiscard]] const Reaction& operator[](const std::string_view& id) const;
bool operator==(const ReactionSet& other) const;
bool operator!=(const ReactionSet& other) const;
[[nodiscard]] uint64_t hash(uint64_t seed = 0) const;
auto begin() { return m_reactions.begin(); }
auto begin() const { return m_reactions.cbegin(); }
auto end() { return m_reactions.end(); }
auto end() const { return m_reactions.cend(); }
protected:
quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
std::vector<std::unique_ptr<Reaction>> m_reactions;
std::string m_id;
std::unordered_map<std::string, Reaction*> m_reactionNameMap; ///< Maps reaction IDs to Reaction objects for quick lookup
};
/**
* @struct REACLIBReaction
* @brief Represents a single nuclear reaction from the JINA REACLIB database.
* @details This struct is designed to be constructed at compile time (constexpr) from
* the data parsed by the Python generation script. It stores all necessary
* information to identify a reaction and calculate its rate.
*/
class REACLIBReaction final : public Reaction {
public:
/**
* @brief Constructs a REACLIBReaction object at compile time.
* @param id A unique string identifier generated by the Python script.
* @param peName
* @param chapter The REACLIB chapter number, defining the reaction structure.
* @param reactants A vector of strings with the names of the reactant species.
* @param products A vector of strings with the names of the product species.
* @param qValue The Q-value of the reaction in MeV.
* @param label The source label for the rate data (e.g., "wc12w", "st08").
* @param sets A vector of RateFitSet, containing the fitting coefficients for the rate.
* @param reverse A boolean indicating if the reaction is reversed (default is false).
*/
REACLIBReaction(
const std::string_view id,
const std::string_view peName,
const int chapter,
const std::vector<fourdst::atomic::Species> &reactants,
const std::vector<fourdst::atomic::Species> &products,
const double qValue,
const std::string_view label,
const REACLIBRateCoefficientSet &sets,
const bool reverse = false);
[[nodiscard]] std::unique_ptr<Reaction> clone() const override;
[[nodiscard]] double calculate_rate(const double T9) const override;
[[nodiscard]] CppAD::AD<double> calculate_rate(const CppAD::AD<double> T9) const override;
template <typename GeneralScalarType>
[[nodiscard]] GeneralScalarType calculate_rate(const GeneralScalarType T9) const {
const GeneralScalarType T913 = CppAD::pow(T9, 1.0/3.0);
const GeneralScalarType rateExponent = m_rateCoefficients.a0 +
m_rateCoefficients.a1 / T9 +
m_rateCoefficients.a2 / T913 +
m_rateCoefficients.a3 * T913 +
m_rateCoefficients.a4 * T9 +
m_rateCoefficients.a5 * CppAD::pow(T9, 5.0/3.0) +
m_rateCoefficients.a6 * CppAD::log(T9);
return CppAD::exp(rateExponent);
}
[[nodiscard]] std::string_view peName() const override { return m_peName; }
[[nodiscard]] int chapter() const { return m_chapter; }
[[nodiscard]] std::string_view sourceLabel() const { return m_sourceLabel; }
[[nodiscard]] const REACLIBRateCoefficientSet& rateCoefficients() const { return m_rateCoefficients; }
friend std::ostream& operator<<(std::ostream& os, const REACLIBReaction& reaction);
private:
std::string m_peName; ///< Name of the reaction in (projectile, ejectile) notation (e.g. p(p, g)d)
int m_chapter; ///< Chapter number from the REACLIB database, defining the reaction structure.
std::string m_sourceLabel; ///< Source label for the rate data, indicating the origin of the rate coefficients (e.g., "wc12w", "st08").
REACLIBRateCoefficientSet m_rateCoefficients;
};
class REACLIBReactionSet final : public ReactionSet {
public:
explicit REACLIBReactionSet(std::vector<REACLIBReaction>);
std::unordered_set<std::string> peNames() const;
friend std::ostream& operator<<(std::ostream& os, const REACLIBReactionSet& set);
};
class REACLIBLogicalReaction final : public Reaction {
public:
explicit REACLIBLogicalReaction(const std::vector<REACLIBReaction> &reactions);
explicit REACLIBLogicalReaction(const REACLIBReaction &reaction);
void add_reaction(const REACLIBReaction& reaction);
[[nodiscard]] std::unique_ptr<Reaction> clone() const override;
[[nodiscard]] std::string_view peName() const override { return m_id; };
[[nodiscard]] size_t size() const { return m_rates.size(); }
[[nodiscard]] std::vector<std::string> sources() const { return m_sources; }
[[nodiscard]] double calculate_rate(const double T9) const override;
[[nodiscard]] CppAD::AD<double> calculate_rate(const CppAD::AD<double> T9) const override;
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);
// ReSharper disable once CppUseStructuredBinding
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]] int chapter() const { return m_chapter; }
auto begin() { return m_rates.begin(); }
auto begin() const { return m_rates.cbegin(); }
auto end() { return m_rates.end(); }
auto end() const { return m_rates.cend(); }
private:
int m_chapter;
std::vector<std::string> m_sources;
std::vector<REACLIBRateCoefficientSet> m_rates;
};
class REACLIBLogicalReactionSet final : public ReactionSet {
public:
REACLIBLogicalReactionSet() = delete;
explicit REACLIBLogicalReactionSet(const REACLIBReactionSet& reactionSet);
[[nodiscard]] std::unordered_set<std::string> peNames() const;
private:
std::unordered_set<std::string> m_peNames;
};
}

256
src/network/include/gridfire/solver/solver.h Normal file
View File

@@ -0,0 +1,256 @@
#pragma once
#include "gridfire/engine/engine_graph.h"
#include "gridfire/engine/engine_abstract.h"
#include "gridfire/network.h"
#include "fourdst/logging/logging.h"
#include "fourdst/config/config.h"
#include "quill/Logger.h"
#include "Eigen/Dense"
#include <vector>
namespace gridfire::solver {
struct dynamicQSESpeciesIndices {
std::vector<size_t> dynamicSpeciesIndices; // Slow species that are not in QSE
std::vector<size_t> QSESpeciesIndices; // Fast species that are in QSE
};
template <typename EngineT>
class NetworkSolverStrategy {
public:
explicit NetworkSolverStrategy(EngineT& engine) : m_engine(engine) {};
virtual ~NetworkSolverStrategy() = default;
virtual NetOut evaluate(const NetIn& netIn) = 0;
protected:
EngineT& m_engine;
};
using DynamicNetworkSolverStrategy = NetworkSolverStrategy<DynamicEngine>;
using StaticNetworkSolverStrategy = NetworkSolverStrategy<Engine>;
class QSENetworkSolver final : public DynamicNetworkSolverStrategy {
public:
using DynamicNetworkSolverStrategy::DynamicNetworkSolverStrategy;
NetOut evaluate(const NetIn& netIn) override;
private: // methods
dynamicQSESpeciesIndices packSpeciesTypeIndexVectors(
const std::vector<double>& Y,
const double T9,
const double rho
) const;
Eigen::VectorXd calculateSteadyStateAbundances(
const std::vector<double>& Y,
const double T9,
const double rho,
const dynamicQSESpeciesIndices& indices
) const;
NetOut initializeNetworkWithShortIgnition(
const NetIn& netIn
) const;
private: // Nested functors for ODE integration
struct RHSFunctor {
DynamicEngine& m_engine;
const std::vector<size_t>& m_dynamicSpeciesIndices;
const std::vector<size_t>& m_QSESpeciesIndices;
const Eigen::VectorXd& m_Y_QSE;
const double m_T9;
const double m_rho;
RHSFunctor(
DynamicEngine& engine,
const std::vector<size_t>& dynamicSpeciesIndices,
const std::vector<size_t>& QSESpeciesIndices,
const Eigen::VectorXd& Y_QSE,
const double T9,
const double rho
) :
m_engine(engine),
m_dynamicSpeciesIndices(dynamicSpeciesIndices),
m_QSESpeciesIndices(QSESpeciesIndices),
m_Y_QSE(Y_QSE),
m_T9(T9),
m_rho(rho) {}
void operator()(
const boost::numeric::ublas::vector<double>& YDynamic,
boost::numeric::ublas::vector<double>& dYdtDynamic,
double t
) const;
};
struct JacobianFunctor {
DynamicEngine& m_engine;
const std::vector<size_t>& m_dynamicSpeciesIndices;
const std::vector<size_t>& m_QSESpeciesIndices;
const double m_T9;
const double m_rho;
JacobianFunctor(
DynamicEngine& engine,
const std::vector<size_t>& dynamicSpeciesIndices,
const std::vector<size_t>& QSESpeciesIndices,
const double T9,
const double rho
) :
m_engine(engine),
m_dynamicSpeciesIndices(dynamicSpeciesIndices),
m_QSESpeciesIndices(QSESpeciesIndices),
m_T9(T9),
m_rho(rho) {}
void operator()(
const boost::numeric::ublas::vector<double>& YDynamic,
boost::numeric::ublas::matrix<double>& JDynamic,
double t,
boost::numeric::ublas::vector<double>& dfdt
) const;
};
template<typename T>
struct EigenFunctor {
using InputType = Eigen::Matrix<T, Eigen::Dynamic, 1>;
using OutputType = Eigen::Matrix<T, Eigen::Dynamic, 1>;
using JacobianType = Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>;
DynamicEngine& m_engine;
const std::vector<double>& m_YDynamic;
const std::vector<size_t>& m_dynamicSpeciesIndices;
const std::vector<size_t>& m_QSESpeciesIndices;
const double m_T9;
const double m_rho;
EigenFunctor(
DynamicEngine& engine,
const std::vector<double>& YDynamic,
const std::vector<size_t>& dynamicSpeciesIndices,
const std::vector<size_t>& QSESpeciesIndices,
const double T9,
const double rho
) :
m_engine(engine),
m_YDynamic(YDynamic),
m_dynamicSpeciesIndices(dynamicSpeciesIndices),
m_QSESpeciesIndices(QSESpeciesIndices),
m_T9(T9),
m_rho(rho) {}
int operator()(const InputType& v_QSE, OutputType& f_QSE) const;
int df(const InputType& v_QSE, JacobianType& J_QSE) const;
};
private:
quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
fourdst::config::Config& m_config = fourdst::config::Config::getInstance();
};
class DirectNetworkSolver final : public DynamicNetworkSolverStrategy {
public:
using DynamicNetworkSolverStrategy::DynamicNetworkSolverStrategy;
NetOut evaluate(const NetIn& netIn) override;
private:
struct RHSFunctor {
DynamicEngine& m_engine;
const double m_T9;
const double m_rho;
const size_t m_numSpecies;
RHSFunctor(
DynamicEngine& engine,
const double T9,
const double rho
) :
m_engine(engine),
m_T9(T9),
m_rho(rho),
m_numSpecies(engine.getNetworkSpecies().size()) {}
void operator()(
const boost::numeric::ublas::vector<double>& Y,
boost::numeric::ublas::vector<double>& dYdt,
double t
) const;
};
struct JacobianFunctor {
DynamicEngine& m_engine;
const double m_T9;
const double m_rho;
const size_t m_numSpecies;
JacobianFunctor(
DynamicEngine& engine,
const double T9,
const double rho
) :
m_engine(engine),
m_T9(T9),
m_rho(rho),
m_numSpecies(engine.getNetworkSpecies().size()) {}
void operator()(
const boost::numeric::ublas::vector<double>& Y,
boost::numeric::ublas::matrix<double>& J,
double t,
boost::numeric::ublas::vector<double>& dfdt
) const;
};
private:
quill::Logger* m_logger = fourdst::logging::LogManager::getInstance().getLogger("log");
fourdst::config::Config& m_config = fourdst::config::Config::getInstance();
};
template<typename T>
int QSENetworkSolver::EigenFunctor<T>::operator()(const InputType &v_QSE, OutputType &f_QSE) const {
std::vector<double> YFull(m_engine.getNetworkSpecies().size(), 0.0);
Eigen::VectorXd Y_QSE(v_QSE.array().exp());
for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) {
YFull[m_dynamicSpeciesIndices[i]] = m_YDynamic[i];
}
for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
YFull[m_QSESpeciesIndices[i]] = Y_QSE(i);
}
const auto [full_dYdt, specificEnergyGenerationRate] = m_engine.calculateRHSAndEnergy(YFull, m_T9, m_rho);
f_QSE.resize(m_QSESpeciesIndices.size());
for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
f_QSE(i) = full_dYdt[m_QSESpeciesIndices[i]];
}
return 0; // Success
}
template <typename T>
int QSENetworkSolver::EigenFunctor<T>::df(const InputType& v_QSE, JacobianType& J_QSE) const {
std::vector<double> YFull(m_engine.getNetworkSpecies().size(), 0.0);
Eigen::VectorXd Y_QSE(v_QSE.array().exp());
for (size_t i = 0; i < m_dynamicSpeciesIndices.size(); ++i) {
YFull[m_dynamicSpeciesIndices[i]] = m_YDynamic[i];
}
for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
YFull[m_QSESpeciesIndices[i]] = Y_QSE(i);
}
m_engine.generateJacobianMatrix(YFull, m_T9, m_rho);
Eigen::MatrixXd J_orig(m_QSESpeciesIndices.size(), m_QSESpeciesIndices.size());
for (size_t i = 0; i < m_QSESpeciesIndices.size(); ++i) {
for (size_t j = 0; j < m_QSESpeciesIndices.size(); ++j) {
J_orig(i, j) = m_engine.getJacobianMatrixEntry(m_QSESpeciesIndices[i], m_QSESpeciesIndices[j]);
}
}
J_QSE = J_orig;
for (long j = 0; j < J_QSE.cols(); ++j) {
J_QSE.col(j) *= Y_QSE(j); // Chain rule for log space
}
return 0; // Success
}
}