From c94740f08ffa04d6277dbfa66c54ac65b8eeb246 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Sun, 26 Oct 2025 15:15:03 -0400 Subject: [PATCH] fix(CVODE_solver_strategy): solved a bug wherein mass fractions were not being properly computed from molar abundances. --- .../gridfire/utils/general_composition.h | 56 ++++++++++++++++++- src/lib/engine/views/engine_multiscale.cpp | 8 +-- .../strategies/CVODE_solver_strategy.cpp | 40 +++++++++---- 3 files changed, 89 insertions(+), 15 deletions(-) diff --git a/src/include/gridfire/utils/general_composition.h b/src/include/gridfire/utils/general_composition.h index d4791e9a..e222d645 100644 --- a/src/include/gridfire/utils/general_composition.h +++ b/src/include/gridfire/utils/general_composition.h @@ -3,7 +3,7 @@ #include "fourdst/composition/atomicSpecies.h" namespace gridfire::utils { - inline double massFractionFromMolarAbundance ( + inline double massFractionFromMolarAbundanceAndComposition ( const fourdst::composition::Composition& composition, const fourdst::atomic::Species& species, const double Yi @@ -18,4 +18,58 @@ namespace gridfire::utils { } return (species.mass() * Yi) / sum; }; + + /** + * @brief Convert a vector of molar abundances into a vector of mass fractions + * @param molarAbundances Vector of molar abundances + * @param molarMasses Vector of molar masses + * + * @note The vectors molarAbundances and molarMasses must be parallel. This function does not provide any checks + * to ensure that the correct molar mass is being used with the correct molar abundance. + * @return A vector of molar masses such that each molar mass < 1 and the sum of all is = 1 + */ + inline std::vector massFractionFromMolarAbundanceAndMolarMass ( + const std::vector& molarAbundances, + const std::vector& molarMasses + ) noexcept { + assert(molarMasses.size() == molarAbundances.size()); + assert(!molarMasses.empty()); + + double totalMass = 0; + std::vector masses; + masses.reserve(molarMasses.size()); + for (const auto [m, Y] : std::views::zip(molarMasses, molarAbundances)) { + const double mass = m * Y; + totalMass += mass; + masses.push_back(mass); + } + + assert(totalMass > 0); + + std::vector massFractions; + massFractions.reserve(masses.size()); + std::ranges::transform( + masses, + std::back_inserter(massFractions), + [&totalMass](const double speciesMass) { + const double Xi = speciesMass / totalMass; + if (std::abs(Xi) < 1e-16 && Xi < 0) { + return 0.0; + } + return Xi; + }); + + return massFractions; + } + + inline std::vector molarMassVectorFromComposition( + const fourdst::composition::Composition& composition + ) { + std::vector molarMassVector; + molarMassVector.reserve(composition.getRegisteredSymbols().size()); + for (const auto &entry: composition | std::views::values) { + molarMassVector.push_back(entry.isotope().mass()); + } + return molarMassVector; + } } \ No newline at end of file diff --git a/src/lib/engine/views/engine_multiscale.cpp b/src/lib/engine/views/engine_multiscale.cpp index b410edeb..d322aa5c 100644 --- a/src/lib/engine/views/engine_multiscale.cpp +++ b/src/lib/engine/views/engine_multiscale.cpp @@ -291,7 +291,7 @@ namespace gridfire { for (const auto& species : m_algebraic_species) { const double Yi = m_algebraic_abundances.at(species); - double Xi = utils::massFractionFromMolarAbundance(comp_mutable, species, Yi); + double Xi = utils::massFractionFromMolarAbundanceAndComposition(comp_mutable, species, Yi); comp_mutable.setMassFraction(species, Xi); // Convert Yi (mol/g) to Xi (mass fraction) if (!comp_mutable.finalize(false)) { LOG_ERROR(m_logger, "Failed to finalize composition after setting algebraic species abundance for species '{}'.", species.name()); @@ -1279,7 +1279,7 @@ namespace gridfire { Y_final_qse(i) ); // double Xi = Y_final_qse(i) * species.mass(); // Convert from molar abundance to mass fraction - double Xi = utils::massFractionFromMolarAbundance(normalized_composition, species, Y_final_qse(i)); + double Xi = utils::massFractionFromMolarAbundanceAndComposition(normalized_composition, species, Y_final_qse(i)); if (!outputComposition.hasSpecies(species)) { outputComposition.registerSpecies(species); } @@ -1505,7 +1505,7 @@ namespace gridfire { } auto index = static_cast(m_qse_solve_species_index_map.at(species)); const double molarAbundance = y_qse[index]; - double massFraction = utils::massFractionFromMolarAbundance(m_initial_comp, species, molarAbundance); + double massFraction = utils::massFractionFromMolarAbundanceAndComposition(m_initial_comp, species, molarAbundance); comp_trial.setMassFraction(species, massFraction); } @@ -1572,7 +1572,7 @@ namespace gridfire { comp_trial.registerSpecies(species); } const double molarAbundance = y_qse[static_cast(m_qse_solve_species_index_map.at(species))]; - double massFraction = utils::massFractionFromMolarAbundance(m_initial_comp, species, molarAbundance); + double massFraction = utils::massFractionFromMolarAbundanceAndComposition(m_initial_comp, species, molarAbundance); comp_trial.setMassFraction(species, massFraction); } diff --git a/src/lib/solver/strategies/CVODE_solver_strategy.cpp b/src/lib/solver/strategies/CVODE_solver_strategy.cpp index 666680c2..2ff9f20e 100644 --- a/src/lib/solver/strategies/CVODE_solver_strategy.cpp +++ b/src/lib/solver/strategies/CVODE_solver_strategy.cpp @@ -20,7 +20,9 @@ #include "gridfire/engine/engine_graph.h" #include "gridfire/solver/strategies/triggers/engine_partitioning_trigger.h" #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 { @@ -136,22 +138,26 @@ 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...."); auto trigger = trigger::solver::CVODE::makeEnginePartitioningTrigger(1e12, 1e10, 1, true, 10); + LOG_TRACE_L2(m_logger, "Engine update trigger built!"); + const double T9 = netIn.temperature / 1e9; // Convert temperature from Kelvin to T9 (T9 = T / 1e9) 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..."); fourdst::composition::Composition equilibratedComposition = m_engine.update(netIn); - std::cout << "Equilibrium d molar abundances: " << equilibratedComposition.getMolarAbundance(fourdst::atomic::H_2) << std::endl; - std::cout << "Equilibrium d mass fraction: " << equilibratedComposition.getMassFraction(fourdst::atomic::H_2) << std::endl; - std::cout << "EXITED AT EXPECTED TESTING POINT" << std::endl; - exit(0); + LOG_TRACE_L2(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"); m_cvode_mem = CVodeCreate(CV_BDF, m_sun_ctx); check_cvode_flag(m_cvode_mem == nullptr ? -1 : 0, "CVodeCreate"); @@ -160,6 +166,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!"); double current_time = 0; // ReSharper disable once CppTooWideScope @@ -167,6 +174,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"); while (current_time < netIn.tMax) { user_data.T9 = T9; user_data.rho = netIn.density; @@ -286,6 +294,7 @@ namespace gridfire::solver { } } + LOG_TRACE_L2(m_logger, "CVODE iteration complete"); sunrealtype* y_data = N_VGetArrayPointer(m_Y); accumulated_energy += y_data[numSpecies]; @@ -313,6 +322,7 @@ namespace gridfire::solver { throw std::runtime_error("Failed to finalize output composition after CVODE integration."); } + LOG_TRACE_L2(m_logger, "Constructing output data..."); NetOut netOut; netOut.composition = outputComposition; netOut.energy = accumulated_energy; @@ -326,6 +336,8 @@ namespace gridfire::solver { netOut.dEps_dT = dEps_dT; netOut.dEps_dRho = dEps_dRho; + LOG_TRACE_L2(m_logger, "Output data built!"); + LOG_TRACE_L2(m_logger, "Solver evaluation complete!."); return netOut; } @@ -425,16 +437,15 @@ namespace gridfire::solver { for (const auto& species : m_engine.getNetworkSpecies()) { symbols.emplace_back(species.name()); } - std::vector X; - X.reserve(numSpecies); + std::vector M; + M.reserve(numSpecies); for (size_t i = 0; i < numSpecies; ++i) { const double molarMass = m_engine.getNetworkSpecies()[i].mass(); - X.push_back(y_vec[i] * molarMass); // Convert from molar abundance to mass fraction + M.push_back(molarMass); } + std::vector X = utils::massFractionFromMolarAbundanceAndMolarMass(y_vec, M); fourdst::composition::Composition composition(symbols, X); - std::ranges::replace_if(y_vec, [](const double val) { return val < 0.0; }, 0.0); - const auto result = m_engine.calculateRHSAndEnergy(composition, data->T9, data->rho); if (!result) { throw exceptions::StaleEngineTrigger({data->T9, data->rho, y_vec, t, m_num_steps, y_data[numSpecies]}); @@ -443,8 +454,17 @@ 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) { - ydot_data[i] = dydt.at(m_engine.getNetworkSpecies()[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 }