Merge pull request #68 from tboudreaux/refactor/useComposition

Refactor Network and EOS modules to use composition module

assets/static/atomic/meson.build

@@ -1,3 +1,4 @@
 species_weight_dep = declare_dependency(
     include_directories: include_directories('include'),
 )
+message('✅ SERiF species_weight dependency declared')

src/composition/meson.build

@@ -29,4 +29,4 @@ composition_dep = declare_dependency(
 )
 
 # Make headers accessible
-install_headers(composition_headers, subdir : '4DSSE/composition')
+install_headers(composition_headers, subdir : 'SERiF/composition')

src/composition/private/composition.cpp

@@ -25,6 +25,7 @@
 #include <unordered_map>
 #include <vector>
 #include <array>
+#include <ranges>
 
 #include <utility>
 
@@ -32,9 +33,11 @@
 
 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);
     }
 
@@ -172,6 +175,29 @@ Composition::Composition(const std::vector<std::string>& symbols, const std::vec
         finalize();
     }
 
+    Composition::Composition(const Composition &composition) {
+        m_finalized = composition.m_finalized;
+        m_specificNumberDensity = composition.m_specificNumberDensity;
+        m_meanParticleMass = composition.m_meanParticleMass;
+        m_massFracMode = composition.m_massFracMode;
+        m_registeredSymbols = composition.m_registeredSymbols;
+        m_compositions = composition.m_compositions;
+    }
+
+    Composition& Composition::operator=(const Composition &other) {
+        if (this != &other) {
+            m_finalized               = other.m_finalized;
+            m_specificNumberDensity   = other.m_specificNumberDensity;
+            m_meanParticleMass        = other.m_meanParticleMass;
+            m_massFracMode            = other.m_massFracMode;
+            m_registeredSymbols       = other.m_registeredSymbols;
+            m_compositions            = other.m_compositions;
+            // note: m_config remains bound to the same singleton, so we skip it
+        }
+        return *this;
+
+    }
+
     void Composition::registerSymbol(const std::string& symbol, bool massFracMode) {
         if (!isValidSymbol(symbol)) {
             LOG_ERROR(m_logger, "Invalid symbol: {}", symbol);
@@ -229,12 +255,12 @@ 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) {
-    return chemSpecies::species.count(symbol) > 0;
+        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.");
         }
@@ -250,7 +276,7 @@ double Composition::setMassFraction(const std::string& symbol, const double& mas
         }
 
         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;
@@ -307,7 +333,7 @@ std::vector<double> Composition::setNumberFraction(const std::vector<std::string
         return old_number_fractions;
     }
 
-bool Composition::finalize(bool norm) {
+    bool Composition::finalize(const bool norm) {
         bool finalized = false;
         if (m_massFracMode) {
             finalized = finalizeMassFracMode(norm);
@@ -406,9 +432,9 @@ Composition Composition::mix(const Composition& other, double fraction) const {
 
         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);
@@ -423,7 +449,7 @@ double Composition::getMassFraction(const std::string& symbol) const {
             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.");
         }
@@ -436,7 +462,7 @@ double Composition::getMassFraction(const std::string& symbol) 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) {
+        for (const auto &symbol: m_compositions | std::views::keys) {
             mass_fractions[symbol] = getMassFraction(symbol);
         }
         return mass_fractions;
@@ -448,7 +474,7 @@ double Composition::getNumberFraction(const std::string& symbol) const {
             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.");
         }
@@ -461,7 +487,7 @@ double Composition::getNumberFraction(const std::string& symbol) 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) {
+        for (const auto &symbol: m_compositions | std::views::keys) {
             number_fractions[symbol] = getNumberFraction(symbol);
         }
         return number_fractions;
@@ -472,7 +498,7 @@ std::pair<CompositionEntry, GlobalComposition> Composition::getComposition(const
             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,18 +513,43 @@ std::pair<std::unordered_map<std::string, CompositionEntry>, GlobalComposition>
         return {m_compositions, {m_specificNumberDensity, m_meanParticleMass}};
     }
 
-Composition Composition::subset(const std::vector<std::string>& symbols, std::string method) const {
+    double Composition::getMeanParticleMass() const {
-    std::array<std::string, 2> methods = {"norm", "none"};
+        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;
+    }
 
-    if (std::find(methods.begin(), methods.end(), method) == methods.end()) {
+    double Composition::getMeanAtomicNumber() const {
-        std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'.";
+        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 zSum = 0.0;
+
+        // Loop through all registered species in the composition.
+        for (const auto &val: m_compositions | std::views::values) {
+            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 {
+        const std::array<std::string, 2> methods = {"norm", "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 {
@@ -507,7 +558,7 @@ Composition Composition::subset(const std::vector<std::string>& symbols, std::st
             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.");
@@ -516,14 +567,14 @@ Composition Composition::subset(const std::vector<std::string>& symbols, std::st
         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 {
@@ -537,6 +588,53 @@ void Composition::setCompositionMode(bool massFracMode) {
         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;
     }

src/composition/public/composition.h

@@ -18,8 +18,7 @@
 //   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 //
 // *********************************************************************** */
-#ifndef COMPOSITION_H
+#pragma once
-#define COMPOSITION_H
 
 #include <iostream>
 #include <string>
@@ -34,6 +33,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.
      */
@@ -68,7 +81,7 @@ namespace serif::composition {
          * @brief Constructs a CompositionEntry with the given symbol and mode.
          * @param symbol The chemical symbol of the species.
          * @param massFracMode True if mass fraction mode, false if number fraction mode.
-         * @example
+         * *Example Usage:*
          * @code
          * CompositionEntry entry("H", true);
          * @endcode
@@ -185,13 +198,13 @@ namespace serif::composition {
      * - The only exception to the finalize rule is if the compositon was constructed with symbols and mass fractions at instantiation time. In this case, the composition is automatically finalized.  
      * however, this means that the composition passed to the constructor must be valid.
      * 
-     * @example Constructing a finalized composition with symbols and mass fractions:
+     * *Example Usage:* Constructing a finalized composition with symbols and mass fractions:
      * @code
      * std::vector<std::string> symbols = {"H", "He"};
      * std::vector<double> mass_fractions = {0.7, 0.3};
      * Composition comp(symbols, mass_fractions);
      * @endcode
-     * @example Constructing a composition with symbols and finalizing it later:
+     * *Example Usage:* Constructing a composition with symbols and finalizing it later:
      * @code
      * std::vector<std::string> symbols = {"H", "He"};
      * Composition comp(symbols);
@@ -220,7 +233,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.
@@ -271,31 +284,31 @@ namespace serif::composition {
         /**
          * @brief Constructs a Composition with the given symbols.
          * @param symbols The symbols to initialize the composition with.
-         * @example
+         * *Example Usage:*
          * @code
          * std::vector<std::string> symbols = {"H", "O"};
          * 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.
          * @param symbols The symbols to initialize the composition with.
-         * @example
+         * *Example Usage:*
          * @code
          * std::set<std::string> symbols = {"H", "O"};
          * 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.
          * @param symbols The symbols to initialize the composition with.
          * @param mass_fractions The mass fractions corresponding to the symbols.
          * @param massFracMode True if mass fraction mode, false if number fraction mode.
-         * @example
+         * *Example Usage:*
          * @code
          * std::vector<std::string> symbols = {"H", "O"};
          * std::vector<double> mass_fractions = {0.1, 0.9};
@@ -304,11 +317,19 @@ 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.
          * @param massFracMode True if mass fraction mode, false if number fraction mode.
-         * @example
+         * *Example Usage:*
          * @code
          * Composition comp;
          * comp.registerSymbol("H");
@@ -320,7 +341,7 @@ namespace serif::composition {
          * @brief Registers multiple new symbols.
          * @param symbols The symbols to register.
          * @param massFracMode True if mass fraction mode, false if number fraction mode.
-         * @example
+         * *Example Usage:*
          * @code
          * std::vector<std::string> symbols = {"H", "O"};
          * Composition comp;
@@ -333,14 +354,14 @@ 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.
          * @param symbol The symbol to set the mass fraction for.
          * @param mass_fraction The mass fraction to set.
          * @return The mass fraction that was set.
-         * @example
+         * *Example Usage:*
          * @code
          * Composition comp;
          * comp.setMassFraction("H", 0.1);
@@ -353,7 +374,7 @@ namespace serif::composition {
          * @param symbols The symbols to set the mass fraction for.
          * @param mass_fractions The mass fractions corresponding to the symbols.
          * @return A vector of mass fractions that were set.
-         * @example
+         * *Example Usage:*
          * @code
          * std::vector<std::string> symbols = {"H", "O"};
          * std::vector<double> mass_fractions = {0.1, 0.9};
@@ -390,40 +411,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.
@@ -446,6 +479,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.
@@ -464,5 +507,3 @@ namespace serif::composition {
 
     };
 }; // namespace serif::composition
-
-#endif // COMPOSITION_H