Merge pull request #68 from tboudreaux/refactor/useComposition
Refactor Network and EOS modules to use composition module
This commit is contained in:
GitHub
@@ -1,3 +1,4 @@
|
||||
species_weight_dep = declare_dependency(
|
||||
include_directories: include_directories('include'),
|
||||
)
|
||||
)
|
||||
message('✅ SERiF species_weight dependency declared')
|
||||
@@ -106,7 +106,7 @@ $(function(){initNavTree('annotated.html',''); initResizable(true); });
|
||||
<div class="levels">[detail level <span onclick="javascript:dynsection.toggleLevel(1);">1</span><span onclick="javascript:dynsection.toggleLevel(2);">2</span><span onclick="javascript:dynsection.toggleLevel(3);">3</span><span onclick="javascript:dynsection.toggleLevel(4);">4</span>]</div><table class="directory">
|
||||
<tr id="row_0_" class="even"><td class="entry"><span style="width:0px;display:inline-block;"> </span><span id="arr_0_" class="arrow" onclick="dynsection.toggleFolder('0_')">▼</span><span class="icona"><span class="icon">N</span></span><a class="el" href="namespaceserif.html" target="_self">serif</a></td><td class="desc"></td></tr>
|
||||
<tr id="row_0_0_" class="odd"><td class="entry"><span style="width:16px;display:inline-block;"> </span><span id="arr_0_0_" class="arrow" onclick="dynsection.toggleFolder('0_0_')">▼</span><span class="icona"><span class="icon">N</span></span><a class="el" href="namespaceserif_1_1composition.html" target="_self">composition</a></td><td class="desc"></td></tr>
|
||||
<tr id="row_0_0_0_" class="even"><td class="entry"><span style="width:48px;display:inline-block;"> </span><span class="icona"><span class="icon">C</span></span><a class="el" href="classserif_1_1composition_1_1_composition.html" target="_self">Composition</a></td><td class="desc"></td></tr>
|
||||
<tr id="row_0_0_0_" class="even"><td class="entry"><span style="width:48px;display:inline-block;"> </span><span class="icona"><span class="icon">C</span></span><a class="el" href="classserif_1_1composition_1_1_composition.html" target="_self">Composition</a></td><td class="desc">Manages the composition of elements </td></tr>
|
||||
<tr id="row_0_0_1_" class="odd"><td class="entry"><span style="width:48px;display:inline-block;"> </span><span class="icona"><span class="icon">C</span></span><a class="el" href="structserif_1_1composition_1_1_composition_entry.html" target="_self">CompositionEntry</a></td><td class="desc">Represents an entry in the composition with a symbol and mass fraction </td></tr>
|
||||
<tr id="row_0_0_2_" class="even"><td class="entry"><span style="width:48px;display:inline-block;"> </span><span class="icona"><span class="icon">C</span></span><a class="el" href="structserif_1_1composition_1_1_global_composition.html" target="_self">GlobalComposition</a></td><td class="desc">Represents the global composition of a system. This tends to be used after finalize and is primarily for internal use </td></tr>
|
||||
<tr id="row_0_1_" class="odd"><td class="entry"><span style="width:16px;display:inline-block;"> </span><span id="arr_0_1_" class="arrow" onclick="dynsection.toggleFolder('0_1_')">▼</span><span class="icona"><span class="icon">N</span></span><a class="el" href="namespaceserif_1_1constant.html" target="_self">constant</a></td><td class="desc"></td></tr>
|
||||
|
||||
@@ -10,12 +10,10 @@
|
||||
<body>
|
||||
<a href="4_d_s_t_a_r_types_8h.html"/>
|
||||
<a href="4_d_s_t_a_r_types_8h_source.html"/>
|
||||
<a href="_2_users_2tboudreaux_2_programming_2_s_e_ri_f_2src_2composition_2public_2composition_8h-example.html"/>
|
||||
<a href="_2_users_2tboudreaux_2_programming_2_s_e_ri_f_2src_2network_2public_2approx8_8h-example.html"/>
|
||||
<a href="_2_users_2tboudreaux_2_programming_2_s_e_ri_f_2src_2polytrope_2solver_2public_2poly_solver_8h-example.html"/>
|
||||
<a href="____init_____8py.html"/>
|
||||
<a href="____init_____8py_source.html"/>
|
||||
<a href="_constructing-example.html"/>
|
||||
<a href="_e_o_sio_8cpp.html"/>
|
||||
<a href="_e_o_sio_8cpp_source.html"/>
|
||||
<a href="_e_o_sio_8h.html"/>
|
||||
|
||||
@@ -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';
|
||||
|
||||
@@ -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]
|
||||
};
|
||||
|
||||
@@ -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}}{}
|
||||
|
||||
@@ -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 ---
|
||||
|
||||
@@ -29,4 +29,4 @@ composition_dep = declare_dependency(
|
||||
)
|
||||
|
||||
# Make headers accessible
|
||||
install_headers(composition_headers, subdir : '4DSSE/composition')
|
||||
install_headers(composition_headers, subdir : 'SERiF/composition')
|
||||
File diff suppressed because it is too large
Load Diff
@@ -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 <iostream>
|
||||
#include <string>
|
||||
@@ -34,6 +33,20 @@
|
||||
#include "atomicSpecies.h"
|
||||
|
||||
namespace serif::composition {
|
||||
struct CanonicalComposition {
|
||||
double X = 0.0; ///< Mass fraction of Hydrogen.
|
||||
double Y = 0.0; ///< Mass fraction of Helium.
|
||||
double Z = 0.0; ///< Mass fraction of Metals.
|
||||
|
||||
friend std::ostream& operator<<(std::ostream& os, const CanonicalComposition& composition) {
|
||||
os << "<CanonicalComposition: "
|
||||
<< "X = " << composition.X << ", "
|
||||
<< "Y = " << composition.Y << ", "
|
||||
<< "Z = " << composition.Z << ">";
|
||||
return os;
|
||||
}
|
||||
};
|
||||
|
||||
/**
|
||||
* @brief Represents the global composition of a system. This tends to be used after finalize and is primarily for internal use.
|
||||
*/
|
||||
@@ -68,7 +81,7 @@ namespace serif::composition {
|
||||
* @brief Constructs a CompositionEntry with the given symbol and mode.
|
||||
* @param symbol The chemical symbol of the species.
|
||||
* @param massFracMode True if mass fraction mode, false if number fraction mode.
|
||||
* @example
|
||||
* *Example Usage:*
|
||||
* @code
|
||||
* CompositionEntry entry("H", true);
|
||||
* @endcode
|
||||
@@ -185,13 +198,13 @@ namespace serif::composition {
|
||||
* - The only exception to the finalize rule is if the compositon was constructed with symbols and mass fractions at instantiation time. In this case, the composition is automatically finalized.
|
||||
* however, this means that the composition passed to the constructor must be valid.
|
||||
*
|
||||
* @example Constructing a finalized composition with symbols and mass fractions:
|
||||
* *Example Usage:* Constructing a finalized composition with symbols and mass fractions:
|
||||
* @code
|
||||
* std::vector<std::string> symbols = {"H", "He"};
|
||||
* std::vector<double> mass_fractions = {0.7, 0.3};
|
||||
* Composition comp(symbols, mass_fractions);
|
||||
* @endcode
|
||||
* @example Constructing a composition with symbols and finalizing it later:
|
||||
* *Example Usage:* Constructing a composition with symbols and finalizing it later:
|
||||
* @code
|
||||
* std::vector<std::string> symbols = {"H", "He"};
|
||||
* Composition comp(symbols);
|
||||
@@ -220,7 +233,7 @@ namespace serif::composition {
|
||||
* @param symbol The symbol to check.
|
||||
* @return True if the symbol is valid, false otherwise.
|
||||
*/
|
||||
bool isValidSymbol(const std::string& symbol) const;
|
||||
static bool isValidSymbol(const std::string& symbol);
|
||||
|
||||
/**
|
||||
* @brief Checks if the given mass fractions are valid.
|
||||
@@ -271,31 +284,31 @@ namespace serif::composition {
|
||||
/**
|
||||
* @brief Constructs a Composition with the given symbols.
|
||||
* @param symbols The symbols to initialize the composition with.
|
||||
* @example
|
||||
* *Example Usage:*
|
||||
* @code
|
||||
* std::vector<std::string> symbols = {"H", "O"};
|
||||
* Composition comp(symbols);
|
||||
* @endcode
|
||||
*/
|
||||
Composition(const std::vector<std::string>& symbols);
|
||||
explicit Composition(const std::vector<std::string>& symbols);
|
||||
|
||||
/**
|
||||
* @brief Constructs a Composition with the given symbols as a set.
|
||||
* @param symbols The symbols to initialize the composition with.
|
||||
* @example
|
||||
* *Example Usage:*
|
||||
* @code
|
||||
* std::set<std::string> symbols = {"H", "O"};
|
||||
* Composition comp(symbols);
|
||||
* @endcode
|
||||
*/
|
||||
Composition(const std::set<std::string>& symbols);
|
||||
explicit Composition(const std::set<std::string>& symbols);
|
||||
|
||||
/**
|
||||
* @brief Constructs a Composition with the given symbols and mass fractions.
|
||||
* @param symbols The symbols to initialize the composition with.
|
||||
* @param mass_fractions The mass fractions corresponding to the symbols.
|
||||
* @param massFracMode True if mass fraction mode, false if number fraction mode.
|
||||
* @example
|
||||
* *Example Usage:*
|
||||
* @code
|
||||
* std::vector<std::string> symbols = {"H", "O"};
|
||||
* std::vector<double> mass_fractions = {0.1, 0.9};
|
||||
@@ -304,11 +317,19 @@ namespace serif::composition {
|
||||
*/
|
||||
Composition(const std::vector<std::string>& symbols, const std::vector<double>& mass_fractions, bool massFracMode=true);
|
||||
|
||||
/**
|
||||
* @brief Constructs a Composition from another Composition.
|
||||
* @param composition The Composition to copy.
|
||||
*/
|
||||
Composition(const Composition& composition);
|
||||
|
||||
Composition& operator=(Composition const& other);
|
||||
|
||||
/**
|
||||
* @brief Registers a new symbol.
|
||||
* @param symbol The symbol to register.
|
||||
* @param massFracMode True if mass fraction mode, false if number fraction mode.
|
||||
* @example
|
||||
* *Example Usage:*
|
||||
* @code
|
||||
* Composition comp;
|
||||
* comp.registerSymbol("H");
|
||||
@@ -320,7 +341,7 @@ namespace serif::composition {
|
||||
* @brief Registers multiple new symbols.
|
||||
* @param symbols The symbols to register.
|
||||
* @param massFracMode True if mass fraction mode, false if number fraction mode.
|
||||
* @example
|
||||
* *Example Usage:*
|
||||
* @code
|
||||
* std::vector<std::string> symbols = {"H", "O"};
|
||||
* Composition comp;
|
||||
@@ -333,14 +354,14 @@ namespace serif::composition {
|
||||
* @brief Gets the registered symbols.
|
||||
* @return A set of registered symbols.
|
||||
*/
|
||||
std::set<std::string> getRegisteredSymbols() const;
|
||||
[[nodiscard]] std::set<std::string> getRegisteredSymbols() const;
|
||||
|
||||
/**
|
||||
* @brief Sets the mass fraction for a given symbol.
|
||||
* @param symbol The symbol to set the mass fraction for.
|
||||
* @param mass_fraction The mass fraction to set.
|
||||
* @return The mass fraction that was set.
|
||||
* @example
|
||||
* *Example Usage:*
|
||||
* @code
|
||||
* Composition comp;
|
||||
* comp.setMassFraction("H", 0.1);
|
||||
@@ -353,7 +374,7 @@ namespace serif::composition {
|
||||
* @param symbols The symbols to set the mass fraction for.
|
||||
* @param mass_fractions The mass fractions corresponding to the symbols.
|
||||
* @return A vector of mass fractions that were set.
|
||||
* @example
|
||||
* *Example Usage:*
|
||||
* @code
|
||||
* std::vector<std::string> symbols = {"H", "O"};
|
||||
* std::vector<double> mass_fractions = {0.1, 0.9};
|
||||
@@ -390,40 +411,52 @@ namespace serif::composition {
|
||||
* @brief Gets the mass fractions of all compositions.
|
||||
* @return An unordered map of compositions with their mass fractions.
|
||||
*/
|
||||
std::unordered_map<std::string, double> getMassFraction() const;
|
||||
[[nodiscard]] std::unordered_map<std::string, double> getMassFraction() const;
|
||||
|
||||
/**
|
||||
* @brief Gets the mass fraction for a given symbol.
|
||||
* @param symbol The symbol to get the mass fraction for.
|
||||
* @return The mass fraction for the given symbol.
|
||||
*/
|
||||
double getMassFraction(const std::string& symbol) const;
|
||||
[[nodiscard]] double getMassFraction(const std::string& symbol) const;
|
||||
|
||||
/**
|
||||
* @brief Gets the number fraction for a given symbol.
|
||||
* @param symbol The symbol to get the number fraction for.
|
||||
* @return The number fraction for the given symbol.
|
||||
*/
|
||||
double getNumberFraction(const std::string& symbol) const;
|
||||
[[nodiscard]] double getNumberFraction(const std::string& symbol) const;
|
||||
|
||||
/**
|
||||
* @brief Gets the number fractions of all compositions.
|
||||
* @return An unordered map of compositions with their number fractions.
|
||||
*/
|
||||
std::unordered_map<std::string, double> getNumberFraction() const;
|
||||
[[nodiscard]] std::unordered_map<std::string, double> getNumberFraction() const;
|
||||
|
||||
/**
|
||||
* @brief Gets the composition entry and global composition for a given symbol.
|
||||
* @param symbol The symbol to get the composition for.
|
||||
* @return A pair containing the CompositionEntry and GlobalComposition for the given symbol.
|
||||
*/
|
||||
std::pair<CompositionEntry, GlobalComposition> getComposition(const std::string& symbol) const;
|
||||
[[nodiscard]] std::pair<CompositionEntry, GlobalComposition> getComposition(const std::string& symbol) const;
|
||||
|
||||
/**
|
||||
* @brief Gets all composition entries and the global composition.
|
||||
* @return A pair containing an unordered map of CompositionEntries and the GlobalComposition.
|
||||
*/
|
||||
std::pair<std::unordered_map<std::string, CompositionEntry>, GlobalComposition> getComposition() const;
|
||||
[[nodiscard]] std::pair<std::unordered_map<std::string, CompositionEntry>, GlobalComposition> getComposition() const;
|
||||
|
||||
/**
|
||||
* @brief Compute the mean particle mass of the composition.
|
||||
* @return Mean particle mass in g.
|
||||
*/
|
||||
[[nodiscard]] double getMeanParticleMass() const;
|
||||
|
||||
/**
|
||||
* @brief Compute the mean atomic mass number of the composition.
|
||||
* @return Mean atomic mass number.
|
||||
*/
|
||||
[[nodiscard]] double getMeanAtomicNumber() const;
|
||||
|
||||
/**
|
||||
* @brief Gets a subset of the composition.
|
||||
@@ -446,6 +479,16 @@ namespace serif::composition {
|
||||
*/
|
||||
void setCompositionMode(bool massFracMode);
|
||||
|
||||
/**
|
||||
* @brief Gets the current canonical composition (X, Y, Z).
|
||||
* @param harsh If true, this will throw an error if X-Y != Z where Z is computed as the sum of all other elements.
|
||||
* @return True if mass fraction mode, false if number fraction mode.
|
||||
*
|
||||
* @throws std::runtime_error if the composition is not finalized or if the canonical composition cannot be computed.
|
||||
* @throws std::runtime_error if harsh is true and the canonical composition is not valid.
|
||||
*/
|
||||
[[nodiscard]] CanonicalComposition getCanonicalComposition(bool harsh=false) const;
|
||||
|
||||
/**
|
||||
* @brief Overloaded output stream operator for Composition.
|
||||
* @param os The output stream.
|
||||
@@ -463,6 +506,4 @@ namespace serif::composition {
|
||||
Composition operator+(const Composition& other) const;
|
||||
|
||||
};
|
||||
}; // namespace serif::composition
|
||||
|
||||
#endif // COMPOSITION_H
|
||||
}; // namespace serif::composition
|
||||
@@ -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,
|
||||
|
||||
64
src/eos/private/EOS.cpp
Normal file
64
src/eos/private/EOS.cpp
Normal file
@@ -0,0 +1,64 @@
|
||||
#include "EOS.h"
|
||||
#include "EOSio.h"
|
||||
#include "helm.h"
|
||||
#include <string>
|
||||
|
||||
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<std::unique_ptr<helmholtz::HELMTable>>(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;
|
||||
}
|
||||
}
|
||||
@@ -28,25 +28,32 @@
|
||||
#include <string>
|
||||
|
||||
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();
|
||||
}
|
||||
}
|
||||
|
||||
@@ -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<std::string>("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;
|
||||
|
||||
266
src/eos/public/EOS.h
Normal file
266
src/eos/public/EOS.h
Normal file
@@ -0,0 +1,266 @@
|
||||
#pragma once
|
||||
|
||||
#include "EOSio.h"
|
||||
#include "helm.h"
|
||||
#include <string>
|
||||
#include "composition.h"
|
||||
#include <iomanip>
|
||||
|
||||
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 << "<EOSInput: "
|
||||
<< "Density: " << input.density << " g/cm^3, "
|
||||
<< "Temperature: " << input.temperature << " K, "
|
||||
<< "Composition: " << input.composition << ">";
|
||||
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) << "<EOSParameter (" << param.name << "): " << param.total << " (gas: " << param.gas
|
||||
<< ", radiation: " << param.radiation << ") "
|
||||
<< "d/dRho: " << param.dDensity << ", d/dT: " << param.dTemperature
|
||||
<< ", d/dAbar: " << param.dMeanAtomicMassNumber
|
||||
<< ", d/dZbar: " << param.dMeanAtomicNumber << ">";
|
||||
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 <iostream>
|
||||
*
|
||||
* 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<std::pair<int, int>, double>& mass_fractions);
|
||||
* // std::map<std::pair<int, int>, 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.
|
||||
};
|
||||
}
|
||||
@@ -22,15 +22,20 @@
|
||||
#include <string>
|
||||
#include <variant>
|
||||
#include <memory>
|
||||
#include <unordered_map>
|
||||
#include "helm.h"
|
||||
|
||||
namespace serif::eos {
|
||||
enum EOSFormat {
|
||||
HELM, ///< Helmholtz EOS format.
|
||||
};
|
||||
|
||||
static inline std::unordered_map<EOSFormat, std::string> FormatStringLookup = {
|
||||
{HELM, "Helmholtz"}
|
||||
};
|
||||
// EOS table format includes
|
||||
|
||||
using EOSTable = std::variant<
|
||||
std::unique_ptr<serif::eos::helmholtz::HELMTable>
|
||||
>;
|
||||
using EOSTable = std::variant<std::unique_ptr<serif::eos::helmholtz::HELMTable>>;
|
||||
|
||||
/**
|
||||
* @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; }
|
||||
};
|
||||
}
|
||||
|
||||
|
||||
@@ -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);
|
||||
|
||||
}
|
||||
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -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<Net::niso; j++) {
|
||||
for (int i=0; i<Net::niso; i++) {
|
||||
J(Net::iener,j) += J(i,j)*Net::mion[i];
|
||||
// energy accounting
|
||||
for (int j=0; j<Approx8Net::nIso; j++) {
|
||||
for (int i=0; i<Approx8Net::nIso; i++) {
|
||||
J(Approx8Net::iEnergy,j) += J(i,j)*Approx8Net::mIon[i];
|
||||
}
|
||||
J(Net::iener,j) *= -avo*clight*clight;
|
||||
J(Approx8Net::iEnergy,j) *= -avo*clight*clight;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void ODE::operator() ( const vector_type &y, vector_type &dydt, double /* t */) {
|
||||
serif::constant::Constants& constants = serif::constant::Constants::getInstance();
|
||||
void ODE::operator() ( const vector_type &y, vector_type &dydt, double /* t */) const {
|
||||
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
|
||||
double T=y[Net::itemp];
|
||||
double den=y[Net::iden];
|
||||
vec7 T9=get_T9_array(T);
|
||||
const double T = y[Approx8Net::iTemp];
|
||||
const double den = y[Approx8Net::iDensity];
|
||||
const vec7 T9=get_T9_array(T);
|
||||
|
||||
// rates
|
||||
double rpp=den*pp_rate(T9);
|
||||
double r33=den*he3he3_rate(T9);
|
||||
double r34=den*he3he4_rate(T9);
|
||||
double r3a=den*den*triple_alpha_rate(T9);
|
||||
double rc12p=den*c12p_rate(T9);
|
||||
double rc12a=den*c12a_rate(T9);
|
||||
double rn14p=den*n14p_rate(T9);
|
||||
double rn14a=n14a_rate(T9);
|
||||
double ro16p=den*o16p_rate(T9);
|
||||
double ro16a=den*o16a_rate(T9);
|
||||
double rne20a=den*ne20a_rate(T9);
|
||||
double r1212=den*c12c12_rate(T9);
|
||||
double r1216=den*c12o16_rate(T9);
|
||||
const double rpp=den*pp_rate(T9);
|
||||
const double r33=den*he3he3_rate(T9);
|
||||
const double r34=den*he3he4_rate(T9);
|
||||
const double r3a=den*den*triple_alpha_rate(T9);
|
||||
const double rc12p=den*c12p_rate(T9);
|
||||
const double rc12a=den*c12a_rate(T9);
|
||||
const double rn14p=den*n14p_rate(T9);
|
||||
const double rn14a=n14a_rate(T9);
|
||||
const double ro16p=den*o16p_rate(T9);
|
||||
const double ro16a=den*o16a_rate(T9);
|
||||
const double rne20a=den*ne20a_rate(T9);
|
||||
const double r1212=den*c12c12_rate(T9);
|
||||
const double r1216=den*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];
|
||||
|
||||
dydt[Net::ih1] = -1.5*yh1*yh1*rpp;
|
||||
dydt[Net::ih1] += yhe3*yhe3*r33;
|
||||
dydt[Net::ih1] += -yhe3*yhe4*r34;
|
||||
dydt[Net::ih1] += -2*yh1*yc12*rc12p;
|
||||
dydt[Net::ih1] += -2*yh1*yn14*rn14p;
|
||||
dydt[Net::ih1] += -2*yh1*yo16*ro16p;
|
||||
|
||||
dydt[Net::ihe3] = 0.5*yh1*yh1*rpp;
|
||||
dydt[Net::ihe3] += -yhe3*yhe3*r33;
|
||||
dydt[Net::ihe3] += -yhe3*yhe4*r34;
|
||||
|
||||
dydt[Net::ihe4] = 0.5*yhe3*yhe3*r33;
|
||||
dydt[Net::ihe4] += yhe3*yhe4*r34;
|
||||
dydt[Net::ihe4] += -yhe4*yc12*rc12a;
|
||||
dydt[Net::ihe4] += yh1*yn14*afrac*rn14p;
|
||||
dydt[Net::ihe4] += yh1*yo16*ro16p;
|
||||
dydt[Net::ihe4] += -0.5*yhe4*yhe4*yhe4*r3a;
|
||||
dydt[Net::ihe4] += -yhe4*yo16*ro16a;
|
||||
dydt[Net::ihe4] += 0.5*yc12*yc12*r1212;
|
||||
dydt[Net::ihe4] += yc12*yo16*r1216;
|
||||
dydt[Net::ihe4] += -yhe4*yne20*rne20a;
|
||||
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];
|
||||
|
||||
dydt[Net::ic12] = (1./6.)*yhe4*yhe4*yhe4*r3a;
|
||||
dydt[Net::ic12] += -yhe4*yc12*rc12a;
|
||||
dydt[Net::ic12] += -yh1*yc12*rc12p;
|
||||
dydt[Net::ic12] += yh1*yn14*afrac*rn14p;
|
||||
dydt[Net::ic12] += -yc12*yc12*r1212;
|
||||
dydt[Net::ic12] += -yc12*yo16*r1216;
|
||||
dydt[Approx8Net::ih1] = -1.5*yh1*yh1*rpp;
|
||||
dydt[Approx8Net::ih1] += yhe3*yhe3*r33;
|
||||
dydt[Approx8Net::ih1] += -yhe3*yhe4*r34;
|
||||
dydt[Approx8Net::ih1] += -2*yh1*yc12*rc12p;
|
||||
dydt[Approx8Net::ih1] += -2*yh1*yn14*rn14p;
|
||||
dydt[Approx8Net::ih1] += -2*yh1*yo16*ro16p;
|
||||
|
||||
dydt[Net::in14] = yh1*yc12*rc12p;
|
||||
dydt[Net::in14] += -yh1*yn14*rn14p;
|
||||
dydt[Net::in14] += yh1*yo16*ro16p;
|
||||
dydt[Net::in14] += -yhe4*yn14*rn14a;
|
||||
|
||||
dydt[Net::io16] = yhe4*yc12*rc12a;
|
||||
dydt[Net::io16] += yh1*yn14*pfrac*rn14p;
|
||||
dydt[Net::io16] += -yh1*yo16*ro16p;
|
||||
dydt[Net::io16] += -yc12*yo16*r1216;
|
||||
dydt[Net::io16] += -yhe4*yo16*ro16a;
|
||||
dydt[Approx8Net::ihe3] = 0.5*yh1*yh1*rpp;
|
||||
dydt[Approx8Net::ihe3] += -yhe3*yhe3*r33;
|
||||
dydt[Approx8Net::ihe3] += -yhe3*yhe4*r34;
|
||||
|
||||
dydt[Net::ine20] = 0.5*yc12*yc12*r1212;
|
||||
dydt[Net::ine20] += yhe4*yn14*rn14a;
|
||||
dydt[Net::ine20] += yhe4*yo16*ro16a;
|
||||
dydt[Net::ine20] += -yhe4*yne20*rne20a;
|
||||
dydt[Approx8Net::ihe4] = 0.5*yhe3*yhe3*r33;
|
||||
dydt[Approx8Net::ihe4] += yhe3*yhe4*r34;
|
||||
dydt[Approx8Net::ihe4] += -yhe4*yc12*rc12a;
|
||||
dydt[Approx8Net::ihe4] += yh1*yn14*aFrac*rn14p;
|
||||
dydt[Approx8Net::ihe4] += yh1*yo16*ro16p;
|
||||
dydt[Approx8Net::ihe4] += -0.5*yhe4*yhe4*yhe4*r3a;
|
||||
dydt[Approx8Net::ihe4] += -yhe4*yo16*ro16a;
|
||||
dydt[Approx8Net::ihe4] += 0.5*yc12*yc12*r1212;
|
||||
dydt[Approx8Net::ihe4] += yc12*yo16*r1216;
|
||||
dydt[Approx8Net::ihe4] += -yhe4*yne20*rne20a;
|
||||
|
||||
dydt[Net::img24] = yc12*yo16*r1216;
|
||||
dydt[Net::img24] += yhe4*yne20*rne20a;
|
||||
dydt[Approx8Net::ic12] = (1./6.)*yhe4*yhe4*yhe4*r3a;
|
||||
dydt[Approx8Net::ic12] += -yhe4*yc12*rc12a;
|
||||
dydt[Approx8Net::ic12] += -yh1*yc12*rc12p;
|
||||
dydt[Approx8Net::ic12] += yh1*yn14*aFrac*rn14p;
|
||||
dydt[Approx8Net::ic12] += -yc12*yc12*r1212;
|
||||
dydt[Approx8Net::ic12] += -yc12*yo16*r1216;
|
||||
|
||||
dydt[Net::itemp] = 0.;
|
||||
dydt[Net::iden] = 0.;
|
||||
dydt[Approx8Net::in14] = yh1*yc12*rc12p;
|
||||
dydt[Approx8Net::in14] += -yh1*yn14*rn14p;
|
||||
dydt[Approx8Net::in14] += yh1*yo16*ro16p;
|
||||
dydt[Approx8Net::in14] += -yhe4*yn14*rn14a;
|
||||
|
||||
dydt[Approx8Net::io16] = yhe4*yc12*rc12a;
|
||||
dydt[Approx8Net::io16] += yh1*yn14*pFrac*rn14p;
|
||||
dydt[Approx8Net::io16] += -yh1*yo16*ro16p;
|
||||
dydt[Approx8Net::io16] += -yc12*yo16*r1216;
|
||||
dydt[Approx8Net::io16] += -yhe4*yo16*ro16a;
|
||||
|
||||
dydt[Approx8Net::ine20] = 0.5*yc12*yc12*r1212;
|
||||
dydt[Approx8Net::ine20] += yhe4*yn14*rn14a;
|
||||
dydt[Approx8Net::ine20] += yhe4*yo16*ro16a;
|
||||
dydt[Approx8Net::ine20] += -yhe4*yne20*rne20a;
|
||||
|
||||
dydt[Approx8Net::img24] = yc12*yo16*r1216;
|
||||
dydt[Approx8Net::img24] += yhe4*yne20*rne20a;
|
||||
|
||||
dydt[Approx8Net::iTemp] = 0.;
|
||||
dydt[Approx8Net::iDensity] = 0.;
|
||||
|
||||
// energy accounting
|
||||
double enuc = 0.;
|
||||
for (int i=0; i<Net::niso; i++) {
|
||||
enuc += Net::mion[i]*dydt[i];
|
||||
double eNuc = 0.;
|
||||
for (int i=0; i<Approx8Net::nIso; i++) {
|
||||
eNuc += Approx8Net::mIon[i]*dydt[i];
|
||||
}
|
||||
dydt[Net::iener] = -enuc*avo*clight*clight;
|
||||
}
|
||||
dydt[Approx8Net::iEnergy] = -eNuc*avo*clight*clight;
|
||||
}
|
||||
|
||||
Approx8Network::Approx8Network() : Network(APPROX8) {}
|
||||
|
||||
NetOut Approx8Network::evaluate(const NetIn &netIn) {
|
||||
m_y = convert_netIn(netIn);
|
||||
m_tmax = netIn.tmax;
|
||||
m_tMax = netIn.tMax;
|
||||
m_dt0 = netIn.dt0;
|
||||
|
||||
const double stiff_abs_tol = m_config.get<double>("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<double> 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<std::string> 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;
|
||||
|
||||
@@ -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;
|
||||
}
|
||||
}
|
||||
|
||||
@@ -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<double,7> 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<double,7> 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<int,niso> aion = {
|
||||
static constexpr std::array<int,nIso> aIon = {
|
||||
1,
|
||||
3,
|
||||
4,
|
||||
@@ -89,7 +90,7 @@ namespace serif::network::approx8{
|
||||
24
|
||||
};
|
||||
|
||||
static constexpr std::array<double,niso> mion = {
|
||||
static constexpr std::array<double,nIso> 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
|
||||
|
||||
@@ -25,15 +25,27 @@
|
||||
#include "probe.h"
|
||||
#include "config.h"
|
||||
#include "quill/Logger.h"
|
||||
#include "composition.h"
|
||||
#include <unordered_map>
|
||||
|
||||
namespace serif::network {
|
||||
|
||||
enum NetworkFormat {
|
||||
APPROX8, ///< Approx8 nuclear reaction network format.
|
||||
UNKNOWN,
|
||||
};
|
||||
|
||||
static inline std::unordered_map<NetworkFormat, std::string> 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<double> 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<double> 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
|
||||
|
||||
@@ -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 "<EOSio(filename='" + eos.getFilename() + "', format='" + eos.getFormat() + "')>";
|
||||
return "<EOSio(filename='" + eos.getFilename() + "', format='" + eos.getFormatName() + "')>";
|
||||
});
|
||||
|
||||
py::class_<serif::eos::EOSTable>(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_<serif::eos::helmholtz::EOS>(eos_submodule, "EOS")
|
||||
py::class_<serif::eos::helmholtz::HELMEOSOutput>(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 "<EOS (output from helmholtz eos)>";
|
||||
});
|
||||
|
||||
py::class_<serif::eos::helmholtz::EOSInput>(eos_submodule, "EOSInput")
|
||||
py::class_<serif::eos::helmholtz::HELMEOSInput>(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 "<EOSInput(T=" + std::to_string(input.T) +
|
||||
.def_readwrite("T", &serif::eos::helmholtz::HELMEOSInput::T)
|
||||
.def_readwrite("rho", &serif::eos::helmholtz::HELMEOSInput::rho)
|
||||
.def_readwrite("abar", &serif::eos::helmholtz::HELMEOSInput::abar)
|
||||
.def_readwrite("zbar", &serif::eos::helmholtz::HELMEOSInput::zbar)
|
||||
.def("__repr__", [](const serif::eos::helmholtz::HELMEOSInput &input) {
|
||||
return "<HELMEOSInput(T=" + std::to_string(input.T) +
|
||||
", rho=" + std::to_string(input.rho) +
|
||||
", abar=" + std::to_string(input.abar) +
|
||||
", zbar=" + std::to_string(input.zbar) + ")>";
|
||||
|
||||
@@ -1,5 +1,4 @@
|
||||
#include <gtest/gtest.h>
|
||||
#include <iostream>
|
||||
#include <memory>
|
||||
#include <string>
|
||||
#include <sstream>
|
||||
@@ -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<std::unique_ptr<serif::eos::EOSio>>(rm.getResource("eos:helm"));
|
||||
auto& table = eos->getTable();
|
||||
auto& helmTable = *std::get<std::unique_ptr<serif::eos::helmholtz::HELMTable>>(table);
|
||||
const auto& table = eos->getTable();
|
||||
const auto& helmTable = *std::get<std::unique_ptr<serif::eos::helmholtz::HELMTable>>(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<nel; i++) {
|
||||
asum += xmass[i]/aion[i];
|
||||
zsum += xmass[i]*zion[i]/aion[i];
|
||||
aSum += xMass[i]/aIon[i];
|
||||
zSum += xMass[i]*zIon[i]/aIon[i];
|
||||
}
|
||||
eos1.abar = 1.0/asum;
|
||||
eos1.zbar = eos1.abar*zsum;
|
||||
eos1.abar = 1.0/aSum;
|
||||
eos1.zbar = eos1.abar*zSum;
|
||||
|
||||
serif::resource::ResourceManager& rm = serif::resource::ResourceManager::getInstance();
|
||||
auto& eos = std::get<std::unique_ptr<serif::eos::EOSio>>(rm.getResource("eos:helm"));
|
||||
auto& table = eos->getTable();
|
||||
auto& helmTable = *std::get<std::unique_ptr<serif::eos::helmholtz::HELMTable>>(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<std::unique_ptr<serif::eos::EOSio>>(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);
|
||||
}
|
||||
|
||||
@@ -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
|
||||
)
|
||||
|
||||
|
||||
@@ -4,6 +4,7 @@
|
||||
#include "approx8.h"
|
||||
#include "config.h"
|
||||
#include "network.h"
|
||||
#include "composition.h"
|
||||
|
||||
#include <vector>
|
||||
|
||||
@@ -32,19 +33,31 @@ TEST_F(approx8Test, evaluate) {
|
||||
serif::network::NetIn netIn;
|
||||
|
||||
std::vector<double> comp = {0.708, 2.94e-5, 0.276, 0.003, 0.0011, 9.62e-3, 1.62e-3, 5.16e-4};
|
||||
std::vector<std::string> 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);
|
||||
}
|
||||
|
||||
@@ -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
|
||||
|
||||
@@ -48,7 +48,7 @@ TEST_F(resourceManagerTest, getResource) {
|
||||
const serif::resource::types::Resource &r = rm.getResource(name);
|
||||
// BREAKPOINT();
|
||||
const auto &eos = std::get<std::unique_ptr<serif::eos::EOSio>>(r);
|
||||
EXPECT_EQ("helm", eos->getFormat());
|
||||
EXPECT_EQ("Helmholtz", eos->getFormatName());
|
||||
serif::eos::EOSTable &table = eos->getTable();
|
||||
|
||||
// -- Extract the Helm table from the EOSTable
|
||||
|
||||
Reference in New Issue
Block a user