From 3b8a0a1f33c7c6d443a91d53da46dd078f5f3c35 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 14 Oct 2025 13:37:48 -0400 Subject: [PATCH] fix(engine_multiscale): resolved a major species index ordering bug All jacobian calculations were broken because the indexing used to record the AD tape was broken (see not parallel to) the indexing used by the composition object. A fix for this was to sort the network species by mass. However, more generally we should introduce a mechanism to ensure these two indexed sets always remain parallel --- src/include/gridfire/engine/engine_graph.h | 11 +- .../gridfire/engine/procedures/priming.h | 11 +- .../gridfire/engine/views/engine_multiscale.h | 4 +- src/include/gridfire/reaction/reaction.h | 12 + src/lib/engine/engine_graph.cpp | 28 +- src/lib/engine/procedures/construction.cpp | 21 +- src/lib/engine/procedures/priming.cpp | 284 +++++++----------- src/lib/engine/views/engine_multiscale.cpp | 124 ++++++-- .../strategies/CVODE_solver_strategy.cpp | 2 + tests/graphnet_sandbox/main.cpp | 11 +- 10 files changed, 276 insertions(+), 232 deletions(-) diff --git a/src/include/gridfire/engine/engine_graph.h b/src/include/gridfire/engine/engine_graph.h index 71ce37e5..de8ef769 100644 --- a/src/include/gridfire/engine/engine_graph.h +++ b/src/include/gridfire/engine/engine_graph.h @@ -515,6 +515,7 @@ namespace gridfire { */ void rebuild(const fourdst::composition::Composition& comp, const BuildDepthType depth) override; + private: struct PrecomputedReaction { // Forward cacheing @@ -983,7 +984,6 @@ namespace gridfire { const T mue, const std::function(const fourdst::atomic::Species &)>& speciesIDLookup ) const { - // --- Pre-setup (flags to control conditionals in an AD safe / branch aware manner) --- // ----- Constants for AD safe calculations --- const T zero = static_cast(0.0); @@ -992,10 +992,10 @@ namespace gridfire { const T k_reaction = reaction.calculate_rate(T9, rho, Ye, mue, Y, m_indexToSpeciesMap); // --- Cound the number of each reactant species to account for species multiplicity --- - std::unordered_map reactant_counts; + std::unordered_map reactant_counts; reactant_counts.reserve(reaction.reactants().size()); for (const auto& reactant : reaction.reactants()) { - reactant_counts[std::string(reactant.name())]++; + reactant_counts[reactant]++; } const int totalReactants = static_cast(reaction.reactants().size()); @@ -1003,10 +1003,9 @@ namespace gridfire { auto molar_concentration_product = static_cast(1.0); // --- Loop through each unique reactant species and calculate the molar concentration for that species then multiply that into the accumulator --- - for (const auto& [species_name, count] : reactant_counts) { + for (const auto& [species, count] : reactant_counts) { // --- Resolve species to molar abundance --- - // TODO: We need some way to handle the case when a species in the reaction is not part of the composition being tracked - const std::optional species_index = speciesIDLookup(m_networkSpeciesMap.at(species_name)); + const std::optional species_index = speciesIDLookup(species); if (!species_index.has_value()) { return static_cast(0.0); // If any reactant is not present, the reaction cannot proceed } diff --git a/src/include/gridfire/engine/procedures/priming.h b/src/include/gridfire/engine/procedures/priming.h index 19a686da..a1efcea3 100644 --- a/src/include/gridfire/engine/procedures/priming.h +++ b/src/include/gridfire/engine/procedures/priming.h @@ -19,6 +19,7 @@ namespace gridfire { * * @param netIn Input network data containing initial composition, temperature, and density. * @param engine DynamicEngine used to build and evaluate the reaction network. + * @param ignoredReactionTypes Types of reactions to ignore during priming (e.g., weak reactions). * @pre netIn.composition defines species and their mass fractions; engine is constructed with a valid network. * @post engine.networkReactions restored to its initial state; returned report contains primedComposition, * massFractionChanges for each species, success flag, and status code. @@ -26,7 +27,8 @@ namespace gridfire { */ PrimingReport primeNetwork( const NetIn& netIn, - DynamicEngine& engine + DynamicEngine& engine, + const std::optional>& ignoredReactionTypes ); /** @@ -40,6 +42,7 @@ namespace gridfire { * @param composition Current composition providing abundances for all species. * @param T9 Temperature in units of 10^9 K. * @param rho Density of the medium. + * @param reactionTypesToIgnore types of reactions to ignore during calculation. * @pre Y.size() matches engine.getNetworkReactions().size() mapping species order. * @post Returned rate constant is non-negative. * @return Sum of absolute stoichiometry-weighted destruction flows for the species. @@ -49,7 +52,8 @@ namespace gridfire { const fourdst::atomic::Species& species, const fourdst::composition::Composition& composition, double T9, - double rho + double rho, const std::optional> & + reactionTypesToIgnore ); /** @@ -63,6 +67,7 @@ namespace gridfire { * @param composition Composition object containing current abundances. * @param T9 Temperature in units of 10^9 K. * @param rho Density of the medium. + * @param reactionTypesToIgnore types of reactions to ignore during calculation. * @pre Y.size() matches engine.getNetworkReactions().size() mapping species order. * @post Returned creation rate is non-negative. * @return Sum of stoichiometry-weighted creation flows for the species. @@ -72,6 +77,6 @@ namespace gridfire { const fourdst::atomic::Species& species, const fourdst::composition::Composition& composition, double T9, - double rho + double rho, const std::optional> &reactionTypesToIgnore ); } \ No newline at end of file diff --git a/src/include/gridfire/engine/views/engine_multiscale.h b/src/include/gridfire/engine/views/engine_multiscale.h index 7c949660..9f428891 100644 --- a/src/include/gridfire/engine/views/engine_multiscale.h +++ b/src/include/gridfire/engine/views/engine_multiscale.h @@ -193,7 +193,7 @@ namespace gridfire { * It must be a `GraphEngine` and not a more general `DynamicEngine` * because this view relies on its specific implementation details. */ - explicit MultiscalePartitioningEngineView(GraphEngine& baseEngine); + explicit MultiscalePartitioningEngineView(DynamicEngine& baseEngine); /** * @brief Gets the list of species in the network. @@ -977,7 +977,7 @@ namespace gridfire { /** * @brief The base engine to which this view delegates calculations. */ - GraphEngine& m_baseEngine; + DynamicEngine& m_baseEngine; /** * @brief The list of identified equilibrium groups. */ diff --git a/src/include/gridfire/reaction/reaction.h b/src/include/gridfire/reaction/reaction.h index 18305ea1..e7c63fd2 100644 --- a/src/include/gridfire/reaction/reaction.h +++ b/src/include/gridfire/reaction/reaction.h @@ -28,6 +28,18 @@ namespace gridfire::reaction { REACLIB, LOGICAL_REACLIB, }; + + static std::unordered_map ReactionTypeNames = { + {ReactionType::WEAK, "weak"}, + {ReactionType::REACLIB, "reaclib"}, + {ReactionType::LOGICAL_REACLIB, "logical_reaclib"}, + }; + + static std::unordered_map ReactionPhysicalTypeNames = { + {ReactionType::WEAK, "Weak"}, + {ReactionType::REACLIB, "Strong"}, + {ReactionType::LOGICAL_REACLIB, "Strong"}, + }; /** * @struct RateCoefficientSet * @brief Holds the seven coefficients for the REACLIB rate equation. diff --git a/src/lib/engine/engine_graph.cpp b/src/lib/engine/engine_graph.cpp index 431a171e..755cf54c 100644 --- a/src/lib/engine/engine_graph.cpp +++ b/src/lib/engine/engine_graph.cpp @@ -192,7 +192,9 @@ namespace gridfire { throw std::runtime_error("Species not found in global atomic species database: " + std::string(name)); } } - + std::ranges::sort(m_networkSpecies, [](const fourdst::atomic::Species& a, const fourdst::atomic::Species& b) -> bool { + return a.mass() < b.mass(); // Otherwise, sort by mass + }); } void GraphEngine::populateReactionIDMap() { @@ -508,12 +510,21 @@ namespace gridfire { composition.setMassFraction(symbol, 0.0); } } - composition.finalize(true); + const bool didFinalize = composition.finalize(true); + if (!didFinalize) { + LOG_ERROR(m_logger, "Failed to finalize composition during priming. Check input mass fractions for validity."); + throw std::runtime_error("Failed to finalize composition during network priming."); + } fullNetIn.composition = composition; fullNetIn.temperature = netIn.temperature; fullNetIn.density = netIn.density; - auto primingReport = primeNetwork(fullNetIn, *this); + std::optional> reactionTypesToIgnore = std::nullopt; + if (!m_useReverseReactions) { + reactionTypesToIgnore = {reaction::ReactionType::WEAK}; + } + + auto primingReport = primeNetwork(fullNetIn, *this, reactionTypesToIgnore); return primingReport; } @@ -737,7 +748,6 @@ namespace gridfire { const double T9, const double rho ) const { - LOG_TRACE_L1_LIMIT_EVERY_N(1000, m_logger, "Generating jacobian matrix for T9={}, rho={}..", T9, rho); const size_t numSpecies = m_networkSpecies.size(); @@ -760,6 +770,7 @@ namespace gridfire { const double value = dotY[i * (numSpecies + 2) + j]; if (std::abs(value) > MIN_JACOBIAN_THRESHOLD || i == j) { // Always keep diagonal elements to avoid pathological stiffness m_jacobianMatrix(i, j) = value; + } } } @@ -1023,7 +1034,11 @@ namespace gridfire { baseUpdatedComposition.setMassFraction(species, 0.0); } } - baseUpdatedComposition.finalize(false); + const bool didFinalize = baseUpdatedComposition.finalize(false); + if (!didFinalize) { + LOG_ERROR(m_logger, "Failed to finalize composition during update. Check input mass fractions for validity."); + throw std::runtime_error("Failed to finalize composition during network update."); + } return baseUpdatedComposition; } @@ -1079,10 +1094,11 @@ namespace gridfire { adYe, adMue, [&](const fourdst::atomic::Species& querySpecies) -> size_t { - return m_speciesToIndexMap.at(querySpecies); // TODO: This is bad, needs to be fixed + return m_speciesToIndexMap.at(querySpecies); } ); + // Extract the raw vector from the associative map std::vector> dydt_vec; dydt_vec.reserve(dydt.size()); diff --git a/src/lib/engine/procedures/construction.cpp b/src/lib/engine/procedures/construction.cpp index 544d6c72..1c90fcee 100644 --- a/src/lib/engine/procedures/construction.cpp +++ b/src/lib/engine/procedures/construction.cpp @@ -120,6 +120,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; LOG_TRACE_L1(logger, "Collecting reactions for layer {} with {} remaining reactions. Currently there are {} available species", layer, remainingReactions.size(), availableSpecies.size()); std::vector reactionsForNextPass; std::unordered_set newProductsThisLayer; @@ -138,6 +141,12 @@ namespace gridfire { if (allReactantsAvailable) { collectedReactionPtrs.push_back(reaction); + if (reaction->type() == reaction::ReactionType::WEAK) { + collectedWeak++; + } else { + collectedStrong++; + } + collectedThisLayer++; newReactionsAdded = true; for (const auto& product : reaction->products()) { @@ -153,14 +162,18 @@ namespace gridfire { break; } + size_t oldProductCount = availableSpecies.size(); + availableSpecies.insert(newProductsThisLayer.begin(), newProductsThisLayer.end()); + size_t newProductCount = availableSpecies.size() - oldProductCount; LOG_TRACE_L1( logger, - "Layer {}: Collected {} new reactions. New products this layer: {}", + "Layer {}: Collected {} new reactions ({} strong, {} weak). New products this layer: {}", layer, - collectedReactionPtrs.size() - collectedReactionPtrs.size(), - newProductsThisLayer.size() + collectedThisLayer, + collectedStrong, + collectedWeak, + newProductCount ); - availableSpecies.insert(newProductsThisLayer.begin(), newProductsThisLayer.end()); remainingReactions = std::move(reactionsForNextPass); } diff --git a/src/lib/engine/procedures/priming.cpp b/src/lib/engine/procedures/priming.cpp index 93ae9664..1436ce73 100644 --- a/src/lib/engine/procedures/priming.cpp +++ b/src/lib/engine/procedures/priming.cpp @@ -10,6 +10,46 @@ #include "quill/Logger.h" #include "quill/LogMacros.h" +namespace { + // Create a dummy wrapper Composition to measure the unrestricted flow rate of species + class UnrestrictedComposition final : public fourdst::composition::Composition { + private: + const fourdst::atomic::Species& m_unrestrictedSpecies; + public: + explicit UnrestrictedComposition(const Composition& baseComposition, const fourdst::atomic::Species& unrestrictedSpecies) : + Composition(baseComposition), + m_unrestrictedSpecies(unrestrictedSpecies) {} + + double getMolarAbundance(const fourdst::atomic::Species &species) const override { + if (species == m_unrestrictedSpecies) { + return 1.0; // Set to a high value to simulate unrestricted abundance + } + return Composition::getMolarAbundance(species); + } + + double getMolarAbundance(const std::string &symbol) const override { + if (symbol == m_unrestrictedSpecies.name()) { + return 1.0; // Set to a high value to simulate unrestricted abundance + } + return Composition::getMolarAbundance(symbol); + } + }; + + bool isReactionIgnorable( + const gridfire::reaction::Reaction& reaction, + const std::optional>& reactionsTypesToIgnore + ) { + if (reactionsTypesToIgnore.has_value()) { + for (const auto& type : reactionsTypesToIgnore.value()) { + if (reaction.type() == type) { + return true; + } + } + } + return false; + } +} + namespace gridfire { using fourdst::composition::Composition; using fourdst::atomic::Species; @@ -19,11 +59,14 @@ namespace gridfire { const Species& species, const Composition &comp, const double T9, - const double rho + const double rho, + const std::optional> &reactionsTypesToIgnore ) { const reaction::Reaction* dominateReaction = nullptr; double maxFlow = -1.0; for (const auto& reaction : engine.getNetworkReactions()) { + if (isReactionIgnorable(*reaction, reactionsTypesToIgnore)) continue; + if (reaction->contains(species) && reaction->stoichiometry(species) > 0) { const double flow = engine.calculateMolarReactionFlow(*reaction, comp, T9, rho); if (flow > maxFlow) { @@ -63,10 +106,15 @@ namespace gridfire { * * @param netIn Input network data containing initial composition, temperature, and density. * @param engine DynamicEngine used to build and evaluate the reaction network. + * @param ignoredReactionTypes Types of reactions to ignore during priming (e.g., weak reactions). * @return PrimingReport encapsulating the results of the priming operation, including the new * robustly primed composition. */ - PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) { + PrimingReport primeNetwork( + const NetIn& netIn, + DynamicEngine& engine, + const std::optional>& ignoredReactionTypes + ) { auto logger = fourdst::logging::LogManager::getInstance().getLogger("log"); // --- Initial Setup --- @@ -77,6 +125,12 @@ namespace gridfire { speciesToPrime.push_back(entry.isotope()); } } + + // sort primingSpecies by mass number, lightest to heaviest. This ensures we prime in a physically sensible order. + std::ranges::sort(speciesToPrime, [](const Species& a, const Species& b) { + return a.mass() < b.mass(); + }); + LOG_DEBUG(logger, "Priming {} species in the network.", speciesToPrime.size()); // If no species need priming, return immediately. @@ -125,7 +179,13 @@ namespace gridfire { tempComp.registerSymbol(std::string(sp.name())); tempComp.setMassFraction(sp, std::max(0.0, mf)); } - tempComp.finalize(true); + bool didFinalize = tempComp.finalize(true); + if (!didFinalize) { + LOG_ERROR(logger, "Failed to finalize temporary composition during priming."); + report.success = false; + report.status = PrimingReportStatus::FAILED_TO_FINALIZE_COMPOSITION; + continue; + } NetworkPrimingEngineView primer(primingSpecies, engine); if (primer.getNetworkReactions().size() == 0) { @@ -135,14 +195,37 @@ namespace gridfire { continue; } - const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, tempComp, T9, rho); + const double destructionRateConstant = calculateDestructionRateConstant( + primer, + primingSpecies, + tempComp, + T9, + rho, + ignoredReactionTypes + ); if (destructionRateConstant > 1e-99) { - const double creationRate = calculateCreationRate(primer, primingSpecies, tempComp, T9, rho); + const double creationRate = calculateCreationRate( + primer, + primingSpecies, + tempComp, + T9, + rho, + ignoredReactionTypes + ); double equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass(); + + // ReSharper disable once CppDFAUnusedValue if (std::isnan(equilibriumMassFraction)) equilibriumMassFraction = 0.0; - if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, tempComp, T9, rho)) { + if (const reaction::Reaction* dominantChannel = findDominantCreationChannel( + primer, + primingSpecies, + tempComp, + T9, + rho, + ignoredReactionTypes) + ) { // Store the request instead of applying it immediately. requests.push_back({primingSpecies, equilibriumMassFraction, dominantChannel->reactants()}); } else { @@ -242,174 +325,22 @@ namespace gridfire { return report; } - - - // PrimingReport primeNetwork(const NetIn& netIn, DynamicEngine& engine) { - // auto logger = LogManager::getInstance().getLogger("log"); - // - // std::vector speciesToPrime; - // for (const auto &entry: netIn.composition | std::views::values) { - // std::cout << "Checking species: " << entry.isotope().name() << " with mass fraction: " << entry.mass_fraction() << std::endl; - // if (entry.mass_fraction() == 0.0) { - // speciesToPrime.push_back(entry.isotope()); - // } - // } - // LOG_DEBUG(logger, "Priming {} species in the network.", speciesToPrime.size()); - // - // PrimingReport report; - // if (speciesToPrime.empty()) { - // report.primedComposition = netIn.composition; - // report.success = true; - // report.status = PrimingReportStatus::NO_SPECIES_TO_PRIME; - // return report; - // } - // - // const double T9 = netIn.temperature / 1e9; - // const double rho = netIn.density; - // const auto initialReactionSet = engine.getNetworkReactions(); - // - // report.status = PrimingReportStatus::FULL_SUCCESS; - // report.success = true; - // - // // --- 1: pack composition into internal map --- - // std::unordered_map currentMassFractions; - // for (const auto& entry : netIn.composition | std::views::values) { - // currentMassFractions[entry.isotope()] = entry.mass_fraction(); - // } - // for (const auto& entry : speciesToPrime) { - // currentMassFractions[entry] = 0.0; // Initialize priming species with 0 mass fraction - // } - // - // std::unordered_map totalMassFractionChanges; - // - // engine.rebuild(netIn.composition, NetworkBuildDepth::Full); - // - // for (const auto& primingSpecies : speciesToPrime) { - // LOG_TRACE_L3(logger, "Priming species: {}", primingSpecies.name()); - // - // // Create a temporary composition from the current internal state for the primer - // Composition tempComp; - // for(const auto& [sp, mf] : currentMassFractions) { - // tempComp.registerSymbol(std::string(sp.name())); - // if (mf < 0.0 && std::abs(mf) < 1e-16) { - // tempComp.setMassFraction(sp, 0.0); // Avoid negative mass fractions - // } else { - // tempComp.setMassFraction(sp, mf); - // } - // } - // tempComp.finalize(true); // Finalize with normalization - // - // 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()); - // report.success = false; - // report.status = PrimingReportStatus::FAILED_TO_FIND_PRIMING_REACTIONS; - // continue; - // } - // - // const auto Y = primer.mapNetInToMolarAbundanceVector(tempNetIn); - // const double destructionRateConstant = calculateDestructionRateConstant(primer, primingSpecies, Y, T9, rho); - // - // if (destructionRateConstant > 1e-99) { - // double equilibriumMassFraction = 0.0; - // const double creationRate = calculateCreationRate(primer, primingSpecies, Y, T9, rho); - // equilibriumMassFraction = (creationRate / destructionRateConstant) * primingSpecies.mass(); - // if (std::isnan(equilibriumMassFraction)) { - // LOG_WARNING(logger, "Equilibrium mass fraction for {} is NaN. Setting to 0.0. This is likely not an issue. It probably originates from all reactions leading to creation and destruction being frozen out. In that case 0.0 should be a good approximation. Hint: This happens often when the network temperature is very the low. ", primingSpecies.name()); - // equilibriumMassFraction = 0.0; - // } - // LOG_TRACE_L3(logger, "Found equilibrium for {}: X_eq = {:.4e}", primingSpecies.name(), equilibriumMassFraction); - // - // if (const reaction::Reaction* dominantChannel = findDominantCreationChannel(primer, primingSpecies, Y, T9, rho)) { - // LOG_TRACE_L3(logger, "Dominant creation channel for {}: {}", primingSpecies.name(), dominantChannel->id()); - // - // double totalReactantMass = 0.0; - // for (const auto& reactant : dominantChannel->reactants()) { - // totalReactantMass += reactant.mass(); - // } - // - // double scalingFactor = 1.0; - // for (const auto& reactant : dominantChannel->reactants()) { - // const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass); - // double availableMass = 0.0; - // if (currentMassFractions.contains(reactant)) { - // availableMass = currentMassFractions.at(reactant); - // } - // if (massToSubtract > availableMass && availableMass > 0) { - // scalingFactor = std::min(scalingFactor, availableMass / massToSubtract); - // } - // } - // - // if (scalingFactor < 1.0) { - // LOG_WARNING(logger, "Priming for {} was limited by reactant availability. Scaling transfer by {:.4e}", primingSpecies.name(), scalingFactor); - // equilibriumMassFraction *= scalingFactor; - // } - // - // // Update the internal mass fraction map and accumulate total changes - // totalMassFractionChanges[primingSpecies] += equilibriumMassFraction; - // currentMassFractions.at(primingSpecies) += equilibriumMassFraction; - // - // for (const auto& reactant : dominantChannel->reactants()) { - // const double massToSubtract = equilibriumMassFraction * (reactant.mass() / totalReactantMass); - // std::cout << "[Priming: " << primingSpecies.name() << ", Channel: " << dominantChannel->id() << "] Subtracting " << massToSubtract << " from reactant " << reactant.name() << std::endl; - // totalMassFractionChanges[reactant] -= massToSubtract; - // currentMassFractions[reactant] -= massToSubtract; - // } - // } else { - // LOG_ERROR(logger, "Failed to find dominant creation channel for {}.", primingSpecies.name()); - // report.status = PrimingReportStatus::FAILED_TO_FIND_CREATION_CHANNEL; - // totalMassFractionChanges[primingSpecies] += 1e-40; - // currentMassFractions.at(primingSpecies) += 1e-40; - // } - // } else { - // LOG_WARNING(logger, "No destruction channel found for {}. Using fallback abundance.", primingSpecies.name()); - // totalMassFractionChanges.at(primingSpecies) += 1e-40; - // currentMassFractions.at(primingSpecies) += 1e-40; - // report.status = PrimingReportStatus::BASE_NETWORK_TOO_SHALLOW; - // } - // } - // - // // --- Final Composition Construction --- - // std::vector final_symbols; - // std::vector final_mass_fractions; - // for(const auto& [species, mass_fraction] : currentMassFractions) { - // final_symbols.emplace_back(species.name()); - // if (mass_fraction < 0.0 && std::abs(mass_fraction) < 1e-16) { - // final_mass_fractions.push_back(0.0); // Avoid negative mass fractions - // } else { - // final_mass_fractions.push_back(mass_fraction); - // } - // } - // - // // Create the final composition object from the pre-normalized mass fractions - // Composition primedComposition(final_symbols, final_mass_fractions, true); - // - // report.primedComposition = primedComposition; - // for (const auto& [species, change] : totalMassFractionChanges) { - // report.massFractionChanges.emplace_back(species, change); - // } - // - // engine.setNetworkReactions(initialReactionSet); - // return report; - // } - double calculateDestructionRateConstant( - const DynamicEngine& engine, - const Species& species, - const Composition& comp, - const double T9, - const double rho - ) { - //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. + const DynamicEngine& engine, + const Species& species, + const Composition& composition, + const double T9, + const double rho, + const std::optional> &reactionTypesToIgnore + ) { + const UnrestrictedComposition unrestrictedComp(composition, species); // Create a composition that simulates an enormous source abundance of the target species (getMolarAbundance(species) always returns 1.0) 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, comp, T9, rho); + if (isReactionIgnorable(*reaction, reactionTypesToIgnore)) continue; + + const int stoichiometry = reaction->stoichiometry(species); + if (stoichiometry < 0) { + destructionRateConstant += std::abs(stoichiometry) * engine.calculateMolarReactionFlow(*reaction, unrestrictedComp, T9, rho); } } return destructionRateConstant; @@ -418,17 +349,18 @@ namespace gridfire { double calculateCreationRate( const DynamicEngine& engine, const Species& species, - const Composition& comp, + const Composition& composition, const double T9, - const double rho + const double rho, + const std::optional> &reactionTypesToIgnore ) { double creationRate = 0.0; for (const auto& reaction: engine.getNetworkReactions()) { + if (isReactionIgnorable(*reaction, reactionTypesToIgnore)) continue; + const int stoichiometry = reaction->stoichiometry(species); if (stoichiometry > 0) { - if (engine.calculateMolarReactionFlow(*reaction, comp, T9, rho) > 0.0) { - } - creationRate += stoichiometry * engine.calculateMolarReactionFlow(*reaction, comp, T9, rho); + creationRate += stoichiometry * engine.calculateMolarReactionFlow(*reaction, composition, T9, rho); } } return creationRate; diff --git a/src/lib/engine/views/engine_multiscale.cpp b/src/lib/engine/views/engine_multiscale.cpp index 6b7721d8..90372ebc 100644 --- a/src/lib/engine/views/engine_multiscale.cpp +++ b/src/lib/engine/views/engine_multiscale.cpp @@ -22,7 +22,7 @@ namespace { // guarantee that this function will return the same ordering as the canonical vector representation) std::vector packCompositionToVector( const fourdst::composition::Composition& composition, - const gridfire::GraphEngine& engine + const gridfire::DynamicEngine& engine ) { std::vector Y(engine.getNetworkSpecies().size(), 0.0); const auto& allSpecies = engine.getNetworkSpecies(); @@ -172,7 +172,7 @@ namespace gridfire { using fourdst::atomic::Species; MultiscalePartitioningEngineView::MultiscalePartitioningEngineView( - GraphEngine& baseEngine + DynamicEngine& baseEngine ) : m_baseEngine(baseEngine) {} const std::vector & MultiscalePartitioningEngineView::getNetworkSpecies() const { @@ -772,7 +772,12 @@ namespace gridfire { } } - qseComposition.finalize(true); + bool didFinalize = qseComposition.finalize(true); + if (!didFinalize) { + LOG_ERROR(m_logger, "Failed to finalize composition after solving QSE abundances."); + m_logger->flush_log(); + throw std::runtime_error("Failed to finalize composition after solving QSE abundances."); + } return qseComposition; } @@ -799,10 +804,10 @@ namespace gridfire { const double rho ) const { LOG_TRACE_L1(m_logger, "Partitioning by timescale..."); - const auto result= m_baseEngine.getSpeciesDestructionTimescales(comp, T9, rho); + const auto destructionTimescale= m_baseEngine.getSpeciesDestructionTimescales(comp, T9, rho); const auto netTimescale = m_baseEngine.getSpeciesTimescales(comp, T9, rho); - if (!result) { + if (!destructionTimescale) { LOG_ERROR(m_logger, "Failed to get species destruction timescales due to stale engine state"); m_logger->flush_log(); throw exceptions::StaleEngineError("Failed to get species destruction timescales due to stale engine state"); @@ -812,26 +817,30 @@ namespace gridfire { m_logger->flush_log(); throw exceptions::StaleEngineError("Failed to get net species timescales due to stale engine state"); } - const std::unordered_map& all_timescales = result.value(); + const std::unordered_map& destruction_timescales = destructionTimescale.value(); const std::unordered_map& 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); + } + const auto& all_species = m_baseEngine.getNetworkSpecies(); - std::vector> sorted_timescales; - for (const auto & all_specie : all_species) { - double timescale = all_timescales.at(all_specie); - double net_timescale = net_timescales.at(all_specie); - if (std::isfinite(timescale) && timescale > 0) { - LOG_TRACE_L3(m_logger, "Species {} has finite destruction timescale: destruction: {} s, net: {} s", all_specie.name(), timescale, net_timescale); - sorted_timescales.emplace_back(timescale, all_specie); + std::vector> 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); + 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", all_specie.name(), timescale, net_timescale); + LOG_TRACE_L3(m_logger, "Species {} has infinite or negative destruction timescale: destruction: {} s, net: {} s", species.name(), destruction_timescale, net_timescale); } } std::ranges::sort( - sorted_timescales, + sorted_destruction_timescales, [](const auto& a, const auto& b) { return a.first > b.first; @@ -839,31 +848,45 @@ namespace gridfire { ); std::vector> final_pools; - if (sorted_timescales.empty()) { + if (sorted_destruction_timescales.empty()) { return final_pools; } constexpr double ABSOLUTE_QSE_TIMESCALE_THRESHOLD = 3.156e7; // Absolute threshold for QSE timescale (1 yr) constexpr double MIN_GAP_THRESHOLD = 2.0; // Require a 2 order of magnitude gap + constexpr double MIN_MOLAR_ABUNDANCE_THRESHOLD = 1e-10; // Minimum abundance threshold to consider a species for QSE. Any species above this will always be considered dynamic. - LOG_TRACE_L1(m_logger, "Found {} species with finite timescales.", sorted_timescales.size()); + LOG_TRACE_L1(m_logger, "Found {} species with finite timescales.", sorted_destruction_timescales.size()); LOG_TRACE_L1(m_logger, "Absolute QSE timescale threshold: {} seconds ({} years).", ABSOLUTE_QSE_TIMESCALE_THRESHOLD, ABSOLUTE_QSE_TIMESCALE_THRESHOLD / 3.156e7); LOG_TRACE_L1(m_logger, "Minimum gap threshold: {} orders of magnitude.", MIN_GAP_THRESHOLD); + LOG_TRACE_L1(m_logger, "Minimum molar abundance threshold: {}.", MIN_MOLAR_ABUNDANCE_THRESHOLD); std::vector dynamic_pool_species; std::vector> fast_candidates; // 1. First Pass: Absolute Timescale Cutoff - for (const auto& [timescale, species] : sorted_timescales) { - if (timescale > ABSOLUTE_QSE_TIMESCALE_THRESHOLD) { + for (const auto& [destruction_timescale, species] : sorted_destruction_timescales) { + if (species == n_1) { + LOG_TRACE_L3(m_logger, "Skipping neutron (n) from timescale analysis. Neutrons are always considered dynamic due to their extremely high connectivity."); + dynamic_pool_species.push_back(species); + continue; + } + if (destruction_timescale > ABSOLUTE_QSE_TIMESCALE_THRESHOLD) { LOG_TRACE_L3(m_logger, "Species {} with timescale {} is considered dynamic (slower than qse timescale threshold).", - species.name(), timescale); + species.name(), destruction_timescale); dynamic_pool_species.push_back(species); } else { - LOG_TRACE_L3(m_logger, "Species {} with timescale {} is a candidate fast species (faster than qse timescale threshold).", - species.name(), timescale); - fast_candidates.emplace_back(timescale, species); + const double Yi = comp.getMolarAbundance(species); + if (Yi > MIN_MOLAR_ABUNDANCE_THRESHOLD) { + LOG_TRACE_L3(m_logger, "Species {} with abundance {} is considered dynamic (above minimum abundance threshold of {}).", + species.name(), Yi, MIN_MOLAR_ABUNDANCE_THRESHOLD); + dynamic_pool_species.push_back(species); + continue; + } + LOG_TRACE_L3(m_logger, "Species {} with timescale {} and molar abundance {} is a candidate fast species (faster than qse timescale threshold and less than the molar abundance threshold).", + species.name(), destruction_timescale, Yi); + fast_candidates.emplace_back(destruction_timescale, species); } } @@ -939,11 +962,13 @@ namespace gridfire { const fourdst::composition::Composition &comp, const double T9, const double rho ) const { - constexpr double FLUX_RATIO_THRESHOLD = 5; std::vector validated_groups; std::vector invalidated_groups; 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 algebraic_group_members( group.algebraic_species.begin(), group.algebraic_species.end() @@ -954,9 +979,14 @@ namespace gridfire { group.seed_species.end() ); + // Values for measuring the flux coupling vs leakage 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) { @@ -967,7 +997,11 @@ namespace gridfire { for (const auto& reactant : reaction->reactants()) { if (algebraic_group_members.contains(reactant)) { 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; } + } bool has_internal_algebraic_product = false; @@ -975,6 +1009,9 @@ namespace gridfire { for (const auto& product : reaction->products()) { if (algebraic_group_members.contains(product)) { 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; } } @@ -1016,12 +1053,17 @@ namespace gridfire { leakage_flux += flow * leakage_fraction; coupling_flux += flow * coupling_fraction; + + } - if ((coupling_flux / leakage_flux ) > FLUX_RATIO_THRESHOLD) { + bool group_is_coupled = (coupling_flux / leakage_flux) > FLUX_RATIO_THRESHOLD; + bool group_is_balanced = std::abs(std::log(creationFlux) - std::log(destructionFlux)) < LOG_FLOW_RATIO_THRESHOLD; + + if (group_is_coupled) { LOG_TRACE_L1( m_logger, - "Group containing {} is in equilibrium due to high coupling flux threshold: leakage flux = {}, coupling flux = {}, ratio = {} (Threshold: {})", + "Group containing {} is in equilibrium due to high coupling flux and balanced creation and destruction rate: , ", [&]() -> std::string { std::stringstream ss; int count = 0; @@ -1037,14 +1079,18 @@ namespace gridfire { leakage_flux, coupling_flux, coupling_flux / leakage_flux, - FLUX_RATIO_THRESHOLD + FLUX_RATIO_THRESHOLD, + std::log10(creationFlux), + std::log10(destructionFlux), + std::abs(std::log10(creationFlux) - std::log10(destructionFlux)), + LOG_FLOW_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: leakage flux = {}, coupling flux = {}, ratio = {} (Threshold: {})", + "Group containing {} is NOT in equilibrium: , ", [&]() -> std::string { std::stringstream ss; int count = 0; @@ -1060,7 +1106,11 @@ namespace gridfire { leakage_flux, coupling_flux, coupling_flux / leakage_flux, - FLUX_RATIO_THRESHOLD + FLUX_RATIO_THRESHOLD, + std::log10(creationFlux), + std::log10(destructionFlux), + std::abs(std::log10(creationFlux) - std::log10(destructionFlux)), + LOG_FLOW_RATIO_THRESHOLD ); invalidated_groups.emplace_back(group); invalidated_groups.back().is_in_equilibrium = false; @@ -1156,13 +1206,19 @@ namespace gridfire { ); //TODO: Check this conversion double Xi = Y_final_qse(i) * species.mass(); // Convert from molar abundance to mass fraction - if (!outputComposition.contains(species)) { + if (!outputComposition.hasSpecies(species)) { outputComposition.registerSpecies(species); } outputComposition.setMassFraction(species, Xi); i++; } } + bool didFinalize = outputComposition.finalize(false); + if (!didFinalize) { + LOG_ERROR(m_logger, "Failed to finalize composition after solving QSE abundances."); + m_logger->flush_log(); + throw std::runtime_error("Failed to finalize composition after solving QSE abundances."); + } return outputComposition; } @@ -1374,7 +1430,10 @@ namespace gridfire { comp_trial.registerSpecies(species); } const double molarAbundance = y_qse[static_cast(m_qse_solve_species_index_map.at(species))]; - const double massFraction = molarAbundance * species.mass(); + double massFraction = molarAbundance * species.mass(); + if (massFraction < 0 && std::abs(massFraction) < 1e-20) { // if there is a larger negative mass fraction, let the composition module throw an error + massFraction = 0.0; // Avoid setting minuscule negative mass fractions due to numerical noise + } comp_trial.setMassFraction(species, massFraction); } @@ -1415,7 +1474,7 @@ namespace gridfire { 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__); + const std::string msg = std::format("Failed to finalize composition (comp_trial) in {} at line {}", __FILE__, __LINE__); throw std::runtime_error(msg); } @@ -1432,6 +1491,7 @@ namespace gridfire { colSpecies ); colID += 1; + LOG_TRACE_L3(m_view->m_logger, "Jacobian[{}, {}] (d(dY({}))/dY({})) = {}", rowID, colID - 1, rowSpecies.name(), colSpecies.name(), J_qse(rowID, colID - 1)); } rowID += 1; } diff --git a/src/lib/solver/strategies/CVODE_solver_strategy.cpp b/src/lib/solver/strategies/CVODE_solver_strategy.cpp index efeddfc6..a1e23b66 100644 --- a/src/lib/solver/strategies/CVODE_solver_strategy.cpp +++ b/src/lib/solver/strategies/CVODE_solver_strategy.cpp @@ -144,6 +144,8 @@ namespace gridfire::solver { const auto relTol = m_config.get("gridfire:solver:CVODESolverStrategy:relTol", 1.0e-8); fourdst::composition::Composition equilibratedComposition = m_engine.update(netIn); + std::cout << "EXITED AT EXPECTED TESTING POINT" << std::endl; + exit(0); size_t numSpecies = m_engine.getNetworkSpecies().size(); uint64_t N = numSpecies + 1; diff --git a/tests/graphnet_sandbox/main.cpp b/tests/graphnet_sandbox/main.cpp index 939b889f..b32f232d 100644 --- a/tests/graphnet_sandbox/main.cpp +++ b/tests/graphnet_sandbox/main.cpp @@ -23,6 +23,7 @@ #include #include +#include "gridfire/engine/views/engine_defined.h" #include "gridfire/partition/composite/partition_composite.h" static std::terminate_handler g_previousHandler = nullptr; @@ -81,7 +82,7 @@ int main(int argc, char* argv[]){ g_previousHandler = std::set_terminate(quill_terminate_handler); quill::Logger* logger = fourdst::logging::LogManager::getInstance().getLogger("log"); - logger->set_log_level(quill::LogLevel::TraceL2); + logger->set_log_level(quill::LogLevel::TraceL3); LOG_INFO(logger, "Starting Adaptive Engine View Example..."); using namespace gridfire; @@ -92,7 +93,11 @@ int main(int argc, char* argv[]){ fourdst::composition::Composition composition; composition.registerSymbol(symbols, true); composition.setMassFraction(symbols, comp); - composition.finalize(true); + bool didFinalize = composition.finalize(true); + if (!didFinalize) { + std::cerr << "Failed to finalize initial composition." << std::endl; + return 1; + } using partition::BasePartitionType; const auto partitionFunction = partition::CompositePartitionFunction({ BasePartitionType::RauscherThielemann, @@ -109,9 +114,9 @@ int main(int argc, char* argv[]){ netIn.dt0 = 1e-12; GraphEngine ReaclibEngine(composition, partitionFunction, NetworkBuildDepth::SecondOrder); - ReaclibEngine.setUseReverseReactions(false); ReaclibEngine.setPrecomputation(false); + // DefinedEngineView ppEngine({"p(p,e+)d", "d(p,g)he3", "he3(he3,2p)he4"}, ReaclibEngine); MultiscalePartitioningEngineView partitioningView(ReaclibEngine); AdaptiveEngineView adaptiveView(partitioningView);