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. diff --git a/src/eos/meson.build b/src/eos/meson.build index 00e285c..43d66e9 100644 --- a/src/eos/meson.build +++ b/src/eos/meson.build @@ -1,15 +1,18 @@ # Define the library eos_sources = files( 'private/helm.cpp', - 'private/EOSio.cpp' + 'private/EOSio.cpp', + 'private/EOS.cpp' ) eos_headers = files( 'public/helm.h', - 'public/EOSio.h' + 'public/EOSio.h', + 'public/EOS.h' ) dependencies = [ + composition_dep, const_dep, quill_dep, probe_dep, diff --git a/src/eos/private/EOS.cpp b/src/eos/private/EOS.cpp new file mode 100644 index 0000000..ba3593d --- /dev/null +++ b/src/eos/private/EOS.cpp @@ -0,0 +1,64 @@ +#include "EOS.h" +#include "EOSio.h" +#include "helm.h" +#include + +namespace serif::eos { + EOS::EOS(const EOSio& reader) : m_reader(reader) {} + EOS::EOS(const std::string& filename, const EOSFormat format) : m_reader(EOSio(filename, format)) {} + + EOSOutput EOS::get(const EOSInput& in) { + EOSOutput output; + if (getFormat() == EOSFormat::HELM) { + helmholtz::HELMEOSInput q; + q.T = in.temperature; // Temperature in K + q.rho = in.density; // Density in g/cm^3 + + q.abar = in.composition.getMeanParticleMass(); // Mean atomic mass in g + q.zbar = in.composition.getMeanAtomicNumber(); // Mean atomic number (dimensionless) + + helmholtz::HELMEOSOutput tempOutput; + tempOutput = helmholtz::get_helm_EOS(q, *std::get>(m_reader.getTable())); + + output.electronFraction = tempOutput.ye; + output.electronChemicalPotential = tempOutput.etaele; + output.neutronExcessFraction = tempOutput.xnefer; + + // --- Pressure Variables --- + output.pressure.total = tempOutput.ptot; + output.pressure.gas = tempOutput.pgas; + output.pressure.radiation = tempOutput.ptot; + output.pressure.dDensity = tempOutput.dpresdd; + output.pressure.dTemperature = tempOutput.dpresdt; + output.pressure.dMeanAtomicMassNumber = tempOutput.dpresda; + output.pressure.dMeanAtomicNumber = tempOutput.dpresdz; + + // --- Energy Variables --- + output.energy.total = tempOutput.etot; + output.energy.gas = tempOutput.egas; + output.energy.radiation = tempOutput.erad; + output.energy.dDensity = tempOutput.denerdd; + output.energy.dTemperature = tempOutput.denerdt; + output.energy.dMeanAtomicMassNumber = tempOutput.denerda; + output.energy.dMeanAtomicNumber = tempOutput.denerdz; + + // --- Entropy Variables --- + output.entropy.total = tempOutput.stot; + output.entropy.gas = tempOutput.sgas; + output.entropy.radiation = tempOutput.srad; + output.entropy.dDensity = tempOutput.dentrdd; + output.entropy.dTemperature = tempOutput.dentrdt; + output.entropy.dMeanAtomicMassNumber = tempOutput.dentrda; + output.entropy.dMeanAtomicNumber = tempOutput.dentrdz; + } + return output; + } + + EOSFormat EOS::getFormat() const { + return m_reader.getFormat(); + } + + const EOSio& EOS::getReader() const { + return m_reader; + } +} \ No newline at end of file diff --git a/src/eos/private/EOSio.cpp b/src/eos/private/EOSio.cpp index 03a6b43..6945e1a 100644 --- a/src/eos/private/EOSio.cpp +++ b/src/eos/private/EOSio.cpp @@ -28,25 +28,32 @@ #include namespace serif::eos { - EOSio::EOSio(const std::string &filename) : m_filename(filename) { + EOSio::EOSio(const std::string &filename, const EOSFormat format) : m_filename(filename), m_format(format){ load(); } - std::string EOSio::getFormat() const { + EOSio::EOSio(const EOSio &other) { + m_filename = other.m_filename; + m_format = other.m_format; + load(); + } + + EOSFormat EOSio::getFormat() const { return m_format; } + std::string EOSio::getFormatName() const { + return FormatStringLookup.at(m_format); + } + + EOSTable& EOSio::getTable() { return m_table; } void EOSio::load() { - // Load the EOS table from the file - // For now, just set the format to HELM - - m_format = "helm"; - if (m_format == "helm") { + if (m_format == EOSFormat::HELM) { loadHelm(); } } diff --git a/src/eos/private/helm.cpp b/src/eos/private/helm.cpp index 1c18b05..e2df9de 100644 --- a/src/eos/private/helm.cpp +++ b/src/eos/private/helm.cpp @@ -240,7 +240,7 @@ namespace serif::eos::helmholtz { ion, radiation, electron-positron and Coulomb interaction and returns the calculated quantities in the input ***/ - serif::eos::helmholtz::EOS get_helm_EOS(serif::eos::helmholtz::EOSInput &q, const serif::eos::helmholtz::HELMTable &table) { + serif::eos::helmholtz::HELMEOSOutput get_helm_EOS(serif::eos::helmholtz::HELMEOSInput &q, const serif::eos::helmholtz::HELMTable &table) { serif::config::Config& config = serif::config::Config::getInstance(); auto logFile = config.get("EOS:Helm:LogFile", "log"); serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance(); @@ -829,7 +829,7 @@ namespace serif::eos::helmholtz { double csound = clight * sqrt(gamma1/z); // package in q: - serif::eos::helmholtz::EOS eos; + serif::eos::helmholtz::HELMEOSOutput eos; eos.ptot = ptot; eos.etot = etot; eos.stot = stot; eos.pgas = pgas; eos.egas = egas; eos.sgas = sgas; eos.prad = prad; eos.erad = erad; eos.srad = srad; diff --git a/src/eos/public/EOS.h b/src/eos/public/EOS.h new file mode 100644 index 0000000..b2ecd22 --- /dev/null +++ b/src/eos/public/EOS.h @@ -0,0 +1,230 @@ +#pragma once + +#include "EOSio.h" +#include "helm.h" +#include +#include "composition.h" + +namespace serif::eos { + + /** + * @brief Input parameters for an EOS calculation. + * + * This struct holds the necessary physical conditions (composition, density, temperature) + * required to query the Equation of State. + */ + struct EOSInput { + serif::composition::Composition composition; ///< The composition of the system. + double density; ///< The density of the system in cgs (g/cm^3). + double temperature; ///< The temperature of the system in cgs (K). + }; + + /** + * @brief Represents a thermodynamic parameter and its derivatives. + * + * This struct stores a specific thermodynamic quantity (e.g., pressure, energy, entropy), + * its breakdown into gas and radiation components, and its partial derivatives + * with respect to density, temperature, mean atomic mass number, and mean atomic number. + * All values are in cgs units unless otherwise specified. + */ + struct EOSParameter { + double total; ///< Total value of the parameter (gas + radiation) (cgs). + double gas; ///< Gas contribution to the parameter (cgs). + double radiation; ///< Radiation contribution to the parameter (cgs). + + double dDensity; ///< Derivative of the total parameter with respect to density (cgs units / (g/cm^3)). + double dTemperature; ///< Derivative of the total parameter with respect to temperature (cgs units / K). + double dMeanAtomicMassNumber; ///< Derivative of the total parameter with respect to mean atomic mass number (Abar) (cgs units / (g/mol)). + double dMeanAtomicNumber; ///< Derivative of the total parameter with respect to mean atomic number (Zbar) (cgs units / dimensionless). + + std::string name; ///< Name of the parameter (e.g., "Pressure", "Energy", "Entropy"). + }; + + /** + * @brief Output from an EOS calculation. + * + * This struct contains various thermodynamic quantities and their derivatives + * calculated by the EOS for a given set of input conditions. + * It includes fundamental properties like electron fraction and chemical potential, + * as well as detailed breakdowns of pressure, energy, and entropy. + * Additionally, it provides methods to calculate derived quantities like + * susceptibilities, sound speed, adiabatic gradients, and specific heats. + */ + struct EOSOutput { + double electronFraction{}; ///< Electron fraction (ye), dimensionless. + double electronChemicalPotential{}; ///< Electron chemical potential (eta_e) in cgs (erg/g). + double neutronExcessFraction{}; ///< Neutron excess fraction (xnefer), dimensionless. + + EOSParameter pressure; ///< Pressure output data, including total, gas, radiation, and derivatives. + EOSParameter energy; ///< Internal energy output data, including total, gas, radiation, and derivatives. + EOSParameter entropy; ///< Entropy output data, including total, gas, radiation, and derivatives. + + /** + * @brief Calculates the temperature susceptibility (chi_T). + * @return Temperature susceptibility, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, chi_T = (d ln P / d ln T)_rho. + */ + double chiTemperature(); + /** + * @brief Calculates the density susceptibility (chi_rho). + * @return Density susceptibility, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, chi_rho = (d ln P / d ln rho)_T. + */ + double chiRho(); + /** + * @brief Calculates the adiabatic sound speed. + * @return Sound speed in cgs (cm/s). + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, c_s^2 = gamma1 * P / rho. + */ + double soundSpeed(); + /** + * @brief Calculates the adiabatic temperature gradient (nabla_ad). + * @return Adiabatic gradient, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, nabla_ad = (P * chi_T) / (rho * T * c_p * chi_rho). + */ + double adiabaticGradient(); + /** + * @brief Calculates the first adiabatic index (Gamma1). + * @return First adiabatic index, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, Gamma1 = (d ln P / d ln rho)_S. + */ + double gamma1(); + /** + * @brief Calculates the second adiabatic index (Gamma2). + * @return Second adiabatic index, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, Gamma2 / (Gamma2 - 1) = (d ln P / d ln T)_S. + */ + double gamma2(); + /** + * @brief Calculates the third adiabatic index (Gamma3). + * @return Third adiabatic index, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, Gamma3 - 1 = (d ln T / d ln rho)_S. + */ + double gamma3(); + /** + * @brief Calculates the specific heat capacity at constant volume (c_v). + * @return Specific heat capacity at constant volume in cgs (erg/K/g). + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, c_v = (dE/dT)_rho. + */ + double specificHeatCapacityAtConstantVolume(); + /** + * @brief Calculates the specific heat capacity at constant pressure (c_p). + * @return Specific heat capacity at constant pressure in cgs (erg/K/g). + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, c_p = c_v + (T / rho^2) * ( (d P / d T)_rho^2 / (d P / d rho)_T ). + */ + double specificHeatCapacityAtConstantPressure(); + + /** + * @brief Returns the format of the EOS data used to generate this output. + * @return The EOSFormat enum value (currently only EOSFormat::HELM). + */ + EOSFormat EOSFormat() const; + }; + + /** + * @class EOS + * @brief Main class for accessing Equation of State data. + * + * This class provides an interface to an underlying EOS table (e.g., Helmholtz EOS). + * It handles loading the EOS data and provides a method to retrieve thermodynamic + * properties for given physical conditions. + * + * @example + * @code + * #include "EOS.h" + * #include "composition.h" // For serif::composition::Composition + * #include + * + * int main() { + * try { + * // Initialize EOS from a Helmholtz table file + * serif::eos::EOS helmEOS("path/to/helm_table.dat", serif::eos::EOSFormat::HELM); + * + * // Define input conditions + * serif::eos::EOSInput input; + * input.density = 1.0e6; // g/cm^3 + * input.temperature = 1.0e7; // K + * // Assuming a simple composition (e.g., pure Helium-4 for demonstration) + * // In a real scenario, initialize Composition properly. + * // For example, if Composition has a constructor like: + * // Composition(const std::map, double>& mass_fractions); + * // std::map, double> he4_mass_fraction = {{{2, 4}, 1.0}}; + * // input.composition = serif::composition::Composition(he4_mass_fraction); + * // For now, let's assume Composition can be default constructed or set up simply: + * input.composition.addSpecies(2, 4, 1.0); // Z=2, A=4 (He-4), mass fraction 1.0 + * + * // Get EOS output + * serif::eos::EOSOutput output = helmEOS.get(input); + * + * // Access results + * std::cout << "Pressure (total): " << output.pressure.total << " dyne/cm^2" << std::endl; + * std::cout << "Energy (total): " << output.energy.total << " erg/g" << std::endl; + * std::cout << "Entropy (total): " << output.entropy.total << " erg/K/g" << std::endl; + * std::cout << "Electron fraction: " << output.electronFraction << std::endl; + * + * // Example of accessing derivatives + * std::cout << "dP/dRho: " << output.pressure.dDensity << std::endl; + * std::cout << "dE/dT: " << output.energy.dTemperature << std::endl; + * + * } catch (const std::exception& e) { + * std::cerr << "An error occurred: " << e.what() << std::endl; + * return 1; + * } + * return 0; + * } + * @endcode + */ + class EOS { + public: + /** + * @brief Constructs an EOS object by loading data from a file. + * @param filename The path to the EOS data file. + * @param format The format of the EOS data file (e.g., EOSFormat::HELM). + * @throw std::runtime_error If the file cannot be opened or read, or if the format is unsupported. + */ + explicit EOS(const std::string& filename, EOSFormat format=EOSFormat::HELM); + /** + * @brief Constructs an EOS object from an existing EOSio reader. + * @param reader An EOSio object that has already loaded the EOS data. + */ + explicit EOS(const EOSio& reader); + /** + * @brief Default destructor. + */ + ~EOS() = default; + + /** + * @brief Retrieves thermodynamic properties for the given input conditions. + * @param in An EOSInput struct containing the density, temperature, and composition. + * @return An EOSOutput struct containing the calculated thermodynamic properties. + * @throw std::runtime_error If the underlying EOS calculation fails (e.g., out of table bounds for Helmholtz). + * + * This method queries the loaded EOS table (e.g., Helmholtz) using the provided + * density, temperature, and composition (mean atomic mass Abar, mean atomic number Zbar). + * It populates and returns an EOSOutput struct with various thermodynamic quantities + * such as pressure, energy, entropy, their derivatives, electron fraction, etc. + */ + [[nodiscard]] EOSOutput get(const EOSInput& in); + /** + * @brief Gets the format of the loaded EOS data. + * @return The EOSFormat enum value. + */ + [[nodiscard]] EOSFormat getFormat() const; + /** + * @brief Gets a constant reference to the internal EOSio reader. + * @return A const reference to the EOSio object. + */ + [[nodiscard]] const EOSio& getReader() const; + private: + EOSio m_reader; ///< The EOS I/O handler responsible for reading and storing EOS table data. + }; +} diff --git a/src/eos/public/EOSio.h b/src/eos/public/EOSio.h index 1882c07..fff0b64 100644 --- a/src/eos/public/EOSio.h +++ b/src/eos/public/EOSio.h @@ -22,15 +22,20 @@ #include #include #include +#include #include "helm.h" namespace serif::eos { + enum EOSFormat { + HELM, ///< Helmholtz EOS format. + }; + static inline std::unordered_map FormatStringLookup = { + {HELM, "Helmholtz"} + }; // EOS table format includes - using EOSTable = std::variant< - std::unique_ptr - >; + using EOSTable = std::variant>; /** * @class EOSio @@ -41,16 +46,19 @@ namespace serif::eos { * * Example usage: * @code - * EosIO eosIO("path/to/file"); - * std::string format = eosIO.getFormat(); - * EOSTable& table = eosIO.getTable(); + * EOSio eosReader("path/to/file"); + * std::string format = eosReader.getFormatName(); + * EOSTable& table = eosReader.getTable(); * @endcode + * + * @note The default format used for reading tables is HELM + * @note Currently only the HELM format is implemented */ class EOSio { private: std::string m_filename; ///< The filename of the EOS table. bool m_loaded = false; ///< Flag indicating if the table is loaded. - std::string m_format; ///< The format of the EOS table. + EOSFormat m_format; EOSTable m_table; ///< The EOS table data. /** @@ -66,8 +74,15 @@ namespace serif::eos { /** * @brief Constructs an EosIO object with the given filename. * @param filename The filename of the EOS table. + * @param format The EOS file format (currently only HELM) */ - explicit EOSio(const std::string &filename); + explicit EOSio(const std::string &filename, EOSFormat format = EOSFormat::HELM); + + /** + * @brief Explicit copy constructor + * @param other The EOSio to be copied + */ + EOSio(const EOSio& other); /** * @brief Default destructor. @@ -75,10 +90,12 @@ namespace serif::eos { ~EOSio() = default; /** - * @brief Gets the format of the EOS table. + * @brief Gets the format name (as a string) of the EOS table. * @return The format of the EOS table as a string. */ - [[nodiscard]] std::string getFormat() const; + [[nodiscard]] std::string getFormatName() const; + + [[nodiscard]] EOSFormat getFormat() const; /** * @brief Gets the EOS table. @@ -87,6 +104,8 @@ namespace serif::eos { [[nodiscard]] EOSTable& getTable(); [[nodiscard]] std::string getFilename() const { return m_filename; } + + bool isLoaded() const { return m_loaded; } }; } diff --git a/src/eos/public/helm.h b/src/eos/public/helm.h index 501db49..a20ba49 100644 --- a/src/eos/public/helm.h +++ b/src/eos/public/helm.h @@ -169,18 +169,18 @@ namespace serif::eos::helmholtz { }; /** - * @struct EOSInput + * @struct HELMEOSInput * @brief Structure to hold the input parameters for the EOS calculation. */ - struct EOSInput + struct HELMEOSInput { double T; ///< Temperature. double rho; ///< Density. double abar; ///< Mean atomic mass. double zbar; ///< Mean atomic number. - friend std::ostream& operator<<(std::ostream& os, const helmholtz::EOSInput& eosInput) { - os << "EOSInput Data:\n"; + friend std::ostream& operator<<(std::ostream& os, const helmholtz::HELMEOSInput& eosInput) { + os << "HELMEOSInput Data:\n"; os << " Temperature: " << eosInput.T << "\n"; os << " Density: " << eosInput.rho << "\n"; os << " Mean Atomic Mass: " << eosInput.abar << "\n"; @@ -194,7 +194,7 @@ namespace serif::eos::helmholtz { * @struct EOS * @brief Structure to hold the output parameters and derivatives of the EOS calculation. */ - struct EOS + struct HELMEOSOutput { // output double ye, etaele, xnefer; // @@ -212,7 +212,7 @@ namespace serif::eos::helmholtz { double gamma1, gamma2, gamma3, cV, cP; // derived quantities double dse, dpe, dsp; // Maxwell relations - friend std::ostream& operator<<(std::ostream& os, const helmholtz::EOS& eos) { + friend std::ostream& operator<<(std::ostream& os, const helmholtz::HELMEOSOutput& eos) { os << "EOS Data:\n" << std::setw(20) << std::left; os << " Electron Fraction: " << std::format("{0:24.16e}",eos.ye) << "\n"; os << " Electron Chemical Potential: " << std::format("{0:24.16e}",eos.etaele) << "\n"; @@ -379,7 +379,7 @@ namespace serif::eos::helmholtz { * @param table HELMTable structure containing the table data. * @return EOS structure containing the calculated quantities. */ - EOS get_helm_EOS(EOSInput &q, const HELMTable &table); + HELMEOSOutput get_helm_EOS(HELMEOSInput &q, const HELMTable &table); } diff --git a/src/python/eos/bindings.cpp b/src/python/eos/bindings.cpp index e8c6314..ca4b87c 100644 --- a/src/python/eos/bindings.cpp +++ b/src/python/eos/bindings.cpp @@ -34,7 +34,7 @@ void register_eos_bindings(pybind11::module &eos_submodule) { }, py::return_value_policy::reference_internal, // IMPORTANT: Keep this policy! "Get the EOS table data.") .def("__repr__", [](const serif::eos::EOSio &eos) { - return ""; + return ""; }); py::class_(eos_submodule, "EOSTable"); @@ -88,64 +88,64 @@ void register_eos_bindings(pybind11::module &eos_submodule) { ); }, py::return_value_policy::reference_internal); // Keep parent 'table' alive - py::class_(eos_submodule, "EOS") + py::class_(eos_submodule, "EOS") .def(py::init<>()) - .def_readonly("ye", &serif::eos::helmholtz::EOS::ye) - .def_readonly("etaele", &serif::eos::helmholtz::EOS::etaele) - .def_readonly("xnefer", &serif::eos::helmholtz::EOS::xnefer) + .def_readonly("ye", &serif::eos::helmholtz::HELMEOSOutput::ye) + .def_readonly("etaele", &serif::eos::helmholtz::HELMEOSOutput::etaele) + .def_readonly("xnefer", &serif::eos::helmholtz::HELMEOSOutput::xnefer) - .def_readonly("ptot", &serif::eos::helmholtz::EOS::ptot) - .def_readonly("pgas", &serif::eos::helmholtz::EOS::pgas) - .def_readonly("prad", &serif::eos::helmholtz::EOS::prad) + .def_readonly("ptot", &serif::eos::helmholtz::HELMEOSOutput::ptot) + .def_readonly("pgas", &serif::eos::helmholtz::HELMEOSOutput::pgas) + .def_readonly("prad", &serif::eos::helmholtz::HELMEOSOutput::prad) - .def_readonly("etot", &serif::eos::helmholtz::EOS::etot) - .def_readonly("egas", &serif::eos::helmholtz::EOS::egas) - .def_readonly("erad", &serif::eos::helmholtz::EOS::erad) + .def_readonly("etot", &serif::eos::helmholtz::HELMEOSOutput::etot) + .def_readonly("egas", &serif::eos::helmholtz::HELMEOSOutput::egas) + .def_readonly("erad", &serif::eos::helmholtz::HELMEOSOutput::erad) - .def_readonly("stot", &serif::eos::helmholtz::EOS::stot) - .def_readonly("sgas", &serif::eos::helmholtz::EOS::sgas) - .def_readonly("srad", &serif::eos::helmholtz::EOS::srad) + .def_readonly("stot", &serif::eos::helmholtz::HELMEOSOutput::stot) + .def_readonly("sgas", &serif::eos::helmholtz::HELMEOSOutput::sgas) + .def_readonly("srad", &serif::eos::helmholtz::HELMEOSOutput::srad) - .def_readonly("dpresdd", &serif::eos::helmholtz::EOS::dpresdd) - .def_readonly("dpresdt", &serif::eos::helmholtz::EOS::dpresdt) - .def_readonly("dpresda", &serif::eos::helmholtz::EOS::dpresda) - .def_readonly("dpresdz", &serif::eos::helmholtz::EOS::dpresdz) - // TODO: Finish adding all the derivatives to the bound class - .def_readonly("dentrdd", &serif::eos::helmholtz::EOS::dentrdd) - .def_readonly("dentrdt", &serif::eos::helmholtz::EOS::dentrdt) - .def_readonly("dentrda", &serif::eos::helmholtz::EOS::dentrda) - .def_readonly("dentrdz", &serif::eos::helmholtz::EOS::dentrdz) + .def_readonly("dpresdd", &serif::eos::helmholtz::HELMEOSOutput::dpresdd) + .def_readonly("dpresdt", &serif::eos::helmholtz::HELMEOSOutput::dpresdt) + .def_readonly("dpresda", &serif::eos::helmholtz::HELMEOSOutput::dpresda) + .def_readonly("dpresdz", &serif::eos::helmholtz::HELMEOSOutput::dpresdz) - .def_readonly("denerdd", &serif::eos::helmholtz::EOS::denerdd) - .def_readonly("denerdt", &serif::eos::helmholtz::EOS::denerdt) - .def_readonly("denerda", &serif::eos::helmholtz::EOS::denerda) - .def_readonly("denerdz", &serif::eos::helmholtz::EOS::denerdz) + .def_readonly("dentrdd", &serif::eos::helmholtz::HELMEOSOutput::dentrdd) + .def_readonly("dentrdt", &serif::eos::helmholtz::HELMEOSOutput::dentrdt) + .def_readonly("dentrda", &serif::eos::helmholtz::HELMEOSOutput::dentrda) + .def_readonly("dentrdz", &serif::eos::helmholtz::HELMEOSOutput::dentrdz) - .def_readonly("chiT", &serif::eos::helmholtz::EOS::chiT) - .def_readonly("chiRho", &serif::eos::helmholtz::EOS::chiRho) - .def_readonly("csound", &serif::eos::helmholtz::EOS::csound) - .def_readonly("grad_ad", &serif::eos::helmholtz::EOS::grad_ad) - .def_readonly("gamma1", &serif::eos::helmholtz::EOS::gamma1) - .def_readonly("gamma2", &serif::eos::helmholtz::EOS::gamma2) - .def_readonly("gamma3", &serif::eos::helmholtz::EOS::gamma3) - .def_readonly("cV", &serif::eos::helmholtz::EOS::cV) - .def_readonly("cP", &serif::eos::helmholtz::EOS::cP) - .def_readonly("dse", &serif::eos::helmholtz::EOS::dse) - .def_readonly("dpe", &serif::eos::helmholtz::EOS::dpe) - .def_readonly("dsp", &serif::eos::helmholtz::EOS::dsp) + .def_readonly("denerdd", &serif::eos::helmholtz::HELMEOSOutput::denerdd) + .def_readonly("denerdt", &serif::eos::helmholtz::HELMEOSOutput::denerdt) + .def_readonly("denerda", &serif::eos::helmholtz::HELMEOSOutput::denerda) + .def_readonly("denerdz", &serif::eos::helmholtz::HELMEOSOutput::denerdz) - .def("__repr__", [](const serif::eos::helmholtz::EOS &eos) { + .def_readonly("chiT", &serif::eos::helmholtz::HELMEOSOutput::chiT) + .def_readonly("chiRho", &serif::eos::helmholtz::HELMEOSOutput::chiRho) + .def_readonly("csound", &serif::eos::helmholtz::HELMEOSOutput::csound) + .def_readonly("grad_ad", &serif::eos::helmholtz::HELMEOSOutput::grad_ad) + .def_readonly("gamma1", &serif::eos::helmholtz::HELMEOSOutput::gamma1) + .def_readonly("gamma2", &serif::eos::helmholtz::HELMEOSOutput::gamma2) + .def_readonly("gamma3", &serif::eos::helmholtz::HELMEOSOutput::gamma3) + .def_readonly("cV", &serif::eos::helmholtz::HELMEOSOutput::cV) + .def_readonly("cP", &serif::eos::helmholtz::HELMEOSOutput::cP) + .def_readonly("dse", &serif::eos::helmholtz::HELMEOSOutput::dse) + .def_readonly("dpe", &serif::eos::helmholtz::HELMEOSOutput::dpe) + .def_readonly("dsp", &serif::eos::helmholtz::HELMEOSOutput::dsp) + + .def("__repr__", [](const serif::eos::helmholtz::HELMEOSOutput &eos) { return ""; }); - py::class_(eos_submodule, "EOSInput") + py::class_(eos_submodule, "HELMEOSInput") .def(py::init<>()) - .def_readwrite("T", &serif::eos::helmholtz::EOSInput::T) - .def_readwrite("rho", &serif::eos::helmholtz::EOSInput::rho) - .def_readwrite("abar", &serif::eos::helmholtz::EOSInput::abar) - .def_readwrite("zbar", &serif::eos::helmholtz::EOSInput::zbar) - .def("__repr__", [](const serif::eos::helmholtz::EOSInput &input) { - return "";