diff --git a/src/network/private/netgraph.cpp b/src/network/private/netgraph.cpp index da6496b0..a855dfdb 100644 --- a/src/network/private/netgraph.cpp +++ b/src/network/private/netgraph.cpp @@ -1,63 +1,54 @@ #include "netgraph.h" #include "atomicSpecies.h" -#include "reaclib.h" -#include "network.h" -#include "species.h" #include "const.h" +#include "network.h" +#include "reaclib.h" +#include "species.h" #include "quill/LogMacros.h" -#include -#include -#include -#include -#include -#include #include #include +#include +#include +#include +#include +#include +#include -#include "const.h" +#include +#include namespace serif::network { - template - GraphNetwork::GraphNetwork( + GraphNetwork::GraphNetwork( const serif::composition::Composition &composition ): Network(REACLIB), - m_reactions(build_reaclib_nuclear_network(composition)), - m_initialComposition(composition), - m_currentComposition(composition), - m_cullingThreshold(0), - m_T9(0) { + m_reactions(build_reaclib_nuclear_network(composition)) { syncInternalMaps(); } - template - GraphNetwork::GraphNetwork( + GraphNetwork::GraphNetwork( const serif::composition::Composition &composition, const double cullingThreshold, const double T9 ): Network(REACLIB), - m_reactions(build_reaclib_nuclear_network(composition, cullingThreshold, T9)), - m_initialComposition(composition), - m_currentComposition(composition), - m_cullingThreshold(cullingThreshold), - m_T9(T9) { + m_reactions(build_reaclib_nuclear_network(composition, cullingThreshold, T9)) { syncInternalMaps(); } - template - void GraphNetwork::syncInternalMaps() { + void GraphNetwork::syncInternalMaps() { collectNetworkSpecies(); populateReactionIDMap(); populateSpeciesToIndexMap(); + reserveJacobianMatrix(); + recordADTape(); } // --- Network Graph Construction Methods --- - template - void GraphNetwork::collectNetworkSpecies() { + void GraphNetwork::collectNetworkSpecies() { m_networkSpecies.clear(); m_networkSpeciesMap.clear(); @@ -85,8 +76,7 @@ namespace serif::network { } - template - void GraphNetwork::populateReactionIDMap() { + void GraphNetwork::populateReactionIDMap() { LOG_INFO(m_logger, "Populating reaction ID map for REACLIB graph network (serif::network::GraphNetwork)..."); m_reactionIDMap.clear(); for (const auto& reaction: m_reactions) { @@ -102,14 +92,17 @@ namespace serif::network { } } - void GraphNetwork::buildNetworkGraph() { - LOG_INFO(m_logger, "Building network graph..."); - m_reactions = build_reaclib_nuclear_network(m_initialComposition, m_cullingThreshold, m_T9); - syncInternalMaps(); - LOG_INFO(m_logger, "Network graph built with {} reactions and {} species.", m_reactions.size(), m_networkSpecies.size()); + void GraphNetwork::reserveJacobianMatrix() { + // The implementation of this function (and others) constrains this nuclear network to a constant temperature and density during + // each evaluation. + size_t numSpecies = m_networkSpecies.size(); + m_jacobianMatrix.clear(); + m_jacobianMatrix.resize(numSpecies, numSpecies, false); // Sparse matrix, no initial values + LOG_INFO(m_logger, "Jacobian matrix resized to {} rows and {} columns.", + m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); } - // --- Basic Accessors and Queries --- + // --- Basic Accessors and Queries --- const std::vector& GraphNetwork::getNetworkSpecies() const { // Returns a constant reference to the vector of unique species in the network. LOG_DEBUG(m_logger, "Providing access to network species vector. Size: {}.", m_networkSpecies.size()); @@ -137,7 +130,7 @@ namespace serif::network { for (const auto& reactant : reaction.reactants()) { auto it = m_networkSpeciesMap.find(reactant.name()); if (it != m_networkSpeciesMap.end()) { - stoichiometry[it->second]--; // Copy Species by value (PERF: Future performance improvements by using pointers or references) + stoichiometry[it->second]--; // Copy Species by value (PERF: Future performance improvements by using pointers or references (std::reference_wrapper) or something like that) } else { LOG_WARNING(m_logger, "Reactant species '{}' in reaction '{}' not found in network species map during stoichiometry calculation.", reactant.name(), reaction.id()); @@ -148,7 +141,7 @@ namespace serif::network { for (const auto& product : reaction.products()) { auto it = m_networkSpeciesMap.find(product.name()); if (it != m_networkSpeciesMap.end()) { - stoichiometry[it->second]++; // Copy Species by value (PERF: Future performance improvements by using pointers or references) + stoichiometry[it->second]++; // Copy Species by value (PERF: Future performance improvements by using pointers or references (std::reference_wrapper) or something like that) } else { LOG_WARNING(m_logger, "Product species '{}' in reaction '{}' not found in network species map during stoichiometry calculation.", product.name(), reaction.id()); @@ -214,6 +207,21 @@ namespace serif::network { return true; // All reactions passed the conservation check } + void GraphNetwork::validateComposition(const serif::composition::Composition &composition, double culling, double T9) { + + // Check if the requested network has already been cached. + // PERF: Rebuilding this should be pretty fast but it may be a good point of optimization in the future. + const reaclib::REACLIBReactionSet validationReactionSet = build_reaclib_nuclear_network(composition, culling, T9); + + + // This allows for dynamic network modification while retaining caching for networks which are very similar. + if (validationReactionSet != m_reactions) { + LOG_INFO(m_logger, "Reaction set not cached. Rebuilding the reaction set for T9={} and culling={}.", T9, culling); + m_reactions = validationReactionSet; + syncInternalMaps(); // Re-sync internal maps after updating reactions. Note this will also retrace the AD tape. + } + } + // --- Generate Stoichiometry Matrix --- void GraphNetwork::generateStoichiometryMatrix() { LOG_INFO(m_logger, "Generating stoichiometry matrix..."); @@ -258,100 +266,191 @@ namespace serif::network { m_stoichiometryMatrix.nnz()); // Assuming nnz() exists for compressed_matrix } - void GraphNetwork::generateJacobianMatrix() { - LOG_INFO(m_logger, "Generating jacobian matrix..."); + void GraphNetwork::generateJacobianMatrix(const std::vector &Y, const double T9, + const double rho) { - // Task 1: Set dimensions and initialize the matrix - size_t numSpecies = m_networkSpecies.size(); - m_stoichiometryMatrix.resize(numSpecies, numSpecies, false); - - LOG_INFO(m_logger, "Jacobian matrix initialized with dimensions: {} rows (species) x {} columns (species).", - numSpecies, numSpecies); - } - - std::vector GraphNetwork::calculateRHS(const std::vector& Y, const double T9, const double rho) const { - std::vector dotY(m_networkSpecies.size(), 0.0); - const size_t numReactions = m_reactions.size(); + LOG_INFO(m_logger, "Generating jacobian matrix for T9={}, rho={}..", T9, rho); const size_t numSpecies = m_networkSpecies.size(); - if (m_stoichiometryMatrix.size1() != numSpecies || m_stoichiometryMatrix.size2() != numReactions) { - LOG_ERROR(m_logger, "Stoichiometry matrix dimensions do not match network species and reactions sizes."); - throw std::runtime_error("Stoichiometry matrix dimensions mismatch."); + // 1. Pack the input variables into a vector for CppAD + std::vector adInput(numSpecies + 2, 0.0); // +2 for T9 and rho + for (size_t i = 0; i < numSpecies; ++i) { + adInput[i] = Y[i]; } + adInput[numSpecies] = T9; // T9 + adInput[numSpecies + 1] = rho; // rho - for (size_t reactionIndex = 0; reactionIndex < numReactions; ++reactionIndex) { - const auto& currentReaction = m_reactions[reactionIndex]; - const double reactionRatePerVolumePerTime = calculateReactionRate(currentReaction, Y, T9, rho); + // 2. Calculate the full jacobian + const std::vector dotY = m_rhsADFun.Jacobian(adInput); - // dY_i/dt = ∑ υ_ij * R_j * mass_in_grams - // TODO: make sure the unit conversion is correct after calculate reaction rate. - - for (size_t speciesIndex = 0; speciesIndex < numSpecies; ++speciesIndex) { - const int nu_ij = m_stoichiometryMatrix(speciesIndex, reactionIndex); - if (nu_ij != 0) { - const double speciesAtomicMassAMU = m_networkSpecies[speciesIndex].mass(); - dotY[speciesIndex] += (nu_ij * reactionRatePerVolumePerTime * speciesAtomicMassAMU)/rho; + // 3. Pack jacobian vector into sparse matrix + m_jacobianMatrix.clear(); + for (size_t i = 0; i < numSpecies; ++i) { + for (size_t j = 0; j < numSpecies; ++j) { + const double value = dotY[i * (numSpecies + 2) + j]; + if (std::abs(value) > MIN_JACOBIAN_THRESHOLD) { + m_jacobianMatrix(i, j) = value; } } } - return dotY; + LOG_INFO(m_logger, "Jacobian matrix generated with dimensions: {} rows x {} columns.", m_jacobianMatrix.size1(), m_jacobianMatrix.size2()); } - double GraphNetwork::calculateReactionRate(const reaclib::REACLIBReaction &reaction, const std::vector &Y, - const double T9, const double rho) const { - const auto &constants = serif::constant::Constants::getInstance(); + void GraphNetwork::detectStiff(const NetIn &netIn, const double T9, const double numSpecies, const boost::numeric::ublas::vector& Y) { + // --- Heuristic for automatic stiffness detection --- + const std::vector initial_y_stl(Y.begin(), Y.begin() + numSpecies); // Copy only the species abundances, not the specific energy rate + const auto derivatives = calculateAllDerivatives(initial_y_stl, T9, netIn.density); + const std::vector& initial_dotY = derivatives.dydt; - const auto u = constants.get("u"); // Atomic mass unit in g/mol - const double k_reaction = reaction.calculate_rate(T9); // PERF: Consider precomputing all of these and putting them into an O(1) lookup table. + double min_timescale = std::numeric_limits::max(); + double max_timescale = 0.0; - double reactant_product = 1.0; - - std::unordered_map reactant_counts; - reactant_counts.reserve(reaction.reactants().size()); - for (const auto& reactant : reaction.reactants()) { - reactant_counts[std::string(reactant.name())]++; - } - - for (const auto& [species_name, count] : reactant_counts) { - constexpr double minThreshold = 1e-18; - 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 double Yi = Y[species_index]; - double Ai = m_networkSpecies[species_index].a(); - - if (Yi < minThreshold) { - return 0.0; // If any reactant is below a threshold, return zero rate. - } - - double atomicMassAMU = m_networkSpecies[species_index].mass(); - - // Convert to number density - double ni; - const double denominator = atomicMassAMU * u.value; - if (denominator > minThreshold) { - ni = (Yi * rho) / (denominator); - } else { - ni = 0.0; - } - - reactant_product *= ni; - - // Apply factorial correction for identical reactions - if (count > 1) { - reactant_product /= static_cast(std::tgamma(count + 1)); // Gamma function for factorial + for (size_t i = 0; i < numSpecies; ++i) { + if (std::abs(initial_dotY[i]) > 0) { + const double timescale = std::abs(Y(i) / initial_dotY[i]); + if (timescale > max_timescale) {max_timescale = timescale;} + if (timescale < min_timescale) {min_timescale = timescale;} } } - const double Na = constants.get("N_a").value; // Avogadro's number in mol^-1 - const double molarCorrectionFactor = std::pow(Na, reaction.reactants().size() - 1); - return (reactant_product * k_reaction) / molarCorrectionFactor; // reaction rate in per volume per time (particles/cm^3/s) + + const double stiffnessRatio = max_timescale / min_timescale; + + // TODO: Pull this out into a configuration option or a more sophisticated heuristic. + constexpr double stiffnessThreshold = 1.0e6; // This is a heuristic threshold, can be tuned based on network characteristics. + + LOG_INFO(m_logger, "Stiffness ratio is {} (max timescale: {}, min timescale: {}).", stiffnessRatio, max_timescale, min_timescale); + if (stiffnessRatio > stiffnessThreshold) { + LOG_INFO(m_logger, "Network is detected to be stiff. Using stiff ODE solver."); + m_stiff = true; + } else { + LOG_INFO(m_logger, "Network is detected to be non-stiff. Using non-stiff ODE solver."); + m_stiff = false; + } } NetOut GraphNetwork::evaluate(const NetIn &netIn) { - return Network::evaluate(netIn); + namespace ublas = boost::numeric::ublas; + namespace odeint = boost::numeric::odeint; + + const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) + validateComposition(netIn.composition, netIn.culling, T9); + + const double numSpecies = m_networkSpecies.size(); + constexpr double abs_tol = 1.0e-8; + constexpr double rel_tol = 1.0e-8; + + int stepCount = 0; + + // TODO: Pull these out into configuration options + + ODETerm rhs_functor(*this, T9, netIn.density); + + + ublas::vector Y(numSpecies + 1); + for (size_t i = 0; i < numSpecies; ++i) { + const auto& species = m_networkSpecies[i]; + // Get the mass fraction for this specific species from the input object + Y(i) = netIn.composition.getMassFraction(std::string(species.name())); + } + Y(numSpecies) = 0; // initial specific energy rate, will be updated later + + detectStiff(netIn, T9, numSpecies, Y); + + if (m_stiff) { + JacobianTerm jacobian_functor(*this, T9, netIn.density); + LOG_INFO(m_logger, "Making use of stiff ODE solver for network evaluation."); + auto stepper = odeint::make_controlled>(abs_tol, rel_tol); + stepCount = odeint::integrate_adaptive( + stepper, + std::make_pair(rhs_functor, jacobian_functor), + Y, + 0.0, // Start time + netIn.tMax, + netIn.dt0 + ); + + } else { + LOG_INFO(m_logger, "Making use of ODE solver (non-stiff) for network evaluation."); + using state_type = ublas::vector; + auto stepper = odeint::make_controlled>(abs_tol, rel_tol); + stepCount = odeint::integrate_adaptive( + stepper, + rhs_functor, + Y, + 0.0, // Start time + netIn.tMax, + netIn.dt0 + ); + + + } + + double sumY = 0.0; + for (int i = 0; i < numSpecies; ++i) { sumY += Y(i); } + for (int i = 0; i < numSpecies; ++i) { Y(i) /= sumY; } + + // --- Marshall output variables --- + // PERF: Im sure this step could be tuned to avoid so many copies, that is a job for another day + std::vector speciesNames; + speciesNames.reserve(numSpecies); + for (const auto& species : m_networkSpecies) { + speciesNames.push_back(std::string(species.name())); + } + + std::vector finalAbundances(Y.begin(), Y.begin() + numSpecies); + serif::composition::Composition outputComposition(speciesNames, finalAbundances); + outputComposition.finalize(true); + + NetOut netOut; + netOut.composition = outputComposition; + netOut.num_steps = stepCount; + netOut.energy = Y(numSpecies); // The last element in Y is the specific energy rate + + return netOut; + + } + + void GraphNetwork::recordADTape() { + LOG_INFO(m_logger, "Recording AD tape for the RHS calculation..."); + + // Task 1: Set dimensions and initialize the matrix + const size_t numSpecies = m_networkSpecies.size(); + if (numSpecies == 0) { + LOG_ERROR(m_logger, "Cannot record AD tape: No species in the network."); + throw std::runtime_error("Cannot record AD tape: No species in the network."); + } + const size_t numADInputs = numSpecies + 2; // Note here that by not letting T9 and rho be independent variables, we are constraining the network to a constant temperature and density during each evaluation. + + // --- CppAD Tape Recording --- + // 1. Declare independent variable (adY) + // We also initialize the dummy variable for tape recording (these tell CppAD what the derivative chain looks like). + // Their numeric values are irrelevant except for in so far as they avoid numerical instabilities. + + // Distribute total mass fraction uniformly between species in the dummy variable space + const auto uniformMassFraction = static_cast>(1.0 / numSpecies); + std::vector> adInput(numADInputs, uniformMassFraction); + adInput[numSpecies] = 1.0; // Dummy T9 + adInput[numSpecies + 1] = 1.0; // Dummy rho + + // 3. Declare independent variables (what CppAD will differentiate wrt.) + // This also beings the tape recording process. + CppAD::Independent(adInput); + + std::vector> adY(numSpecies); + for(size_t i = 0; i < numSpecies; ++i) { + adY[i] = adInput[i]; + } + const CppAD::AD adT9 = adInput[numSpecies]; + const CppAD::AD adRho = adInput[numSpecies + 1]; + + + // 5. Call the actual templated function + // We let T9 and rho be constant, so we pass them as fixed values. + auto derivatives = calculateAllDerivatives>(adY, adT9, adRho); + + m_rhsADFun.Dependent(adInput, derivatives.dydt); + + LOG_INFO(m_logger, "AD tape recorded successfully for the RHS calculation. Number of independent variables: {}.", + adInput.size()); } } diff --git a/src/network/private/network.cpp b/src/network/private/network.cpp index 75b177a8..fba44f5d 100644 --- a/src/network/private/network.cpp +++ b/src/network/private/network.cpp @@ -33,7 +33,8 @@ namespace serif::network { m_config(serif::config::Config::getInstance()), m_logManager(serif::probe::LogManager::getInstance()), m_logger(m_logManager.getLogger("log")), - m_format(format) { + m_format(format), + m_constants(serif::constant::Constants::getInstance()){ if (format == NetworkFormat::UNKNOWN) { LOG_ERROR(m_logger, "nuclearNetwork::Network::Network() called with UNKNOWN format"); throw std::runtime_error("nuclearNetwork::Network::Network() called with UNKNOWN format"); @@ -73,6 +74,7 @@ namespace serif::network { serif::network::reaclib::REACLIBReactionSet build_reaclib_nuclear_network(const serif::composition::Composition &composition) { using namespace serif::network::reaclib; REACLIBReactionSet reactions; + auto logger = serif::probe::LogManager::getInstance().getLogger("log"); if (!s_initialized) { LOG_INFO(serif::probe::LogManager::getInstance().getLogger("log"), "REACLIB reactions not initialized. Calling initializeAllReaclibReactions()..."); @@ -89,6 +91,7 @@ namespace serif::network { } } if (gotReaction) { + LOG_INFO(logger, "Adding reaction {} to REACLIB reaction set.", reaction.id()); reactions.add_reaction(reaction); } } diff --git a/src/network/public/approx8.h b/src/network/public/approx8.h index a82d6fd3..b47bef47 100644 --- a/src/network/public/approx8.h +++ b/src/network/public/approx8.h @@ -306,13 +306,13 @@ namespace serif::network::approx8{ * @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); + 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 { return m_stiff; } + bool isStiff() const override { return m_stiff; } private: vector_type m_y; double m_tMax = 0; diff --git a/src/network/public/netgraph.h b/src/network/public/netgraph.h index 3ac3dcae..88252cf7 100644 --- a/src/network/public/netgraph.h +++ b/src/network/public/netgraph.h @@ -1,16 +1,20 @@ #pragma once -#include "network.h" -#include "reaclib.h" #include "atomicSpecies.h" #include "composition.h" +#include "network.h" +#include "reaclib.h" +#include #include #include -#include #include +#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. @@ -18,54 +22,248 @@ // REACLIBReactions are quite large data structures, so this could be a performance bottleneck. namespace serif::network { - template + // Define a concept to check if a type is one of our allowed scalar types + template + concept IsArithmeticOrAD = std::is_same_v || std::is_same_v>; + + 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 GraphNetwork final : public Network { public: explicit GraphNetwork(const serif::composition::Composition &composition); - explicit GraphNetwork(const serif::composition::Composition &composition, const double cullingThreshold, const double T9); + explicit GraphNetwork(const serif::composition::Composition &composition, + const double cullingThreshold, const double T9); NetOut evaluate(const NetIn &netIn) override; [[nodiscard]] const std::vector& getNetworkSpecies() const; [[nodiscard]] const reaclib::REACLIBReactionSet& getNetworkReactions() const; - [[nodiscard]] std::unordered_map getNetReactionStoichiometry(const reaclib::REACLIBReaction& reaction) const; + [[nodiscard]] std::unordered_map getNetReactionStoichiometry( + const reaclib::REACLIBReaction &reaction) const; [[nodiscard]] bool involvesSpecies(const serif::atomic::Species& species) const; - std::vector> detectCycles() const; - - std::vector topologicalSortReactions() const; + [[deprecated("not implimented")]] std::vector> detectCycles() const; + [[deprecated("not implimented")]] std::vector topologicalSortReactions() const; private: reaclib::REACLIBReactionSet m_reactions; ///< Set of REACLIB reactions for the network. - serif::composition::Composition m_initialComposition; - serif::composition::Composition m_currentComposition; - GeneralScalarType m_cullingThreshold; ///< Threshold for culling reactions. - GeneralScalarType m_T9; ///< Temperature in T9 units to evaluate culling. std::vector m_networkSpecies; ///< The species in the network. std::unordered_map m_networkSpeciesMap; std::unordered_map m_reactionIDMap; ///< Map of reaction IDs to REACLIB reactions. boost::numeric::ublas::compressed_matrix m_stoichiometryMatrix; ///< Stoichiometry matrix for the network. - boost::numeric::ublas::compressed_matrix m_jacobianMatrix; ///< Jacobian matrix for the network. + boost::numeric::ublas::compressed_matrix m_jacobianMatrix; ///< Jacobian matrix for the network. std::unordered_map m_speciesToIndexMap; ///< Map of species to their index in the stoichiometry matrix. + CppAD::ADFun m_rhsADFun; ///< AD tape + + struct ODETerm { + GraphNetwork& m_network; ///< Reference to the network + const double m_T9; + const double m_rho; + const size_t m_numSpecies; + + ODETerm(GraphNetwork& network, const double T9, double rho) : + m_network(network), m_T9(T9), m_rho(rho), m_numSpecies(m_network.m_networkSpecies.size()) {} + + void operator()(const boost::numeric::ublas::vector&Y, boost::numeric::ublas::vector& dYdt, double /* t */) { + const std::vector y(Y.begin(), Y.begin() + m_numSpecies); + + StepDerivatives derivatives = m_network.calculateAllDerivatives(y, m_T9, m_rho); + + dYdt.resize(m_numSpecies + 1); + std::copy(derivatives.dydt.begin(), derivatives.dydt.end(), dYdt.begin()); + dYdt(m_numSpecies) = derivatives.specificEnergyRate; // Last element is the specific energy rate + } + }; + + struct JacobianTerm { + GraphNetwork& m_network; + const double m_T9; + const double m_rho; + const size_t m_numSpecies; + + JacobianTerm(GraphNetwork& network, const double T9, const double rho) : + m_network(network), m_T9(T9), m_rho(rho), m_numSpecies(m_network.m_networkSpecies.size()) {} + + void operator()(const boost::numeric::ublas::vector& Y, boost::numeric::ublas::matrix& J, double /* t */, boost::numeric::ublas::vector& /* dfdt */) { + const std::vector 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; + } + } + } + }; + + template + struct StepDerivatives { + std::vector dydt; + T specificEnergyRate = T(0.0); + }; + private: void syncInternalMaps(); void collectNetworkSpecies(); void populateReactionIDMap(); void populateSpeciesToIndexMap(); - void buildNetworkGraph(); + void reserveJacobianMatrix(); + void recordADTape(); + + // --- Validation methods --- bool validateConservation() const; + void validateComposition(const serif::composition::Composition &composition, double culling, double T9); + + // --- Simple Derivative Calculations --- + template + StepDerivatives calculateAllDerivatives(const std::vector& Y, T T9, T rho) const; + // --- Generate the system matrices --- void generateStoichiometryMatrix(); - void generateJacobianMatrix(); + void generateJacobianMatrix(const std::vector& Y, const double T9, const double rho); - std::vector calculateRHS(const std::vector& Y, const GeneralScalarType T9, const GeneralScalarType rho) const; - GeneralScalarType calculateReactionRate(const reaclib::REACLIBReaction &reaction, const std::vector& Y, const GeneralScalarType T9, const GeneralScalarType rho) const; + template + std::vector calculateRHS(const std::vector &Y, const GeneralScalarType T9, + const GeneralScalarType rho) const; + template + GeneralScalarType calculateReactionRate(const reaclib::REACLIBReaction &reaction, + const std::vector &Y, const GeneralScalarType T9, + const GeneralScalarType rho) const; - void pruneNetworkGraph(); + void detectStiff(const NetIn &netIn, double T9, double numSpecies, const boost::numeric::ublas::vector& Y); }; + + + template + typename GraphNetwork::template StepDerivatives GraphNetwork::calculateAllDerivatives( + const std::vector &Y, T T9, T rho) const { + StepDerivatives result; + result.dydt.resize(m_networkSpecies.size(), static_cast(0.0)); + + if (rho < MIN_DENSITY_THRESHOLD) { + return result; // Return zero for everything if density is too low + } + + const T u = static_cast(m_constants.get("u").value); // Atomic mass unit in grams + const T MeV_to_erg = static_cast(m_constants.get("MeV_to_erg").value); + + T volumetricEnergyRate = static_cast(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(m_stoichiometryMatrix(speciesIndex, reactionIndex)); + if (nu_ij != 0) { + const T speciesAtomicMassAMU = static_cast(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(reaction.qValue()) * MeV_to_erg; + volumetricEnergyRate += reactionRate * q_value_ergs; + } + + result.specificEnergyRate = volumetricEnergyRate / rho; + + return result; + } + + + template + GeneralScalarType GraphNetwork::calculateReactionRate(const reaclib::REACLIBReaction &reaction, const std::vector &Y, + const GeneralScalarType T9, const GeneralScalarType rho) const { + const auto &constants = serif::constant::Constants::getInstance(); + + const auto u = constants.get("u"); // Atomic mass unit in g/mol + const GeneralScalarType uValue = static_cast(u.value); // Convert to double for calculations + const GeneralScalarType k_reaction = reaction.calculate_rate(T9); + + GeneralScalarType reactant_product_or_particle_densities = static_cast(1.0); + + std::unordered_map reactant_counts; + reactant_counts.reserve(reaction.reactants().size()); + for (const auto& reactant : reaction.reactants()) { + reactant_counts[std::string(reactant.name())]++; + } + + const GeneralScalarType minAbundanceThreshold = static_cast(MIN_ABUNDANCE_THRESHOLD); + const GeneralScalarType minDensityThreshold = static_cast(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(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(0.0); // If any reactant is below a threshold, return zero rate. + } + + GeneralScalarType atomicMassAMU = static_cast(m_networkSpecies[species_index].mass()); + + // Convert to number density + GeneralScalarType ni; + const GeneralScalarType denominator = atomicMassAMU * uValue; + if (denominator > minDensityThreshold) { + ni = (Yi * rho) / (denominator); + } else { + ni = static_cast(0.0); + } + + reactant_product_or_particle_densities *= ni; + + // Apply factorial correction for identical reactions + if (count > 1) { + reactant_product_or_particle_densities /= static_cast(std::tgamma(static_cast(count + 1))); // Gamma function for factorial + } + } + const GeneralScalarType Na = static_cast(constants.get("N_a").value); // Avogadro's number in mol^-1 + const int numReactants = static_cast(reaction.reactants().size()); + GeneralScalarType molarCorrectionFactor = static_cast(1.0); // No correction needed for single reactant reactions + if (numReactants > 1) { + molarCorrectionFactor = CppAD::pow(Na, static_cast(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 + std::vector GraphNetwork::calculateRHS( + const std::vector& Y, + const GeneralScalarType T9, + const GeneralScalarType rho) const { + + auto derivatives = calculateAllDerivatives(Y, T9, rho); + return derivatives.dydt; + } + }; \ No newline at end of file diff --git a/src/network/public/network.h b/src/network/public/network.h index 4c9ca58f..608bc4ca 100644 --- a/src/network/public/network.h +++ b/src/network/public/network.h @@ -29,6 +29,8 @@ #include "reaclib.h" #include +#include "const.h" + namespace serif::network { enum NetworkFormat { @@ -67,6 +69,7 @@ namespace serif::network { 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) }; /** @@ -87,6 +90,11 @@ namespace serif::network { serif::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; + } }; /** @@ -119,12 +127,18 @@ namespace serif::network { */ virtual NetOut evaluate(const NetIn &netIn); + virtual bool isStiff() const { return m_stiff; } + virtual void setStiff(const bool stiff) { m_stiff = stiff; } + protected: serif::config::Config& m_config; ///< Configuration instance serif::probe::LogManager& m_logManager; ///< Log manager instance quill::Logger* m_logger; ///< Logger instance NetworkFormat m_format; ///< Format of the network + serif::constant::Constants& m_constants; + + bool m_stiff = false; ///< Flag indicating if the network is stiff }; class ReaclibNetwork final : public Network { diff --git a/tests/network/approx8Test.cpp b/tests/network/approx8Test.cpp index fc722948..fafe9c99 100644 --- a/tests/network/approx8Test.cpp +++ b/tests/network/approx8Test.cpp @@ -6,7 +6,7 @@ #include "network.h" #include "composition.h" #include "reaclib.h" -#include "reactions.h" +#include "netgraph.h" #include @@ -74,8 +74,19 @@ TEST_F(approx8Test, reaclib) { serif::composition::Composition composition; composition.registerSymbol(symbols, true); composition.setMassFraction(symbols, comp); - bool isFinalized = composition.finalize(true); + composition.finalize(true); - reaclib::REACLIBReactionSet reactionSet = build_reaclib_nuclear_network(composition); - std::cout << reactionSet.size() << std::endl; + NetIn netIn; + netIn.composition = composition; + netIn.temperature = 1e7; + netIn.density = 1e2; + netIn.energy = 0.0; + + netIn.tMax = 3.15e17; + netIn.dt0 = 1e12; + + GraphNetwork network(composition); + NetOut netOut; + netOut = network.evaluate(netIn); + std::cout << netOut << std::endl; }