feat(eos): EOS now uses composition module

EOS code updated to make use of composition module. This is part of the larger change to change all composition handling to use the composition module. Note that the current implimentation is a bit hacky since it simply copies back and forth to the alredy used HELMEOSInput and HELMEOSOutput structs. In furture this can be more tightly connected to avoid extra copies.

src/composition/private/composition.cpp

@@ -25,6 +25,7 @@
 #include <unordered_map>
 #include <vector>
 #include <array>
+#include <ranges>
 
 #include <utility>
 
@@ -32,13 +33,15 @@
 
 namespace serif::composition {
 
-CompositionEntry::CompositionEntry() : m_symbol("H-1"), m_isotope(chemSpecies::species.at("H-1")), m_initialized(false) {}
+    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) {
+    CompositionEntry::CompositionEntry(const std::string& symbol, const bool massFracMode) : m_symbol(symbol), m_isotope(chemSpecies::species.at(symbol)), m_massFracMode(massFracMode) {
         setSpecies(symbol);
-}
+    }
 
-CompositionEntry::CompositionEntry(const CompositionEntry& entry) : 
+    CompositionEntry::CompositionEntry(const CompositionEntry& entry) :
         m_symbol(entry.m_symbol),
         m_isotope(entry.m_isotope),
         m_massFracMode(entry.m_massFracMode),
@@ -47,7 +50,7 @@ CompositionEntry::CompositionEntry(const CompositionEntry& entry) :
         m_relAbundance(entry.m_relAbundance),
         m_initialized(entry.m_initialized) {}
 
