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

feat(weak): major weak rate progress

Browse Source 
Major weak rate progress which includes: A refactor of many of the public interfaces for GridFire Engines to use composition objects as opposed to raw abundance vectors. This helps prevent index mismatch errors. Further, the weak reaction class has been expanded with the majority of an implimentation, including an atomic_base derived class to allow for proper auto diff tracking of the interpolated table results. Some additional changes are that the version of fourdst and libcomposition have been bumped to versions with smarter caching of intermediate vectors and a few bug fixes.
This commit is contained in:
Emily Boudreaux
2025-10-07 15:16:03 -04:00
parent 4f1c260444
commit 8a0b5b2c36
53 changed files with 2310 additions and 1759 deletions
Download Patch File Download Diff File Expand all files Collapse all files

39
src/lib/engine/diagnostics/dynamic_engine_diagnostics.cpp
View File

@@ -5,7 +5,6 @@
#include <vector>
#include <string>
#include <iostream>
#include <iomanip>
#include <algorithm>
namespace gridfire::diagnostics {
@@ -73,7 +72,7 @@ namespace gridfire::diagnostics {
void inspect_species_balance(
const DynamicEngine& engine,
const std::string& species_name,
const std::vector<double>& Y_full,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) {
@@ -89,15 +88,15 @@ namespace gridfire::diagnostics {
const int stoichiometry = reaction->stoichiometry(species_obj);
if (stoichiometry == 0) continue;
const double flow = engine.calculateMolarReactionFlow(*reaction, Y_full, T9, rho);
const double flow = engine.calculateMolarReactionFlow(*reaction, comp, T9, rho);
if (stoichiometry > 0) {
creation_ids.push_back(std::string(reaction->id()));
creation_ids.emplace_back(reaction->id());
creation_stoichiometry.push_back(stoichiometry);
creation_flows.push_back(flow);
total_creation_flow += stoichiometry * flow;
} else {
destruction_ids.push_back(std::string(reaction->id()));
destruction_ids.emplace_back(reaction->id());
destruction_stoichiometry.push_back(stoichiometry);
destruction_flows.push_back(flow);
total_destruction_flow += std::abs(stoichiometry) * flow;
@@ -129,38 +128,38 @@ namespace gridfire::diagnostics {
void inspect_jacobian_stiffness(
const DynamicEngine& engine,
const std::vector<double>& Y_full,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) {
engine.generateJacobianMatrix(Y_full, T9, rho);
engine.generateJacobianMatrix(comp, T9, rho);
const auto& species_list = engine.getNetworkSpecies();
double max_diag = 0.0;
double max_off_diag = 0.0;
int max_diag_idx = -1;
int max_off_diag_i = -1, max_off_diag_j = -1;
std::optional<fourdst::atomic::Species> max_diag_species = std::nullopt;
std::optional<std::pair<fourdst::atomic::Species, fourdst::atomic::Species>> max_off_diag_species = std::nullopt;
for (size_t i = 0; i < species_list.size(); ++i) {
for (size_t j = 0; j < species_list.size(); ++j) {
const double val = std::abs(engine.getJacobianMatrixEntry(i, j));
if (i == j) {
if (val > max_diag) { max_diag = val; max_diag_idx = i; }
for (const auto& rowSpecies : species_list) {
for (const auto& colSpecies : species_list) {
const double val = std::abs(engine.getJacobianMatrixEntry(rowSpecies, colSpecies));
if (rowSpecies == colSpecies) {
if (val > max_diag) { max_diag = val; max_diag_species = colSpecies; }
} else {
if (val > max_off_diag) { max_off_diag = val; max_off_diag_i = i; max_off_diag_j = j; }
if (val > max_off_diag) { max_off_diag = val; max_off_diag_species = {rowSpecies, colSpecies};}
}
}
}
std::cout << "\n--- Jacobian Stiffness Report ---" << std::endl;
if (max_diag_idx != -1) {
if (max_diag_species.has_value()) {
std::cout << " Largest Diagonal Element (d(dYi/dt)/dYi): " << std::scientific << max_diag
<< " for species " << species_list[max_diag_idx].name() << std::endl;
<< " for species " << max_diag_species->name() << std::endl;
}
if (max_off_diag_i != -1) {
if (max_off_diag_species.has_value()) {
std::cout << " Largest Off-Diagonal Element (d(dYi/dt)/dYj): " << std::scientific << max_off_diag
<< " for d(" << species_list[max_off_diag_i].name()
<< ")/d(" << species_list[max_off_diag_j].name() << ")" << std::endl;
<< " for d(" << max_off_diag_species->first.name()
<< ")/d(" << max_off_diag_species->second.name() << ")" << std::endl;
}
std::cout << "---------------------------------" << std::endl;
}

201
src/lib/engine/engine_graph.cpp
View File

@@ -12,7 +12,6 @@
#include "quill/LogMacros.h"
#include <cstdint>
#include <iostream>
#include <set>
#include <stdexcept>
#include <string>
@@ -40,7 +39,8 @@ namespace gridfire {
const fourdst::composition::Composition &composition,
const partition::PartitionFunction& partitionFunction,
const BuildDepthType buildDepth) :
m_reactions(build_reaclib_nuclear_network(composition, buildDepth, false)),
m_weakRateInterpolator(rates::weak::UNIFIED_WEAK_DATA),
m_reactions(build_reaclib_nuclear_network(composition, m_weakRateInterpolator, buildDepth, false)),
m_depth(buildDepth),
m_partitionFunction(partitionFunction.clone())
{
@@ -50,15 +50,19 @@ namespace gridfire {
GraphEngine::GraphEngine(
const reaction::ReactionSet &reactions
) :
m_reactions(reactions) {
m_weakRateInterpolator(rates::weak::UNIFIED_WEAK_DATA),
m_reactions(reactions)
{
syncInternalMaps();
}
std::expected<StepDerivatives<double>, expectations::StaleEngineError> GraphEngine::calculateRHSAndEnergy(
const std::vector<double> &Y,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
const double Ye = comp.getElectronAbundance();
const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9;
if (m_usePrecomputation) {
std::vector<double> bare_rates;
std::vector<double> bare_reverse_rates;
@@ -66,33 +70,35 @@ namespace gridfire {
bare_reverse_rates.reserve(m_reactions.size());
// TODO: Add cache to this
for (const auto& reaction: m_reactions) {
bare_rates.push_back(reaction->calculate_rate(T9, rho, Y));
bare_reverse_rates.push_back(calculateReverseRate(*reaction, T9, rho, Y));
bare_rates.push_back(reaction->calculate_rate(T9, rho, Ye, mue, comp.getMolarAbundanceVector(), m_indexToSpeciesMap));
bare_reverse_rates.push_back(calculateReverseRate(*reaction, T9, rho, comp));
}
// --- The public facing interface can always use the precomputed version since taping is done internally ---
return calculateAllDerivativesUsingPrecomputation(Y, bare_rates, bare_reverse_rates, T9, rho);
return calculateAllDerivativesUsingPrecomputation(comp, bare_rates, bare_reverse_rates, T9, rho);
} else {
return calculateAllDerivatives<double>(Y, T9, rho);
return calculateAllDerivatives<double>(comp.getMolarAbundanceVector(), T9, rho, Ye, mue);
}
}
EnergyDerivatives GraphEngine::calculateEpsDerivatives(
const std::vector<double> &Y,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
const size_t numSpecies = m_networkSpecies.size();
const size_t numADInputs = numSpecies + 2; // +2 for T9 and rho
if (Y.size() != numSpecies) {
if (comp.getRegisteredSpecies().size() != numSpecies) {
LOG_ERROR(m_logger, "Input abundance vector size ({}) does not match number of species in the network ({}).",
Y.size(), numSpecies);
comp.getRegisteredSpecies().size(), numSpecies);
throw std::invalid_argument("Input abundance vector size does not match number of species in the network.");
}
std::vector<double> x(numADInputs);
const std::vector<double> Y = comp.getMolarAbundanceVector();
for (size_t i = 0; i < numSpecies; ++i) {
x[i] = Y[i];
}
@@ -148,6 +154,7 @@ namespace gridfire {
// --- Network Graph Construction Methods ---
void GraphEngine::collectNetworkSpecies() {
//TODO: Ensure consistent ordering in the m_networkSpecies vector so that it is ordered by species mass.
m_networkSpecies.clear();
m_networkSpeciesMap.clear();
@@ -190,6 +197,10 @@ namespace gridfire {
for (size_t i = 0; i < m_networkSpecies.size(); ++i) {
m_speciesToIndexMap.insert({m_networkSpecies[i], i});
}
m_indexToSpeciesMap.clear();
for (size_t i = 0; i < m_networkSpecies.size(); ++i) {
m_indexToSpeciesMap.insert({i, m_networkSpecies[i]});
}
}
void GraphEngine::reserveJacobianMatrix() const {
@@ -287,13 +298,18 @@ namespace gridfire {
const reaction::Reaction &reaction,
const double T9,
const double rho,
const std::vector<double> &Y
const fourdst::composition::Composition &comp
) const {
if (!m_useReverseReactions) {
LOG_TRACE_L3_LIMIT_EVERY_N(std::numeric_limits<int>::max(), m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id());
return 0.0; // If reverse reactions are not used, return 0.0
}
const double temp = T9 * 1e9; // Convert T9 to Kelvin
const double Ye = comp.getElectronAbundance();
// TODO: This is a dummy value for the electron chemical potential. We eventually need to replace this with an EOS call.
const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9;
// In debug builds we check the units on kB to ensure it is in erg/K. This is removed in release builds to avoid overhead. (Note assert is a no-op in release builds)
assert(Constants::getInstance().get("kB").unit == "erg / K");
@@ -301,7 +317,7 @@ namespace gridfire {
const double kBMeV = m_constants.kB * 624151; // Convert kB to MeV/K NOTE: This relies on the fact that m_constants.kB is in erg/K!
const double expFactor = std::exp(-reaction.qValue() / (kBMeV * temp));
double reverseRate = 0.0;
const double forwardRate = reaction.calculate_rate(T9, rho, Y);
const double forwardRate = reaction.calculate_rate(T9, rho, Ye, mue, comp.getMolarAbundanceVector(), m_indexToSpeciesMap);
if (reaction.reactants().size() == 2 && reaction.products().size() == 2) {
reverseRate = calculateReverseRateTwoBody(reaction, T9, forwardRate, expFactor);
@@ -393,14 +409,17 @@ namespace gridfire {
const reaction::Reaction &reaction,
const double T9,
const double rho,
const std::vector<double> &Y,
const fourdst::composition::Composition& comp,
const double reverseRate
) const {
if (!m_useReverseReactions) {
LOG_TRACE_L3_LIMIT_EVERY_N(std::numeric_limits<int>::max(), m_logger, "Reverse reactions are disabled. Returning 0.0 for reverse rate of reaction '{}'.", reaction.id());
return 0.0; // If reverse reactions are not used, return 0.0
}
const double d_log_kFwd = reaction.calculate_forward_rate_log_derivative(T9, rho, Y);
double Ye = comp.getElectronAbundance();
// TODO: This is a dummy value for the electron chemical potential. We eventually need to replace this with an EOS call.
double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9;
const double d_log_kFwd = reaction.calculate_forward_rate_log_derivative(T9, rho, Ye, mue, comp);
auto log_deriv_pf_op = [&](double acc, const auto& species) {
const double g = m_partitionFunction->evaluate(species.z(), species.a(), T9);
@@ -486,7 +505,7 @@ namespace gridfire {
void GraphEngine::rebuild(const fourdst::composition::Composition& comp, const BuildDepthType depth) {
if (depth != m_depth) {
m_depth = depth;
m_reactions = build_reaclib_nuclear_network(comp, m_depth, false);
m_reactions = build_reaclib_nuclear_network(comp, m_weakRateInterpolator, m_depth, false);
syncInternalMaps(); // Resync internal maps after changing the depth
} else {
LOG_DEBUG(m_logger, "Rebuild requested with the same depth. No changes made to the network.");
@@ -494,7 +513,7 @@ namespace gridfire {
}
StepDerivatives<double> GraphEngine::calculateAllDerivativesUsingPrecomputation(
const std::vector<double> &Y_in,
const fourdst::composition::Composition& comp,
const std::vector<double> &bare_rates,
const std::vector<double> &bare_reverse_rates,
const double T9,
@@ -504,33 +523,27 @@ namespace gridfire {
const std::vector<double> screeningFactors = m_screeningModel->calculateScreeningFactors(
m_reactions,
m_networkSpecies,
Y_in,
comp.getMolarAbundanceVector(),
T9,
rho
);
// TODO: Fix up the precomputation to use the new comp in interface as opposed to a raw vector of molar abundances
// This will require carefully checking the way the precomputation is stashed.
// --- Optimized loop ---
std::vector<double> molarReactionFlows;
molarReactionFlows.reserve(m_precomputedReactions.size());
for (const auto& precomp : m_precomputedReactions) {
double forwardAbundanceProduct = 1.0;
// bool below_threshold = false;
for (size_t i = 0; i < precomp.unique_reactant_indices.size(); ++i) {
const size_t reactantIndex = precomp.unique_reactant_indices[i];
const fourdst::atomic::Species& reactant = m_networkSpecies[reactantIndex];
const int power = precomp.reactant_powers[i];
// const double abundance = Y_in[reactantIndex];
// if (abundance < MIN_ABUNDANCE_THRESHOLD) {
// below_threshold = true;
// break;
// }
forwardAbundanceProduct *= std::pow(Y_in[reactantIndex], power);
forwardAbundanceProduct *= std::pow(comp.getMolarAbundance(reactant), power);
}
// if (below_threshold) {
// molarReactionFlows.push_back(0.0);
// continue; // Skip this reaction if any reactant is below the abundance threshold
// }
const double bare_rate = bare_rates[precomp.reaction_index];
const double screeningFactor = screeningFactors[precomp.reaction_index];
@@ -549,7 +562,9 @@ namespace gridfire {
const double bare_reverse_rate = bare_reverse_rates[precomp.reaction_index];
double reverseAbundanceProduct = 1.0;
for (size_t i = 0; i < precomp.unique_product_indices.size(); ++i) {
reverseAbundanceProduct *= std::pow(Y_in[precomp.unique_product_indices[i]], precomp.product_powers[i]);
const size_t productIndex = precomp.unique_product_indices[i];
const fourdst::atomic::Species& product = m_networkSpecies[productIndex];
reverseAbundanceProduct *= std::pow(comp.getMolarAbundance(product), precomp.product_powers[i]);
}
reverseMolarReactionFlow = screeningFactor *
bare_reverse_rate *
@@ -564,25 +579,28 @@ namespace gridfire {
// --- Assemble molar abundance derivatives ---
StepDerivatives<double> result;
result.dydt.assign(m_networkSpecies.size(), 0.0); // Initialize derivatives to zero
for (const auto& species: m_networkSpecies) {
result.dydt[species] = 0.0; // Initialize the change in abundance for each network species to 0
}
for (size_t j = 0; j < m_precomputedReactions.size(); ++j) {
const auto& precomp = m_precomputedReactions[j];
const double R_j = molarReactionFlows[j];
for (size_t i = 0; i < precomp.affected_species_indices.size(); ++i) {
const size_t speciesIndex = precomp.affected_species_indices[i];
const fourdst::atomic::Species& species = m_networkSpecies[speciesIndex];
const int stoichiometricCoefficient = precomp.stoichiometric_coefficients[i];
// Update the derivative for this species
result.dydt[speciesIndex] += static_cast<double>(stoichiometricCoefficient) * R_j;
result.dydt.at(species) += static_cast<double>(stoichiometricCoefficient) * R_j;
}
}
// --- Calculate the nuclear energy generation rate ---
double massProductionRate = 0.0; // [mol][s^-1]
for (size_t i = 0; i < m_networkSpecies.size(); ++i) {
const auto& species = m_networkSpecies[i];
massProductionRate += result.dydt[i] * species.mass() * m_constants.u;
for (const auto & species : m_networkSpecies) {
massProductionRate += result.dydt[species] * species.mass() * m_constants.u;
}
result.nuclearEnergyGenerationRate = -massProductionRate * m_constants.Na * m_constants.c * m_constants.c; // [erg][s^-1][g^-1]
return result;
@@ -631,21 +649,22 @@ namespace gridfire {
m_stoichiometryMatrix.nnz()); // Assuming nnz() exists for compressed_matrix
}
StepDerivatives<double> GraphEngine::calculateAllDerivatives(
const std::vector<double> &Y_in,
const double T9,
const double rho
) const {
return calculateAllDerivatives<double>(Y_in, T9, rho);
}
StepDerivatives<ADDouble> GraphEngine::calculateAllDerivatives(
const std::vector<ADDouble> &Y_in,
const ADDouble &T9,
const ADDouble &rho
) const {
return calculateAllDerivatives<ADDouble>(Y_in, T9, rho);
}
// StepDerivatives<double> GraphEngine::calculateAllDerivatives(
// const std::vector<double> &Y,
// const double T9,
// const double rho
// ) const {
// return calculateAllDerivatives<double>(Y, T9, rho);
// }
//
// StepDerivatives<ADDouble> GraphEngine::calculateAllDerivatives(
// const std::vector<ADDouble> &Y,
// ADDouble T9,
// ADDouble rho
// ) const {
// //TODO: Sort out why this template is not being found
// return calculateAllDerivatives<ADDouble>(Y, T9, rho);
// }
void GraphEngine::setScreeningModel(const screening::ScreeningType model) {
m_screeningModel = screening::selectScreeningModel(model);
@@ -670,15 +689,21 @@ namespace gridfire {
double GraphEngine::calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const std::vector<double> &Y,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
return calculateMolarReactionFlow<double>(reaction, Y, T9, rho);
const double Ye = comp.getElectronAbundance();
// TODO: This is a dummy placeholder which must be replaced with an EOS call
const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9;
return calculateMolarReactionFlow<double>(reaction, comp.getMolarAbundanceVector(), T9, rho, Ye, mue);
}
void GraphEngine::generateJacobianMatrix(
const std::vector<double> &Y_dynamic,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
@@ -688,6 +713,7 @@ namespace gridfire {
// 1. Pack the input variables into a vector for CppAD
std::vector<double> adInput(numSpecies + 2, 0.0); // +2 for T9 and rho
const std::vector<double>& Y_dynamic = comp.getMolarAbundanceVector();
for (size_t i = 0; i < numSpecies; ++i) {
adInput[i] = std::max(Y_dynamic[i], 1e-99); // regularize the jacobian...
}
@@ -711,7 +737,7 @@ namespace gridfire {
}
void GraphEngine::generateJacobianMatrix(
const std::vector<double> &Y_dynamic,
const fourdst::composition::Composition &comp,
const double T9,
const double rho,
const SparsityPattern &sparsityPattern
@@ -719,6 +745,7 @@ namespace gridfire {
// --- Pack the input variables into a vector for CppAD ---
const size_t numSpecies = m_networkSpecies.size();
std::vector<double> x(numSpecies + 2, 0.0);
const std::vector<double>& Y_dynamic = comp.getMolarAbundanceVector();
for (size_t i = 0; i < numSpecies; ++i) {
x[i] = Y_dynamic[i];
}
@@ -767,8 +794,13 @@ namespace gridfire {
}
}
double GraphEngine::getJacobianMatrixEntry(const int i, const int j) const {
// LOG_TRACE_L3(m_logger, "Getting jacobian matrix entry for {},{} = {}", i, j, m_jacobianMatrix(i, j));
double GraphEngine::getJacobianMatrixEntry(
const fourdst::atomic::Species& rowSpecies,
const fourdst::atomic::Species& colSpecies
) const {
//PERF: There may be some way to make this more efficient
const size_t i = getSpeciesIndex(rowSpecies);
const size_t j = getSpeciesIndex(colSpecies);
return m_jacobianMatrix(i, j);
}
@@ -779,10 +811,10 @@ namespace gridfire {
}
int GraphEngine::getStoichiometryMatrixEntry(
const int speciesIndex,
const int reactionIndex
const fourdst::atomic::Species& species,
const reaction::Reaction &reaction
) const {
return m_stoichiometryMatrix(speciesIndex, reactionIndex);
return reaction.stoichiometry(species);
}
void GraphEngine::exportToDot(const std::string &filename) const {
@@ -871,18 +903,21 @@ namespace gridfire {
}
std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> GraphEngine::getSpeciesTimescales(
const std::vector<double> &Y,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
auto [dydt, _] = calculateAllDerivatives<double>(Y, T9, rho);
const double Ye = comp.getElectronAbundance();
// TODO: This is a dummy placeholder which must be replaced with an EOS call
const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9;
auto [dydt, _] = calculateAllDerivatives<double>(comp.getMolarAbundanceVector(), T9, rho, Ye, mue);
std::unordered_map<fourdst::atomic::Species, double> speciesTimescales;
speciesTimescales.reserve(m_networkSpecies.size());
for (size_t i = 0; i < m_networkSpecies.size(); ++i) {
for (const auto& species : m_networkSpecies) {
double timescale = std::numeric_limits<double>::infinity();
const auto species = m_networkSpecies[i];
if (std::abs(dydt[i]) > 0.0) {
timescale = std::abs(Y[i] / dydt[i]);
if (std::abs(dydt.at(species)) > 0.0) {
timescale = std::abs(comp.getMolarAbundance(species) / dydt.at(species));
}
speciesTimescales.emplace(species, timescale);
}
@@ -890,24 +925,28 @@ namespace gridfire {
}
std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError> GraphEngine::getSpeciesDestructionTimescales(
const std::vector<double> &Y,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
auto [dydt, _] = calculateAllDerivatives<double>(Y, T9, rho);
const double Ye = comp.getElectronAbundance();
// TODO: This is a dummy placeholder which must be replaced with an EOS call
const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9;
auto [dydt, _] = calculateAllDerivatives<double>(comp.getMolarAbundanceVector(), T9, rho, Ye, mue);
std::unordered_map<fourdst::atomic::Species, double> speciesDestructionTimescales;
speciesDestructionTimescales.reserve(m_networkSpecies.size());
for (const auto& species : m_networkSpecies) {
double netDestructionFlow = 0.0;
for (const auto& reaction : m_reactions) {
if (reaction->stoichiometry(species) < 0) {
const auto flow = calculateMolarReactionFlow<double>(*reaction, Y, T9, rho);
const auto flow = calculateMolarReactionFlow<double>(*reaction, comp.getMolarAbundanceVector(), T9, rho, Ye, mue);
netDestructionFlow += flow;
}
}
double timescale = std::numeric_limits<double>::infinity();
if (netDestructionFlow != 0.0) {
timescale = Y[getSpeciesIndex(species)] / netDestructionFlow;
timescale = comp.getMolarAbundance(species) / netDestructionFlow;
}
speciesDestructionTimescales.emplace(species, timescale);
}
@@ -964,12 +1003,21 @@ namespace gridfire {
const CppAD::AD<double> adT9 = adInput[numSpecies];
const CppAD::AD<double> adRho = adInput[numSpecies + 1];
// Dummy values for Ye and mue to let taping happen
const CppAD::AD<double> adYe = 1.0;
const CppAD::AD<double> adMue = 1.0;
// 5. Call the actual templated function
// We let T9 and rho be constant, so we pass them as fixed values.
auto [dydt, nuclearEnergyGenerationRate] = calculateAllDerivatives<CppAD::AD<double>>(adY, adT9, adRho);
auto [dydt, nuclearEnergyGenerationRate] = calculateAllDerivatives<CppAD::AD<double>>(adY, adT9, adRho, adYe, adMue);
m_rhsADFun.Dependent(adInput, dydt);
// Extract the raw vector from the associative map
std::vector<CppAD::AD<double>> dydt_vec;
dydt_vec.reserve(dydt.size());
std::ranges::transform(dydt, std::back_inserter(dydt_vec),[](const auto& kv) { return kv.second; });
m_rhsADFun.Dependent(adInput, dydt_vec);
LOG_TRACE_L1(m_logger, "AD tape recorded successfully for the RHS calculation. Number of independent variables: {}.",
adInput.size());
@@ -996,16 +1044,19 @@ namespace gridfire {
CppAD::Independent(adInput);
const std::vector<CppAD::AD<double>> adY(adInput.begin(), adInput.begin() + numSpecies);
const std::vector<CppAD::AD<double>> adY(adInput.begin(), adInput.begin() + static_cast<long>(numSpecies));
const CppAD::AD<double> adT9 = adInput[numSpecies];
const CppAD::AD<double> adRho = adInput[numSpecies + 1];
StepDerivatives<CppAD::AD<double>> derivatives = calculateAllDerivatives<CppAD::AD<double>>(adY, adT9, adRho);
// Dummy values for Ye and mue to let taping happen
const CppAD::AD<double> adYe = 1.0;
const CppAD::AD<double> adMue = 1.0;
auto [dydt, nuclearEnergyGenerationRate] = calculateAllDerivatives<CppAD::AD<double>>(adY, adT9, adRho, adYe, adMue);
std::vector<CppAD::AD<double>> adOutput(1);
adOutput[0] = derivatives.nuclearEnergyGenerationRate;
adOutput[0] = nuclearEnergyGenerationRate;
m_epsADFun.Dependent(adInput, adOutput);
// m_epsADFun.optimize();
LOG_TRACE_L1(m_logger, "AD tape recorded successfully for the EPS calculation. Number of independent variables: {}.", adInput.size());

106
src/lib/engine/procedures/construction.cpp
View File

@@ -1,7 +1,11 @@
#include "gridfire/engine/procedures/construction.h"
#include "gridfire/reaction/weak/weak_interpolator.h"
#include "gridfire/reaction/weak/weak.h"
#include "gridfire/reaction/weak/weak_types.h"
#include <ranges>
#include <stdexcept>
#include <memory>
#include "gridfire/reaction/reaction.h"
#include "gridfire/reaction/reaclib.h"
@@ -20,11 +24,11 @@ namespace gridfire {
using fourdst::atomic::Species;
ReactionSet build_reaclib_nuclear_network(
const Composition &composition,
BuildDepthType maxLayers,
bool reverse
) {
ReactionSet build_nuclear_network(
const Composition& composition,
const rates::weak::WeakRateInterpolator& weak_interpolator,
BuildDepthType maxLayers,
bool reverse_reaclib) {
int depth;
if (std::holds_alternative<NetworkBuildDepth>(maxLayers)) {
depth = static_cast<int>(std::get<NetworkBuildDepth>(maxLayers));
@@ -37,31 +41,71 @@ namespace gridfire {
throw std::logic_error("Network build depth is set to 0. No reactions will be collected.");
}
const auto& allReactions = reaclib::get_all_reaclib_reactions();
std::vector<Reaction*> remainingReactions;
for (const auto& reaction : allReactions) {
if (reaction->is_reverse() == reverse) {
remainingReactions.push_back(reaction.get());
// --- Step 1: Create a single master pool that owns all possible reactions ---
ReactionSet master_reaction_pool;
// Clone all relevant REACLIB reactions into the master pool
const auto& allReaclibReactions = reaclib::get_all_reaclib_reactions();
for (const auto& reaction : allReaclibReactions) {
if (reaction->is_reverse() == reverse_reaclib) {
master_reaction_pool.add_reaction(reaction->clone());
}
}
if (depth == static_cast<int>(NetworkBuildDepth::Full)) {
LOG_INFO(logger, "Building full nuclear network with a total of {} reactions.", allReactions.size());
const ReactionSet reactionSet(remainingReactions);
return reactionSet;
for (const auto& parent_species: weak_interpolator.available_isotopes()) {
master_reaction_pool.add_reaction(
std::make_unique<rates::weak::WeakReaction>(
parent_species,
rates::weak::WeakReactionType::BETA_PLUS_DECAY,
weak_interpolator
)
);
master_reaction_pool.add_reaction(
std::make_unique<rates::weak::WeakReaction>(
parent_species,
rates::weak::WeakReactionType::BETA_MINUS_DECAY,
weak_interpolator
)
);
master_reaction_pool.add_reaction(
std::make_unique<rates::weak::WeakReaction>(
parent_species,
rates::weak::WeakReactionType::ELECTRON_CAPTURE,
weak_interpolator
)
);
master_reaction_pool.add_reaction(
std::make_unique<rates::weak::WeakReaction>(
parent_species,
rates::weak::WeakReactionType::POSITRON_CAPTURE,
weak_interpolator
)
);
}
// --- Step 2: Use non-owning raw pointers for the fast build algorithm ---
std::vector<Reaction*> remainingReactions;
remainingReactions.reserve(master_reaction_pool.size());
for(const auto& reaction : master_reaction_pool) {
remainingReactions.push_back(reaction.get());
}
if (depth == static_cast<int>(NetworkBuildDepth::Full)) {
LOG_INFO(logger, "Building full nuclear network with a total of {} reactions.", remainingReactions.size());
return master_reaction_pool; // No need to build layer by layer if we want the full network
}
// --- Step 3: Execute the layered network build using observing pointers ---
std::unordered_set<Species> availableSpecies;
for (const auto &entry: composition | std::views::values) {
for (const auto& entry : composition | std::views::values) {
if (entry.mass_fraction() > 0.0) {
availableSpecies.insert(entry.isotope());
}
}
std::vector<Reaction*> collectedReactions;
std::vector<Reaction*> collectedReactionPtrs;
LOG_INFO(logger, "Starting network construction with {} available species.", availableSpecies.size());
for (int layer = 0; layer < depth && !remainingReactions.empty(); ++layer) {
LOG_TRACE_L1(logger, "Collecting reactions for layer {} with {} remaining reactions. Currently there are {} available species", layer, remainingReactions.size(), availableSpecies.size());
std::vector<Reaction*> reactionsForNextPass;
@@ -70,7 +114,7 @@ namespace gridfire {
reactionsForNextPass.reserve(remainingReactions.size());
for (const auto &reaction : remainingReactions) {
for (Reaction* reaction : remainingReactions) {
bool allReactantsAvailable = true;
for (const auto& reactant : reaction->reactants()) {
if (!availableSpecies.contains(reactant)) {
@@ -80,7 +124,7 @@ namespace gridfire {
}
if (allReactantsAvailable) {
collectedReactions.push_back(reaction);
collectedReactionPtrs.push_back(reaction);
newReactionsAdded = true;
for (const auto& product : reaction->products()) {
@@ -92,19 +136,29 @@ namespace gridfire {
}
if (!newReactionsAdded) {
LOG_INFO(logger, "No new reactions added in layer {}. Stopping network construction with {} reactions collected.", layer, collectedReactions.size());
LOG_INFO(logger, "No new reactions added in layer {}. Stopping network construction.", layer);
break;
}
LOG_TRACE_L1(logger, "Layer {}: Collected {} reactions. New products this layer: {}", layer, collectedReactions.size(), newProductsThisLayer.size());
LOG_TRACE_L1(logger, "Layer {}: Collected {} new reactions. New products this layer: {}", collectedReactionPtrs.size() - collectedReactionPtrs.size(), newProductsThisLayer.size());
availableSpecies.insert(newProductsThisLayer.begin(), newProductsThisLayer.end());
remainingReactions = std::move(reactionsForNextPass);
}
LOG_INFO(logger, "Network construction completed with {} reactions collected.", collectedReactions.size());
const ReactionSet reactionSet(collectedReactions);
return reactionSet;
// --- Step 4: Construct the final ReactionSet by moving ownership ---
LOG_INFO(logger, "Network construction completed. Assembling final set of {} reactions.", collectedReactionPtrs.size());
ReactionSet finalReactionSet;
std::unordered_set<Reaction*> collectedPtrSet(collectedReactionPtrs.begin(), collectedReactionPtrs.end());
for (auto& reaction_ptr : master_reaction_pool) {
if (reaction_ptr && collectedPtrSet.contains(reaction_ptr.get())) {
finalReactionSet.add_reaction(std::move(reaction_ptr));
}
}
return finalReactionSet;
}
}

32
src/lib/engine/procedures/priming.cpp
View File

@@ -17,7 +17,7 @@ namespace gridfire {
const reaction::Reaction* findDominantCreationChannel (
const DynamicEngine& engine,
const Species& species,
const std::vector<double>& Y,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) {
@@ -25,7 +25,7 @@ namespace gridfire {
double maxFlow = -1.0;
for (const auto& reaction : engine.getNetworkReactions()) {
if (reaction->contains(species) && reaction->stoichiometry(species) > 0) {
const double flow = engine.calculateMolarReactionFlow(*reaction, Y, T9, rho);
const double flow = engine.calculateMolarReactionFlow(*reaction, comp, T9, rho);
if (flow > maxFlow) {
maxFlow = flow;
dominateReaction = reaction.get();
@@ -67,7 +67,7 @@ namespace gridfire {
* robustly primed composition.
*/
PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) {
auto logger = LogManager::getInstance().getLogger("log");
auto logger = fourdst::logging::LogManager::getInstance().getLogger("log");
// --- Initial Setup ---
// Identify all species with zero initial mass fraction that need to be primed.
@@ -127,9 +127,6 @@ namespace gridfire {
}
tempComp.finalize(true);
NetIn tempNetIn = netIn;
tempNetIn.composition = tempComp;
NetworkPrimingEngineView primer(primingSpecies, engine);
if (primer.getNetworkReactions().size() == 0) {
LOG_ERROR(logger, "No priming reactions found for species {}.", primingSpecies.name());
@@ -138,15 +135,14 @@ namespace gridfire {
continue;
}
const auto Y = primer.mapNetInToMolarAbundanceVector(tempNetIn);
const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho);
const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, tempComp, T9, rho);
if (destructionRateConstant > 1e-99) {
const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho);
const double creationRate = calculateCreationRate(primer, primingSpecies, tempComp, T9, rho);
double equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass();
if (std::isnan(equilibriumMassFraction)) equilibriumMassFraction = 0.0;
if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho)) {
if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, tempComp, T9, rho)) {
// Store the request instead of applying it immediately.
requests.push_back({primingSpecies, equilibriumMassFraction, dominantChannel->reactants()});
} else {
@@ -403,19 +399,17 @@ namespace gridfire {
double calculateDestructionRateConstant(
const DynamicEngine& engine,
const fourdst::atomic::Species& species,
const std::vector<double>& Y,
const Species& species,
const Composition& comp,
const double T9,
const double rho
) {
const size_t speciesIndex = engine.getSpeciesIndex(species);
std::vector<double> Y_scaled(Y.begin(), Y.end());
Y_scaled[speciesIndex] = 1.0; // Set the abundance of the species to 1.0 for rate constant calculation
//TODO: previously (when using raw vectors) I set y[speciesIndex] = 1.0 to let there be enough so that just the destruction rate could be found (without bottlenecks from abundance) we will need to do a similar thing here.
double destructionRateConstant = 0.0;
for (const auto& reaction: engine.getNetworkReactions()) {
if (reaction->contains_reactant(species)) {
const int stoichiometry = reaction->stoichiometry(species);
destructionRateConstant += std::abs(stoichiometry) * engine.calculateMolarReactionFlow(*reaction, Y_scaled, T9, rho);
destructionRateConstant += std::abs(stoichiometry) * engine.calculateMolarReactionFlow(*reaction, comp, T9, rho);
}
}
return destructionRateConstant;
@@ -424,7 +418,7 @@ namespace gridfire {
double calculateCreationRate(
const DynamicEngine& engine,
const Species& species,
const std::vector<double>& Y,
const Composition& comp,
const double T9,
const double rho
) {
@@ -432,9 +426,9 @@ namespace gridfire {
for (const auto& reaction: engine.getNetworkReactions()) {
const int stoichiometry = reaction->stoichiometry(species);
if (stoichiometry > 0) {
if (engine.calculateMolarReactionFlow(*reaction, Y, T9, rho) > 0.0) {
if (engine.calculateMolarReactionFlow(*reaction, comp, T9, rho) > 0.0) {
}
creationRate += stoichiometry * engine.calculateMolarReactionFlow(*reaction, Y, T9, rho);
creationRate += stoichiometry * engine.calculateMolarReactionFlow(*reaction, comp, T9, rho);
}
}
return creationRate;

119
src/lib/engine/views/engine_adaptive.cpp
View File

@@ -95,8 +95,7 @@ namespace gridfire {
LOG_TRACE_L1(m_logger, "Updating AdaptiveEngineView with new network input...");
std::vector<double> Y_Full;
std::vector<ReactionFlow> allFlows = calculateAllReactionFlows(updatedNetIn, Y_Full);
auto [allFlows, composition] = calculateAllReactionFlows(updatedNetIn);
double maxFlow = 0.0;
@@ -110,11 +109,11 @@ namespace gridfire {
const std::unordered_set<Species> reachableSpecies = findReachableSpecies(updatedNetIn);
LOG_DEBUG(m_logger, "Found {} reachable species in adaptive engine view.", reachableSpecies.size());
const std::vector<const reaction::Reaction*> finalReactions = cullReactionsByFlow(allFlows, reachableSpecies, Y_Full, maxFlow);
const std::vector<const reaction::Reaction*> finalReactions = cullReactionsByFlow(allFlows, reachableSpecies, composition, maxFlow);
finalizeActiveSet(finalReactions);
auto [rescuedReactions, rescuedSpecies] = rescueEdgeSpeciesDestructionChannel(Y_Full, netIn.temperature/1e9, netIn.density, m_activeSpecies, m_activeReactions);
auto [rescuedReactions, rescuedSpecies] = rescueEdgeSpeciesDestructionChannel(composition, netIn.temperature/1e9, netIn.density, m_activeSpecies, m_activeReactions);
for (const auto& reactionPtr : rescuedReactions) {
m_activeReactions.add_reaction(*reactionPtr);
@@ -145,59 +144,46 @@ namespace gridfire {
}
std::expected<StepDerivatives<double>, expectations::StaleEngineError> AdaptiveEngineView::calculateRHSAndEnergy(
const std::vector<double> &Y_culled,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
validateState();
const auto Y_full = mapCulledToFull(Y_culled);
auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho);
// TODO: Think about if I need to reach in and adjust the composition to zero out inactive species.
auto result = m_baseEngine.calculateRHSAndEnergy(comp, T9, rho);
if (!result) {
return std::unexpected{result.error()};
}
const auto [dydt, nuclearEnergyGenerationRate] = result.value();
StepDerivatives<double> culledResults;
culledResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate;
culledResults.dydt = mapFullToCulled(dydt);
return culledResults;
return result.value();
}
EnergyDerivatives AdaptiveEngineView::calculateEpsDerivatives(
const std::vector<double> &Y_culled,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
validateState();
const auto Y_full = mapCulledToFull(Y_culled);
return m_baseEngine.calculateEpsDerivatives(Y_full, T9, rho);
return m_baseEngine.calculateEpsDerivatives(comp, T9, rho);
}
void AdaptiveEngineView::generateJacobianMatrix(
const std::vector<double> &Y_dynamic,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
validateState();
const auto Y_full = mapCulledToFull(Y_dynamic);
m_baseEngine.generateJacobianMatrix(Y_full, T9, rho);
m_baseEngine.generateJacobianMatrix(comp, T9, rho);
}
double AdaptiveEngineView::getJacobianMatrixEntry(
const int i_culled,
const int j_culled
const Species &rowSpecies,
const Species &colSpecies
) const {
validateState();
const size_t i_full = mapCulledToFullSpeciesIndex(i_culled);
const size_t j_full = mapCulledToFullSpeciesIndex(j_culled);
return m_baseEngine.getJacobianMatrixEntry(static_cast<int>(i_full), static_cast<int>(j_full));
return m_baseEngine.getJacobianMatrixEntry(rowSpecies, colSpecies);
}
void AdaptiveEngineView::generateStoichiometryMatrix() {
@@ -206,18 +192,16 @@ namespace gridfire {
}
int AdaptiveEngineView::getStoichiometryMatrixEntry(
const int speciesIndex_culled,
const int reactionIndex_culled
const Species &species,
const reaction::Reaction& reaction
) const {
validateState();
const size_t speciesIndex_full = mapCulledToFullSpeciesIndex(speciesIndex_culled);
const size_t reactionIndex_full = mapCulledToFullReactionIndex(reactionIndex_culled);
return m_baseEngine.getStoichiometryMatrixEntry(static_cast<int>(speciesIndex_full), static_cast<int>(reactionIndex_full));
return m_baseEngine.getStoichiometryMatrixEntry(species, reaction);
}
double AdaptiveEngineView::calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const std::vector<double> &Y_culled,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
@@ -227,9 +211,8 @@ namespace gridfire {
m_logger -> flush_log();
throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id()));
}
const auto Y = mapCulledToFull(Y_culled);
return m_baseEngine.calculateMolarReactionFlow(reaction, Y, T9, rho);
return m_baseEngine.calculateMolarReactionFlow(reaction, comp, T9, rho);
}
const reaction::ReactionSet & AdaptiveEngineView::getNetworkReactions() const {
@@ -242,13 +225,12 @@ namespace gridfire {
}
std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError> AdaptiveEngineView::getSpeciesTimescales(
const std::vector<double> &Y_culled,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
validateState();
const auto Y_full = mapCulledToFull(Y_culled);
const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
const auto result = m_baseEngine.getSpeciesTimescales(comp, T9, rho);
if (!result) {
return std::unexpected{result.error()};
@@ -270,15 +252,13 @@ namespace gridfire {
std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError>
AdaptiveEngineView::getSpeciesDestructionTimescales(
const std::vector<double> &Y,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
validateState();
const std::vector<double> Y_full = mapCulledToFull(Y);
const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho);
const auto result = m_baseEngine.getSpeciesDestructionTimescales(comp, T9, rho);
if (!result) {
return std::unexpected{result.error()};
}
@@ -344,7 +324,7 @@ namespace gridfire {
}
size_t AdaptiveEngineView::mapCulledToFullSpeciesIndex(size_t culledSpeciesIndex) const {
if (culledSpeciesIndex < 0 || culledSpeciesIndex >= m_speciesIndexMap.size()) {
if (culledSpeciesIndex >= m_speciesIndexMap.size()) {
LOG_ERROR(m_logger, "Culled index {} is out of bounds for species index map of size {}.", culledSpeciesIndex, m_speciesIndexMap.size());
m_logger->flush_log();
throw std::out_of_range("Culled index " + std::to_string(culledSpeciesIndex) + " is out of bounds for species index map of size " + std::to_string(m_speciesIndexMap.size()) + ".");
@@ -353,7 +333,7 @@ namespace gridfire {
}
size_t AdaptiveEngineView::mapCulledToFullReactionIndex(size_t culledReactionIndex) const {
if (culledReactionIndex < 0 || culledReactionIndex >= m_reactionIndexMap.size()) {
if (culledReactionIndex >= m_reactionIndexMap.size()) {
LOG_ERROR(m_logger, "Culled index {} is out of bounds for reaction index map of size {}.", culledReactionIndex, m_reactionIndexMap.size());
m_logger->flush_log();
throw std::out_of_range("Culled index " + std::to_string(culledReactionIndex) + " is out of bounds for reaction index map of size " + std::to_string(m_reactionIndexMap.size()) + ".");
@@ -369,21 +349,17 @@ namespace gridfire {
}
}
// TODO: Change this to use a return value instead of an output parameter.
std::vector<AdaptiveEngineView::ReactionFlow> AdaptiveEngineView::calculateAllReactionFlows(
const NetIn &netIn,
std::vector<double> &out_Y_Full
std::pair<std::vector<AdaptiveEngineView::ReactionFlow>, fourdst::composition::Composition> AdaptiveEngineView::calculateAllReactionFlows(
const NetIn &netIn
) const {
const auto& fullSpeciesList = m_baseEngine.getNetworkSpecies();
out_Y_Full.clear();
out_Y_Full.reserve(fullSpeciesList.size());
fourdst::composition::Composition composition = netIn.composition;
for (const auto& species: fullSpeciesList) {
if (netIn.composition.contains(species)) {
out_Y_Full.push_back(netIn.composition.getMolarAbundance(std::string(species.name())));
} else {
if (!netIn.composition.contains(species)) {
LOG_TRACE_L2(m_logger, "Species '{}' not found in composition. Setting abundance to 0.0.", species.name());
out_Y_Full.push_back(0.0);
composition.registerSpecies(species);
composition.setMassFraction(species, 0.0);
}
}
@@ -394,11 +370,11 @@ namespace gridfire {
const auto& fullReactionSet = m_baseEngine.getNetworkReactions();
reactionFlows.reserve(fullReactionSet.size());
for (const auto& reaction : fullReactionSet) {
const double flow = m_baseEngine.calculateMolarReactionFlow(*reaction, out_Y_Full, T9, rho);
const double flow = m_baseEngine.calculateMolarReactionFlow(*reaction, composition, T9, rho);
reactionFlows.push_back({reaction.get(), flow});
LOG_TRACE_L1(m_logger, "Reaction '{}' has flow rate: {:0.3E} [mol/s/g]", reaction->id(), flow);
}
return reactionFlows;
return {reactionFlows, composition};
}
std::unordered_set<Species> AdaptiveEngineView::findReachableSpecies(
@@ -447,7 +423,7 @@ namespace gridfire {
std::vector<const reaction::Reaction *> AdaptiveEngineView::cullReactionsByFlow(
const std::vector<ReactionFlow> &allFlows,
const std::unordered_set<fourdst::atomic::Species> &reachableSpecies,
const std::vector<double> &Y_full,
const fourdst::composition::Composition &comp,
const double maxFlow
) const {
LOG_TRACE_L1(m_logger, "Culling reactions based on flow rates...");
@@ -464,9 +440,7 @@ namespace gridfire {
bool zero_flow_due_to_reachable_reactants = false;
if (flowRate < 1e-99 && flowRate > 0.0) {
for (const auto& reactant: reactionPtr->reactants()) {
const auto it = std::ranges::find(m_baseEngine.getNetworkSpecies(), reactant);
const size_t index = std::distance(m_baseEngine.getNetworkSpecies().begin(), it);
if (Y_full[index] < 1e-99 && reachableSpecies.contains(reactant)) {
if (comp.getMolarAbundance(reactant) < 1e-99 && reachableSpecies.contains(reactant)) {
LOG_TRACE_L1(m_logger, "Maintaining reaction '{}' with low flow ({:0.3E} [mol/s/g]) due to reachable reactant '{}'.", reactionPtr->id(), flowRate, reactant.name());
zero_flow_due_to_reachable_reactants = true;
break;
@@ -488,13 +462,13 @@ namespace gridfire {
}
AdaptiveEngineView::RescueSet AdaptiveEngineView::rescueEdgeSpeciesDestructionChannel(
const std::vector<double> &Y_full,
const fourdst::composition::Composition &comp,
const double T9,
const double rho,
const std::vector<Species> &activeSpecies,
const reaction::ReactionSet &activeReactions
) const {
const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
const auto result = m_baseEngine.getSpeciesTimescales(comp, T9, rho);
if (!result) {
LOG_ERROR(m_logger, "Failed to get species timescales due to stale engine state.");
throw exceptions::StaleEngineError("Failed to get species timescales");
@@ -565,8 +539,23 @@ namespace gridfire {
allOtherReactantsAreAvailable = false;
}
}
if (allOtherReactantsAreAvailable && speciesToCheckIsConsumed) {
double rate = reaction->calculate_rate(T9, rho, Y_full);
if (allOtherReactantsAreAvailable) {
std::vector<double> Y = comp.getMolarAbundanceVector();
const double Ye = comp.getElectronAbundance();
// TODO: This is a dummy placeholder which must be replaced with an EOS call
const double mue = 5.0e-3 * std::pow(rho * Ye, 1.0 / 3.0) + 0.5 * T9;
std::unordered_map<Species, double> speciesMassMap;
for (const auto &entry: comp | std::views::values) {
speciesMassMap[entry.isotope()] = entry.isotope().mass();
}
std::unordered_map<size_t, Species> speciesIndexMap;
for (const auto& entry: comp | std::views::values) {
size_t distance = std::distance(speciesMassMap.begin(), speciesMassMap.find(entry.isotope()));
speciesIndexMap.emplace(distance, entry.isotope());
}
double rate = reaction->calculate_rate(T9, rho, Ye, mue, Y, speciesIndexMap);
if (rate > maxSpeciesConsumptionRate) {
maxSpeciesConsumptionRate = rate;
reactionsToRescue[species] = reaction.get();

57
src/lib/engine/views/engine_defined.cpp
View File

@@ -28,58 +28,48 @@ namespace gridfire {
}
std::expected<StepDerivatives<double>, expectations::StaleEngineError> DefinedEngineView::calculateRHSAndEnergy(
const std::vector<double> &Y_defined,
const fourdst::composition::Composition& comp,
const double T9,
const double rho
) const {
validateNetworkState();
const auto Y_full = mapViewToFull(Y_defined);
const auto result = m_baseEngine.calculateRHSAndEnergy(Y_full, T9, rho);
const auto result = m_baseEngine.calculateRHSAndEnergy(comp, T9, rho);
if (!result) {
return std::unexpected{result.error()};
}
const auto [dydt, nuclearEnergyGenerationRate] = result.value();
StepDerivatives<double> definedResults;
definedResults.nuclearEnergyGenerationRate = nuclearEnergyGenerationRate;
definedResults.dydt = mapFullToView(dydt);
return definedResults;
return result.value();
}
EnergyDerivatives DefinedEngineView::calculateEpsDerivatives(
const std::vector<double> &Y_dynamic,
const fourdst::composition::Composition& comp,
const double T9,
const double rho
) const {
validateNetworkState();
const auto Y_full = mapViewToFull(Y_dynamic);
return m_baseEngine.calculateEpsDerivatives(Y_full, T9, rho);
return m_baseEngine.calculateEpsDerivatives(comp, T9, rho);
}
void DefinedEngineView::generateJacobianMatrix(
const std::vector<double> &Y_dynamic,
const fourdst::composition::Composition& comp,
const double T9,
const double rho
) const {
validateNetworkState();
const auto Y_full = mapViewToFull(Y_dynamic);
m_baseEngine.generateJacobianMatrix(Y_full, T9, rho);
m_baseEngine.generateJacobianMatrix(comp, T9, rho);
}
double DefinedEngineView::getJacobianMatrixEntry(
const int i_defined,
const int j_defined
const Species& rowSpecies,
const Species& colSpecies
) const {
validateNetworkState();
const size_t i_full = mapViewToFullSpeciesIndex(i_defined);
const size_t j_full = mapViewToFullSpeciesIndex(j_defined);
return m_baseEngine.getJacobianMatrixEntry(static_cast<int>(i_full), static_cast<int>(j_full));
return m_baseEngine.getJacobianMatrixEntry(rowSpecies, colSpecies);
}
void DefinedEngineView::generateStoichiometryMatrix() {
@@ -89,19 +79,17 @@ namespace gridfire {
}
int DefinedEngineView::getStoichiometryMatrixEntry(
const int speciesIndex_defined,
const int reactionIndex_defined
const Species& species,
const reaction::Reaction& reaction
) const {
validateNetworkState();
const size_t i_full = mapViewToFullSpeciesIndex(speciesIndex_defined);
const size_t j_full = mapViewToFullReactionIndex(reactionIndex_defined);
return m_baseEngine.getStoichiometryMatrixEntry(static_cast<int>(i_full), static_cast<int>(j_full));
return m_baseEngine.getStoichiometryMatrixEntry(species, reaction);
}
double DefinedEngineView::calculateMolarReactionFlow(
const reaction::Reaction &reaction,
const std::vector<double> &Y_defined,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
@@ -112,8 +100,7 @@ namespace gridfire {
m_logger -> flush_log();
throw std::runtime_error("Reaction not found in active reactions: " + std::string(reaction.id()));
}
const auto Y_full = mapViewToFull(Y_defined);
return m_baseEngine.calculateMolarReactionFlow(reaction, Y_full, T9, rho);
return m_baseEngine.calculateMolarReactionFlow(reaction, comp, T9, rho);
}
const reaction::ReactionSet & DefinedEngineView::getNetworkReactions() const {
@@ -131,14 +118,13 @@ namespace gridfire {
}
std::expected<std::unordered_map<Species, double>, expectations::StaleEngineError> DefinedEngineView::getSpeciesTimescales(
const std::vector<double> &Y_defined,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
validateNetworkState();
const auto Y_full = mapViewToFull(Y_defined);
const auto result = m_baseEngine.getSpeciesTimescales(Y_full, T9, rho);
const auto result = m_baseEngine.getSpeciesTimescales(comp, T9, rho);
if (!result) {
return std::unexpected{result.error()};
}
@@ -155,14 +141,13 @@ namespace gridfire {
std::expected<std::unordered_map<fourdst::atomic::Species, double>, expectations::StaleEngineError>
DefinedEngineView::getSpeciesDestructionTimescales(
const std::vector<double> &Y_defined,
const fourdst::composition::Composition &comp,
const double T9,
const double rho
) const {
validateNetworkState();
const auto Y_full = mapViewToFull(Y_defined);
const auto result = m_baseEngine.getSpeciesDestructionTimescales(Y_full, T9, rho);
const auto result = m_baseEngine.getSpeciesDestructionTimescales(comp, T9, rho);
if (!result) {
return std::unexpected{result.error()};
@@ -304,7 +289,7 @@ namespace gridfire {
}
size_t DefinedEngineView::mapViewToFullSpeciesIndex(size_t culledSpeciesIndex) const {
if (culledSpeciesIndex < 0 || culledSpeciesIndex >= m_speciesIndexMap.size()) {
if (culledSpeciesIndex >= m_speciesIndexMap.size()) {
LOG_ERROR(m_logger, "Defined index {} is out of bounds for species index map of size {}.", culledSpeciesIndex, m_speciesIndexMap.size());
m_logger->flush_log();
throw std::out_of_range("Defined index " + std::to_string(culledSpeciesIndex) + " is out of bounds for species index map of size " + std::to_string(m_speciesIndexMap.size()) + ".");
@@ -313,7 +298,7 @@ namespace gridfire {
}
size_t DefinedEngineView::mapViewToFullReactionIndex(size_t culledReactionIndex) const {
if (culledReactionIndex < 0 || culledReactionIndex >= m_reactionIndexMap.size()) {
if (culledReactionIndex >= m_reactionIndexMap.size()) {
LOG_ERROR(m_logger, "Defined index {} is out of bounds for reaction index map of size {}.", culledReactionIndex, m_reactionIndexMap.size());
m_logger->flush_log();
throw std::out_of_range("Defined index " + std::to_string(culledReactionIndex) + " is out of bounds for reaction index map of size " + std::to_string(m_reactionIndexMap.size()) + ".");

681
src/lib/engine/views/engine_multiscale.cpp
View File

File diff suppressed because it is too large Load Diff

2
src/lib/io/network_file.cpp
View File

@@ -39,7 +39,7 @@ namespace gridfire::io {
}
SimpleReactionListFileParser::SimpleReactionListFileParser() {}
SimpleReactionListFileParser::SimpleReactionListFileParser() = default;
ParsedNetworkData SimpleReactionListFileParser::parse(const std::string& filename) const {
LOG_TRACE_L1(m_logger, "Parsing simple reaction list file: {}", filename);

3
src/lib/network.cpp
View File

@@ -19,7 +19,6 @@
//
// *********************************************************************** */
#include "gridfire/network.h"
#include "gridfire/reaction/reaclib.h"
#include "gridfire/reaction/reaction.h"
#include <ranges>
@@ -74,7 +73,7 @@ namespace gridfire {
}
const auto ritr = std::find_if(str.rbegin(), std::string::const_reverse_iterator(startIt),
[](const unsigned char ch){ return !std::isspace(ch); });
return std::string(startIt, ritr.base());
return {startIt, ritr.base()};
}
}

1
src/lib/partition/partition_rauscher_thielemann.cpp
View File

@@ -7,7 +7,6 @@
#include <stdexcept>
#include <algorithm>
#include <vector>
#include <array>
#include <iostream>

56
src/lib/reaction/reaction.cpp
View File

@@ -45,15 +45,35 @@ namespace gridfire::reaction {
m_rateCoefficients(sets),
m_reverse(reverse) {}
double ReaclibReaction::calculate_rate(const double T9, const double rho, const std::vector<double>& Y) const {
double ReaclibReaction::calculate_rate(
const double T9,
const double rho,
double Ye,
double mue,
const std::vector<double> &Y,
const std::unordered_map<size_t, Species>& index_to_species_map
) const {
return calculate_rate<double>(T9);
}
CppAD::AD<double> ReaclibReaction::calculate_rate(const CppAD::AD<double> T9, const CppAD::AD<double> rho, const std::vector<CppAD::AD<double>>& Y) const {
CppAD::AD<double> ReaclibReaction::calculate_rate(
const CppAD::AD<double> T9,
const CppAD::AD<double> rho,
CppAD::AD<double> Ye,
CppAD::AD<double> mue,
const std::vector<CppAD::AD<double>>& Y,
const std::unordered_map<size_t, Species>& index_to_species_map
) const {
return calculate_rate<CppAD::AD<double>>(T9);
}
double ReaclibReaction::calculate_forward_rate_log_derivative(const double T9, const double rho, const std::vector<double>& Y) const {
double ReaclibReaction::calculate_forward_rate_log_derivative(
const double T9,
const double rho,
double Ye,
double mue,
const fourdst::composition::Composition& comp
) const {
constexpr double r_p13 = 1.0 / 3.0;
constexpr double r_p53 = 5.0 / 3.0;
constexpr double r_p23 = 2.0 / 3.0;
@@ -80,21 +100,11 @@ namespace gridfire::reaction {
bool ReaclibReaction::contains_reactant(const Species& species) const {
for (const auto& reactant : m_reactants) {
if (reactant == species) {
return true;
}
}
return false;
return std::ranges::any_of(m_reactants, [&](const Species& r) { return r == species; });
}
bool ReaclibReaction::contains_product(const Species& species) const {
for (const auto& product : m_products) {
if (product == species) {
return true;
}
}
return false;
return std::ranges::any_of(m_products, [&](const Species& p) { return p == species; });
}
std::unordered_set<Species> ReaclibReaction::all_species() const {
@@ -224,11 +234,20 @@ namespace gridfire::reaction {
m_rates.push_back(reaction.rateCoefficients());
}
double LogicalReaclibReaction::calculate_rate(const double T9, const double rho, const std::vector<double>& Y) const {
double LogicalReaclibReaction::calculate_rate(
const double T9,
const double rho,
double Ye,
double mue, const std::vector<double> &Y, const std::unordered_map<size_t, Species>& index_to_species_map
) const {
return calculate_rate<double>(T9);
}
double LogicalReaclibReaction::calculate_forward_rate_log_derivative(const double T9, const double rho, const std::vector<double>& Y) const {
double LogicalReaclibReaction::calculate_forward_rate_log_derivative(
const double T9,
const double rho,
double Ye, double mue, const fourdst::composition::Composition& comp
) const {
constexpr double r_p13 = 1.0 / 3.0;
constexpr double r_p53 = 5.0 / 3.0;
constexpr double r_p23 = 2.0 / 3.0;
@@ -286,7 +305,8 @@ namespace gridfire::reaction {
CppAD::AD<double> LogicalReaclibReaction::calculate_rate(
const CppAD::AD<double> T9,
const CppAD::AD<double> rho,
const std::vector<CppAD::AD<double>>& Y
CppAD::AD<double> Ye,
CppAD::AD<double> mue, const std::vector<CppAD::AD<double>>& Y, const std::unordered_map<size_t, Species>& index_to_species_map
) const {
return calculate_rate<CppAD::AD<double>>(T9);
}

450
src/lib/reaction/weak/weak.cpp
View File

@@ -7,20 +7,81 @@
#include <ranges>
#include <unordered_map>
#include <expected>
#include <vector>
#include <format>
#define GRIDFIRE_WEAK_REACTION_LIB_SENTINEL -60.0
#include "gridfire/reaction/weak/weak_interpolator.h"
#include "xxhash64.h"
namespace {
fourdst::atomic::Species resolve_weak_product(
const gridfire::rates::weak::WeakReactionType type,
const fourdst::atomic::Species& reactant
) {
using namespace fourdst::atomic;
using namespace gridfire::rates::weak;
std::optional<Species> product; // Use optional so that we can start in a valid "null" state
switch (type) {
case WeakReactionType::BETA_MINUS_DECAY:
product = az_to_species(reactant.a(), reactant.z() + 1);
return product.value();
case WeakReactionType::BETA_PLUS_DECAY:
product = az_to_species(reactant.a(), reactant.z() - 1);
return product.value();
case WeakReactionType::ELECTRON_CAPTURE:
product = az_to_species(reactant.a(), reactant.z() - 1);
return product.value();
case WeakReactionType::POSITRON_CAPTURE:
product = az_to_species(reactant.a(), reactant.z() + 1);
break;
}
if (!product.has_value()) {
throw std::runtime_error("Failed to resolve weak reaction product for reactant: " + std::string(reactant.name()));
}
return product.value();
}
std::string resolve_weak_id(
const gridfire::rates::weak::WeakReactionType type,
const fourdst::atomic::Species& reactant,
const fourdst::atomic::Species& product
) {
using namespace gridfire::rates::weak;
std::string id;
switch (type) {
case WeakReactionType::BETA_MINUS_DECAY:
id = std::format("{}(,ν|)e-,{}", reactant.name(), product.name());
break;
case WeakReactionType::BETA_PLUS_DECAY:
id = std::format("{}(,ν)e+,{}", reactant.name(), product.name());
break;
case WeakReactionType::ELECTRON_CAPTURE:
id = std::format("{}(e-,ν){}", reactant.name(), product.name());
break;
case WeakReactionType::POSITRON_CAPTURE:
id = std::format("{}(e+,ν|){}", reactant.name(), product.name());
break;
}
return id;
}
}
namespace gridfire::rates::weak {
WeakReactionMap::WeakReactionMap() {
using namespace fourdst::atomic;
// ReSharper disable once CppUseStructuredBinding
for (const auto& weak_reaction_record : UNIFIED_WEAK_DATA) {
Species species = az_to_species(weak_reaction_record.A, weak_reaction_record.Z);
if (weak_reaction_record.log_beta_minus > GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) {
m_weak_network[species].push_back(
WeakReaction{
WeakReactionEntry{
WeakReactionType::BETA_MINUS_DECAY,
weak_reaction_record.t9,
weak_reaction_record.log_rhoye,
@@ -32,7 +93,7 @@ namespace gridfire::rates::weak {
}
if (weak_reaction_record.log_beta_plus > GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) {
m_weak_network[species].push_back(
WeakReaction{
WeakReactionEntry{
WeakReactionType::BETA_PLUS_DECAY,
weak_reaction_record.t9,
weak_reaction_record.log_rhoye,
@@ -44,7 +105,7 @@ namespace gridfire::rates::weak {
}
if (weak_reaction_record.log_electron_capture > GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) {
m_weak_network[species].push_back(
WeakReaction{
WeakReactionEntry{
WeakReactionType::ELECTRON_CAPTURE,
weak_reaction_record.t9,
weak_reaction_record.log_rhoye,
@@ -56,7 +117,7 @@ namespace gridfire::rates::weak {
}
if (weak_reaction_record.log_positron_capture > GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) {
m_weak_network[species].push_back(
WeakReaction{
WeakReactionEntry{
WeakReactionType::POSITRON_CAPTURE,
weak_reaction_record.t9,
weak_reaction_record.log_rhoye,
@@ -69,26 +130,391 @@ namespace gridfire::rates::weak {
}
}
std::vector<WeakReaction> WeakReactionMap::get_all_reactions() const {
std::vector<WeakReaction> reactions;
std::vector<WeakReactionEntry> WeakReactionMap::get_all_reactions() const {
std::vector<WeakReactionEntry> reactions;
for (const auto &species_reactions: m_weak_network | std::views::values) {
reactions.insert(reactions.end(), species_reactions.begin(), species_reactions.end());
}
return reactions;
}
std::expected<std::vector<WeakReaction>, bool> WeakReactionMap::get_species_reactions(const fourdst::atomic::Species &species) const {
std::expected<std::vector<WeakReactionEntry>, WeakMapError> WeakReactionMap::get_species_reactions(
const fourdst::atomic::Species &species) const {
if (m_weak_network.contains(species)) {
return m_weak_network.at(species);
}
return std::unexpected(false);
return std::unexpected(WeakMapError::SPECIES_NOT_FOUND);
}
std::expected<std::vector<WeakReaction>, bool> WeakReactionMap::get_species_reactions(const std::string &species_name) const {
fourdst::atomic::Species species = fourdst::atomic::species.at(species_name);
std::expected<std::vector<WeakReactionEntry>, WeakMapError> WeakReactionMap::get_species_reactions(
const std::string &species_name) const {
const fourdst::atomic::Species species = fourdst::atomic::species.at(species_name);
if (m_weak_network.contains(species)) {
return m_weak_network.at(species);
}
return std::unexpected(false);
return std::unexpected(WeakMapError::SPECIES_NOT_FOUND);
}
WeakReaction::WeakReaction(
const fourdst::atomic::Species &species,
const WeakReactionType type,
const WeakRateInterpolator &interpolator
) :
m_reactant(species),
m_product(resolve_weak_product(type, species)),
m_reactant_a(species.a()),
m_reactant_z(species.z()),
m_product_a(m_product.a()),
m_product_z(m_product.z()),
m_id(resolve_weak_id(type, species, m_product)),
m_type(type),
m_interpolator(interpolator),
m_atomic(m_interpolator, m_reactant_a, m_reactant_z, m_type) {}
double WeakReaction::calculate_rate(
const double T9,
const double rho,
const double Ye,
const double mue,
const std::vector<double> &Y,
const std::unordered_map<size_t, fourdst::atomic::Species>& index_to_species_map
) const {
return calculate_rate<double>(T9, rho, Ye, mue, Y, index_to_species_map);
}
CppAD::AD<double> WeakReaction::calculate_rate(
CppAD::AD<double> T9,
CppAD::AD<double> rho,
CppAD::AD<double> Ye,
CppAD::AD<double> mue, const std::vector<CppAD::AD<double>> &Y, const std::unordered_map<size_t,fourdst::atomic::Species>& index_to_species_map
) const {
return static_cast<CppAD::AD<double>>(0.0);
}
bool WeakReaction::contains(const fourdst::atomic::Species &species) const {
return contains_reactant(species) || contains_product(species);
}
bool WeakReaction::contains_reactant(const fourdst::atomic::Species& species) const {
if (m_reactant == species) {
return true;
}
return false;
}
bool WeakReaction::contains_product(const fourdst::atomic::Species &species) const {
if (m_product == species) {
return true;
}
return false;
}
std::unordered_set<fourdst::atomic::Species> WeakReaction::all_species() const {
return {m_reactant, m_product};
}
std::unordered_set<fourdst::atomic::Species> WeakReaction::reactant_species() const {
return {m_reactant};
}
std::unordered_set<fourdst::atomic::Species> WeakReaction::product_species() const {
return {m_product};
}
int WeakReaction::stoichiometry(const fourdst::atomic::Species &species) const {
if (species == m_reactant) {
return -1;
}
if (species == m_product) {
return 1;
}
return 0;
}
std::unordered_map<fourdst::atomic::Species, int> WeakReaction::stoichiometry() const {
return {
{m_reactant, -1},
{m_product, 1}
};
}
uint64_t WeakReaction::hash(const uint64_t seed) const {
const std::string reaction_string = std::format(
"{}:{}({})",
m_reactant.name(),
m_product.name(),
static_cast<int>(m_type)
);
return XXHash64::hash(reaction_string.data(), reaction_string.size(), seed);
}
double WeakReaction::qValue() const {
// We ignore neutrino mass as it is negligible compared to other masses here.
double Q_MeV = 0.0;
const double parentMass_u = m_reactant.mass();
const double daughterMass_u = m_product.mass();
const double electronMass_MeV = m_constants.electronMassMeV;
const double nuclearMassDiff_MeV = (parentMass_u - daughterMass_u) * m_constants.u_to_MeV;
switch (m_type) {
case WeakReactionType::BETA_PLUS_DECAY:
Q_MeV = nuclearMassDiff_MeV - 2.0 * electronMass_MeV;
break;
case WeakReactionType::POSITRON_CAPTURE:
Q_MeV = nuclearMassDiff_MeV + 2.0 * electronMass_MeV;
break;
case WeakReactionType::BETA_MINUS_DECAY:
Q_MeV = nuclearMassDiff_MeV;
break;
case WeakReactionType::ELECTRON_CAPTURE:
Q_MeV = nuclearMassDiff_MeV;
break;
}
return Q_MeV;
}
double WeakReaction::calculate_energy_generation_rate(
const double T9,
const double rho,
const double Ye,
const double mue,
const std::vector<double> &Y,
const std::unordered_map<size_t, fourdst::atomic::Species> &index_to_species_map
) const {
std::expected<WeakRatePayload, InterpolationError> rates = m_interpolator.get_rates(
static_cast<uint16_t>(m_reactant_a),
static_cast<uint8_t>(m_reactant_z),
T9,
std::log10(rho * Ye),
mue
);
if (!rates.has_value()) {
const InterpolationErrorType type = rates.error().type;
const std::string msg = std::format(
"Failed to interpolate weak rate for (A={}, Z={}) at T9={}, log10(rho*Ye)={}, mu_e={} with error: {}",
m_reactant.name(), m_reactant_a, m_reactant_z, T9, std::log10(rho * Ye), mue, InterpolationErrorTypeMap.at(type)
);
throw std::runtime_error(msg);
}
double logNeutrinoLossRate = 0.0;
if (m_type == WeakReactionType::BETA_PLUS_DECAY || m_type == WeakReactionType::ELECTRON_CAPTURE) {
logNeutrinoLossRate = rates->log_antineutrino_loss_bd;
} else if (m_type == WeakReactionType::BETA_MINUS_DECAY || m_type == WeakReactionType::POSITRON_CAPTURE) {
logNeutrinoLossRate = rates->log_neutrino_loss_ec;
}
const double neutrinoLossRate = std::pow(10, logNeutrinoLossRate);
const double EDeposited_MeV = qValue() - neutrinoLossRate;
// We reimplement this logic here instead of calling calculate_rate() to avoid
// doing the interpolation twice (since the payload has already been interpolated).
const double logRate = get_log_rate_from_payload(rates.value());
double lambda = 0.0;
if (logRate > GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) {
lambda = std::pow(10, logRate);
}
return lambda * EDeposited_MeV; // returns in MeV / s
}
CppAD::AD<double> WeakReaction::calculate_energy_generation_rate(
const CppAD::AD<double> &T9,
const CppAD::AD<double> &rho,
const CppAD::AD<double> &Ye,
const CppAD::AD<double> &mue,
const std::vector<CppAD::AD<double>> &Y,
const std::unordered_map<size_t, fourdst::atomic::Species> &index_to_species_map
) const {
const CppAD::AD<double> log_rhoYe = CppAD::log10(rho * Ye);
std::vector<CppAD::AD<double>> ax = {T9, log_rhoYe, mue};
std::vector<CppAD::AD<double>> ay(1);
m_atomic(ax, ay); // TODO: Sort out why this isn't working and checkline 222 in weak.h where a similar line is
//TODO: think about how to get out neutrino loss in a autodiff safe way. This may mean I need to add an extra output to the atomic base
// so that I can get out both the rate and the neutrino loss rate. This will also mean that the sparsity pattern will need to
// be updated to account for the extra output.
CppAD::AD<double> rateConstant = ay[0];
}
std::unique_ptr<reaction::Reaction> WeakReaction::clone() const {
std::unique_ptr<reaction::Reaction> reaction_ptr = std::make_unique<WeakReaction>(
m_reactant,
m_type,
m_interpolator
);
return reaction_ptr;
}
double WeakReaction::get_log_rate_from_payload(const WeakRatePayload &payload) const {
double logRate = 0.0;
switch (m_type) {
case WeakReactionType::BETA_MINUS_DECAY:
logRate = payload.log_beta_minus;
break;
case WeakReactionType::BETA_PLUS_DECAY:
logRate = payload.log_beta_plus;
break;
case WeakReactionType::ELECTRON_CAPTURE:
logRate = payload.log_electron_capture;
break;
case WeakReactionType::POSITRON_CAPTURE:
logRate = payload.log_positron_capture;
break;
}
return logRate;
}
bool WeakReaction::AtomicWeakRate::forward (
const size_t p,
const size_t q,
const CppAD::vector<bool> &vx,
CppAD::vector<bool> &vy,
const CppAD::vector<double> &tx,
CppAD::vector<double> &ty
) {
// Doing this explicitly (only allowing p == 0) makes forward mode AD impossible for now.
if (p != 0) {
return false;
}
const double T9 = tx[0];
const double log10_rhoye = tx[1];
const double mu_e = tx[2];
const std::expected<WeakRatePayload, InterpolationError> result = m_interpolator.get_rates(
static_cast<uint16_t>(m_a),
static_cast<uint8_t>(m_z),
T9,
log10_rhoye,
mu_e
);
if (!result.has_value()) {
const InterpolationErrorType type = result.error().type;
std::string msg = std::format(
"Failed to interpolate weak rate for (A={}, Z={}) at T9={}, log10(rho*Ye)={}, mu_e={} with error: {}",
m_a, m_z, T9, log10_rhoye, mu_e, InterpolationErrorTypeMap.at(type)
);
}
switch (m_type) {
case WeakReactionType::BETA_MINUS_DECAY:
ty[0] = std::pow(10, result.value().log_beta_minus);
break;
case WeakReactionType::BETA_PLUS_DECAY:
ty[0] = std::pow(10, result.value().log_beta_plus);
break;
case WeakReactionType::ELECTRON_CAPTURE:
ty[0] = std::pow(10, result.value().log_electron_capture);
break;
case WeakReactionType::POSITRON_CAPTURE:
ty[0] = std::pow(10, result.value().log_positron_capture);
break;
}
if (vx.size() > 0) {
vy[0] = vx[0] || vx[1] || vx[2]; // Sets the output sparsity pattern
}
return true;
}
bool WeakReaction::AtomicWeakRate::reverse(
size_t q,
const CppAD::vector<double> &tx,
const CppAD::vector<double> &ty,
CppAD::vector<double> &px,
const CppAD::vector<double> &py
) {
const double T9 = tx[0];
const double log10_rhoye = tx[1];
const double mu_e = tx[2];
const std::expected<WeakRateDerivatives, InterpolationError> result = m_interpolator.get_rate_derivatives(
static_cast<uint16_t>(m_a),
static_cast<uint8_t>(m_z),
T9,
log10_rhoye,
mu_e
);
if (!result.has_value()) {
const InterpolationErrorType type = result.error().type;
const std::string msg = std::format(
"Failed to interpolate weak rate derivatives for (A={}, Z={}) at T9={}, log10(rho*Ye)={}, mu_e={} with error: {}",
m_a, m_z, T9, log10_rhoye, mu_e, InterpolationErrorTypeMap.at(type)
);
throw std::runtime_error(msg);
}
WeakRateDerivatives derivatives = result.value();
double dT9 = 0.0;
double dRho = 0.0;
double dMuE = 0.0;
switch (m_type) {
case WeakReactionType::BETA_MINUS_DECAY:
dT9 = py[0] * derivatives.d_log_beta_minus[0];
dRho = py[0] * derivatives.d_log_beta_minus[1];
dMuE = py[0] * derivatives.d_log_beta_minus[2];
break;
case WeakReactionType::BETA_PLUS_DECAY:
dT9 = py[0] * derivatives.d_log_beta_plus[0];
dRho = py[0] * derivatives.d_log_beta_plus[1];
dMuE = py[0] * derivatives.d_log_beta_plus[2];
break;
case WeakReactionType::ELECTRON_CAPTURE:
dT9 = py[0] * derivatives.d_log_electron_capture[0];
dRho = py[0] * derivatives.d_log_electron_capture[1];
dMuE = py[0] * derivatives.d_log_electron_capture[2];
break;
case WeakReactionType::POSITRON_CAPTURE:
dT9 = py[0] * derivatives.d_log_positron_capture[0];
dRho = py[0] * derivatives.d_log_positron_capture[1];
dMuE = py[0] * derivatives.d_log_positron_capture[2];
break;
}
px[0] = py[0] * dT9; // d(rate)/dT9
px[1] = py[0] * dRho; // d(rate)/dlogRhoYe
px[2] = py[0] * dMuE; // d(rate)/dMuE
return true;
}
bool WeakReaction::AtomicWeakRate::for_sparse_jac(
size_t q,
const CppAD::vector<std::set<size_t> > &r,
CppAD::vector<std::set<size_t> > &s
) {
s[0] = r[0];
s[0].insert(r[1].begin(), r[1].end());
s[0].insert(r[2].begin(), r[2].end());
return true;
}
bool WeakReaction::AtomicWeakRate::rev_sparse_jac(
size_t q,
const CppAD::vector<std::set<size_t> > &rt,
CppAD::vector<std::set<size_t> > &st
) {
// What this is saying is that each of the three input variables (T9, rho, Ye)
// all only affect the output variable (the rate) since there is only
// one output variable.
st[0] = rt[0];
st[1] = rt[0];
st[2] = rt[0];
return true;
}
}

308
src/lib/reaction/weak/weak_interpolator.cpp Normal file
View File

@@ -0,0 +1,308 @@
#include "gridfire/reaction/weak/weak_interpolator.h"
#include "gridfire/reaction/reaction.h"
#include "gridfire/reaction/weak/weak.h"
#include <algorithm>
#include <map>
#include <set>
#include <unordered_map>
#include <vector>
#include <optional>
#include <expected>
#include <ranges>
#include "fourdst/composition/species.h"
namespace gridfire::rates::weak {
WeakRateInterpolator::WeakRateInterpolator(const RowDataTable &raw_data) {
std::map<uint32_t, std::vector<const RateDataRow*>> grouped_rows;
for (const auto& row : raw_data) {
grouped_rows[pack_isotope_id(row.A, row.Z)].push_back(&row);
}
for (auto const& [isotope_id, rows] : grouped_rows) {
IsotopeGrid grid;
std::set<float> unique_t9, unique_rhoYe, unique_mue;
for (const auto* row : rows) {
unique_t9.emplace(row->t9);
unique_rhoYe.emplace(row->log_rhoye);
unique_mue.emplace(row->mu_e);
}
grid.t9_axis.reserve(unique_t9.size());
grid.rhoYe_axis.reserve(unique_rhoYe.size());
grid.mue_axis.reserve(unique_mue.size());
grid.t9_axis.insert(grid.t9_axis.begin(), unique_t9.begin(), unique_t9.end());
grid.rhoYe_axis.insert(grid.rhoYe_axis.begin(), unique_rhoYe.begin(), unique_rhoYe.end());
grid.mue_axis.insert(grid.mue_axis.begin(), unique_mue.begin(), unique_mue.end());
std::ranges::sort(grid.t9_axis);
std::ranges::sort(grid.rhoYe_axis);
std::ranges::sort(grid.mue_axis);
const size_t nt9 = grid.t9_axis.size();
const size_t nrhoYe = grid.rhoYe_axis.size();
const size_t nmue = grid.mue_axis.size();
grid.data.resize(nt9 * nrhoYe * nmue);
// Reverse map for quick index lookup
std::unordered_map<float, size_t> t9_map, rhoYe_map, mue_map;
for (size_t i = 0; i < nt9; i++) { t9_map[grid.t9_axis[i]] = i; }
for (size_t j = 0; j < nrhoYe; j++) { rhoYe_map[grid.rhoYe_axis[j]] = j; }
for (size_t k = 0; k < nmue; k++) { mue_map[grid.mue_axis[k]] = k; }
for (const auto* row: rows) {
size_t i_t9 = t9_map.at(row->t9);
size_t j_rhoYe = rhoYe_map.at(row->log_rhoye);
size_t k_mue = mue_map.at(row->mu_e);
size_t index = (i_t9 * nrhoYe + j_rhoYe) * nmue + k_mue;
grid.data[index] = WeakRatePayload{
row->log_beta_plus,
row->log_electron_capture,
row->log_neutrino_loss_ec,
row->log_beta_minus,
row->log_positron_capture,
row->log_antineutrino_loss_bd
};
}
m_rate_table[isotope_id] = std::move(grid);
}
}
std::vector<fourdst::atomic::Species> WeakRateInterpolator::available_isotopes() const {
std::vector<fourdst::atomic::Species> isotopes;
for (const auto &packed_id: m_rate_table | std::views::keys) {
const uint16_t A = static_cast<uint16_t>(packed_id >> 8);
const uint8_t Z = static_cast<uint8_t>(packed_id & 0xFF);
try {
fourdst::atomic::Species species = fourdst::atomic::az_to_species(A, Z);
isotopes.push_back(species);
} catch (const std::exception& e) {
throw std::runtime_error("Error converting A=" + std::to_string(A) + ", Z=" + std::to_string(Z) + " to Species: " + e.what());
}
}
return isotopes;
}
std::expected<WeakRatePayload, InterpolationError> WeakRateInterpolator::get_rates(
const uint16_t A,
const uint8_t Z,
const double t9,
const double log_rhoYe,
const double mu_e
) const {
const auto it = m_rate_table.find(pack_isotope_id(A, Z));
if (it == m_rate_table.end()) {
return std::unexpected(InterpolationError{InterpolationErrorType::UNKNOWN_SPECIES_ERROR});
}
const auto&[t9_axis, rhoYe_axis, mue_axis, data] = it->second;
// Now find the bracketing indices for t9, log_rhoYe, and mu_e
auto find_lower_index = [](const std::vector<double>& axis, const double value) -> std::optional<size_t> {
const auto upperBoundIterator = std::ranges::upper_bound(axis, value);
if (upperBoundIterator == axis.begin() || upperBoundIterator == axis.end()) {
return std::nullopt; // Out of bounds
}
return std::distance(axis.begin(), upperBoundIterator) - 1;
};
const auto i_t9_opt = find_lower_index(t9_axis, t9);
const auto j_rhoYe_opt = find_lower_index(rhoYe_axis, log_rhoYe);
const auto k_mue_opt = find_lower_index(mue_axis, mu_e);
if (!i_t9_opt || !j_rhoYe_opt || !k_mue_opt) {
std::unordered_map<TableAxes, BoundsErrorInfo> boundsInfo;
if (!i_t9_opt) {
boundsInfo[TableAxes::T9] = BoundsErrorInfo{
TableAxes::T9,
t9_axis.front(),
t9_axis.back(),
t9
};
}
if (!j_rhoYe_opt) {
boundsInfo[TableAxes::LOG_RHOYE] = BoundsErrorInfo{
TableAxes::LOG_RHOYE,
rhoYe_axis.front(),
rhoYe_axis.back(),
log_rhoYe
};
}
if (!k_mue_opt) {
boundsInfo[TableAxes::MUE] = BoundsErrorInfo{
TableAxes::MUE,
mue_axis.front(),
mue_axis.back(),
mu_e
};
}
return std::unexpected(
InterpolationError{
InterpolationErrorType::BOUNDS_ERROR,
boundsInfo
}
);
}
const size_t i = i_t9_opt.value();
const size_t j = j_rhoYe_opt.value();
const size_t k = k_mue_opt.value();
// Coordinates of the bounding cube
const double t1 = t9_axis[i];
const double t2 = t9_axis[i + 1];
const double r1 = rhoYe_axis[j];
const double r2 = rhoYe_axis[j + 1];
const double m1 = mue_axis[k];
const double m2 = mue_axis[k + 1];
const double td = (t9 - t1) / (t2 - t1);
const double rd = (log_rhoYe - r1) / (r2 - r1);
const double md = (mu_e - m1) / (m2 - m1);
auto lerp = [](const double v0, const double v1, const double t) {
return v0 * (1 - t) + v1 * t;
};
auto interpolationField = [&](auto field_accessor) {
const size_t nrhoYe = rhoYe_axis.size();
const size_t nmue = mue_axis.size();
auto get_val = [&](const size_t i_t, const size_t j_r, const size_t k_m) {
return field_accessor(data[(i_t * nrhoYe + j_r) * nmue + k_m]);
};
const double c000 = get_val(i, j, k);
const double c001 = get_val(i, j, k + 1);
const double c010 = get_val(i, j + 1, k);
const double c011 = get_val(i, j + 1, k + 1);
const double c100 = get_val(i + 1, j, k);
const double c101 = get_val(i + 1, j, k + 1);
const double c110 = get_val(i + 1, j + 1, k);
const double c111 = get_val(i + 1, j + 1, k + 1);
const double c00 = lerp(c000, c001, md);
const double c01 = lerp(c010, c011, md);
const double c10 = lerp(c100, c101, md);
const double c11 = lerp(c110, c111, md);
const double c0 = lerp(c00, c01, rd);
const double c1 = lerp(c10, c11, rd);
return lerp(c0, c1, td);
};
WeakRatePayload result;
result.log_beta_plus = interpolationField([](const WeakRatePayload& p) { return p.log_beta_plus; });
result.log_electron_capture = interpolationField([](const WeakRatePayload& p) { return p.log_electron_capture; });
result.log_neutrino_loss_ec = interpolationField([](const WeakRatePayload& p) { return p.log_neutrino_loss_ec; });
result.log_beta_minus = interpolationField([](const WeakRatePayload& p) { return p.log_beta_minus; });
result.log_positron_capture = interpolationField([](const WeakRatePayload& p) { return p.log_positron_capture; });
result.log_antineutrino_loss_bd = interpolationField([](const WeakRatePayload& p) { return p.log_antineutrino_loss_bd; });
return result;
}
std::expected<WeakRateDerivatives, InterpolationError> WeakRateInterpolator::get_rate_derivatives(
uint16_t A,
uint8_t Z,
double t9,
double log_rhoYe,
double mu_e
) const {
WeakRateDerivatives result;
constexpr double eps = 1e-6; // Small perturbation for finite difference
// Perturbations for finite difference
const double h_t9 = (t9 > 1e-9) ? t9 * eps : eps;
const auto payload_plus_t9 = get_rates(A, Z, t9 + h_t9, log_rhoYe, mu_e);
const auto payload_minus_t9 = get_rates(A, Z, t9 - h_t9, log_rhoYe, mu_e);
const double h_rhoYe = (std::abs(log_rhoYe) > 1e-9) ? std::abs(log_rhoYe) * eps : eps;
const auto payload_plus_rhoYe = get_rates(A, Z, t9, log_rhoYe + h_rhoYe, mu_e);
const auto payload_minus_rhoYe = get_rates(A, Z, t9, log_rhoYe - h_rhoYe, mu_e);
const double h_mue = (std::abs(mu_e) > 1e-9) ? std::abs(mu_e) * eps : eps;
const auto payload_plus_mue = get_rates(A, Z, t9, log_rhoYe, mu_e + h_mue);
const auto payload_minus_mue = get_rates(A, Z, t9, log_rhoYe, mu_e - h_mue);
if (!payload_plus_t9 || !payload_minus_t9 || !payload_plus_rhoYe || !payload_minus_rhoYe || !payload_plus_mue || !payload_minus_mue) {
const auto it = m_rate_table.find(pack_isotope_id(A, Z));
if (it == m_rate_table.end()) {
return std::unexpected(InterpolationError{InterpolationErrorType::UNKNOWN_SPECIES_ERROR});
}
const IsotopeGrid& grid = it->second;
InterpolationError error;
std::unordered_map<TableAxes, BoundsErrorInfo> boundsInfo;
if (!payload_minus_t9 || !payload_plus_t9) {
boundsInfo[TableAxes::T9] = BoundsErrorInfo{
TableAxes::T9,
grid.t9_axis.front(),
grid.t9_axis.back(),
t9
};
}
if (!payload_minus_rhoYe || !payload_plus_rhoYe) {
boundsInfo[TableAxes::LOG_RHOYE] = BoundsErrorInfo{
TableAxes::LOG_RHOYE,
grid.rhoYe_axis.front(),
grid.rhoYe_axis.back(),
log_rhoYe
};
}
if (!payload_minus_mue || !payload_plus_mue) {
boundsInfo[TableAxes::MUE] = BoundsErrorInfo{
TableAxes::MUE,
grid.mue_axis.front(),
grid.mue_axis.back(),
mu_e
};
}
error.type = InterpolationErrorType::BOUNDS_ERROR;
error.boundsErrorInfo = boundsInfo;
return std::unexpected(error);
}
// Derivatives wrt. T9
const double t9_denominator = 2 * h_t9;
result.d_log_beta_plus[0] = (payload_plus_t9->log_beta_plus - payload_minus_t9->log_beta_plus) / t9_denominator;
result.d_log_beta_minus[0] = (payload_plus_t9->log_beta_minus - payload_minus_t9->log_beta_minus) / t9_denominator;
result.d_log_electron_capture[0] = (payload_plus_t9->log_electron_capture - payload_minus_t9->log_electron_capture) / t9_denominator;
result.d_log_neutrino_loss_ec[0] = (payload_plus_t9->log_neutrino_loss_ec - payload_minus_t9->log_neutrino_loss_ec) / t9_denominator;
result.d_log_positron_capture[0] = (payload_plus_t9->log_positron_capture - payload_minus_t9->log_positron_capture) / t9_denominator;
result.d_log_antineutrino_loss_bd[0] = (payload_plus_t9->log_antineutrino_loss_bd - payload_minus_t9->log_antineutrino_loss_bd) / t9_denominator;
// Derivatives wrt. logRhoYe
const double rhoYe_denominator = 2 * h_rhoYe;
result.d_log_beta_plus[1] = (payload_plus_rhoYe->log_beta_plus - payload_minus_rhoYe->log_beta_plus) / rhoYe_denominator;
result.d_log_beta_minus[1] = (payload_plus_rhoYe->log_beta_minus - payload_minus_rhoYe->log_beta_minus) / rhoYe_denominator;
result.d_log_electron_capture[1] = (payload_plus_rhoYe->log_electron_capture - payload_minus_rhoYe->log_electron_capture) / rhoYe_denominator;
result.d_log_neutrino_loss_ec[1] = (payload_plus_rhoYe->log_neutrino_loss_ec - payload_minus_rhoYe->log_neutrino_loss_ec) / rhoYe_denominator;
result.d_log_positron_capture[1] = (payload_plus_rhoYe->log_positron_capture - payload_minus_rhoYe->log_positron_capture) / rhoYe_denominator;
result.d_log_antineutrino_loss_bd[1] = (payload_plus_rhoYe->log_antineutrino_loss_bd - payload_minus_rhoYe->log_antineutrino_loss_bd) / rhoYe_denominator;
// Derivatives wrt. MuE
const double mue_denominator = 2 * h_mue;
result.d_log_beta_plus[2] = (payload_plus_mue->log_beta_plus - payload_minus_mue->log_beta_plus) / mue_denominator;
result.d_log_beta_minus[2] = (payload_plus_mue->log_beta_minus - payload_minus_mue->log_beta_minus) / mue_denominator;
result.d_log_electron_capture[2] = (payload_plus_mue->log_electron_capture - payload_minus_mue->log_electron_capture) / mue_denominator;
result.d_log_neutrino_loss_ec[2] = (payload_plus_mue->log_neutrino_loss_ec - payload_minus_mue->log_neutrino_loss_ec) / mue_denominator;
result.d_log_positron_capture[2] = (payload_plus_mue->log_positron_capture - payload_minus_mue->log_positron_capture) / mue_denominator;
result.d_log_antineutrino_loss_bd[2] = (payload_plus_mue->log_antineutrino_loss_bd - payload_minus_mue->log_antineutrino_loss_bd) / mue_denominator;
return result;
}
uint32_t WeakRateInterpolator::pack_isotope_id(const uint16_t A, const uint8_t Z) {
return (static_cast<uint32_t>(A) << 8) | static_cast<uint32_t>(Z);
}
}

2
src/lib/screening/screening_bare.cpp
View File

@@ -12,7 +12,7 @@ namespace gridfire::screening {
std::vector<ADDouble> BareScreeningModel::calculateScreeningFactors(
const reaction::ReactionSet &reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<ADDouble> &Y,
const std::vector<screening::ADDouble> &Y,
const ADDouble T9,
const ADDouble rho
) const {

2
src/lib/screening/screening_weak.cpp
View File

@@ -12,7 +12,7 @@ namespace gridfire::screening {
std::vector<ADDouble> WeakScreeningModel::calculateScreeningFactors(
const reaction::ReactionSet &reactions,
const std::vector<fourdst::atomic::Species>& species,
const std::vector<ADDouble> &Y,
const std::vector<screening::ADDouble> &Y,
const ADDouble T9,
const ADDouble rho
) const {

356
src/lib/solver/solver.cpp
View File

@@ -1,356 +0,0 @@
#include "gridfire/solver/solver.h"
#include "gridfire/engine/engine_graph.h"
#include "gridfire/network.h"
#include "gridfire/exceptions/error_engine.h"
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/composition/composition.h"
#include "fourdst/config/config.h"
#include "fourdst/plugin/plugin.h"
#include "gridfire/interfaces/solver/solver_interfaces.h"
#include "unsupported/Eigen/NonLinearOptimization"
#include <boost/numeric/odeint.hpp>
#include <vector>
#include <string>
#include <stdexcept>
#include <iomanip>
#include "quill/LogMacros.h"
namespace gridfire::solver {
NetOut DirectNetworkSolver::evaluate(const NetIn &netIn) {
namespace ublas = boost::numeric::ublas;
namespace odeint = boost::numeric::odeint;
using fourdst::composition::Composition;
const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9)
const auto absTol = m_config.get<double>("gridfire:solver:DirectNetworkSolver:absTol", 1.0e-8);
const auto relTol = m_config.get<double>("gridfire:solver:DirectNetworkSolver:relTol", 1.0e-8);
Composition equilibratedComposition = m_engine.update(netIn);
size_t numSpecies = m_engine.getNetworkSpecies().size();
ublas::vector<double> Y(numSpecies + 1);
RHSManager manager(m_engine, T9, netIn.density, m_callback, m_engine.getNetworkSpecies());
JacobianFunctor jacobianFunctor(m_engine, T9, netIn.density);
auto populateY = [&](const Composition& comp) {
const size_t numSpeciesInternal = m_engine.getNetworkSpecies().size();
Y.resize(numSpeciesInternal + 1);
for (size_t i = 0; i < numSpeciesInternal; i++) {
const auto& species = m_engine.getNetworkSpecies()[i];
if (!comp.contains(species)) {
double lim = std::numeric_limits<double>::min();
LOG_DEBUG(m_logger, "Species '{}' not found in composition. Setting abundance to {:0.3E}.", species.name(), lim);
Y(i) = lim; // Species not in the composition, set to zero
} else {
Y(i) = comp.getMolarAbundance(species);
}
}
// TODO: a good starting point to make the temperature, density, and energy self consistent would be to turn this into an accumulator
Y(numSpeciesInternal) = 0.0; // Specific energy rate, initialized to zero
};
// This is a quick debug that can be turned on. For solar code input parameters (T~1.5e7K, ρ~1.5e3 g/cm^3) this should be near 8e-17
// std::cout << "D/H: " << equilibratedComposition.getMolarAbundance("H-2") / equilibratedComposition.getMolarAbundance("H-1") << std::endl;
populateY(equilibratedComposition);
const auto stepper = odeint::make_controlled<odeint::rosenbrock4<double>>(absTol, relTol);
double current_time = 0.0;
double current_initial_timestep = netIn.dt0;
double accumulated_energy = 0.0;
// size_t total_update_stages_triggered = 0;
while (current_time < netIn.tMax) {
try {
odeint::integrate_adaptive(
stepper,
std::make_pair(manager, jacobianFunctor),
Y,
current_time,
netIn.tMax,
current_initial_timestep,
[&](const auto& state, double t) {
current_time = t;
manager.observe(state, t);
}
);
current_time = netIn.tMax;
} catch (const exceptions::StaleEngineTrigger &e) {
LOG_INFO(m_logger, "Catching StaleEngineTrigger at t = {:0.3E} with T9 = {:0.3E}, rho = {:0.3E}. Triggering update stage (last stage took {} steps).", current_time, T9, netIn.density, e.totalSteps());
exceptions::StaleEngineTrigger::state staleState = e.getState();
accumulated_energy += e.energy(); // Add the specific energy rate to the accumulated energy
// total_update_stages_triggered++;
Composition temp_comp;
std::vector<double> mass_fractions;
size_t num_species_at_stop = e.numSpecies();
if (num_species_at_stop != m_engine.getNetworkSpecies().size()) {
throw std::runtime_error(
"StaleEngineError state has a different number of species than the engine. This should not happen."
);
}
mass_fractions.reserve(num_species_at_stop);
for (size_t i = 0; i < num_species_at_stop; ++i) {
const auto& species = m_engine.getNetworkSpecies()[i];
temp_comp.registerSpecies(species);
mass_fractions.push_back(e.getMolarAbundance(i) * species.mass()); // Convert from molar abundance to mass fraction
}
temp_comp.setMassFraction(m_engine.getNetworkSpecies(), mass_fractions);
temp_comp.finalize(true);
NetIn netInTemp = netIn;
netInTemp.temperature = e.temperature();
netInTemp.density = e.density();
netInTemp.composition = temp_comp;
Composition currentComposition = m_engine.update(netInTemp);
populateY(currentComposition);
Y(Y.size() - 1) = e.energy(); // Set the specific energy rate from the stale state
numSpecies = m_engine.getNetworkSpecies().size();
// current_initial_timestep = 0.001 * manager.m_last_step_time; // set the new timestep to the last successful timestep before repartitioning
}
}
accumulated_energy += Y(Y.size() - 1); // Add the specific energy rate to the accumulated energy
std::vector<double> finalMassFractions(numSpecies);
for (size_t i = 0; i < numSpecies; ++i) {
const double molarMass = m_engine.getNetworkSpecies()[i].mass();
finalMassFractions[i] = Y(i) * molarMass; // Convert from molar abundance to mass fraction
if (finalMassFractions[i] < MIN_ABUNDANCE_THRESHOLD) {
finalMassFractions[i] = 0.0;
}
}
std::vector<std::string> speciesNames;
speciesNames.reserve(numSpecies);
for (const auto& species : m_engine.getNetworkSpecies()) {
speciesNames.emplace_back(species.name());
}
Composition outputComposition(speciesNames);
outputComposition.setMassFraction(speciesNames, finalMassFractions);
outputComposition.finalize(true);
NetOut netOut;
netOut.composition = outputComposition;
netOut.energy = accumulated_energy; // Specific energy rate
netOut.num_steps = manager.m_num_steps;
auto [dEps_dT, dEps_dRho] = m_engine.calculateEpsDerivatives(
std::vector<double>(Y.begin(), Y.begin() + numSpecies), // TODO: This narrowing should probably be solved. Its possible unforeseen bugs will arise from this
T9,
netIn.density
);
netOut.dEps_dT = dEps_dT;
netOut.dEps_dRho = dEps_dRho;
return netOut;
}
void DirectNetworkSolver::set_callback(const std::any& callback) {
if (!callback.has_value()) {
m_callback = {};
return;
}
using FunctionPtrType = void (*)(const TimestepContext&);
if (callback.type() == typeid(TimestepCallback)) {
m_callback = std::any_cast<TimestepCallback>(callback);
}
else if (callback.type() == typeid(FunctionPtrType)) {
auto func_ptr = std::any_cast<FunctionPtrType>(callback);
m_callback = func_ptr;
}
else {
throw std::invalid_argument("Unsupported type passed to set_callback. "
"Provide a std::function or a matching function pointer.");
}
}
std::vector<std::tuple<std::string, std::string>> DirectNetworkSolver::describe_callback_context() const {
const TimestepContext context(
0.0, // time
boost::numeric::ublas::vector<double>(), // state
0.0, // dt
0.0, // cached_time
0.0, // last_observed_time
0.0, // last_step_time
0.0, // T9
0.0, // rho
std::nullopt, // cached_result
0, // num_steps
m_engine, // engine,
{}
);
return context.describe();
}
void DirectNetworkSolver::RHSManager::operator()(
const boost::numeric::ublas::vector<double> &Y,
boost::numeric::ublas::vector<double> &dYdt,
const double t
) const {
const size_t numSpecies = m_engine.getNetworkSpecies().size();
if (t != m_cached_time || !m_cached_result.has_value() || m_cached_result.value().dydt.size() != numSpecies + 1) {
compute_and_cache(Y, t);
}
const auto&[dydt, nuclearEnergyGenerationRate] = m_cached_result.value();
dYdt.resize(numSpecies + 1);
std::ranges::copy(dydt, dYdt.begin());
dYdt(numSpecies) = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate
}
void DirectNetworkSolver::RHSManager::observe(
const boost::numeric::ublas::vector<double> &state,
const double t
) const {
double dt = t - m_last_observed_time;
compute_and_cache(state, t);
LOG_INFO(
m_logger,
"(Step {}) Observed state at t = {:0.3E} (dt = {:0.3E})",
m_num_steps,
t,
dt
);
std::ostringstream oss;
oss << std::scientific << std::setprecision(3);
oss << "(Step: " << std::setw(10) << m_num_steps << ") t = " << t << " (dt = " << dt << ", eps_nuc: " << state(state.size() - 1) << " [erg])\n";
std::cout << oss.str();
fourdst::plugin::manager::PluginManager &pluginManager = fourdst::plugin::manager::PluginManager::getInstance();
if (pluginManager.has("gridfire/solver")) {
auto* plugin = pluginManager.get<SolverPluginInterface>("gridfire/solver");
plugin -> log_time(t, dt);
}
// Callback logic
if (m_callback) {
LOG_TRACE_L1(m_logger, "Calling user callback function at t = {:0.3E} with dt = {:0.3E}", t, dt);
const TimestepContext context(
t,
state,
dt,
m_cached_time,
m_last_observed_time,
m_last_step_time,
m_T9,
m_rho,
m_cached_result,
m_num_steps,
m_engine,
m_networkSpecies
);
m_callback(context);
LOG_TRACE_L1(m_logger, "User callback function completed at t = {:0.3E} with dt = {:0.3E}", t, dt);
}
m_last_observed_time = t;
m_last_step_time = dt;
}
void DirectNetworkSolver::RHSManager::compute_and_cache(
const boost::numeric::ublas::vector<double> &state,
double t
) const {
std::vector<double> y_vec(state.begin(), state.end() - 1);
std::ranges::replace_if(
y_vec,
[](const double yi){
return yi < 0.0;
},
0.0 // Avoid negative abundances
);
const auto result = m_engine.calculateRHSAndEnergy(y_vec, m_T9, m_rho);
if (!result) {
LOG_INFO(m_logger,
"Triggering update stage due to stale engine detected at t = {:0.3E} with T9 = {:0.3E}, rho = {:0.3E}, y_vec (size: {})",
t, m_T9, m_rho, y_vec.size());
throw exceptions::StaleEngineTrigger({m_T9, m_rho, y_vec, t, m_num_steps, state(state.size() - 1)});
}
m_cached_result = result.value();
m_cached_time = t;
m_num_steps++;
}
void DirectNetworkSolver::JacobianFunctor::operator()(
const boost::numeric::ublas::vector<double> &Y,
boost::numeric::ublas::matrix<double> &J,
double t,
boost::numeric::ublas::vector<double> &dfdt
) const {
size_t numSpecies = m_engine.getNetworkSpecies().size();
J.resize(numSpecies+1, numSpecies+1);
J.clear();
for (size_t i = 0; i < numSpecies; ++i) {
for (size_t j = 0; j < numSpecies; ++j) {
J(i, j) = m_engine.getJacobianMatrixEntry(i, j);
}
}
}
DirectNetworkSolver::TimestepContext::TimestepContext(
const double t,
const boost::numeric::ublas::vector<double> &state,
const double dt,
const double cached_time,
const double last_observed_time,
const double last_step_time,
const double t9,
const double rho,
const std::optional<StepDerivatives<double>> &cached_result,
const int num_steps,
const DynamicEngine &engine,
const std::vector<fourdst::atomic::Species> &networkSpecies
)
: t(t),
state(state),
dt(dt),
cached_time(cached_time),
last_observed_time(last_observed_time),
last_step_time(last_step_time),
T9(t9),
rho(rho),
cached_result(cached_result),
num_steps(num_steps),
engine(engine),
networkSpecies(networkSpecies) {}
std::vector<std::tuple<std::string, std::string>> DirectNetworkSolver::TimestepContext::describe() const {
return {
{"time", "double"},
{"state", "boost::numeric::ublas::vector<double>&"},
{"dt", "double"},
{"cached_time", "double"},
{"last_observed_time", "double"},
{"last_step_time", "double"},
{"T9", "double"},
{"rho", "double"},
{"cached_result", "std::optional<StepDerivatives<double>>&"},
{"num_steps", "int"},
{"engine", "DynamicEngine&"},
{"networkSpecies", "std::vector<fourdst::atomic::Species>&"}
};
}
}

49
src/lib/solver/strategies/CVODE_solver_strategy.cpp
View File

@@ -17,6 +17,7 @@
#include <algorithm>
#include "fourdst/composition/exceptions/exceptions_composition.h"
#include "gridfire/engine/engine_graph.h"
#include "gridfire/solver/strategies/triggers/engine_partitioning_trigger.h"
#include "gridfire/trigger/procedures/trigger_pprint.h"
@@ -71,7 +72,7 @@ namespace {
#elif SUNDIALS_HAVE_PTHREADS
N_Vector vec = N_VNew_Pthreads(size, sun_ctx);
#else
N_Vector vec = N_VNew_Serial(size, sun_ctx);
N_Vector vec = N_VNew_Serial(static_cast<long long>(size), sun_ctx);
#endif
check_cvode_flag(vec == nullptr ? -1 : 0, "N_VNew");
return vec;
@@ -87,7 +88,7 @@ namespace gridfire::solver {
const double last_step_time,
const double t9,
const double rho,
const int num_steps,
const size_t num_steps,
const DynamicEngine &engine,
const std::vector<fourdst::atomic::Species> &networkSpecies
) :
@@ -157,6 +158,7 @@ namespace gridfire::solver {
user_data.engine = &m_engine;
double current_time = 0;
// ReSharper disable once CppTooWideScope
[[maybe_unused]] double last_callback_time = 0;
m_num_steps = 0;
double accumulated_energy = 0.0;
@@ -169,7 +171,7 @@ namespace gridfire::solver {
check_cvode_flag(CVodeSetUserData(m_cvode_mem, &user_data), "CVodeSetUserData");
int flag = -1;
int flag{};
if (m_stdout_logging_enabled) {
flag = CVode(m_cvode_mem, netIn.tMax, m_Y, &current_time, CV_ONE_STEP);
} else {
@@ -304,10 +306,8 @@ namespace gridfire::solver {
netOut.energy = accumulated_energy;
check_cvode_flag(CVodeGetNumSteps(m_cvode_mem, reinterpret_cast<long int *>(&netOut.num_steps)), "CVodeGetNumSteps");
outputComposition.setCompositionMode(false); // set to number fraction mode
std::vector<double> Y = outputComposition.getNumberFractionVector(); // TODO need to ensure that the canonical vector representation is used throughout the code to make sure tracking does not get messed up
auto [dEps_dT, dEps_dRho] = m_engine.calculateEpsDerivatives(
std::vector<double>(Y.begin(), Y.begin() + numSpecies), // TODO: This narrowing should probably be solved. Its possible unforeseen bugs will arise from this
outputComposition,
T9,
netIn.density
);
@@ -374,9 +374,11 @@ namespace gridfire::solver {
for (size_t j = 0; j < numSpecies; ++j) {
for (size_t i = 0; i < numSpecies; ++i) {
const fourdst::atomic::Species& species_j = engine->getNetworkSpecies()[j];
const fourdst::atomic::Species& species_i = engine->getNetworkSpecies()[i];
// J(i,j) = d(f_i)/d(y_j)
// Column-major order format for SUNDenseMatrix: J_data[j*N + i]
J_data[j * N + i] = engine->getJacobianMatrixEntry(i, j);
J_data[j * N + i] = engine->getJacobianMatrixEntry(species_i, species_j);
}
}
@@ -398,11 +400,30 @@ namespace gridfire::solver {
const size_t numSpecies = m_engine.getNetworkSpecies().size();
sunrealtype* y_data = N_VGetArrayPointer(y);
// PERF: The trade off of ensured index consistency is some performance here. If this becomes a bottleneck we can revisit.
// The specific trade off is that we have decided to enforce that all interfaces accept composition objects rather
// than raw vectors of molar abundances. This then lets any method lookup the species by name rather than relying on
// the index in the vector being consistent. The trade off is that sometimes we need to construct a composition object
// which, at the moment, requires a somewhat expensive set of operations. Perhaps in the future we could enforce
// some consistent memory layout for the composition object to make this cheeper. That optimization would need to be
// done in the libcomposition library though...
std::vector<double> y_vec(y_data, y_data + numSpecies);
std::vector<std::string> symbols;
symbols.reserve(numSpecies);
for (const auto& species : m_engine.getNetworkSpecies()) {
symbols.emplace_back(species.name());
}
std::vector<double> X;
X.reserve(numSpecies);
for (size_t i = 0; i < numSpecies; ++i) {
const double molarMass = m_engine.getNetworkSpecies()[i].mass();
X.push_back(y_vec[i] * molarMass); // Convert from molar abundance to mass fraction
}
fourdst::composition::Composition composition(symbols, X);
std::ranges::replace_if(y_vec, [](const double val) { return val < 0.0; }, 0.0);
const auto result = m_engine.calculateRHSAndEnergy(y_vec, data->T9, data->rho);
const auto result = m_engine.calculateRHSAndEnergy(composition, data->T9, data->rho);
if (!result) {
throw exceptions::StaleEngineTrigger({data->T9, data->rho, y_vec, t, m_num_steps, y_data[numSpecies]});
}
@@ -411,7 +432,7 @@ namespace gridfire::solver {
const auto& [dydt, nuclearEnergyGenerationRate] = result.value();
for (size_t i = 0; i < numSpecies; ++i) {
ydot_data[i] = dydt[i];
ydot_data[i] = dydt.at(m_engine.getNetworkSpecies()[i]);
}
ydot_data[numSpecies] = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate
}
@@ -513,6 +534,7 @@ namespace gridfire::solver {
std::vector<double> Y_full(y_data, y_data + num_components - 1);
std::ranges::replace_if(
Y_full,
[](const double val) {
@@ -528,8 +550,9 @@ namespace gridfire::solver {
const double err_ratio = std::abs(y_err_data[i]) / weight;
err_ratios[i] = err_ratio;
speciesNames.push_back(std::string(user_data.networkSpecies->at(i).name()));
speciesNames.emplace_back(user_data.networkSpecies->at(i).name());
}
fourdst::composition::Composition composition(speciesNames, Y_full);
if (err_ratios.empty()) {
return;
@@ -565,9 +588,9 @@ namespace gridfire::solver {
columns.push_back(std::make_unique<utils::Column<double>>("Error Ratio", sorted_err_ratios));
std::cout << utils::format_table("Species Error Ratios", columns) << std::endl;
diagnostics::inspect_jacobian_stiffness(*user_data.engine, Y_full, user_data.T9, user_data.rho);
diagnostics::inspect_species_balance(*user_data.engine, "N-14", Y_full, user_data.T9, user_data.rho);
diagnostics::inspect_species_balance(*user_data.engine, "n-1", Y_full, user_data.T9, user_data.rho);
diagnostics::inspect_jacobian_stiffness(*user_data.engine, composition, user_data.T9, user_data.rho);
diagnostics::inspect_species_balance(*user_data.engine, "N-14", composition, user_data.T9, user_data.rho);
diagnostics::inspect_species_balance(*user_data.engine, "n-1", composition, user_data.T9, user_data.rho);
}
}

21
src/lib/solver/strategies/triggers/engine_partitioning_trigger.cpp
View File

@@ -94,13 +94,12 @@ namespace gridfire::trigger::solver::CVODE {
}
bool OffDiagonalTrigger::check(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const {
const size_t numSpecies = ctx.engine.getNetworkSpecies().size();
for (int row = 0; row < numSpecies; ++row) {
for (int col = 0; col < numSpecies; ++col) {
double DRowDCol = std::abs(ctx.engine.getJacobianMatrixEntry(row, col));
if (row != col && DRowDCol > m_threshold) {
for (const auto& rowSpecies : ctx.engine.getNetworkSpecies()) {
for (const auto& colSpecies : ctx.engine.getNetworkSpecies()) {
double DRowDCol = std::abs(ctx.engine.getJacobianMatrixEntry(rowSpecies, colSpecies));
if (rowSpecies != colSpecies && DRowDCol > m_threshold) {
m_hits++;
LOG_TRACE_L2(m_logger, "OffDiagonalTrigger triggered at t = {} due to entry ({}, {}) = {}", ctx.t, row, col, DRowDCol);
LOG_TRACE_L2(m_logger, "OffDiagonalTrigger triggered at t = {} due to entry ({}, {}) = {}", ctx.t, rowSpecies.name(), colSpecies.name(), DRowDCol);
return true;
}
}
@@ -173,7 +172,7 @@ namespace gridfire::trigger::solver::CVODE {
}
bool TimestepCollapseTrigger::check(const gridfire::solver::CVODESolverStrategy::TimestepContext &ctx) const {
if (m_timestep_window.size() < 1) {
if (m_timestep_window.empty()) {
m_misses++;
return false;
}
@@ -181,7 +180,7 @@ namespace gridfire::trigger::solver::CVODE {
for (const auto& dt : m_timestep_window) {
averageTimestep += dt;
}
averageTimestep /= m_timestep_window.size();
averageTimestep /= static_cast<double>(m_timestep_window.size());
if (m_relative && (std::abs(ctx.dt - averageTimestep) / averageTimestep) >= m_threshold) {
m_hits++;
LOG_TRACE_L2(m_logger, "TimestepCollapseTrigger triggered at t = {} due to relative growth: dt = {}, average dt = {}, threshold = {}", ctx.t, ctx.dt, averageTimestep, m_threshold);
@@ -250,6 +249,12 @@ namespace gridfire::trigger::solver::CVODE {
using ctx_t = gridfire::solver::CVODESolverStrategy::TimestepContext;
// Create the individual conditions that can trigger a repartitioning
// The current trigger logic is as follows
// 1. Trigger every 1000th time that the simulation time exceeds the simulationTimeInterval
// 2. OR if any off-diagonal Jacobian entry exceeds the offDiagonalThreshold
// 3. OR every 10th time that the timestep growth exceeds the timestepGrowthThreshold (relative or absolute)
// TODO: This logic likely needs to be revisited; however, for now it is easy enough to change and test and it works reasonably well
auto simulationTimeTrigger = std::make_unique<EveryNthTrigger<ctx_t>>(std::make_unique<SimulationTimeTrigger>(simulationTimeInterval), 1000);
auto offDiagTrigger = std::make_unique<OffDiagonalTrigger>(offDiagonalThreshold);
auto timestepGrowthTrigger = std::make_unique<EveryNthTrigger<ctx_t>>(std::make_unique<TimestepCollapseTrigger>(timestepGrowthThreshold, timestepGrowthRelative, timestepGrowthWindowSize), 10);

4
src/lib/utils/logging.cpp
View File

@@ -11,11 +11,11 @@
std::string gridfire::utils::formatNuclearTimescaleLogString(
const DynamicEngine& engine,
std::vector<double> const& Y,
const fourdst::composition::Composition& composition,
const double T9,
const double rho
) {
auto const& result = engine.getSpeciesTimescales(Y, T9, rho);
auto const& result = engine.getSpeciesTimescales(composition, T9, rho);
if (!result) {
std::ostringstream ss;
ss << "Failed to get species timescales: " << result.error();