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

feat(Added-ability-to-get-electron-abundance-and-fixed-some-conversion-bugs): Now Ye can be retrived directly from the composition object. Further a bug which prevented proper conversion to and from number or mass frac modes without messing up the numbers has been resolved

Browse Source
This commit is contained in:
Emily Boudreaux
2025-10-06 14:29:33 -04:00
parent 9d0d877e0f
commit 5fba54c1a2
3 changed files with 746 additions and 216 deletions
Download Patch File Download Diff File Expand all files Collapse all files

63
src/composition/include/fourdst/composition/composition.h
View File

@@ -88,6 +88,9 @@ namespace fourdst::composition {
double m_massFraction = 0.0; ///< The mass fraction of the species. Valid only if `m_massFracMode` is true.
double m_numberFraction = 0.0; ///< The number fraction (mole fraction) of the species. Valid only if `m_massFracMode` is false.
double m_relAbundance = 0.0; ///< The relative abundance, used internally for conversions. For mass fraction mode, this is X_i / A_i; for number fraction mode, it's n_i * A_i.
double m_molesPerMass = 0.0;
double m_cachedNumberFraction = 0.0; ///< Cached number fraction for conversions when in mass fraction mode.
bool m_initialized = false; ///< True if the composition entry has been initialized with a valid species.
@@ -137,13 +140,6 @@ namespace fourdst::composition {
*/
[[nodiscard]] double mass_fraction() const;
/**
* @brief Gets the mass fraction, converting from number fraction if necessary.
* @param meanMolarMass The mean molar mass of the entire composition, required for conversion.
* @return The mass fraction of the species.
*/
[[nodiscard]] double mass_fraction(double meanMolarMass) const;
/**
* @brief Gets the number fraction of the species.
* @pre The entry must be in number fraction mode.
@@ -154,10 +150,10 @@ namespace fourdst::composition {
/**
* @brief Gets the number fraction, converting from mass fraction if necessary.
* @param totalMoles The total moles per unit mass (specific number density) of the entire composition.
* @param totalMolesPerMass The total moles per unit mass (specific number density) of the entire composition.
* @return The number fraction of the species.
*/
[[nodiscard]] double number_fraction(double totalMoles) const;
[[nodiscard]] double number_fraction(double totalMolesPerMass) const;
/**
* @brief Gets the relative abundance of the species.
@@ -202,10 +198,10 @@ namespace fourdst::composition {
/**
* @brief Switches the mode to number fraction mode.
* @param totalMoles The total moles per unit mass (specific number density) of the composition.
* @param totalMolesPerMass The total moles per unit mass (specific number density) of the composition.
* @return True if the mode was successfully set, false otherwise.
*/
bool setNumberFracMode(double totalMoles);
bool setNumberFracMode(double totalMolesPerMass);
/**
* @brief Overloaded output stream operator for CompositionEntry.
@@ -257,8 +253,36 @@ namespace fourdst::composition {
*/
class Composition {
private:
fourdst::config::Config& m_config = fourdst::config::Config::getInstance();
fourdst::logging::LogManager& m_logManager = fourdst::logging::LogManager::getInstance();
struct CompositionCache {
std::optional<GlobalComposition> globalComp; ///< Cached global composition data.
std::optional<CanonicalComposition> canonicalComp; ///< Cached canonical composition data.
std::optional<std::vector<double>> massFractions; ///< Cached vector of mass fractions.
std::optional<std::vector<double>> numberFractions; ///< Cached vector of number fractions.
std::optional<std::vector<double>> molarAbundances; ///< Cached vector of molar abundances.
std::optional<std::vector<atomic::Species>> sortedSpecies; ///< Cached vector of sorted species (by mass).
std::optional<std::vector<std::string>> sortedSymbols; ///< Cached vector of sorted species (by mass).
std::optional<double> Ye; ///< Cached electron abundance.
void clear() {
globalComp = std::nullopt;
canonicalComp = std::nullopt;
massFractions = std::nullopt;
numberFractions = std::nullopt;
molarAbundances = std::nullopt;
sortedSymbols = std::nullopt;
sortedSpecies = std::nullopt;
Ye = std::nullopt;
}
[[nodiscard]] bool is_clear() const {
return !globalComp.has_value() && !canonicalComp.has_value() && !massFractions.has_value() &&
!numberFractions.has_value() && !molarAbundances.has_value() && !sortedSymbols.has_value() &&
!Ye.has_value() && !sortedSpecies.has_value();
}
};
private:
config::Config& m_config = config::Config::getInstance();
logging::LogManager& m_logManager = logging::LogManager::getInstance();
quill::Logger* m_logger = m_logManager.getLogger("log");
bool m_finalized = false; ///< True if the composition is finalized.
@@ -269,6 +293,9 @@ namespace fourdst::composition {
std::set<std::string> m_registeredSymbols; ///< The registered symbols.
std::unordered_map<std::string, CompositionEntry> m_compositions; ///< The compositions.
mutable CompositionCache m_cache; ///< Cache for computed properties to avoid redundant calculations.
/**
* @brief Checks if the given symbol is valid by checking against the global species database.
* @param symbol The symbol to check.
@@ -686,6 +713,14 @@ namespace fourdst::composition {
*/
[[nodiscard]] double getMeanAtomicNumber() const;
/**
* @brief Compute the electron abundance of the composition.
* @details Ye is defined as the sum over all species of (Z_i * X_i / A_i), where Z_i is the atomic number, X_i is the mass fraction, and A_i is the atomic mass of species i.
* @return Ye (electron abundance).
* @pre The composition must be finalized.
*/
[[nodiscard]] double getElectronAbundance() const;
/**
* @brief Creates a new Composition object containing a subset of species from this one.
* @param symbols The symbols to include in the subset.
@@ -711,7 +746,7 @@ namespace fourdst::composition {
* @return True if the isotope is in the composition, false otherwise.
* @throws exceptions::CompositionNotFinalizedError if the composition is not finalized.
*/
[[nodiscard]] bool contains(const fourdst::atomic::Species& isotope) const;
[[nodiscard]] bool contains(const atomic::Species& isotope) const;
/**
* @brief Sets the composition mode (mass fraction vs. number fraction).

522
src/composition/lib/composition.cpp
View File

@@ -2,7 +2,7 @@
//
// Copyright (C) 2025 -- The 4D-STAR Collaboration
// File Author: Emily Boudreaux
// Last Modified: March 26, 2025
// Last Modified: October 6, 2025
//
// 4DSSE is free software; you can use it and/or modify
// it under the terms and restrictions the GNU General Library Public
@@ -33,11 +33,17 @@
#include "fourdst/composition/atomicSpecies.h"
#include "fourdst/composition/species.h"
#include "fourdst/composition/composition.h"
#include <numeric>
#include "fourdst/composition/exceptions/exceptions_composition.h"
namespace {
template<typename A, typename B>
std::vector<A> sortVectorBy(std::vector<A> toSort, const std::vector<B>& by) {
std::vector<A> sortVectorBy(
std::vector<A> toSort,
const std::vector<B>& by
) {
std::vector<std::size_t> indices(by.size());
for (size_t i = 0; i < indices.size(); i++) {
indices[i] = i;
@@ -63,20 +69,19 @@ namespace fourdst::composition {
CompositionEntry::CompositionEntry() :
m_symbol("H-1"),
m_isotope(fourdst::atomic::species.at("H-1")),
m_initialized(false) {} // Note: Default entry is uninitialized, must be explicitly set.
m_initialized(false) {}
CompositionEntry::CompositionEntry(const std::string& symbol, const bool massFracMode) : m_symbol(symbol), m_isotope(fourdst::atomic::species.at(symbol)), m_massFracMode(massFracMode) {
CompositionEntry::CompositionEntry(
const std::string& symbol,
const bool massFracMode
) :
m_symbol(symbol),
m_isotope(atomic::species.at(symbol)),
m_massFracMode(massFracMode) {
setSpecies(symbol);
}
CompositionEntry::CompositionEntry(const CompositionEntry& entry) :
m_symbol(entry.m_symbol),
m_isotope(entry.m_isotope),
m_massFracMode(entry.m_massFracMode),
m_massFraction(entry.m_massFraction),
m_numberFraction(entry.m_numberFraction),
m_relAbundance(entry.m_relAbundance),
m_initialized(entry.m_initialized) {}
CompositionEntry::CompositionEntry(const CompositionEntry& entry) = default;
void CompositionEntry::setSpecies(const std::string& symbol) {
if (m_initialized) {
@@ -86,7 +91,7 @@ namespace fourdst::composition {
throw exceptions::InvalidSpeciesSymbolError("Invalid symbol.");
}
m_symbol = symbol;
m_isotope = fourdst::atomic::species.at(symbol);
m_isotope = atomic::species.at(symbol);
m_initialized = true;
}
@@ -98,74 +103,80 @@ namespace fourdst::composition {
if (!m_massFracMode) {
throw exceptions::CompositionModeError("Composition entry is in number fraction mode.");
}
return m_massFraction;
// X_i = (moles_i / mass_total) * (mass_i / moles_i) = m_molesPerMass * A_i
return m_molesPerMass * m_isotope.mass();
}
double CompositionEntry::mass_fraction(const double meanMolarMass) const {
if (m_massFracMode) {
return m_massFraction;
}
// Convert from number fraction to mass fraction using: X_i = n_i * A_i / <A>
// where m_relAbundance is n_i * A_i and meanMolarMass is <A>.
return m_relAbundance / meanMolarMass;
}
double CompositionEntry::number_fraction() const {
if (m_massFracMode) {
throw exceptions::CompositionModeError("Composition entry is in mass fraction mode.");
}
return m_numberFraction;
// In number fraction mode, the value is cached during the mode switch.
return m_cachedNumberFraction;
}
double CompositionEntry::number_fraction(const double totalMoles) const {
if (m_massFracMode) {
// Convert from mass fraction to number fraction using: n_i = (X_i / A_i) / sum(X_j / A_j)
// where m_relAbundance is X_i / A_i and totalMoles is sum(X_j / A_j).
return m_relAbundance / totalMoles;
}
return m_numberFraction;
double CompositionEntry::number_fraction(
const double totalMolesPerMass
) const {
// n_i = (moles_i / mass_total) / (moles_total / mass_total)
if (totalMolesPerMass == 0.0) return 0.0;
return m_molesPerMass / totalMolesPerMass;
}
double CompositionEntry::rel_abundance() const {
return m_relAbundance;
return m_molesPerMass;
}
fourdst::atomic::Species CompositionEntry::isotope() const {
atomic::Species CompositionEntry::isotope() const {
return m_isotope;
}
void CompositionEntry::setMassFraction(const double mass_fraction) {
void CompositionEntry::setMassFraction(
const double mass_fraction
) {
if (!m_massFracMode) {
throw exceptions::CompositionModeError("Composition entry is in number fraction mode.");
}
m_massFraction = mass_fraction;
m_relAbundance = m_massFraction / m_isotope.mass();
// Set the invariant from the given mass fraction
if (m_isotope.mass() == 0.0) {
m_molesPerMass = 0.0;
} else {
m_molesPerMass = mass_fraction / m_isotope.mass();
}
}
void CompositionEntry::setNumberFraction(const double number_fraction) {
void CompositionEntry::setNumberFraction(
const double number_fraction
) {
if (m_massFracMode) {
throw exceptions::CompositionModeError("Composition entry is in mass fraction mode.");
}
m_numberFraction = number_fraction;
m_relAbundance = m_numberFraction * m_isotope.mass();
// In number fraction mode, we only cache the value. The invariant
// m_molesPerMass cannot be calculated until finalize() provides global context.
m_cachedNumberFraction = number_fraction;
}
bool CompositionEntry::setMassFracMode(const double meanParticleMass) {
bool CompositionEntry::setMassFracMode(
[[maybe_unused]] const double meanMolarMass
) {
if (m_massFracMode) {
return false;
}
m_massFracMode = true;
m_massFraction = m_relAbundance / meanParticleMass;
// The invariant m_molesPerMass does not change when switching mode.
// The cached number fraction is now stale, but that's okay.
return true;
}
bool CompositionEntry::setNumberFracMode(const double specificNumberDensity) {
bool CompositionEntry::setNumberFracMode(
const double totalMolesPerMass
) {
if (!m_massFracMode) {
return false;
}
m_massFracMode = false;
m_numberFraction = m_relAbundance / specificNumberDensity;
// Calculate and cache the number fraction for the new mode.
m_cachedNumberFraction = number_fraction(totalMolesPerMass);
return true;
}
@@ -173,19 +184,27 @@ namespace fourdst::composition {
return m_massFracMode;
}
Composition::Composition(const std::vector<std::string>& symbols) {
Composition::Composition(
const std::vector<std::string>& symbols
) {
for (const auto& symbol : symbols) {
registerSymbol(symbol);
}
}
Composition::Composition(const std::set<std::string>& symbols) {
Composition::Composition(
const std::set<std::string>& symbols
) {
for (const auto& symbol : symbols) {
registerSymbol(symbol);
}
}
Composition::Composition(const std::vector<std::string>& symbols, const std::vector<double>& fractions, const bool massFracMode) : m_massFracMode(massFracMode) {
Composition::Composition(
const std::vector<std::string>& symbols,
const std::vector<double>& fractions,
const bool massFracMode
) : m_massFracMode(massFracMode) {
if (symbols.size() != fractions.size()) {
LOG_CRITICAL(m_logger, "The number of symbols and fractions must be equal (got {} symbols and {} fractions).", symbols.size(), fractions.size());
throw exceptions::InvalidCompositionError("The number of symbols and fractions must be equal. Got " + std::to_string(symbols.size()) + " symbols and " + std::to_string(fractions.size()) + " fractions.");
@@ -194,7 +213,7 @@ namespace fourdst::composition {
validateComposition(fractions);
for (const auto &symbol : symbols) {
registerSymbol(symbol);
registerSymbol(symbol, m_massFracMode);
}
for (size_t i = 0; i < symbols.size(); ++i) {
@@ -224,25 +243,25 @@ namespace fourdst::composition {
m_massFracMode = other.m_massFracMode;
m_registeredSymbols = other.m_registeredSymbols;
m_compositions = other.m_compositions;
// note: m_config remains bound to the same singleton, so we skip it
}
return *this;
}
void Composition::registerSymbol(const std::string& symbol, bool massFracMode) {
void Composition::registerSymbol(
const std::string& symbol,
const bool massFracMode
) {
if (!isValidSymbol(symbol)) {
LOG_ERROR(m_logger, "Invalid symbol: {}", symbol);
throw exceptions::InvalidSymbolError("Invalid symbol: " + symbol);
}
// If no symbols have been registered allow mode to be set
if (m_registeredSymbols.empty()) {
m_massFracMode = massFracMode;
} else {
if (m_massFracMode != massFracMode) {
LOG_ERROR(m_logger, "Composition is in mass fraction mode. Cannot register symbol ({}) in number fraction mode.", symbol);
throw exceptions::CompositionModeError("Composition is in mass fraction mode. Cannot register symbol (" + symbol + ") in number fraction mode.");
LOG_ERROR(m_logger, "Composition is in {} fraction mode. Cannot register symbol ({}) in {} fraction mode.", m_massFracMode ? "mass" : "number", symbol, massFracMode ? "mass" : "number");
throw exceptions::CompositionModeError("Composition mode mismatch.");
}
}
@@ -252,22 +271,31 @@ namespace fourdst::composition {
}
m_registeredSymbols.insert(symbol);
const CompositionEntry entry(symbol, m_massFracMode);
m_compositions[symbol] = entry;
m_compositions[symbol] = CompositionEntry(symbol, m_massFracMode);
m_finalized = false;
LOG_TRACE_L3(m_logger, "Registered symbol: {}", symbol);
}
void Composition::registerSymbol(const std::vector<std::string>& symbols, bool massFracMode) {
void Composition::registerSymbol(
const std::vector<std::string>& symbols,
const bool massFracMode
) {
for (const auto& symbol : symbols) {
registerSymbol(symbol, massFracMode);
}
}
void Composition::registerSpecies(const fourdst::atomic::Species &species, bool massFracMode) {
registerSymbol(std::string(species.name()));
void Composition::registerSpecies(
const atomic::Species &species,
const bool massFracMode
) {
registerSymbol(std::string(species.name()), massFracMode);
}
void Composition::registerSpecies(const std::vector<fourdst::atomic::Species> &species, bool massFracMode) {
void Composition::registerSpecies(
const std::vector<atomic::Species> &species,
const bool massFracMode
) {
for (const auto& s : species) {
registerSpecies(s, massFracMode);
}
@@ -277,14 +305,20 @@ namespace fourdst::composition {
return m_registeredSymbols;
}
std::set<fourdst::atomic::Species> Composition::getRegisteredSpecies() const {
std::set<fourdst::atomic::Species> result;
std::set<atomic::Species> Composition::getRegisteredSpecies() const {
std::set<atomic::Species> result;
for (const auto& entry : m_compositions | std::views::values) {
result.insert(entry.isotope());
}
return result;
}
bool Composition::isValidSymbol(
const std::string& symbol
) {
return atomic::species.contains(symbol);
}
void Composition::validateComposition(const std::vector<double>& fractions) const {
if (!isValidComposition(fractions)) {
LOG_ERROR(m_logger, "Invalid composition.");
@@ -293,51 +327,37 @@ namespace fourdst::composition {
}
bool Composition::isValidComposition(const std::vector<double>& fractions) const {
double sum = 0.0;
for (const auto& fraction : fractions) {
sum += fraction;
}
double sum = std::accumulate(fractions.begin(), fractions.end(), 0.0);
if (sum < 0.999999 || sum > 1.000001) {
LOG_ERROR(m_logger, "The sum of fractions must be equal to 1 (expected 1, got {}).", sum);
return false;
}
return true;
}
bool Composition::isValidSymbol(const std::string& symbol) {
return fourdst::atomic::species.contains(symbol);
}
double Composition::setMassFraction(const std::string& symbol, const double& mass_fraction) {
if (!m_registeredSymbols.contains(symbol)) {
LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol);
throw exceptions::UnregisteredSymbolError("Symbol (" + symbol + ") is not registered.");
}
if (!m_massFracMode) {
LOG_ERROR(m_logger, "Composition is in number fraction mode.");
throw exceptions::CompositionModeError("Composition is in number fraction mode.");
}
if (mass_fraction < 0.0 || mass_fraction > 1.0) {
LOG_ERROR(m_logger, "Mass fraction must be between 0 and 1 for symbol {}. Currently it is {}.", symbol, mass_fraction);
throw exceptions::InvalidCompositionError("Mass fraction must be between 0 and 1 for symbol " + symbol + ". Currently it is " + std::to_string(mass_fraction) + ".");
throw exceptions::InvalidCompositionError("Mass fraction must be between 0 and 1.");
}
m_finalized = false;
const double old_mass_fraction = m_compositions.at(symbol).mass_fraction();
m_compositions.at(symbol).setMassFraction(mass_fraction);
return old_mass_fraction;
}
std::vector<double> Composition::setMassFraction(const std::vector<std::string>& symbols, const std::vector<double>& mass_fractions) {
if (symbols.size() != mass_fractions.size()) {
LOG_ERROR(m_logger, "The number of symbols and mass fractions must be equal (currently {} symbols and {} mass fractions).", symbols.size(), mass_fractions.size());
throw exceptions::InvalidCompositionError("The number of symbols and mass fractions must be equal (currently " + std::to_string(symbols.size()) + " symbols and " + std::to_string(mass_fractions.size()) + " mass fractions).");
throw exceptions::InvalidCompositionError("The number of symbols and mass fractions must be equal.");
}
std::vector<double> old_mass_fractions;
old_mass_fractions.reserve(symbols.size());
for (size_t i = 0; i < symbols.size(); ++i) {
@@ -346,49 +366,35 @@ namespace fourdst::composition {
return old_mass_fractions;
}
double Composition::setMassFraction(const fourdst::atomic::Species &species, const double &mass_fraction) {
return setMassFraction(std::string(species.name()), mass_fraction);
}
std::vector<double> Composition::setMassFraction(const std::vector<fourdst::atomic::Species> &species,
const std::vector<double> &mass_fractions) {
std::vector<double> old_mass_fractions;
old_mass_fractions.reserve(species.size());
for (const auto& spec : species) {
old_mass_fractions.push_back(setMassFraction(spec, mass_fractions[&spec - &species[0]]));
}
return old_mass_fractions;
}
double Composition::setNumberFraction(const std::string& symbol, const double& number_fraction) {
double Composition::setNumberFraction(
const std::string& symbol,
const double& number_fraction
) {
if (!m_registeredSymbols.contains(symbol)) {
LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol);
throw exceptions::UnregisteredSymbolError("Symbol (" + symbol + ") is not registered.");
}
if (m_massFracMode) {
LOG_ERROR(m_logger, "Composition is in mass fraction mode, should be in number fraction mode to call setNumberFraction. Hint: The mode can be switched by first finalizing and then calling setCompositionMode(false).");
throw exceptions::CompositionModeError("Composition is in mass fraction mode, should be in number fraction mode to call setNumberFraction. Hint: The mode can be switched by first finalizing and then calling setCompositionMode(false).");
LOG_ERROR(m_logger, "Composition is in mass fraction mode.");
throw exceptions::CompositionModeError("Composition is in mass fraction mode.");
}
if (number_fraction < 0.0 || number_fraction > 1.0) {
LOG_ERROR(m_logger, "Number fraction must be between 0 and 1 for symbol {}. Currently it is {}.", symbol, number_fraction);
throw exceptions::InvalidCompositionError("Number fraction must be between 0 and 1 for symbol " + symbol + ". Currently it is " + std::to_string(number_fraction) + ".");
throw exceptions::InvalidCompositionError("Number fraction must be between 0 and 1.");
}
m_finalized = false;
const double old_number_fraction = m_compositions.at(symbol).number_fraction();
m_compositions.at(symbol).setNumberFraction(number_fraction);
return old_number_fraction;
}
std::vector<double> Composition::setNumberFraction(const std::vector<std::string>& symbols, const std::vector<double>& number_fractions) {
std::vector<double> Composition::setNumberFraction(
const std::vector<std::string>& symbols,
const std::vector<double>& number_fractions
) {
if (symbols.size() != number_fractions.size()) {
LOG_ERROR(m_logger, "The number of symbols and number fractions must be equal. (Currently {} symbols and {} number fractions).", symbols.size(), number_fractions.size());
throw exceptions::InvalidCompositionError("The number of symbols and number fractions must be equal. (Currently " + std::to_string(symbols.size()) + " symbols and " + std::to_string(number_fractions.size()) + " number fractions).");
throw exceptions::InvalidCompositionError("The number of symbols and number fractions must be equal.");
}
std::vector<double> old_number_fractions;
old_number_fractions.reserve(symbols.size());
for (size_t i = 0; i < symbols.size(); ++i) {
@@ -397,101 +403,125 @@ namespace fourdst::composition {
return old_number_fractions;
}
double Composition::setNumberFraction(const fourdst::atomic::Species &species, const double &number_fraction) {
double Composition::setMassFraction(
const atomic::Species &species,
const double &mass_fraction
) {
return setMassFraction(std::string(species.name()), mass_fraction);
}
std::vector<double> Composition::setMassFraction(
const std::vector<atomic::Species> &species,
const std::vector<double> &mass_fractions
) {
std::vector<std::string> symbols;
symbols.reserve(species.size());
for(const auto& s : species) symbols.push_back(std::string(s.name()));
return setMassFraction(symbols, mass_fractions);
}
double Composition::setNumberFraction(
const atomic::Species &species,
const double &number_fraction
) {
return setNumberFraction(std::string(species.name()), number_fraction);
}
std::vector<double> Composition::setNumberFraction(const std::vector<fourdst::atomic::Species> &species,
const std::vector<double> &number_fractions) {
std::vector<double> old_number_fractions;
old_number_fractions.reserve(species.size());
for (const auto& spec : species) {
old_number_fractions.push_back(setNumberFraction(spec, number_fractions[&spec - &species[0]]));
}
return old_number_fractions;
std::vector<double> Composition::setNumberFraction(
const std::vector<atomic::Species> &species,
const std::vector<double> &number_fractions
) {
std::vector<std::string> symbols;
symbols.reserve(species.size());
for(const auto& s : species) symbols.push_back(std::string(s.name()));
return setNumberFraction(symbols, number_fractions);
}
bool Composition::finalize(const bool norm) {
bool finalized = false;
if (m_massFracMode) {
finalized = finalizeMassFracMode(norm);
} else {
finalized = finalizeNumberFracMode(norm);
}
if (finalized) {
m_finalized = true;
}
return finalized;
m_specificNumberDensity = 0.0;
m_meanParticleMass = 0.0;
m_finalized = m_massFracMode ? finalizeMassFracMode(norm) : finalizeNumberFracMode(norm);
m_cache.clear();
return m_finalized;
}
bool Composition::finalizeMassFracMode(bool norm) {
bool Composition::finalizeMassFracMode(const bool norm) {
std::vector<double> mass_fractions;
mass_fractions.reserve(m_compositions.size());
for (const auto &entry: m_compositions | std::views::values) {
mass_fractions.push_back(entry.mass_fraction());
}
if (norm) {
double sum = 0.0;
for (const auto& mass_fraction : mass_fractions) {
sum += mass_fraction;
}
for (double & mass_fraction : mass_fractions) {
mass_fraction /= sum;
}
double sum = std::accumulate(mass_fractions.begin(), mass_fractions.end(), 0.0);
if (norm && sum > 0) {
for (auto& [symbol, entry] : m_compositions) {
setMassFraction(symbol, entry.mass_fraction() / sum);
}
// Recalculate fractions vector after normalization for validation
mass_fractions.clear();
for (const auto &entry: m_compositions | std::views::values) {
mass_fractions.push_back(entry.mass_fraction());
}
}
try {
validateComposition(mass_fractions);
} catch ([[maybe_unused]] const exceptions::InvalidCompositionError& e) {
double massSum = 0.0;
for (const auto &entry: m_compositions | std::views::values) {
massSum += entry.mass_fraction();
}
LOG_ERROR(m_logger, "Composition is invalid (Total mass {}).", massSum);
m_finalized = false;
LOG_ERROR(m_logger, "Composition is invalid after mass frac finalization (Total mass {}).", sum);
return false;
}
// After validation, calculate the specific number density (total moles per unit mass).
for (const auto &entry: m_compositions | std::views::values) {
m_specificNumberDensity += entry.rel_abundance();
m_specificNumberDensity += entry.rel_abundance(); // rel_abundance is now consistently moles/mass
}
if (m_specificNumberDensity > 0) {
m_meanParticleMass = 1.0 / m_specificNumberDensity;
}
m_meanParticleMass = 1.0/m_specificNumberDensity;
return true;
}
bool Composition::finalizeNumberFracMode(bool norm) {
bool Composition::finalizeNumberFracMode(const bool norm) {
std::vector<double> number_fractions;
number_fractions.reserve(m_compositions.size());
for (const auto &entry: m_compositions | std::views::values) {
number_fractions.push_back(entry.number_fraction());
}
if (norm) {
double sum = 0.0;
for (const auto& number_fraction : number_fractions) {
sum += number_fraction;
}
double sum = std::accumulate(number_fractions.begin(), number_fractions.end(), 0.0);
if (norm && sum > 0) {
for (auto& [symbol, entry] : m_compositions) {
setNumberFraction(symbol, entry.number_fraction() / sum);
}
// Recalculate fractions vector after normalization for validation
number_fractions.clear();
for (const auto &entry: m_compositions | std::views::values) {
number_fractions.push_back(entry.number_fraction());
}
}
try {
validateComposition(number_fractions);
} catch ([[maybe_unused]] const std::runtime_error& e) {
double numberSum = 0.0;
for (const auto &entry: m_compositions | std::views::values) {
numberSum += entry.number_fraction();
}
LOG_ERROR(m_logger, "Composition is invalid (Total number {}).", numberSum);
m_finalized = false;
} catch ([[maybe_unused]] const exceptions::InvalidCompositionError& e) {
LOG_ERROR(m_logger, "Composition is invalid after number frac finalization (Total number frac {}).", sum);
return false;
}
// After validation, calculate the mean particle mass.
// Calculate mean particle mass <A> = sum(n_i * A_i)
for (const auto &entry: m_compositions | std::views::values) {
m_meanParticleMass += entry.rel_abundance();
m_meanParticleMass += entry.number_fraction() * entry.isotope().mass();
}
for (auto &entry: m_compositions | std::views::values) {
const double X_i = (m_meanParticleMass > 0) ? (entry.number_fraction() * entry.isotope().mass() / m_meanParticleMass) : 0.0;
entry.m_massFracMode = true;
entry.setMassFraction(X_i);
entry.m_massFracMode = false;
}
if (m_meanParticleMass > 0) {
m_specificNumberDensity = 1.0 / m_meanParticleMass;
}
m_specificNumberDensity = 1.0/m_meanParticleMass;
return true;
}
@@ -547,12 +577,14 @@ namespace fourdst::composition {
}
if (m_massFracMode) {
return m_compositions.at(symbol).mass_fraction();
} else {
return m_compositions.at(symbol).mass_fraction(m_meanParticleMass);
}
return m_compositions.at(symbol).mass_fraction();
}
double Composition::getMassFraction(const fourdst::atomic::Species &species) const {
double Composition::getMassFraction(
const atomic::Species &species
) const {
return getMassFraction(std::string(species.name()));
}
@@ -565,7 +597,9 @@ namespace fourdst::composition {
}
double Composition::getNumberFraction(const std::string& symbol) const {
double Composition::getNumberFraction(
const std::string& symbol
) const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize().");
@@ -576,12 +610,13 @@ namespace fourdst::composition {
}
if (!m_massFracMode) {
return m_compositions.at(symbol).number_fraction();
} else {
return m_compositions.at(symbol).number_fraction(m_specificNumberDensity);
}
return m_compositions.at(symbol).number_fraction(m_specificNumberDensity);
}
double Composition::getNumberFraction(const fourdst::atomic::Species &species) const {
double Composition::getNumberFraction(
const atomic::Species &species
) const {
return getNumberFraction(std::string(species.name()));
}
@@ -593,7 +628,9 @@ namespace fourdst::composition {
return number_fractions;
}
double Composition::getMolarAbundance(const std::string &symbol) const {
double Composition::getMolarAbundance(
const std::string &symbol
) const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize().");
@@ -606,11 +643,15 @@ namespace fourdst::composition {
}
double Composition::getMolarAbundance(const fourdst::atomic::Species &species) const {
double Composition::getMolarAbundance(
const atomic::Species &species
) const {
return getMolarAbundance(std::string(species.name()));
}
std::pair<CompositionEntry, GlobalComposition> Composition::getComposition(const std::string& symbol) const {
std::pair<CompositionEntry, GlobalComposition> Composition::getComposition(
const std::string& symbol
) const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize().");
@@ -623,7 +664,8 @@ namespace fourdst::composition {
}
std::pair<CompositionEntry, GlobalComposition> Composition::getComposition(
const fourdst::atomic::Species &species) const {
const atomic::Species &species
) const {
return getComposition(std::string(species.name()));
}
@@ -651,21 +693,39 @@ namespace fourdst::composition {
double zSum = 0.0;
// Loop through all registered species in the composition.
for (const auto &val: m_compositions | std::views::values) {
// Sum of (X_i * Z_i / A_i)
zSum += (val.mass_fraction() * val.m_isotope.z())/val.m_isotope.a();
}
// Calculate mean atomic number <Z> = <A> * sum(X_i * Z_i / A_i)
// <Z> = <A> * sum(X_i * Z_i / A_i)
const double mean_A = m_meanParticleMass * zSum;
return mean_A;
}
Composition Composition::subset(const std::vector<std::string>& symbols, const std::string& method) const {
const std::array<std::string, 2> methods = {"norm", "none"};
double Composition::getElectronAbundance() const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition must be finalized before getting the electron abundance. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition not finalized. Cannot retrieve electron abundance. Hint: Consider running .finalize().");
}
if (std::ranges::find(methods, method) == methods.end()) {
if (m_cache.Ye.has_value()) {
return m_cache.Ye.value();
}
double Ye = 0.0;
for (const auto &val: m_compositions | std::views::values) {
Ye += (val.mass_fraction() * val.m_isotope.z())/val.m_isotope.a();
}
m_cache.Ye = Ye;
return Ye;
}
Composition Composition::subset(
const std::vector<std::string>& symbols,
const std::string& method
) const {
if (const std::array<std::string, 2> methods = {"norm", "none"}; std::ranges::find(methods, method) == methods.end()) {
const std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'.";
LOG_ERROR(m_logger, "Invalid method: {}. Valid methods are norm and none.", method);
throw exceptions::InvalidMixingMode(errorMessage);
@@ -676,14 +736,12 @@ namespace fourdst::composition {
if (!m_compositions.contains(symbol)) {
LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol);
throw exceptions::UnregisteredSymbolError("Symbol " + symbol + " is not in the composition.");
} else {
subsetComposition.registerSymbol(symbol);
}
subsetComposition.registerSymbol(symbol);
subsetComposition.setMassFraction(symbol, m_compositions.at(symbol).mass_fraction());
}
if (method == "norm") {
const bool isNorm = subsetComposition.finalize(true);
if (!isNorm) {
if (const bool isNorm = subsetComposition.finalize(true); !isNorm) {
LOG_ERROR(m_logger, "Subset composition is invalid. (Unable to finalize with normalization).");
throw exceptions::FailedToFinalizeCompositionError("Subset composition is invalid. (Unable to finalize with normalization).");
}
@@ -691,13 +749,15 @@ namespace fourdst::composition {
return subsetComposition;
}
void Composition::setCompositionMode(const bool massFracMode) {
void Composition::setCompositionMode(
const bool massFracMode
) {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized. Mode cannot be set unless composition is finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Mode cannot be set unless composition is finalized. Hint: Consider running .finalize().");
}
bool okay = true;
bool okay;
for (auto &entry: m_compositions | std::views::values) {
if (massFracMode) {
okay = entry.setMassFracMode(m_meanParticleMass);
@@ -712,11 +772,16 @@ namespace fourdst::composition {
m_massFracMode = massFracMode;
}
CanonicalComposition Composition::getCanonicalComposition(bool harsh) const {
CanonicalComposition Composition::getCanonicalComposition(
const bool harsh
) const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize().");
}
if (m_cache.canonicalComp.has_value()) {
return m_cache.canonicalComp.value(); // Short circuit if we have cached the canonical composition
}
CanonicalComposition canonicalComposition;
const std::array<std::string, 7> canonicalH = {
"H-1", "H-2", "H-3", "H-4", "H-5", "H-6", "H-7"
@@ -737,6 +802,7 @@ namespace fourdst::composition {
for (const auto& symbol : getRegisteredSymbols()) {
const bool isHSymbol = std::ranges::find(canonicalH, symbol) != std::end(canonicalH);
// ReSharper disable once CppTooWideScopeInitStatement
const bool isHeSymbol = std::ranges::find(canonicalHe, symbol) != std::end(canonicalHe);
if (isHSymbol || isHeSymbol) {
@@ -746,6 +812,7 @@ namespace fourdst::composition {
canonicalComposition.Z += getMassFraction(symbol);
}
// ReSharper disable once CppTooWideScopeInitStatement
const double Z = 1.0 - (canonicalComposition.X + canonicalComposition.Y);
if (std::abs(Z - canonicalComposition.Z) > 1e-6) {
if (!harsh) {
@@ -756,6 +823,7 @@ namespace fourdst::composition {
throw std::runtime_error("Validation composition Z (X-Y = " + std::to_string(Z) + ") is different than canonical composition Z (" + std::to_string(canonicalComposition.Z) + ") (∑a_i where a_i != H/He).");
}
}
m_cache.canonicalComp = canonicalComposition;
return canonicalComposition;
}
@@ -764,6 +832,9 @@ namespace fourdst::composition {
LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize().");
}
if (m_cache.massFractions.has_value()) {
return m_cache.massFractions.value(); // Short circuit if we have cached the mass fractions
}
std::vector<double> massFractionVector;
std::vector<double> speciesMass;
@@ -776,7 +847,9 @@ namespace fourdst::composition {
speciesMass.push_back(entry.isotope().mass());
}
return sortVectorBy(massFractionVector, speciesMass);
std::vector<double> massFractions = sortVectorBy(massFractionVector, speciesMass);
m_cache.massFractions = massFractions; // Cache the result
return massFractions;
}
@@ -785,6 +858,9 @@ namespace fourdst::composition {
LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize().");
}
if (m_cache.numberFractions.has_value()) {
return m_cache.numberFractions.value(); // Short circuit if we have cached the number fractions
}
std::vector<double> numberFractionVector;
std::vector<double> speciesMass;
@@ -797,7 +873,9 @@ namespace fourdst::composition {
speciesMass.push_back(entry.isotope().mass());
}
return sortVectorBy(numberFractionVector, speciesMass);
std::vector<double> numberFractions = sortVectorBy(numberFractionVector, speciesMass);
m_cache.numberFractions = numberFractions; // Cache the result
return numberFractions;
}
std::vector<double> Composition::getMolarAbundanceVector() const {
@@ -805,6 +883,9 @@ namespace fourdst::composition {
LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize().");
}
if (m_cache.molarAbundances.has_value()) {
return m_cache.molarAbundances.value(); // Short circuit if we have cached the molar abundances
}
std::vector<double> molarAbundanceVector;
std::vector<double> speciesMass;
@@ -817,10 +898,15 @@ namespace fourdst::composition {
speciesMass.push_back(entry.isotope().mass());
}
return sortVectorBy(molarAbundanceVector, speciesMass);
std::vector<double> molarAbundances = sortVectorBy(molarAbundanceVector, speciesMass);
m_cache.molarAbundances = molarAbundances; // Cache the result
return molarAbundances;
}
size_t Composition::getSpeciesIndex(const std::string &symbol) const {
size_t Composition::getSpeciesIndex(
const std::string &symbol
) const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize().");
@@ -829,6 +915,16 @@ namespace fourdst::composition {
LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol);
throw exceptions::UnregisteredSymbolError("Symbol " + symbol + " is not in the composition.");
}
if (m_cache.sortedSymbols.has_value()) {
return std::distance(
m_cache.sortedSymbols->begin(),
std::ranges::find(
m_cache.sortedSymbols.value().begin(),
m_cache.sortedSymbols.value().end(),
symbol
)
);
}
std::vector<std::string> symbols;
std::vector<double> speciesMass;
@@ -842,10 +938,13 @@ namespace fourdst::composition {
}
std::vector<std::string> sortedSymbols = sortVectorBy(symbols, speciesMass);
m_cache.sortedSymbols = sortedSymbols;
return std::distance(sortedSymbols.begin(), std::ranges::find(sortedSymbols, symbol));
}
size_t Composition::getSpeciesIndex(const atomic::Species &species) const {
size_t Composition::getSpeciesIndex(
const atomic::Species &species
) const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize().");
@@ -854,6 +953,16 @@ namespace fourdst::composition {
LOG_ERROR(m_logger, "Species {} is not in the composition.", species.name());
throw exceptions::UnregisteredSymbolError("Species " + std::string(species.name()) + " is not in the composition.");
}
if (m_cache.sortedSpecies.has_value()) {
return std::distance(
m_cache.sortedSpecies->begin(),
std::ranges::find(
m_cache.sortedSpecies.value().begin(),
m_cache.sortedSpecies.value().end(),
species
)
);
}
std::vector<atomic::Species> speciesVector;
std::vector<double> speciesMass;
@@ -867,10 +976,13 @@ namespace fourdst::composition {
}
std::vector<atomic::Species> sortedSpecies = sortVectorBy(speciesVector, speciesMass);
m_cache.sortedSpecies = sortedSpecies;
return std::distance(sortedSpecies.begin(), std::ranges::find(sortedSpecies, species));
}
atomic::Species Composition::getSpeciesAtIndex(size_t index) const {
atomic::Species Composition::getSpeciesAtIndex(
size_t index
) const {
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize().");
@@ -879,9 +991,8 @@ namespace fourdst::composition {
LOG_ERROR(m_logger, "Index {} is out of bounds for composition of size {}.", index, m_compositions.size());
throw std::out_of_range("Index " + std::to_string(index) + " is out of bounds for composition of size " + std::to_string(m_compositions.size()) + ".");
}
if (index < 0) {
LOG_ERROR(m_logger, "Index {} is negative. Cannot get species at negative index.", index);
throw std::out_of_range("Index " + std::to_string(index) + " is negative. Cannot get species at negative index.");
if (m_cache.sortedSpecies.has_value()) {
return m_cache.sortedSpecies.value().at(index);
}
std::vector<atomic::Species> speciesVector;
@@ -899,18 +1010,21 @@ namespace fourdst::composition {
return sortedSymbols.at(index);
}
bool Composition::hasSymbol(const std::string& symbol) const {
bool Composition::hasSymbol(
const std::string& symbol
) const {
return m_compositions.contains(symbol);
}
bool Composition::contains(const fourdst::atomic::Species &isotope) const {
bool Composition::contains(
const atomic::Species &isotope
) const {
// Check if the isotope's symbol is in the composition
if (!m_finalized) {
LOG_ERROR(m_logger, "Composition has not been finalized. Hint: Consider running .finalize().");
throw exceptions::CompositionNotFinalizedError("Composition has not been finalized. Hint: Consider running .finalize().");
}
const auto symbol = static_cast<std::string>(isotope.name());
if (m_compositions.contains(symbol)) {
if (const auto symbol = static_cast<std::string>(isotope.name()); m_compositions.contains(symbol)) {
return true;
}
return false;
@@ -918,23 +1032,34 @@ namespace fourdst::composition {
/// OVERLOADS
Composition Composition::operator+(const Composition& other) const {
Composition Composition::operator+(
const Composition& other
) const {
return mix(other, 0.5);
}
std::ostream& operator<<(std::ostream& os, const GlobalComposition& comp) {
std::ostream& operator<<(
std::ostream& os,
const GlobalComposition& comp
) {
os << "Global Composition: \n";
os << "\tSpecific Number Density: " << comp.specificNumberDensity << "\n";
os << "\tMean Particle Mass: " << comp.meanParticleMass << "\n";
return os;
}
std::ostream& operator<<(std::ostream& os, const CompositionEntry& entry) {
std::ostream& operator<<(
std::ostream& os,
const CompositionEntry& entry
) {
os << "<" << entry.m_symbol << " : m_frac = " << entry.mass_fraction() << ">";
return os;
}
std::ostream& operator<<(std::ostream& os, const Composition& composition) {
std::ostream& operator<<(
std::ostream& os,
const Composition& composition
) {
os << "Composition(finalized: " << (composition.m_finalized ? "true" : "false") << ", " ;
size_t count = 0;
for (const auto &entry: composition.m_compositions | std::views::values) {
@@ -942,6 +1067,7 @@ namespace fourdst::composition {
if (count < composition.m_compositions.size() - 1) {
os << ", ";
}
count++;
}
os << ")";
return os;

377
tests/composition/compositionTest.cpp
View File

@@ -116,10 +116,10 @@ TEST_F(compositionTest, registerSymbol) {
EXPECT_THROW(comp.registerSymbol("He-21"), fourdst::composition::exceptions::InvalidSymbolError);
std::set<std::string> registeredSymbols = comp.getRegisteredSymbols();
EXPECT_TRUE(registeredSymbols.find("H-1") != registeredSymbols.end());
EXPECT_TRUE(registeredSymbols.find("He-4") != registeredSymbols.end());
EXPECT_TRUE(registeredSymbols.find("H-19") == registeredSymbols.end());
EXPECT_TRUE(registeredSymbols.find("He-21") == registeredSymbols.end());
EXPECT_TRUE(registeredSymbols.contains("H-1"));
EXPECT_TRUE(registeredSymbols.contains("He-4"));
EXPECT_TRUE(!registeredSymbols.contains("H-19"));
EXPECT_TRUE(!registeredSymbols.contains("He-21"));
}
/**
@@ -433,3 +433,372 @@ TEST_F(compositionTest, getSpeciesFromAZ) {
EXPECT_EQ(fourdst::atomic::O_12, fourdst::atomic::az_to_species(12, 8));
}
TEST_F(compositionTest, constructorWithSymbolsVectorAndSet) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
std::vector<std::string> vs = {"H-1", "He-4"};
std::set<std::string> ss = {"H-1", "He-4"};
fourdst::composition::Composition compVec(vs);
EXPECT_TRUE(compVec.hasSymbol("H-1"));
EXPECT_TRUE(compVec.hasSymbol("He-4"));
fourdst::composition::Composition compSet(ss);
EXPECT_TRUE(compSet.hasSymbol("H-1"));
EXPECT_TRUE(compSet.hasSymbol("He-4"));
}
TEST_F(compositionTest, constructorWithSymbolsAndFractionsMassAndNumber) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
using fourdst::composition::Composition;
using fourdst::atomic::species;
// Mass-fraction constructor
std::vector<std::string> symM = {"H-1", "He-4"};
std::vector<double> fracM = {0.6, 0.4};
Composition compM(symM, fracM, true);
EXPECT_NEAR(compM.getMassFraction("H-1"), 0.6, 1e-12);
EXPECT_NEAR(compM.getMassFraction("He-4"), 0.4, 1e-12);
// Mean particle mass and specific number density are reciprocals
double sn = 0.6/species.at("H-1").mass() + 0.4/species.at("He-4").mass();
double mp = 1.0/sn;
EXPECT_NEAR(compM.getMeanParticleMass(), mp, 1e-12);
// Number-fraction constructor
std::vector<std::string> symN = {"H-1", "He-4"};
std::vector<double> fracN = {0.9, 0.1};
Composition compN(symN, fracN, false);
EXPECT_NEAR(compN.getNumberFraction("H-1"), 0.9, 1e-12);
EXPECT_NEAR(compN.getNumberFraction("He-4"), 0.1, 1e-12);
double meanA = 0.9*species.at("H-1").mass() + 0.1*species.at("He-4").mass();
EXPECT_NEAR(compN.getMeanParticleMass(), meanA, 1e-12);
// Check converted mass fractions X_i = n_i * A_i / <A>
double xH = 0.9*species.at("H-1").mass()/meanA;
double xHe = 0.1*species.at("He-4").mass()/meanA;
compN.setCompositionMode(true);
EXPECT_NEAR(compN.getMassFraction("H-1"), xH, 1e-12);
EXPECT_NEAR(compN.getMassFraction("He-4"), xHe, 1e-12);
}
TEST_F(compositionTest, registerSymbolVectorAndSingleSpeciesAndModeMismatch) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
using fourdst::composition::Composition;
Composition comp;
comp.registerSymbol(std::vector<std::string>{"H-1", "He-4"});
EXPECT_TRUE(comp.hasSymbol("H-1"));
EXPECT_TRUE(comp.hasSymbol("He-4"));
// Register by Species
Composition comp2;
comp2.registerSpecies(fourdst::atomic::H_1);
comp2.registerSpecies(fourdst::atomic::He_4);
EXPECT_TRUE(comp2.hasSymbol("H-1"));
EXPECT_TRUE(comp2.hasSymbol("He-4"));
// Mode mismatch should throw
Composition comp3;
comp3.registerSymbol("H-1", true); // mass mode
EXPECT_THROW(comp3.registerSymbol("He-4", false), fourdst::composition::exceptions::CompositionModeError);
}
TEST_F(compositionTest, setMassFractionBySpeciesAndVector) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
using fourdst::composition::Composition;
using fourdst::atomic::H_1;
using fourdst::atomic::He_4;
Composition comp;
comp.registerSymbol("H-1");
comp.registerSymbol("He-4");
// Single species overload
double old = comp.setMassFraction(H_1, 0.7);
EXPECT_NEAR(old, 0.0, 1e-15);
old = comp.setMassFraction(He_4, 0.3);
EXPECT_NEAR(old, 0.0, 1e-15);
// Vector overload
std::vector<fourdst::atomic::Species> sp = {H_1, He_4};
std::vector<double> xs = {0.6, 0.4};
auto olds = comp.setMassFraction(sp, xs);
ASSERT_EQ(olds.size(), 2u);
EXPECT_NEAR(olds[0], 0.7, 1e-12);
EXPECT_NEAR(olds[1], 0.3, 1e-12);
EXPECT_TRUE(comp.finalize());
EXPECT_NEAR(comp.getMassFraction("H-1"), 0.6, 1e-12);
EXPECT_NEAR(comp.getMassFraction(He_4), 0.4, 1e-12);
}
TEST_F(compositionTest, setNumberFractionOverloads) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
using fourdst::composition::Composition;
using fourdst::atomic::H_1;
using fourdst::atomic::He_4;
Composition comp;
comp.registerSymbol("H-1", false);
comp.registerSymbol("He-4", false);
// Single symbol
double old = comp.setNumberFraction("H-1", 0.8);
EXPECT_NEAR(old, 0.0, 1e-15);
// Vector of symbols
auto oldv = comp.setNumberFraction(std::vector<std::string>{"H-1", "He-4"}, std::vector<double>{0.75, 0.25});
ASSERT_EQ(oldv.size(), 2u);
EXPECT_NEAR(oldv[0], 0.8, 1e-12);
EXPECT_NEAR(oldv[1], 0.0, 1e-12);
// Species and vector<Species>
old = comp.setNumberFraction(H_1, 0.7);
EXPECT_NEAR(old, 0.75, 1e-12);
auto oldsv = comp.setNumberFraction(std::vector<fourdst::atomic::Species>{H_1, He_4}, std::vector<double>{0.6, 0.4});
ASSERT_EQ(oldsv.size(), 2u);
EXPECT_NEAR(oldsv[0], 0.7, 1e-12);
EXPECT_NEAR(oldsv[1], 0.25, 1e-12);
EXPECT_TRUE(comp.finalize());
EXPECT_NEAR(comp.getNumberFraction("H-1"), 0.6, 1e-12);
EXPECT_NEAR(comp.getNumberFraction(He_4), 0.4, 1e-12);
}
TEST_F(compositionTest, mixErrorCases) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
using fourdst::composition::Composition;
Composition a; a.registerSymbol("H-1"); a.registerSymbol("He-4"); a.setMassFraction("H-1", 0.6); a.setMassFraction("He-4", 0.4); a.finalize();
Composition b; b.registerSymbol("H-1"); b.registerSymbol("He-4"); b.setMassFraction("H-1", 0.5); b.setMassFraction("He-4", 0.5);
// Not finalized second comp
EXPECT_THROW(static_cast<void>(a.mix(b, 0.5)), fourdst::composition::exceptions::CompositionNotFinalizedError);
b.finalize();
// Invalid fraction
EXPECT_THROW(static_cast<void>(a.mix(b, -0.1)), fourdst::composition::exceptions::InvalidCompositionError);
EXPECT_THROW(static_cast<void>(a.mix(b, 1.1)), fourdst::composition::exceptions::InvalidCompositionError);
}
TEST_F(compositionTest, getMassAndNumberFractionMapsAndSpeciesOverloads) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
using fourdst::composition::Composition;
Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4);
ASSERT_TRUE(comp.finalize());
auto m = comp.getMassFraction();
ASSERT_EQ(m.size(), 2u);
EXPECT_NEAR(m.at("H-1"), 0.6, 1e-12);
EXPECT_NEAR(m.at("He-4"), 0.4, 1e-12);
EXPECT_NEAR(comp.getMassFraction(fourdst::atomic::H_1), 0.6, 1e-12);
EXPECT_NEAR(comp.getMolarAbundance(fourdst::atomic::H_1), m.at("H-1")/fourdst::atomic::H_1.mass(), 1e-12);
// Switch to number-fraction mode and verify number maps
comp.setCompositionMode(false);
// Must re-finalize after modifications (mode switch itself keeps values consistent but not finalized status changed? setCompositionMode requires to be finalized; here we just switched modes)
// Set specific number fractions and finalize
comp.setNumberFraction("H-1", 0.7);
comp.setNumberFraction("He-4", 0.3);
ASSERT_TRUE(comp.finalize());
auto n = comp.getNumberFraction();
ASSERT_EQ(n.size(), 2u);
EXPECT_NEAR(n.at("H-1"), 0.7, 1e-12);
EXPECT_NEAR(n.at("He-4"), 0.3, 1e-12);
EXPECT_NEAR(comp.getNumberFraction(fourdst::atomic::He_4), 0.3, 1e-12);
}
TEST_F(compositionTest, meanAtomicNumberAndElectronAbundance) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
using fourdst::atomic::species;
using fourdst::composition::Composition;
Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4);
ASSERT_TRUE(comp.finalize());
// Compute expected Ye = sum(X_i * Z_i / A_i)
constexpr double xH = 0.6, xHe = 0.4;
const double aH = species.at("H-1").a();
const double aHe = species.at("He-4").a();
const double zH = species.at("H-1").z();
const double zHe = species.at("He-4").z();
const double expectedYe = xH*zH/aH + xHe*zHe/aHe;
EXPECT_NEAR(comp.getElectronAbundance(), expectedYe, 1e-12);
// <Z> = <A> * sum(X_i * Z_i / A_i)
const double expectedZ = comp.getMeanParticleMass() * expectedYe;
EXPECT_NEAR(comp.getMeanAtomicNumber(), expectedZ, 1e-12);
}
TEST_F(compositionTest, canonicalCompositionAndCaching) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
using fourdst::composition::Composition;
Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4);
ASSERT_TRUE(comp.finalize());
auto canon1 = comp.getCanonicalComposition();
EXPECT_NEAR(canon1.X, 0.6, 1e-12);
EXPECT_NEAR(canon1.Y, 0.4, 1e-12);
EXPECT_NEAR(canon1.Z, 0.0, 1e-12);
// Call again to exercise caching code path
auto canon2 = comp.getCanonicalComposition();
EXPECT_NEAR(canon2.X, 0.6, 1e-12);
EXPECT_NEAR(canon2.Y, 0.4, 1e-12);
EXPECT_NEAR(canon2.Z, 0.0, 1e-12);
// Add a metal and re-check
Composition comp2; comp2.registerSymbol("H-1"); comp2.registerSymbol("He-4"); comp2.registerSymbol("O-16");
comp2.setMassFraction("H-1", 0.6); comp2.setMassFraction("He-4", 0.35); comp2.setMassFraction("O-16", 0.05);
ASSERT_TRUE(comp2.finalize());
auto canon3 = comp2.getCanonicalComposition(true);
EXPECT_NEAR(canon3.X, 0.6, 1e-12);
EXPECT_NEAR(canon3.Y, 0.35, 1e-12);
EXPECT_NEAR(canon3.Z, 0.05, 1e-12);
}
TEST_F(compositionTest, vectorsAndIndexingAndSpeciesAtIndex) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
using fourdst::composition::Composition;
using fourdst::atomic::species;
Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); comp.registerSymbol("O-16");
comp.setMassFraction("H-1", 0.5); comp.setMassFraction("He-4", 0.3); comp.setMassFraction("O-16", 0.2);
ASSERT_TRUE(comp.finalize());
// Mass fraction vector sorted by mass: H-1, He-4, O-16
auto mv = comp.getMassFractionVector();
ASSERT_EQ(mv.size(), 3u);
EXPECT_NEAR(mv[0], 0.5, 1e-12);
EXPECT_NEAR(mv[1], 0.3, 1e-12);
EXPECT_NEAR(mv[2], 0.2, 1e-12);
// Species indices ordering
size_t iH = comp.getSpeciesIndex("H-1");
size_t iHe = comp.getSpeciesIndex("He-4");
size_t iO = comp.getSpeciesIndex("O-16");
EXPECT_LT(iH, iHe);
EXPECT_LT(iHe, iO);
EXPECT_EQ(comp.getSpeciesIndex(fourdst::atomic::H_1), iH);
EXPECT_EQ(comp.getSpeciesIndex(species.at("He-4")), iHe);
// Species at index
auto sAtHe = comp.getSpeciesAtIndex(iHe);
EXPECT_EQ(std::string(sAtHe.name()), std::string("He-4"));
// Number fraction vector after switching modes
comp.setCompositionMode(false);
// Tweak number fractions and finalize
// Compute expected number fractions from original mass fractions first
double denom = 0.5/species.at("H-1").mass() + 0.3/species.at("He-4").mass() + 0.2/species.at("O-16").mass();
double nH_exp = (0.5/species.at("H-1").mass())/denom;
double nHe_exp = (0.3/species.at("He-4").mass())/denom;
double nO_exp = (0.2/species.at("O-16").mass())/denom;
auto nv0 = comp.getNumberFractionVector();
ASSERT_EQ(nv0.size(), 3u);
EXPECT_NEAR(nv0[iH], nH_exp, 1e-12);
EXPECT_NEAR(nv0[iHe], nHe_exp, 1e-12);
EXPECT_NEAR(nv0[iO], nO_exp, 1e-12);
// Molar abundance vector X_i/A_i in mass mode; switch back to mass mode to verify
comp.setCompositionMode(true);
comp.finalize(true);
auto av = comp.getMolarAbundanceVector();
ASSERT_EQ(av.size(), 3u);
EXPECT_NEAR(av[iH], 0.5/species.at("H-1").mass(), 1e-12);
EXPECT_NEAR(av[iHe], 0.3/species.at("He-4").mass(), 1e-12);
EXPECT_NEAR(av[iO], 0.2/species.at("O-16").mass(), 1e-12);
}
TEST_F(compositionTest, containsAndPreFinalizationGuards) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
fourdst::composition::Composition comp;
comp.registerSymbol("H-1"); comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4);
// contains should throw before finalize
EXPECT_THROW(static_cast<void>(comp.contains(fourdst::atomic::H_1)), fourdst::composition::exceptions::CompositionNotFinalizedError);
ASSERT_TRUE(comp.finalize());
EXPECT_TRUE(comp.contains(fourdst::atomic::H_1));
EXPECT_FALSE(comp.contains(fourdst::atomic::Li_6));
}
TEST_F(compositionTest, subsetNoneMethodAndNormalizationFlow) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
fourdst::composition::Composition comp;
comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); comp.registerSymbol("O-16");
comp.setMassFraction("H-1", 0.5); comp.setMassFraction("He-4", 0.3); comp.setMassFraction("O-16", 0.2);
ASSERT_TRUE(comp.finalize());
fourdst::composition::Composition sub = comp.subset(std::vector<std::string>{"H-1", "He-4"}, "none");
// Not normalized: finalize without normalization should fail
EXPECT_FALSE(sub.finalize(false));
// With normalization, it should succeed and scale to sum to 1
EXPECT_TRUE(sub.finalize(true));
double sum = sub.getMassFraction("H-1") + sub.getMassFraction("He-4");
EXPECT_NEAR(sum, 1.0, 1e-12);
EXPECT_NEAR(sub.getMassFraction("H-1"), 0.5/(0.5+0.3), 1e-12);
EXPECT_NEAR(sub.getMassFraction("He-4"), 0.3/(0.5+0.3), 1e-12);
}
TEST_F(compositionTest, copyConstructorAndAssignmentIndependence) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
using fourdst::composition::Composition;
Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4);
ASSERT_TRUE(comp.finalize());
Composition copy(comp); // copy ctor
EXPECT_NEAR(copy.getMassFraction("H-1"), 0.6, 1e-12);
EXPECT_NEAR(copy.getMassFraction("He-4"), 0.4, 1e-12);
Composition assigned; assigned = comp; // assignment
EXPECT_NEAR(assigned.getMassFraction("H-1"), 0.6, 1e-12);
EXPECT_NEAR(assigned.getMassFraction("He-4"), 0.4, 1e-12);
// Modify original and ensure copies do not change
comp.setMassFraction("H-1", 0.7); comp.setMassFraction("He-4", 0.3); ASSERT_TRUE(comp.finalize());
EXPECT_NEAR(copy.getMassFraction("H-1"), 0.6, 1e-12);
EXPECT_NEAR(assigned.getMassFraction("He-4"), 0.4, 1e-12);
}
TEST_F(compositionTest, getCompositionBySpeciesAndAllEntries) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
using fourdst::composition::Composition;
Composition comp; comp.registerSymbol("H-1"); comp.registerSymbol("He-4");
comp.setMassFraction("H-1", 0.6); comp.setMassFraction("He-4", 0.4);
ASSERT_TRUE(comp.finalize());
auto pairBySpec = comp.getComposition(fourdst::atomic::H_1);
EXPECT_NEAR(pairBySpec.first.mass_fraction(), 0.6, 1e-12);
EXPECT_NEAR(pairBySpec.second.meanParticleMass, comp.getMeanParticleMass(), 1e-15);
auto all = comp.getComposition();
ASSERT_EQ(all.first.size(), 2u);
EXPECT_NEAR(all.first.at("H-1").mass_fraction(), 0.6, 1e-12);
EXPECT_NEAR(all.first.at("He-4").mass_fraction(), 0.4, 1e-12);
EXPECT_NEAR(all.second.meanParticleMass, comp.getMeanParticleMass(), 1e-15);
}
TEST_F(compositionTest, iterationBeginEndAndIndexOutOfRange) {
fourdst::config::Config::getInstance().loadConfig(EXAMPLE_FILENAME);
fourdst::composition::Composition comp;
comp.registerSymbol("H-1"); comp.registerSymbol("He-4"); comp.registerSymbol("O-16");
comp.setMassFraction("H-1", 0.5); comp.setMassFraction("He-4", 0.3); comp.setMassFraction("O-16", 0.2);
ASSERT_TRUE(comp.finalize());
// Iterate and count entries
size_t count = 0;
for (auto it = comp.begin(); it != comp.end(); ++it) {
count++;
}
EXPECT_EQ(count, 3u);
// Out-of-range access
EXPECT_THROW(static_cast<void>(comp.getSpeciesAtIndex(100)), std::out_of_range);
}