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);