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

feat(reverse-rates): fixed unit error in reverse rate calculation

Browse Source
This commit is contained in:
Emily Boudreaux
2025-07-29 07:49:01 -04:00
parent ba9b3e2392
commit 2d9216600f
2 changed files with 7 additions and 59 deletions
Download Patch File Download Diff File Expand all files Collapse all files

30
src/include/gridfire/engine/engine_graph.h
View File

@@ -36,7 +36,6 @@
// REACLIBReactions are quite large data structures, so this could be a performance bottleneck.
// static bool isF17 = false;
namespace gridfire {
static bool s_debug = false; // Global debug flag for the GraphEngine
/**
* @brief Alias for CppAD AD type for double precision.
*
@@ -548,24 +547,6 @@ namespace gridfire {
*/
[[nodiscard]] bool validateConservation() const;
/**
* @brief Validates the composition against the current reaction set.
*
* @param composition The composition to validate.
* @param culling The culling threshold to use.
* @param T9 The temperature to use.
*
* This method validates the composition against the current reaction set.
* If the composition is not compatible with the reaction set, the
* reaction set is rebuilt from the composition.
*/
void validateComposition(
const fourdst::composition::Composition &composition,
double culling,
double T9
);
[[nodiscard]] StepDerivatives<double> calculateAllDerivativesUsingPrecomputation(
const std::vector<double> &Y_in,
@@ -675,14 +656,6 @@ namespace gridfire {
if (!m_useReverseReactions) {
return static_cast<T>(0.0); // If reverse reactions are not used, return zero
}
s_debug = false;
if (reaction.peName() == "p(p,e+)d" || reaction.peName() =="d(d,n)he3" || reaction.peName() == "c12(p,g)n13") {
if constexpr (!std::is_same_v<T, ADDouble>) {
s_debug = true;
std::cout << "Calculating reverse molar flow for reaction: " << reaction.peName() << std::endl;
std::cout << "\tT9: " << T9 << ", rho: " << rho << std::endl;
}
}
T reverseMolarFlow = static_cast<T>(0.0);
if (reaction.qValue() != 0.0) {
@@ -702,9 +675,6 @@ namespace gridfire {
}
} else {
// A,B If not calling with an AD type, calculate the reverse rate directly
if (s_debug) {
std::cout << "\tUsing double overload\n";
}
reverseRateConstant = calculateReverseRate(reaction, T9);
}

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

@@ -239,26 +239,6 @@ namespace gridfire {
return true; // All reactions passed the conservation check
}
void GraphEngine::validateComposition(const fourdst::composition::Composition &composition, double culling, double T9) {
// Check if the requested network has already been cached.
// PERF: Rebuilding this should be pretty fast but it may be a good point of optimization in the future.
const reaction::LogicalReactionSet validationReactionSet = build_reaclib_nuclear_network(composition, false);
// TODO: need some more robust method here to
// A. Build the basic network from the composition's species with non zero mass fractions.
// B. rebuild a new composition from the reaction set's reactants + products (with the mass fractions from the things that are only products set to 0)
// C. Rebuild the reaction set from the new composition
// D. Cull reactions where all reactants have mass fractions below the culling threshold.
// E. Be careful about maintaining caching through all of this
// This allows for dynamic network modification while retaining caching for networks which are very similar.
if (validationReactionSet != m_reactions) {
LOG_DEBUG(m_logger, "Reaction set not cached. Rebuilding the reaction set for T9={} and culling={}.", T9, culling);
m_reactions = validationReactionSet;
syncInternalMaps(); // Re-sync internal maps after updating reactions. Note this will also retrace the AD tape.
}
}
double GraphEngine::calculateReverseRate(
const reaction::Reaction &reaction,
const double T9
@@ -267,15 +247,13 @@ namespace gridfire {
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 expFactor = std::exp(-reaction.qValue() / (m_constants.kB * T9));
if (s_debug) {
std::cout << "\texpFactor = exp(-qValue/(kB * T9))\n";
std::cout << "\texpFactor: " << expFactor << " for reaction: " << reaction.peName() << std::endl;
std::cout << "\tQ-value: " << reaction.qValue() << " at T9: " << T9 << std::endl;
std::cout << "\tT9: " << T9 << std::endl;
std::cout << "\tkB * T9: " << m_constants.kB * T9 << std::endl;
std::cout << "\tqValue/(kB * T9): " << reaction.qValue() / (m_constants.kB * T9) << std::endl;
}
const double temp = T9 * 1e9; // Convert T9 to Kelvin
// 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");
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);