From 39f9c1856e125553352cfa1aa5cbf94cec47253c Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 16 Jun 2025 13:15:34 -0400 Subject: [PATCH 1/8] docs(composition): @example -> *Example Usage:* --- src/composition/public/composition.h | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/composition/public/composition.h b/src/composition/public/composition.h index 7226393..e61a5c6 100644 --- a/src/composition/public/composition.h +++ b/src/composition/public/composition.h @@ -68,7 +68,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 +185,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 symbols = {"H", "He"}; * std::vector 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 symbols = {"H", "He"}; * Composition comp(symbols); @@ -271,7 +271,7 @@ 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 symbols = {"H", "O"}; * Composition comp(symbols); @@ -282,7 +282,7 @@ namespace serif::composition { /** * @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 symbols = {"H", "O"}; * Composition comp(symbols); @@ -295,7 +295,7 @@ namespace serif::composition { * @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 symbols = {"H", "O"}; * std::vector mass_fractions = {0.1, 0.9}; @@ -308,7 +308,7 @@ namespace serif::composition { * @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 +320,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 symbols = {"H", "O"}; * Composition comp; @@ -340,7 +340,7 @@ namespace serif::composition { * @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 +353,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 symbols = {"H", "O"}; * std::vector mass_fractions = {0.1, 0.9}; From 1a102368680de5eb7d0d777c6e95df4b7c441c06 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 16 Jun 2025 13:15:42 -0400 Subject: [PATCH 2/8] docs(docs): rebuilt docs --- docs/html/annotated.html | 2 +- docs/html/doxygen_crawl.html | 2 -- docs/html/navtreedata.js | 10 +++++----- docs/html/navtreeindex0.js | 10 +++++----- docs/latex/annotated.tex | 2 +- docs/latex/refman.tex | 2 -- 6 files changed, 12 insertions(+), 16 deletions(-) diff --git a/docs/html/annotated.html b/docs/html/annotated.html index 7aad003..79a63f8 100644 --- a/docs/html/annotated.html +++ b/docs/html/annotated.html @@ -106,7 +106,7 @@ $(function(){initNavTree('annotated.html',''); initResizable(true); });
[detail level 1234]
- + diff --git a/docs/html/doxygen_crawl.html b/docs/html/doxygen_crawl.html index 2caf3b1..f6a2ef5 100644 --- a/docs/html/doxygen_crawl.html +++ b/docs/html/doxygen_crawl.html @@ -10,12 +10,10 @@ - - diff --git a/docs/html/navtreedata.js b/docs/html/navtreedata.js index 39dd5ae..3af9c40 100644 --- a/docs/html/navtreedata.js +++ b/docs/html/navtreedata.js @@ -64,11 +64,11 @@ var NAVTREE = var NAVTREEINDEX = [ "4_d_s_t_a_r_types_8h.html", -"classserif_1_1polytrope_1_1_g_m_r_e_s_inverter.html#a698ccfa71e87e381b6c0a4d7a081a776", -"classserif_1_1polytrope_1_1poly_m_f_e_m_utils_1_1_nonlinear_power_integrator.html#ac636b3bdb45524d65d2a3aa1b6e43c7b", -"dir_daf60af48849bf958020a18d83ad56ea.html", -"namespaceserif_1_1polytrope_1_1poly_m_f_e_m_utils.html#ac43f5dda296efc47fbdbd351f2f4bf00", -"structserif_1_1eos_1_1helmholtz_1_1_h_e_l_m_table.html" +"classserif_1_1polytrope_1_1_g_m_r_e_s_inverter.html#a78c5540e756ad46201097cc83d90d356", +"classserif_1_1probe_1_1_log_manager.html", +"dir_e296cd0c311faef12afe540130b9e3be.html", +"namespaceserif_1_1polytrope_1_1polycoeff.html", +"structserif_1_1eos_1_1helmholtz_1_1_h_e_l_m_table.html#a0333b97d0f0498b86718cc20cb812e2c" ]; var SYNCONMSG = 'click to disable panel synchronization'; diff --git a/docs/html/navtreeindex0.js b/docs/html/navtreeindex0.js index b94e174..93f336f 100644 --- a/docs/html/navtreeindex0.js +++ b/docs/html/navtreeindex0.js @@ -2,12 +2,10 @@ var NAVTREEINDEX0 = { "4_d_s_t_a_r_types_8h.html":[4,0,1,11,0,0], "4_d_s_t_a_r_types_8h_source.html":[4,0,1,11,0,0], -"_2_users_2tboudreaux_2_programming_2_s_e_ri_f_2src_2composition_2public_2composition_8h-example.html":[5,0], -"_2_users_2tboudreaux_2_programming_2_s_e_ri_f_2src_2network_2public_2approx8_8h-example.html":[5,2], -"_2_users_2tboudreaux_2_programming_2_s_e_ri_f_2src_2polytrope_2solver_2public_2poly_solver_8h-example.html":[5,3], +"_2_users_2tboudreaux_2_programming_2_s_e_ri_f_2src_2network_2public_2approx8_8h-example.html":[5,0], +"_2_users_2tboudreaux_2_programming_2_s_e_ri_f_2src_2polytrope_2solver_2public_2poly_solver_8h-example.html":[5,1], "____init_____8py.html":[4,0,1,9,6,0], "____init_____8py_source.html":[4,0,1,9,6,0], -"_constructing-example.html":[5,1], "_e_o_sio_8cpp.html":[4,0,1,3,0,0], "_e_o_sio_8cpp_source.html":[4,0,1,3,0,0], "_e_o_sio_8h.html":[4,0,1,3,1,0], @@ -249,5 +247,7 @@ var NAVTREEINDEX0 = "classserif_1_1polytrope_1_1_g_m_r_e_s_inverter.html#a26e52db2e259ceb0efc74a188a0626df":[3,0,0,5,3,0], "classserif_1_1polytrope_1_1_g_m_r_e_s_inverter.html#a56eb2c60f01e499c905ccf1f9c1766dc":[2,0,3,6,4,4], "classserif_1_1polytrope_1_1_g_m_r_e_s_inverter.html#a56eb2c60f01e499c905ccf1f9c1766dc":[3,0,0,5,3,4], -"classserif_1_1polytrope_1_1_g_m_r_e_s_inverter.html#a698ccfa71e87e381b6c0a4d7a081a776":[2,0,3,6,4,2] +"classserif_1_1polytrope_1_1_g_m_r_e_s_inverter.html#a698ccfa71e87e381b6c0a4d7a081a776":[2,0,3,6,4,2], +"classserif_1_1polytrope_1_1_g_m_r_e_s_inverter.html#a698ccfa71e87e381b6c0a4d7a081a776":[3,0,0,5,3,2], +"classserif_1_1polytrope_1_1_g_m_r_e_s_inverter.html#a78c5540e756ad46201097cc83d90d356":[2,0,3,6,4,1] }; diff --git a/docs/latex/annotated.tex b/docs/latex/annotated.tex index a74b130..e15013a 100644 --- a/docs/latex/annotated.tex +++ b/docs/latex/annotated.tex @@ -2,7 +2,7 @@ Here are the classes, structs, unions and interfaces with brief descriptions\+:\begin{DoxyCompactList} \item\contentsline{section}{\mbox{\hyperlink{classserif_1_1network_1_1approx8_1_1_approx8_network}{serif\+::network\+::approx8\+::\+Approx8\+Network}} \\*Class for the Approx8 nuclear reaction network }{\pageref{classserif_1_1network_1_1approx8_1_1_approx8_network}}{} \item\contentsline{section}{\mbox{\hyperlink{classapprox8_test}{approx8\+Test}} }{\pageref{classapprox8_test}}{} -\item\contentsline{section}{\mbox{\hyperlink{classserif_1_1composition_1_1_composition}{serif\+::composition\+::\+Composition}} }{\pageref{classserif_1_1composition_1_1_composition}}{} +\item\contentsline{section}{\mbox{\hyperlink{classserif_1_1composition_1_1_composition}{serif\+::composition\+::\+Composition}} \\*Manages the composition of elements }{\pageref{classserif_1_1composition_1_1_composition}}{} \item\contentsline{section}{\mbox{\hyperlink{structserif_1_1composition_1_1_composition_entry}{serif\+::composition\+::\+Composition\+Entry}} \\*Represents an entry in the composition with a symbol and mass fraction }{\pageref{structserif_1_1composition_1_1_composition_entry}}{} \item\contentsline{section}{\mbox{\hyperlink{classcomposition_test}{composition\+Test}} \\*Test suite for the composition class }{\pageref{classcomposition_test}}{} \item\contentsline{section}{\mbox{\hyperlink{classconfig_test}{config\+Test}} \\*Test suite for the Config class }{\pageref{classconfig_test}}{} diff --git a/docs/latex/refman.tex b/docs/latex/refman.tex index 9e63bb4..b837d2c 100644 --- a/docs/latex/refman.tex +++ b/docs/latex/refman.tex @@ -447,8 +447,6 @@ \input{resource_manager_test_8cpp} \input{resource_manager_test_8cpp_source} \chapter{Examples} -\input{_2_users_2tboudreaux_2_programming_2_s_e_ri_f_2src_2composition_2public_2composition_8h-example} -\input{_constructing-example} \input{_2_users_2tboudreaux_2_programming_2_s_e_ri_f_2src_2network_2public_2approx8_8h-example} \input{_2_users_2tboudreaux_2_programming_2_s_e_ri_f_2src_2polytrope_2solver_2public_2poly_solver_8h-example} %--- End generated contents --- From cf25c54cda14a0bb8575719f388b2f6b62d37de1 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 16 Jun 2025 15:00:33 -0400 Subject: [PATCH 3/8] 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 | 1018 ++++++++++++----------- src/composition/public/composition.h | 38 +- src/eos/meson.build | 7 +- src/eos/private/EOS.cpp | 64 ++ src/eos/private/EOSio.cpp | 21 +- src/eos/private/helm.cpp | 4 +- src/eos/public/EOS.h | 230 +++++ src/eos/public/EOSio.h | 39 +- src/eos/public/helm.h | 14 +- src/python/eos/bindings.cpp | 94 +-- 10 files changed, 967 insertions(+), 562 deletions(-) create mode 100644 src/eos/private/EOS.cpp create mode 100644 src/eos/public/EOS.h diff --git a/src/composition/private/composition.cpp b/src/composition/private/composition.cpp index 6af4a88..9543836 100644 --- a/src/composition/private/composition.cpp +++ b/src/composition/private/composition.cpp @@ -25,6 +25,7 @@ #include #include #include +#include #include @@ -32,539 +33,600 @@ namespace serif::composition { -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) { - 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) {} - -void CompositionEntry::setSpecies(const std::string& symbol) { - if (m_initialized) { - throw std::runtime_error("Composition entry is already initialized."); + CompositionEntry::CompositionEntry() : m_symbol("H-1"), m_isotope(chemSpecies::species.at("H-1")), + m_initialized(false) { } - if (chemSpecies::species.count(symbol) == 0) { - throw std::runtime_error("Invalid symbol."); + + CompositionEntry::CompositionEntry(const std::string& symbol, const bool massFracMode) : m_symbol(symbol), m_isotope(chemSpecies::species.at(symbol)), m_massFracMode(massFracMode) { + setSpecies(symbol); } - m_symbol = symbol; - m_isotope = chemSpecies::species.at(symbol); - m_initialized = true; -} -std::string CompositionEntry::symbol() const { - return m_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) {} -double CompositionEntry::mass_fraction() const { - if (!m_massFracMode) { - throw std::runtime_error("Composition entry is in number fraction mode."); + void CompositionEntry::setSpecies(const std::string& symbol) { + if (m_initialized) { + throw std::runtime_error("Composition entry is already initialized."); + } + if (chemSpecies::species.count(symbol) == 0) { + throw std::runtime_error("Invalid symbol."); + } + m_symbol = symbol; + m_isotope = chemSpecies::species.at(symbol); + m_initialized = true; } - return m_massFraction; -} -double CompositionEntry::mass_fraction(double meanMolarMass) const { - if (m_massFracMode) { + std::string CompositionEntry::symbol() const { + return m_symbol; + } + + double CompositionEntry::mass_fraction() const { + if (!m_massFracMode) { + throw std::runtime_error("Composition entry is in number fraction mode."); + } return m_massFraction; } - return m_relAbundance / meanMolarMass; -} - -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 { - if (m_massFracMode) { - return m_relAbundance / totalMoles; - } - return m_numberFraction; -} - -double CompositionEntry::rel_abundance() const { - return m_relAbundance; -} - -chemSpecies::Species CompositionEntry::isotope() const { - return m_isotope; -} - -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) { - 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) { - if (m_massFracMode) { - return false; - } - m_massFracMode = true; - m_massFraction = m_relAbundance / meanParticleMass; - return true; -} - -bool CompositionEntry::setNumberFracMode(double specificNumberDensity) { - if (!m_massFracMode) { - return false; - } - m_massFracMode = false; - m_numberFraction = m_relAbundance / specificNumberDensity; - return true; -} - -bool CompositionEntry::getMassFracMode() const { - return m_massFracMode; -} - -Composition::Composition(const std::vector& symbols) { - for (const auto& symbol : symbols) { - registerSymbol(symbol); - } -} - -Composition::Composition(const std::set& symbols) { - for (const auto& symbol : symbols) { - registerSymbol(symbol); - } -} - -Composition::Composition(const std::vector& symbols, const std::vector& 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."); - } - - validateComposition(fractions); - - for (const auto &symbol : symbols) { - registerSymbol(symbol); - } - - for (size_t i = 0; i < symbols.size(); ++i) { + double CompositionEntry::mass_fraction(double meanMolarMass) const { if (m_massFracMode) { - setMassFraction(symbols[i], fractions[i]); + return m_massFraction; + } + return m_relAbundance / meanMolarMass; + } + + + 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 { + if (m_massFracMode) { + return m_relAbundance / totalMoles; + } + return m_numberFraction; + } + + double CompositionEntry::rel_abundance() const { + return m_relAbundance; + } + + chemSpecies::Species CompositionEntry::isotope() const { + return m_isotope; + } + + 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) { + 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) { + if (m_massFracMode) { + return false; + } + m_massFracMode = true; + m_massFraction = m_relAbundance / meanParticleMass; + return true; + } + + bool CompositionEntry::setNumberFracMode(double specificNumberDensity) { + if (!m_massFracMode) { + return false; + } + m_massFracMode = false; + m_numberFraction = m_relAbundance / specificNumberDensity; + return true; + } + + bool CompositionEntry::getMassFracMode() const { + return m_massFracMode; + } + + Composition::Composition(const std::vector& symbols) { + for (const auto& symbol : symbols) { + registerSymbol(symbol); + } + } + + Composition::Composition(const std::set& symbols) { + for (const auto& symbol : symbols) { + registerSymbol(symbol); + } + } + + Composition::Composition(const std::vector& symbols, const std::vector& 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."); + } + + validateComposition(fractions); + + for (const auto &symbol : symbols) { + registerSymbol(symbol); + } + + for (size_t i = 0; i < symbols.size(); ++i) { + if (m_massFracMode) { + setMassFraction(symbols[i], fractions[i]); + } else { + setNumberFraction(symbols[i], fractions[i]); + } + } + 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); + throw std::runtime_error("Invalid symbol."); + } + + // If no symbols have been registered allow mode to be set + if (m_registeredSymbols.size() == 0) { + m_massFracMode = massFracMode; } else { - setNumberFraction(symbols[i], fractions[i]); + if (m_massFracMode != massFracMode) { + LOG_ERROR(m_logger, "Composition is in mass fraction mode. Cannot register symbol in number fraction mode."); + throw std::runtime_error("Composition is in mass fraction mode. Cannot register symbol in number fraction mode."); + } } - } - finalize(); -} -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."); + if (m_registeredSymbols.find(symbol) != m_registeredSymbols.end()) { + LOG_WARNING(m_logger, "Symbol {} is already registered.", symbol); + return; + } + + m_registeredSymbols.insert(symbol); + CompositionEntry entry(symbol, m_massFracMode); + m_compositions[symbol] = entry; + LOG_INFO(m_logger, "Registered symbol: {}", symbol); } - // If no symbols have been registered allow mode to be set - if (m_registeredSymbols.size() == 0) { - 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."); - throw std::runtime_error("Composition is in mass fraction mode. Cannot register symbol in number fraction mode."); + void Composition::registerSymbol(const std::vector& symbols, bool massFracMode) { + for (const auto& symbol : symbols) { + registerSymbol(symbol, massFracMode); } } - if (m_registeredSymbols.find(symbol) != m_registeredSymbols.end()) { - LOG_WARNING(m_logger, "Symbol {} is already registered.", symbol); - return; + std::set Composition::getRegisteredSymbols() const { + return m_registeredSymbols; } - m_registeredSymbols.insert(symbol); - CompositionEntry entry(symbol, m_massFracMode); - m_compositions[symbol] = entry; - LOG_INFO(m_logger, "Registered symbol: {}", symbol); -} - -void Composition::registerSymbol(const std::vector& symbols, bool massFracMode) { - for (const auto& symbol : symbols) { - registerSymbol(symbol, massFracMode); - } -} - -std::set Composition::getRegisteredSymbols() const { - return m_registeredSymbols; -} - -void Composition::validateComposition(const std::vector& fractions) const { - if (!isValidComposition(fractions)) { - LOG_ERROR(m_logger, "Invalid composition."); - throw std::runtime_error("Invalid composition."); - } -} - -bool Composition::isValidComposition(const std::vector& fractions) const { - double sum = 0.0; - for (const auto& fraction : fractions) { - sum += fraction; - } - if (sum < 0.999999 || sum > 1.000001) { - LOG_ERROR(m_logger, "The sum of fractions must be equal to 1."); - return false; + void Composition::validateComposition(const std::vector& fractions) const { + if (!isValidComposition(fractions)) { + LOG_ERROR(m_logger, "Invalid composition."); + throw std::runtime_error("Invalid composition."); + } } - return true; -} - -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) { - 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."); - } - - if (!m_massFracMode) { - LOG_ERROR(m_logger, "Composition is in number fraction mode."); - throw std::runtime_error("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 std::runtime_error("Mass fraction must be between 0 and 1."); - } - - m_finalized = false; - double old_mass_fraction = m_compositions.at(symbol).mass_fraction(); - m_compositions.at(symbol).setMassFraction(mass_fraction); - - return old_mass_fraction; -} - -std::vector Composition::setMassFraction(const std::vector& symbols, const std::vector& 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."); - } - - std::vector old_mass_fractions; - old_mass_fractions.reserve(symbols.size()); - for (size_t i = 0; i < symbols.size(); ++i) { - 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) { - 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."); - } - - if (m_massFracMode) { - LOG_ERROR(m_logger, "Composition is in mass fraction mode."); - throw std::runtime_error("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 std::runtime_error("Number fraction must be between 0 and 1."); - } - - m_finalized = false; - double old_number_fraction = m_compositions.at(symbol).number_fraction(); - m_compositions.at(symbol).setNumberFraction(number_fraction); - - return old_number_fraction; -} - -std::vector Composition::setNumberFraction(const std::vector& symbols, const std::vector& 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."); - } - - std::vector old_number_fractions; - old_number_fractions.reserve(symbols.size()); - for (size_t i = 0; i < symbols.size(); ++i) { - old_number_fractions.push_back(setNumberFraction(symbols[i], number_fractions[i])); - } - return old_number_fractions; -} - -bool Composition::finalize(bool norm) { - bool finalized = false; - if (m_massFracMode) { - finalized = finalizeMassFracMode(norm); - } else { - finalized = finalizeNumberFracMode(norm); - } - if (finalized) { - m_finalized = true; - } - return finalized; -} - -bool Composition::finalizeMassFracMode(bool norm) { - std::vector mass_fractions; - mass_fractions.reserve(m_compositions.size()); - for (const auto& [_, entry] : m_compositions) { - mass_fractions.push_back(entry.mass_fraction()); - } - if (norm) { + bool Composition::isValidComposition(const std::vector& fractions) const { double sum = 0.0; - for (const auto& mass_fraction : mass_fractions) { - sum += mass_fraction; + for (const auto& fraction : fractions) { + sum += fraction; } - for (int i = 0; i < mass_fractions.size(); ++i) { - mass_fractions[i] /= sum; - } - for (auto& [symbol, entry] : m_compositions) { - setMassFraction(symbol, entry.mass_fraction() / sum); + if (sum < 0.999999 || sum > 1.000001) { + LOG_ERROR(m_logger, "The sum of fractions must be equal to 1."); + return false; } + + return true; } - try { - validateComposition(mass_fractions); - } catch (const std::runtime_error& e) { - double massSum = 0.0; - for (const auto& [_, entry] : m_compositions) { - massSum += entry.mass_fraction(); + + 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) { + 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."); } - LOG_ERROR(m_logger, "Composition is invalid (Total mass {}).", massSum); + + if (!m_massFracMode) { + LOG_ERROR(m_logger, "Composition is in number fraction mode."); + throw std::runtime_error("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 std::runtime_error("Mass fraction must be between 0 and 1."); + } + m_finalized = false; - return false; - } - for (const auto& [_, entry] : m_compositions) { - m_specificNumberDensity += entry.rel_abundance(); - } - m_meanParticleMass = 1.0/m_specificNumberDensity; - return true; -} + double old_mass_fraction = m_compositions.at(symbol).mass_fraction(); + m_compositions.at(symbol).setMassFraction(mass_fraction); -bool Composition::finalizeNumberFracMode(bool norm) { - std::vector number_fractions; - number_fractions.reserve(m_compositions.size()); - for (const auto& [_, entry] : m_compositions) { - number_fractions.push_back(entry.number_fraction()); + return old_mass_fraction; } - if (norm) { - double sum = 0.0; - for (const auto& number_fraction : number_fractions) { - sum += number_fraction; + + std::vector Composition::setMassFraction(const std::vector& symbols, const std::vector& 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."); } - for (auto& [symbol, entry] : m_compositions) { - setNumberFraction(symbol, entry.number_fraction() / sum); + + std::vector old_mass_fractions; + old_mass_fractions.reserve(symbols.size()); + for (size_t i = 0; i < symbols.size(); ++i) { + old_mass_fractions.push_back(setMassFraction(symbols[i], mass_fractions[i])); } + return old_mass_fractions; } - try { - validateComposition(number_fractions); - } catch (const std::runtime_error& e) { - double numberSum = 0.0; - for (const auto& [_, entry] : m_compositions) { - numberSum += entry.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."); } - LOG_ERROR(m_logger, "Composition is invalid (Total number {}).", numberSum); + + if (m_massFracMode) { + LOG_ERROR(m_logger, "Composition is in mass fraction mode."); + throw std::runtime_error("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 std::runtime_error("Number fraction must be between 0 and 1."); + } + m_finalized = false; - return false; - } - for (const auto& [_, entry] : m_compositions) { - m_meanParticleMass += entry.rel_abundance(); - } - m_specificNumberDensity = 1.0/m_meanParticleMass; - return true; -} + double old_number_fraction = m_compositions.at(symbol).number_fraction(); + m_compositions.at(symbol).setNumberFraction(number_fraction); -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())."); + return old_number_fraction; } - if (fraction < 0.0 || fraction > 1.0) { - LOG_ERROR(m_logger, "Fraction must be between 0 and 1."); - throw std::runtime_error("Fraction must be between 0 and 1."); + std::vector Composition::setNumberFraction(const std::vector& symbols, const std::vector& 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."); + } + + std::vector old_number_fractions; + old_number_fractions.reserve(symbols.size()); + for (size_t i = 0; i < symbols.size(); ++i) { + old_number_fractions.push_back(setNumberFraction(symbols[i], number_fractions[i])); + } + return old_number_fractions; } - std::set mixedSymbols = other.getRegisteredSymbols(); - // Get the union of the two sets - mixedSymbols.insert(m_registeredSymbols.begin(), m_registeredSymbols.end()); - - Composition mixedComposition(mixedSymbols); - for (const auto& symbol : mixedSymbols) { - double thisMassFrac, otherMassFrac = 0.0; - - thisMassFrac = hasSymbol(symbol) ? getMassFraction(symbol) : 0.0; - otherMassFrac = other.hasSymbol(symbol) ? other.getMassFraction(symbol) : 0.0; - - double massFraction = fraction * thisMassFrac + otherMassFrac * (1-fraction); - mixedComposition.setMassFraction(symbol, massFraction); - } - mixedComposition.finalize(); - return mixedComposition; -} - -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())."); - } - if (m_compositions.count(symbol) == 0) { - LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); - throw std::runtime_error("Symbol is not in the composition."); - } - if (m_massFracMode) { - return m_compositions.at(symbol).mass_fraction(); - } else { - return m_compositions.at(symbol).mass_fraction(m_meanParticleMass); - } -} - -std::unordered_map Composition::getMassFraction() const { - std::unordered_map 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 { - if (!m_finalized) { - 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) { - LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); - throw std::runtime_error("Symbol is not in the composition."); - } - if (!m_massFracMode) { - return m_compositions.at(symbol).number_fraction(); - } else { - return m_compositions.at(symbol).number_fraction(m_specificNumberDensity); - } -} - -std::unordered_map Composition::getNumberFraction() const { - std::unordered_map number_fractions; - for (const auto& [symbol, entry] : m_compositions) { - number_fractions[symbol] = getNumberFraction(symbol); - } - return number_fractions; -} - -std::pair 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())."); - } - if (m_compositions.count(symbol) == 0) { - LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); - throw std::runtime_error("Symbol is not in the composition."); - } - return {m_compositions.at(symbol), {m_specificNumberDensity, m_meanParticleMass}}; -} - -std::pair, 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& symbols, std::string method) const { - std::array methods = {"norm", "none"}; - - if (std::find(methods.begin(), methods.end(), method) == methods.end()) { - std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'."; - LOG_ERROR(m_logger, "Invalid method: {}. Valid methods are norm and none.", method); - throw std::runtime_error(errorMessage); + bool Composition::finalize(bool norm) { + bool finalized = false; + if (m_massFracMode) { + finalized = finalizeMassFracMode(norm); + } else { + finalized = finalizeNumberFracMode(norm); + } + if (finalized) { + m_finalized = true; + } + return finalized; } - Composition subsetComposition; - for (const auto& symbol : symbols) { + bool Composition::finalizeMassFracMode(bool norm) { + std::vector mass_fractions; + mass_fractions.reserve(m_compositions.size()); + for (const auto& [_, entry] : m_compositions) { + mass_fractions.push_back(entry.mass_fraction()); + } + if (norm) { + double sum = 0.0; + for (const auto& mass_fraction : mass_fractions) { + sum += mass_fraction; + } + for (int i = 0; i < mass_fractions.size(); ++i) { + mass_fractions[i] /= sum; + } + for (auto& [symbol, entry] : m_compositions) { + setMassFraction(symbol, entry.mass_fraction() / sum); + } + } + try { + validateComposition(mass_fractions); + } catch (const std::runtime_error& e) { + double massSum = 0.0; + for (const auto& [_, entry] : m_compositions) { + massSum += entry.mass_fraction(); + } + LOG_ERROR(m_logger, "Composition is invalid (Total mass {}).", massSum); + m_finalized = false; + return false; + } + for (const auto& [_, entry] : m_compositions) { + m_specificNumberDensity += entry.rel_abundance(); + } + m_meanParticleMass = 1.0/m_specificNumberDensity; + return true; + } + + bool Composition::finalizeNumberFracMode(bool norm) { + std::vector number_fractions; + number_fractions.reserve(m_compositions.size()); + for (const auto& [_, entry] : m_compositions) { + number_fractions.push_back(entry.number_fraction()); + } + if (norm) { + double sum = 0.0; + for (const auto& number_fraction : number_fractions) { + sum += number_fraction; + } + for (auto& [symbol, entry] : m_compositions) { + setNumberFraction(symbol, entry.number_fraction() / sum); + } + } + try { + validateComposition(number_fractions); + } catch (const std::runtime_error& e) { + double numberSum = 0.0; + for (const auto& [_, entry] : m_compositions) { + numberSum += entry.number_fraction(); + } + LOG_ERROR(m_logger, "Composition is invalid (Total number {}).", numberSum); + m_finalized = false; + return false; + } + for (const auto& [_, entry] : m_compositions) { + m_meanParticleMass += entry.rel_abundance(); + } + m_specificNumberDensity = 1.0/m_meanParticleMass; + return true; + } + + 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())."); + } + + if (fraction < 0.0 || fraction > 1.0) { + LOG_ERROR(m_logger, "Fraction must be between 0 and 1."); + throw std::runtime_error("Fraction must be between 0 and 1."); + } + + std::set mixedSymbols = other.getRegisteredSymbols(); + // Get the union of the two sets + mixedSymbols.insert(m_registeredSymbols.begin(), m_registeredSymbols.end()); + + Composition mixedComposition(mixedSymbols); + for (const auto& symbol : mixedSymbols) { + double thisMassFrac, otherMassFrac = 0.0; + + thisMassFrac = hasSymbol(symbol) ? getMassFraction(symbol) : 0.0; + otherMassFrac = other.hasSymbol(symbol) ? other.getMassFraction(symbol) : 0.0; + + double massFraction = fraction * thisMassFrac + otherMassFrac * (1-fraction); + mixedComposition.setMassFraction(symbol, massFraction); + } + mixedComposition.finalize(); + return mixedComposition; + } + + 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())."); + } if (m_compositions.count(symbol) == 0) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); + } + if (m_massFracMode) { + return m_compositions.at(symbol).mass_fraction(); } else { - subsetComposition.registerSymbol(symbol); - } - subsetComposition.setMassFraction(symbol, m_compositions.at(symbol).mass_fraction()); - } - if (method == "norm") { - bool isNorm = subsetComposition.finalize(true); - if (!isNorm) { - LOG_ERROR(m_logger, "Subset composition is invalid."); - throw std::runtime_error("Subset composition is invalid."); + return m_compositions.at(symbol).mass_fraction(m_meanParticleMass); } } - return subsetComposition; -} -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."); + std::unordered_map Composition::getMassFraction() const { + std::unordered_map mass_fractions; + for (const auto& [symbol, entry] : m_compositions) { + mass_fractions[symbol] = getMassFraction(symbol); + } + return mass_fractions; } - bool okay = true; - for (auto& [_, entry] : m_compositions) { - if (massFracMode) { - okay = entry.setMassFracMode(m_meanParticleMass); + + 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())."); + } + if (m_compositions.count(symbol) == 0) { + LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); + throw std::runtime_error("Symbol is not in the composition."); + } + if (!m_massFracMode) { + return m_compositions.at(symbol).number_fraction(); } else { - okay = entry.setNumberFracMode(m_specificNumberDensity); - } - if (!okay) { - LOG_ERROR(m_logger, "Composition mode could not be set."); - throw std::runtime_error("Composition mode could not be set due to an unknown error."); + return m_compositions.at(symbol).number_fraction(m_specificNumberDensity); } } - m_massFracMode = massFracMode; -} -bool Composition::hasSymbol(const std::string& symbol) const { - return m_compositions.count(symbol) > 0; -} - -/// OVERLOADS - -Composition Composition::operator+(const Composition& other) const { - return mix(other, 0.5); -} - -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) { - os << "<" << entry.m_symbol << " : m_frac = " << entry.mass_fraction() << ">"; - return os; -} - -std::ostream& operator<<(std::ostream& os, const Composition& composition) { - os << "Composition: \n"; - for (const auto& [symbol, entry] : composition.m_compositions) { - os << entry << "\n"; + std::unordered_map Composition::getNumberFraction() const { + std::unordered_map number_fractions; + for (const auto& [symbol, entry] : m_compositions) { + number_fractions[symbol] = getNumberFraction(symbol); + } + return number_fractions; + } + + std::pair 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())."); + } + if (m_compositions.count(symbol) == 0) { + LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); + throw std::runtime_error("Symbol is not in the composition."); + } + return {m_compositions.at(symbol), {m_specificNumberDensity, m_meanParticleMass}}; + } + + std::pair, 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}}; + } + + 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& symbols, std::string method) const { + std::array methods = {"norm", "none"}; + + if (std::find(methods.begin(), methods.end(), method) == methods.end()) { + std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'."; + 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) { + LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); + throw std::runtime_error("Symbol is not in the composition."); + } else { + subsetComposition.registerSymbol(symbol); + } + subsetComposition.setMassFraction(symbol, m_compositions.at(symbol).mass_fraction()); + } + if (method == "norm") { + bool isNorm = subsetComposition.finalize(true); + if (!isNorm) { + LOG_ERROR(m_logger, "Subset composition is invalid."); + throw std::runtime_error("Subset composition is invalid."); + } + } + return subsetComposition; + } + + 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."); + } + + bool okay = true; + for (auto& [_, entry] : m_compositions) { + if (massFracMode) { + okay = entry.setMassFracMode(m_meanParticleMass); + } else { + okay = entry.setNumberFracMode(m_specificNumberDensity); + } + if (!okay) { + LOG_ERROR(m_logger, "Composition mode could not be set."); + throw std::runtime_error("Composition mode could not be set due to an unknown error."); + } + } + m_massFracMode = massFracMode; + } + + bool Composition::hasSymbol(const std::string& symbol) const { + return m_compositions.count(symbol) > 0; + } + + /// OVERLOADS + + Composition Composition::operator+(const Composition& other) const { + return mix(other, 0.5); + } + + 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) { + os << "<" << entry.m_symbol << " : m_frac = " << entry.mass_fraction() << ">"; + return os; + } + + 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; } - return os; -} } // namespace serif::composition \ No newline at end of file diff --git a/src/composition/public/composition.h b/src/composition/public/composition.h index e61a5c6..457f96d 100644 --- a/src/composition/public/composition.h +++ b/src/composition/public/composition.h @@ -277,7 +277,7 @@ namespace serif::composition { * Composition comp(symbols); * @endcode */ - Composition(const std::vector& symbols); + explicit Composition(const std::vector& 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& symbols); + explicit Composition(const std::set& symbols); /** * @brief Constructs a Composition with the given symbols and mass fractions. @@ -304,6 +304,14 @@ namespace serif::composition { */ Composition(const std::vector& symbols, const std::vector& 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 getRegisteredSymbols() const; + [[nodiscard]] std::set 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 getMassFraction() const; + [[nodiscard]] std::unordered_map 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 getNumberFraction() const; + [[nodiscard]] std::unordered_map 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 getComposition(const std::string& symbol) const; + [[nodiscard]] std::pair 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, GlobalComposition> getComposition() const; + [[nodiscard]] std::pair, 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. diff --git a/src/eos/meson.build b/src/eos/meson.build index 00e285c..43d66e9 100644 --- a/src/eos/meson.build +++ b/src/eos/meson.build @@ -1,15 +1,18 @@ # Define the library eos_sources = files( 'private/helm.cpp', - 'private/EOSio.cpp' + 'private/EOSio.cpp', + 'private/EOS.cpp' ) eos_headers = files( 'public/helm.h', - 'public/EOSio.h' + 'public/EOSio.h', + 'public/EOS.h' ) dependencies = [ + composition_dep, const_dep, quill_dep, probe_dep, diff --git a/src/eos/private/EOS.cpp b/src/eos/private/EOS.cpp new file mode 100644 index 0000000..ba3593d --- /dev/null +++ b/src/eos/private/EOS.cpp @@ -0,0 +1,64 @@ +#include "EOS.h" +#include "EOSio.h" +#include "helm.h" +#include + +namespace serif::eos { + EOS::EOS(const EOSio& reader) : m_reader(reader) {} + EOS::EOS(const std::string& filename, const EOSFormat format) : m_reader(EOSio(filename, format)) {} + + EOSOutput EOS::get(const EOSInput& in) { + EOSOutput output; + if (getFormat() == EOSFormat::HELM) { + helmholtz::HELMEOSInput q; + q.T = in.temperature; // Temperature in K + q.rho = in.density; // Density in g/cm^3 + + q.abar = in.composition.getMeanParticleMass(); // Mean atomic mass in g + q.zbar = in.composition.getMeanAtomicNumber(); // Mean atomic number (dimensionless) + + helmholtz::HELMEOSOutput tempOutput; + tempOutput = helmholtz::get_helm_EOS(q, *std::get>(m_reader.getTable())); + + output.electronFraction = tempOutput.ye; + output.electronChemicalPotential = tempOutput.etaele; + output.neutronExcessFraction = tempOutput.xnefer; + + // --- Pressure Variables --- + output.pressure.total = tempOutput.ptot; + output.pressure.gas = tempOutput.pgas; + output.pressure.radiation = tempOutput.ptot; + output.pressure.dDensity = tempOutput.dpresdd; + output.pressure.dTemperature = tempOutput.dpresdt; + output.pressure.dMeanAtomicMassNumber = tempOutput.dpresda; + output.pressure.dMeanAtomicNumber = tempOutput.dpresdz; + + // --- Energy Variables --- + output.energy.total = tempOutput.etot; + output.energy.gas = tempOutput.egas; + output.energy.radiation = tempOutput.erad; + output.energy.dDensity = tempOutput.denerdd; + output.energy.dTemperature = tempOutput.denerdt; + output.energy.dMeanAtomicMassNumber = tempOutput.denerda; + output.energy.dMeanAtomicNumber = tempOutput.denerdz; + + // --- Entropy Variables --- + output.entropy.total = tempOutput.stot; + output.entropy.gas = tempOutput.sgas; + output.entropy.radiation = tempOutput.srad; + output.entropy.dDensity = tempOutput.dentrdd; + output.entropy.dTemperature = tempOutput.dentrdt; + output.entropy.dMeanAtomicMassNumber = tempOutput.dentrda; + output.entropy.dMeanAtomicNumber = tempOutput.dentrdz; + } + return output; + } + + EOSFormat EOS::getFormat() const { + return m_reader.getFormat(); + } + + const EOSio& EOS::getReader() const { + return m_reader; + } +} \ No newline at end of file diff --git a/src/eos/private/EOSio.cpp b/src/eos/private/EOSio.cpp index 03a6b43..6945e1a 100644 --- a/src/eos/private/EOSio.cpp +++ b/src/eos/private/EOSio.cpp @@ -28,25 +28,32 @@ #include namespace serif::eos { - EOSio::EOSio(const std::string &filename) : m_filename(filename) { + EOSio::EOSio(const std::string &filename, const EOSFormat format) : m_filename(filename), m_format(format){ load(); } - std::string EOSio::getFormat() const { + EOSio::EOSio(const EOSio &other) { + m_filename = other.m_filename; + m_format = other.m_format; + load(); + } + + EOSFormat EOSio::getFormat() const { return m_format; } + std::string EOSio::getFormatName() const { + return FormatStringLookup.at(m_format); + } + + EOSTable& EOSio::getTable() { return m_table; } void EOSio::load() { - // Load the EOS table from the file - // For now, just set the format to HELM - - m_format = "helm"; - if (m_format == "helm") { + if (m_format == EOSFormat::HELM) { loadHelm(); } } diff --git a/src/eos/private/helm.cpp b/src/eos/private/helm.cpp index 1c18b05..e2df9de 100644 --- a/src/eos/private/helm.cpp +++ b/src/eos/private/helm.cpp @@ -240,7 +240,7 @@ namespace serif::eos::helmholtz { ion, radiation, electron-positron and Coulomb interaction and returns the calculated quantities in the input ***/ - serif::eos::helmholtz::EOS get_helm_EOS(serif::eos::helmholtz::EOSInput &q, const serif::eos::helmholtz::HELMTable &table) { + serif::eos::helmholtz::HELMEOSOutput get_helm_EOS(serif::eos::helmholtz::HELMEOSInput &q, const serif::eos::helmholtz::HELMTable &table) { serif::config::Config& config = serif::config::Config::getInstance(); auto logFile = config.get("EOS:Helm:LogFile", "log"); serif::probe::LogManager& logManager = serif::probe::LogManager::getInstance(); @@ -829,7 +829,7 @@ namespace serif::eos::helmholtz { double csound = clight * sqrt(gamma1/z); // package in q: - serif::eos::helmholtz::EOS eos; + serif::eos::helmholtz::HELMEOSOutput eos; eos.ptot = ptot; eos.etot = etot; eos.stot = stot; eos.pgas = pgas; eos.egas = egas; eos.sgas = sgas; eos.prad = prad; eos.erad = erad; eos.srad = srad; diff --git a/src/eos/public/EOS.h b/src/eos/public/EOS.h new file mode 100644 index 0000000..b2ecd22 --- /dev/null +++ b/src/eos/public/EOS.h @@ -0,0 +1,230 @@ +#pragma once + +#include "EOSio.h" +#include "helm.h" +#include +#include "composition.h" + +namespace serif::eos { + + /** + * @brief Input parameters for an EOS calculation. + * + * This struct holds the necessary physical conditions (composition, density, temperature) + * required to query the Equation of State. + */ + struct EOSInput { + serif::composition::Composition composition; ///< The composition of the system. + double density; ///< The density of the system in cgs (g/cm^3). + double temperature; ///< The temperature of the system in cgs (K). + }; + + /** + * @brief Represents a thermodynamic parameter and its derivatives. + * + * This struct stores a specific thermodynamic quantity (e.g., pressure, energy, entropy), + * its breakdown into gas and radiation components, and its partial derivatives + * with respect to density, temperature, mean atomic mass number, and mean atomic number. + * All values are in cgs units unless otherwise specified. + */ + struct EOSParameter { + double total; ///< Total value of the parameter (gas + radiation) (cgs). + double gas; ///< Gas contribution to the parameter (cgs). + double radiation; ///< Radiation contribution to the parameter (cgs). + + double dDensity; ///< Derivative of the total parameter with respect to density (cgs units / (g/cm^3)). + double dTemperature; ///< Derivative of the total parameter with respect to temperature (cgs units / K). + double dMeanAtomicMassNumber; ///< Derivative of the total parameter with respect to mean atomic mass number (Abar) (cgs units / (g/mol)). + double dMeanAtomicNumber; ///< Derivative of the total parameter with respect to mean atomic number (Zbar) (cgs units / dimensionless). + + std::string name; ///< Name of the parameter (e.g., "Pressure", "Energy", "Entropy"). + }; + + /** + * @brief Output from an EOS calculation. + * + * This struct contains various thermodynamic quantities and their derivatives + * calculated by the EOS for a given set of input conditions. + * It includes fundamental properties like electron fraction and chemical potential, + * as well as detailed breakdowns of pressure, energy, and entropy. + * Additionally, it provides methods to calculate derived quantities like + * susceptibilities, sound speed, adiabatic gradients, and specific heats. + */ + struct EOSOutput { + double electronFraction{}; ///< Electron fraction (ye), dimensionless. + double electronChemicalPotential{}; ///< Electron chemical potential (eta_e) in cgs (erg/g). + double neutronExcessFraction{}; ///< Neutron excess fraction (xnefer), dimensionless. + + EOSParameter pressure; ///< Pressure output data, including total, gas, radiation, and derivatives. + EOSParameter energy; ///< Internal energy output data, including total, gas, radiation, and derivatives. + EOSParameter entropy; ///< Entropy output data, including total, gas, radiation, and derivatives. + + /** + * @brief Calculates the temperature susceptibility (chi_T). + * @return Temperature susceptibility, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, chi_T = (d ln P / d ln T)_rho. + */ + double chiTemperature(); + /** + * @brief Calculates the density susceptibility (chi_rho). + * @return Density susceptibility, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, chi_rho = (d ln P / d ln rho)_T. + */ + double chiRho(); + /** + * @brief Calculates the adiabatic sound speed. + * @return Sound speed in cgs (cm/s). + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, c_s^2 = gamma1 * P / rho. + */ + double soundSpeed(); + /** + * @brief Calculates the adiabatic temperature gradient (nabla_ad). + * @return Adiabatic gradient, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, nabla_ad = (P * chi_T) / (rho * T * c_p * chi_rho). + */ + double adiabaticGradient(); + /** + * @brief Calculates the first adiabatic index (Gamma1). + * @return First adiabatic index, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, Gamma1 = (d ln P / d ln rho)_S. + */ + double gamma1(); + /** + * @brief Calculates the second adiabatic index (Gamma2). + * @return Second adiabatic index, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, Gamma2 / (Gamma2 - 1) = (d ln P / d ln T)_S. + */ + double gamma2(); + /** + * @brief Calculates the third adiabatic index (Gamma3). + * @return Third adiabatic index, dimensionless. + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, Gamma3 - 1 = (d ln T / d ln rho)_S. + */ + double gamma3(); + /** + * @brief Calculates the specific heat capacity at constant volume (c_v). + * @return Specific heat capacity at constant volume in cgs (erg/K/g). + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, c_v = (dE/dT)_rho. + */ + double specificHeatCapacityAtConstantVolume(); + /** + * @brief Calculates the specific heat capacity at constant pressure (c_p). + * @return Specific heat capacity at constant pressure in cgs (erg/K/g). + * @note Placeholder: Actual calculation needs to be implemented based on available EOS derivatives. + * Typically, c_p = c_v + (T / rho^2) * ( (d P / d T)_rho^2 / (d P / d rho)_T ). + */ + double specificHeatCapacityAtConstantPressure(); + + /** + * @brief Returns the format of the EOS data used to generate this output. + * @return The EOSFormat enum value (currently only EOSFormat::HELM). + */ + EOSFormat EOSFormat() const; + }; + + /** + * @class EOS + * @brief Main class for accessing Equation of State data. + * + * This class provides an interface to an underlying EOS table (e.g., Helmholtz EOS). + * It handles loading the EOS data and provides a method to retrieve thermodynamic + * properties for given physical conditions. + * + * @example + * @code + * #include "EOS.h" + * #include "composition.h" // For serif::composition::Composition + * #include + * + * int main() { + * try { + * // Initialize EOS from a Helmholtz table file + * serif::eos::EOS helmEOS("path/to/helm_table.dat", serif::eos::EOSFormat::HELM); + * + * // Define input conditions + * serif::eos::EOSInput input; + * input.density = 1.0e6; // g/cm^3 + * input.temperature = 1.0e7; // K + * // Assuming a simple composition (e.g., pure Helium-4 for demonstration) + * // In a real scenario, initialize Composition properly. + * // For example, if Composition has a constructor like: + * // Composition(const std::map, double>& mass_fractions); + * // std::map, double> he4_mass_fraction = {{{2, 4}, 1.0}}; + * // input.composition = serif::composition::Composition(he4_mass_fraction); + * // For now, let's assume Composition can be default constructed or set up simply: + * input.composition.addSpecies(2, 4, 1.0); // Z=2, A=4 (He-4), mass fraction 1.0 + * + * // Get EOS output + * serif::eos::EOSOutput output = helmEOS.get(input); + * + * // Access results + * std::cout << "Pressure (total): " << output.pressure.total << " dyne/cm^2" << std::endl; + * std::cout << "Energy (total): " << output.energy.total << " erg/g" << std::endl; + * std::cout << "Entropy (total): " << output.entropy.total << " erg/K/g" << std::endl; + * std::cout << "Electron fraction: " << output.electronFraction << std::endl; + * + * // Example of accessing derivatives + * std::cout << "dP/dRho: " << output.pressure.dDensity << std::endl; + * std::cout << "dE/dT: " << output.energy.dTemperature << std::endl; + * + * } catch (const std::exception& e) { + * std::cerr << "An error occurred: " << e.what() << std::endl; + * return 1; + * } + * return 0; + * } + * @endcode + */ + class EOS { + public: + /** + * @brief Constructs an EOS object by loading data from a file. + * @param filename The path to the EOS data file. + * @param format The format of the EOS data file (e.g., EOSFormat::HELM). + * @throw std::runtime_error If the file cannot be opened or read, or if the format is unsupported. + */ + explicit EOS(const std::string& filename, EOSFormat format=EOSFormat::HELM); + /** + * @brief Constructs an EOS object from an existing EOSio reader. + * @param reader An EOSio object that has already loaded the EOS data. + */ + explicit EOS(const EOSio& reader); + /** + * @brief Default destructor. + */ + ~EOS() = default; + + /** + * @brief Retrieves thermodynamic properties for the given input conditions. + * @param in An EOSInput struct containing the density, temperature, and composition. + * @return An EOSOutput struct containing the calculated thermodynamic properties. + * @throw std::runtime_error If the underlying EOS calculation fails (e.g., out of table bounds for Helmholtz). + * + * This method queries the loaded EOS table (e.g., Helmholtz) using the provided + * density, temperature, and composition (mean atomic mass Abar, mean atomic number Zbar). + * It populates and returns an EOSOutput struct with various thermodynamic quantities + * such as pressure, energy, entropy, their derivatives, electron fraction, etc. + */ + [[nodiscard]] EOSOutput get(const EOSInput& in); + /** + * @brief Gets the format of the loaded EOS data. + * @return The EOSFormat enum value. + */ + [[nodiscard]] EOSFormat getFormat() const; + /** + * @brief Gets a constant reference to the internal EOSio reader. + * @return A const reference to the EOSio object. + */ + [[nodiscard]] const EOSio& getReader() const; + private: + EOSio m_reader; ///< The EOS I/O handler responsible for reading and storing EOS table data. + }; +} diff --git a/src/eos/public/EOSio.h b/src/eos/public/EOSio.h index 1882c07..fff0b64 100644 --- a/src/eos/public/EOSio.h +++ b/src/eos/public/EOSio.h @@ -22,15 +22,20 @@ #include #include #include +#include #include "helm.h" namespace serif::eos { + enum EOSFormat { + HELM, ///< Helmholtz EOS format. + }; + static inline std::unordered_map FormatStringLookup = { + {HELM, "Helmholtz"} + }; // EOS table format includes - using EOSTable = std::variant< - std::unique_ptr - >; + using EOSTable = std::variant>; /** * @class EOSio @@ -41,16 +46,19 @@ namespace serif::eos { * * Example usage: * @code - * EosIO eosIO("path/to/file"); - * std::string format = eosIO.getFormat(); - * EOSTable& table = eosIO.getTable(); + * EOSio eosReader("path/to/file"); + * std::string format = eosReader.getFormatName(); + * EOSTable& table = eosReader.getTable(); * @endcode + * + * @note The default format used for reading tables is HELM + * @note Currently only the HELM format is implemented */ class EOSio { private: std::string m_filename; ///< The filename of the EOS table. bool m_loaded = false; ///< Flag indicating if the table is loaded. - std::string m_format; ///< The format of the EOS table. + EOSFormat m_format; EOSTable m_table; ///< The EOS table data. /** @@ -66,8 +74,15 @@ namespace serif::eos { /** * @brief Constructs an EosIO object with the given filename. * @param filename The filename of the EOS table. + * @param format The EOS file format (currently only HELM) */ - explicit EOSio(const std::string &filename); + explicit EOSio(const std::string &filename, EOSFormat format = EOSFormat::HELM); + + /** + * @brief Explicit copy constructor + * @param other The EOSio to be copied + */ + EOSio(const EOSio& other); /** * @brief Default destructor. @@ -75,10 +90,12 @@ namespace serif::eos { ~EOSio() = default; /** - * @brief Gets the format of the EOS table. + * @brief Gets the format name (as a string) of the EOS table. * @return The format of the EOS table as a string. */ - [[nodiscard]] std::string getFormat() const; + [[nodiscard]] std::string getFormatName() const; + + [[nodiscard]] EOSFormat getFormat() const; /** * @brief Gets the EOS table. @@ -87,6 +104,8 @@ namespace serif::eos { [[nodiscard]] EOSTable& getTable(); [[nodiscard]] std::string getFilename() const { return m_filename; } + + bool isLoaded() const { return m_loaded; } }; } diff --git a/src/eos/public/helm.h b/src/eos/public/helm.h index 501db49..a20ba49 100644 --- a/src/eos/public/helm.h +++ b/src/eos/public/helm.h @@ -169,18 +169,18 @@ namespace serif::eos::helmholtz { }; /** - * @struct EOSInput + * @struct HELMEOSInput * @brief Structure to hold the input parameters for the EOS calculation. */ - struct EOSInput + struct HELMEOSInput { double T; ///< Temperature. double rho; ///< Density. double abar; ///< Mean atomic mass. double zbar; ///< Mean atomic number. - friend std::ostream& operator<<(std::ostream& os, const helmholtz::EOSInput& eosInput) { - os << "EOSInput Data:\n"; + friend std::ostream& operator<<(std::ostream& os, const helmholtz::HELMEOSInput& eosInput) { + os << "HELMEOSInput Data:\n"; os << " Temperature: " << eosInput.T << "\n"; os << " Density: " << eosInput.rho << "\n"; os << " Mean Atomic Mass: " << eosInput.abar << "\n"; @@ -194,7 +194,7 @@ namespace serif::eos::helmholtz { * @struct EOS * @brief Structure to hold the output parameters and derivatives of the EOS calculation. */ - struct EOS + struct HELMEOSOutput { // output double ye, etaele, xnefer; // @@ -212,7 +212,7 @@ namespace serif::eos::helmholtz { double gamma1, gamma2, gamma3, cV, cP; // derived quantities double dse, dpe, dsp; // Maxwell relations - friend std::ostream& operator<<(std::ostream& os, const helmholtz::EOS& eos) { + friend std::ostream& operator<<(std::ostream& os, const helmholtz::HELMEOSOutput& eos) { os << "EOS Data:\n" << std::setw(20) << std::left; os << " Electron Fraction: " << std::format("{0:24.16e}",eos.ye) << "\n"; os << " Electron Chemical Potential: " << std::format("{0:24.16e}",eos.etaele) << "\n"; @@ -379,7 +379,7 @@ namespace serif::eos::helmholtz { * @param table HELMTable structure containing the table data. * @return EOS structure containing the calculated quantities. */ - EOS get_helm_EOS(EOSInput &q, const HELMTable &table); + HELMEOSOutput get_helm_EOS(HELMEOSInput &q, const HELMTable &table); } diff --git a/src/python/eos/bindings.cpp b/src/python/eos/bindings.cpp index e8c6314..ca4b87c 100644 --- a/src/python/eos/bindings.cpp +++ b/src/python/eos/bindings.cpp @@ -34,7 +34,7 @@ void register_eos_bindings(pybind11::module &eos_submodule) { }, py::return_value_policy::reference_internal, // IMPORTANT: Keep this policy! "Get the EOS table data.") .def("__repr__", [](const serif::eos::EOSio &eos) { - return ""; + return ""; }); py::class_(eos_submodule, "EOSTable"); @@ -88,64 +88,64 @@ void register_eos_bindings(pybind11::module &eos_submodule) { ); }, py::return_value_policy::reference_internal); // Keep parent 'table' alive - py::class_(eos_submodule, "EOS") + py::class_(eos_submodule, "EOS") .def(py::init<>()) - .def_readonly("ye", &serif::eos::helmholtz::EOS::ye) - .def_readonly("etaele", &serif::eos::helmholtz::EOS::etaele) - .def_readonly("xnefer", &serif::eos::helmholtz::EOS::xnefer) + .def_readonly("ye", &serif::eos::helmholtz::HELMEOSOutput::ye) + .def_readonly("etaele", &serif::eos::helmholtz::HELMEOSOutput::etaele) + .def_readonly("xnefer", &serif::eos::helmholtz::HELMEOSOutput::xnefer) - .def_readonly("ptot", &serif::eos::helmholtz::EOS::ptot) - .def_readonly("pgas", &serif::eos::helmholtz::EOS::pgas) - .def_readonly("prad", &serif::eos::helmholtz::EOS::prad) + .def_readonly("ptot", &serif::eos::helmholtz::HELMEOSOutput::ptot) + .def_readonly("pgas", &serif::eos::helmholtz::HELMEOSOutput::pgas) + .def_readonly("prad", &serif::eos::helmholtz::HELMEOSOutput::prad) - .def_readonly("etot", &serif::eos::helmholtz::EOS::etot) - .def_readonly("egas", &serif::eos::helmholtz::EOS::egas) - .def_readonly("erad", &serif::eos::helmholtz::EOS::erad) + .def_readonly("etot", &serif::eos::helmholtz::HELMEOSOutput::etot) + .def_readonly("egas", &serif::eos::helmholtz::HELMEOSOutput::egas) + .def_readonly("erad", &serif::eos::helmholtz::HELMEOSOutput::erad) - .def_readonly("stot", &serif::eos::helmholtz::EOS::stot) - .def_readonly("sgas", &serif::eos::helmholtz::EOS::sgas) - .def_readonly("srad", &serif::eos::helmholtz::EOS::srad) + .def_readonly("stot", &serif::eos::helmholtz::HELMEOSOutput::stot) + .def_readonly("sgas", &serif::eos::helmholtz::HELMEOSOutput::sgas) + .def_readonly("srad", &serif::eos::helmholtz::HELMEOSOutput::srad) - .def_readonly("dpresdd", &serif::eos::helmholtz::EOS::dpresdd) - .def_readonly("dpresdt", &serif::eos::helmholtz::EOS::dpresdt) - .def_readonly("dpresda", &serif::eos::helmholtz::EOS::dpresda) - .def_readonly("dpresdz", &serif::eos::helmholtz::EOS::dpresdz) - // TODO: Finish adding all the derivatives to the bound class - .def_readonly("dentrdd", &serif::eos::helmholtz::EOS::dentrdd) - .def_readonly("dentrdt", &serif::eos::helmholtz::EOS::dentrdt) - .def_readonly("dentrda", &serif::eos::helmholtz::EOS::dentrda) - .def_readonly("dentrdz", &serif::eos::helmholtz::EOS::dentrdz) + .def_readonly("dpresdd", &serif::eos::helmholtz::HELMEOSOutput::dpresdd) + .def_readonly("dpresdt", &serif::eos::helmholtz::HELMEOSOutput::dpresdt) + .def_readonly("dpresda", &serif::eos::helmholtz::HELMEOSOutput::dpresda) + .def_readonly("dpresdz", &serif::eos::helmholtz::HELMEOSOutput::dpresdz) - .def_readonly("denerdd", &serif::eos::helmholtz::EOS::denerdd) - .def_readonly("denerdt", &serif::eos::helmholtz::EOS::denerdt) - .def_readonly("denerda", &serif::eos::helmholtz::EOS::denerda) - .def_readonly("denerdz", &serif::eos::helmholtz::EOS::denerdz) + .def_readonly("dentrdd", &serif::eos::helmholtz::HELMEOSOutput::dentrdd) + .def_readonly("dentrdt", &serif::eos::helmholtz::HELMEOSOutput::dentrdt) + .def_readonly("dentrda", &serif::eos::helmholtz::HELMEOSOutput::dentrda) + .def_readonly("dentrdz", &serif::eos::helmholtz::HELMEOSOutput::dentrdz) - .def_readonly("chiT", &serif::eos::helmholtz::EOS::chiT) - .def_readonly("chiRho", &serif::eos::helmholtz::EOS::chiRho) - .def_readonly("csound", &serif::eos::helmholtz::EOS::csound) - .def_readonly("grad_ad", &serif::eos::helmholtz::EOS::grad_ad) - .def_readonly("gamma1", &serif::eos::helmholtz::EOS::gamma1) - .def_readonly("gamma2", &serif::eos::helmholtz::EOS::gamma2) - .def_readonly("gamma3", &serif::eos::helmholtz::EOS::gamma3) - .def_readonly("cV", &serif::eos::helmholtz::EOS::cV) - .def_readonly("cP", &serif::eos::helmholtz::EOS::cP) - .def_readonly("dse", &serif::eos::helmholtz::EOS::dse) - .def_readonly("dpe", &serif::eos::helmholtz::EOS::dpe) - .def_readonly("dsp", &serif::eos::helmholtz::EOS::dsp) + .def_readonly("denerdd", &serif::eos::helmholtz::HELMEOSOutput::denerdd) + .def_readonly("denerdt", &serif::eos::helmholtz::HELMEOSOutput::denerdt) + .def_readonly("denerda", &serif::eos::helmholtz::HELMEOSOutput::denerda) + .def_readonly("denerdz", &serif::eos::helmholtz::HELMEOSOutput::denerdz) - .def("__repr__", [](const serif::eos::helmholtz::EOS &eos) { + .def_readonly("chiT", &serif::eos::helmholtz::HELMEOSOutput::chiT) + .def_readonly("chiRho", &serif::eos::helmholtz::HELMEOSOutput::chiRho) + .def_readonly("csound", &serif::eos::helmholtz::HELMEOSOutput::csound) + .def_readonly("grad_ad", &serif::eos::helmholtz::HELMEOSOutput::grad_ad) + .def_readonly("gamma1", &serif::eos::helmholtz::HELMEOSOutput::gamma1) + .def_readonly("gamma2", &serif::eos::helmholtz::HELMEOSOutput::gamma2) + .def_readonly("gamma3", &serif::eos::helmholtz::HELMEOSOutput::gamma3) + .def_readonly("cV", &serif::eos::helmholtz::HELMEOSOutput::cV) + .def_readonly("cP", &serif::eos::helmholtz::HELMEOSOutput::cP) + .def_readonly("dse", &serif::eos::helmholtz::HELMEOSOutput::dse) + .def_readonly("dpe", &serif::eos::helmholtz::HELMEOSOutput::dpe) + .def_readonly("dsp", &serif::eos::helmholtz::HELMEOSOutput::dsp) + + .def("__repr__", [](const serif::eos::helmholtz::HELMEOSOutput &eos) { return ""; }); - py::class_(eos_submodule, "EOSInput") + py::class_(eos_submodule, "HELMEOSInput") .def(py::init<>()) - .def_readwrite("T", &serif::eos::helmholtz::EOSInput::T) - .def_readwrite("rho", &serif::eos::helmholtz::EOSInput::rho) - .def_readwrite("abar", &serif::eos::helmholtz::EOSInput::abar) - .def_readwrite("zbar", &serif::eos::helmholtz::EOSInput::zbar) - .def("__repr__", [](const serif::eos::helmholtz::EOSInput &input) { - return ""; From 9100af3fc562269439715293a68d2afebcd3c6ce Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Mon, 16 Jun 2025 15:01:05 -0400 Subject: [PATCH 4/8] test(tests/eos): added test for the new composition version of EOS note that this test is currently minimal and should be expanded --- tests/eos/eosTest.cpp | 28 ++++++++++++++++++++++++-- tests/eos/meson.build | 12 ++++++++++- tests/resource/resourceManagerTest.cpp | 2 +- 3 files changed, 38 insertions(+), 4 deletions(-) diff --git a/tests/eos/eosTest.cpp b/tests/eos/eosTest.cpp index b7a670b..a4c12ca 100644 --- a/tests/eos/eosTest.cpp +++ b/tests/eos/eosTest.cpp @@ -7,6 +7,8 @@ #include "helm.h" #include "resourceManager.h" #include "config.h" +#include "composition.h" +#include "EOS.h" /** * @file constTest.cpp @@ -39,7 +41,7 @@ TEST_F(eosTest, get_helm_EOS) { const int nel=3; double xmass[nel], aion[nel], zion[nel]; - serif::eos::helmholtz::EOSInput eos1; + serif::eos::helmholtz::HELMEOSInput eos1; xmass[0] = 0.75; aion[0] = 1.0; zion[0] = 1.0; xmass[1] = 0.23; aion[1] = 4.0; zion[1] = 2.0; @@ -61,7 +63,7 @@ TEST_F(eosTest, get_helm_EOS) { auto& eos = std::get>(rm.getResource("eos:helm")); auto& table = eos->getTable(); auto& helmTable = *std::get>(table); - serif::eos::helmholtz::EOS helmEos = get_helm_EOS(eos1, helmTable); + serif::eos::helmholtz::HELMEOSOutput helmEos = get_helm_EOS(eos1, helmTable); const double absErr = 1e-12; @@ -85,3 +87,25 @@ TEST_F(eosTest, get_helm_EOS) { EXPECT_NEAR( helmEos.dpe, 0, absErr); EXPECT_NEAR( helmEos.dsp, 0, absErr); } + +TEST_F(eosTest, eos_using_composition) { + serif::composition::Composition composition; + composition.registerSymbol({"H-1", "He-4", "C-12"}, true); + composition.setMassFraction("H-1", 0.75); + composition.setMassFraction("He-4", 0.23); + composition.setMassFraction("C-12", 0.02); + composition.finalize(); + + serif::resource::ResourceManager& rm = serif::resource::ResourceManager::getInstance(); + auto& EOSio = std::get>(rm.getResource("eos:helm")); + + serif::eos::EOS EOS(*EOSio); + + serif::eos::EOSInput eosInput; + eosInput.temperature = 1.0e8; // Temperature in K + eosInput.density = 1.0e6; // Density in g/cm^3 + eosInput.composition = composition; // Set the composition + + serif::eos::EOSOutput eosOutput; + EXPECT_NO_THROW(eosOutput = EOS.get(eosInput)); +} diff --git a/tests/eos/meson.build b/tests/eos/meson.build index d8154e6..61b86e3 100644 --- a/tests/eos/meson.build +++ b/tests/eos/meson.build @@ -3,6 +3,16 @@ test_sources = [ 'eosTest.cpp', ] +dependencies = [ + gtest_dep, + eos_dep, + gtest_main, + resourceManager_dep, + config_dep, + composition_dep, +] + + foreach test_file : test_sources exe_name = test_file.split('.')[0] message('Building test: ' + exe_name) @@ -11,7 +21,7 @@ foreach test_file : test_sources test_exe = executable( exe_name, test_file, - dependencies: [gtest_dep, eos_dep, gtest_main, resourceManager_dep, config_dep], + dependencies: dependencies, install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly ) diff --git a/tests/resource/resourceManagerTest.cpp b/tests/resource/resourceManagerTest.cpp index 2936653..4238689 100644 --- a/tests/resource/resourceManagerTest.cpp +++ b/tests/resource/resourceManagerTest.cpp @@ -48,7 +48,7 @@ TEST_F(resourceManagerTest, getResource) { const serif::resource::types::Resource &r = rm.getResource(name); // BREAKPOINT(); const auto &eos = std::get>(r); - EXPECT_EQ("helm", eos->getFormat()); + EXPECT_EQ("Helmholtz", eos->getFormatName()); serif::eos::EOSTable &table = eos->getTable(); // -- Extract the Helm table from the EOSTable From 3961c745e3660212243d3124952ca2fd798d2c9f Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 17 Jun 2025 08:12:41 -0400 Subject: [PATCH 5/8] fix(eos): fixed calculation of mean atomic number --- src/composition/private/composition.cpp | 100 ++++++++++++++++-------- src/composition/public/composition.h | 26 +++++- src/eos/private/EOS.cpp | 2 + src/eos/public/EOS.h | 42 +++++++++- 4 files changed, 134 insertions(+), 36 deletions(-) diff --git a/src/composition/private/composition.cpp b/src/composition/private/composition.cpp index 9543836..98a3668 100644 --- a/src/composition/private/composition.cpp +++ b/src/composition/private/composition.cpp @@ -255,12 +255,12 @@ namespace serif::composition { return true; } - bool Composition::isValidSymbol(const std::string& symbol) const { - return chemSpecies::species.count(symbol) > 0; + bool Composition::isValidSymbol(const std::string& symbol) { + return chemSpecies::species.contains(symbol); } double Composition::setMassFraction(const std::string& symbol, const double& mass_fraction) { - if (m_registeredSymbols.find(symbol) == m_registeredSymbols.end()) { + if (!m_registeredSymbols.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol); throw std::runtime_error("Symbol is not registered."); } @@ -276,7 +276,7 @@ namespace serif::composition { } m_finalized = false; - double old_mass_fraction = m_compositions.at(symbol).mass_fraction(); + const double old_mass_fraction = m_compositions.at(symbol).mass_fraction(); m_compositions.at(symbol).setMassFraction(mass_fraction); return old_mass_fraction; @@ -432,9 +432,9 @@ namespace serif::composition { Composition mixedComposition(mixedSymbols); for (const auto& symbol : mixedSymbols) { - double thisMassFrac, otherMassFrac = 0.0; + double otherMassFrac = 0.0; - thisMassFrac = hasSymbol(symbol) ? getMassFraction(symbol) : 0.0; + const double thisMassFrac = hasSymbol(symbol) ? getMassFraction(symbol) : 0.0; otherMassFrac = other.hasSymbol(symbol) ? other.getMassFraction(symbol) : 0.0; double massFraction = fraction * thisMassFrac + otherMassFrac * (1-fraction); @@ -449,7 +449,7 @@ namespace serif::composition { LOG_ERROR(m_logger, "Composition has not been finalized."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); } - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } @@ -462,7 +462,7 @@ namespace serif::composition { std::unordered_map Composition::getMassFraction() const { std::unordered_map mass_fractions; - for (const auto& [symbol, entry] : m_compositions) { + for (const auto &symbol: m_compositions | std::views::keys) { mass_fractions[symbol] = getMassFraction(symbol); } return mass_fractions; @@ -474,7 +474,7 @@ namespace serif::composition { LOG_ERROR(m_logger, "Composition has not been finalized."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); } - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } @@ -487,7 +487,7 @@ namespace serif::composition { std::unordered_map Composition::getNumberFraction() const { std::unordered_map number_fractions; - for (const auto& [symbol, entry] : m_compositions) { + for (const auto &symbol: m_compositions | std::views::keys) { number_fractions[symbol] = getNumberFraction(symbol); } return number_fractions; @@ -498,7 +498,7 @@ namespace serif::composition { LOG_ERROR(m_logger, "Composition has not been finalized."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); } - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } @@ -527,40 +527,29 @@ namespace serif::composition { throw std::runtime_error("Composition not finalized. Cannot retrieve mean atomic mass number."); } - double mean_A = 0.0; + double zSum = 0.0; // Loop through all registered species in the composition. for (const auto &val: m_compositions | std::views::values) { - const CompositionEntry& entry = val; - const chemSpecies::Species& species = entry.isotope(); - - const double mass_fraction = entry.mass_fraction(); - const double particle_mass_g = species.mass(); - const int mass_number = species.a(); - - // Avoid division by zero, though a valid species should have a positive mass. - if (particle_mass_g > 0) { - // Calculate the number fraction for this species. - const double number_fraction = (mass_fraction / particle_mass_g) * m_meanParticleMass; - mean_A += number_fraction * mass_number; - } + zSum += (val.mass_fraction() * val.m_isotope.z())/val.m_isotope.a(); } + const double mean_A = m_meanParticleMass * zSum; return mean_A; } Composition Composition::subset(const std::vector& symbols, std::string method) const { - std::array methods = {"norm", "none"}; + const std::array methods = {"norm", "none"}; - if (std::find(methods.begin(), methods.end(), method) == methods.end()) { - std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'."; + if (std::ranges::find(methods, method) == methods.end()) { + const std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'."; LOG_ERROR(m_logger, "Invalid method: {}. Valid methods are norm and none.", method); throw std::runtime_error(errorMessage); } Composition subsetComposition; for (const auto& symbol : symbols) { - if (m_compositions.count(symbol) == 0) { + if (!m_compositions.contains(symbol)) { LOG_ERROR(m_logger, "Symbol {} is not in the composition.", symbol); throw std::runtime_error("Symbol is not in the composition."); } else { @@ -569,7 +558,7 @@ namespace serif::composition { subsetComposition.setMassFraction(symbol, m_compositions.at(symbol).mass_fraction()); } if (method == "norm") { - bool isNorm = subsetComposition.finalize(true); + const bool isNorm = subsetComposition.finalize(true); if (!isNorm) { LOG_ERROR(m_logger, "Subset composition is invalid."); throw std::runtime_error("Subset composition is invalid."); @@ -578,14 +567,14 @@ namespace serif::composition { return subsetComposition; } - void Composition::setCompositionMode(bool massFracMode) { + void Composition::setCompositionMode(const bool massFracMode) { if (!m_finalized) { LOG_ERROR(m_logger, "Composition has not been finalized. Mode cannot be set unless composition is finalized."); throw std::runtime_error("Composition has not been finalized (Consider running .finalize()). The mode cannot be set unless the composition is finalized."); } bool okay = true; - for (auto& [_, entry] : m_compositions) { + for (auto &entry: m_compositions | std::views::values) { if (massFracMode) { okay = entry.setMassFracMode(m_meanParticleMass); } else { @@ -599,6 +588,53 @@ namespace serif::composition { m_massFracMode = massFracMode; } + CanonicalComposition Composition::getCanonicalComposition(bool harsh) const { + if (!m_finalized) { + LOG_ERROR(m_logger, "Composition has not been finalized."); + throw std::runtime_error("Composition has not been finalized (Consider running .finalize())."); + } + CanonicalComposition canonicalComposition; + constexpr std::array canonicalH = { + "H-1", "H-2", "H-3", "H-4", "H-5", "H-6", "H-7" + }; + constexpr std::array canonicalHe = { + "He-3", "He-4", "He-5", "He-6", "He-7", "He-8", "He-9", "He-10" + }; + for (const auto& symbol : canonicalH) { + if (hasSymbol(symbol)) { + canonicalComposition.X += getMassFraction(symbol); + } + } + for (const auto& symbol : canonicalHe) { + if (hasSymbol(symbol)) { + canonicalComposition.Y += getMassFraction(symbol); + } + } + + for (const auto& symbol : getRegisteredSymbols()) { + const bool isHSymbol = std::ranges::find(canonicalH, symbol) != std::end(canonicalH); + const bool isHeSymbol = std::ranges::find(canonicalHe, symbol) != std::end(canonicalHe); + + if (isHSymbol || isHeSymbol) { + continue; // Skip canonical H and He symbols + } + + canonicalComposition.Z += getMassFraction(symbol); + } + + const double Z = 1.0 - (canonicalComposition.X + canonicalComposition.Y); + if (std::abs(Z - canonicalComposition.Z) > 1e-6) { + if (!harsh) { + LOG_WARNING(m_logger, "Validation composition Z (X-Y = {}) is different than canonical composition Z ({}) (∑a_i where a_i != H/He).", Z, canonicalComposition.Z); + } + else { + LOG_ERROR(m_logger, "Validation composition Z (X-Y = {}) is different than canonical composition Z ({}) (∑a_i where a_i != H/He).", Z, canonicalComposition.Z); + throw std::runtime_error("Validation composition Z (X-Y = " + std::to_string(Z) + ") is different than canonical composition Z (" + std::to_string(canonicalComposition.Z) + ") (∑a_i where a_i != H/He)."); + } + } + return canonicalComposition; + } + bool Composition::hasSymbol(const std::string& symbol) const { return m_compositions.count(symbol) > 0; } diff --git a/src/composition/public/composition.h b/src/composition/public/composition.h index 457f96d..861cc47 100644 --- a/src/composition/public/composition.h +++ b/src/composition/public/composition.h @@ -34,6 +34,20 @@ #include "atomicSpecies.h" namespace serif::composition { + struct CanonicalComposition { + double X = 0.0; ///< Mass fraction of Hydrogen. + double Y = 0.0; ///< Mass fraction of Helium. + double Z = 0.0; ///< Mass fraction of Metals. + + friend std::ostream& operator<<(std::ostream& os, const CanonicalComposition& composition) { + os << ""; + return os; + } + }; + /** * @brief Represents the global composition of a system. This tends to be used after finalize and is primarily for internal use. */ @@ -220,7 +234,7 @@ namespace serif::composition { * @param symbol The symbol to check. * @return True if the symbol is valid, false otherwise. */ - bool isValidSymbol(const std::string& symbol) const; + static bool isValidSymbol(const std::string& symbol); /** * @brief Checks if the given mass fractions are valid. @@ -466,6 +480,16 @@ namespace serif::composition { */ void setCompositionMode(bool massFracMode); + /** + * @brief Gets the current canonical composition (X, Y, Z). + * @param harsh If true, this will throw an error if X-Y != Z where Z is computed as the sum of all other elements. + * @return True if mass fraction mode, false if number fraction mode. + * + * @throws std::runtime_error if the composition is not finalized or if the canonical composition cannot be computed. + * @throws std::runtime_error if harsh is true and the canonical composition is not valid. + */ + [[nodiscard]] CanonicalComposition getCanonicalComposition(bool harsh=false) const; + /** * @brief Overloaded output stream operator for Composition. * @param os The output stream. diff --git a/src/eos/private/EOS.cpp b/src/eos/private/EOS.cpp index ba3593d..0ec82f9 100644 --- a/src/eos/private/EOS.cpp +++ b/src/eos/private/EOS.cpp @@ -17,6 +17,8 @@ namespace serif::eos { q.abar = in.composition.getMeanParticleMass(); // Mean atomic mass in g q.zbar = in.composition.getMeanAtomicNumber(); // Mean atomic number (dimensionless) + std::cout << "(EOS) abar: " << q.abar << ", zbar: " << q.zbar << std::endl; + helmholtz::HELMEOSOutput tempOutput; tempOutput = helmholtz::get_helm_EOS(q, *std::get>(m_reader.getTable())); diff --git a/src/eos/public/EOS.h b/src/eos/public/EOS.h index b2ecd22..5baf725 100644 --- a/src/eos/public/EOS.h +++ b/src/eos/public/EOS.h @@ -4,6 +4,7 @@ #include "helm.h" #include #include "composition.h" +#include namespace serif::eos { @@ -17,6 +18,13 @@ namespace serif::eos { serif::composition::Composition composition; ///< The composition of the system. double density; ///< The density of the system in cgs (g/cm^3). double temperature; ///< The temperature of the system in cgs (K). + friend std::ostream& operator<<(std::ostream& os, const EOSInput& input) { + os << ""; + return os; + } }; /** @@ -28,6 +36,12 @@ namespace serif::eos { * All values are in cgs units unless otherwise specified. */ struct EOSParameter { + explicit EOSParameter(std::string name_) + : total(), gas(), radiation(), + dDensity(), dTemperature(), + dMeanAtomicMassNumber(), + dMeanAtomicNumber(), + name(std::move(name_)) {} double total; ///< Total value of the parameter (gas + radiation) (cgs). double gas; ///< Gas contribution to the parameter (cgs). double radiation; ///< Radiation contribution to the parameter (cgs). @@ -38,6 +52,17 @@ namespace serif::eos { double dMeanAtomicNumber; ///< Derivative of the total parameter with respect to mean atomic number (Zbar) (cgs units / dimensionless). std::string name; ///< Name of the parameter (e.g., "Pressure", "Energy", "Entropy"). + + friend std::ostream& operator<<(std::ostream& os, const EOSParameter& param) { + os << std::setprecision(3) << ""; + return os; + } + + }; /** @@ -55,9 +80,9 @@ namespace serif::eos { double electronChemicalPotential{}; ///< Electron chemical potential (eta_e) in cgs (erg/g). double neutronExcessFraction{}; ///< Neutron excess fraction (xnefer), dimensionless. - EOSParameter pressure; ///< Pressure output data, including total, gas, radiation, and derivatives. - EOSParameter energy; ///< Internal energy output data, including total, gas, radiation, and derivatives. - EOSParameter entropy; ///< Entropy output data, including total, gas, radiation, and derivatives. + EOSParameter pressure {"pressure"}; ///< Pressure output data, including total, gas, radiation, and derivatives. + EOSParameter energy {"energy"}; ///< Internal energy output data, including total, gas, radiation, and derivatives. + EOSParameter entropy {"entropy"}; ///< Entropy output data, including total, gas, radiation, and derivatives. /** * @brief Calculates the temperature susceptibility (chi_T). @@ -128,6 +153,17 @@ namespace serif::eos { * @return The EOSFormat enum value (currently only EOSFormat::HELM). */ EOSFormat EOSFormat() const; + + friend std::ostream& operator<<(std::ostream& os, const EOSOutput& output) { + os << "EOSOutput:\n" + << "\tElectron Fraction: " << output.electronFraction << "\n" + << "\tElectron Chemical Potential: " << output.electronChemicalPotential << "\n" + << "\tNeutron Excess Fraction: " << output.neutronExcessFraction << "\n\t" + << output.pressure << "\n\t" + << output.energy << "\n\t" + << output.entropy; + return os; + } }; /** From bc548ed311e9913584f038219a9496b2ee8b4984 Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 17 Jun 2025 08:18:06 -0400 Subject: [PATCH 6/8] test(tests/eos): added test for composition module mode --- src/eos/private/EOS.cpp | 2 -- src/eos/public/EOS.h | 2 +- tests/eos/eosTest.cpp | 25 ++++++++++++++++++++++--- 3 files changed, 23 insertions(+), 6 deletions(-) diff --git a/src/eos/private/EOS.cpp b/src/eos/private/EOS.cpp index 0ec82f9..ba3593d 100644 --- a/src/eos/private/EOS.cpp +++ b/src/eos/private/EOS.cpp @@ -17,8 +17,6 @@ namespace serif::eos { q.abar = in.composition.getMeanParticleMass(); // Mean atomic mass in g q.zbar = in.composition.getMeanAtomicNumber(); // Mean atomic number (dimensionless) - std::cout << "(EOS) abar: " << q.abar << ", zbar: " << q.zbar << std::endl; - helmholtz::HELMEOSOutput tempOutput; tempOutput = helmholtz::get_helm_EOS(q, *std::get>(m_reader.getTable())); diff --git a/src/eos/public/EOS.h b/src/eos/public/EOS.h index 5baf725..207d746 100644 --- a/src/eos/public/EOS.h +++ b/src/eos/public/EOS.h @@ -54,7 +54,7 @@ namespace serif::eos { std::string name; ///< Name of the parameter (e.g., "Pressure", "Energy", "Entropy"). friend std::ostream& operator<<(std::ostream& os, const EOSParameter& param) { - os << std::setprecision(3) << ">(table); serif::eos::helmholtz::HELMEOSOutput helmEos = get_helm_EOS(eos1, helmTable); - const double absErr = 1e-12; + constexpr double absErr = 1e-12; //Check composition info EXPECT_NEAR( helmEos.ye, 8.75e-01, absErr); @@ -90,10 +90,26 @@ TEST_F(eosTest, get_helm_EOS) { TEST_F(eosTest, eos_using_composition) { serif::composition::Composition composition; - composition.registerSymbol({"H-1", "He-4", "C-12"}, true); + composition.registerSymbol({ + "H-1", + "He-4", + "C-12", + "O-16", + "Ne-20", + "Fe-56", + "N-14", + "Si-28", + "Mg-24" + }, true); composition.setMassFraction("H-1", 0.75); composition.setMassFraction("He-4", 0.23); - composition.setMassFraction("C-12", 0.02); + composition.setMassFraction("C-12", 0.0044); + composition.setMassFraction("O-16", 0.0096); + composition.setMassFraction("Ne-20", 0.002); + composition.setMassFraction("Fe-56", 0.0018); + composition.setMassFraction("N-14", 0.001); + composition.setMassFraction("Si-28", 0.0008); + composition.setMassFraction("Mg-24", 0.0004); composition.finalize(); serif::resource::ResourceManager& rm = serif::resource::ResourceManager::getInstance(); @@ -108,4 +124,7 @@ TEST_F(eosTest, eos_using_composition) { serif::eos::EOSOutput eosOutput; EXPECT_NO_THROW(eosOutput = EOS.get(eosInput)); + eosOutput = EOS.get(eosInput); + constexpr double absErr = 1e-8; + EXPECT_NEAR(eosOutput.pressure.total, 6.9548533046915791E+22, absErr); } From 3c14fd54bd1b22c1f3a5626fc5111b2d036a79bd Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 17 Jun 2025 08:19:50 -0400 Subject: [PATCH 7/8] test(tests/eos): improved const and constexpr correctness --- tests/eos/eosTest.cpp | 32 +++++++++++++++----------------- 1 file changed, 15 insertions(+), 17 deletions(-) diff --git a/tests/eos/eosTest.cpp b/tests/eos/eosTest.cpp index d6e39d0..d98a25c 100644 --- a/tests/eos/eosTest.cpp +++ b/tests/eos/eosTest.cpp @@ -1,5 +1,4 @@ #include -#include #include #include #include @@ -11,7 +10,7 @@ #include "EOS.h" /** - * @file constTest.cpp + * @file eosTest.cpp * @brief Unit tests for the const class. */ @@ -28,36 +27,35 @@ std::string TEST_CONFIG = std::string(getenv("MESON_SOURCE_ROOT")) + "/tests/tes TEST_F(eosTest, read_helm_table) { serif::config::Config::getInstance().loadConfig(TEST_CONFIG); - serif::resource::ResourceManager& rm = serif::resource::ResourceManager::getInstance(); + const serif::resource::ResourceManager& rm = serif::resource::ResourceManager::getInstance(); auto& eos = std::get>(rm.getResource("eos:helm")); - auto& table = eos->getTable(); - auto& helmTable = *std::get>(table); + const auto& table = eos->getTable(); + const auto& helmTable = *std::get>(table); std::stringstream ss; ss << helmTable; EXPECT_EQ(ss.str(), "HELMTable Data:\n imax: 541, jmax: 201\n Temperature Range: [1000, 1e+13]\n Density Range: [1e-12, 1e+15]\n"); } TEST_F(eosTest, get_helm_EOS) { - - const int nel=3; - double xmass[nel], aion[nel], zion[nel]; + constexpr int nel=3; + double xMass[nel], aIon[nel], zIon[nel]; serif::eos::helmholtz::HELMEOSInput eos1; - xmass[0] = 0.75; aion[0] = 1.0; zion[0] = 1.0; - xmass[1] = 0.23; aion[1] = 4.0; zion[1] = 2.0; - xmass[2] = 0.02; aion[2] = 12.0; zion[2] = 6.0; + xMass[0] = 0.75; aIon[0] = 1.0; zIon[0] = 1.0; + xMass[1] = 0.23; aIon[1] = 4.0; zIon[1] = 2.0; + xMass[2] = 0.02; aIon[2] = 12.0; zIon[2] = 6.0; eos1.T = 1.0e8; eos1.rho = 1.0e6; - double asum = 0.0; - double zsum = 0.0; + double aSum = 0.0; + double zSum = 0.0; for (int i=0; i>(rm.getResource("eos:helm")); From 0cfb020a8c81e6f65506a5c176990658eee529bb Mon Sep 17 00:00:00 2001 From: Emily Boudreaux Date: Tue, 17 Jun 2025 09:43:43 -0400 Subject: [PATCH 8/8] refactor(network): updated network and network::approx8 to use composition module This is a very basic wrapper implimentation currently. This is sufficient to lock the interface down so that other code can target it. However, internally there is just a "convert" function. Eventually we should rework the code itself to use the composition module more directly. --- assets/static/atomic/meson.build | 3 +- src/composition/meson.build | 2 +- src/composition/private/composition.cpp | 2 +- src/composition/public/composition.h | 7 +- src/network/meson.build | 4 +- src/network/private/approx8.cpp | 467 ++++++++++++------------ src/network/private/network.cpp | 42 ++- src/network/public/approx8.h | 92 ++--- src/network/public/network.h | 31 +- tests/network/approx8Test.cpp | 23 +- tests/network/meson.build | 2 +- 11 files changed, 370 insertions(+), 305 deletions(-) diff --git a/assets/static/atomic/meson.build b/assets/static/atomic/meson.build index 9e024aa..b5b79b7 100644 --- a/assets/static/atomic/meson.build +++ b/assets/static/atomic/meson.build @@ -1,3 +1,4 @@ species_weight_dep = declare_dependency( include_directories: include_directories('include'), -) \ No newline at end of file +) +message('✅ SERiF species_weight dependency declared') \ No newline at end of file diff --git a/src/composition/meson.build b/src/composition/meson.build index fe928d8..d97cbec 100644 --- a/src/composition/meson.build +++ b/src/composition/meson.build @@ -29,4 +29,4 @@ composition_dep = declare_dependency( ) # Make headers accessible -install_headers(composition_headers, subdir : '4DSSE/composition') \ No newline at end of file +install_headers(composition_headers, subdir : 'SERiF/composition') \ No newline at end of file diff --git a/src/composition/private/composition.cpp b/src/composition/private/composition.cpp index 98a3668..aad84dc 100644 --- a/src/composition/private/composition.cpp +++ b/src/composition/private/composition.cpp @@ -333,7 +333,7 @@ namespace serif::composition { return old_number_fractions; } - bool Composition::finalize(bool norm) { + bool Composition::finalize(const bool norm) { bool finalized = false; if (m_massFracMode) { finalized = finalizeMassFracMode(norm); diff --git a/src/composition/public/composition.h b/src/composition/public/composition.h index 861cc47..3dcf8a0 100644 --- a/src/composition/public/composition.h +++ b/src/composition/public/composition.h @@ -18,8 +18,7 @@ // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // // *********************************************************************** */ -#ifndef COMPOSITION_H -#define COMPOSITION_H +#pragma once #include #include @@ -507,6 +506,4 @@ namespace serif::composition { Composition operator+(const Composition& other) const; }; -}; // namespace serif::composition - -#endif // COMPOSITION_H \ No newline at end of file +}; // namespace serif::composition \ No newline at end of file diff --git a/src/network/meson.build b/src/network/meson.build index 41abf3f..97d410e 100644 --- a/src/network/meson.build +++ b/src/network/meson.build @@ -15,7 +15,9 @@ dependencies = [ quill_dep, mfem_dep, config_dep, - probe_dep + probe_dep, + species_weight_dep, + composition_dep, ] # Define the libnetwork library so it can be linked against by other parts of the build system diff --git a/src/network/private/approx8.cpp b/src/network/private/approx8.cpp index c485765..f46bb7e 100644 --- a/src/network/private/approx8.cpp +++ b/src/network/private/approx8.cpp @@ -31,7 +31,7 @@ #include "approx8.h" #include "network.h" -/* Nuclear reaction network in cgs units based on Frank Timmes' "aprox8". +/* Nuclear reaction network in cgs units based on Frank Timmes' "approx8". At this time it does neither screening nor neutrino losses. It includes the following 8 isotopes: @@ -78,7 +78,7 @@ namespace serif::network::approx8{ using namespace boost::numeric::odeint; //helper functions - // a function to multilpy two arrays and then sum the resulting elements: sum(a*b) + // a function to multiply two arrays and then sum the resulting elements: sum(a*b) double sum_product( const vec7 &a, const vec7 &b){ if (a.size() != b.size()) { throw std::runtime_error("Error: array size mismatch in sum_product"); @@ -96,8 +96,8 @@ namespace serif::network::approx8{ // this function returns an array of the T9 terms in that order, where T9 is the temperatures in GigaKelvin vec7 get_T9_array(const double &T) { vec7 arr; - double T9=1e-9*T; - double T913=pow(T9,1./3.); + const double T9=1e-9*T; + const double T913=pow(T9,1./3.); arr[0]=1; arr[1]=1/T9; @@ -117,125 +117,125 @@ namespace serif::network::approx8{ // p + p -> d; this, like some of the other rates, this is a composite of multiple fits double pp_rate(const vec7 &T9) { - vec7 a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170}; - vec7 a2 = { -4.364990e+1,-2.460640e-3,-2.750700,-4.248770e-1,1.598700e-2,-6.908750e-4,-2.076250e-1}; + constexpr vec7 a1 = {-34.78630, 0,-3.511930, 3.100860, -0.1983140, 1.262510e-2, -1.025170}; + constexpr vec7 a2 = { -4.364990e+1,-2.460640e-3,-2.750700,-4.248770e-1,1.598700e-2,-6.908750e-4,-2.076250e-1}; return rate_fit(T9,a1) + rate_fit(T9,a2); } // p + d -> he3 double dp_rate(const vec7 &T9) { - vec7 a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670}; - vec7 a2 = {8.935250, 0, -3.720800, 0.1986540, 0, 0, 0.3333330}; + constexpr vec7 a1 = {7.528980, 0, -3.720800, 0.8717820, 0, 0,-0.6666670}; + constexpr vec7 a2 = {8.935250, 0, -3.720800, 0.1986540, 0, 0, 0.3333330}; return rate_fit(T9,a1) + rate_fit(T9,a2); } // he3 + he3 -> he4 + 2p double he3he3_rate(const vec7 &T9){ - vec7 a = {2.477880e+01,0,-12.27700,-0.1036990,-6.499670e-02,1.681910e-02,-6.666670e-01}; + constexpr vec7 a = {2.477880e+01,0,-12.27700,-0.1036990,-6.499670e-02,1.681910e-02,-6.666670e-01}; return rate_fit(T9,a); } // he3(he3,2p)he4 double he3he4_rate(const vec7 &T9){ - vec7 a1 = {1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01}; - vec7 a2 = {1.770750e+01,0.000000e+00,-1.282710e+01,-3.812600e+00,9.422850e-02,-3.010180e-03,1.333330e+00}; + constexpr vec7 a1 = {1.560990e+01,0.000000e+00,-1.282710e+01,-3.082250e-02,-6.546850e-01,8.963310e-02,-6.666670e-01}; + constexpr vec7 a2 = {1.770750e+01,0.000000e+00,-1.282710e+01,-3.812600e+00,9.422850e-02,-3.010180e-03,1.333330e+00}; return rate_fit(T9,a1) + rate_fit(T9,a2); } // he4 + he4 + he4 -> c12 double triple_alpha_rate(const vec7 &T9){ - vec7 a1 = {-9.710520e-01,0.000000e+00,-3.706000e+01,2.934930e+01,-1.155070e+02,-1.000000e+01,-1.333330e+00}; - vec7 a2 = {-1.178840e+01,-1.024460e+00,-2.357000e+01,2.048860e+01,-1.298820e+01,-2.000000e+01,-2.166670e+00}; - vec7 a3 = {-2.435050e+01,-4.126560e+00,-1.349000e+01,2.142590e+01,-1.347690e+00,8.798160e-02,-1.316530e+01}; + constexpr vec7 a1 = {-9.710520e-01,0.000000e+00,-3.706000e+01,2.934930e+01,-1.155070e+02,-1.000000e+01,-1.333330e+00}; + constexpr vec7 a2 = {-1.178840e+01,-1.024460e+00,-2.357000e+01,2.048860e+01,-1.298820e+01,-2.000000e+01,-2.166670e+00}; + constexpr vec7 a3 = {-2.435050e+01,-4.126560e+00,-1.349000e+01,2.142590e+01,-1.347690e+00,8.798160e-02,-1.316530e+01}; return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3); } // c12 + p -> n13 double c12p_rate(const vec7 &T9){ - vec7 a1={1.714820e+01,0.000000e+00,-1.369200e+01,-2.308810e-01,4.443620e+00,-3.158980e+00,-6.666670e-01}; - vec7 a2={1.754280e+01,-3.778490e+00,-5.107350e+00,-2.241110e+00,1.488830e-01,0.000000e+00,-1.500000e+00}; + constexpr vec7 a1={1.714820e+01,0.000000e+00,-1.369200e+01,-2.308810e-01,4.443620e+00,-3.158980e+00,-6.666670e-01}; + constexpr vec7 a2={1.754280e+01,-3.778490e+00,-5.107350e+00,-2.241110e+00,1.488830e-01,0.000000e+00,-1.500000e+00}; return rate_fit(T9,a1) + rate_fit(T9,a2); } // c12 + he4 -> o16 double c12a_rate(const vec7 &T9){ - vec7 a1={6.965260e+01,-1.392540e+00,5.891280e+01,-1.482730e+02,9.083240e+00,-5.410410e-01,7.035540e+01}; - vec7 a2={2.546340e+02,-1.840970e+00,1.034110e+02,-4.205670e+02,6.408740e+01,-1.246240e+01,1.373030e+02}; + constexpr vec7 a1={6.965260e+01,-1.392540e+00,5.891280e+01,-1.482730e+02,9.083240e+00,-5.410410e-01,7.035540e+01}; + constexpr vec7 a2={2.546340e+02,-1.840970e+00,1.034110e+02,-4.205670e+02,6.408740e+01,-1.246240e+01,1.373030e+02}; return rate_fit(T9,a1) + rate_fit(T9,a2); } // n14(p,g)o15 - o15 + p -> c12 + he4 double n14p_rate(const vec7 &T9){ - vec7 a1={1.701000e+01,0.000000e+00,-1.519300e+01,-1.619540e-01,-7.521230e+00,-9.875650e-01,-6.666670e-01}; - vec7 a2={2.011690e+01,0.000000e+00,-1.519300e+01,-4.639750e+00,9.734580e+00,-9.550510e+00,3.333330e-01}; - vec7 a3={7.654440e+00,-2.998000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec7 a4={6.735780e+00,-4.891000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,6.820000e-02}; + constexpr vec7 a1={1.701000e+01,0.000000e+00,-1.519300e+01,-1.619540e-01,-7.521230e+00,-9.875650e-01,-6.666670e-01}; + constexpr vec7 a2={2.011690e+01,0.000000e+00,-1.519300e+01,-4.639750e+00,9.734580e+00,-9.550510e+00,3.333330e-01}; + constexpr vec7 a3={7.654440e+00,-2.998000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + constexpr vec7 a4={6.735780e+00,-4.891000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,6.820000e-02}; return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4); } // n14(a,g)f18 assumed to go on to ne20 double n14a_rate(const vec7 &T9){ - vec7 a1={2.153390e+01,0.000000e+00,-3.625040e+01,0.000000e+00,0.000000e+00,-5.000000e+00,-6.666670e-01}; - vec7 a2={1.968380e-01,-5.160340e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec7 a3={1.389950e+01,-1.096560e+01,-5.622700e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + constexpr vec7 a1={2.153390e+01,0.000000e+00,-3.625040e+01,0.000000e+00,0.000000e+00,-5.000000e+00,-6.666670e-01}; + constexpr vec7 a2={1.968380e-01,-5.160340e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + constexpr vec7 a3={1.389950e+01,-1.096560e+01,-5.622700e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3); } // n15(p,a)c12 (CNO I) double n15pa_rate(const vec7 &T9){ - vec7 a1 = {2.747640e+01,0.000000e+00,-1.525300e+01,1.593180e+00,2.447900e+00,-2.197080e+00,-6.666670e-01}; - vec7 a2 = {-4.873470e+00,-2.021170e+00,0.000000e+00,3.084970e+01,-8.504330e+00,-1.544260e+00,-1.500000e+00}; - vec7 a3 = {2.089720e+01,-7.406000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec7 a4 = {-6.575220e+00,-1.163800e+00,0.000000e+00,2.271050e+01,-2.907070e+00,2.057540e-01,-1.500000e+00}; + constexpr vec7 a1 = {2.747640e+01,0.000000e+00,-1.525300e+01,1.593180e+00,2.447900e+00,-2.197080e+00,-6.666670e-01}; + constexpr vec7 a2 = {-4.873470e+00,-2.021170e+00,0.000000e+00,3.084970e+01,-8.504330e+00,-1.544260e+00,-1.500000e+00}; + constexpr vec7 a3 = {2.089720e+01,-7.406000e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + constexpr vec7 a4 = {-6.575220e+00,-1.163800e+00,0.000000e+00,2.271050e+01,-2.907070e+00,2.057540e-01,-1.500000e+00}; return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4); } // n15(p,g)o16 (CNO II) double n15pg_rate(const vec7 &T9){ - vec7 a1 = {2.001760e+01,0.000000e+00,-1.524000e+01,3.349260e-01,4.590880e+00,-4.784680e+00,-6.666670e-01}; - vec7 a2 = {6.590560e+00,-2.923150e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec7 a3 = {1.454440e+01,-1.022950e+01,0.000000e+00,0.000000e+00,4.590370e-02,0.000000e+00,-1.500000e+00}; + constexpr vec7 a1 = {2.001760e+01,0.000000e+00,-1.524000e+01,3.349260e-01,4.590880e+00,-4.784680e+00,-6.666670e-01}; + constexpr vec7 a2 = {6.590560e+00,-2.923150e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + constexpr vec7 a3 = {1.454440e+01,-1.022950e+01,0.000000e+00,0.000000e+00,4.590370e-02,0.000000e+00,-1.500000e+00}; return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3); } double n15pg_frac(const vec7 &T9){ - double f1=n15pg_rate(T9); - double f2=n15pa_rate(T9); + const double f1=n15pg_rate(T9); + const double f2=n15pa_rate(T9); return f1/(f1+f2); } // o16(p,g)f17 then f17 -> o17(p,a)n14 double o16p_rate(const vec7 &T9){ - vec7 a={1.909040e+01,0.000000e+00,-1.669600e+01,-1.162520e+00,2.677030e-01,-3.384110e-02,-6.666670e-01}; + constexpr vec7 a={1.909040e+01,0.000000e+00,-1.669600e+01,-1.162520e+00,2.677030e-01,-3.384110e-02,-6.666670e-01}; return rate_fit(T9,a); } // o16(a,g)ne20 double o16a_rate(const vec7 &T9){ - vec7 a1={2.390300e+01,0.000000e+00,-3.972620e+01,-2.107990e-01,4.428790e-01,-7.977530e-02,-6.666670e-01}; - vec7 a2={3.885710e+00,-1.035850e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec7 a3={9.508480e+00,-1.276430e+01,0.000000e+00,-3.659250e+00,7.142240e-01,-1.075080e-03,-1.500000e+00}; + constexpr vec7 a1={2.390300e+01,0.000000e+00,-3.972620e+01,-2.107990e-01,4.428790e-01,-7.977530e-02,-6.666670e-01}; + constexpr vec7 a2={3.885710e+00,-1.035850e+01,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + constexpr vec7 a3={9.508480e+00,-1.276430e+01,0.000000e+00,-3.659250e+00,7.142240e-01,-1.075080e-03,-1.500000e+00}; return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3); } // ne20(a,g)mg24 double ne20a_rate(const vec7 &T9){ - vec7 a1={2.450580e+01,0.000000e+00,-4.625250e+01,5.589010e+00,7.618430e+00,-3.683000e+00,-6.666670e-01}; - vec7 a2={-3.870550e+01,-2.506050e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec7 a3={1.983070e+00,-9.220260e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; - vec7 a4={-8.798270e+00,-1.278090e+01,0.000000e+00,1.692290e+01,-2.573250e+00,2.089970e-01,-1.500000e+00}; + constexpr vec7 a1={2.450580e+01,0.000000e+00,-4.625250e+01,5.589010e+00,7.618430e+00,-3.683000e+00,-6.666670e-01}; + constexpr vec7 a2={-3.870550e+01,-2.506050e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + constexpr vec7 a3={1.983070e+00,-9.220260e+00,0.000000e+00,0.000000e+00,0.000000e+00,0.000000e+00,-1.500000e+00}; + constexpr vec7 a4={-8.798270e+00,-1.278090e+01,0.000000e+00,1.692290e+01,-2.573250e+00,2.089970e-01,-1.500000e+00}; return rate_fit(T9,a1) + rate_fit(T9,a2) + rate_fit(T9,a3) + rate_fit(T9,a4); } // c12(c12,a)ne20 double c12c12_rate(const vec7 &T9){ - vec7 a={6.128630e+01,0.000000e+00,-8.416500e+01,-1.566270e+00,-7.360840e-02,-7.279700e-02,-6.666670e-01}; + constexpr vec7 a={6.128630e+01,0.000000e+00,-8.416500e+01,-1.566270e+00,-7.360840e-02,-7.279700e-02,-6.666670e-01}; return rate_fit(T9,a); } // c12(o16,a)mg24 double c12o16_rate(const vec7 &T9){ - vec7 a={4.853410e+01,3.720400e-01,-1.334130e+02,5.015720e+01,-3.159870e+00,1.782510e-02,-2.370270e+01}; + constexpr vec7 a={4.853410e+01,3.720400e-01,-1.334130e+02,5.015720e+01,-3.159870e+00,1.782510e-02,-2.370270e+01}; return rate_fit(T9,a); } @@ -244,209 +244,211 @@ namespace serif::network::approx8{ // a Jacobian matrix for implicit solvers - void Jacobian::operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) { + void Jacobian::operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const { serif::constant::Constants& constants = serif::constant::Constants::getInstance(); const double avo = constants.get("N_a").value; const double clight = constants.get("c").value; // EOS - vec7 T9=get_T9_array(y[Net::itemp]); + const vec7 T9=get_T9_array(y[Approx8Net::iTemp]); // evaluate rates once per call - double rpp=pp_rate(T9); - double r33=he3he3_rate(T9); - double r34=he3he4_rate(T9); - double r3a=triple_alpha_rate(T9); - double rc12p=c12p_rate(T9); - double rc12a=c12a_rate(T9); - double rn14p=n14p_rate(T9); - double rn14a=n14a_rate(T9); - double ro16p=o16p_rate(T9); - double ro16a=o16a_rate(T9); - double rne20a=ne20a_rate(T9); - double r1212=c12c12_rate(T9); - double r1216=c12o16_rate(T9); + const double rpp=pp_rate(T9); + const double r33=he3he3_rate(T9); + const double r34=he3he4_rate(T9); + const double r3a=triple_alpha_rate(T9); + const double rc12p=c12p_rate(T9); + const double rc12a=c12a_rate(T9); + const double rn14p=n14p_rate(T9); + const double rn14a=n14a_rate(T9); + const double ro16p=o16p_rate(T9); + const double ro16a=o16a_rate(T9); + const double rne20a=ne20a_rate(T9); + const double r1212=c12c12_rate(T9); + const double r1216=c12o16_rate(T9); - double pfrac=n15pg_frac(T9); - double afrac=1-pfrac; + const double pFrac=n15pg_frac(T9); + const double aFrac=1-pFrac; - double yh1 = y[Net::ih1]; - double yhe3 = y[Net::ihe3]; - double yhe4 = y[Net::ihe4]; - double yc12 = y[Net::ic12]; - double yn14 = y[Net::in14]; - double yo16 = y[Net::io16]; - double yne20 = y[Net::ine20]; + const double yh1 = y[Approx8Net::ih1]; + const double yhe3 = y[Approx8Net::ihe3]; + const double yhe4 = y[Approx8Net::ihe4]; + const double yc12 = y[Approx8Net::ic12]; + const double yn14 = y[Approx8Net::in14]; + const double yo16 = y[Approx8Net::io16]; + const double yne20 = y[Approx8Net::ine20]; // zero all elements to begin - for (int i=0; i < Net::nvar; i++) { - for (int j=0; j < Net::nvar; j++) { - J(i,j)=0.0; + for (int i=0; i < Approx8Net::nVar; i++) { + for (int j=0; j < Approx8Net::nVar; j++) { + J(i,j)=0.0; } } - + // h1 jacobian elements - J(Net::ih1,Net::ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p; - J(Net::ih1,Net::ihe3) = 2*yhe3*r33 - yhe4*r34; - J(Net::ih1,Net::ihe4) = -yhe3*r34; - J(Net::ih1,Net::ic12) = -2*yh1*rc12p; - J(Net::ih1,Net::in14) = -2*yh1*rn14p; - J(Net::ih1,Net::io16) = -2*yh1*ro16p; + J(Approx8Net::ih1,Approx8Net::ih1) = -3*yh1*rpp - 2*yc12*rc12p -2*yn14*rn14p -2*yo16*ro16p; + J(Approx8Net::ih1,Approx8Net::ihe3) = 2*yhe3*r33 - yhe4*r34; + J(Approx8Net::ih1,Approx8Net::ihe4) = -yhe3*r34; + J(Approx8Net::ih1,Approx8Net::ic12) = -2*yh1*rc12p; + J(Approx8Net::ih1,Approx8Net::in14) = -2*yh1*rn14p; + J(Approx8Net::ih1,Approx8Net::io16) = -2*yh1*ro16p; // he3 jacobian elements - J(Net::ihe3,Net::ih1) = yh1*rpp; - J(Net::ihe3,Net::ihe3) = -2*yhe3*r33 - yhe4*r34; - J(Net::ihe3,Net::ihe4) = -yhe3*r34; - + J(Approx8Net::ihe3,Approx8Net::ih1) = yh1*rpp; + J(Approx8Net::ihe3,Approx8Net::ihe3) = -2*yhe3*r33 - yhe4*r34; + J(Approx8Net::ihe3,Approx8Net::ihe4) = -yhe3*r34; + // he4 jacobian elements - J(Net::ihe4,Net::ih1) = yn14*afrac*rn14p + yo16*ro16p; - J(Net::ihe4,Net::ihe3) = yhe3*r33 - yhe4*r34; - J(Net::ihe4,Net::ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a; - J(Net::ihe4,Net::ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216; - J(Net::ihe4,Net::in14) = yh1*afrac*rn14p - 1.5*yhe4*rn14a; - J(Net::ihe4,Net::io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216; - J(Net::ihe4,Net::ine20) = -yhe4*rne20a; + J(Approx8Net::ihe4,Approx8Net::ih1) = yn14*aFrac*rn14p + yo16*ro16p; + J(Approx8Net::ihe4,Approx8Net::ihe3) = yhe3*r33 - yhe4*r34; + J(Approx8Net::ihe4,Approx8Net::ihe4) = yhe3*r34 - 1.5*yhe4*yhe4*r3a - yc12*rc12a - 1.5*yn14*rn14a - yo16*ro16a - yne20*rne20a; + J(Approx8Net::ihe4,Approx8Net::ic12) = -yhe4*rc12a + yc12*r1212 + yo16*r1216; + J(Approx8Net::ihe4,Approx8Net::in14) = yh1*aFrac*rn14p - 1.5*yhe4*rn14a; + J(Approx8Net::ihe4,Approx8Net::io16) = yh1*ro16p - yhe4*ro16a + yc12*r1216; + J(Approx8Net::ihe4,Approx8Net::ine20) = -yhe4*rne20a; // c12 jacobian elements - J(Net::ic12,Net::ih1) = -yc12*rc12p + yn14*afrac*rn14p; - J(Net::ic12,Net::ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a; - J(Net::ic12,Net::ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212; - J(Net::ic12,Net::in14) = yh1*yn14*afrac*rn14p; - J(Net::ic12,Net::io16) = -yc12*r1216; + J(Approx8Net::ic12,Approx8Net::ih1) = -yc12*rc12p + yn14*aFrac*rn14p; + J(Approx8Net::ic12,Approx8Net::ihe4) = 0.5*yhe4*yhe4*r3a - yhe4*rc12a; + J(Approx8Net::ic12,Approx8Net::ic12) = -yh1*rc12p - yhe4*rc12a - yo16*r1216 - 2*yc12*r1212; + J(Approx8Net::ic12,Approx8Net::in14) = yh1*yn14*aFrac*rn14p; + J(Approx8Net::ic12,Approx8Net::io16) = -yc12*r1216; // n14 jacobian elements - J(Net::in14,Net::ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p; - J(Net::in14,Net::ihe4) = -yn14*rn14a; - J(Net::in14,Net::ic12) = yh1*rc12p; - J(Net::in14,Net::in14) = -yh1*rn14p - yhe4*rn14a; - J(Net::in14,Net::io16) = yo16*ro16p; - + J(Approx8Net::in14,Approx8Net::ih1) = yc12*rc12p - yn14*rn14p + yo16*ro16p; + J(Approx8Net::in14,Approx8Net::ihe4) = -yn14*rn14a; + J(Approx8Net::in14,Approx8Net::ic12) = yh1*rc12p; + J(Approx8Net::in14,Approx8Net::in14) = -yh1*rn14p - yhe4*rn14a; + J(Approx8Net::in14,Approx8Net::io16) = yo16*ro16p; + // o16 jacobian elements - J(Net::io16,Net::ih1) = yn14*pfrac*rn14p - yo16*ro16p; - J(Net::io16,Net::ihe4) = yc12*rc12a - yo16*ro16a; - J(Net::io16,Net::ic12) = yhe4*rc12a - yo16*r1216; - J(Net::io16,Net::in14) = yh1*pfrac*rn14p; - J(Net::io16,Net::io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a; + J(Approx8Net::io16,Approx8Net::ih1) = yn14*pFrac*rn14p - yo16*ro16p; + J(Approx8Net::io16,Approx8Net::ihe4) = yc12*rc12a - yo16*ro16a; + J(Approx8Net::io16,Approx8Net::ic12) = yhe4*rc12a - yo16*r1216; + J(Approx8Net::io16,Approx8Net::in14) = yh1*pFrac*rn14p; + J(Approx8Net::io16,Approx8Net::io16) = yh1*ro16p - yc12*r1216 -yhe4*ro16a; // ne20 jacobian elements - J(Net::ine20,Net::ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a; - J(Net::ine20,Net::ic12) = yc12*r1212; - J(Net::ine20,Net::in14) = yhe4*rn14a; - J(Net::ine20,Net::io16) = yo16*ro16a; - J(Net::ine20,Net::ine20) = -yhe4*rne20a; + J(Approx8Net::ine20,Approx8Net::ihe4) = yn14*rn14a + yo16*ro16a - yne20*rne20a; + J(Approx8Net::ine20,Approx8Net::ic12) = yc12*r1212; + J(Approx8Net::ine20,Approx8Net::in14) = yhe4*rn14a; + J(Approx8Net::ine20,Approx8Net::io16) = yo16*ro16a; + J(Approx8Net::ine20,Approx8Net::ine20) = -yhe4*rne20a; // mg24 jacobian elements - J(Net::img24,Net::ihe4) = yne20*rne20a; - J(Net::img24,Net::ic12) = yo16*r1216; - J(Net::img24,Net::io16) = yc12*r1216; - J(Net::img24,Net::ine20) = yhe4*rne20a; + J(Approx8Net::img24,Approx8Net::ihe4) = yne20*rne20a; + J(Approx8Net::img24,Approx8Net::ic12) = yo16*r1216; + J(Approx8Net::img24,Approx8Net::io16) = yc12*r1216; + J(Approx8Net::img24,Approx8Net::ine20) = yhe4*rne20a; - // energy accountings - for (int j=0; j("Network:Approx8:Stiff:AbsTol", 1.0e-6); @@ -463,7 +465,7 @@ namespace serif::network::approx8{ std::make_pair(ODE(), Jacobian()), m_y, 0.0, - m_tmax, + m_tMax, m_dt0 ); @@ -474,31 +476,33 @@ namespace serif::network::approx8{ ODE(), m_y, 0.0, - m_tmax, + m_tMax, m_dt0 ); } - double ysum = 0.0; - for (int i = 0; i < Net::niso; i++) { - m_y[i] *= Net::aion[i]; - ysum += m_y[i]; + double ySum = 0.0; + for (int i = 0; i < Approx8Net::nIso; i++) { + m_y[i] *= Approx8Net::aIon[i]; + ySum += m_y[i]; } - for (int i = 0; i < Net::niso; i++) { - m_y[i] /= ysum; + for (int i = 0; i < Approx8Net::nIso; i++) { + m_y[i] /= ySum; } NetOut netOut; std::vector outComposition; - outComposition.reserve(Net::nvar); + outComposition.reserve(Approx8Net::nVar); - for (int i = 0; i < Net::niso; i++) { + for (int i = 0; i < Approx8Net::nIso; i++) { outComposition.push_back(m_y[i]); } - netOut.energy = m_y[Net::iener]; - netOut.composition = outComposition; + netOut.energy = m_y[Approx8Net::iEnergy]; netOut.num_steps = num_steps; + const std::vector symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"}; + netOut.composition = serif::composition::Composition(symbols, outComposition); + return netOut; } @@ -507,31 +511,26 @@ namespace serif::network::approx8{ } vector_type Approx8Network::convert_netIn(const NetIn &netIn) { - if (netIn.composition.size() != Net::niso) { - LOG_ERROR(m_logger, "Error: composition size mismatch in convert_netIn"); - throw std::runtime_error("Error: composition size mismatch in convert_netIn"); - } + vector_type y(Approx8Net::nVar, 0.0); + y[Approx8Net::ih1] = netIn.composition.getMassFraction("H-1"); + y[Approx8Net::ihe3] = netIn.composition.getMassFraction("He-3"); + y[Approx8Net::ihe4] = netIn.composition.getMassFraction("He-4"); + y[Approx8Net::ic12] = netIn.composition.getMassFraction("C-12"); + y[Approx8Net::in14] = netIn.composition.getMassFraction("N-14"); + y[Approx8Net::io16] = netIn.composition.getMassFraction("O-16"); + y[Approx8Net::ine20] = netIn.composition.getMassFraction("Ne-20"); + y[Approx8Net::img24] = netIn.composition.getMassFraction("Mg-24"); + y[Approx8Net::iTemp] = netIn.temperature; + y[Approx8Net::iDensity] = netIn.density; + y[Approx8Net::iEnergy] = netIn.energy; - vector_type y(Net::nvar, 0.0); - y[Net::ih1] = netIn.composition[0]; - y[Net::ihe3] = netIn.composition[1]; - y[Net::ihe4] = netIn.composition[2]; - y[Net::ic12] = netIn.composition[3]; - y[Net::in14] = netIn.composition[4]; - y[Net::io16] = netIn.composition[5]; - y[Net::ine20] = netIn.composition[6]; - y[Net::img24] = netIn.composition[7]; - y[Net::itemp] = netIn.temperature; - y[Net::iden] = netIn.density; - y[Net::iener] = netIn.energy; - - double ysum = 0.0; - for (int i = 0; i < Net::niso; i++) { - y[i] /= Net::aion[i]; - ysum += y[i]; + double ySum = 0.0; + for (int i = 0; i < Approx8Net::nIso; i++) { + y[i] /= Approx8Net::aIon[i]; + ySum += y[i]; } - for (int i = 0; i < Net::niso; i++) { - y[i] /= ysum; + for (int i = 0; i < Approx8Net::nIso; i++) { + y[i] /= ySum; } return y; diff --git a/src/network/private/network.cpp b/src/network/private/network.cpp index d3b3b25..b0e3e2d 100644 --- a/src/network/private/network.cpp +++ b/src/network/private/network.cpp @@ -19,18 +19,50 @@ // // *********************************************************************** */ #include "network.h" + +#include "approx8.h" #include "probe.h" #include "quill/LogMacros.h" namespace serif::network { - Network::Network() : + Network::Network(const NetworkFormat format) : m_config(serif::config::Config::getInstance()), m_logManager(serif::probe::LogManager::getInstance()), - m_logger(m_logManager.getLogger("log")) { + m_logger(m_logManager.getLogger("log")), + m_format(format) { + if (format == NetworkFormat::UNKNOWN) { + LOG_ERROR(m_logger, "nuclearNetwork::Network::Network() called with UNKNOWN format"); + throw std::runtime_error("nuclearNetwork::Network::Network() called with UNKNOWN format"); + } } + + NetworkFormat Network::getFormat() const { + return m_format; + } + + NetworkFormat Network::setFormat(const NetworkFormat format) { + const NetworkFormat oldFormat = m_format; + m_format = format; + return oldFormat; + } + NetOut Network::evaluate(const NetIn &netIn) { - // You can throw an exception here or log a warning if it should never be used. - LOG_ERROR(m_logger, "nuclearNetwork::Network::evaluate() is not implemented"); - throw std::runtime_error("nuclearNetwork::Network::evaluate() is not implemented"); + NetOut netOut; + switch (m_format) { + case APPROX8: { + approx8::Approx8Network network; + netOut = network.evaluate(netIn); + break; + } + case UNKNOWN: { + LOG_ERROR(m_logger, "Network format {} is not implemented.", FormatStringLookup.at(m_format)); + throw std::runtime_error("Network format not implemented."); + } + default: { + LOG_ERROR(m_logger, "Unknown network format."); + throw std::runtime_error("Unknown network format."); + } + } + return netOut; } } diff --git a/src/network/public/approx8.h b/src/network/public/approx8.h index 711d2a1..634225e 100644 --- a/src/network/public/approx8.h +++ b/src/network/public/approx8.h @@ -31,54 +31,55 @@ * @brief Header file for the Approx8 nuclear reaction network. * * This file contains the definitions and declarations for the Approx8 nuclear reaction network. - * The network is based on Frank Timmes' "aprox8" and includes 8 isotopes and various nuclear reactions. + * The network is based on Frank Timmes' "approx8" and includes 8 isotopes and various nuclear reactions. * The rates are evaluated using a fitting function with coefficients from reaclib.jinaweb.org. */ -/** - * @typedef vector_type - * @brief Alias for a vector of doubles using Boost uBLAS. - */ -typedef boost::numeric::ublas::vector< double > vector_type; - -/** - * @typedef matrix_type - * @brief Alias for a matrix of doubles using Boost uBLAS. - */ -typedef boost::numeric::ublas::matrix< double > matrix_type; - -/** - * @typedef vec7 - * @brief Alias for a std::array of 7 doubles. - */ -typedef std::array vec7; namespace serif::network::approx8{ + /** + * @typedef vector_type + * @brief Alias for a vector of doubles using Boost uBLAS. + */ + typedef boost::numeric::ublas::vector< double > vector_type; + + /** + * @typedef matrix_type + * @brief Alias for a matrix of doubles using Boost uBLAS. + */ + typedef boost::numeric::ublas::matrix< double > matrix_type; + + /** + * @typedef vec7 + * @brief Alias for a std::array of 7 doubles. + */ + typedef std::array vec7; + using namespace boost::numeric::odeint; /** - * @struct Net + * @struct Approx8Net * @brief Contains constants and arrays related to the nuclear network. */ - struct Net{ - const static int ih1=0; - const static int ihe3=1; - const static int ihe4=2; - const static int ic12=3; - const static int in14=4; - const static int io16=5; - const static int ine20=6; - const static int img24=7; + struct Approx8Net{ + static constexpr int ih1=0; + static constexpr int ihe3=1; + static constexpr int ihe4=2; + static constexpr int ic12=3; + static constexpr int in14=4; + static constexpr int io16=5; + static constexpr int ine20=6; + static constexpr int img24=7; - const static int itemp=img24+1; - const static int iden =itemp+1; - const static int iener=iden+1; + static constexpr int iTemp=img24+1; + static constexpr int iDensity =iTemp+1; + static constexpr int iEnergy=iDensity+1; - const static int niso=img24+1; // number of isotopes - const static int nvar=iener+1; // number of variables + static constexpr int nIso=img24+1; // number of isotopes + static constexpr int nVar=iEnergy+1; // number of variables - static constexpr std::array aion = { + static constexpr std::array aIon = { 1, 3, 4, @@ -89,7 +90,7 @@ namespace serif::network::approx8{ 24 }; - static constexpr std::array mion = { + static constexpr std::array mIon = { 1.67262164e-24, 5.00641157e-24, 6.64465545e-24, @@ -270,10 +271,9 @@ namespace serif::network::approx8{ * @brief Calculates the Jacobian matrix. * @param y State vector. * @param J Jacobian matrix. - * @param t Time. * @param dfdt Derivative of the state vector. */ - void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ); + void operator() ( const vector_type &y, matrix_type &J, double /* t */, vector_type &dfdt ) const; }; /** @@ -285,23 +285,24 @@ namespace serif::network::approx8{ * @brief Calculates the derivatives of the state vector. * @param y State vector. * @param dydt Derivative of the state vector. - * @param t Time. */ - void operator() ( const vector_type &y, vector_type &dydt, double /* t */); + void operator() ( const vector_type &y, vector_type &dydt, double /* t */) const; }; /** * @class Approx8Network * @brief Class for the Approx8 nuclear reaction network. */ - class Approx8Network : public Network { + class Approx8Network final : public Network { public: + Approx8Network(); + /** * @brief Evaluates the nuclear network. * @param netIn Input parameters for the network. * @return Output results from the network. */ - virtual NetOut evaluate(const NetIn &netIn); + NetOut evaluate(const NetIn &netIn) override; /** * @brief Sets whether the solver should use a stiff method. @@ -313,11 +314,11 @@ namespace serif::network::approx8{ * @brief Checks if the solver is using a stiff method. * @return Boolean indicating if a stiff method is being used. */ - bool isStiff() { return m_stiff; } + bool isStiff() const { return m_stiff; } private: vector_type m_y; - double m_tmax; - double m_dt0; + double m_tMax = 0; + double m_dt0 = 0; bool m_stiff = false; /** @@ -325,7 +326,8 @@ namespace serif::network::approx8{ * @param netIn Input parameters for the network. * @return Internal state vector. */ - vector_type convert_netIn(const NetIn &netIn); + static vector_type convert_netIn(const NetIn &netIn); }; + } // namespace nnApprox8 diff --git a/src/network/public/network.h b/src/network/public/network.h index 0ae4677..0d7692b 100644 --- a/src/network/public/network.h +++ b/src/network/public/network.h @@ -25,15 +25,27 @@ #include "probe.h" #include "config.h" #include "quill/Logger.h" +#include "composition.h" +#include namespace serif::network { + enum NetworkFormat { + APPROX8, ///< Approx8 nuclear reaction network format. + UNKNOWN, + }; + + static inline std::unordered_map FormatStringLookup = { + {APPROX8, "Approx8"}, + {UNKNOWN, "Unknown"} + }; + /** * @struct NetIn * @brief Input structure for the network evaluation. - * + * * This structure holds the input parameters required for the network evaluation. - * + * * Example usage: * @code * nuclearNetwork::NetIn netIn; @@ -46,8 +58,8 @@ namespace serif::network { * @endcode */ struct NetIn { - std::vector composition; ///< Composition of the network - double tmax; ///< Maximum time + serif::composition::Composition composition; ///< Composition of the network + double tMax; ///< Maximum time double dt0; ///< Initial time step double temperature; ///< Temperature in Kelvin double density; ///< Density in g/cm^3 @@ -69,7 +81,7 @@ namespace serif::network { * @endcode */ struct NetOut { - std::vector composition; ///< Composition of the network after evaluation + serif::composition::Composition composition; ///< Composition of the network after evaluation int num_steps; ///< Number of steps taken in the evaluation double energy; ///< Energy in ergs after evaluation }; @@ -90,9 +102,12 @@ namespace serif::network { */ class Network { public: - Network(); + explicit Network(const NetworkFormat format = NetworkFormat::APPROX8); virtual ~Network() = default; + NetworkFormat getFormat() const; + NetworkFormat setFormat(const NetworkFormat format); + /** * @brief Evaluate the network based on the input parameters. * @@ -105,6 +120,10 @@ namespace serif::network { serif::config::Config& m_config; ///< Configuration instance serif::probe::LogManager& m_logManager; ///< Log manager instance quill::Logger* m_logger; ///< Logger instance + + NetworkFormat m_format; ///< Format of the network }; + + } // namespace nuclearNetwork diff --git a/tests/network/approx8Test.cpp b/tests/network/approx8Test.cpp index 2bfe7fa..899eaf6 100644 --- a/tests/network/approx8Test.cpp +++ b/tests/network/approx8Test.cpp @@ -4,6 +4,7 @@ #include "approx8.h" #include "config.h" #include "network.h" +#include "composition.h" #include @@ -32,19 +33,31 @@ TEST_F(approx8Test, evaluate) { serif::network::NetIn netIn; std::vector comp = {0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4}; + std::vector symbols = {"H-1", "He-3", "He-4", "C-12", "N-14", "O-16", "Ne-20", "Mg-24"}; - netIn.composition = comp; + serif::composition::Composition composition; + composition.registerSymbol(symbols, true); + composition.setMassFraction(symbols, comp); + bool isFinalized = composition.finalize(true); + EXPECT_TRUE(isFinalized); + + netIn.composition = composition; netIn.temperature = 1e7; netIn.density = 1e2; netIn.energy = 0.0; - netIn.tmax = 3.15e17; + netIn.tMax = 3.15e17; netIn.dt0 = 1e12; serif::network::NetOut netOut; EXPECT_NO_THROW(netOut = network.evaluate(netIn)); - EXPECT_DOUBLE_EQ(netOut.composition[serif::network::approx8::Net::ih1], 0.50166260916650918); - EXPECT_DOUBLE_EQ(netOut.composition[serif::network::approx8::Net::ihe4],0.48172270591286032); - EXPECT_DOUBLE_EQ(netOut.energy, 1.6433049870528356e+18); + double energyFraction = netOut.energy / 1.6433051127589775E+18; + double H1MassFraction = netOut.composition.getMassFraction("H-1")/ 0.50166262445895604; + double He4MassFraction = netOut.composition.getMassFraction("He-4") / 0.48172273720971226; + + double relError = 1e-8; + EXPECT_NEAR(H1MassFraction, 1.0, relError); + EXPECT_NEAR(He4MassFraction, 1.0, relError); + EXPECT_NEAR(energyFraction, 1.0, relError); } diff --git a/tests/network/meson.build b/tests/network/meson.build index 7b84063..1e6ccf3 100644 --- a/tests/network/meson.build +++ b/tests/network/meson.build @@ -11,7 +11,7 @@ foreach test_file : test_sources test_exe = executable( exe_name, test_file, - dependencies: [gtest_dep, gtest_main, network_dep], + dependencies: [gtest_dep, gtest_main, network_dep, species_weight_dep, composition_dep], include_directories: include_directories('../../src/network/public'), link_with: libnetwork, # Link the dobj library install_rpath: '@loader_path/../../src' # Ensure runtime library path resolves correctly
 ▼Nserif
 ▼Ncomposition
 CComposition
 CCompositionManages the composition of elements
 CCompositionEntryRepresents an entry in the composition with a symbol and mass fraction
 CGlobalCompositionRepresents the global composition of a system. This tends to be used after finalize and is primarily for internal use
 ▼Nconstant