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

fix(weakRates): major progress in resolving bugs

Browse Source 
bigs were introduced by the interface change from accepting raw molar abundance vectors to using the composition vector. This commit resolves many of these, including preformant ways to report that a species is not present in the composition and unified index lookups using composition object tooling.

BREAKING CHANGE:
This commit is contained in:
Emily Boudreaux
2025-10-10 09:12:40 -04:00
parent 13e2ea9ffa
commit 2f1077c02d
21 changed files with 17953 additions and 375 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -79,7 +79,19 @@ namespace gridfire {
// --- The public facing interface can always use the precomputed version since taping is done internally ---
return calculateAllDerivativesUsingPrecomputation(comp, bare_rates, bare_reverse_rates, T9, rho);
} else {
return calculateAllDerivatives<double>(comp.getMolarAbundanceVector(), T9, rho, Ye, mue);
return calculateAllDerivatives<double>(
comp.getMolarAbundanceVector(),
T9,
rho,
Ye,
mue,
[&comp](const fourdst::atomic::Species& species) -> std::optional<size_t> {
if (comp.contains(species)) {
return comp.getSpeciesIndex(species); // Return the index of the species in the composition
}
return std::nullopt; // Species not found in the composition
}
);
}
}
@@ -305,11 +317,18 @@ namespace gridfire {
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;
// Reverse reactions are only relevant for strong reactions (at least during the vast majority of stellar evolution)
// So here we just let these be dummy values since we know
// 1. The reaction should always be strong
// 2. The strong reaction rate is independent of Ye and mue
//
// In development builds the assert below will confirm this
constexpr double Ye = 0.0;
constexpr double mue = 0.0;
// It is a logic error to call this function on a weak reaction
assert(reaction.type() != gridfire::reaction::ReactionType::WEAK);
// 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");
@@ -317,7 +336,8 @@ 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, Ye, mue, comp.getMolarAbundanceVector(), m_indexToSpeciesMap);
// We also let Y be an empy vector since the strong reaction rate is independent of Y
const double forwardRate = reaction.calculate_rate(T9, rho, Ye, mue, {}, m_indexToSpeciesMap);
if (reaction.reactants().size() == 2 && reaction.products().size() == 2) {
reverseRate = calculateReverseRateTwoBody(reaction, T9, forwardRate, expFactor);
@@ -696,10 +716,20 @@ namespace gridfire {
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);
return calculateMolarReactionFlow<double>(
reaction,
comp.getMolarAbundanceVector(),
T9,
rho,
Ye,
0.0,
[&comp](const fourdst::atomic::Species& species) -> std::optional<size_t> {
if (comp.contains(species)) { // Species present in the composition
return comp.getSpeciesIndex(species);
}
return std::nullopt; // Species not present
}
);
}
void GraphEngine::generateJacobianMatrix(
@@ -908,10 +938,20 @@ namespace gridfire {
const double rho
) const {
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);
auto [dydt, _] = calculateAllDerivatives<double>(
comp.getMolarAbundanceVector(),
T9,
rho,
Ye,
0.0,
[&comp](const fourdst::atomic::Species& species) -> std::optional<size_t> {
if (comp.contains(species)) { // Species present in the composition
return comp.getSpeciesIndex(species);
}
return std::nullopt; // Species not present
}
);
std::unordered_map<fourdst::atomic::Species, double> speciesTimescales;
speciesTimescales.reserve(m_networkSpecies.size());
for (const auto& species : m_networkSpecies) {
@@ -930,17 +970,39 @@ namespace gridfire {
const double rho
) const {
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;
const std::vector<double>& Y = comp.getMolarAbundanceVector();
auto speciesLookup = [&comp](const fourdst::atomic::Species& species) -> std::optional<size_t> {
if (comp.contains(species)) { // Species present in the composition
return comp.getSpeciesIndex(species);
}
return std::nullopt; // Species not present
};
auto [dydt, _] = calculateAllDerivatives<double>(
Y,
T9,
rho,
Ye,
0.0,
speciesLookup
);
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, comp.getMolarAbundanceVector(), T9, rho, Ye, mue);
const auto flow = calculateMolarReactionFlow<double>(
*reaction,
Y,
T9,
rho,
Ye,
0.0,
speciesLookup
);
netDestructionFlow += flow;
}
}
@@ -979,7 +1041,7 @@ namespace gridfire {
m_logger->flush_log();
throw std::runtime_error("Cannot record AD tape: No species in the network.");
}
const size_t numADInputs = numSpecies + 2; // Note here that by not letting T9 and rho be independent variables, we are constraining the network to a constant temperature and density during each evaluation.
const size_t numADInputs = numSpecies + 2; // Y + T9 + rho
// --- CppAD Tape Recording ---
// 1. Declare independent variable (adY)
@@ -1004,13 +1066,22 @@ namespace gridfire {
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;
const CppAD::AD<double> adYe = 1e6;
const CppAD::AD<double> adMue = 10.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, adYe, adMue);
auto [dydt, nuclearEnergyGenerationRate] = calculateAllDerivatives<CppAD::AD<double>>(
adY,
adT9,
adRho,
adYe,
adMue,
[&](const fourdst::atomic::Species& querySpecies) -> size_t {
return m_speciesToIndexMap.at(querySpecies); // TODO: This is bad, needs to be fixed
}
);
// Extract the raw vector from the associative map
std::vector<CppAD::AD<double>> dydt_vec;
@@ -1049,10 +1120,19 @@ namespace gridfire {
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> adYe = 1e6;
const CppAD::AD<double> adMue = 1.0;
auto [dydt, nuclearEnergyGenerationRate] = calculateAllDerivatives<CppAD::AD<double>>(adY, adT9, adRho, adYe, adMue);
auto [dydt, nuclearEnergyGenerationRate] = calculateAllDerivatives<CppAD::AD<double>>(
adY,
adT9,
adRho,
adYe,
adMue,
[&](const fourdst::atomic::Species& querySpecies) -> size_t {
return m_speciesToIndexMap.at(querySpecies); // TODO: This is bad, needs to be fixed
}
);
std::vector<CppAD::AD<double>> adOutput(1);
adOutput[0] = nuclearEnergyGenerationRate;
@@ -1157,8 +1237,10 @@ namespace gridfire {
if ( p != 0) { return false; }
const double T9 = tx[0];
// TODO: Handle rho and Y
const double reverseRate = m_engine.calculateReverseRate(m_reaction, T9, 0, {});
// This is an interesting problem because the reverse rate should only ever be computed for strong reactions
// Which do not depend on rho or Y. However, the signature requires them...
// For now, we just pass dummy values for rho and Y
const double reverseRate = m_engine.calculateReverseRate(m_reaction, T9, 0.0, {});
// std::cout << m_reaction.peName() << " reverseRate: " << reverseRate << " at T9: " << T9 << "\n";
ty[0] = reverseRate; // Store the reverse rate in the output vector
@@ -1178,7 +1260,9 @@ namespace gridfire {
const double T9 = tx[0];
const double reverseRate = ty[0];
// TODO: Handle rho and Y
// This is an interesting problem because the reverse rate should only ever be computed for strong reactions
// Which do not depend on rho or Y. However, the signature requires them...
// For now, we just pass dummy values for rho and Y
const double derivative = m_engine.calculateReverseRateTwoBodyDerivative(m_reaction, T9, 0, {}, reverseRate);
px[0] = py[0] * derivative; // Return the derivative of the reverse rate with respect to T9

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

@@ -25,10 +25,11 @@ namespace gridfire {
ReactionSet build_nuclear_network(
const Composition& composition,
const rates::weak::WeakRateInterpolator& weak_interpolator,
BuildDepthType maxLayers,
bool reverse_reaclib) {
const Composition& composition,
const rates::weak::WeakRateInterpolator& weakInterpolator,
BuildDepthType maxLayers,
bool reverse
) {
int depth;
if (std::holds_alternative<NetworkBuildDepth>(maxLayers)) {
depth = static_cast<int>(std::get<NetworkBuildDepth>(maxLayers));
@@ -47,40 +48,52 @@ namespace gridfire {
// 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) {
if (reaction->is_reverse() == reverse) {
master_reaction_pool.add_reaction(reaction->clone());
}
}
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
)
for (const auto& parent_species: weakInterpolator.available_isotopes()) {
std::expected<Species, fourdst::atomic::SpeciesErrorType> upProduct = fourdst::atomic::az_to_species(
parent_species.a(),
parent_species.z() + 1
);
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
)
std::expected<Species, fourdst::atomic::SpeciesErrorType> downProduct = fourdst::atomic::az_to_species(
parent_species.a(),
parent_species.z() - 1
);
if (downProduct.has_value()) { // Only add the reaction if the Species map contains the product
master_reaction_pool.add_reaction(
std::make_unique<rates::weak::WeakReaction>(
parent_species,
rates::weak::WeakReactionType::BETA_PLUS_DECAY,
weakInterpolator
)
);
master_reaction_pool.add_reaction(
std::make_unique<rates::weak::WeakReaction>(
parent_species,
rates::weak::WeakReactionType::ELECTRON_CAPTURE,
weakInterpolator
)
);
}
if (upProduct.has_value()) { // Only add the reaction if the Species map contains the product
master_reaction_pool.add_reaction(
std::make_unique<rates::weak::WeakReaction>(
parent_species,
rates::weak::WeakReactionType::BETA_MINUS_DECAY,
weakInterpolator
)
);
master_reaction_pool.add_reaction(
std::make_unique<rates::weak::WeakReaction>(
parent_species,
rates::weak::WeakReactionType::POSITRON_CAPTURE,
weakInterpolator
)
);
}
}
// --- Step 2: Use non-owning raw pointers for the fast build algorithm ---
@@ -140,7 +153,13 @@ namespace gridfire {
break;
}
LOG_TRACE_L1(logger, "Layer {}: Collected {} new reactions. New products this layer: {}", collectedReactionPtrs.size() - collectedReactionPtrs.size(), newProductsThisLayer.size());
LOG_TRACE_L1(
logger,
"Layer {}: Collected {} new reactions. New products this layer: {}",
layer,
collectedReactionPtrs.size() - collectedReactionPtrs.size(),
newProductsThisLayer.size()
);
availableSpecies.insert(newProductsThisLayer.begin(), newProductsThisLayer.end());
remainingReactions = std::move(reactionsForNextPass);
}

2
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 fourdst::composition::Composition &comp,
const Composition &comp,
const double T9,
const double rho
) {

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

@@ -17,7 +17,13 @@
namespace {
using namespace fourdst::atomic;
std::vector<double> packCompositionToVector(const fourdst::composition::Composition& composition, const gridfire::GraphEngine& engine) {
//TODO: Replace all calls to this function with composition.getMolarAbundanceVector() so that
// we don't have to keep this function around. (Cant do this right now because there is not a
// guarantee that this function will return the same ordering as the canonical vector representation)
std::vector<double> packCompositionToVector(
const fourdst::composition::Composition& composition,
const gridfire::GraphEngine& engine
) {
std::vector<double> Y(engine.getNetworkSpecies().size(), 0.0);
const auto& allSpecies = engine.getNetworkSpecies();
for (size_t i = 0; i < allSpecies.size(); ++i) {
@@ -780,6 +786,7 @@ namespace gridfire {
LOG_TRACE_L1(m_logger, "Partitioning by timescale...");
const auto result= m_baseEngine.getSpeciesDestructionTimescales(comp, T9, rho);
const auto netTimescale = m_baseEngine.getSpeciesTimescales(comp, T9, rho);
if (!result) {
LOG_ERROR(m_logger, "Failed to get species destruction timescales due to stale engine state");
m_logger->flush_log();
@@ -792,6 +799,14 @@ namespace gridfire {
}
const std::unordered_map<Species, double>& all_timescales = result.value();
const std::unordered_map<Species, double>& net_timescales = netTimescale.value();
for (const auto& [species, destructionTimescale] : all_timescales) {
std::cout << "Species: " << species.name()
<< ", Destruction Timescale: " << (std::isfinite(destructionTimescale) ? std::to_string(destructionTimescale) : "inf")
<< " s, Net Timescale: " << (std::isfinite(net_timescales.at(species)) ? std::to_string(net_timescales.at(species)) : "inf")
<< " s\n";
}
const auto& all_species = m_baseEngine.getNetworkSpecies();
std::vector<std::pair<double, Species>> sorted_timescales;
@@ -1075,6 +1090,11 @@ namespace gridfire {
normalized_composition.setMassFraction(species, 0.0);
}
}
bool normCompFinalizedOkay = normalized_composition.finalize(true);
if (!normCompFinalizedOkay) {
LOG_ERROR(m_logger, "Failed to finalize composition before QSE solve.");
throw std::runtime_error("Failed to finalize composition before QSE solve.");
}
Eigen::VectorXd Y_scale(algebraic_species.size());
Eigen::VectorXd v_initial(algebraic_species.size());
@@ -1330,10 +1350,18 @@ namespace gridfire {
Eigen::VectorXd y_qse = m_Y_scale.array() * v_qse.array().sinh(); // Convert to physical abundances using asinh scaling
for (const auto& species: m_qse_solve_species) {
if (!comp_trial.contains(species)) {
if (!comp_trial.hasSymbol(std::string(species.name()))) {
comp_trial.registerSpecies(species);
}
comp_trial.setMassFraction(species, y_qse[static_cast<long>(m_qse_solve_species_index_map.at(species))]);
const double molarAbundance = y_qse[static_cast<long>(m_qse_solve_species_index_map.at(species))];
const double massFraction = molarAbundance * species.mass();
comp_trial.setMassFraction(species, massFraction);
}
const bool didFinalize = comp_trial.finalize(false);
if (!didFinalize) {
std::string msg = std::format("Failed to finalize composition (comp_trial) in {} at line {}", __FILE__, __LINE__);
throw std::runtime_error(msg);
}
const auto result = m_view->getBaseEngine().calculateRHSAndEnergy(comp_trial, m_T9, m_rho);
@@ -1357,10 +1385,18 @@ namespace gridfire {
Eigen::VectorXd y_qse = m_Y_scale.array() * v_qse.array().sinh(); // Convert to physical abundances using asinh scaling
for (const auto& species: m_qse_solve_species) {
if (!comp_trial.contains(species)) {
if (!comp_trial.hasSymbol(std::string(species.name()))) {
comp_trial.registerSpecies(species);
}
comp_trial.setMassFraction(species, y_qse[static_cast<long>(m_qse_solve_species_index_map.at(species))]);
const double molarAbundance = y_qse[static_cast<long>(m_qse_solve_species_index_map.at(species))];
const double massFraction = molarAbundance * species.mass();
comp_trial.setMassFraction(species, massFraction);
}
const bool didFinalize = comp_trial.finalize(false);
if (!didFinalize) {
std::string msg = std::format("Failed to finalize composition (comp_trial) in {} at line {}", __FILE__, __LINE__);
throw std::runtime_error(msg);
}
m_view->getBaseEngine().generateJacobianMatrix(comp_trial, m_T9, m_rho);
@@ -1368,14 +1404,16 @@ namespace gridfire {
const long N = static_cast<long>(m_qse_solve_species.size());
J_qse.resize(N, N);
long rowID = 0;
long colID = 0;
for (const auto& rowSpecies : m_qse_solve_species) {
long colID = 0;
for (const auto& colSpecies: m_qse_solve_species) {
J_qse(rowID++, colID++) = m_view->getBaseEngine().getJacobianMatrixEntry(
J_qse(rowID, colID) = m_view->getBaseEngine().getJacobianMatrixEntry(
rowSpecies,
colSpecies
);
colID += 1;
}
rowID += 1;
}
// Chain rule for asinh scaling:

4
src/lib/partition/partition_ground.cpp
View File

@@ -19,7 +19,7 @@ namespace gridfire::partition {
const int a,
const double T9
) const {
LOG_TRACE_L2(m_logger, "Evaluating ground state partition function for Z={} A={} T9={}", z, a, T9);
LOG_TRACE_L3(m_logger, "Evaluating ground state partition function for Z={} A={} T9={}", z, a, T9);
const int key = make_key(z, a);
const double spin = m_ground_state_spin.at(key);
return (2.0 * spin) + 1.0;
@@ -30,7 +30,7 @@ namespace gridfire::partition {
const int a,
const double T9
) const {
LOG_TRACE_L2(m_logger, "Evaluating derivative of ground state partition function for Z={} A={} T9={}", z, a, T9);
LOG_TRACE_L3(m_logger, "Evaluating derivative of ground state partition function for Z={} A={} T9={}", z, a, T9);
return 0.0;
}

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

@@ -44,21 +44,21 @@ namespace gridfire::partition {
const int a,
const double T9
) const {
LOG_TRACE_L2(m_logger, "Evaluating Rauscher-Thielemann partition function for Z={} A={} T9={}", z, a, T9);
LOG_TRACE_L3(m_logger, "Evaluating Rauscher-Thielemann partition function for Z={} A={} T9={}", z, a, T9);
const auto [bound, data, upperIndex, lowerIndex] = find(z, a, T9);
switch (bound) {
case FRONT: {
LOG_TRACE_L2(m_logger, "Using FRONT bound for Z={} A={} T9={}", z, a, T9);
LOG_TRACE_L3(m_logger, "Using FRONT bound for Z={} A={} T9={}", z, a, T9);
return data.normalized_g_values.front() * (2.0 * data.ground_state_spin + 1.0);
}
case BACK: {
LOG_TRACE_L2(m_logger, "Using BACK bound for Z={} A={} T9={}", z, a, T9);
LOG_TRACE_L3(m_logger, "Using BACK bound for Z={} A={} T9={}", z, a, T9);
return data.normalized_g_values.back() * (2.0 * data.ground_state_spin + 1.0);
}
case MIDDLE: {
LOG_TRACE_L2(m_logger, "Using MIDDLE bound for Z={} A={} T9={}", z, a, T9);
LOG_TRACE_L3(m_logger, "Using MIDDLE bound for Z={} A={} T9={}", z, a, T9);
}
}
@@ -79,10 +79,10 @@ namespace gridfire::partition {
const int a,
const double T9
) const {
LOG_TRACE_L2(m_logger, "Evaluating derivative of Rauscher-Thielemann partition function for Z={} A={} T9={}", z, a, T9);
LOG_TRACE_L3(m_logger, "Evaluating derivative of Rauscher-Thielemann partition function for Z={} A={} T9={}", z, a, T9);
const auto [bound, data, upperIndex, lowerIndex] = find(z, a, T9);
if (bound == FRONT || bound == BACK) {
LOG_TRACE_L2(m_logger, "Derivative is zero for Z={} A={} T9={} (bound: {})", z, a, T9, bound == FRONT ? "FRONT" : "BACK");
LOG_TRACE_L3(m_logger, "Derivative is zero for Z={} A={} T9={} (bound: {})", z, a, T9, bound == FRONT ? "FRONT" : "BACK");
return 0.0; // Derivative is zero at the boundaries
}
const auto [T9_high, G_norm_high, T9_low, G_norm_low] = get_interpolation_points(

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

@@ -16,6 +16,11 @@
namespace {
std::unordered_map<fourdst::atomic::SpeciesErrorType, std::string> SpeciesErrorTypeMap = {
{fourdst::atomic::SpeciesErrorType::ELEMENT_SYMBOL_NOT_FOUND, "Element symbol not found (Z out of range)"},
{fourdst::atomic::SpeciesErrorType::SPECIES_SYMBOL_NOT_FOUND, "Species symbol not found ((A,Z) out of range)"}
};
fourdst::atomic::Species resolve_weak_product(
const gridfire::rates::weak::WeakReactionType type,
const fourdst::atomic::Species& reactant
@@ -23,25 +28,36 @@ namespace {
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
int zMod = 0;
switch (type) {
case WeakReactionType::BETA_MINUS_DECAY:
product = az_to_species(reactant.a(), reactant.z() + 1);
return product.value();
zMod = 1;
break;
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();
zMod = -1;
break;
case WeakReactionType::POSITRON_CAPTURE:
product = az_to_species(reactant.a(), reactant.z() + 1);
zMod = 1;
break;
case WeakReactionType::ELECTRON_CAPTURE:
zMod = -1;
break;
}
if (!product.has_value()) {
throw std::runtime_error("Failed to resolve weak reaction product for reactant: " + std::string(reactant.name()));
std::expected<Species, SpeciesErrorType> product = az_to_species(reactant.a(), reactant.z() + zMod);
if (product.has_value()) {
return product.value();
}
return product.value();
const std::string msg = std::format(
"Failed to resolve weak reaction product (A: {}, Z: {}) for reactant {} (looked up A: {}, Z: {}): {}",
reactant.a(),
reactant.z(),
reactant.name(),
reactant.a(),
reactant.z() + zMod,
SpeciesErrorTypeMap.at(product.error())
);
throw std::runtime_error(msg);
}
std::string resolve_weak_id(
@@ -77,7 +93,16 @@ namespace gridfire::rates::weak {
// 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);
std::expected<Species, SpeciesErrorType> species_result = az_to_species(weak_reaction_record.A, weak_reaction_record.Z);
if (!species_result.has_value()) {
const SpeciesErrorType type = species_result.error();
const std::string msg = std::format(
"Failed to load weak reaction data for (A={}, Z={}) with error: {}",
weak_reaction_record.A, weak_reaction_record.Z, SpeciesErrorTypeMap.at(type)
);
throw std::runtime_error(msg);
}
const Species& species = species_result.value();
if (weak_reaction_record.log_beta_minus > GRIDFIRE_WEAK_REACTION_LIB_SENTINEL) {
m_weak_network[species].push_back(
@@ -148,7 +173,7 @@ namespace gridfire::rates::weak {
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);
const fourdst::atomic::Species& species = fourdst::atomic::species.at(species_name);
if (m_weak_network.contains(species)) {
return m_weak_network.at(species);
}
@@ -284,9 +309,7 @@ namespace gridfire::rates::weak {
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::BETA_MINUS_DECAY: // Same as electron capture so we can simply fall through
case WeakReactionType::ELECTRON_CAPTURE:
Q_MeV = nuclearMassDiff_MeV;
break;
@@ -306,8 +329,7 @@ namespace gridfire::rates::weak {
static_cast<uint16_t>(m_reactant_a),
static_cast<uint8_t>(m_reactant_z),
T9,
std::log10(rho * Ye),
mue
std::log10(rho * Ye)
);
if (!rates.has_value()) {
@@ -374,8 +396,7 @@ namespace gridfire::rates::weak {
static_cast<uint16_t>(m_reactant_a),
static_cast<uint8_t>(m_reactant_z),
T9,
log_rhoYe,
mue
log_rhoYe
);
if (!rates.has_value()) {
const InterpolationErrorType type = rates.error().type;
@@ -386,9 +407,22 @@ namespace gridfire::rates::weak {
throw std::runtime_error(msg);
}
// TODO: Finish implementing this (just need a switch statement)
return 0.0;
double logRate = 0.0;
switch (m_type) {
case WeakReactionType::BETA_MINUS_DECAY:
logRate = rates->d_log_beta_minus[0];
break;
case WeakReactionType::BETA_PLUS_DECAY:
logRate = rates->d_log_beta_plus[0];
break;
case WeakReactionType::ELECTRON_CAPTURE:
logRate = rates->d_log_electron_capture[0];
break;
case WeakReactionType::POSITRON_CAPTURE:
logRate = rates->d_log_positron_capture[0];
break;
}
return logRate;
}
reaction::ReactionType WeakReaction::type() const {
@@ -437,9 +471,7 @@ namespace gridfire::rates::weak {
case WeakReactionType::BETA_MINUS_DECAY:
logNeutrinoLoss = payload.log_antineutrino_loss_bd;
break;
case WeakReactionType::BETA_PLUS_DECAY:
logNeutrinoLoss = payload.log_neutrino_loss_ec;
break;
case WeakReactionType::BETA_PLUS_DECAY: // Same as electron capture so we can simply fall through
case WeakReactionType::ELECTRON_CAPTURE:
logNeutrinoLoss = payload.log_neutrino_loss_ec;
break;
@@ -450,6 +482,7 @@ namespace gridfire::rates::weak {
return logNeutrinoLoss;
}
// Note that the input vector tx is of size 3: [T9, log10(rho*Ye), mu_e]
bool WeakReaction::AtomicWeakRate::forward (
const size_t p,
const size_t q,
@@ -464,20 +497,18 @@ namespace gridfire::rates::weak {
}
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
log10_rhoye
);
if (!result.has_value()) {
const InterpolationErrorType type = result.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_a, m_z, T9, log10_rhoye, mu_e, InterpolationErrorTypeMap.at(type)
"Failed to interpolate weak rate for (A={}, Z={}) at T9={}, log10(ρ Ye)={}, with error: {}",
m_a, m_z, T9, log10_rhoye, InterpolationErrorTypeMap.at(type)
);
throw std::runtime_error(msg);
}
@@ -501,7 +532,7 @@ namespace gridfire::rates::weak {
}
if (vx.size() > 0) { // Set up the sparsity pattern. This is saying that all input variables affect the output variable.
const bool any_input_varies = vx[0] || vx[1] || vx[2];
const bool any_input_varies = vx[0] || vx[1];
vy[0] = any_input_varies;
vy[1] = any_input_varies;
}
@@ -517,7 +548,6 @@ namespace gridfire::rates::weak {
) {
const double T9 = tx[0];
const double log10_rhoye = tx[1];
const double mu_e = tx[2];
const double forwardPassRate = ty[0]; // This is the rate from the forward pass.
const double forwardPassNeutrinoLossRate = ty[1]; // This is the neutrino loss rate from the forward pass.
@@ -526,55 +556,47 @@ namespace gridfire::rates::weak {
static_cast<uint16_t>(m_a),
static_cast<uint8_t>(m_z),
T9,
log10_rhoye,
mu_e
log10_rhoye
);
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)
"Failed to interpolate weak rate derivatives for (A={}, Z={}) at T9={}, log10(ρ Ye)={}, with error: {}",
m_a, m_z, T9, log10_rhoye, InterpolationErrorTypeMap.at(type)
);
throw std::runtime_error(msg);
}
// ReSharper disable once CppUseStructuredBinding
const WeakRateDerivatives derivatives = result.value();
std::array<double, 3> dLogRate; // d(rate)/dT9, d(rate)/dlogRhoYe, d(rate)/dMuE
std::array<double, 3> dLogNuLoss; // d(nu loss)/dT9, d(nu loss)/dlogRhoYe, d(nu loss)/dMuE
std::array<double, 2> dLogRate{}; // d(rate)/dT9, d(rate)/dlogRhoYe
std::array<double, 2> dLogNuLoss{}; // d(nu loss)/dT9, d(nu loss)/dlogRhoYe
switch (m_type) {
case WeakReactionType::BETA_MINUS_DECAY:
dLogRate[0] = derivatives.d_log_beta_minus[0];
dLogRate[1] = derivatives.d_log_beta_minus[1];
dLogRate[2] = derivatives.d_log_beta_minus[2];
dLogNuLoss[0] = derivatives.d_log_antineutrino_loss_bd[0];
dLogNuLoss[1] = derivatives.d_log_antineutrino_loss_bd[1];
dLogNuLoss[2] = derivatives.d_log_antineutrino_loss_bd[2];
break;
case WeakReactionType::BETA_PLUS_DECAY:
dLogRate[0] = derivatives.d_log_beta_plus[0];
dLogRate[1] = derivatives.d_log_beta_plus[1];
dLogRate[2] = derivatives.d_log_beta_plus[2];
dLogNuLoss[0] = derivatives.d_log_neutrino_loss_ec[0];
dLogNuLoss[1] = derivatives.d_log_neutrino_loss_ec[1];
dLogNuLoss[2] = derivatives.d_log_neutrino_loss_ec[2];
break;
case WeakReactionType::ELECTRON_CAPTURE:
dLogRate[0] = derivatives.d_log_electron_capture[0];
dLogRate[1] = derivatives.d_log_electron_capture[1];
dLogRate[2] = derivatives.d_log_electron_capture[2];
dLogNuLoss[0] = derivatives.d_log_neutrino_loss_ec[0];
dLogNuLoss[1] = derivatives.d_log_neutrino_loss_ec[1];
dLogNuLoss[2] = derivatives.d_log_neutrino_loss_ec[2];
break;
case WeakReactionType::POSITRON_CAPTURE:
dLogRate[0] = derivatives.d_log_positron_capture[0];
dLogRate[1] = derivatives.d_log_positron_capture[1];
dLogRate[2] = derivatives.d_log_positron_capture[2];
dLogNuLoss[0] = derivatives.d_log_antineutrino_loss_bd[0];
dLogNuLoss[1] = derivatives.d_log_antineutrino_loss_bd[1];
dLogNuLoss[2] = derivatives.d_log_antineutrino_loss_bd[2];
break;
}
@@ -583,12 +605,10 @@ namespace gridfire::rates::weak {
// Contributions from the reaction rate (output 0)
px[0] = py[0] * forwardPassRate * ln10 * dLogRate[0];
px[1] = py[0] * forwardPassRate * ln10 * dLogRate[1];
px[2] = py[0] * forwardPassRate * ln10 * dLogRate[2];
// Contributions from the neutrino loss rate (output 1)
px[0] += py[1] * forwardPassNeutrinoLossRate * ln10 * dLogNuLoss[0];
px[1] += py[1] * forwardPassNeutrinoLossRate * ln10 * dLogNuLoss[1];
px[2] += py[1] * forwardPassNeutrinoLossRate * ln10 * dLogNuLoss[2];
return true;
@@ -602,7 +622,6 @@ namespace gridfire::rates::weak {
std::set<size_t> all_input_deps;
all_input_deps.insert(r[0].begin(), r[0].end());
all_input_deps.insert(r[1].begin(), r[1].end());
all_input_deps.insert(r[2].begin(), r[2].end());
// What this is saying is that both output variables depend on all input variables.
s[0] = all_input_deps;
@@ -623,7 +642,6 @@ namespace gridfire::rates::weak {
st[0] = all_output_deps;
st[1] = all_output_deps;
st[2] = all_output_deps;
return true;
}

256
src/lib/reaction/weak/weak_interpolator.cpp
View File

@@ -13,54 +13,61 @@
#include "fourdst/composition/species.h"
#include "quill/LogMacros.h"
namespace gridfire::rates::weak {
WeakRateInterpolator::WeakRateInterpolator(const RowDataTable &raw_data) {
// Group all raw data rows by their isotope ID.
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);
}
// Process each isotope's data to build a simple 2D grid.
for (auto const& [isotope_id, rows] : grouped_rows) {
IsotopeGrid grid;
std::set<float> unique_t9, unique_rhoYe, unique_mue;
// Establish the T9 and log(rho*Ye) axes
std::set<float> unique_t9, unique_rhoYe;
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);
grid.t9_axis.assign(unique_t9.begin(), unique_t9.end());
grid.rhoYe_axis.assign(unique_rhoYe.begin(), unique_rhoYe.end());
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);
grid.data.resize(nt9 * nrhoYe);
// Reverse map for quick index lookup
std::unordered_map<float, size_t> t9_map, rhoYe_map, mue_map;
// Create reverse maps for efficient index lookups.
std::unordered_map<double, size_t> t9_map, rhoYe_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; }
// Use a set to detect duplicate (T9, rhoYe) pairs, which would be a data error.
std::set<std::pair<float, float>> seen_coords;
// Populate the 2D grid.
for (const auto* row: rows) {
if (auto [it, inserted] = seen_coords.insert({row->t9, row->log_rhoye}); !inserted) {
auto A = static_cast<uint16_t>(isotope_id >> 8);
auto Z = static_cast<uint8_t>(isotope_id & 0xFF);
std::string msg = std::format(
"Duplicate data point for isotope (A={}, Z={}) at (T9={}, log10(rho*Ye)={}) in weak rate table. This indicates corrupted or malformed input data and should be taken as an unrecoverable error.",
A, Z, row->t9, row->log_rhoye
);
LOG_ERROR(m_logger, "{}", msg);
throw std::runtime_error(msg);
}
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;
size_t index = i_t9 * nrhoYe + j_rhoYe;
grid.data[index] = WeakRatePayload{
row->log_beta_plus,
row->log_electron_capture,
@@ -74,35 +81,41 @@ namespace gridfire::rates::weak {
}
}
std::vector<fourdst::atomic::Species> WeakRateInterpolator::available_isotopes() const {
std::vector<fourdst::atomic::Species> isotopes;
using namespace fourdst::atomic;
std::vector<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());
const auto A = static_cast<uint16_t>(packed_id >> 8);
const auto Z = static_cast<uint8_t>(packed_id & 0xFF);
std::expected<Species, SpeciesErrorType> result = az_to_species(A, Z);
if (!result.has_value()) {
std::string msg = "Could not convert A=" + std::to_string(A) + ", Z=" + std::to_string(Z) + " to Species: ";
msg += (result.error() == SpeciesErrorType::ELEMENT_SYMBOL_NOT_FOUND) ? "Unknown element (Z out of range)." : "Invalid isotope (A < Z or A out of range).";
LOG_TRACE_L3(m_logger, "{}", msg);
} else {
isotopes.emplace_back(result.value());
}
}
return isotopes;
}
std::expected<WeakRatePayload, InterpolationError> WeakRateInterpolator::get_rates(
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 double log_rhoYe
) 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;
const auto& grid = it->second;
const auto& t9_axis = grid.t9_axis;
const auto& rhoYe_axis = grid.rhoYe_axis;
// Now find the bracketing indices for t9, log_rhoYe, and mu_e
// Find bracketing indices for the 2D (t9, rhoYe) grid
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()) {
@@ -113,161 +126,79 @@ namespace gridfire::rates::weak {
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) {
// Handle bounds errors for the 2D grid
if (!i_t9_opt || !j_rhoYe_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 (!i_t9_opt.has_value()) {
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 (!j_rhoYe_opt.has_value()) {
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
}
);
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();
const size_t nrhoYe = rhoYe_axis.size();
// 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];
// Get the four corner payloads for the bilinear interpolation
const auto& p00 = grid.data[(i * nrhoYe) + j];
const auto& p01 = grid.data[(i * nrhoYe) + j + 1];
const auto& p10 = grid.data[((i + 1) * nrhoYe) + j];
const auto& p11 = grid.data[((i + 1) * nrhoYe) + j + 1];
const double td = (t9 - t1) / (t2 - t1);
const double rd = (log_rhoYe - r1) / (r2 - r1);
const double md = (mu_e - m1) / (m2 - m1);
// Fractional distances for the 2D bilinear interpolation
const double td = (t9 - t9_axis[i]) / (t9_axis[i + 1] - t9_axis[i]);
const double rd = (log_rhoYe - rhoYe_axis[j]) / (rhoYe_axis[j + 1] - rhoYe_axis[j]);
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]);
// Helper lambda to linearly interpolate between two full payloads
auto lerp_payload = [](const WeakRatePayload& p0, const WeakRatePayload& p1, double t) {
return WeakRatePayload{
.log_beta_plus = std::lerp(p0.log_beta_plus, p1.log_beta_plus, t),
.log_electron_capture = std::lerp(p0.log_electron_capture, p1.log_electron_capture, t),
.log_neutrino_loss_ec = std::lerp(p0.log_neutrino_loss_ec, p1.log_neutrino_loss_ec, t),
.log_beta_minus = std::lerp(p0.log_beta_minus, p1.log_beta_minus, t),
.log_positron_capture = std::lerp(p0.log_positron_capture, p1.log_positron_capture, t),
.log_antineutrino_loss_bd = std::lerp(p0.log_antineutrino_loss_bd, p1.log_antineutrino_loss_bd, t),
};
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;
// Perform the bilinear interpolation
const WeakRatePayload p0 = lerp_payload(p00, p01, rd);
const WeakRatePayload p1 = lerp_payload(p10, p11, rd);
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;
return lerp_payload(p0, p1, td);
}
std::expected<WeakRateDerivatives, InterpolationError> WeakRateInterpolator::get_rate_derivatives(
uint16_t A,
uint8_t Z,
double t9,
double log_rhoYe,
double mu_e
double log_rhoYe
) const {
WeakRateDerivatives result;
WeakRateDerivatives result{};
//TODO: Make this perturbation scale aware
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 auto payload_plus_t9 = get_rates(A, Z, t9 + h_t9, log_rhoYe);
const auto payload_minus_t9 = get_rates(A, Z, t9 - h_t9, log_rhoYe);
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 auto payload_plus_rhoYe = get_rates(A, Z, t9, log_rhoYe + h_rhoYe);
const auto payload_minus_rhoYe = get_rates(A, Z, t9, log_rhoYe - h_rhoYe);
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);
if (!payload_plus_t9 || !payload_minus_t9 || !payload_plus_rhoYe || !payload_minus_rhoYe) {
// Determine which perturbation failed and return a consolidated error
auto first_error = !payload_plus_t9 ? payload_plus_t9.error() :
!payload_minus_t9 ? payload_minus_t9.error() :
!payload_plus_rhoYe ? payload_plus_rhoYe.error() :
payload_minus_rhoYe.error();
return std::unexpected(first_error);
}
// Derivatives wrt. T9
@@ -288,15 +219,6 @@ namespace gridfire::rates::weak {
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;
}

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

@@ -393,8 +393,8 @@ namespace gridfire::solver {
void CVODESolverStrategy::calculate_rhs(
const sunrealtype t,
const N_Vector y,
const N_Vector ydot,
N_Vector y,
N_Vector ydot,
const CVODEUserData *data
) const {
const size_t numSpecies = m_engine.getNetworkSpecies().size();