-void CompositionEntry::setSpecies(const std::string& symbol) {
+    void CompositionEntry::setSpecies(const std::string& symbol) {
         if (m_initialized) {
             throw std::runtime_error("Composition entry is already initialized.");
         }
@@ -57,100 +60,100 @@ void CompositionEntry::setSpecies(const std::string& symbol) {
         m_symbol = symbol;
         m_isotope = chemSpecies::species.at(symbol);
         m_initialized = true;
-}
+    }
 
-std::string CompositionEntry::symbol() const {
+    std::string CompositionEntry::symbol() const {
         return m_symbol;
-}
+    }
 
-double CompositionEntry::mass_fraction() const {
+    double CompositionEntry::mass_fraction() const {
         if (!m_massFracMode) {
             throw std::runtime_error("Composition entry is in number fraction mode.");
         }
         return m_massFraction;
-}
+    }
 
-double CompositionEntry::mass_fraction(double meanMolarMass) const {
+    double CompositionEntry::mass_fraction(double meanMolarMass) const {
         if (m_massFracMode) {
             return m_massFraction;
         }
         return m_relAbundance / meanMolarMass;
-}
+    }
 
 
-double CompositionEntry::number_fraction() const {
+    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 {
+    double CompositionEntry::number_fraction(double totalMoles) const {
         if (m_massFracMode) {
             return m_relAbundance / totalMoles;
         }
         return m_numberFraction;
-}
+    }
 
-double CompositionEntry::rel_abundance() const {
+    double CompositionEntry::rel_abundance() const {
         return m_relAbundance;
-}
+    }
 
-chemSpecies::Species CompositionEntry::isotope() const {
+    chemSpecies::Species CompositionEntry::isotope() const {
         return m_isotope;
-}
+    }
 
-void CompositionEntry::setMassFraction(double mass_fraction) {
+    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) {
+    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) {
+    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) {
+    bool CompositionEntry::setNumberFracMode(double specificNumberDensity) {
         if (!m_massFracMode) {
             return false;
         }
         m_massFracMode = false;
         m_numberFraction = m_relAbundance / specificNumberDensity;
         return true;
-}
+    }
 
-bool CompositionEntry::getMassFracMode() const {
+    bool CompositionEntry::getMassFracMode() const {
         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, bool massFracMode) : m_massFracMode(massFracMode) {
+    Composition::Composition(const std::vector<std::string>& symbols, const std::vector<double>& 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.");
@@ -170,9 +173,32 @@ Composition::Composition(const std::vector<std::string>& symbols, const std::vec
             }
         }
         finalize();
-}
+    }
 
-void Composition::registerSymbol(const std::string& symbol, bool massFracMode) {
+    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.");
@@ -197,26 +223,26 @@ void Composition::registerSymbol(const std::string& symbol, bool massFracMode) {
         CompositionEntry entry(symbol, m_massFracMode);
         m_compositions[symbol] = entry;
         LOG_INFO(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, bool massFracMode) {
         for (const auto& symbol : symbols) {
             registerSymbol(symbol, massFracMode);
         }
-}
+    }
 
-std::set<std::string> Composition::getRegisteredSymbols() const {
+    std::set<std::string> Composition::getRegisteredSymbols() const {
         return m_registeredSymbols;
-}
+    }
 
-void Composition::validateComposition(const std::vector<double>& fractions) const {
+    void Composition::validateComposition(const std::vector<double>& fractions) const {
         if (!isValidComposition(fractions)) {
             LOG_ERROR(m_logger, "Invalid composition.");
             throw std::runtime_error("Invalid composition.");
         }
-}
+    }
 
-bool Composition::isValidComposition(const std::vector<double>& fractions) const {
+    bool Composition::isValidComposition(const std::vector<double>& fractions) const {
         double sum = 0.0;
         for (const auto& fraction : fractions) {
             sum += fraction;
@@ -227,13 +253,13 @@ bool Composition::isValidComposition(const std::vector<double>& fractions) const
         }
 
         return true;
-}
+    }
 
-bool Composition::isValidSymbol(const std::string& symbol) const {
+    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) {
+    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.");
@@ -254,9 +280,9 @@ double Composition::setMassFraction(const std::string& symbol, const double& mas
         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) {
+    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.");
             throw std::runtime_error("The number of symbols and mass fractions must be equal.");
@@ -268,9 +294,9 @@ std::vector<double> Composition::setMassFraction(const std::vector<std::string>&
             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) {
+    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.");
@@ -291,9 +317,9 @@ double Composition::setNumberFraction(const std::string& symbol, const double& n
         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.");
             throw std::runtime_error("The number of symbols and number fractions must be equal.");
@@ -305,9 +331,9 @@ std::vector<double> Composition::setNumberFraction(const std::vector<std::string
             old_number_fractions.push_back(setNumberFraction(symbols[i], number_fractions[i]));
         }
         return old_number_fractions;
-}
+    }
 
-bool Composition::finalize(bool norm) {
+    bool Composition::finalize(bool norm) {
         bool finalized = false;
         if (m_massFracMode) {
             finalized = finalizeMassFracMode(norm);
@@ -318,9 +344,9 @@ bool Composition::finalize(bool norm) {
             m_finalized = true;
         }
         return finalized;
-}
+    }
 
-bool Composition::finalizeMassFracMode(bool norm) {
+    bool Composition::finalizeMassFracMode(bool norm) {
         std::vector<double> mass_fractions;
         mass_fractions.reserve(m_compositions.size());
         for (const auto& [_, entry] : m_compositions) {
@@ -354,9 +380,9 @@ bool Composition::finalizeMassFracMode(bool norm) {
         }
         m_meanParticleMass = 1.0/m_specificNumberDensity;
         return true;
-}
+    }
 
-bool Composition::finalizeNumberFracMode(bool norm) {
+    bool Composition::finalizeNumberFracMode(bool norm) {
         std::vector<double> number_fractions;
         number_fractions.reserve(m_compositions.size());
         for (const auto& [_, entry] : m_compositions) {
@@ -387,9 +413,9 @@ bool Composition::finalizeNumberFracMode(bool norm) {
         }
         m_specificNumberDensity = 1.0/m_meanParticleMass;
         return true;
-}
+    }
 
-Composition Composition::mix(const Composition& other, double fraction) const {
+    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()).");
@@ -416,9 +442,9 @@ Composition Composition::mix(const Composition& other, double fraction) const {
         }
         mixedComposition.finalize();
         return mixedComposition;
-}
+    }
 
-double Composition::getMassFraction(const std::string& symbol) const {
+    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()).");
@@ -432,18 +458,18 @@ double Composition::getMassFraction(const std::string& symbol) const {
         } else {
             return m_compositions.at(symbol).mass_fraction(m_meanParticleMass);
         }
-}
+    }
 
-std::unordered_map<std::string, double> Composition::getMassFraction() const {
+    std::unordered_map<std::string, double> Composition::getMassFraction() const {
         std::unordered_map<std::string, double> 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 {
+    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()).");
@@ -457,17 +483,17 @@ double Composition::getNumberFraction(const std::string& symbol) const {
         } else {
             return m_compositions.at(symbol).number_fraction(m_specificNumberDensity);
         }
-}
+    }
 
-std::unordered_map<std::string, double> Composition::getNumberFraction() const {
+    std::unordered_map<std::string, double> Composition::getNumberFraction() const {
         std::unordered_map<std::string, double> number_fractions;
         for (const auto& [symbol, entry] : m_compositions) {
             number_fractions[symbol] = getNumberFraction(symbol);
         }
         return number_fractions;
-}
+    }
 
-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.");
             throw std::runtime_error("Composition has not been finalized (Consider running .finalize()).");
@@ -477,17 +503,53 @@ std::pair<CompositionEntry, GlobalComposition> Composition::getComposition(const
             throw std::runtime_error("Symbol is not in the composition.");
         }
         return {m_compositions.at(symbol), {m_specificNumberDensity, m_meanParticleMass}};
-}
+    }
 
-std::pair<std::unordered_map<std::string, CompositionEntry>, GlobalComposition> Composition::getComposition() const {
+    std::pair<std::unordered_map<std::string, CompositionEntry>, 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<std::string>& symbols, std::string method) const {
+    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<std::string>& symbols, std::string method) const {
         std::array<std::string, 2> methods = {"norm", "none"};
 
         if (std::find(methods.begin(), methods.end(), method) == methods.end()) {
@@ -514,9 +576,9 @@ Composition Composition::subset(const std::vector<std::string>& symbols, std::st
             }
         }
         return subsetComposition;
-}
+    }
 
-void Composition::setCompositionMode(bool massFracMode) {
+    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.");
@@ -535,36 +597,36 @@ void Composition::setCompositionMode(bool massFracMode) {
             }
         }
         m_massFracMode = massFracMode;
-}
+    }
 
