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/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 --- 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 6af4a88..aad84dc 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,636 @@ 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) { + return chemSpecies::species.contains(symbol); + } + + double Composition::setMassFraction(const std::string& symbol, const double& mass_fraction) { + if (!m_registeredSymbols.contains(symbol)) { + LOG_ERROR(m_logger, "Symbol {} is not registered.", symbol); + throw 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; -} + const 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(const 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) { - if (m_compositions.count(symbol) == 0) { + 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 otherMassFrac = 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); + 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.contains(symbol)) { 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: m_compositions | std::views::keys) { + 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.contains(symbol)) { + 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: m_compositions | std::views::keys) { + 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.contains(symbol)) { + 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 zSum = 0.0; + + // Loop through all registered species in the composition. + for (const auto &val: m_compositions | std::views::values) { + zSum += (val.mass_fraction() * val.m_isotope.z())/val.m_isotope.a(); + } + + const double mean_A = m_meanParticleMass * zSum; + return mean_A; + } + + Composition Composition::subset(const std::vector& symbols, std::string method) const { + const std::array methods = {"norm", "none"}; + + if (std::ranges::find(methods, method) == methods.end()) { + const std::string errorMessage = "Invalid method: " + method + ". Valid methods are 'norm' and 'none'."; + LOG_ERROR(m_logger, "Invalid method: {}. Valid methods are norm and none.", method); + throw std::runtime_error(errorMessage); + } + + Composition subsetComposition; + for (const auto& symbol : symbols) { + if (!m_compositions.contains(symbol)) { + 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") { + const 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(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 | std::views::values) { + 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; + } + + 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; + } + + /// 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 7226393..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 @@ -34,6 +33,20 @@ #include "atomicSpecies.h" namespace serif::composition { + struct CanonicalComposition { + double X = 0.0; ///< Mass fraction of Hydrogen. + double Y = 0.0; ///< Mass fraction of Helium. + double Z = 0.0; ///< Mass fraction of Metals. + + friend std::ostream& operator<<(std::ostream& os, const CanonicalComposition& composition) { + os << ""; + return os; + } + }; + /** * @brief Represents the global composition of a system. This tends to be used after finalize and is primarily for internal use. */ @@ -68,7 +81,7 @@ namespace serif::composition { * @brief Constructs a CompositionEntry with the given symbol and mode. * @param symbol The chemical symbol of the species. * @param massFracMode True if mass fraction mode, false if number fraction mode. - * @example + * *Example Usage:* * @code * CompositionEntry entry("H", true); * @endcode @@ -185,13 +198,13 @@ namespace serif::composition { * - The only exception to the finalize rule is if the compositon was constructed with symbols and mass fractions at instantiation time. In this case, the composition is automatically finalized. * however, this means that the composition passed to the constructor must be valid. * - * @example Constructing a finalized composition with symbols and mass fractions: + * *Example Usage:* Constructing a finalized composition with symbols and mass fractions: * @code * std::vector 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); @@ -220,7 +233,7 @@ namespace serif::composition { * @param symbol The symbol to check. * @return True if the symbol is valid, false otherwise. */ - bool isValidSymbol(const std::string& symbol) const; + static bool isValidSymbol(const std::string& symbol); /** * @brief Checks if the given mass fractions are valid. @@ -271,31 +284,31 @@ namespace serif::composition { /** * @brief Constructs a Composition with the given symbols. * @param symbols The symbols to initialize the composition with. - * @example + * *Example Usage:* * @code * std::vector symbols = {"H", "O"}; * 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. * @param symbols The symbols to initialize the composition with. - * @example + * *Example Usage:* * @code * std::set symbols = {"H", "O"}; * 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. * @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}; @@ -304,11 +317,19 @@ 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. * @param massFracMode True if mass fraction mode, false if number fraction mode. - * @example + * *Example Usage:* * @code * Composition comp; * comp.registerSymbol("H"); @@ -320,7 +341,7 @@ namespace serif::composition { * @brief Registers multiple new symbols. * @param symbols The symbols to register. * @param massFracMode True if mass fraction mode, false if number fraction mode. - * @example + * *Example Usage:* * @code * std::vector symbols = {"H", "O"}; * Composition comp; @@ -333,14 +354,14 @@ namespace serif::composition { * @brief Gets the registered symbols. * @return A set of registered symbols. */ - std::set getRegisteredSymbols() const; + [[nodiscard]] std::set getRegisteredSymbols() const; /** * @brief Sets the mass fraction for a given symbol. * @param symbol The symbol to set the mass fraction for. * @param mass_fraction The mass fraction to set. * @return The mass fraction that was set. - * @example + * *Example Usage:* * @code * Composition comp; * comp.setMassFraction("H", 0.1); @@ -353,7 +374,7 @@ namespace serif::composition { * @param symbols The symbols to set the mass fraction for. * @param mass_fractions The mass fractions corresponding to the symbols. * @return A vector of mass fractions that were set. - * @example + * *Example Usage:* * @code * std::vector symbols = {"H", "O"}; * std::vector mass_fractions = {0.1, 0.9}; @@ -390,40 +411,52 @@ namespace serif::composition { * @brief Gets the mass fractions of all compositions. * @return An unordered map of compositions with their mass fractions. */ - std::unordered_map 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. @@ -446,6 +479,16 @@ namespace serif::composition { */ void setCompositionMode(bool massFracMode); + /** + * @brief Gets the current canonical composition (X, Y, Z). + * @param harsh If true, this will throw an error if X-Y != Z where Z is computed as the sum of all other elements. + * @return True if mass fraction mode, false if number fraction mode. + * + * @throws std::runtime_error if the composition is not finalized or if the canonical composition cannot be computed. + * @throws std::runtime_error if harsh is true and the canonical composition is not valid. + */ + [[nodiscard]] CanonicalComposition getCanonicalComposition(bool harsh=false) const; + /** * @brief Overloaded output stream operator for Composition. * @param os The output stream. @@ -463,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/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..207d746 --- /dev/null +++ b/src/eos/public/EOS.h @@ -0,0 +1,266 @@ +#pragma once + +#include "EOSio.h" +#include "helm.h" +#include +#include "composition.h" +#include + +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). + friend std::ostream& operator<<(std::ostream& os, const EOSInput& input) { + os << ""; + return os; + } + }; + + /** + * @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 { + 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). + + 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"). + + friend std::ostream& operator<<(std::ostream& os, const EOSParameter& param) { + os << std::setprecision(5) << ""; + return os; + } + + + }; + + /** + * @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"}; ///< 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). + * @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; + + 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; + } + }; + + /** + * @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/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/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 ""; diff --git a/tests/eos/eosTest.cpp b/tests/eos/eosTest.cpp index b7a670b..d98a25c 100644 --- a/tests/eos/eosTest.cpp +++ b/tests/eos/eosTest.cpp @@ -1,5 +1,4 @@ #include -#include #include #include #include @@ -7,9 +6,11 @@ #include "helm.h" #include "resourceManager.h" #include "config.h" +#include "composition.h" +#include "EOS.h" /** - * @file constTest.cpp + * @file eosTest.cpp * @brief Unit tests for the const class. */ @@ -26,44 +27,43 @@ 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) { + constexpr int nel=3; + double xMass[nel], aIon[nel], zIon[nel]; + serif::eos::helmholtz::HELMEOSInput eos1; - const int nel=3; - double xmass[nel], aion[nel], zion[nel]; - serif::eos::helmholtz::EOSInput 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")); 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; + constexpr double absErr = 1e-12; //Check composition info EXPECT_NEAR( helmEos.ye, 8.75e-01, absErr); @@ -85,3 +85,44 @@ 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", + "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.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(); + 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)); + eosOutput = EOS.get(eosInput); + constexpr double absErr = 1e-8; + EXPECT_NEAR(eosOutput.pressure.total, 6.9548533046915791E+22, absErr); +} 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/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 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
 ▼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