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

feat(GraphEngine): More robust reaction type selection at network construction

Browse Source 
Added new ways to select exactly what types of reactions (strong, beta+, beta-, electron capture, positron capture, or any combination thereof) can be turned on at network construction time. There are a few quality of life masks added as well such as weak which addes all weak type reactions, and all which adds weak + strong reactions. The default is to just add strong reactions for now.
This commit is contained in:
Emily Boudreaux
2025-11-03 08:38:21 -05:00
parent 56f9342052
commit 86d7a283a1
13 changed files with 276 additions and 150 deletions
Download Patch File Download Diff File Expand all files Collapse all files

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

@@ -40,9 +40,16 @@ namespace gridfire {
GraphEngine::GraphEngine(
const fourdst::composition::Composition &composition,
const partition::PartitionFunction& partitionFunction,
const BuildDepthType buildDepth) :
const BuildDepthType buildDepth
) : GraphEngine(composition, partitionFunction, buildDepth, NetworkConstructionFlags::DEFAULT){}
GraphEngine::GraphEngine(
const fourdst::composition::Composition &composition,
const partition::PartitionFunction &partitionFunction,
const BuildDepthType buildDepth,
const NetworkConstructionFlags reactionTypes ) :
m_weakRateInterpolator(rates::weak::UNIFIED_WEAK_DATA),
m_reactions(build_nuclear_network(composition, m_weakRateInterpolator, buildDepth, false)),
m_reactions(build_nuclear_network(composition, m_weakRateInterpolator, buildDepth, reactionTypes)),
m_depth(buildDepth),
m_partitionFunction(partitionFunction.clone())
{
@@ -576,7 +583,7 @@ namespace gridfire {
void GraphEngine::rebuild(const fourdst::composition::Composition& comp, const BuildDepthType depth) {
if (depth != m_depth) {
m_depth = depth;
m_reactions = build_nuclear_network(comp, m_weakRateInterpolator, m_depth, false);
m_reactions = build_nuclear_network(comp, m_weakRateInterpolator, m_depth);
m_jacobianMatrixState = JacobianMatrixState::STALE;
syncInternalMaps(); // Resync internal maps after changing the depth
} else {

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

@@ -37,6 +37,86 @@ namespace {
return true;
}
gridfire::reaction::ReactionSet register_weak_reactions(
const std::optional<gridfire::rates::weak::WeakRateInterpolator> &weakInterpolator,
const gridfire::NetworkConstructionFlags reactionTypes
) {
gridfire::reaction::ReactionSet weak_reaction_pool;
assert(weakInterpolator.has_value());
if (!has_flag(reactionTypes, gridfire::NetworkConstructionFlags::WEAK)) {
return weak_reaction_pool;
}
for (const auto& parent_species: weakInterpolator->available_isotopes()) {
std::expected<fourdst::atomic::Species, fourdst::atomic::SpeciesErrorType> upProduct = fourdst::atomic::az_to_species(
parent_species.a(),
parent_species.z() + 1
);
std::expected<fourdst::atomic::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
if (has_flag(reactionTypes, gridfire::NetworkConstructionFlags::BETA_PLUS)) {
weak_reaction_pool.add_reaction(
std::make_unique<gridfire::rates::weak::WeakReaction>(
parent_species,
gridfire::rates::weak::WeakReactionType::BETA_PLUS_DECAY,
*weakInterpolator
)
);
}
if (has_flag(reactionTypes, gridfire::NetworkConstructionFlags::ELECTRON_CAPTURE)) {
weak_reaction_pool.add_reaction(
std::make_unique<gridfire::rates::weak::WeakReaction>(
parent_species,
gridfire::rates::weak::WeakReactionType::ELECTRON_CAPTURE,
*weakInterpolator
)
);
}
}
if (upProduct.has_value()) { // Only add the reaction if the Species map contains the product
if (has_flag(reactionTypes, gridfire::NetworkConstructionFlags::BETA_MINUS)) {
weak_reaction_pool.add_reaction(
std::make_unique<gridfire::rates::weak::WeakReaction>(
parent_species,
gridfire::rates::weak::WeakReactionType::BETA_MINUS_DECAY,
*weakInterpolator
)
);
}
if (has_flag(reactionTypes, gridfire::NetworkConstructionFlags::POSITRON_CAPTURE)) {
weak_reaction_pool.add_reaction(
std::make_unique<gridfire::rates::weak::WeakReaction>(
parent_species,
gridfire::rates::weak::WeakReactionType::POSITRON_CAPTURE,
*weakInterpolator
)
);
}
}
}
return weak_reaction_pool;
}
gridfire::reaction::ReactionSet register_strong_reactions(
const gridfire::NetworkConstructionFlags reaction_types
) {
gridfire::reaction::ReactionSet strong_reaction_pool;
if (has_flag(reaction_types, gridfire::NetworkConstructionFlags::STRONG)) {
const auto& allReaclibReactions = gridfire::reaclib::get_all_reaclib_reactions();
for (const auto& reaction : allReaclibReactions) {
if (!reaction->is_reverse() && !reaclib_reaction_is_weak(*reaction)) { // Only add reactions of the correct direction and which are not weak. Weak reactions are handled from the WRL separately which provides much higher quality weak reactions than reaclib does
strong_reaction_pool.add_reaction(reaction->clone());
}
}
}
return strong_reaction_pool;
}
}
namespace gridfire {
@@ -45,20 +125,26 @@ namespace gridfire {
using fourdst::composition::Composition;
using fourdst::atomic::Species;
ReactionSet build_nuclear_network(
const Composition& composition,
const rates::weak::WeakRateInterpolator& weakInterpolator,
const std::optional<rates::weak::WeakRateInterpolator> &weakInterpolator,
BuildDepthType maxLayers,
bool reverse
NetworkConstructionFlags ReactionTypes
) {
auto logger = fourdst::logging::LogManager::getInstance().getLogger("log");
LOG_INFO(logger, "Constructing network topology from reaction types : {}", NetworkConstructionFlagsToString(ReactionTypes));
if (ReactionTypes == NetworkConstructionFlags::NONE) {
LOG_ERROR(logger, "Reaction types is set to NONE. No reactions will be collected");
throw std::logic_error("Reaction types is set to NONE. No reactions will be collected");
}
int depth;
if (std::holds_alternative<NetworkBuildDepth>(maxLayers)) {
depth = static_cast<int>(std::get<NetworkBuildDepth>(maxLayers));
} else {
depth = std::get<int>(maxLayers);
}
auto logger = fourdst::logging::LogManager::getInstance().getLogger("log");
if (depth == 0) {
LOG_ERROR(logger, "Network build depth is set to 0. No reactions will be collected.");
throw std::logic_error("Network build depth is set to 0. No reactions will be collected.");
@@ -68,56 +154,10 @@ namespace gridfire {
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_reaction_is_weak(*reaction)) { // Only add reactions of the correct direction and which are not weak. Weak reactions are handled from the WRL separately which provides much higher quality weak reactions than reaclib does
master_reaction_pool.add_reaction(reaction->clone());
}
}
master_reaction_pool.extend(register_strong_reactions(ReactionTypes));
// --- Clone all possible weak reactions into the master reaction pool ---
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
);
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
)
);
}
} // TODO: Remove comments, weak reactions have been disabled for testing
// --- Clone all possible weak reactions into the master reaction pool if the construction function is told to use weak reactions ---
master_reaction_pool.extend(register_weak_reactions(weakInterpolator, ReactionTypes));
// --- Step 2: Use non-owning raw pointers for the fast build algorithm ---
std::vector<Reaction*> remainingReactions;
@@ -143,9 +183,9 @@ namespace gridfire {
LOG_INFO(logger, "Starting network construction with {} available species.", availableSpecies.size());
for (int layer = 0; layer < depth && !remainingReactions.empty(); ++layer) {
size_t collectedThisLayer = 0;
size_t collectedStrong = 0;
size_t collectedWeak = 0;
[[maybe_unused]] size_t collectedThisLayer = 0;
[[maybe_unused]] size_t collectedStrong = 0;
[[maybe_unused]] size_t collectedWeak = 0;
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;
std::unordered_set<Species> newProductsThisLayer;
@@ -185,9 +225,9 @@ namespace gridfire {
break;
}
size_t oldProductCount = availableSpecies.size();
[[maybe_unused]] size_t oldProductCount = availableSpecies.size();
availableSpecies.insert(newProductsThisLayer.begin(), newProductsThisLayer.end());
size_t newProductCount = availableSpecies.size() - oldProductCount;
[[maybe_unused]] size_t newProductCount = availableSpecies.size() - oldProductCount;
LOG_TRACE_L1(
logger,
"Layer {}: Collected {} new reactions ({} strong, {} weak). New products this layer: {}",

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

@@ -936,24 +936,30 @@ namespace gridfire {
throw exceptions::StaleEngineError("Failed to get net species timescales due to stale engine state");
}
const std::unordered_map<Species, double>& destruction_timescales = destructionTimescale.value();
const std::unordered_map<Species, double>& net_timescales = netTimescale.value();
[[maybe_unused]] const std::unordered_map<Species, double>& net_timescales = netTimescale.value();
for (const auto& [species, destruction_timescale] : destruction_timescales) {
LOG_TRACE_L3(m_logger, "For {} destruction timescale is {} s", species.name(), destruction_timescale);
}
LOG_TRACE_L3(
m_logger,
"{}",
[&]() -> std::string {
std::stringstream ss;
for (const auto& [species, destruction_timescale] : destruction_timescales) {
ss << std::format("For {} destruction timescale is {}s\n", species.name(), destruction_timescale);
}
return ss.str();
}()
);
const auto& all_species = m_baseEngine.getNetworkSpecies();
std::vector<std::pair<double, Species>> sorted_destruction_timescales;
for (const auto & species : all_species) {
double destruction_timescale = destruction_timescales.at(species);
double net_timescale = net_timescales.at(species);
if (std::isfinite(destruction_timescale) && destruction_timescale > 0) {
LOG_TRACE_L3(m_logger, "Species {} has finite destruction timescale: destruction: {} s, net: {} s", species.name(), destruction_timescale, net_timescale);
LOG_TRACE_L3(m_logger, "Species {} has finite destruction timescale: destruction: {} s, net: {} s", species.name(), destruction_timescale, net_timescales.at(species));
sorted_destruction_timescales.emplace_back(destruction_timescale, species);
} else {
LOG_TRACE_L3(m_logger, "Species {} has infinite or negative destruction timescale: destruction: {} s, net: {} s", species.name(), destruction_timescale, net_timescale);
LOG_TRACE_L3(m_logger, "Species {} has infinite or negative destruction timescale: destruction: {} s, net: {} s", species.name(), destruction_timescale, net_timescales.at(species));
}
}
@@ -1093,7 +1099,6 @@ namespace gridfire {
validated_groups.reserve(candidate_groups.size());
for (auto& group : candidate_groups) {
constexpr double FLUX_RATIO_THRESHOLD = 5;
constexpr double LOG_FLOW_RATIO_THRESHOLD = 2;
const std::unordered_set<Species> algebraic_group_members(
group.algebraic_species.begin(),
@@ -1109,10 +1114,6 @@ namespace gridfire {
double coupling_flux = 0.0;
double leakage_flux = 0.0;
// Values for validating if the group could physically be in equilibrium
double creationFlux = 0.0;
double destructionFlux = 0.0;
for (const auto& reaction: m_baseEngine.getNetworkReactions()) {
const double flow = std::abs(m_baseEngine.calculateMolarReactionFlow(*reaction, comp, T9, rho));
if (flow == 0.0) {
@@ -1125,7 +1126,6 @@ namespace gridfire {
has_internal_algebraic_reactant = true;
LOG_TRACE_L3(m_logger, "Adjusting destruction flux (+= {} mol g^-1 s^-1) for QSEGroup due to reactant {} from reaction {}",
flow, reactant.name(), reaction->id());
destructionFlux += flow;
}
}
@@ -1137,7 +1137,6 @@ namespace gridfire {
has_internal_algebraic_product = true;
LOG_TRACE_L3(m_logger, "Adjusting creation flux (+= {} mol g^-1 s^-1) for QSEGroup due to product {} from reaction {}",
flow, product.name(), reaction->id());
creationFlux += flow;
}
}
@@ -1186,7 +1185,7 @@ namespace gridfire {
if (coupling_flux / leakage_flux > FLUX_RATIO_THRESHOLD) {
LOG_TRACE_L1(
m_logger,
"Group containing {} is in equilibrium due to high coupling flux and balanced creation and destruction rate: <coupling: leakage flux = {}, coupling flux = {}, ratio = {} (Threshold: {})>, <creation: creation flux = {}, destruction flux = {}, ratio = {} order of mag (Threshold: {} order of mag)>",
"Group containing {} is in equilibrium due to high coupling flux and balanced creation and destruction rate: <coupling: leakage flux = {}, coupling flux = {}, ratio = {} (Threshold: {})>",
[&]() -> std::string {
std::stringstream ss;
int count = 0;
@@ -1202,18 +1201,14 @@ namespace gridfire {
leakage_flux,
coupling_flux,
coupling_flux / leakage_flux,
FLUX_RATIO_THRESHOLD,
std::log10(creationFlux),
std::log10(destructionFlux),
std::abs(std::log10(creationFlux) - std::log10(destructionFlux)),
LOG_FLOW_RATIO_THRESHOLD
FLUX_RATIO_THRESHOLD
);
validated_groups.emplace_back(group);
validated_groups.back().is_in_equilibrium = true;
} else {
LOG_TRACE_L1(
m_logger,
"Group containing {} is NOT in equilibrium: <coupling: leakage flux = {}, coupling flux = {}, ratio = {} (Threshold: {})>, <creation: creation flux = {}, destruction flux = {}, ratio = {} order of mag (Threshold: {} order of mag)>",
"Group containing {} is NOT in equilibrium: <coupling: leakage flux = {}, coupling flux = {}, ratio = {} (Threshold: {})>",
[&]() -> std::string {
std::stringstream ss;
int count = 0;
@@ -1229,11 +1224,7 @@ namespace gridfire {
leakage_flux,
coupling_flux,
coupling_flux / leakage_flux,
FLUX_RATIO_THRESHOLD,
std::log10(creationFlux),
std::log10(destructionFlux),
std::abs(std::log10(creationFlux) - std::log10(destructionFlux)),
LOG_FLOW_RATIO_THRESHOLD
FLUX_RATIO_THRESHOLD
);
invalidated_groups.emplace_back(group);
invalidated_groups.back().is_in_equilibrium = false;

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

@@ -397,6 +397,12 @@ namespace gridfire::reaction {
m_reactionNameMap.emplace(std::move(reaction_id), new_index);
}
void ReactionSet::extend(const ReactionSet &other) {
for (const auto& reaction : other.m_reactions) {
add_reaction(*reaction);
}
}
void ReactionSet::remove_reaction(const Reaction& reaction) {
const auto reaction_id = std::string(reaction.id());
if (!m_reactionNameMap.contains(reaction_id)) {

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

@@ -195,7 +195,8 @@ namespace gridfire::solver {
std::rethrow_exception(std::make_exception_ptr(*user_data.captured_exception));
}
// log_step_diagnostics(user_data, false);
// log_step_diagnostics(user_data, true);
// exit(0);
check_cvode_flag(flag, "CVode");
long int n_steps;