-bool Composition::hasSymbol(const std::string& symbol) const {
+    bool Composition::hasSymbol(const std::string& symbol) const {
         return m_compositions.count(symbol) > 0;
-}
+    }
 
-/// OVERLOADS
+    /// 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: \n";
         for (const auto& [symbol, entry] : composition.m_compositions) {
             os << entry << "\n";
         }
         return os;
-}
+    }
 
 } // namespace serif::composition

src/composition/public/composition.h

@@ -277,7 +277,7 @@ namespace serif::composition {
          * Composition comp(symbols);
          * @endcode
          */
-        Composition(const std::vector<std::string>& symbols);
+        explicit Composition(const std::vector<std::string>& 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<std::string>& symbols);
+        explicit Composition(const std::set<std::string>& symbols);
 
         /**
          * @brief Constructs a Composition with the given symbols and mass fractions.
@@ -304,6 +304,14 @@ namespace serif::composition {
          */
         Composition(const std::vector<std::string>& symbols, const std::vector<double>& 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<std::string> getRegisteredSymbols() const;
+        [[nodiscard]] std::set<std::string> 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<std::string, double> getMassFraction() const;
+        [[nodiscard]] std::unordered_map<std::string, double> 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<std::string, double> getNumberFraction() const;
+        [[nodiscard]] std::unordered_map<std::string, double> 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<CompositionEntry, GlobalComposition> getComposition(const std::string& symbol) const;
+        [[nodiscard]] std::pair<CompositionEntry, GlobalComposition> 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<std::unordered_map<std::string, CompositionEntry>, GlobalComposition> getComposition() const;
+        [[nodiscard]] std::pair<std::unordered_map<std::string, CompositionEntry>, 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.