diff --git a/src/include/gridfire/reaction/reaction.h b/src/include/gridfire/reaction/reaction.h index 18cef474..cb0401a0 100644 --- a/src/include/gridfire/reaction/reaction.h +++ b/src/include/gridfire/reaction/reaction.h @@ -1,5 +1,6 @@ #pragma once +#include #include #include "fourdst/composition/atomicSpecies.h" @@ -613,12 +614,22 @@ namespace gridfire::reaction { class LogicalReaclibReaction final : public ReaclibReaction { public: /** - * @brief Constructs a LogicalReaction from a vector of `Reaction` objects. + * @brief Constructs a LogicalReaction from a vector of `Reaction` objects. Implicitly assumes that the + * logical reaction is for a forward (i.e. not reverse) reaction. * @param reactions A vector of reactions that represent the same logical process. * @throws std::runtime_error if the provided reactions have inconsistent Q-values. */ explicit LogicalReaclibReaction(const std::vector &reactions); + /** + * @breif Constructs a LogicalReaction from a vector of `Reaction` objects and allows the user + * to specify if the logical set is for a reverse reaction explicitly + * @param reactions A vector of reactions that represent the same logical process + * @param reverse A flag to control if this logical reaction is reverse or not + * @returns std::runtime_error if the provided reactions have inconsistent Q-values. + */ + explicit LogicalReaclibReaction(const std::vector &reactions, bool reverse); + /** * @brief Adds another `Reaction` source to this logical reaction. * @param reaction The reaction to add. @@ -719,7 +730,7 @@ namespace gridfire::reaction { const T T953 = CppAD::pow(T9, 5.0/3.0); const T logT9 = CppAD::log(T9); // ReSharper disable once CppUseStructuredBinding - for (const auto& rate : m_rates) { + for (const auto& [rate, source] : std::views::zip(m_rates, m_sources)) { const T exponent = rate.a0 + rate.a1 / T9 + rate.a2 / T913 + diff --git a/src/lib/engine/engine_graph.cpp b/src/lib/engine/engine_graph.cpp index 1cc3194f..14777ff2 100644 --- a/src/lib/engine/engine_graph.cpp +++ b/src/lib/engine/engine_graph.cpp @@ -854,9 +854,18 @@ namespace gridfire { // --- Pack the input variables into a vector for CppAD --- const size_t numSpecies = m_networkSpecies.size(); std::vector x(numSpecies + 2, 0.0); - const std::vector& Y_dynamic = comp.getMolarAbundanceVector(); - for (size_t i = 0; i < numSpecies; ++i) { - x[i] = Y_dynamic[i]; + // const std::vector& Y_dynamic = comp.getMolarAbundanceVector(); + // for (size_t i = 0; i < numSpecies; ++i) { + // x[i] = Y_dynamic[i]; + // } + size_t i = 0; + for (const auto& species: m_networkSpecies) { + double Yi = 0.0; + if (comp.contains(species)) { + Yi = comp.getMolarAbundance(species); + } + x[i] = Yi; + i++; } x[numSpecies] = T9; x[numSpecies + 1] = rho; diff --git a/src/lib/engine/views/engine_multiscale.cpp b/src/lib/engine/views/engine_multiscale.cpp index d322aa5c..8f1ccd9d 100644 --- a/src/lib/engine/views/engine_multiscale.cpp +++ b/src/lib/engine/views/engine_multiscale.cpp @@ -225,9 +225,21 @@ namespace gridfire { const bool activeSpeciesIsSubset = std::ranges::any_of(activeSpecies, [&](const auto& species) -> bool { return !involvesSpecies(species); }); - if (activeSpeciesIsSubset) { - LOG_CRITICAL(m_logger, "Active species contains species not in the network partition. Cannot generate Jacobian matrix for active species."); - throw std::runtime_error("Active species contains species not in the network partition. Cannot generate Jacobian matrix for active species."); + if (!activeSpeciesIsSubset) { + std::string msg = std::format( + "Active species set contains species ({}) not present in network partition. Cannot generate jacobian matrix due to this.", + [&]() -> std::string { + std::stringstream ss; + for (const auto& species : activeSpecies) { + if (!this->involvesSpecies(species)) { + ss << species << " "; + } + } + return ss.str(); + }() + ); + LOG_CRITICAL(m_logger, "{}", msg); + throw std::runtime_error(msg); } std::vector dynamicActiveSpeciesIntersection; diff --git a/src/lib/reaction/reaclib.cpp b/src/lib/reaction/reaclib.cpp index 94bb50c9..d899d8c2 100644 --- a/src/lib/reaction/reaclib.cpp +++ b/src/lib/reaction/reaclib.cpp @@ -90,7 +90,6 @@ namespace gridfire::reaclib { for (size_t i = 0; i < num_reactions; ++i) { const auto&[chapter, qValue, coeffs, reverse, label, rpName, reactants_str, products_str] = records[i]; - // The char arrays from the binary are not guaranteed to be null-terminated // if the string fills the entire buffer. We create null-terminated string_views. const std::string_view label_sv(label, strnlen(label, sizeof(label))); diff --git a/src/lib/reaction/reaction.cpp b/src/lib/reaction/reaction.cpp index fa1d578a..8dad84f8 100644 --- a/src/lib/reaction/reaction.cpp +++ b/src/lib/reaction/reaction.cpp @@ -182,22 +182,29 @@ namespace gridfire::reaction { } - LogicalReaclibReaction::LogicalReaclibReaction(const std::vector& reactants) : - ReaclibReaction( - safe_check_reactant_id(reactants), // Use this first to check if the reactants array is empty and safely exit if so - reactants.front().peName(), - reactants.front().chapter(), - reactants.front().reactants(), - reactants.front().products(), - reactants.front().qValue(), - reactants.front().sourceLabel(), - reactants.front().rateCoefficients(), - reactants.front().is_reverse()) { + LogicalReaclibReaction::LogicalReaclibReaction( + const std::vector& reactions + ) : LogicalReaclibReaction(reactions, false) {} - m_sources.reserve(reactants.size()); - m_rates.reserve(reactants.size()); - for (const auto& reaction : reactants) { - if (std::abs(std::abs(reaction.qValue()) - std::abs(m_qValue)) > 1e-6) { + LogicalReaclibReaction::LogicalReaclibReaction( + const std::vector &reactions, + const bool reverse + ) : + ReaclibReaction( + safe_check_reactant_id(reactions), + reactions.front().peName(), + reactions.front().chapter(), + reactions.front().reactants(), + reactions.front().products(), + reactions.front().qValue(), + reactions.front().sourceLabel(), + reactions.front().rateCoefficients(), + reverse) + { + m_sources.reserve(reactions.size()); + m_rates.reserve(reactions.size()); + for (const auto& reaction : reactions) { + if (std::abs(reaction.qValue() - m_qValue) > 1e-6) { LOG_ERROR( m_logger, "LogicalReaclibReaction constructed with reactions having different Q-values. Expected {} got {}.", @@ -212,6 +219,7 @@ namespace gridfire::reaction { } } + void LogicalReaclibReaction::add_reaction(const ReaclibReaction& reaction) { if (reaction.peName() != m_id) { LOG_ERROR(m_logger, "Cannot add reaction with different peName to LogicalReaclibReaction. Expected {} got {}.", m_id, reaction.peName()); @@ -521,7 +529,8 @@ namespace gridfire::reaction { switch (reaction_ptr->type()) { case ReactionType::REACLIB: { const auto& reaclib_cast_reaction = static_cast(*reaction_ptr); // NOLINT(*-pro-type-static-cast-downcast) - groupedReaclibReactions[std::string(reaclib_cast_reaction.peName())].push_back(reaclib_cast_reaction); + std::string rid = std::format("{}{}", reaclib_cast_reaction.peName(), (reaclib_cast_reaction.is_reverse() ? "_r" : "")); + groupedReaclibReactions[rid].push_back(reaclib_cast_reaction); break; } case ReactionType::LOGICAL_REACLIB: { diff --git a/src/lib/reaction/weak/weak.cpp b/src/lib/reaction/weak/weak.cpp index e47e28c4..0d0713f2 100644 --- a/src/lib/reaction/weak/weak.cpp +++ b/src/lib/reaction/weak/weak.cpp @@ -21,6 +21,13 @@ namespace { {fourdst::atomic::SpeciesErrorType::SPECIES_SYMBOL_NOT_FOUND, "Species symbol not found ((A,Z) out of range)"} }; + std::string normalize_species_id(const fourdst::atomic::Species& species) { + auto result = std::string(species.name()); + std::ranges::transform(result, result.begin(), ::tolower); + std::erase(result, '-'); + return result; + } + fourdst::atomic::Species resolve_weak_product( const gridfire::rates::weak::WeakReactionType type, const fourdst::atomic::Species& reactant @@ -68,18 +75,20 @@ namespace { using namespace gridfire::rates::weak; std::string id; + std::string reactant_id = normalize_species_id(reactant); + std::string product_id = normalize_species_id(product); switch (type) { case WeakReactionType::BETA_MINUS_DECAY: - id = std::format("{}(,ν|)e-,{}", reactant.name(), product.name()); + id = std::format("{}(,ν|)e-,{}", reactant_id, product_id); break; case WeakReactionType::BETA_PLUS_DECAY: - id = std::format("{}(,ν)e+,{}", reactant.name(), product.name()); + id = std::format("{}(,ν)e+,{}", reactant_id, product_id); break; case WeakReactionType::ELECTRON_CAPTURE: - id = std::format("{}(e-,ν){}", reactant.name(), product.name()); + id = std::format("{}(e-,ν){}", reactant_id, product_id); break; case WeakReactionType::POSITRON_CAPTURE: - id = std::format("{}(e+,ν|){}", reactant.name(), product.name()); + id = std::format("{}(e+,ν|){}", reactant_id, product_id); break; } return id; diff --git a/src/lib/solver/strategies/CVODE_solver_strategy.cpp b/src/lib/solver/strategies/CVODE_solver_strategy.cpp index 2ff9f20e..65b6639c 100644 --- a/src/lib/solver/strategies/CVODE_solver_strategy.cpp +++ b/src/lib/solver/strategies/CVODE_solver_strategy.cpp @@ -22,8 +22,6 @@ #include "gridfire/trigger/procedures/trigger_pprint.h" #include "gridfire/utils/general_composition.h" -static size_t s_call_counter = 0; - namespace { std::unordered_map cvode_ret_code_map { {0, "CV_SUCCESS: The solver succeeded."}, @@ -57,7 +55,7 @@ namespace { {-28, "CV_VECTOROP_ERR: A vector operation failed."}, {-29, "CV_PROJ_MEM_NULL: The projection memory structure is NULL."}, {-30, "CV_PROJFUNC_FAIL: The projection function failed in an unrecoverable manner."}, - {-31, "CV_REPTD_PROJFUNC_ERR: THe projection function has repeated recoverable errors."} + {-31, "CV_REPTD_PROJFUNC_ERR: The projection function has repeated recoverable errors."} }; void check_cvode_flag(const int flag, const std::string& func_name) { if (flag < 0) { @@ -74,7 +72,7 @@ namespace { #elif SUNDIALS_HAVE_PTHREADS N_Vector vec = N_VNew_Pthreads(size, sun_ctx); #else - N_Vector vec = N_VNew_Serial(static_cast(size), sun_ctx); + const N_Vector vec = N_VNew_Serial(static_cast(size), sun_ctx); #endif check_cvode_flag(vec == nullptr ? -1 : 0, "N_VNew"); return vec; @@ -129,7 +127,7 @@ namespace gridfire::solver { } CVODESolverStrategy::~CVODESolverStrategy() { - std::cout << "Cleaning up CVODE resources..." << std::endl; + LOG_TRACE_L1(m_logger, "Cleaning up CVODE resources..."); cleanup_cvode_resources(true); if (m_sun_ctx) { @@ -138,10 +136,10 @@ namespace gridfire::solver { } NetOut CVODESolverStrategy::evaluate(const NetIn& netIn) { - LOG_TRACE_L2(m_logger, "Starting solver evaluation with T9: {} and rho: {}", netIn.temperature/1e9, netIn.density); - LOG_TRACE_L2(m_logger, "Building engine update trigger...."); + LOG_TRACE_L1(m_logger, "Starting solver evaluation with T9: {} and rho: {}", netIn.temperature/1e9, netIn.density); + LOG_TRACE_L1(m_logger, "Building engine update trigger...."); auto trigger = trigger::solver::CVODE::makeEnginePartitioningTrigger(1e12, 1e10, 1, true, 10); - LOG_TRACE_L2(m_logger, "Engine update trigger built!"); + LOG_TRACE_L1(m_logger, "Engine update trigger built!"); const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) @@ -149,15 +147,15 @@ namespace gridfire::solver { const auto absTol = m_config.get("gridfire:solver:CVODESolverStrategy:absTol", 1.0e-8); const auto relTol = m_config.get("gridfire:solver:CVODESolverStrategy:relTol", 1.0e-8); - LOG_TRACE_L2(m_logger, "Starting engine update chain..."); + LOG_TRACE_L1(m_logger, "Starting engine update chain..."); fourdst::composition::Composition equilibratedComposition = m_engine.update(netIn); - LOG_TRACE_L2(m_logger, "Engine updated and equilibrated composition found!"); + LOG_TRACE_L1(m_logger, "Engine updated and equilibrated composition found!"); size_t numSpecies = m_engine.getNetworkSpecies().size(); uint64_t N = numSpecies + 1; - LOG_TRACE_L2(m_logger, "Number of species: {}", N); - LOG_TRACE_L2(m_logger, "Initializing CVODE resources"); + LOG_TRACE_L1(m_logger, "Number of species: {}", N); + LOG_TRACE_L1(m_logger, "Initializing CVODE resources"); m_cvode_mem = CVodeCreate(CV_BDF, m_sun_ctx); check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate"); @@ -166,7 +164,7 @@ namespace gridfire::solver { CVODEUserData user_data; user_data.solver_instance = this; user_data.engine = &m_engine; - LOG_TRACE_L2(m_logger, "CVODE resources successfully initialized!"); + LOG_TRACE_L1(m_logger, "CVODE resources successfully initialized!"); double current_time = 0; // ReSharper disable once CppTooWideScope @@ -174,7 +172,7 @@ namespace gridfire::solver { m_num_steps = 0; double accumulated_energy = 0.0; int total_update_stages_triggered = 0; - LOG_TRACE_L2(m_logger, "Starting CVODE iteration"); + LOG_TRACE_L1(m_logger, "Starting CVODE iteration"); while (current_time < netIn.tMax) { user_data.T9 = T9; user_data.rho = netIn.density; @@ -194,6 +192,7 @@ namespace gridfire::solver { std::rethrow_exception(std::make_exception_ptr(*user_data.captured_exception)); } + log_step_diagnostics(user_data); check_cvode_flag(flag, "CVode"); long int n_steps; @@ -359,9 +358,9 @@ namespace gridfire::solver { } int CVODESolverStrategy::cvode_rhs_wrapper( - sunrealtype t, - N_Vector y, - N_Vector ydot, + const sunrealtype t, + const N_Vector y, + const N_Vector ydot, void *user_data ) { auto* data = static_cast(user_data); @@ -454,17 +453,9 @@ namespace gridfire::solver { sunrealtype* ydot_data = N_VGetArrayPointer(ydot); const auto& [dydt, nuclearEnergyGenerationRate] = result.value(); - std::cout << "=========================\n"; for (size_t i = 0; i < numSpecies; ++i) { fourdst::atomic::Species species = m_engine.getNetworkSpecies()[i]; ydot_data[i] = dydt.at(species); - std::cout << "\t" << species << " dydt = " << dydt.at(species) << "\n"; - } - std::cout << "\tε = " << nuclearEnergyGenerationRate << std::endl; - s_call_counter++; - std::cout << "\tMethod called " << s_call_counter << " times" << std::endl; - if (s_call_counter == 31) { - std::cout << "Last reliable step\n"; } ydot_data[numSpecies] = nuclearEnergyGenerationRate; // Set the last element to the specific energy rate } @@ -575,6 +566,8 @@ namespace gridfire::solver { 0.0 ); + std::vector M; + M.reserve(num_components); for (size_t i = 0; i < num_components - 1; i++) { const double weight = relTol * std::abs(y_data[i]) + absTol; if (weight == 0.0) continue; // Skip components with zero weight @@ -583,8 +576,10 @@ namespace gridfire::solver { err_ratios[i] = err_ratio; speciesNames.emplace_back(user_data.networkSpecies->at(i).name()); + M.push_back(user_data.networkSpecies->at(i).mass()); } - fourdst::composition::Composition composition(speciesNames, Y_full); + std::vector X = utils::massFractionFromMolarAbundanceAndMolarMass(Y_full, M); + fourdst::composition::Composition composition(speciesNames, X); if (err_ratios.empty()) { return; @@ -620,9 +615,11 @@ namespace gridfire::solver { columns.push_back(std::make_unique>("Error Ratio", sorted_err_ratios)); std::cout << utils::format_table("Species Error Ratios", columns) << std::endl; + diagnostics::inspect_jacobian_stiffness(*user_data.engine, composition, user_data.T9, user_data.rho); - diagnostics::inspect_species_balance(*user_data.engine, "N-14", composition, user_data.T9, user_data.rho); - diagnostics::inspect_species_balance(*user_data.engine, "n-1", composition, user_data.T9, user_data.rho); + for (const auto& species : sorted_speciesNames) { + diagnostics::inspect_species_balance(*user_data.engine, species, composition, user_data.T9, user_data.rho); + } } }