From 5fba54c1a22c20acd73fb4d5945e28d7f65ae453 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 6 Oct 2025 14:29:33 -0400 Subject: [PATCH] feat(Added-ability-to-get-electron-abundance-and-fixed-some-conversion-bugs): Now Ye can be retrived directly from the composition object. Further a bug which prevented proper conversion to and from number or mass frac modes without messing up the numbers has been resolved --- .../include/fourdst/composition/composition.h | 63 ++- src/composition/lib/composition.cpp | 522 +++++++++++------- tests/composition/compositionTest.cpp | 377 ++++++++++++- 3 files changed, 746 insertions(+), 216 deletions(-) diff --git a/src/composition/include/fourdst/composition/composition.h b/src/composition/include/fourdst/composition/composition.h index b85437d..69af37a 100644 --- a/src/composition/include/fourdst/composition/composition.h +++ b/src/composition/include/fourdst/composition/composition.h @@ -88,6 +88,9 @@ namespace fourdst::composition { double m_massFraction = 0.0; ///< The mass fraction of the species. Valid only if `m_massFracMode` is true. double m_numberFraction = 0.0; ///< The number fraction (mole fraction) of the species. Valid only if `m_massFracMode` is false. double m_relAbundance = 0.0; ///< The relative abundance, used internally for conversions. For mass fraction mode, this is X_i / A_i; for number fraction mode, it's n_i * A_i. + double m_molesPerMass = 0.0; + double m_cachedNumberFraction = 0.0; ///< Cached number fraction for conversions when in mass fraction mode. + bool m_initialized = false; ///< True if the composition entry has been initialized with a valid species. @@ -137,13 +140,6 @@ namespace fourdst::composition { */ [[nodiscard]] double mass_fraction() const; - /** - * @brief Gets the mass fraction, converting from number fraction if necessary. - * @param meanMolarMass The mean molar mass of the entire composition, required for conversion. - * @return The mass fraction of the species. - */ - [[nodiscard]] double mass_fraction(double meanMolarMass) const; - /** * @brief Gets the number fraction of the species. * @pre The entry must be in number fraction mode. @@ -154,10 +150,10 @@ namespace fourdst::composition { /** * @brief Gets the number fraction, converting from mass fraction if necessary. - * @param totalMoles The total moles per unit mass (specific number density) of the entire composition. + * @param totalMolesPerMass The total moles per unit mass (specific number density) of the entire composition. * @return The number fraction of the species. */ - [[nodiscard]] double number_fraction(double totalMoles) const; + [[nodiscard]] double number_fraction(double totalMolesPerMass) const; /** * @brief Gets the relative abundance of the species. @@ -202,10 +198,10 @@ namespace fourdst::composition { /** * @brief Switches the mode to number fraction mode. - * @param totalMoles The total moles per unit mass (specific number density) of the composition. + * @param totalMolesPerMass The total moles per unit mass (specific number density) of the composition. * @return True if the mode was successfully set, false otherwise. */ - bool setNumberFracMode(double totalMoles); + bool setNumberFracMode(double totalMolesPerMass); /** * @brief Overloaded output stream operator for CompositionEntry. @@ -257,8 +253,36 @@ namespace fourdst::composition { */ class Composition { private: - fourdst::config::Config& m_config = fourdst::config::Config::getInstance(); - fourdst::logging::LogManager& m_logManager = fourdst::logging::LogManager::getInstance(); + struct CompositionCache { + std::optional globalComp; ///< Cached global composition data. + std::optional canonicalComp; ///< Cached canonical composition data. + std::optional> massFractions; ///< Cached vector of mass fractions. + std::optional> numberFractions; ///< Cached vector of number fractions. + std::optional> molarAbundances; ///< Cached vector of molar abundances. + std::optional> sortedSpecies; ///< Cached vector of sorted species (by mass). + std::optional> sortedSymbols; ///< Cached vector of sorted species (by mass). + std::optional Ye; ///< Cached electron abundance. + + void clear() { + globalComp = std::nullopt; + canonicalComp = std::nullopt; + massFractions = std::nullopt; + numberFractions = std::nullopt; + molarAbundances = std::nullopt; + sortedSymbols = std::nullopt; + sortedSpecies = std::nullopt; + Ye = std::nullopt; + } + + [[nodiscard]] bool is_clear() const { + return !globalComp.has_value() && !canonicalComp.has_value() && !massFractions.has_value() && + !numberFractions.has_value() && !molarAbundances.has_value() && !sortedSymbols.has_value() && + !Ye.has_value() && !sortedSpecies.has_value(); + } + }; + private: + config::Config& m_config = config::Config::getInstance(); + logging::LogManager& m_logManager = logging::LogManager::getInstance(); quill::Logger* m_logger = m_logManager.getLogger("log"); bool m_finalized = false; ///< True if the composition is finalized. @@ -269,6 +293,9 @@ namespace fourdst::composition { std::set m_registeredSymbols; ///< The registered symbols. std::unordered_map m_compositions; ///< The compositions. + mutable CompositionCache m_cache; ///< Cache for computed properties to avoid redundant calculations. + + /** * @brief Checks if the given symbol is valid by checking against the global species database. * @param symbol The symbol to check. @@ -686,6 +713,14 @@ namespace fourdst::composition { */ [[nodiscard]] double getMeanAtomicNumber() const; + /** + * @brief Compute the electron abundance of the composition. + * @details Ye is defined as the sum over all species of (Z_i * X_i / A_i), where Z_i is the atomic number, X_i is the mass fraction, and A_i is the atomic mass of species i. + * @return Ye (electron abundance). + * @pre The composition must be finalized. + */ + [[nodiscard]] double getElectronAbundance() const; + /** * @brief Creates a new Composition object containing a subset of species from this one. * @param symbols The symbols to include in the subset. @@ -711,7 +746,7 @@ namespace fourdst::composition { * @return True if the isotope is in the composition, false otherwise. * @throws exceptions::CompositionNotFinalizedError if the composition is not finalized. */ - [[nodiscard]] bool contains(const fourdst::atomic::Species& isotope) const; + [[nodiscard]] bool contains(const atomic::Species& isotope) const; /** * @brief Sets the composition mode (mass fraction vs. number fraction). diff --git a/src/composition/lib/composition.cpp b/src/composition/lib/composition.cpp index 4e68031..9574035 100644 --- a/src/composition/lib/composition.cpp +++ b/src/composition/lib/composition.cpp @@ -2,7 +2,7 @@ // // Copyright (C) 2025 -- The 4D-STAR Collaboration // File Author: Emily Boudreaux -// Last Modified: March 26, 2025 +// Last Modified: October 6, 2025 // // 4DSSE is free software; you can use it and/or modify // it under the terms and restrictions the GNU General Library Public @@ -33,11 +33,17 @@ #include "fourdst/composition/atomicSpecies.h" #include "fourdst/composition/species.h" #include "fourdst/composition/composition.h" + +#include + #include "fourdst/composition/exceptions/exceptions_composition.h" namespace { template - std::vector sortVectorBy(std::vector toSort, const std::vector& by) { + std::vector sortVectorBy( + std::vector toSort, + const std::vector& by + ) { std::vector indices(by.size()); for (size_t i = 0; i < indices.size(); i++) { indices[i] = i; @@ -63,20 +69,19 @@ namespace fourdst::composition { CompositionEntry::CompositionEntry() : m_symbol("H-1"), m_isotope(fourdst::atomic::species.at("H-1")), - m_initialized(false) {} // Note: Default entry is uninitialized, must be explicitly set. + m_initialized(false) {} - CompositionEntry::CompositionEntry(const std::string& symbol, const bool massFracMode) : m_symbol(symbol), m_isotope(fourdst::atomic::species.at(symbol)), m_massFracMode(massFracMode) { + CompositionEntry::CompositionEntry( + const std::string& symbol, + const bool massFracMode + ) : + m_symbol(symbol), + m_isotope(atomic::species.at(symbol)), + m_massFracMode(massFracMode) { setSpecies(symbol); } - CompositionEntry::CompositionEntry(const CompositionEntry& entry) : - m_symbol(entry.m_symbol), - m_isotope(entry.m_isotope), - m_massFracMode(entry.m_massFracMode), - m_massFraction(entry.m_massFraction), - m_numberFraction(entry.m_numberFraction), - m_relAbundance(entry.m_relAbundance), - m_initialized(entry.m_initialized) {} + CompositionEntry::CompositionEntry(const CompositionEntry& entry) = default; void CompositionEntry::setSpecies(const std::string& symbol) { if (m_initialized) { @@ -86,7 +91,7 @@ namespace fourdst::composition { throw exceptions::InvalidSpeciesSymbolError("Invalid symbol."); } m_symbol = symbol; - m_isotope = fourdst::atomic::species.at(symbol); + m_isotope = atomic::species.at(symbol); m_initialized = true; } @@ -98,74 +103,80 @@ namespace fourdst::composition { if (!m_massFracMode) { throw exceptions::CompositionModeError("Composition entry is in number fraction mode."); } - return m_massFraction; + // X_i = (moles_i / mass_total) * (mass_i / moles_i) = m_molesPerMass * A_i + return m_molesPerMass * m_isotope.mass(); } - double CompositionEntry::mass_fraction(const double meanMolarMass) const { - if (m_massFracMode) { - return m_massFraction; - } - // Convert from number fraction to mass fraction using: X_i = n_i * A_i / - // where m_relAbundance is n_i * A_i and meanMolarMass is . - return m_relAbundance / meanMolarMass; - } - - double CompositionEntry::number_fraction() const { if (m_massFracMode) { throw exceptions::CompositionModeError("Composition entry is in mass fraction mode."); } - return m_numberFraction; + // In number fraction mode, the value is cached during the mode switch. + return m_cachedNumberFraction; } - double CompositionEntry::number_fraction(const double totalMoles) const { - if (m_massFracMode) { - // Convert from mass fraction to number fraction using: n_i = (X_i / A_i) / sum(X_j / A_j) - // where m_relAbundance is X_i / A_i and totalMoles is sum(X_j / A_j). - return m_relAbundance / totalMoles; - } - return m_numberFraction; + double CompositionEntry::number_fraction( + const double totalMolesPerMass + ) const { + // n_i = (moles_i / mass_total) / (moles_total / mass_total) + if (totalMolesPerMass == 0.0) return 0.0; + return m_molesPerMass / totalMolesPerMass; } double CompositionEntry::rel_abundance() const { - return m_relAbundance; + return m_molesPerMass; } - fourdst::atomic::Species CompositionEntry::isotope() const { + atomic::Species CompositionEntry::isotope() const { return m_isotope; } - void CompositionEntry::setMassFraction(const double mass_fraction) { + void CompositionEntry::setMassFraction( + const double mass_fraction + ) { if (!m_massFracMode) { throw exceptions::CompositionModeError("Composition entry is in number fraction mode."); } - m_massFraction = mass_fraction; - m_relAbundance = m_massFraction / m_isotope.mass(); + // Set the invariant from the given mass fraction + if (m_isotope.mass() == 0.0) { + m_molesPerMass = 0.0; + } else { + m_molesPerMass = mass_fraction / m_isotope.mass(); + } } - void CompositionEntry::setNumberFraction(const double number_fraction) { + void CompositionEntry::setNumberFraction( + const double number_fraction + ) { if (m_massFracMode) { throw exceptions::CompositionModeError("Composition entry is in mass fraction mode."); } - m_numberFraction = number_fraction; - m_relAbundance = m_numberFraction * m_isotope.mass(); + // In number fraction mode, we only cache the value. The invariant + // m_molesPerMass cannot be calculated until finalize() provides global context. + m_cachedNumberFraction = number_fraction; } - bool CompositionEntry::setMassFracMode(const double meanParticleMass) { + bool CompositionEntry::setMassFracMode( + [[maybe_unused]] const double meanMolarMass + ) { if (m_massFracMode) { return false; } m_massFracMode = true; - m_massFraction = m_relAbundance / meanParticleMass; + // The invariant m_molesPerMass does not change when switching mode. + // The cached number fraction is now stale, but that's okay. return true; } - bool CompositionEntry::setNumberFracMode(const double specificNumberDensity) { + bool CompositionEntry::setNumberFracMode( + const double totalMolesPerMass + ) { if (!m_massFracMode) { return false; } m_massFracMode = false; - m_numberFraction = m_relAbundance / specificNumberDensity; + // Calculate and cache the number fraction for the new mode. + m_cachedNumberFraction = number_fraction(totalMolesPerMass); return true; } @@ -173,19 +184,27 @@ namespace fourdst::composition { return m_massFracMode; } - Composition::Composition(const std::vector& symbols) { + Composition::Composition( + const std::vector& symbols + ) { for (const auto& symbol : symbols) { registerSymbol(symbol); } } - Composition::Composition(const std::set& symbols) { + Composition::Composition( + const std::set& symbols + ) { for (const auto& symbol : symbols) { registerSymbol(symbol); } } - Composition::Composition(const std::vector& symbols, const std::vector& fractions, const bool massFracMode) : m_massFracMode(massFracMode) { + Composition::Composition( + const std::vector& symbols, + const std::vector& fractions, + const bool massFracMode + ) : m_massFracMode(massFracMode) { if (symbols.size() != fractions.size()) { LOG_CRITICAL(m_logger, "The number of symbols and fractions must be equal (got {} symbols and {} fractions).", symbols.size(), fractions.size()); throw exceptions::InvalidCompositionError("The number of symbols and fractions must be equal. Got " + std::to_string(symbols.size()) + " symbols and " + std::to_string(fractions.size()) + " fractions."); @@ -194,7 +213,7 @@ namespace fourdst::composition { validateComposition(fractions); for (const auto &symbol : symbols) { - registerSymbol(symbol); + registerSymbol(symbol, m_massFracMode); } for (size_t i = 0; i < symbols.size(); ++i) { @@ -224,25 +243,25 @@ namespace fourdst::composition { m_massFracMode = other.m_massFracMode; m_registeredSymbols = other.m_registeredSymbols; m_compositions = other.m_compositions; - // note: m_config remains bound to the same singleton, so we skip it } return *this; - } - void Composition::registerSymbol(const std::string& symbol, bool massFracMode) { + void Composition::registerSymbol( + const std::string& symbol, + const bool massFracMode + ) { if (!isValidSymbol(symbol)) { LOG_ERROR(m_logger, "Invalid symbol: {}", symbol); throw exceptions::InvalidSymbolError("Invalid symbol: " + symbol); } - // If no symbols have been registered allow mode to be set if (m_registeredSymbols.empty()) { m_massFracMode = massFracMode; } else { if (m_massFracMode != massFracMode) { - LOG_ERROR(m_logger, "Composition is in mass fraction mode. Cannot register symbol ({}) in number fraction mode.", symbol); - throw exceptions::CompositionModeError("Composition is in mass fraction mode. Cannot register symbol (" + symbol + ") in number fraction mode."); + LOG_ERROR(m_logger, "Composition is in {} fraction mode. Cannot register symbol ({}) in {} fraction mode.", m_massFracMode ? "mass" : "number", symbol, massFracMode ? "mass" : "number"); + throw exceptions::CompositionModeError("Composition mode mismatch."); } } @@ -252,22 +271,31 @@ namespace fourdst::composition { } m_registeredSymbols.insert(symbol); - const CompositionEntry entry(symbol, m_massFracMode); - m_compositions[symbol] = entry; + m_compositions[symbol] = CompositionEntry(symbol, m_massFracMode); + m_finalized = false; LOG_TRACE_L3(m_logger, "Registered symbol: {}", symbol); } - void Composition::registerSymbol(const std::vector& symbols, bool massFracMode) { + void Composition::registerSymbol( + const std::vector& symbols, + const bool massFracMode + ) { for (const auto& symbol : symbols) { registerSymbol(symbol, massFracMode); } } - void Composition::registerSpecies(const fourdst::atomic::Species &species, bool massFracMode) { - registerSymbol(std::string(species.name())); + void Composition::registerSpecies( + const atomic::Species &species, + const bool massFracMode + ) { + registerSymbol(std::string(species.name()), massFracMode); } - void Composition::registerSpecies(const std::vector &species, bool massFracMode) { + void Composition::registerSpecies( + const std::vector &species, + const bool massFracMode + ) { for (const auto& s : species) { registerSpecies(s, massFracMode); } @@ -277,14 +305,20 @@ namespace fourdst::composition { return m_registeredSymbols; } - std::set Composition::getRegisteredSpecies() const { - std::set result; + std::set Composition::getRegisteredSpecies() const { + std::set result; for (const auto& entry : m_compositions | std::views::values) { result.insert(entry.isotope()); } return result; } + bool Composition::isValidSymbol( + const std::string& symbol + ) { + return atomic::species.contains(symbol); + } + void Composition::validateComposition(const std::vector& fractions) const { if (!isValidComposition(fractions)) { LOG_ERROR(m_logger, "Invalid composition."); @@ -293,51 +327,37 @@ namespace fourdst::composition { } bool Composition::isValidComposition(const std::vector& fractions) const { - double sum = 0.0; - for (const auto& fraction : fractions) { - sum += fraction; - } + double sum = std::accumulate(fractions.begin(), fractions.end(), 0.0); if (sum < 0.999999 || sum > 1.000001) { LOG_ERROR(m_logger, "The sum of fractions must be equal to 1 (expected 1, got {}).", sum); return false; } - return true; } - bool Composition::isValidSymbol(const std::string& symbol) { - return fourdst::atomic::species.contains(symbol); - } - double Composition::setMassFraction(const std::string& symbol, const double& mass_fraction) { if (!m_registeredSymbols.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol); throw exceptions::UnregisteredSymbolError("Symbol (" + symbol + ") is not registered."); } - if (!m_massFracMode) { LOG_ERROR(m_logger, "Composition is in number fraction mode."); throw exceptions::CompositionModeError("Composition is in number fraction mode."); } - if (mass_fraction < 0.0 || mass_fraction > 1.0) { LOG_ERROR(m_logger, "Mass fraction must be between 0 and 1 for symbol {}. Currently it is {}.", symbol, mass_fraction); - throw exceptions::InvalidCompositionError("Mass fraction must be between 0 and 1 for symbol " + symbol + ". Currently it is " + std::to_string(mass_fraction) + "."); + throw exceptions::InvalidCompositionError("Mass fraction must be between 0 and 1."); } - m_finalized = false; const double old_mass_fraction = m_compositions.at(symbol).mass_fraction(); m_compositions.at(symbol).setMassFraction(mass_fraction); - return old_mass_fraction; } std::vector Composition::setMassFraction(const std::vector& symbols, const std::vector& mass_fractions) { if (symbols.size() != mass_fractions.size()) { - LOG_ERROR(m_logger, "The number of symbols and mass fractions must be equal (currently {} symbols and {} mass fractions).", symbols.size(), mass_fractions.size()); - throw exceptions::InvalidCompositionError("The number of symbols and mass fractions must be equal (currently " + std::to_string(symbols.size()) + " symbols and " + std::to_string(mass_fractions.size()) + " mass fractions)."); + throw exceptions::InvalidCompositionError("The number of symbols and mass fractions must be equal."); } - std::vector old_mass_fractions; old_mass_fractions.reserve(symbols.size()); for (size_t i = 0; i < symbols.size(); ++i) { @@ -346,49 +366,35 @@ namespace fourdst::composition { return old_mass_fractions; } - double Composition::setMassFraction(const fourdst::atomic::Species &species, const double &mass_fraction) { - return setMassFraction(std::string(species.name()), mass_fraction); - } - - std::vector Composition::setMassFraction(const std::vector &species, - const std::vector &mass_fractions) { - std::vector old_mass_fractions; - old_mass_fractions.reserve(species.size()); - for (const auto& spec : species) { - old_mass_fractions.push_back(setMassFraction(spec, mass_fractions[&spec - &species[0]])); - } - return old_mass_fractions; - } - - double Composition::setNumberFraction(const std::string& symbol, const double& number_fraction) { + double Composition::setNumberFraction( + const std::string& symbol, + const double& number_fraction + ) { if (!m_registeredSymbols.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol); throw exceptions::UnregisteredSymbolError("Symbol (" + symbol + ") is not registered."); } - if (m_massFracMode) { - LOG_ERROR(m_logger, "Composition is in mass fraction mode, should be in number fraction mode to call setNumberFraction. Hint: The mode can be switched by first finalizing and then calling setCompositionMode(false)."); - throw exceptions::CompositionModeError("Composition is in mass fraction mode, should be in number fraction mode to call setNumberFraction. Hint: The mode can be switched by first finalizing and then calling setCompositionMode(false)."); + LOG_ERROR(m_logger, "Composition is in mass fraction mode."); + throw exceptions::CompositionModeError("Composition is in mass fraction mode."); } - if (number_fraction < 0.0 || number_fraction > 1.0) { LOG_ERROR(m_logger, "Number fraction must be between 0 and 1 for symbol {}. Currently it is {}.", symbol, number_fraction); - throw exceptions::InvalidCompositionError("Number fraction must be between 0 and 1 for symbol " + symbol + ". Currently it is " + std::to_string(number_fraction) + "."); + throw exceptions::InvalidCompositionError("Number fraction must be between 0 and 1."); } - m_finalized = false; const double old_number_fraction = m_compositions.at(symbol).number_fraction(); m_compositions.at(symbol).setNumberFraction(number_fraction); - return old_number_fraction; } - std::vector Composition::setNumberFraction(const std::vector& symbols, const std::vector& number_fractions) { + std::vector Composition::setNumberFraction( + const std::vector& symbols, + const std::vector& number_fractions + ) { if (symbols.size() != number_fractions.size()) { - LOG_ERROR(m_logger, "The number of symbols and number fractions must be equal. (Currently {} symbols and {} number fractions).", symbols.size(), number_fractions.size()); - throw exceptions::InvalidCompositionError("The number of symbols and number fractions must be equal. (Currently " + std::to_string(symbols.size()) + " symbols and " + std::to_string(number_fractions.size()) + " number fractions)."); + throw exceptions::InvalidCompositionError("The number of symbols and number fractions must be equal."); } - std::vector old_number_fractions; old_number_fractions.reserve(symbols.size()); for (size_t i = 0; i < symbols.size(); ++i) { @@ -397,101 +403,125 @@ namespace fourdst::composition { return old_number_fractions; } - double Composition::setNumberFraction(const fourdst::atomic::Species &species, const double &number_fraction) { + double Composition::setMassFraction( + const atomic::Species &species, + const double &mass_fraction + ) { + return setMassFraction(std::string(species.name()), mass_fraction); + } + + std::vector Composition::setMassFraction( + const std::vector &species, + const std::vector &mass_fractions + ) { + std::vector symbols; + symbols.reserve(species.size()); + for(const auto& s : species) symbols.push_back(std::string(s.name())); + return setMassFraction(symbols, mass_fractions); + } + + double Composition::setNumberFraction( + const atomic::Species &species, + const double &number_fraction + ) { return setNumberFraction(std::string(species.name()), number_fraction); } - std::vector Composition::setNumberFraction(const std::vector &species, - const std::vector &number_fractions) { - std::vector old_number_fractions; - old_number_fractions.reserve(species.size()); - for (const auto& spec : species) { - old_number_fractions.push_back(setNumberFraction(spec, number_fractions[&spec - &species[0]])); - } - return old_number_fractions; + std::vector Composition::setNumberFraction( + const std::vector &species, + const std::vector &number_fractions + ) { + std::vector symbols; + symbols.reserve(species.size()); + for(const auto& s : species) symbols.push_back(std::string(s.name())); + return setNumberFraction(symbols, number_fractions); } bool Composition::finalize(const bool norm) { - bool finalized = false; - if (m_massFracMode) { - finalized = finalizeMassFracMode(norm); - } else { - finalized = finalizeNumberFracMode(norm); - } - if (finalized) { - m_finalized = true; - } - return finalized; + m_specificNumberDensity = 0.0; + m_meanParticleMass = 0.0; + m_finalized = m_massFracMode ? finalizeMassFracMode(norm) : finalizeNumberFracMode(norm); + m_cache.clear(); + return m_finalized; } - bool Composition::finalizeMassFracMode(bool norm) { + bool Composition::finalizeMassFracMode(const bool norm) { std::vector mass_fractions; mass_fractions.reserve(m_compositions.size()); for (const auto &entry: m_compositions | std::views::values) { mass_fractions.push_back(entry.mass_fraction()); } - if (norm) { - double sum = 0.0; - for (const auto& mass_fraction : mass_fractions) { - sum += mass_fraction; - } - for (double & mass_fraction : mass_fractions) { - mass_fraction /= sum; - } + + double sum = std::accumulate(mass_fractions.begin(), mass_fractions.end(), 0.0); + if (norm && sum > 0) { for (auto& [symbol, entry] : m_compositions) { setMassFraction(symbol, entry.mass_fraction() / sum); } + // Recalculate fractions vector after normalization for validation + mass_fractions.clear(); + for (const auto &entry: m_compositions | std::views::values) { + mass_fractions.push_back(entry.mass_fraction()); + } } + try { validateComposition(mass_fractions); } catch ([[maybe_unused]] const exceptions::InvalidCompositionError& e) { - double massSum = 0.0; - for (const auto &entry: m_compositions | std::views::values) { - massSum += entry.mass_fraction(); - } - LOG_ERROR(m_logger, "Composition is invalid (Total mass {}).", massSum); - m_finalized = false; + LOG_ERROR(m_logger, "Composition is invalid after mass frac finalization (Total mass {}).", sum); return false; } - // After validation, calculate the specific number density (total moles per unit mass). + for (const auto &entry: m_compositions | std::views::values) { - m_specificNumberDensity += entry.rel_abundance(); + m_specificNumberDensity += entry.rel_abundance(); // rel_abundance is now consistently moles/mass + } + + if (m_specificNumberDensity > 0) { + m_meanParticleMass = 1.0 / m_specificNumberDensity; } - m_meanParticleMass = 1.0/m_specificNumberDensity; return true; } - bool Composition::finalizeNumberFracMode(bool norm) { + bool Composition::finalizeNumberFracMode(const bool norm) { std::vector number_fractions; number_fractions.reserve(m_compositions.size()); for (const auto &entry: m_compositions | std::views::values) { number_fractions.push_back(entry.number_fraction()); } - if (norm) { - double sum = 0.0; - for (const auto& number_fraction : number_fractions) { - sum += number_fraction; - } + + double sum = std::accumulate(number_fractions.begin(), number_fractions.end(), 0.0); + if (norm && sum > 0) { for (auto& [symbol, entry] : m_compositions) { setNumberFraction(symbol, entry.number_fraction() / sum); } + // Recalculate fractions vector after normalization for validation + number_fractions.clear(); + for (const auto &entry: m_compositions | std::views::values) { + number_fractions.push_back(entry.number_fraction()); + } } + try { validateComposition(number_fractions); - } catch ([[maybe_unused]] const std::runtime_error& e) { - double numberSum = 0.0; - for (const auto &entry: m_compositions | std::views::values) { - numberSum += entry.number_fraction(); - } - LOG_ERROR(m_logger, "Composition is invalid (Total number {}).", numberSum); - m_finalized = false; + } catch ([[maybe_unused]] const exceptions::InvalidCompositionError& e) { + LOG_ERROR(m_logger, "Composition is invalid after number frac finalization (Total number frac {}).", sum); return false; } - // After validation, calculate the mean particle mass. + + // Calculate mean particle mass = sum(n_i * A_i) for (const auto &entry: m_compositions | std::views::values) { - m_meanParticleMass += entry.rel_abundance(); + m_meanParticleMass += entry.number_fraction() * entry.isotope().mass(); + } + + for (auto &entry: m_compositions | std::views::values) { + const double X_i = (m_meanParticleMass > 0) ? (entry.number_fraction() * entry.isotope().mass() / m_meanParticleMass) : 0.0; + entry.m_massFracMode = true; + entry.setMassFraction(X_i); + entry.m_massFracMode = false; + } + + if (m_meanParticleMass > 0) { + m_specificNumberDensity = 1.0 / m_meanParticleMass; } - m_specificNumberDensity = 1.0/m_meanParticleMass; return true; } @@ -547,12 +577,14 @@ namespace fourdst::composition { } if (m_massFracMode) { return m_compositions.at(symbol).mass_fraction(); - } else { - return m_compositions.at(symbol).mass_fraction(m_meanParticleMass); } + + return m_compositions.at(symbol).mass_fraction(); } - double Composition::getMassFraction(const fourdst::atomic::Species &species) const { + double Composition::getMassFraction( + const atomic::Species &species + ) const { return getMassFraction(std::string(species.name())); } @@ -565,7 +597,9 @@ namespace fourdst::composition { } - double Composition::getNumberFraction(const std::string& symbol) const { + double Composition::getNumberFraction( + const std::string& symbol + ) const { if (!m_finalized) { LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize()."); @@ -576,12 +610,13 @@ namespace fourdst::composition { } if (!m_massFracMode) { return m_compositions.at(symbol).number_fraction(); - } else { - return m_compositions.at(symbol).number_fraction(m_specificNumberDensity); } + return m_compositions.at(symbol).number_fraction(m_specificNumberDensity); } - double Composition::getNumberFraction(const fourdst::atomic::Species &species) const { + double Composition::getNumberFraction( + const atomic::Species &species + ) const { return getNumberFraction(std::string(species.name())); } @@ -593,7 +628,9 @@ namespace fourdst::composition { return number_fractions; } - double Composition::getMolarAbundance(const std::string &symbol) const { + double Composition::getMolarAbundance( + const std::string &symbol + ) const { if (!m_finalized) { LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize()."); @@ -606,11 +643,15 @@ namespace fourdst::composition { } - double Composition::getMolarAbundance(const fourdst::atomic::Species &species) const { + double Composition::getMolarAbundance( + const atomic::Species &species + ) const { return getMolarAbundance(std::string(species.name())); } - std::pair Composition::getComposition(const std::string& symbol) const { + std::pair Composition::getComposition( + const std::string& symbol + ) const { if (!m_finalized) { LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize()."); @@ -623,7 +664,8 @@ namespace fourdst::composition { } std::pair Composition::getComposition( - const fourdst::atomic::Species &species) const { + const atomic::Species &species + ) const { return getComposition(std::string(species.name())); } @@ -651,21 +693,39 @@ namespace fourdst::composition { double zSum = 0.0; - // Loop through all registered species in the composition. for (const auto &val: m_compositions | std::views::values) { // Sum of (X_i * Z_i / A_i) zSum += (val.mass_fraction() * val.m_isotope.z())/val.m_isotope.a(); } - // Calculate mean atomic number = * sum(X_i * Z_i / A_i) + // = * sum(X_i * Z_i / A_i) const double mean_A = m_meanParticleMass * zSum; return mean_A; } - Composition Composition::subset(const std::vector& symbols, const std::string& method) const { - const std::array methods = {"norm", "none"}; + double Composition::getElectronAbundance() const { + if (!m_finalized) { + LOG_ERROR(m_logger, "Composition must be finalized before getting the electron abundance. Hint: Consider running .finalize()."); + throw exceptions::CompositionNotFinalizedError("Composition not finalized. Cannot retrieve electron abundance. Hint: Consider running .finalize()."); + } - if (std::ranges::find(methods, method) == methods.end()) { + if (m_cache.Ye.has_value()) { + return m_cache.Ye.value(); + } + + double Ye = 0.0; + for (const auto &val: m_compositions | std::views::values) { + Ye += (val.mass_fraction() * val.m_isotope.z())/val.m_isotope.a(); + } + m_cache.Ye = Ye; + return Ye; + } + + Composition Composition::subset( + const std::vector& symbols, + const std::string& method + ) const { + if (const std::array methods = {"norm", "none"}; std::ranges::find(methods, method) == methods.end()) { const std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'."; LOG_ERROR(m_logger, "Invalid method: {}. Valid methods are norm and none.", method); throw exceptions::InvalidMixingMode(errorMessage); @@ -676,14 +736,12 @@ namespace fourdst::composition { if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw exceptions::UnregisteredSymbolError("Symbol " + symbol + " is not in the composition."); - } else { - subsetComposition.registerSymbol(symbol); } + subsetComposition.registerSymbol(symbol); subsetComposition.setMassFraction(symbol, m_compositions.at(symbol).mass_fraction()); } if (method == "norm") { - const bool isNorm = subsetComposition.finalize(true); - if (!isNorm) { + if (const bool isNorm = subsetComposition.finalize(true); !isNorm) { LOG_ERROR(m_logger, "Subset composition is invalid. (Unable to finalize with normalization)."); throw exceptions::FailedToFinalizeCompositionError("Subset composition is invalid. (Unable to finalize with normalization)."); } @@ -691,13 +749,15 @@ namespace fourdst::composition { return subsetComposition; } - void Composition::setCompositionMode(const bool massFracMode) { + void Composition::setCompositionMode( + const bool massFracMode + ) { if (!m_finalized) { LOG_ERROR(m_logger, "Composition has not been finalized. Mode cannot be set unless composition is finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Mode cannot be set unless composition is finalized. Hint: Consider running .finalize()."); } - bool okay = true; + bool okay; for (auto &entry: m_compositions | std::views::values) { if (massFracMode) { okay = entry.setMassFracMode(m_meanParticleMass); @@ -712,11 +772,16 @@ namespace fourdst::composition { m_massFracMode = massFracMode; } - CanonicalComposition Composition::getCanonicalComposition(bool harsh) const { + CanonicalComposition Composition::getCanonicalComposition( + const bool harsh + ) const { if (!m_finalized) { LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize()."); } + if (m_cache.canonicalComp.has_value()) { + return m_cache.canonicalComp.value(); // Short circuit if we have cached the canonical composition + } CanonicalComposition canonicalComposition; const std::array canonicalH = { "H-1", "H-2", "H-3", "H-4", "H-5", "H-6", "H-7" @@ -737,6 +802,7 @@ namespace fourdst::composition { for (const auto& symbol : getRegisteredSymbols()) { const bool isHSymbol = std::ranges::find(canonicalH, symbol) != std::end(canonicalH); + // ReSharper disable once CppTooWideScopeInitStatement const bool isHeSymbol = std::ranges::find(canonicalHe, symbol) != std::end(canonicalHe); if (isHSymbol || isHeSymbol) { @@ -746,6 +812,7 @@ namespace fourdst::composition { canonicalComposition.Z += getMassFraction(symbol); } + // ReSharper disable once CppTooWideScopeInitStatement const double Z = 1.0 - (canonicalComposition.X + canonicalComposition.Y); if (std::abs(Z - canonicalComposition.Z) > 1e-6) { if (!harsh) { @@ -756,6 +823,7 @@ namespace fourdst::composition { throw std::runtime_error("Validation composition Z (X-Y = " + std::to_string(Z) + ") is different than canonical composition Z (" + std::to_string(canonicalComposition.Z) + ") (∑a_i where a_i != H/He)."); } } + m_cache.canonicalComp = canonicalComposition; return canonicalComposition; } @@ -764,6 +832,9 @@ namespace fourdst::composition { LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize()."); } + if (m_cache.massFractions.has_value()) { + return m_cache.massFractions.value(); // Short circuit if we have cached the mass fractions + } std::vector massFractionVector; std::vector speciesMass; @@ -776,7 +847,9 @@ namespace fourdst::composition { speciesMass.push_back(entry.isotope().mass()); } - return sortVectorBy(massFractionVector, speciesMass); + std::vector massFractions = sortVectorBy(massFractionVector, speciesMass); + m_cache.massFractions = massFractions; // Cache the result + return massFractions; } @@ -785,6 +858,9 @@ namespace fourdst::composition { LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize()."); } + if (m_cache.numberFractions.has_value()) { + return m_cache.numberFractions.value(); // Short circuit if we have cached the number fractions + } std::vector numberFractionVector; std::vector speciesMass; @@ -797,7 +873,9 @@ namespace fourdst::composition { speciesMass.push_back(entry.isotope().mass()); } - return sortVectorBy(numberFractionVector, speciesMass); + std::vector numberFractions = sortVectorBy(numberFractionVector, speciesMass); + m_cache.numberFractions = numberFractions; // Cache the result + return numberFractions; } std::vector Composition::getMolarAbundanceVector() const { @@ -805,6 +883,9 @@ namespace fourdst::composition { LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize()."); } + if (m_cache.molarAbundances.has_value()) { + return m_cache.molarAbundances.value(); // Short circuit if we have cached the molar abundances + } std::vector molarAbundanceVector; std::vector speciesMass; @@ -817,10 +898,15 @@ namespace fourdst::composition { speciesMass.push_back(entry.isotope().mass()); } - return sortVectorBy(molarAbundanceVector, speciesMass); + std::vector molarAbundances = sortVectorBy(molarAbundanceVector, speciesMass); + m_cache.molarAbundances = molarAbundances; // Cache the result + return molarAbundances; + } - size_t Composition::getSpeciesIndex(const std::string &symbol) const { + size_t Composition::getSpeciesIndex( + const std::string &symbol + ) const { if (!m_finalized) { LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize()."); @@ -829,6 +915,16 @@ namespace fourdst::composition { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw exceptions::UnregisteredSymbolError("Symbol " + symbol + " is not in the composition."); } + if (m_cache.sortedSymbols.has_value()) { + return std::distance( + m_cache.sortedSymbols->begin(), + std::ranges::find( + m_cache.sortedSymbols.value().begin(), + m_cache.sortedSymbols.value().end(), + symbol + ) + ); + } std::vector symbols; std::vector speciesMass; @@ -842,10 +938,13 @@ namespace fourdst::composition { } std::vector sortedSymbols = sortVectorBy(symbols, speciesMass); + m_cache.sortedSymbols = sortedSymbols; return std::distance(sortedSymbols.begin(), std::ranges::find(sortedSymbols, symbol)); } - size_t Composition::getSpeciesIndex(const atomic::Species &species) const { + size_t Composition::getSpeciesIndex( + const atomic::Species &species + ) const { if (!m_finalized) { LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize()."); @@ -854,6 +953,16 @@ namespace fourdst::composition { LOG_ERROR(m_logger, "Species {} is not in the composition.", species.name()); throw exceptions::UnregisteredSymbolError("Species " + std::string(species.name()) + " is not in the composition."); } + if (m_cache.sortedSpecies.has_value()) { + return std::distance( + m_cache.sortedSpecies->begin(), + std::ranges::find( + m_cache.sortedSpecies.value().begin(), + m_cache.sortedSpecies.value().end(), + species + ) + ); + } std::vector speciesVector; std::vector speciesMass; @@ -867,10 +976,13 @@ namespace fourdst::composition { } std::vector sortedSpecies = sortVectorBy(speciesVector, speciesMass); + m_cache.sortedSpecies = sortedSpecies; return std::distance(sortedSpecies.begin(), std::ranges::find(sortedSpecies, species)); } - atomic::Species Composition::getSpeciesAtIndex(size_t index) const { + atomic::Species Composition::getSpeciesAtIndex( + size_t index + ) const { if (!m_finalized) { LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize()."); @@ -879,9 +991,8 @@ namespace fourdst::composition { LOG_ERROR(m_logger, "Index {} is out of bounds for composition of size {}.", index, m_compositions.size()); throw std::out_of_range("Index " + std::to_string(index) + " is out of bounds for composition of size " + std::to_string(m_compositions.size()) + "."); } - if (index < 0) { - LOG_ERROR(m_logger, "Index {} is negative. Cannot get species at negative index.", index); - throw std::out_of_range("Index " + std::to_string(index) + " is negative. Cannot get species at negative index."); + if (m_cache.sortedSpecies.has_value()) { + return m_cache.sortedSpecies.value().at(index); } std::vector speciesVector; @@ -899,18 +1010,21 @@ namespace fourdst::composition { return sortedSymbols.at(index); } - bool Composition::hasSymbol(const std::string& symbol) const { + bool Composition::hasSymbol( + const std::string& symbol + ) const { return m_compositions.contains(symbol); } - bool Composition::contains(const fourdst::atomic::Species &isotope) const { + bool Composition::contains( + const atomic::Species &isotope + ) const { // Check if the isotope's symbol is in the composition if (!m_finalized) { LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize()."); throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize()."); } - const auto symbol = static_cast(isotope.name()); - if (m_compositions.contains(symbol)) { + if (const auto symbol = static_cast(isotope.name()); m_compositions.contains(symbol)) { return true; } return false; @@ -918,23 +1032,34 @@ namespace fourdst::composition { /// OVERLOADS - Composition Composition::operator+(const Composition& other) const { + Composition Composition::operator+( + const Composition& other + ) const { return mix(other, 0.5); } - std::ostream& operator<<(std::ostream& os, const GlobalComposition& comp) { + std::ostream& operator<<( + std::ostream& os, + const GlobalComposition& comp + ) { os << "Global Composition: \n"; os << "\tSpecific Number Density: " << comp.specificNumberDensity << "\n"; os << "\tMean Particle Mass: " << comp.meanParticleMass << "\n"; return os; } - std::ostream& operator<<(std::ostream& os, const CompositionEntry& entry) { + std::ostream& operator<<( + std::ostream& os, + const CompositionEntry& entry + ) { os << "<" << entry.m_symbol << " : m_frac = " << entry.mass_fraction() << ">"; return os; } - std::ostream& operator<<(std::ostream& os, const Composition& composition) { + std::ostream& operator<<( + std::ostream& os, + const Composition& composition + ) { os << "Composition(finalized: " << (composition.m_finalized ? "true" : "false") << ", " ; size_t count = 0; for (const auto &entry: composition.m_compositions | std::views::values) { @@ -942,6 +1067,7 @@ namespace fourdst::composition { if (count < composition.m_compositions.size() - 1) { os << ", "; } + count++; } os << ")"; return os; diff --git a/tests/composition/compositionTest.cpp b/tests/composition/compositionTest.cpp index a2b8a5e..71b179c 100644 --- a/tests/composition/compositionTest.cpp +++ b/tests/composition/compositionTest.cpp @@ -116,10 +116,10 @@ TEST_F(compositionTest, registerSymbol) { EXPECT_THROW(comp.registerSymbol("He-21"), fourdst::composition::exceptions::InvalidSymbolError); std::set registeredSymbols = comp.getRegisteredSymbols(); - EXPECT_TRUE(registeredSymbols.find("H-1") != registeredSymbols.end()); - EXPECT_TRUE(registeredSymbols.find("He-4") != registeredSymbols.end()); - EXPECT_TRUE(registeredSymbols.find("H-19") == registeredSymbols.end()); - EXPECT_TRUE(registeredSymbols.find("He-21") == registeredSymbols.end()); + EXPECT_TRUE(registeredSymbols.contains("H-1")); + EXPECT_TRUE(registeredSymbols.contains("He-4")); + EXPECT_TRUE(!registeredSymbols.contains("H-19")); + EXPECT_TRUE(!registeredSymbols.contains("He-21")); } /** @@ -433,3 +433,372 @@ TEST_F(compositionTest, getSpeciesFromAZ) { EXPECT_EQ(fourdst::atomic::O_12, fourdst::atomic::az_to_species(12, 8)); } +TEST_F(compositionTest, constructorWithSymbolsVectorAndSet) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + std::vector vs = {"H-1", "He-4"}; + std::set ss = {"H-1", "He-4"}; + + fourdst::composition::Composition compVec(vs); + EXPECT_TRUE(compVec.hasSymbol("H-1")); + EXPECT_TRUE(compVec.hasSymbol("He-4")); + + fourdst::composition::Composition compSet(ss); + EXPECT_TRUE(compSet.hasSymbol("H-1")); + EXPECT_TRUE(compSet.hasSymbol("He-4")); +} + +TEST_F(compositionTest, constructorWithSymbolsAndFractionsMassAndNumber) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + using fourdst::composition::Composition; + using fourdst::atomic::species; + + // Mass-fraction constructor + std::vector symM = {"H-1", "He-4"}; + std::vector fracM = {0.6, 0.4}; + Composition compM(symM, fracM, true); + EXPECT_NEAR(compM.getMassFraction("H-1"), 0.6, 1e-12); + EXPECT_NEAR(compM.getMassFraction("He-4"), 0.4, 1e-12); + // Mean particle mass and specific number density are reciprocals + double sn = 0.6/species.at("H-1").mass() + 0.4/species.at("He-4").mass(); + double mp = 1.0/sn; + EXPECT_NEAR(compM.getMeanParticleMass(), mp, 1e-12); + + // Number-fraction constructor + std::vector symN = {"H-1", "He-4"}; + std::vector fracN = {0.9, 0.1}; + Composition compN(symN, fracN, false); + EXPECT_NEAR(compN.getNumberFraction("H-1"), 0.9, 1e-12); + EXPECT_NEAR(compN.getNumberFraction("He-4"), 0.1, 1e-12); + double meanA = 0.9*species.at("H-1").mass() + 0.1*species.at("He-4").mass(); + EXPECT_NEAR(compN.getMeanParticleMass(), meanA, 1e-12); + // Check converted mass fractions X_i = n_i * A_i / + double xH = 0.9*species.at("H-1").mass()/meanA; + double xHe = 0.1*species.at("He-4").mass()/meanA; + compN.setCompositionMode(true); + EXPECT_NEAR(compN.getMassFraction("H-1"), xH, 1e-12); + EXPECT_NEAR(compN.getMassFraction("He-4"), xHe, 1e-12); +} + +TEST_F(compositionTest, registerSymbolVectorAndSingleSpeciesAndModeMismatch) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + using fourdst::composition::Composition; + Composition comp; + comp.registerSymbol(std::vector{"H-1", "He-4"}); + EXPECT_TRUE(comp.hasSymbol("H-1")); + EXPECT_TRUE(comp.hasSymbol("He-4")); + + // Register by Species + Composition comp2; + comp2.registerSpecies(fourdst::atomic::H_1); + comp2.registerSpecies(fourdst::atomic::He_4); + EXPECT_TRUE(comp2.hasSymbol("H-1")); + EXPECT_TRUE(comp2.hasSymbol("He-4")); + + // Mode mismatch should throw + Composition comp3; + comp3.registerSymbol("H-1", true); // mass mode + EXPECT_THROW(comp3.registerSymbol("He-4", false), fourdst::composition::exceptions::CompositionModeError); +} + +TEST_F(compositionTest, setMassFractionBySpeciesAndVector) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + using fourdst::composition::Composition; + using fourdst::atomic::H_1; + using fourdst::atomic::He_4; + + Composition comp; + comp.registerSymbol("H-1"); + comp.registerSymbol("He-4"); + + // Single species overload + double old = comp.setMassFraction(H_1, 0.7); + EXPECT_NEAR(old, 0.0, 1e-15); + old = comp.setMassFraction(He_4, 0.3); + EXPECT_NEAR(old, 0.0, 1e-15); + + // Vector overload + std::vector sp = {H_1, He_4}; + std::vector xs = {0.6, 0.4}; + auto olds = comp.setMassFraction(sp, xs); + ASSERT_EQ(olds.size(), 2u); + EXPECT_NEAR(olds[0], 0.7, 1e-12); + EXPECT_NEAR(olds[1], 0.3, 1e-12); + + EXPECT_TRUE(comp.finalize()); + EXPECT_NEAR(comp.getMassFraction("H-1"), 0.6, 1e-12); + EXPECT_NEAR(comp.getMassFraction(He_4), 0.4, 1e-12); +} + +TEST_F(compositionTest, setNumberFractionOverloads) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + using fourdst::composition::Composition; + using fourdst::atomic::H_1; + using fourdst::atomic::He_4; + + Composition comp; + comp.registerSymbol("H-1", false); + comp.registerSymbol("He-4", false); + + // Single symbol + double old = comp.setNumberFraction("H-1", 0.8); + EXPECT_NEAR(old, 0.0, 1e-15); + // Vector of symbols + auto oldv = comp.setNumberFraction(std::vector{"H-1", "He-4"}, std::vector{0.75, 0.25}); + ASSERT_EQ(oldv.size(), 2u); + EXPECT_NEAR(oldv[0], 0.8, 1e-12); + EXPECT_NEAR(oldv[1], 0.0, 1e-12); + + // Species and vector + old = comp.setNumberFraction(H_1, 0.7); + EXPECT_NEAR(old, 0.75, 1e-12); + auto oldsv = comp.setNumberFraction(std::vector{H_1, He_4}, std::vector{0.6, 0.4}); + ASSERT_EQ(oldsv.size(), 2u); + EXPECT_NEAR(oldsv[0], 0.7, 1e-12); + EXPECT_NEAR(oldsv[1], 0.25, 1e-12); + + EXPECT_TRUE(comp.finalize()); + EXPECT_NEAR(comp.getNumberFraction("H-1"), 0.6, 1e-12); + EXPECT_NEAR(comp.getNumberFraction(He_4), 0.4, 1e-12); +} + +TEST_F(compositionTest, mixErrorCases) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + using fourdst::composition::Composition; + + Composition a; a.registerSymbol("H-1"); a.registerSymbol("He-4"); a.setMassFraction("H-1", 0.6); a.setMassFraction("He-4", 0.4); a.finalize(); + Composition b; b.registerSymbol("H-1"); b.registerSymbol("He-4"); b.setMassFraction("H-1", 0.5); b.setMassFraction("He-4", 0.5); + // Not finalized second comp + EXPECT_THROW(static_cast(a.mix(b, 0.5)), fourdst::composition::exceptions::CompositionNotFinalizedError); + b.finalize(); + // Invalid fraction + EXPECT_THROW(static_cast(a.mix(b, -0.1)), fourdst::composition::exceptions::InvalidCompositionError); + EXPECT_THROW(static_cast(a.mix(b, 1.1)), fourdst::composition::exceptions::InvalidCompositionError); +} + +TEST_F(compositionTest, getMassAndNumberFractionMapsAndSpeciesOverloads) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + using fourdst::composition::Composition; + + Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); + comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4); + ASSERT_TRUE(comp.finalize()); + + auto m = comp.getMassFraction(); + ASSERT_EQ(m.size(), 2u); + EXPECT_NEAR(m.at("H-1"), 0.6, 1e-12); + EXPECT_NEAR(m.at("He-4"), 0.4, 1e-12); + EXPECT_NEAR(comp.getMassFraction(fourdst::atomic::H_1), 0.6, 1e-12); + EXPECT_NEAR(comp.getMolarAbundance(fourdst::atomic::H_1), m.at("H-1")/fourdst::atomic::H_1.mass(), 1e-12); + + // Switch to number-fraction mode and verify number maps + comp.setCompositionMode(false); + // Must re-finalize after modifications (mode switch itself keeps values consistent but not finalized status changed? setCompositionMode requires to be finalized; here we just switched modes) + // Set specific number fractions and finalize + comp.setNumberFraction("H-1", 0.7); + comp.setNumberFraction("He-4", 0.3); + ASSERT_TRUE(comp.finalize()); + + auto n = comp.getNumberFraction(); + ASSERT_EQ(n.size(), 2u); + EXPECT_NEAR(n.at("H-1"), 0.7, 1e-12); + EXPECT_NEAR(n.at("He-4"), 0.3, 1e-12); + EXPECT_NEAR(comp.getNumberFraction(fourdst::atomic::He_4), 0.3, 1e-12); +} + +TEST_F(compositionTest, meanAtomicNumberAndElectronAbundance) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + using fourdst::atomic::species; + using fourdst::composition::Composition; + + Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); + comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4); + ASSERT_TRUE(comp.finalize()); + + // Compute expected Ye = sum(X_i * Z_i / A_i) + constexpr double xH = 0.6, xHe = 0.4; + const double aH = species.at("H-1").a(); + const double aHe = species.at("He-4").a(); + const double zH = species.at("H-1").z(); + const double zHe = species.at("He-4").z(); + const double expectedYe = xH*zH/aH + xHe*zHe/aHe; + + EXPECT_NEAR(comp.getElectronAbundance(), expectedYe, 1e-12); + + // = * sum(X_i * Z_i / A_i) + const double expectedZ = comp.getMeanParticleMass() * expectedYe; + EXPECT_NEAR(comp.getMeanAtomicNumber(), expectedZ, 1e-12); +} + +TEST_F(compositionTest, canonicalCompositionAndCaching) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + using fourdst::composition::Composition; + + Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); + comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4); + ASSERT_TRUE(comp.finalize()); + + auto canon1 = comp.getCanonicalComposition(); + EXPECT_NEAR(canon1.X, 0.6, 1e-12); + EXPECT_NEAR(canon1.Y, 0.4, 1e-12); + EXPECT_NEAR(canon1.Z, 0.0, 1e-12); + + // Call again to exercise caching code path + auto canon2 = comp.getCanonicalComposition(); + EXPECT_NEAR(canon2.X, 0.6, 1e-12); + EXPECT_NEAR(canon2.Y, 0.4, 1e-12); + EXPECT_NEAR(canon2.Z, 0.0, 1e-12); + + // Add a metal and re-check + Composition comp2; comp2.registerSymbol("H-1"); comp2.registerSymbol("He-4"); comp2.registerSymbol("O-16"); + comp2.setMassFraction("H-1", 0.6); comp2.setMassFraction("He-4", 0.35); comp2.setMassFraction("O-16", 0.05); + ASSERT_TRUE(comp2.finalize()); + auto canon3 = comp2.getCanonicalComposition(true); + EXPECT_NEAR(canon3.X, 0.6, 1e-12); + EXPECT_NEAR(canon3.Y, 0.35, 1e-12); + EXPECT_NEAR(canon3.Z, 0.05, 1e-12); +} + +TEST_F(compositionTest, vectorsAndIndexingAndSpeciesAtIndex) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + using fourdst::composition::Composition; + using fourdst::atomic::species; + + Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); comp.registerSymbol("O-16"); + comp.setMassFraction("H-1", 0.5); comp.setMassFraction("He-4", 0.3); comp.setMassFraction("O-16", 0.2); + ASSERT_TRUE(comp.finalize()); + + // Mass fraction vector sorted by mass: H-1, He-4, O-16 + auto mv = comp.getMassFractionVector(); + ASSERT_EQ(mv.size(), 3u); + EXPECT_NEAR(mv[0], 0.5, 1e-12); + EXPECT_NEAR(mv[1], 0.3, 1e-12); + EXPECT_NEAR(mv[2], 0.2, 1e-12); + + // Species indices ordering + size_t iH = comp.getSpeciesIndex("H-1"); + size_t iHe = comp.getSpeciesIndex("He-4"); + size_t iO = comp.getSpeciesIndex("O-16"); + EXPECT_LT(iH, iHe); + EXPECT_LT(iHe, iO); + EXPECT_EQ(comp.getSpeciesIndex(fourdst::atomic::H_1), iH); + EXPECT_EQ(comp.getSpeciesIndex(species.at("He-4")), iHe); + + // Species at index + auto sAtHe = comp.getSpeciesAtIndex(iHe); + EXPECT_EQ(std::string(sAtHe.name()), std::string("He-4")); + + // Number fraction vector after switching modes + comp.setCompositionMode(false); + // Tweak number fractions and finalize + // Compute expected number fractions from original mass fractions first + double denom = 0.5/species.at("H-1").mass() + 0.3/species.at("He-4").mass() + 0.2/species.at("O-16").mass(); + double nH_exp = (0.5/species.at("H-1").mass())/denom; + double nHe_exp = (0.3/species.at("He-4").mass())/denom; + double nO_exp = (0.2/species.at("O-16").mass())/denom; + + auto nv0 = comp.getNumberFractionVector(); + ASSERT_EQ(nv0.size(), 3u); + EXPECT_NEAR(nv0[iH], nH_exp, 1e-12); + EXPECT_NEAR(nv0[iHe], nHe_exp, 1e-12); + EXPECT_NEAR(nv0[iO], nO_exp, 1e-12); + + // Molar abundance vector X_i/A_i in mass mode; switch back to mass mode to verify + comp.setCompositionMode(true); + comp.finalize(true); + auto av = comp.getMolarAbundanceVector(); + ASSERT_EQ(av.size(), 3u); + EXPECT_NEAR(av[iH], 0.5/species.at("H-1").mass(), 1e-12); + EXPECT_NEAR(av[iHe], 0.3/species.at("He-4").mass(), 1e-12); + EXPECT_NEAR(av[iO], 0.2/species.at("O-16").mass(), 1e-12); +} + +TEST_F(compositionTest, containsAndPreFinalizationGuards) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + fourdst::composition::Composition comp; + comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); + comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4); + // contains should throw before finalize + EXPECT_THROW(static_cast(comp.contains(fourdst::atomic::H_1)), fourdst::composition::exceptions::CompositionNotFinalizedError); + + ASSERT_TRUE(comp.finalize()); + EXPECT_TRUE(comp.contains(fourdst::atomic::H_1)); + EXPECT_FALSE(comp.contains(fourdst::atomic::Li_6)); +} + +TEST_F(compositionTest, subsetNoneMethodAndNormalizationFlow) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + fourdst::composition::Composition comp; + comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); comp.registerSymbol("O-16"); + comp.setMassFraction("H-1", 0.5); comp.setMassFraction("He-4", 0.3); comp.setMassFraction("O-16", 0.2); + ASSERT_TRUE(comp.finalize()); + + fourdst::composition::Composition sub = comp.subset(std::vector{"H-1", "He-4"}, "none"); + // Not normalized: finalize without normalization should fail + EXPECT_FALSE(sub.finalize(false)); + // With normalization, it should succeed and scale to sum to 1 + EXPECT_TRUE(sub.finalize(true)); + double sum = sub.getMassFraction("H-1") + sub.getMassFraction("He-4"); + EXPECT_NEAR(sum, 1.0, 1e-12); + EXPECT_NEAR(sub.getMassFraction("H-1"), 0.5/(0.5+0.3), 1e-12); + EXPECT_NEAR(sub.getMassFraction("He-4"), 0.3/(0.5+0.3), 1e-12); +} + +TEST_F(compositionTest, copyConstructorAndAssignmentIndependence) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + using fourdst::composition::Composition; + + Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); + comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4); + ASSERT_TRUE(comp.finalize()); + + Composition copy(comp); // copy ctor + EXPECT_NEAR(copy.getMassFraction("H-1"), 0.6, 1e-12); + EXPECT_NEAR(copy.getMassFraction("He-4"), 0.4, 1e-12); + + Composition assigned; assigned = comp; // assignment + EXPECT_NEAR(assigned.getMassFraction("H-1"), 0.6, 1e-12); + EXPECT_NEAR(assigned.getMassFraction("He-4"), 0.4, 1e-12); + + // Modify original and ensure copies do not change + comp.setMassFraction("H-1", 0.7); comp.setMassFraction("He-4", 0.3); ASSERT_TRUE(comp.finalize()); + EXPECT_NEAR(copy.getMassFraction("H-1"), 0.6, 1e-12); + EXPECT_NEAR(assigned.getMassFraction("He-4"), 0.4, 1e-12); +} + +TEST_F(compositionTest, getCompositionBySpeciesAndAllEntries) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + using fourdst::composition::Composition; + + Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); + comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4); + ASSERT_TRUE(comp.finalize()); + + auto pairBySpec = comp.getComposition(fourdst::atomic::H_1); + EXPECT_NEAR(pairBySpec.first.mass_fraction(), 0.6, 1e-12); + EXPECT_NEAR(pairBySpec.second.meanParticleMass, comp.getMeanParticleMass(), 1e-15); + + auto all = comp.getComposition(); + ASSERT_EQ(all.first.size(), 2u); + EXPECT_NEAR(all.first.at("H-1").mass_fraction(), 0.6, 1e-12); + EXPECT_NEAR(all.first.at("He-4").mass_fraction(), 0.4, 1e-12); + EXPECT_NEAR(all.second.meanParticleMass, comp.getMeanParticleMass(), 1e-15); +} + +TEST_F(compositionTest, iterationBeginEndAndIndexOutOfRange) { + fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME); + fourdst::composition::Composition comp; + comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); comp.registerSymbol("O-16"); + comp.setMassFraction("H-1", 0.5); comp.setMassFraction("He-4", 0.3); comp.setMassFraction("O-16", 0.2); + ASSERT_TRUE(comp.finalize()); + + // Iterate and count entries + size_t count = 0; + for (auto it = comp.begin(); it != comp.end(); ++it) { + count++; + } + EXPECT_EQ(count, 3u); + + // Out-of-range access + EXPECT_THROW(static_cast(comp.getSpeciesAtIndex(100)), std::out_of_range); +} +