From 3961c745e3660212243d3124952ca2fd798d2c9f Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 17 Jun 2025 08:12:41 -0400 Subject: [PATCH] fix(eos): fixed calculation of mean atomic number --- src/composition/private/composition.cpp | 100 ++++++++++++++++-------- src/composition/public/composition.h | 26 +++++- src/eos/private/EOS.cpp | 2 + src/eos/public/EOS.h | 42 +++++++++- 4 files changed, 134 insertions(+), 36 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. diff --git a/src/eos/private/EOS.cpp b/src/eos/private/EOS.cpp index ba3593d..0ec82f9 100644 --- a/src/eos/private/EOS.cpp +++ b/src/eos/private/EOS.cpp @@ -17,6 +17,8 @@ namespace serif::eos { q.abar = in.composition.getMeanParticleMass(); // Mean atomic mass in g q.zbar = in.composition.getMeanAtomicNumber(); // Mean atomic number (dimensionless) + std::cout << "(EOS) abar: " << q.abar << ", zbar: " << q.zbar << std::endl; + helmholtz::HELMEOSOutput tempOutput; tempOutput = helmholtz::get_helm_EOS(q, *std::get>(m_reader.getTable())); diff --git a/src/eos/public/EOS.h b/src/eos/public/EOS.h index b2ecd22..5baf725 100644 --- a/src/eos/public/EOS.h +++ b/src/eos/public/EOS.h @@ -4,6 +4,7 @@ #include "helm.h" #include #include "composition.h" +#include namespace serif::eos { @@ -17,6 +18,13 @@ namespace serif::eos { 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). + friend std::ostream& operator<<(std::ostream& os, const EOSInput& input) { + os << ""; + return os; + } }; /** @@ -28,6 +36,12 @@ namespace serif::eos { * All values are in cgs units unless otherwise specified. */ struct EOSParameter { + explicit EOSParameter(std::string name_) + : total(), gas(), radiation(), + dDensity(), dTemperature(), + dMeanAtomicMassNumber(), + dMeanAtomicNumber(), + name(std::move(name_)) {} 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). @@ -38,6 +52,17 @@ namespace serif::eos { 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"). + + friend std::ostream& operator<<(std::ostream& os, const EOSParameter& param) { + os << std::setprecision(3) << ""; + return os; + } + + }; /** @@ -55,9 +80,9 @@ namespace serif::eos { 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. + EOSParameter pressure {"pressure"}; ///< Pressure output data, including total, gas, radiation, and derivatives. + EOSParameter energy {"energy"}; ///< Internal energy output data, including total, gas, radiation, and derivatives. + EOSParameter entropy {"entropy"}; ///< Entropy output data, including total, gas, radiation, and derivatives. /** * @brief Calculates the temperature susceptibility (chi_T). @@ -128,6 +153,17 @@ namespace serif::eos { * @return The EOSFormat enum value (currently only EOSFormat::HELM). */ EOSFormat EOSFormat() const; + + friend std::ostream& operator<<(std::ostream& os, const EOSOutput& output) { + os << "EOSOutput:\n" + << "\tElectron Fraction: " << output.electronFraction << "\n" + << "\tElectron Chemical Potential: " << output.electronChemicalPotential << "\n" + << "\tNeutron Excess Fraction: " << output.neutronExcessFraction << "\n\t" + << output.pressure << "\n\t" + << output.energy << "\n\t" + << output.entropy; + return os; + } }; /**