From 517a4cf44da059b443ca4a0ce78c2dd8b6dd1b8d Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 16 Jun 2025 13:15:34 -0400 Subject: [PATCH 1/4] docs(composition): @example -> *Example Usage:* --- src/composition/public/composition.h | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/composition/public/composition.h b/src/composition/public/composition.h index 7226393..e61a5c6 100644 --- a/src/composition/public/composition.h +++ b/src/composition/public/composition.h @@ -68,7 +68,7 @@ namespace serif::composition { * @brief Constructs a CompositionEntry with the given symbol and mode. * @param symbol The chemical symbol of the species. * @param massFracMode True if mass fraction mode, false if number fraction mode. - * @example + * *Example Usage:* * @code * CompositionEntry entry("H", true); * @endcode @@ -185,13 +185,13 @@ namespace serif::composition { * - The only exception to the finalize rule is if the compositon was constructed with symbols and mass fractions at instantiation time. In this case, the composition is automatically finalized. * however, this means that the composition passed to the constructor must be valid. * - * @example Constructing a finalized composition with symbols and mass fractions: + * *Example Usage:* Constructing a finalized composition with symbols and mass fractions: * @code * std::vector symbols = {"H", "He"}; * std::vector mass_fractions = {0.7, 0.3}; * Composition comp(symbols, mass_fractions); * @endcode - * @example Constructing a composition with symbols and finalizing it later: + * *Example Usage:* Constructing a composition with symbols and finalizing it later: * @code * std::vector symbols = {"H", "He"}; * Composition comp(symbols); @@ -271,7 +271,7 @@ namespace serif::composition { /** * @brief Constructs a Composition with the given symbols. * @param symbols The symbols to initialize the composition with. - * @example + * *Example Usage:* * @code * std::vector symbols = {"H", "O"}; * Composition comp(symbols); @@ -282,7 +282,7 @@ namespace serif::composition { /** * @brief Constructs a Composition with the given symbols as a set. * @param symbols The symbols to initialize the composition with. - * @example + * *Example Usage:* * @code * std::set symbols = {"H", "O"}; * Composition comp(symbols); @@ -295,7 +295,7 @@ namespace serif::composition { * @param symbols The symbols to initialize the composition with. * @param mass_fractions The mass fractions corresponding to the symbols. * @param massFracMode True if mass fraction mode, false if number fraction mode. - * @example + * *Example Usage:* * @code * std::vector symbols = {"H", "O"}; * std::vector mass_fractions = {0.1, 0.9}; @@ -308,7 +308,7 @@ namespace serif::composition { * @brief Registers a new symbol. * @param symbol The symbol to register. * @param massFracMode True if mass fraction mode, false if number fraction mode. - * @example + * *Example Usage:* * @code * Composition comp; * comp.registerSymbol("H"); @@ -320,7 +320,7 @@ namespace serif::composition { * @brief Registers multiple new symbols. * @param symbols The symbols to register. * @param massFracMode True if mass fraction mode, false if number fraction mode. - * @example + * *Example Usage:* * @code * std::vector symbols = {"H", "O"}; * Composition comp; @@ -340,7 +340,7 @@ namespace serif::composition { * @param symbol The symbol to set the mass fraction for. * @param mass_fraction The mass fraction to set. * @return The mass fraction that was set. - * @example + * *Example Usage:* * @code * Composition comp; * comp.setMassFraction("H", 0.1); @@ -353,7 +353,7 @@ namespace serif::composition { * @param symbols The symbols to set the mass fraction for. * @param mass_fractions The mass fractions corresponding to the symbols. * @return A vector of mass fractions that were set. - * @example + * *Example Usage:* * @code * std::vector symbols = {"H", "O"}; * std::vector mass_fractions = {0.1, 0.9}; From 4f47dc97a8198aba049e5890a43058777afd0179 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 16 Jun 2025 15:00:33 -0400 Subject: [PATCH 2/4] feat(eos): EOS now uses composition module EOS code updated to make use of composition module. This is part of the larger change to change all composition handling to use the composition module. Note that the current implimentation is a bit hacky since it simply copies back and forth to the alredy used HELMEOSInput and HELMEOSOutput structs. In furture this can be more tightly connected to avoid extra copies. --- src/composition/private/composition.cpp | 1018 ++++++++++++----------- src/composition/public/composition.h | 38 +- 2 files changed, 569 insertions(+), 487 deletions(-) diff --git a/src/composition/private/composition.cpp b/src/composition/private/composition.cpp index 6af4a88..9543836 100644 --- a/src/composition/private/composition.cpp +++ b/src/composition/private/composition.cpp @@ -25,6 +25,7 @@ #include #include #include +#include #include @@ -32,539 +33,600 @@ namespace serif::composition { -CompositionEntry::CompositionEntry() : m_symbol("H-1"), m_isotope(chemSpecies::species.at("H-1")), m_initialized(false) {} - -CompositionEntry::CompositionEntry(const std::string& symbol, bool massFracMode) : m_symbol(symbol), m_isotope(chemSpecies::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) {} - -void CompositionEntry::setSpecies(const std::string& symbol) { - if (m_initialized) { - throw std::runtime_error("Composition entry is already initialized."); + CompositionEntry::CompositionEntry() : m_symbol("H-1"), m_isotope(chemSpecies::species.at("H-1")), + m_initialized(false) { } - if (chemSpecies::species.count(symbol) == 0) { - throw std::runtime_error("Invalid symbol."); + + CompositionEntry::CompositionEntry(const std::string& symbol, const bool massFracMode) : m_symbol(symbol), m_isotope(chemSpecies::species.at(symbol)), m_massFracMode(massFracMode) { + setSpecies(symbol); } - m_symbol = symbol; - m_isotope = chemSpecies::species.at(symbol); - m_initialized = true; -} -std::string CompositionEntry::symbol() const { - return m_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) {} -double CompositionEntry::mass_fraction() const { - if (!m_massFracMode) { - throw std::runtime_error("Composition entry is in number fraction mode."); + void CompositionEntry::setSpecies(const std::string& symbol) { + if (m_initialized) { + throw std::runtime_error("Composition entry is already initialized."); + } + if (chemSpecies::species.count(symbol) == 0) { + throw std::runtime_error("Invalid symbol."); + } + m_symbol = symbol; + m_isotope = chemSpecies::species.at(symbol); + m_initialized = true; } - return m_massFraction; -} -double CompositionEntry::mass_fraction(double meanMolarMass) const { - if (m_massFracMode) { + std::string CompositionEntry::symbol() const { + return m_symbol; + } + + double CompositionEntry::mass_fraction() const { + if (!m_massFracMode) { + throw std::runtime_error("Composition entry is in number fraction mode."); + } return m_massFraction; } - return m_relAbundance / meanMolarMass; -} - -double CompositionEntry::number_fraction() const { - if (m_massFracMode) { - throw std::runtime_error("Composition entry is in mass fraction mode."); - } - return m_numberFraction; -} - -double CompositionEntry::number_fraction(double totalMoles) const { - if (m_massFracMode) { - return m_relAbundance / totalMoles; - } - return m_numberFraction; -} - -double CompositionEntry::rel_abundance() const { - return m_relAbundance; -} - -chemSpecies::Species CompositionEntry::isotope() const { - return m_isotope; -} - -void CompositionEntry::setMassFraction(double mass_fraction) { - if (!m_massFracMode) { - throw std::runtime_error("Composition entry is in number fraction mode."); - } - m_massFraction = mass_fraction; - m_relAbundance = m_massFraction / m_isotope.mass(); -} - -void CompositionEntry::setNumberFraction(double number_fraction) { - if (m_massFracMode) { - throw std::runtime_error("Composition entry is in mass fraction mode."); - } - m_numberFraction = number_fraction; - m_relAbundance = m_numberFraction * m_isotope.mass(); -} - -bool CompositionEntry::setMassFracMode(double meanParticleMass) { - if (m_massFracMode) { - return false; - } - m_massFracMode = true; - m_massFraction = m_relAbundance / meanParticleMass; - return true; -} - -bool CompositionEntry::setNumberFracMode(double specificNumberDensity) { - if (!m_massFracMode) { - return false; - } - m_massFracMode = false; - m_numberFraction = m_relAbundance / specificNumberDensity; - return true; -} - -bool CompositionEntry::getMassFracMode() const { - return m_massFracMode; -} - -Composition::Composition(const std::vector& symbols) { - for (const auto& symbol : symbols) { - registerSymbol(symbol); - } -} - -Composition::Composition(const std::set& symbols) { - for (const auto& symbol : symbols) { - registerSymbol(symbol); - } -} - -Composition::Composition(const std::vector& symbols, const std::vector& fractions, bool massFracMode) : m_massFracMode(massFracMode) { - if (symbols.size() != fractions.size()) { - LOG_ERROR(m_logger, "The number of symbols and fractions must be equal."); - throw std::runtime_error("The number of symbols and fractions must be equal."); - } - - validateComposition(fractions); - - for (const auto &symbol : symbols) { - registerSymbol(symbol); - } - - for (size_t i = 0; i < symbols.size(); ++i) { + double CompositionEntry::mass_fraction(double meanMolarMass) const { if (m_massFracMode) { - setMassFraction(symbols[i], fractions[i]); + return m_massFraction; + } + return m_relAbundance / meanMolarMass; + } + + + double CompositionEntry::number_fraction() const { + if (m_massFracMode) { + throw std::runtime_error("Composition entry is in mass fraction mode."); + } + return m_numberFraction; + } + + double CompositionEntry::number_fraction(double totalMoles) const { + if (m_massFracMode) { + return m_relAbundance / totalMoles; + } + return m_numberFraction; + } + + double CompositionEntry::rel_abundance() const { + return m_relAbundance; + } + + chemSpecies::Species CompositionEntry::isotope() const { + return m_isotope; + } + + void CompositionEntry::setMassFraction(double mass_fraction) { + if (!m_massFracMode) { + throw std::runtime_error("Composition entry is in number fraction mode."); + } + m_massFraction = mass_fraction; + m_relAbundance = m_massFraction / m_isotope.mass(); + } + + void CompositionEntry::setNumberFraction(double number_fraction) { + if (m_massFracMode) { + throw std::runtime_error("Composition entry is in mass fraction mode."); + } + m_numberFraction = number_fraction; + m_relAbundance = m_numberFraction * m_isotope.mass(); + } + + bool CompositionEntry::setMassFracMode(double meanParticleMass) { + if (m_massFracMode) { + return false; + } + m_massFracMode = true; + m_massFraction = m_relAbundance / meanParticleMass; + return true; + } + + bool CompositionEntry::setNumberFracMode(double specificNumberDensity) { + if (!m_massFracMode) { + return false; + } + m_massFracMode = false; + m_numberFraction = m_relAbundance / specificNumberDensity; + return true; + } + + bool CompositionEntry::getMassFracMode() const { + return m_massFracMode; + } + + Composition::Composition(const std::vector& symbols) { + for (const auto& symbol : symbols) { + registerSymbol(symbol); + } + } + + Composition::Composition(const std::set& symbols) { + for (const auto& symbol : symbols) { + registerSymbol(symbol); + } + } + + Composition::Composition(const std::vector& symbols, const std::vector& fractions, bool massFracMode) : m_massFracMode(massFracMode) { + if (symbols.size() != fractions.size()) { + LOG_ERROR(m_logger, "The number of symbols and fractions must be equal."); + throw std::runtime_error("The number of symbols and fractions must be equal."); + } + + validateComposition(fractions); + + for (const auto &symbol : symbols) { + registerSymbol(symbol); + } + + for (size_t i = 0; i < symbols.size(); ++i) { + if (m_massFracMode) { + setMassFraction(symbols[i], fractions[i]); + } else { + setNumberFraction(symbols[i], fractions[i]); + } + } + finalize(); + } + + Composition::Composition(const Composition &composition) { + m_finalized = composition.m_finalized; + m_specificNumberDensity = composition.m_specificNumberDensity; + m_meanParticleMass = composition.m_meanParticleMass; + m_massFracMode = composition.m_massFracMode; + m_registeredSymbols = composition.m_registeredSymbols; + m_compositions = composition.m_compositions; + } + + Composition& Composition::operator=(const Composition &other) { + if (this != &other) { + m_finalized = other.m_finalized; + m_specificNumberDensity = other.m_specificNumberDensity; + m_meanParticleMass = other.m_meanParticleMass; + 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) { + if (!isValidSymbol(symbol)) { + LOG_ERROR(m_logger, "Invalid symbol: {}", symbol); + throw std::runtime_error("Invalid symbol."); + } + + // If no symbols have been registered allow mode to be set + if (m_registeredSymbols.size() == 0) { + m_massFracMode = massFracMode; } else { - setNumberFraction(symbols[i], fractions[i]); + if (m_massFracMode != massFracMode) { + LOG_ERROR(m_logger, "Composition is in mass fraction mode. Cannot register symbol in number fraction mode."); + throw std::runtime_error("Composition is in mass fraction mode. Cannot register symbol in number fraction mode."); + } } - } - finalize(); -} -void Composition::registerSymbol(const std::string& symbol, bool massFracMode) { - if (!isValidSymbol(symbol)) { - LOG_ERROR(m_logger, "Invalid symbol: {}", symbol); - throw std::runtime_error("Invalid symbol."); + if (m_registeredSymbols.find(symbol) != m_registeredSymbols.end()) { + LOG_WARNING(m_logger, "Symbol {} is already registered.", symbol); + return; + } + + m_registeredSymbols.insert(symbol); + CompositionEntry entry(symbol, m_massFracMode); + m_compositions[symbol] = entry; + LOG_INFO(m_logger, "Registered symbol: {}", symbol); } - // If no symbols have been registered allow mode to be set - if (m_registeredSymbols.size() == 0) { - 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."); - throw std::runtime_error("Composition is in mass fraction mode. Cannot register symbol in number fraction mode."); + void Composition::registerSymbol(const std::vector& symbols, bool massFracMode) { + for (const auto& symbol : symbols) { + registerSymbol(symbol, massFracMode); } } - if (m_registeredSymbols.find(symbol) != m_registeredSymbols.end()) { - LOG_WARNING(m_logger, "Symbol {} is already registered.", symbol); - return; + std::set Composition::getRegisteredSymbols() const { + return m_registeredSymbols; } - m_registeredSymbols.insert(symbol); - CompositionEntry entry(symbol, m_massFracMode); - m_compositions[symbol] = entry; - LOG_INFO(m_logger, "Registered symbol: {}", symbol); -} - -void Composition::registerSymbol(const std::vector& symbols, bool massFracMode) { - for (const auto& symbol : symbols) { - registerSymbol(symbol, massFracMode); - } -} - -std::set Composition::getRegisteredSymbols() const { - return m_registeredSymbols; -} - -void Composition::validateComposition(const std::vector& fractions) const { - if (!isValidComposition(fractions)) { - LOG_ERROR(m_logger, "Invalid composition."); - throw std::runtime_error("Invalid composition."); - } -} - -bool Composition::isValidComposition(const std::vector& fractions) const { - double sum = 0.0; - for (const auto& fraction : fractions) { - sum += fraction; - } - if (sum < 0.999999 || sum > 1.000001) { - LOG_ERROR(m_logger, "The sum of fractions must be equal to 1."); - return false; + void Composition::validateComposition(const std::vector& fractions) const { + if (!isValidComposition(fractions)) { + LOG_ERROR(m_logger, "Invalid composition."); + throw std::runtime_error("Invalid composition."); + } } - return true; -} - -bool Composition::isValidSymbol(const std::string& symbol) const { - return chemSpecies::species.count(symbol) > 0; -} - -double Composition::setMassFraction(const std::string& symbol, const double& mass_fraction) { - if (m_registeredSymbols.find(symbol) == m_registeredSymbols.end()) { - LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol); - throw std::runtime_error("Symbol is not registered."); - } - - if (!m_massFracMode) { - LOG_ERROR(m_logger, "Composition is in number fraction mode."); - throw std::runtime_error("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 std::runtime_error("Mass fraction must be between 0 and 1."); - } - - m_finalized = false; - 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."); - throw std::runtime_error("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) { - old_mass_fractions.push_back(setMassFraction(symbols[i], mass_fractions[i])); - } - return old_mass_fractions; -} - -double Composition::setNumberFraction(const std::string& symbol, const double& number_fraction) { - if (m_registeredSymbols.find(symbol) == m_registeredSymbols.end()) { - LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol); - throw std::runtime_error("Symbol is not registered."); - } - - if (m_massFracMode) { - LOG_ERROR(m_logger, "Composition is in mass fraction mode."); - throw std::runtime_error("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 std::runtime_error("Number fraction must be between 0 and 1."); - } - - m_finalized = false; - 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) { - if (symbols.size() != number_fractions.size()) { - LOG_ERROR(m_logger, "The number of symbols and number fractions must be equal."); - throw std::runtime_error("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) { - old_number_fractions.push_back(setNumberFraction(symbols[i], number_fractions[i])); - } - return old_number_fractions; -} - -bool Composition::finalize(bool norm) { - bool finalized = false; - if (m_massFracMode) { - finalized = finalizeMassFracMode(norm); - } else { - finalized = finalizeNumberFracMode(norm); - } - if (finalized) { - m_finalized = true; - } - return finalized; -} - -bool Composition::finalizeMassFracMode(bool norm) { - std::vector mass_fractions; - mass_fractions.reserve(m_compositions.size()); - for (const auto& [_, entry] : m_compositions) { - mass_fractions.push_back(entry.mass_fraction()); - } - if (norm) { + bool Composition::isValidComposition(const std::vector& fractions) const { double sum = 0.0; - for (const auto& mass_fraction : mass_fractions) { - sum += mass_fraction; + for (const auto& fraction : fractions) { + sum += fraction; } - for (int i = 0; i < mass_fractions.size(); ++i) { - mass_fractions[i] /= sum; - } - for (auto& [symbol, entry] : m_compositions) { - setMassFraction(symbol, entry.mass_fraction() / sum); + if (sum < 0.999999 || sum > 1.000001) { + LOG_ERROR(m_logger, "The sum of fractions must be equal to 1."); + return false; } + + return true; } - try { - validateComposition(mass_fractions); - } catch (const std::runtime_error& e) { - double massSum = 0.0; - for (const auto& [_, entry] : m_compositions) { - massSum += entry.mass_fraction(); + + bool Composition::isValidSymbol(const std::string& symbol) const { + return chemSpecies::species.count(symbol) > 0; + } + + double Composition::setMassFraction(const std::string& symbol, const double& mass_fraction) { + if (m_registeredSymbols.find(symbol) == m_registeredSymbols.end()) { + LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol); + throw std::runtime_error("Symbol is not registered."); } - LOG_ERROR(m_logger, "Composition is invalid (Total mass {}).", massSum); + + if (!m_massFracMode) { + LOG_ERROR(m_logger, "Composition is in number fraction mode."); + throw std::runtime_error("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 std::runtime_error("Mass fraction must be between 0 and 1."); + } + m_finalized = false; - return false; - } - for (const auto& [_, entry] : m_compositions) { - m_specificNumberDensity += entry.rel_abundance(); - } - m_meanParticleMass = 1.0/m_specificNumberDensity; - return true; -} + double old_mass_fraction = m_compositions.at(symbol).mass_fraction(); + m_compositions.at(symbol).setMassFraction(mass_fraction); -bool Composition::finalizeNumberFracMode(bool norm) { - std::vector number_fractions; - number_fractions.reserve(m_compositions.size()); - for (const auto& [_, entry] : m_compositions) { - number_fractions.push_back(entry.number_fraction()); + return old_mass_fraction; } - if (norm) { - double sum = 0.0; - for (const auto& number_fraction : number_fractions) { - sum += number_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."); + throw std::runtime_error("The number of symbols and mass fractions must be equal."); } - for (auto& [symbol, entry] : m_compositions) { - setNumberFraction(symbol, entry.number_fraction() / sum); + + std::vector old_mass_fractions; + old_mass_fractions.reserve(symbols.size()); + for (size_t i = 0; i < symbols.size(); ++i) { + old_mass_fractions.push_back(setMassFraction(symbols[i], mass_fractions[i])); } + return old_mass_fractions; } - try { - validateComposition(number_fractions); - } catch (const std::runtime_error& e) { - double numberSum = 0.0; - for (const auto& [_, entry] : m_compositions) { - numberSum += entry.number_fraction(); + + double Composition::setNumberFraction(const std::string& symbol, const double& number_fraction) { + if (m_registeredSymbols.find(symbol) == m_registeredSymbols.end()) { + LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol); + throw std::runtime_error("Symbol is not registered."); } - LOG_ERROR(m_logger, "Composition is invalid (Total number {}).", numberSum); + + if (m_massFracMode) { + LOG_ERROR(m_logger, "Composition is in mass fraction mode."); + throw std::runtime_error("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 std::runtime_error("Number fraction must be between 0 and 1."); + } + m_finalized = false; - return false; - } - for (const auto& [_, entry] : m_compositions) { - m_meanParticleMass += entry.rel_abundance(); - } - m_specificNumberDensity = 1.0/m_meanParticleMass; - return true; -} + double old_number_fraction = m_compositions.at(symbol).number_fraction(); + m_compositions.at(symbol).setNumberFraction(number_fraction); -Composition Composition::mix(const Composition& other, double fraction) const { - if (!m_finalized || !other.m_finalized) { - LOG_ERROR(m_logger, "Compositions have not both been finalized."); - throw std::runtime_error("Compositions have not been finalized (Consider running .finalize())."); + return old_number_fraction; } - if (fraction < 0.0 || fraction > 1.0) { - LOG_ERROR(m_logger, "Fraction must be between 0 and 1."); - throw std::runtime_error("Fraction must be between 0 and 1."); + 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."); + throw std::runtime_error("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) { + old_number_fractions.push_back(setNumberFraction(symbols[i], number_fractions[i])); + } + return old_number_fractions; } - std::set mixedSymbols = other.getRegisteredSymbols(); - // Get the union of the two sets - mixedSymbols.insert(m_registeredSymbols.begin(), m_registeredSymbols.end()); - - Composition mixedComposition(mixedSymbols); - for (const auto& symbol : mixedSymbols) { - double thisMassFrac, otherMassFrac = 0.0; - - thisMassFrac = hasSymbol(symbol) ? getMassFraction(symbol) : 0.0; - otherMassFrac = other.hasSymbol(symbol) ? other.getMassFraction(symbol) : 0.0; - - double massFraction = fraction * thisMassFrac + otherMassFrac * (1-fraction); - mixedComposition.setMassFraction(symbol, massFraction); - } - mixedComposition.finalize(); - return mixedComposition; -} - -double Composition::getMassFraction(const std::string& symbol) const { - if (!m_finalized) { - LOG_ERROR(m_logger, "Composition has not been finalized."); - throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); - } - if (m_compositions.count(symbol) == 0) { - LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); - throw std::runtime_error("Symbol is not in the composition."); - } - if (m_massFracMode) { - return m_compositions.at(symbol).mass_fraction(); - } else { - return m_compositions.at(symbol).mass_fraction(m_meanParticleMass); - } -} - -std::unordered_map Composition::getMassFraction() const { - std::unordered_map mass_fractions; - for (const auto& [symbol, entry] : m_compositions) { - mass_fractions[symbol] = getMassFraction(symbol); - } - return mass_fractions; -} - - -double Composition::getNumberFraction(const std::string& symbol) const { - if (!m_finalized) { - LOG_ERROR(m_logger, "Composition has not been finalized."); - throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); - } - if (m_compositions.count(symbol) == 0) { - LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); - throw std::runtime_error("Symbol is not in the composition."); - } - if (!m_massFracMode) { - return m_compositions.at(symbol).number_fraction(); - } else { - return m_compositions.at(symbol).number_fraction(m_specificNumberDensity); - } -} - -std::unordered_map Composition::getNumberFraction() const { - std::unordered_map number_fractions; - for (const auto& [symbol, entry] : m_compositions) { - number_fractions[symbol] = getNumberFraction(symbol); - } - return number_fractions; -} - -std::pair Composition::getComposition(const std::string& symbol) const { - if (!m_finalized) { - LOG_ERROR(m_logger, "Composition has not been finalized."); - throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); - } - if (m_compositions.count(symbol) == 0) { - LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); - throw std::runtime_error("Symbol is not in the composition."); - } - return {m_compositions.at(symbol), {m_specificNumberDensity, m_meanParticleMass}}; -} - -std::pair, GlobalComposition> Composition::getComposition() const { - if (!m_finalized) { - LOG_ERROR(m_logger, "Composition has not been finalized."); - throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); - } - return {m_compositions, {m_specificNumberDensity, m_meanParticleMass}}; -} - -Composition Composition::subset(const std::vector& symbols, std::string method) const { - std::array methods = {"norm", "none"}; - - if (std::find(methods.begin(), methods.end(), method) == methods.end()) { - 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 std::runtime_error(errorMessage); + bool Composition::finalize(bool norm) { + bool finalized = false; + if (m_massFracMode) { + finalized = finalizeMassFracMode(norm); + } else { + finalized = finalizeNumberFracMode(norm); + } + if (finalized) { + m_finalized = true; + } + return finalized; } - Composition subsetComposition; - for (const auto& symbol : symbols) { + bool Composition::finalizeMassFracMode(bool norm) { + std::vector mass_fractions; + mass_fractions.reserve(m_compositions.size()); + for (const auto& [_, entry] : m_compositions) { + mass_fractions.push_back(entry.mass_fraction()); + } + if (norm) { + double sum = 0.0; + for (const auto& mass_fraction : mass_fractions) { + sum += mass_fraction; + } + for (int i = 0; i < mass_fractions.size(); ++i) { + mass_fractions[i] /= sum; + } + for (auto& [symbol, entry] : m_compositions) { + setMassFraction(symbol, entry.mass_fraction() / sum); + } + } + try { + validateComposition(mass_fractions); + } catch (const std::runtime_error& e) { + double massSum = 0.0; + for (const auto& [_, entry] : m_compositions) { + massSum += entry.mass_fraction(); + } + LOG_ERROR(m_logger, "Composition is invalid (Total mass {}).", massSum); + m_finalized = false; + return false; + } + for (const auto& [_, entry] : m_compositions) { + m_specificNumberDensity += entry.rel_abundance(); + } + m_meanParticleMass = 1.0/m_specificNumberDensity; + return true; + } + + bool Composition::finalizeNumberFracMode(bool norm) { + std::vector number_fractions; + number_fractions.reserve(m_compositions.size()); + for (const auto& [_, entry] : m_compositions) { + number_fractions.push_back(entry.number_fraction()); + } + if (norm) { + double sum = 0.0; + for (const auto& number_fraction : number_fractions) { + sum += number_fraction; + } + for (auto& [symbol, entry] : m_compositions) { + setNumberFraction(symbol, entry.number_fraction() / sum); + } + } + try { + validateComposition(number_fractions); + } catch (const std::runtime_error& e) { + double numberSum = 0.0; + for (const auto& [_, entry] : m_compositions) { + numberSum += entry.number_fraction(); + } + LOG_ERROR(m_logger, "Composition is invalid (Total number {}).", numberSum); + m_finalized = false; + return false; + } + for (const auto& [_, entry] : m_compositions) { + m_meanParticleMass += entry.rel_abundance(); + } + m_specificNumberDensity = 1.0/m_meanParticleMass; + return true; + } + + Composition Composition::mix(const Composition& other, double fraction) const { + if (!m_finalized || !other.m_finalized) { + LOG_ERROR(m_logger, "Compositions have not both been finalized."); + throw std::runtime_error("Compositions have not been finalized (Consider running .finalize())."); + } + + if (fraction < 0.0 || fraction > 1.0) { + LOG_ERROR(m_logger, "Fraction must be between 0 and 1."); + throw std::runtime_error("Fraction must be between 0 and 1."); + } + + std::set mixedSymbols = other.getRegisteredSymbols(); + // Get the union of the two sets + mixedSymbols.insert(m_registeredSymbols.begin(), m_registeredSymbols.end()); + + Composition mixedComposition(mixedSymbols); + for (const auto& symbol : mixedSymbols) { + double thisMassFrac, otherMassFrac = 0.0; + + thisMassFrac = hasSymbol(symbol) ? getMassFraction(symbol) : 0.0; + otherMassFrac = other.hasSymbol(symbol) ? other.getMassFraction(symbol) : 0.0; + + double massFraction = fraction * thisMassFrac + otherMassFrac * (1-fraction); + mixedComposition.setMassFraction(symbol, massFraction); + } + mixedComposition.finalize(); + return mixedComposition; + } + + double Composition::getMassFraction(const std::string& symbol) const { + if (!m_finalized) { + LOG_ERROR(m_logger, "Composition has not been finalized."); + throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); + } if (m_compositions.count(symbol) == 0) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); + } + if (m_massFracMode) { + return m_compositions.at(symbol).mass_fraction(); } else { - subsetComposition.registerSymbol(symbol); - } - subsetComposition.setMassFraction(symbol, m_compositions.at(symbol).mass_fraction()); - } - if (method == "norm") { - bool isNorm = subsetComposition.finalize(true); - if (!isNorm) { - LOG_ERROR(m_logger, "Subset composition is invalid."); - throw std::runtime_error("Subset composition is invalid."); + return m_compositions.at(symbol).mass_fraction(m_meanParticleMass); } } - return subsetComposition; -} -void Composition::setCompositionMode(bool massFracMode) { - if (!m_finalized) { - LOG_ERROR(m_logger, "Composition has not been finalized. Mode cannot be set unless composition is finalized."); - throw std::runtime_error("Composition has not been finalized (Consider running .finalize()). The mode cannot be set unless the composition is finalized."); + std::unordered_map Composition::getMassFraction() const { + std::unordered_map mass_fractions; + for (const auto& [symbol, entry] : m_compositions) { + mass_fractions[symbol] = getMassFraction(symbol); + } + return mass_fractions; } - bool okay = true; - for (auto& [_, entry] : m_compositions) { - if (massFracMode) { - okay = entry.setMassFracMode(m_meanParticleMass); + + double Composition::getNumberFraction(const std::string& symbol) const { + if (!m_finalized) { + LOG_ERROR(m_logger, "Composition has not been finalized."); + throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); + } + if (m_compositions.count(symbol) == 0) { + LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); + throw std::runtime_error("Symbol is not in the composition."); + } + if (!m_massFracMode) { + return m_compositions.at(symbol).number_fraction(); } else { - okay = entry.setNumberFracMode(m_specificNumberDensity); - } - if (!okay) { - LOG_ERROR(m_logger, "Composition mode could not be set."); - throw std::runtime_error("Composition mode could not be set due to an unknown error."); + return m_compositions.at(symbol).number_fraction(m_specificNumberDensity); } } - m_massFracMode = massFracMode; -} -bool Composition::hasSymbol(const std::string& symbol) const { - return m_compositions.count(symbol) > 0; -} - -/// OVERLOADS - -Composition Composition::operator+(const Composition& other) const { - return mix(other, 0.5); -} - -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) { - os << "<" << entry.m_symbol << " : m_frac = " << entry.mass_fraction() << ">"; - return os; -} - -std::ostream& operator<<(std::ostream& os, const Composition& composition) { - os << "Composition: \n"; - for (const auto& [symbol, entry] : composition.m_compositions) { - os << entry << "\n"; + std::unordered_map Composition::getNumberFraction() const { + std::unordered_map number_fractions; + for (const auto& [symbol, entry] : m_compositions) { + number_fractions[symbol] = getNumberFraction(symbol); + } + return number_fractions; + } + + std::pair Composition::getComposition(const std::string& symbol) const { + if (!m_finalized) { + LOG_ERROR(m_logger, "Composition has not been finalized."); + throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); + } + if (m_compositions.count(symbol) == 0) { + LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); + throw std::runtime_error("Symbol is not in the composition."); + } + return {m_compositions.at(symbol), {m_specificNumberDensity, m_meanParticleMass}}; + } + + std::pair, GlobalComposition> Composition::getComposition() const { + if (!m_finalized) { + LOG_ERROR(m_logger, "Composition has not been finalized."); + throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); + } + return {m_compositions, {m_specificNumberDensity, m_meanParticleMass}}; + } + + double Composition::getMeanParticleMass() const { + if (!m_finalized) { + LOG_ERROR(m_logger, "Composition has not been finalized."); + throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); + } + return m_meanParticleMass; + } + + double Composition::getMeanAtomicNumber() const { + if (!m_finalized) { + LOG_ERROR(m_logger, "Composition must be finalized before getting the mean atomic mass number."); + throw std::runtime_error("Composition not finalized. Cannot retrieve mean atomic mass number."); + } + + double mean_A = 0.0; + + // Loop through all registered species in the composition. + for (const auto &val: m_compositions | std::views::values) { + const CompositionEntry& entry = val; + const chemSpecies::Species& species = entry.isotope(); + + const double mass_fraction = entry.mass_fraction(); + const double particle_mass_g = species.mass(); + const int mass_number = species.a(); + + // Avoid division by zero, though a valid species should have a positive mass. + if (particle_mass_g > 0) { + // Calculate the number fraction for this species. + const double number_fraction = (mass_fraction / particle_mass_g) * m_meanParticleMass; + mean_A += number_fraction * mass_number; + } + } + + return mean_A; + } + + Composition Composition::subset(const std::vector& symbols, std::string method) const { + std::array methods = {"norm", "none"}; + + if (std::find(methods.begin(), methods.end(), method) == methods.end()) { + 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 std::runtime_error(errorMessage); + } + + Composition subsetComposition; + for (const auto& symbol : symbols) { + if (m_compositions.count(symbol) == 0) { + LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); + throw std::runtime_error("Symbol is not in the composition."); + } else { + subsetComposition.registerSymbol(symbol); + } + subsetComposition.setMassFraction(symbol, m_compositions.at(symbol).mass_fraction()); + } + if (method == "norm") { + bool isNorm = subsetComposition.finalize(true); + if (!isNorm) { + LOG_ERROR(m_logger, "Subset composition is invalid."); + throw std::runtime_error("Subset composition is invalid."); + } + } + return subsetComposition; + } + + void Composition::setCompositionMode(bool massFracMode) { + if (!m_finalized) { + LOG_ERROR(m_logger, "Composition has not been finalized. Mode cannot be set unless composition is finalized."); + throw std::runtime_error("Composition has not been finalized (Consider running .finalize()). The mode cannot be set unless the composition is finalized."); + } + + bool okay = true; + for (auto& [_, entry] : m_compositions) { + if (massFracMode) { + okay = entry.setMassFracMode(m_meanParticleMass); + } else { + okay = entry.setNumberFracMode(m_specificNumberDensity); + } + if (!okay) { + LOG_ERROR(m_logger, "Composition mode could not be set."); + throw std::runtime_error("Composition mode could not be set due to an unknown error."); + } + } + m_massFracMode = massFracMode; + } + + bool Composition::hasSymbol(const std::string& symbol) const { + return m_compositions.count(symbol) > 0; + } + + /// OVERLOADS + + Composition Composition::operator+(const Composition& other) const { + return mix(other, 0.5); + } + + 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) { + os << "<" << entry.m_symbol << " : m_frac = " << entry.mass_fraction() << ">"; + return os; + } + + std::ostream& operator<<(std::ostream& os, const Composition& composition) { + os << "Composition: \n"; + for (const auto& [symbol, entry] : composition.m_compositions) { + os << entry << "\n"; + } + return os; } - return os; -} } // namespace serif::composition \ No newline at end of file diff --git a/src/composition/public/composition.h b/src/composition/public/composition.h index e61a5c6..457f96d 100644 --- a/src/composition/public/composition.h +++ b/src/composition/public/composition.h @@ -277,7 +277,7 @@ namespace serif::composition { * Composition comp(symbols); * @endcode */ - Composition(const std::vector& symbols); + explicit Composition(const std::vector& symbols); /** * @brief Constructs a Composition with the given symbols as a set. @@ -288,7 +288,7 @@ namespace serif::composition { * Composition comp(symbols); * @endcode */ - Composition(const std::set& symbols); + explicit Composition(const std::set& symbols); /** * @brief Constructs a Composition with the given symbols and mass fractions. @@ -304,6 +304,14 @@ namespace serif::composition { */ Composition(const std::vector& symbols, const std::vector& mass_fractions, bool massFracMode=true); + /** + * @brief Constructs a Composition from another Composition. + * @param composition The Composition to copy. + */ + Composition(const Composition& composition); + + Composition& operator=(Composition const& other); + /** * @brief Registers a new symbol. * @param symbol The symbol to register. @@ -333,7 +341,7 @@ namespace serif::composition { * @brief Gets the registered symbols. * @return A set of registered symbols. */ - std::set getRegisteredSymbols() const; + [[nodiscard]] std::set getRegisteredSymbols() const; /** * @brief Sets the mass fraction for a given symbol. @@ -390,40 +398,52 @@ namespace serif::composition { * @brief Gets the mass fractions of all compositions. * @return An unordered map of compositions with their mass fractions. */ - std::unordered_map getMassFraction() const; + [[nodiscard]] std::unordered_map getMassFraction() const; /** * @brief Gets the mass fraction for a given symbol. * @param symbol The symbol to get the mass fraction for. * @return The mass fraction for the given symbol. */ - double getMassFraction(const std::string& symbol) const; + [[nodiscard]] double getMassFraction(const std::string& symbol) const; /** * @brief Gets the number fraction for a given symbol. * @param symbol The symbol to get the number fraction for. * @return The number fraction for the given symbol. */ - double getNumberFraction(const std::string& symbol) const; + [[nodiscard]] double getNumberFraction(const std::string& symbol) const; /** * @brief Gets the number fractions of all compositions. * @return An unordered map of compositions with their number fractions. */ - std::unordered_map getNumberFraction() const; + [[nodiscard]] std::unordered_map getNumberFraction() const; /** * @brief Gets the composition entry and global composition for a given symbol. * @param symbol The symbol to get the composition for. * @return A pair containing the CompositionEntry and GlobalComposition for the given symbol. */ - std::pair getComposition(const std::string& symbol) const; + [[nodiscard]] std::pair getComposition(const std::string& symbol) const; /** * @brief Gets all composition entries and the global composition. * @return A pair containing an unordered map of CompositionEntries and the GlobalComposition. */ - std::pair, GlobalComposition> getComposition() const; + [[nodiscard]] std::pair, GlobalComposition> getComposition() const; + + /** + * @brief Compute the mean particle mass of the composition. + * @return Mean particle mass in g. + */ + [[nodiscard]] double getMeanParticleMass() const; + + /** + * @brief Compute the mean atomic mass number of the composition. + * @return Mean atomic mass number. + */ + [[nodiscard]] double getMeanAtomicNumber() const; /** * @brief Gets a subset of the composition. From 36a3f832f79d47d2c8b0580aa855254d15c99aec Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 17 Jun 2025 08:12:41 -0400 Subject: [PATCH 3/4] fix(eos): fixed calculation of mean atomic number --- src/composition/private/composition.cpp | 100 ++++++++++++++++-------- src/composition/public/composition.h | 26 +++++- 2 files changed, 93 insertions(+), 33 deletions(-) diff --git a/src/composition/private/composition.cpp b/src/composition/private/composition.cpp index 9543836..98a3668 100644 --- a/src/composition/private/composition.cpp +++ b/src/composition/private/composition.cpp @@ -255,12 +255,12 @@ namespace serif::composition { return true; } - bool Composition::isValidSymbol(const std::string& symbol) const { - return chemSpecies::species.count(symbol) > 0; + bool Composition::isValidSymbol(const std::string& symbol) { + return chemSpecies::species.contains(symbol); } double Composition::setMassFraction(const std::string& symbol, const double& mass_fraction) { - if (m_registeredSymbols.find(symbol) == m_registeredSymbols.end()) { + if (!m_registeredSymbols.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol); throw std::runtime_error("Symbol is not registered."); } @@ -276,7 +276,7 @@ namespace serif::composition { } m_finalized = false; - double old_mass_fraction = m_compositions.at(symbol).mass_fraction(); + const double old_mass_fraction = m_compositions.at(symbol).mass_fraction(); m_compositions.at(symbol).setMassFraction(mass_fraction); return old_mass_fraction; @@ -432,9 +432,9 @@ namespace serif::composition { Composition mixedComposition(mixedSymbols); for (const auto& symbol : mixedSymbols) { - double thisMassFrac, otherMassFrac = 0.0; + double otherMassFrac = 0.0; - thisMassFrac = hasSymbol(symbol) ? getMassFraction(symbol) : 0.0; + const double thisMassFrac = hasSymbol(symbol) ? getMassFraction(symbol) : 0.0; otherMassFrac = other.hasSymbol(symbol) ? other.getMassFraction(symbol) : 0.0; double massFraction = fraction * thisMassFrac + otherMassFrac * (1-fraction); @@ -449,7 +449,7 @@ namespace serif::composition { LOG_ERROR(m_logger, "Composition has not been finalized."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); } - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } @@ -462,7 +462,7 @@ namespace serif::composition { std::unordered_map Composition::getMassFraction() const { std::unordered_map mass_fractions; - for (const auto& [symbol, entry] : m_compositions) { + for (const auto &symbol: m_compositions | std::views::keys) { mass_fractions[symbol] = getMassFraction(symbol); } return mass_fractions; @@ -474,7 +474,7 @@ namespace serif::composition { LOG_ERROR(m_logger, "Composition has not been finalized."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); } - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } @@ -487,7 +487,7 @@ namespace serif::composition { std::unordered_map Composition::getNumberFraction() const { std::unordered_map number_fractions; - for (const auto& [symbol, entry] : m_compositions) { + for (const auto &symbol: m_compositions | std::views::keys) { number_fractions[symbol] = getNumberFraction(symbol); } return number_fractions; @@ -498,7 +498,7 @@ namespace serif::composition { LOG_ERROR(m_logger, "Composition has not been finalized."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); } - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } @@ -527,40 +527,29 @@ namespace serif::composition { throw std::runtime_error("Composition not finalized. Cannot retrieve mean atomic mass number."); } - double mean_A = 0.0; + double zSum = 0.0; // Loop through all registered species in the composition. for (const auto &val: m_compositions | std::views::values) { - const CompositionEntry& entry = val; - const chemSpecies::Species& species = entry.isotope(); - - const double mass_fraction = entry.mass_fraction(); - const double particle_mass_g = species.mass(); - const int mass_number = species.a(); - - // Avoid division by zero, though a valid species should have a positive mass. - if (particle_mass_g > 0) { - // Calculate the number fraction for this species. - const double number_fraction = (mass_fraction / particle_mass_g) * m_meanParticleMass; - mean_A += number_fraction * mass_number; - } + zSum += (val.mass_fraction() * val.m_isotope.z())/val.m_isotope.a(); } + const double mean_A = m_meanParticleMass * zSum; return mean_A; } Composition Composition::subset(const std::vector& symbols, std::string method) const { - std::array methods = {"norm", "none"}; + const std::array methods = {"norm", "none"}; - if (std::find(methods.begin(), methods.end(), method) == methods.end()) { - std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'."; + if (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 std::runtime_error(errorMessage); } Composition subsetComposition; for (const auto& symbol : symbols) { - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } else { @@ -569,7 +558,7 @@ namespace serif::composition { subsetComposition.setMassFraction(symbol, m_compositions.at(symbol).mass_fraction()); } if (method == "norm") { - bool isNorm = subsetComposition.finalize(true); + const bool isNorm = subsetComposition.finalize(true); if (!isNorm) { LOG_ERROR(m_logger, "Subset composition is invalid."); throw std::runtime_error("Subset composition is invalid."); @@ -578,14 +567,14 @@ namespace serif::composition { return subsetComposition; } - void Composition::setCompositionMode(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."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize()). The mode cannot be set unless the composition is finalized."); } bool okay = true; - for (auto& [_, entry] : m_compositions) { + for (auto &entry: m_compositions | std::views::values) { if (massFracMode) { okay = entry.setMassFracMode(m_meanParticleMass); } else { @@ -599,6 +588,53 @@ namespace serif::composition { m_massFracMode = massFracMode; } + CanonicalComposition Composition::getCanonicalComposition(bool harsh) const { + if (!m_finalized) { + LOG_ERROR(m_logger, "Composition has not been finalized."); + throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); + } + CanonicalComposition canonicalComposition; + constexpr std::array canonicalH = { + "H-1", "H-2", "H-3", "H-4", "H-5", "H-6", "H-7" + }; + constexpr std::array canonicalHe = { + "He-3", "He-4", "He-5", "He-6", "He-7", "He-8", "He-9", "He-10" + }; + for (const auto& symbol : canonicalH) { + if (hasSymbol(symbol)) { + canonicalComposition.X += getMassFraction(symbol); + } + } + for (const auto& symbol : canonicalHe) { + if (hasSymbol(symbol)) { + canonicalComposition.Y += getMassFraction(symbol); + } + } + + for (const auto& symbol : getRegisteredSymbols()) { + const bool isHSymbol = std::ranges::find(canonicalH, symbol) != std::end(canonicalH); + const bool isHeSymbol = std::ranges::find(canonicalHe, symbol) != std::end(canonicalHe); + + if (isHSymbol || isHeSymbol) { + continue; // Skip canonical H and He symbols + } + + canonicalComposition.Z += getMassFraction(symbol); + } + + const double Z = 1.0 - (canonicalComposition.X + canonicalComposition.Y); + if (std::abs(Z - canonicalComposition.Z) > 1e-6) { + if (!harsh) { + LOG_WARNING(m_logger, "Validation composition Z (X-Y = {}) is different than canonical composition Z ({}) (∑a_i where a_i != H/He).", Z, canonicalComposition.Z); + } + else { + LOG_ERROR(m_logger, "Validation composition Z (X-Y = {}) is different than canonical composition Z ({}) (∑a_i where a_i != H/He).", Z, canonicalComposition.Z); + 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)."); + } + } + return canonicalComposition; + } + bool Composition::hasSymbol(const std::string& symbol) const { return m_compositions.count(symbol) > 0; } diff --git a/src/composition/public/composition.h b/src/composition/public/composition.h index 457f96d..861cc47 100644 --- a/src/composition/public/composition.h +++ b/src/composition/public/composition.h @@ -34,6 +34,20 @@ #include "atomicSpecies.h" namespace serif::composition { + struct CanonicalComposition { + double X = 0.0; ///< Mass fraction of Hydrogen. + double Y = 0.0; ///< Mass fraction of Helium. + double Z = 0.0; ///< Mass fraction of Metals. + + friend std::ostream& operator<<(std::ostream& os, const CanonicalComposition& composition) { + os << ""; + return os; + } + }; + /** * @brief Represents the global composition of a system. This tends to be used after finalize and is primarily for internal use. */ @@ -220,7 +234,7 @@ namespace serif::composition { * @param symbol The symbol to check. * @return True if the symbol is valid, false otherwise. */ - bool isValidSymbol(const std::string& symbol) const; + static bool isValidSymbol(const std::string& symbol); /** * @brief Checks if the given mass fractions are valid. @@ -466,6 +480,16 @@ namespace serif::composition { */ void setCompositionMode(bool massFracMode); + /** + * @brief Gets the current canonical composition (X, Y, Z). + * @param harsh If true, this will throw an error if X-Y != Z where Z is computed as the sum of all other elements. + * @return True if mass fraction mode, false if number fraction mode. + * + * @throws std::runtime_error if the composition is not finalized or if the canonical composition cannot be computed. + * @throws std::runtime_error if harsh is true and the canonical composition is not valid. + */ + [[nodiscard]] CanonicalComposition getCanonicalComposition(bool harsh=false) const; + /** * @brief Overloaded output stream operator for Composition. * @param os The output stream. From 6d517d937e4d431c4d6e04f51d397e0e9dd74255 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 17 Jun 2025 09:43:43 -0400 Subject: [PATCH 4/4] refactor(network): updated network and network::approx8 to use composition module This is a very basic wrapper implimentation currently. This is sufficient to lock the interface down so that other code can target it. However, internally there is just a "convert" function. Eventually we should rework the code itself to use the composition module more directly. --- assets/static/atomic/meson.build | 3 ++- src/composition/meson.build | 2 +- src/composition/private/composition.cpp | 2 +- src/composition/public/composition.h | 7 ++----- 4 files changed, 6 insertions(+), 8 deletions(-) diff --git a/assets/static/atomic/meson.build b/assets/static/atomic/meson.build index 9e024aa..b5b79b7 100644 --- a/assets/static/atomic/meson.build +++ b/assets/static/atomic/meson.build @@ -1,3 +1,4 @@ species_weight_dep = declare_dependency( include_directories: include_directories('include'), -) \ No newline at end of file +) +message('✅ SERiF species_weight dependency declared') \ No newline at end of file diff --git a/src/composition/meson.build b/src/composition/meson.build index fe928d8..d97cbec 100644 --- a/src/composition/meson.build +++ b/src/composition/meson.build @@ -29,4 +29,4 @@ composition_dep = declare_dependency( ) # Make headers accessible -install_headers(composition_headers, subdir : '4DSSE/composition') \ No newline at end of file +install_headers(composition_headers, subdir : 'SERiF/composition') \ No newline at end of file diff --git a/src/composition/private/composition.cpp b/src/composition/private/composition.cpp index 98a3668..aad84dc 100644 --- a/src/composition/private/composition.cpp +++ b/src/composition/private/composition.cpp @@ -333,7 +333,7 @@ namespace serif::composition { return old_number_fractions; } - bool Composition::finalize(bool norm) { + bool Composition::finalize(const bool norm) { bool finalized = false; if (m_massFracMode) { finalized = finalizeMassFracMode(norm); diff --git a/src/composition/public/composition.h b/src/composition/public/composition.h index 861cc47..3dcf8a0 100644 --- a/src/composition/public/composition.h +++ b/src/composition/public/composition.h @@ -18,8 +18,7 @@ // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // // *********************************************************************** */ -#ifndef COMPOSITION_H -#define COMPOSITION_H +#pragma once #include #include @@ -507,6 +506,4 @@ namespace serif::composition { Composition operator+(const Composition& other) const; }; -}; // namespace serif::composition - -#endif // COMPOSITION_H \ No newline at end of file +}; // namespace serif::composition \ No newline at end of file