tboudreaux/libcomposition
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

fix(eos): fixed calculation of mean atomic number

Browse Source
This commit is contained in:
Emily Boudreaux
2025-06-17 08:12:41 -04:00
parent 4f47dc97a8
commit 36a3f832f7
2 changed files with 93 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files

100
src/composition/private/composition.cpp
View File

@@ -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<std::string, double> Composition::getMassFraction() const {
std::unordered_map<std::string, double> 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<std::string, double> Composition::getNumberFraction() const {
std::unordered_map<std::string, double> 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<std::string>& symbols, std::string method) const {
std::array<std::string, 2> methods = {"norm", "none"};
const std::array<std::string, 2> 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<std::string, 7> canonicalH = {
"H-1", "H-2", "H-3", "H-4", "H-5", "H-6", "H-7"
};
constexpr std::array<std::string, 8> 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;
}

26
src/composition/public/composition.h
View File

@@ -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 << "<CanonicalComposition: "
<< "X = " << composition.X << ", "
<< "Y = " << composition.Y << ", "
<< "Z = " << composition.Z << ">";
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.