From 36da29cfa0d46067c7708f4abde7e52f975b4c20 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Thu, 14 Apr 2022 10:38:28 -0600 Subject: [PATCH 01/34] Initial sap/lap ramp pressure functions. --- singularity-eos/eos/modifiers/ramps_eos.hpp | 163 ++++++++++++++++++++ 1 file changed, 163 insertions(+) create mode 100644 singularity-eos/eos/modifiers/ramps_eos.hpp diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp new file mode 100644 index 0000000000..f20e37cee9 --- /dev/null +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -0,0 +1,163 @@ +//------------------------------------------------------------------------------ +// © 2021-2022. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#ifndef _SINGULARITY_EOS_EOS_RAMPS_EOS_ +#define _SINGULARITY_EOS_EOS_RAMPS_EOS_ + +#include "stdio.h" +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +namespace singularity { + +using namespace eos_base; + +template +class SAPRampEOS : public EosBase> { + public: + // move semantics ensures dynamic memory comes along for the ride + SAPRampEOS(T &&t, const Real shift) : t_(std::forward(t)), shift_(shift) {} + SAPRampEOS() = default; + + auto GetOnDevice() { return SAPRampEOS(t_.GetOnDevice(), shift_); } + inline void Finalize() { t_.Finalize(); } + + PORTABLE_INLINE_FUNCTION + Real get_ramp_pressure(Real rho) const { + const Real p_ramp {rho < rholow ? 0.0 : + rho < rhomid ? a*(rho/r0 - 1.0) : + b*(rho/r0 - c)}; + return p_ramp; + } + + PORTABLE_FUNCTION + Real TemperatureFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const { + return t_.TemperatureFromDensityInternalEnergy(rho, sie, lambda); + } + PORTABLE_FUNCTION + Real InternalEnergyFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const { + return t_.InternalEnergyFromDensityTemperature(rho, temperature, lambda); + } + PORTABLE_FUNCTION + Real PressureFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const { + // ramp pressure + const Real p_ramp {get_ramp_pressure(rho)}; + // pressure from eos + const Real p_eos {t_.PressureFromDensityInternalEnergy(rho, sie, lambda)}; + // return max(p_ramp, p_eos) + return p_eos < p_ramp ? p_ramp : p_eos; + } + PORTABLE_FUNCTION + Real SpecificHeatFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const { + return t_.SpecificHeatFromDensityInternalEnergy(rho, sie, lambda); + } + PORTABLE_FUNCTION + Real BulkModulusFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const { + return t_.BulkModulusFromDensityInternalEnergy(rho, sie, lambda); + } + PORTABLE_FUNCTION + Real GruneisenParamFromDensityInternalEnergy(const Real rho, const Real sie, + Real *lambda = nullptr) const { + return t_.GruneisenParamFromDensityInternalEnergy(rho, sie, lambda); + } + PORTABLE_FUNCTION + Real PressureFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const { + const Real p_ramp {get_ramp_pressure(rho)}; + const Real p_eos {t_.PressureFromDensityTemperature(rho, temperature, lambda)}; + return p_eos < p_ramp ? p_ramp : p_eos; + } + PORTABLE_FUNCTION + Real SpecificHeatFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const { + return t_.SpecificHeatFromDensityTemperature(rho, temperature, lambda); + } + PORTABLE_FUNCTION + Real BulkModulusFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const { + return t_.BulkModulusFromDensityTemperature(rho, temperature, lambda); + } + PORTABLE_FUNCTION + Real GruneisenParamFromDensityTemperature(const Real rho, const Real temperature, + Real *lambda = nullptr) const { + return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda); + } + PORTABLE_FUNCTION + void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, + const unsigned long output, Real *lambda = nullptr) const { + Real senergy; + switch (t_.PreferredInput()) { + case thermalqs::density | thermalqs::temperature: + t_.FillEos(rho, temp, energy, press, cv, bmod, output, lambda); + energy = energy + shift_; + break; + case thermalqs::density | thermalqs::specific_internal_energy: + senergy = energy - shift_; + t_.FillEos(rho, temp, senergy, press, cv, bmod, output, lambda); + break; + default: + EOS_ERROR("Didn't find a valid input for SAPRampEOS::FillEOS\n"); + } + } + + PORTABLE_INLINE_FUNCTION + int nlambda() const noexcept { return t_.nlambda(); } + + PORTABLE_FUNCTION + unsigned long PreferredInput() const { return t_.PreferredInput(); } + + PORTABLE_FUNCTION void PrintParams() const { + t_.PrintParams(); + printf("scaling_ratio = %f\n", shift_); + } + PORTABLE_FUNCTION + void DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Real *lambda, Real &rho, Real &sie) const { + t_.DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie); + sie = sie + shift_; + } + + PORTABLE_FUNCTION + void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, + Real &bmod, Real &dpde, Real &dvdt, + Real *lambda = nullptr) const { + t_.ValuesAtReferenceState(rho, temp, sie, press, cv, bmod, dpde, dvdt, lambda); + sie += shift_; + } + + // Vector functions that overload the scalar versions declared here. + SG_ADD_BASE_CLASS_USINGS(SAPRampEOS) + + private: + T t_; + double shift_; +}; + +} // namespace singularity + +#endif // _SINGULARITY_EOS_EOS_RAMPS_EOS_ From c35a76d758510f09d343a33fe1957258372920da Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Wed, 27 Apr 2022 11:08:12 -0600 Subject: [PATCH 02/34] Finished draft implementation of the sap ramp. Naming may need to change since it it taken from lap. --- singularity-eos/eos/modifiers/ramps_eos.hpp | 53 +++++++++++++-------- 1 file changed, 33 insertions(+), 20 deletions(-) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index f20e37cee9..07f20d9a49 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -36,7 +36,14 @@ template class SAPRampEOS : public EosBase> { public: // move semantics ensures dynamic memory comes along for the ride - SAPRampEOS(T &&t, const Real shift) : t_(std::forward(t)), shift_(shift) {} + SAPRampEOS(T &&t, const Real r0, const Real a, const Real b, const Real c) + : t_(std::forward(t)), + r0_(r0), + a_(a), + b_(b), + c_(c), + rmid_(r0*(a-b*c)/(a-b)) + {} SAPRampEOS() = default; auto GetOnDevice() { return SAPRampEOS(t_.GetOnDevice(), shift_); } @@ -44,9 +51,9 @@ class SAPRampEOS : public EosBase> { PORTABLE_INLINE_FUNCTION Real get_ramp_pressure(Real rho) const { - const Real p_ramp {rho < rholow ? 0.0 : - rho < rhomid ? a*(rho/r0 - 1.0) : - b*(rho/r0 - c)}; + const Real p_ramp {rho < r0_ ? 0.0 : + rho < rmid_ ? a_*(rho/r0_ - 1.0) : + b_*(rho/r0_ - c_)}; return p_ramp; } @@ -110,19 +117,23 @@ class SAPRampEOS : public EosBase> { PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, const unsigned long output, Real *lambda = nullptr) const { - Real senergy; - switch (t_.PreferredInput()) { - case thermalqs::density | thermalqs::temperature: - t_.FillEos(rho, temp, energy, press, cv, bmod, output, lambda); - energy = energy + shift_; - break; - case thermalqs::density | thermalqs::specific_internal_energy: - senergy = energy - shift_; - t_.FillEos(rho, temp, senergy, press, cv, bmod, output, lambda); - break; - default: - EOS_ERROR("Didn't find a valid input for SAPRampEOS::FillEOS\n"); + // density must be an input + assert(!(output & thermalqs::density)); + // output passed into internal filleos can't include pressure + const unsigned long ramp_out = output & ~thermalqs::pressure; + // if pressure is output, calculate it first + if (output & thermalqs::pressure) { + // maybe switch case on preferred input and check for not output of the input + // for now sie lookup is prioritized + if (!(output & thermalqs::specific_internal_energy)) { + press = this->PressureFromDensityInternalEnergy(rho, energy, lambda); + } + else if (!(output & thermalqs::temperature)) { + press = this->PressureFromDensityTemperature(rho, temp, lambda); + } } + // call internal filleos + t_.FillEos(rho, temp, energy, press, cv, bmod, ramp_out, lambda); } PORTABLE_INLINE_FUNCTION @@ -133,13 +144,12 @@ class SAPRampEOS : public EosBase> { PORTABLE_FUNCTION void PrintParams() const { t_.PrintParams(); - printf("scaling_ratio = %f\n", shift_); + printf("r0=%e\na=%e\nb=%e\nc=%e\n", r0_, a_, b_, c_); } PORTABLE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, Real *lambda, Real &rho, Real &sie) const { t_.DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie); - sie = sie + shift_; } PORTABLE_FUNCTION @@ -147,7 +157,6 @@ class SAPRampEOS : public EosBase> { Real &bmod, Real &dpde, Real &dvdt, Real *lambda = nullptr) const { t_.ValuesAtReferenceState(rho, temp, sie, press, cv, bmod, dpde, dvdt, lambda); - sie += shift_; } // Vector functions that overload the scalar versions declared here. @@ -155,7 +164,11 @@ class SAPRampEOS : public EosBase> { private: T t_; - double shift_; + Real r0_; + Real a_; + Real b_; + Real c_; + Real rmid_; }; } // namespace singularity From b4641080b47219e7285643886620a5bfdb4d085c Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Wed, 27 Apr 2022 11:11:31 -0600 Subject: [PATCH 03/34] Fix get on device func. --- singularity-eos/eos/modifiers/ramps_eos.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 07f20d9a49..7b4e3ec47f 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -46,7 +46,7 @@ class SAPRampEOS : public EosBase> { {} SAPRampEOS() = default; - auto GetOnDevice() { return SAPRampEOS(t_.GetOnDevice(), shift_); } + auto GetOnDevice() { return SAPRampEOS(t_.GetOnDevice(), r0_, a_, b_, c_); } inline void Finalize() { t_.Finalize(); } PORTABLE_INLINE_FUNCTION From 1e025e69db343c67568742a1a68f65e59e509f3c Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Thu, 28 Apr 2022 18:35:50 -0600 Subject: [PATCH 04/34] Adding test to ramp eos to make sure it is reflexive. --- singularity-eos/eos/CMakeLists.txt | 1 + singularity-eos/eos/eos.hpp | 1 + singularity-eos/eos/modifiers/ramps_eos.hpp | 8 +++++++- singularity-eos/eos/modifiers/shifted_eos.hpp | 2 +- test/eos_unit_tests.cpp | 6 +++++- 5 files changed, 15 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/CMakeLists.txt b/singularity-eos/eos/CMakeLists.txt index e10efa0345..b9e4973f3d 100644 --- a/singularity-eos/eos/CMakeLists.txt +++ b/singularity-eos/eos/CMakeLists.txt @@ -23,6 +23,7 @@ set(EOS_SRCS modifiers/scaled_eos.hpp modifiers/shifted_eos.hpp modifiers/relativistic_eos.hpp + modifiers/ramps_eos.hpp eos_spiner.cpp eos_jwl.cpp eos_builder.cpp diff --git a/singularity-eos/eos/eos.hpp b/singularity-eos/eos/eos.hpp index b73761f5f1..ca590a28a0 100644 --- a/singularity-eos/eos/eos.hpp +++ b/singularity-eos/eos/eos.hpp @@ -42,6 +42,7 @@ #include #include #include +#include #ifdef SINGULARITY_USE_EOSPAC #include diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 7b4e3ec47f..455d197daf 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -43,7 +43,13 @@ class SAPRampEOS : public EosBase> { b_(b), c_(c), rmid_(r0*(a-b*c)/(a-b)) - {} + { + // add input parameter checks to ensure validity of the ramp + assert(r0 >= 0.0); + assert(a > 0.0); + assert(b >= 0); + assert(a != b); + } SAPRampEOS() = default; auto GetOnDevice() { return SAPRampEOS(t_.GetOnDevice(), r0_, a_, b_, c_); } diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index 4315970503..3adb70bb26 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -171,4 +171,4 @@ class ShiftedEOS : public EosBase> { } // namespace singularity -#endif _SINGULARITY_EOS_EOS_SHIFTED_EOS_ +#endif // _SINGULARITY_EOS_EOS_SHIFTED_EOS_ diff --git a/test/eos_unit_tests.cpp b/test/eos_unit_tests.cpp index 0ea4586768..a9b34310dd 100644 --- a/test/eos_unit_tests.cpp +++ b/test/eos_unit_tests.cpp @@ -34,6 +34,7 @@ using singularity::EOS; using singularity::IdealGas; using singularity::ScaledEOS; using singularity::ShiftedEOS; +using singularity::SAPRampEOS; #ifdef SPINER_USE_HDF using singularity::SpinerEOSDependsRhoSie; @@ -80,7 +81,8 @@ inline void array_compare(int num, X &&x, Y &&y, Z &&z, ZT &&ztrue, XN xname, YN } } -inline void compare_two_eoss(const EOS &test_e, const EOS &ref_e) { +template +inline void compare_two_eoss(E1&& test_e, E2&& ref_e) { // compare all individual member functions with 1 as inputs, // this function is meant to catch mis-implementations of // modifiers that can be initialized in such a way as to @@ -240,9 +242,11 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { EOS ig = IdealGas(gm1, Cv); EOS igsh = ScaledEOS(IdealGas(gm1, Cv), 1.0); EOS igsc = ShiftedEOS(IdealGas(gm1, Cv), 0.0); + auto igra = SAPRampEOS(IdealGas(gm1, Cv), 1.e9, 1.0, 1.0, 1.0); THEN("The modified EOS should produce equivalent results") { compare_two_eoss(igsh, ig); compare_two_eoss(igsc, ig); + compare_two_eoss(igra, ig); } } } From 1e02abacc38f3d2c85115e177f6d697413099200 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Fri, 29 Apr 2022 11:01:10 -0600 Subject: [PATCH 05/34] fix assertion that a cannot equal b. --- test/eos_unit_tests.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/eos_unit_tests.cpp b/test/eos_unit_tests.cpp index a9b34310dd..5c3bd0b26b 100644 --- a/test/eos_unit_tests.cpp +++ b/test/eos_unit_tests.cpp @@ -242,7 +242,7 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { EOS ig = IdealGas(gm1, Cv); EOS igsh = ScaledEOS(IdealGas(gm1, Cv), 1.0); EOS igsc = ShiftedEOS(IdealGas(gm1, Cv), 0.0); - auto igra = SAPRampEOS(IdealGas(gm1, Cv), 1.e9, 1.0, 1.0, 1.0); + auto igra = SAPRampEOS(IdealGas(gm1, Cv), 1.e9, 1.0, 2.0, 1.0); THEN("The modified EOS should produce equivalent results") { compare_two_eoss(igsh, ig); compare_two_eoss(igsc, ig); From 5b036544e4c782a89f0441037c1f685f016ac1d3 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Fri, 29 Apr 2022 12:21:35 -0600 Subject: [PATCH 06/34] Fix formatting. --- singularity-eos/eos/eos.hpp | 2 +- singularity-eos/eos/modifiers/ramps_eos.hpp | 32 +++++++++------------ test/eos_unit_tests.cpp | 8 +++--- 3 files changed, 18 insertions(+), 24 deletions(-) diff --git a/singularity-eos/eos/eos.hpp b/singularity-eos/eos/eos.hpp index ca590a28a0..766176b7b4 100644 --- a/singularity-eos/eos/eos.hpp +++ b/singularity-eos/eos/eos.hpp @@ -39,10 +39,10 @@ #include #include #include +#include #include #include #include -#include #ifdef SINGULARITY_USE_EOSPAC #include diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 455d197daf..ad2c587718 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -36,14 +36,9 @@ template class SAPRampEOS : public EosBase> { public: // move semantics ensures dynamic memory comes along for the ride - SAPRampEOS(T &&t, const Real r0, const Real a, const Real b, const Real c) - : t_(std::forward(t)), - r0_(r0), - a_(a), - b_(b), - c_(c), - rmid_(r0*(a-b*c)/(a-b)) - { + SAPRampEOS(T &&t, const Real r0, const Real a, const Real b, const Real c) + : t_(std::forward(t)), r0_(r0), a_(a), b_(b), c_(c), + rmid_(r0 * (a - b * c) / (a - b)) { // add input parameter checks to ensure validity of the ramp assert(r0 >= 0.0); assert(a > 0.0); @@ -57,9 +52,9 @@ class SAPRampEOS : public EosBase> { PORTABLE_INLINE_FUNCTION Real get_ramp_pressure(Real rho) const { - const Real p_ramp {rho < r0_ ? 0.0 : - rho < rmid_ ? a_*(rho/r0_ - 1.0) : - b_*(rho/r0_ - c_)}; + const Real p_ramp{rho < r0_ ? 0.0 + : rho < rmid_ ? a_ * (rho / r0_ - 1.0) + : b_ * (rho / r0_ - c_)}; return p_ramp; } @@ -77,9 +72,9 @@ class SAPRampEOS : public EosBase> { Real PressureFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda = nullptr) const { // ramp pressure - const Real p_ramp {get_ramp_pressure(rho)}; + const Real p_ramp{get_ramp_pressure(rho)}; // pressure from eos - const Real p_eos {t_.PressureFromDensityInternalEnergy(rho, sie, lambda)}; + const Real p_eos{t_.PressureFromDensityInternalEnergy(rho, sie, lambda)}; // return max(p_ramp, p_eos) return p_eos < p_ramp ? p_ramp : p_eos; } @@ -101,8 +96,8 @@ class SAPRampEOS : public EosBase> { PORTABLE_FUNCTION Real PressureFromDensityTemperature(const Real rho, const Real temperature, Real *lambda = nullptr) const { - const Real p_ramp {get_ramp_pressure(rho)}; - const Real p_eos {t_.PressureFromDensityTemperature(rho, temperature, lambda)}; + const Real p_ramp{get_ramp_pressure(rho)}; + const Real p_eos{t_.PressureFromDensityTemperature(rho, temperature, lambda)}; return p_eos < p_ramp ? p_ramp : p_eos; } PORTABLE_FUNCTION @@ -132,10 +127,9 @@ class SAPRampEOS : public EosBase> { // maybe switch case on preferred input and check for not output of the input // for now sie lookup is prioritized if (!(output & thermalqs::specific_internal_energy)) { - press = this->PressureFromDensityInternalEnergy(rho, energy, lambda); - } - else if (!(output & thermalqs::temperature)) { - press = this->PressureFromDensityTemperature(rho, temp, lambda); + press = this->PressureFromDensityInternalEnergy(rho, energy, lambda); + } else if (!(output & thermalqs::temperature)) { + press = this->PressureFromDensityTemperature(rho, temp, lambda); } } // call internal filleos diff --git a/test/eos_unit_tests.cpp b/test/eos_unit_tests.cpp index 5c3bd0b26b..e058080759 100644 --- a/test/eos_unit_tests.cpp +++ b/test/eos_unit_tests.cpp @@ -32,9 +32,9 @@ using singularity::EOS; using singularity::IdealGas; +using singularity::SAPRampEOS; using singularity::ScaledEOS; using singularity::ShiftedEOS; -using singularity::SAPRampEOS; #ifdef SPINER_USE_HDF using singularity::SpinerEOSDependsRhoSie; @@ -81,8 +81,8 @@ inline void array_compare(int num, X &&x, Y &&y, Z &&z, ZT &&ztrue, XN xname, YN } } -template -inline void compare_two_eoss(E1&& test_e, E2&& ref_e) { +template +inline void compare_two_eoss(E1 &&test_e, E2 &&ref_e) { // compare all individual member functions with 1 as inputs, // this function is meant to catch mis-implementations of // modifiers that can be initialized in such a way as to @@ -246,7 +246,7 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { THEN("The modified EOS should produce equivalent results") { compare_two_eoss(igsh, ig); compare_two_eoss(igsc, ig); - compare_two_eoss(igra, ig); + compare_two_eoss(igra, ig); } } } From c049e689718f5137f46d9893c01c8137269f452d Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 17 May 2022 11:02:25 -0600 Subject: [PATCH 07/34] Create functions necessary to initialize a ramped eos. --- singularity-eos/eos/eos.hpp | 8 +++++- singularity-eos/eos/eos_builder.hpp | 33 +++++++++++++++++++++++++ singularity-eos/eos/singularity_eos.cpp | 10 +++++--- 3 files changed, 46 insertions(+), 5 deletions(-) diff --git a/singularity-eos/eos/eos.hpp b/singularity-eos/eos/eos.hpp index c2f27d347e..248203e7a5 100644 --- a/singularity-eos/eos/eos.hpp +++ b/singularity-eos/eos/eos.hpp @@ -1272,8 +1272,14 @@ static constexpr const auto scaled_of_shifted = static constexpr const auto unit_or_rel = transform_variadic_list(partial_eos_list, apply_to_partial); // create combined list -static constexpr const auto combined_list = singularity::detail::concat( +static constexpr const auto combined_list_1 = singularity::detail::concat( full_eos_list, shifted, scaled, scaled_of_shifted, unit_or_rel); +// make a ramped eos of everything +static constexpr const auto ramped_all = + transform_variadic_list(combined_list_1, al{}); +// final combined list +static constexpr const auto combined_list = singularity::detail::concat( + full_eos_list, shifted, scaled, scaled_of_shifted, unit_or_rel, ramped_all); // a function that returns a Variant from a typelist template struct tl_to_Variant_struct { diff --git a/singularity-eos/eos/eos_builder.hpp b/singularity-eos/eos/eos_builder.hpp index 62f3d29aeb..9eecd5372e 100644 --- a/singularity-eos/eos/eos_builder.hpp +++ b/singularity-eos/eos/eos_builder.hpp @@ -83,6 +83,11 @@ EOS makeRelativistic(T &&eos, Real cl) { return RelativisticEOS(std::move(eos), cl); } +template +EOS makeSAPRamp(T &&eos, Real r0, Real a, Real b, Real c) { + return SAPRampEOS(std::forward(eos), r0, a, b, c); +} + template EOS applyShiftAndScale(T &&eos, bool scaled, bool shifted, Real scale, Real shift) { if (shifted && scaled) { @@ -98,6 +103,34 @@ EOS applyShiftAndScale(T &&eos, bool scaled, bool shifted, Real scale, Real shif } return eos; } + +template class W, typename ... ARGS> +EOS applyWrappedShiftAndScale(T&& eos, bool scaled, bool shifted, Real scale, Real shift, ARGS ... args) { + if (shifted && scaled) { + ShiftedEOS a(std::forward(eos), shift); + ScaledEOS> b(std::move(a), scale); + W>> c(std::move(b), args...); + return c; + } + if (shifted) { + ShiftedEOS sh_eos(std::forward(eos), shift); + return W>(std::move(sh_eos), args...); + } + if (scaled) { + ScaledEOS sc_eos(std::forward(eos), scale); + return W>(std::move(sc_eos), args...); + } + return W(std::forward(eos), args...); +} + +template +EOS applyShiftAndScaleAndSAPRamp(T&& eos, bool scaled, bool shifted, bool ramped, Real scale, Real shift, Real r0, Real a, Real b, Real c) { + if (ramped) { + return applyWrappedShiftAndScale(std::forward(eos), scaled, shifted, scale, shift, r0, a, b, c); + } else { + return applyShiftAndScale(std::forward(eos), scaled, shifted, scale, shift); + } +} } // namespace EOSBuilder } // namespace singularity diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index 3f8a7e2eab..ffde43d4ea 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -35,11 +35,13 @@ int init_sg_eos(const int nmat, EOS *&eos) { return 0; } -#define SGAPPLYMOD(A) \ - EOSBuilder::applyShiftAndScale(A, enabled[0] == 1, enabled[1] == 1, vals[0], vals[1]) +#define SGAPPLYMOD(A) \ + EOSBuilder::applyShiftAndScaleAndSAPRamp(A, enabled[0] == 1, enabled[1] == 1, \ + enabled[2] == 1, vals[0], vals[1], \ + vals[2], vals[3], vals[4], vals[5]) -constexpr const int def_en[2] = {0, 0}; -constexpr const double def_v[2] = {0.0, 0.0}; +constexpr const int def_en[3] = {0, 0, 0}; +constexpr const double def_v[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; int init_sg_IdealGas(const int matindex, EOS *eos, const double gm1, const double Cv, int const *const enabled, double const *const vals) { From 239c6d19eae56388aa22db9bde09dfa9e035434f Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Mon, 6 Jun 2022 12:22:42 -0600 Subject: [PATCH 08/34] hook ramp into eos builder --- singularity-eos/eos/eos_builder.cpp | 21 ++++++++++++++++----- singularity-eos/eos/eos_builder.hpp | 2 +- 2 files changed, 17 insertions(+), 6 deletions(-) diff --git a/singularity-eos/eos/eos_builder.cpp b/singularity-eos/eos/eos_builder.cpp index 74f7a765bc..818dd6e4bd 100644 --- a/singularity-eos/eos/eos_builder.cpp +++ b/singularity-eos/eos/eos_builder.cpp @@ -27,7 +27,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par bool shifted = (modifiers.count(EOSModifier::Shifted) > 0); bool relativistic = (modifiers.count(EOSModifier::Relativistic) > 0); bool units = (modifiers.count(EOSModifier::UnitSystem) > 0); - if ((shifted || scaled || relativistic || units) && !(isModifiable(type))) { + bool ramp = (modifiers.count(EOSModifier::SAPRamp) > 0); + if ((shifted || scaled || relativistic || units || ramp) && !(isModifiable(type))) { EOS_ERROR("Modifiers not supported for this EOS"); } if (relativistic && (shifted || scaled || units)) { @@ -53,6 +54,16 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par use_length_time = mpark::get(modifiers[EOSModifier::UnitSystem]["use_length_time"]); } + Real r0 = 0; + Real a = 1; + Real b = 0; + Real c = 0; + if (ramp) { + r0 = mpark::get(modifiers[EOSModifier::SAPRamp]["r0"]); + a = mpark::get(modifiers[EOSModifier::SAPRamp]["a"]); + b = mpark::get(modifiers[EOSModifier::SAPRamp]["b"]); + c = mpark::get(modifiers[EOSModifier::SAPRamp]["c"]); + } Real rho_unit = 1; Real sie_unit = 1; Real temp_unit = 1; @@ -93,7 +104,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(g), cl); } - return applyShiftAndScale(std::move(g), scaled, shifted, scale, shift); + return appyShiftAndScaleAndSapRamp(std::move(g), scaled, shifted, ramped, scale, shift, r0, a, b, c); } #ifdef SPINER_USE_HDF if (type == EOSType::SpinerEOSDependsRhoT || type == EOSType::SpinerEOSDependsRhoSie) { @@ -120,7 +131,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return applyShiftAndScale(std::move(s), scaled, shifted, scale, shift); + return appyShiftAndScaleAndSapRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); } } else { string materialName = mpark::get(base_params["materialName"]); @@ -143,7 +154,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return applyShiftAndScale(std::move(s), scaled, shifted, scale, shift); + return appyShiftAndScaleAndSapRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); } } } @@ -165,7 +176,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return applyShiftAndScale(std::move(s), scaled, shifted, scale, shift); + return appyShiftAndScaleAndSapRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); } #endif if (type == EOSType::Gruneisen) { diff --git a/singularity-eos/eos/eos_builder.hpp b/singularity-eos/eos/eos_builder.hpp index 9eecd5372e..9c6e787e90 100644 --- a/singularity-eos/eos/eos_builder.hpp +++ b/singularity-eos/eos/eos_builder.hpp @@ -43,7 +43,7 @@ enum class EOSType { StellarCollapse #endif }; -enum class EOSModifier { Scaled, Shifted, Relativistic, UnitSystem }; +enum class EOSModifier { Scaled, Shifted, Relativistic, UnitSystem, SAPRamp }; // evil type erasure using param_t = mpark::variant; From cea0db1c4c6eb19b3547c4ed70cf858bfaf8779a Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Mon, 6 Jun 2022 15:58:54 -0600 Subject: [PATCH 09/34] builder for ramp works --- singularity-eos/eos/eos_builder.cpp | 14 +++++++------- singularity-eos/eos/eos_builder.hpp | 1 + test/eos_unit_tests.cpp | 15 +++++++++++++++ 3 files changed, 23 insertions(+), 7 deletions(-) diff --git a/singularity-eos/eos/eos_builder.cpp b/singularity-eos/eos/eos_builder.cpp index 818dd6e4bd..7af694a728 100644 --- a/singularity-eos/eos/eos_builder.cpp +++ b/singularity-eos/eos/eos_builder.cpp @@ -27,8 +27,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par bool shifted = (modifiers.count(EOSModifier::Shifted) > 0); bool relativistic = (modifiers.count(EOSModifier::Relativistic) > 0); bool units = (modifiers.count(EOSModifier::UnitSystem) > 0); - bool ramp = (modifiers.count(EOSModifier::SAPRamp) > 0); - if ((shifted || scaled || relativistic || units || ramp) && !(isModifiable(type))) { + bool ramped = (modifiers.count(EOSModifier::SAPRamp) > 0); + if ((shifted || scaled || relativistic || units || ramped) && !(isModifiable(type))) { EOS_ERROR("Modifiers not supported for this EOS"); } if (relativistic && (shifted || scaled || units)) { @@ -58,7 +58,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par Real a = 1; Real b = 0; Real c = 0; - if (ramp) { + if (ramped) { r0 = mpark::get(modifiers[EOSModifier::SAPRamp]["r0"]); a = mpark::get(modifiers[EOSModifier::SAPRamp]["a"]); b = mpark::get(modifiers[EOSModifier::SAPRamp]["b"]); @@ -104,7 +104,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(g), cl); } - return appyShiftAndScaleAndSapRamp(std::move(g), scaled, shifted, ramped, scale, shift, r0, a, b, c); + return applyShiftAndScaleAndSAPRamp(std::move(g), scaled, shifted, ramped, scale, shift, r0, a, b, c); } #ifdef SPINER_USE_HDF if (type == EOSType::SpinerEOSDependsRhoT || type == EOSType::SpinerEOSDependsRhoSie) { @@ -131,7 +131,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return appyShiftAndScaleAndSapRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); + return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); } } else { string materialName = mpark::get(base_params["materialName"]); @@ -154,7 +154,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return appyShiftAndScaleAndSapRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); + return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); } } } @@ -176,7 +176,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return appyShiftAndScaleAndSapRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); + return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); } #endif if (type == EOSType::Gruneisen) { diff --git a/singularity-eos/eos/eos_builder.hpp b/singularity-eos/eos/eos_builder.hpp index 9c6e787e90..9d19927089 100644 --- a/singularity-eos/eos/eos_builder.hpp +++ b/singularity-eos/eos/eos_builder.hpp @@ -51,6 +51,7 @@ using params_t = std::map; using modifiers_t = std::map; const params_t NO_PARAMS; +// TODO: Extend if needed const std::unordered_set modifiable({EOSType::IdealGas #ifdef SPINER_USE_HDF , diff --git a/test/eos_unit_tests.cpp b/test/eos_unit_tests.cpp index c76cce5f85..62e827289e 100644 --- a/test/eos_unit_tests.cpp +++ b/test/eos_unit_tests.cpp @@ -266,6 +266,21 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { THEN("The shift and scale parameters pass through correctly") { REQUIRE(eos.PressureFromDensityInternalEnergy(rho, sie) == 0.3); } + WHEN("We add a ramp") { + EOSBuilder::params_t ramp_params; + Real r0 = 1; + Real a = 1; + Real b = 0; + Real c = 0; + ramp_params["r0"] = r0; + ramp_params["a"] = a; + ramp_params["b"] = b; + ramp_params["c"] = c; + modifiers[EOSBuilder::EOSModifier::SAPRamp] = ramp_params; + THEN("The EOS is constructed correctly") { + auto eos_ramped = EOSBuilder::buildEOS(type, base_params, modifiers); + } + } } WHEN("We construct a non-modifying modifier") { EOS ig = IdealGas(gm1, Cv); From 63a65959e9fbd2cdaf641a26ea78e87986c7d3df Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Mon, 6 Jun 2022 16:39:39 -0600 Subject: [PATCH 10/34] Adding palpha2ramp function and a possible way of initializing it. Adding more comprehensive ramp tests to unit tests. --- singularity-eos/eos/modifiers/ramps_eos.hpp | 37 +++++++++++++++++++++ singularity-eos/eos/singularity_eos.cpp | 21 +++++++++--- test/eos_unit_tests.cpp | 22 +++++++++++- 3 files changed, 75 insertions(+), 5 deletions(-) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index ad2c587718..5c0d95b82a 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -32,6 +32,43 @@ namespace singularity { using namespace eos_base; +template +void pAlpha2SAPRampParams(const T& eos, + const Real alpha0, const Real Pe, const Real Pc, + Real& r0, Real& a, Real& b, Real& c) { + // get reference conditions + Real rho0, T0, sie0, P0, cv0, bmod0, dpde0, dvdt0; + Real rmid, r1; + eos.ValuesAtReferenceState(rho0, T0, sie0, P0, cv0, bmod0, dpde0, dvdt0); + // calculate r0 + r0 = rho0 / alpha0; + // calculate rmid + auto rmid_func = [&] (const Real x) { + return eos.PressureFromDensityTemperature(alpha0*x, T0); + }; + RootFinding1D::RootCounts co{}; + // get upper bound to density informed by the reference + // bulk modulus + const Real max_exp_arg = std::log(std::numeric_limits::max()*0.99); + const Real exp_arg = std::min(max_exp_arg, (2.0*Pc - P0) / bmod0); + const Real rho_ub = rho0*std::exp(exp_arg); + // finds where rmid_func = Pe + RootFinding1D::findRoot(rmid_func, Pe, rho0, r0, rho_ub, 1.e-12, 1.e-12, rmid, co); + // calculate r1 + auto r1_func = [&] (const Real x) { + return eos.PressureFromDensityTemperature(x, T0); + }; + // finds where r1_func = Pc + RootFinding1D::findRoot(r1_func, Pc, rmid, r0, rho_ub, 1.e-12, 1.e-12, r1, co); + // a + a = r0 * Pe / (rmid - r0); + // b + b = r0 * (Pc - Pe) / (r1 - rmid); + // c + c = (Pc * rmid - Pe * r1) / (r0 * (Pc - Pe) ); + return; +} + template class SAPRampEOS : public EosBase> { public: diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index b76922573d..3ceb6fe0a7 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -35,17 +35,30 @@ int init_sg_eos(const int nmat, EOS *&eos) { return 0; } +// apply everything but ramp in order to possibly calculate the +// SAP ramp parameters from p-alhpa ramp parameters +#define SGAPPLYMODSIMPLE(A) \ + EOSBuilder::applyShiftAndScale(A, enabled[0] == 1, enabled[1] == 1, \ + vals[0], vals[1]) + #define SGAPPLYMOD(A) \ EOSBuilder::applyShiftAndScaleAndSAPRamp(A, enabled[0] == 1, enabled[1] == 1, \ - enabled[2] == 1, vals[0], vals[1], \ - vals[2], vals[3], vals[4], vals[5]) + enabled[2] == 1 || enabled[3] == 1, \ + vals[0], vals[1], vals[2], vals[3], \ + vals[4], vals[5]) -constexpr const int def_en[3] = {0, 0, 0}; +constexpr const int def_en[4] = {0, 0, 0, 0}; constexpr const double def_v[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; int init_sg_IdealGas(const int matindex, EOS *eos, const double gm1, const double Cv, - int const *const enabled, double const *const vals) { + int const *const enabled, double const *const vals_i) { assert(matindex >= 0); + Real vals[6] = {vals_i[0], vals_i[1], vals_i[2], vals_i[3], vals_i[4], vals_i[5]}; + EOS eosi = SGAPPLYMODSIMPLE(IdealGas(gm1, Cv)); + if (enabled[3] == 1) { + singularity::pAlpha2SAPRampParams(eosi, vals_i[2], vals_i[3], vals_i[4], vals[2], + vals[3], vals[4], vals[5]); + } EOS eos_ = SGAPPLYMOD(IdealGas(gm1, Cv)); eos[matindex] = eos_.GetOnDevice(); return 0; diff --git a/test/eos_unit_tests.cpp b/test/eos_unit_tests.cpp index c76cce5f85..df9ca0ee88 100644 --- a/test/eos_unit_tests.cpp +++ b/test/eos_unit_tests.cpp @@ -26,6 +26,7 @@ #include #include +#include // typename demangler #ifdef __GNUG__ @@ -271,7 +272,26 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { EOS ig = IdealGas(gm1, Cv); EOS igsh = ScaledEOS(IdealGas(gm1, Cv), 1.0); EOS igsc = ShiftedEOS(IdealGas(gm1, Cv), 0.0); - auto igra = SAPRampEOS(IdealGas(gm1, Cv), 1.e9, 1.0, 2.0, 1.0); + EOS igra;// = SAPRampEOS(IdealGas(gm1, Cv), 1.e9, 1.0, 2.0, 1.0); + // test out the c interface + int enabled[4] = {0, 0, 1, 0}; + Real vals[6] = {0.0, 0.0, 1.e9, 1.0, 2.0, 1.0}; + Real rho0 = 1.e6 / (gm1*Cv*293.0); + init_sg_IdealGas(0, &igra, gm1, Cv, enabled, vals); + EOS igra2; + enabled[2] = 0; + enabled[3] = 1; + Real Pe = 2.e6, Pc = 1.e7; + Real alpha0 = 0.98; + vals[2] = alpha0; + vals[3] = Pe; + vals[4] = Pc; + vals[5] = 0.0; + Real r1 = Pc / (gm1*Cv*293.0); + Real rmid = Pe / (gm1*Cv*293.0*alpha0); + init_sg_IdealGas(0, &igra2, gm1, Cv, enabled, vals); + REQUIRE(isClose(Pe, igra2.PressureFromDensityTemperature(alpha0*rmid, 293.0), 1.e-12)); + REQUIRE(isClose(Pc, igra2.PressureFromDensityTemperature(r1, 293.0), 1.e-12)); THEN("The modified EOS should produce equivalent results") { compare_two_eoss(igsh, ig); compare_two_eoss(igsc, ig); From 388333bdced97f30090e99bf9de4f055a9b3273a Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Wed, 8 Jun 2022 19:51:23 -0600 Subject: [PATCH 11/34] Adding alternate ramp bulk modulus calculation. --- singularity-eos/eos/modifiers/ramps_eos.hpp | 21 +++++++++++++++++++-- 1 file changed, 19 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 5c0d95b82a..6f6626ad4e 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -95,6 +95,14 @@ class SAPRampEOS : public EosBase> { return p_ramp; } + PORTABLE_INLINE_FUNCTION + Real get_ramp_dpdrho(Real rho) const { + const Real dpdr{rho < r0_ ? 0.0 + : rho < rmid_ ? a_ / r0_ + : b_ / r0_}; + return dpdr; + } + PORTABLE_FUNCTION Real TemperatureFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda = nullptr) const { @@ -123,7 +131,12 @@ class SAPRampEOS : public EosBase> { PORTABLE_FUNCTION Real BulkModulusFromDensityInternalEnergy(const Real rho, const Real sie, Real *lambda = nullptr) const { - return t_.BulkModulusFromDensityInternalEnergy(rho, sie, lambda); + const Real p_ramp{get_ramp_pressure(rho)}; + const Real p_eos{t_.PressureFromDensityInternalEnergy(rho, sie, lambda)}; + + return + p_eos < p_ramp ? rho*get_ramp_dpdrho(rho) + : t_.BulkModulusFromDensityInternalEnergy(rho, sie, lambda); } PORTABLE_FUNCTION Real GruneisenParamFromDensityInternalEnergy(const Real rho, const Real sie, @@ -145,7 +158,11 @@ class SAPRampEOS : public EosBase> { PORTABLE_FUNCTION Real BulkModulusFromDensityTemperature(const Real rho, const Real temperature, Real *lambda = nullptr) const { - return t_.BulkModulusFromDensityTemperature(rho, temperature, lambda); + const Real p_ramp{get_ramp_pressure(rho)}; + const Real p_eos{t_.PressureFromDensityTemperature(rho, temperature, lambda)}; + return + p_eos < p_ramp ? rho*get_ramp_dpdrho(rho) + : t_.BulkModulusFromDensityTemperature(rho, temperature, lambda); } PORTABLE_FUNCTION Real GruneisenParamFromDensityTemperature(const Real rho, const Real temperature, From 895c187ca9719c1d7250f48066db4210b9018ee4 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 14 Jun 2022 09:38:32 -0600 Subject: [PATCH 12/34] Making the values array nonconst such that the init functions can convert p-alpha inputs to ramp inputs. --- singularity-eos/eos/singularity_eos.cpp | 58 ++++++++++++++++++++----- singularity-eos/eos/singularity_eos.f90 | 32 +++++++------- singularity-eos/eos/singularity_eos.hpp | 16 +++---- 3 files changed, 70 insertions(+), 36 deletions(-) diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index 3ceb6fe0a7..e053babdcf 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -47,16 +47,15 @@ int init_sg_eos(const int nmat, EOS *&eos) { vals[0], vals[1], vals[2], vals[3], \ vals[4], vals[5]) -constexpr const int def_en[4] = {0, 0, 0, 0}; -constexpr const double def_v[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; +int def_en[4] = {0, 0, 0, 0}; +double def_v[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; int init_sg_IdealGas(const int matindex, EOS *eos, const double gm1, const double Cv, - int const *const enabled, double const *const vals_i) { + int const *const enabled, double *const vals) { assert(matindex >= 0); - Real vals[6] = {vals_i[0], vals_i[1], vals_i[2], vals_i[3], vals_i[4], vals_i[5]}; EOS eosi = SGAPPLYMODSIMPLE(IdealGas(gm1, Cv)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals_i[2], vals_i[3], vals_i[4], vals[2], + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(IdealGas(gm1, Cv)); @@ -72,8 +71,13 @@ int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const doubl const double s2, const double s3, const double G0, const double b, const double rho0, const double T0, const double P0, const double Cv, int const *const enabled, - double const *const vals) { + double *const vals) { assert(matindex >= 0); + EOS eosi = SGAPPLYMODSIMPLE(Gruneisen(C0, s1, s2, s3, G0, b, rho0, T0, P0, Cv)); + if (enabled[3] == 1) { + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); + } EOS eos_ = SGAPPLYMOD(Gruneisen(C0, s1, s2, s3, G0, b, rho0, T0, P0, Cv)); eos[matindex] = eos_.GetOnDevice(); return 0; @@ -89,8 +93,13 @@ int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const doubl int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, const double R1, const double R2, const double w, const double rho0, - const double Cv, int const *const enabled, double const *const vals) { + const double Cv, int const *const enabled, double *const vals) { assert(matindex >= 0); + EOS eosi = SGAPPLYMODSIMPLE(JWL(A, B, R1, R2, w, rho0, Cv)); + if (enabled[3] == 1) { + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); + } EOS eos_ = SGAPPLYMOD(JWL(A, B, R1, R2, w, rho0, Cv)); eos[matindex] = eos_.GetOnDevice(); return 0; @@ -105,8 +114,13 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, const double pc, const double Cv, const double E0, - int const *const enabled, double const *const vals) { + int const *const enabled, double *const vals) { assert(matindex >= 0); + EOS eosi = SGAPPLYMODSIMPLE(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); + if (enabled[3] == 1) { + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); + } EOS eos_ = SGAPPLYMOD(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); eos[matindex] = eos_.GetOnDevice(); return 0; @@ -123,8 +137,13 @@ int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, const double A, const double B, const double C, const double G0, const double Z, const double alpha, const double Cv0, int const *const enabled, - double const *const vals) { + double *const vals) { assert(matindex >= 0); + EOS eosi = SGAPPLYMODSIMPLE(DavisReactants(rho0, e0, P0, T0, A, B, C, G0, Z, alpha, Cv0)); + if (enabled[3] == 1) { + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); + } EOS eos_ = SGAPPLYMOD(DavisReactants(rho0, e0, P0, T0, A, B, C, G0, Z, alpha, Cv0)); eos[matindex] = eos_.GetOnDevice(); return 0; @@ -142,8 +161,13 @@ int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, #ifdef SPINER_USE_HDF int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename, const int matid, int const *const enabled, - double const *const vals) { + double *const vals) { assert(matindex >= 0); + EOS eosi = SGAPPLYMODSIMPLE(SpinerEOSDependsRhoT(std::string(filename), matid)); + if (enabled[3] == 1) { + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); + } EOS eos_ = SGAPPLYMOD(SpinerEOSDependsRhoT(std::string(filename), matid)); eos[matindex] = eos_.GetOnDevice(); return 0; @@ -156,8 +180,13 @@ int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename int init_sg_SpinerDependsRhoSie(const int matindex, EOS *eos, const char *filename, const int matid, int const *const enabled, - double const *const vals) { + double *const vals) { assert(matindex >= 0); + EOS eosi = SGAPPLYMODSIMPLE(SpinerEOSDependsRhoSie(std::string(filename), matid)); + if (enabled[3] == 1) { + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); + } EOS eos_ = SGAPPLYMOD(SpinerEOSDependsRhoSie(std::string(filename), matid)); eos[matindex] = eos_.GetOnDevice(); return 0; @@ -170,8 +199,13 @@ int init_sg_SpinerDependsRhoSie(const int matindex, EOS *eos, const char *filena #ifdef SINGULARITY_USE_EOSPAC int init_sg_eospac(const int matindex, EOS *eos, const int id, int const *const enabled, - double const *const vals) { + double *const vals) { assert(matindex >= 0); + EOS eosi = SGAPPLYMODSIMPLE(EOSPAC(id)); + if (enabled[3] == 1) { + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); + } EOS eos_ = SGAPPLYMOD(EOSPAC(id)); eos[matindex] = eos_.GetOnDevice(); return 0; diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index 42691455d6..5f539099c6 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -282,8 +282,8 @@ integer function init_sg_IdealGas_f(matindex, eos, gm1, Cv, & integer(c_int), value, intent(in) :: matindex type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: gm1, Cv - integer(kind=c_int), dimension(:), target, intent(in) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(in) :: sg_mods_values + integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values err = init_sg_IdealGas(matindex-1, eos%ptr, gm1, Cv, & c_loc(sg_mods_enabled), c_loc(sg_mods_values)) end function init_sg_IdealGas_f @@ -296,8 +296,8 @@ integer function init_sg_Gruneisen_f(matindex, eos, C0, s1, s2, s3, G0, b,& type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: C0, s1, s2, s3, G0, b, rho0 real(kind=8), value, intent(in) :: T0, P0, Cv - integer(kind=c_int), dimension(:), target, intent(in) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(in) :: sg_mods_values + integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values err = init_sg_Gruneisen(matindex-1, eos%ptr, C0, s1, s2, s3, G0, b, rho0,& T0, P0, Cv, c_loc(sg_mods_enabled), & c_loc(sg_mods_values)) @@ -309,8 +309,8 @@ integer function init_sg_JWL_f(matindex, eos, A, B, R1, R2, w, rho0, Cv, & integer(c_int), value, intent(in) :: matindex type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: A, B, R1, R2, w, rho0, Cv - integer(kind=c_int), dimension(:), target, intent(in) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(in) :: sg_mods_values + integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values err = init_sg_JWL(matindex-1, eos%ptr, A, B, R1, R2, w, rho0, Cv, & c_loc(sg_mods_enabled), c_loc(sg_mods_values)) end function init_sg_JWL_f @@ -322,8 +322,8 @@ integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, & integer(c_int), value, intent(in) :: matindex type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: a, b, k, n, vc, pc, Cv, E0 - integer(kind=c_int), dimension(:), target, intent(in) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(in) :: sg_mods_values + integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values err = init_sg_DavisProducts(matindex-1, eos%ptr, a, b, k, n, vc, pc, Cv, & E0, c_loc(sg_mods_enabled), & c_loc(sg_mods_values)) @@ -337,8 +337,8 @@ integer function init_sg_DavisReactants_f(matindex, eos, rho0, e0, P0, T0, & type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: rho0, e0, P0, T0, A, B, C, G0, Z real(kind=8), value, intent(in) :: alpha, Cv0 - integer(kind=c_int), dimension(:), target, intent(in) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(in) :: sg_mods_values + integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values err = init_sg_DavisReactants(matindex-1, eos%ptr, rho0, e0, P0, T0, A, B, & C, G0, Z, alpha, Cv0, c_loc(sg_mods_enabled), & c_loc(sg_mods_values)) @@ -352,8 +352,8 @@ integer function init_sg_SpinerDependsRhoT_f(matindex, eos, filename, id, & type(sg_eos_ary_t), intent(in) :: eos character(len=*, kind=c_char), intent(in) :: filename integer(c_int), intent(inout) :: id - integer(kind=c_int), dimension(:), target, intent(in) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(in) :: sg_mods_values + integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values err = init_sg_SpinerDependsRhoT(matindex-1, eos%ptr,& trim(filename)//C_NULL_CHAR, id, & c_loc(sg_mods_enabled), & @@ -367,8 +367,8 @@ integer function init_sg_SpinerDependsRhoSie_f(matindex, eos, filename, id, & integer(c_int), value, intent(in) :: matindex, id type(sg_eos_ary_t), intent(in) :: eos character(len=*, kind=c_char), intent(in) :: filename - integer(kind=c_int), dimension(:), target, intent(in) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(in) :: sg_mods_values + integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values err = init_sg_SpinerDependsRhoSie(matindex-1, eos%ptr,& trim(filename)//C_NULL_CHAR, id, & c_loc(sg_mods_enabled), & @@ -380,8 +380,8 @@ integer function init_sg_eospac_f(matindex, eos, id, sg_mods_enabled, & result(err) integer(c_int), value, intent(in) :: matindex, id type(sg_eos_ary_t), intent(in) :: eos - integer(kind=c_int), dimension(:), target, intent(in) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(in) :: sg_mods_values + integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values err = init_sg_eospac(matindex-1, eos%ptr, id, c_loc(sg_mods_enabled), & c_loc(sg_mods_values)) end function init_sg_eospac_f diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index 02cb0ba170..0c8b4ba75e 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/singularity-eos/eos/singularity_eos.hpp @@ -27,44 +27,44 @@ extern "C" { int init_sg_eos(const int nmat, EOS *&eos); int init_sg_IdealGas(const int matindex, EOS *eos, const double gm1, const double Cv, - int const *const enabled, double const *const vals); + int const *const enabled, double *const vals); int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, const double R1, const double R2, const double w, const double rho0, - const double Cv, int const *const enabled, double const *const vals); + const double Cv, int const *const enabled, double *const vals); int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const double s1, const double s2, const double s3, const double G0, const double b, const double rho0, const double T0, const double P0, const double Cv, int const *const enabled, - double const *const vals); + double *const vals); int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, const double pc, const double Cv, const double E0, - int const *const enabled, double const *const vals); + int const *const enabled, double *const vals); int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, const double e0, const double P0, const double T0, const double A, const double B, const double C, const double G0, const double Z, const double alpha, const double Cv0, int const *const enabled, - double const *const vals); + double *const vals); #ifdef SPINER_USE_HDF int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename, const int id, int const *const enabled, - double const *const vals); + double *const vals); int init_sg_SpinerDependsRhoSie(const int matindex, EOS *eos, const char *filename, const int id, int const *const enabled, - double const *const vals); + double *const vals); #endif #ifdef SINGULARITY_USE_EOSPAC // capitalize? eospaceos Eospac Eospaceos EOSPAC EOSPACeos? int init_sg_eospac(const int matindex, EOS *eos, const int id, int const *const enabled, - double const *const vals); + double *const vals); #endif // SINGULARITY_USE_EOSPAC int get_sg_eos( // sizing information From c2c2b62351e6cdb425feef0d8b2d3fca61147684 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 14 Jun 2022 12:50:39 -0600 Subject: [PATCH 13/34] Formatting. --- singularity-eos/eos/eos_builder.cpp | 12 +++-- singularity-eos/eos/eos_builder.hpp | 12 +++-- singularity-eos/eos/modifiers/ramps_eos.hpp | 40 +++++++-------- singularity-eos/eos/singularity_eos.cpp | 54 ++++++++++----------- singularity-eos/eos/singularity_eos.hpp | 6 +-- test/test_eos_unit.cpp | 11 +++-- 6 files changed, 67 insertions(+), 68 deletions(-) diff --git a/singularity-eos/eos/eos_builder.cpp b/singularity-eos/eos/eos_builder.cpp index 7af694a728..f73fe5b28b 100644 --- a/singularity-eos/eos/eos_builder.cpp +++ b/singularity-eos/eos/eos_builder.cpp @@ -104,7 +104,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(g), cl); } - return applyShiftAndScaleAndSAPRamp(std::move(g), scaled, shifted, ramped, scale, shift, r0, a, b, c); + return applyShiftAndScaleAndSAPRamp(std::move(g), scaled, shifted, ramped, scale, + shift, r0, a, b, c); } #ifdef SPINER_USE_HDF if (type == EOSType::SpinerEOSDependsRhoT || type == EOSType::SpinerEOSDependsRhoSie) { @@ -131,7 +132,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); + return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, + shift, r0, a, b, c); } } else { string materialName = mpark::get(base_params["materialName"]); @@ -154,7 +156,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); + return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, + shift, r0, a, b, c); } } } @@ -176,7 +179,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); + return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, + shift, r0, a, b, c); } #endif if (type == EOSType::Gruneisen) { diff --git a/singularity-eos/eos/eos_builder.hpp b/singularity-eos/eos/eos_builder.hpp index 9d19927089..20c598d578 100644 --- a/singularity-eos/eos/eos_builder.hpp +++ b/singularity-eos/eos/eos_builder.hpp @@ -105,8 +105,9 @@ EOS applyShiftAndScale(T &&eos, bool scaled, bool shifted, Real scale, Real shif return eos; } -template class W, typename ... ARGS> -EOS applyWrappedShiftAndScale(T&& eos, bool scaled, bool shifted, Real scale, Real shift, ARGS ... args) { +template class W, typename... ARGS> +EOS applyWrappedShiftAndScale(T &&eos, bool scaled, bool shifted, Real scale, Real shift, + ARGS... args) { if (shifted && scaled) { ShiftedEOS a(std::forward(eos), shift); ScaledEOS> b(std::move(a), scale); @@ -125,9 +126,12 @@ EOS applyWrappedShiftAndScale(T&& eos, bool scaled, bool shifted, Real scale, Re } template -EOS applyShiftAndScaleAndSAPRamp(T&& eos, bool scaled, bool shifted, bool ramped, Real scale, Real shift, Real r0, Real a, Real b, Real c) { +EOS applyShiftAndScaleAndSAPRamp(T &&eos, bool scaled, bool shifted, bool ramped, + Real scale, Real shift, Real r0, Real a, Real b, + Real c) { if (ramped) { - return applyWrappedShiftAndScale(std::forward(eos), scaled, shifted, scale, shift, r0, a, b, c); + return applyWrappedShiftAndScale(std::forward(eos), scaled, shifted, + scale, shift, r0, a, b, c); } else { return applyShiftAndScale(std::forward(eos), scaled, shifted, scale, shift); } diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 6f6626ad4e..81bb637119 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -32,10 +32,9 @@ namespace singularity { using namespace eos_base; -template -void pAlpha2SAPRampParams(const T& eos, - const Real alpha0, const Real Pe, const Real Pc, - Real& r0, Real& a, Real& b, Real& c) { +template +void pAlpha2SAPRampParams(const T &eos, const Real alpha0, const Real Pe, const Real Pc, + Real &r0, Real &a, Real &b, Real &c) { // get reference conditions Real rho0, T0, sie0, P0, cv0, bmod0, dpde0, dvdt0; Real rmid, r1; @@ -43,21 +42,19 @@ void pAlpha2SAPRampParams(const T& eos, // calculate r0 r0 = rho0 / alpha0; // calculate rmid - auto rmid_func = [&] (const Real x) { - return eos.PressureFromDensityTemperature(alpha0*x, T0); + auto rmid_func = [&](const Real x) { + return eos.PressureFromDensityTemperature(alpha0 * x, T0); }; RootFinding1D::RootCounts co{}; // get upper bound to density informed by the reference // bulk modulus - const Real max_exp_arg = std::log(std::numeric_limits::max()*0.99); - const Real exp_arg = std::min(max_exp_arg, (2.0*Pc - P0) / bmod0); - const Real rho_ub = rho0*std::exp(exp_arg); + const Real max_exp_arg = std::log(std::numeric_limits::max() * 0.99); + const Real exp_arg = std::min(max_exp_arg, (2.0 * Pc - P0) / bmod0); + const Real rho_ub = rho0 * std::exp(exp_arg); // finds where rmid_func = Pe RootFinding1D::findRoot(rmid_func, Pe, rho0, r0, rho_ub, 1.e-12, 1.e-12, rmid, co); // calculate r1 - auto r1_func = [&] (const Real x) { - return eos.PressureFromDensityTemperature(x, T0); - }; + auto r1_func = [&](const Real x) { return eos.PressureFromDensityTemperature(x, T0); }; // finds where r1_func = Pc RootFinding1D::findRoot(r1_func, Pc, rmid, r0, rho_ub, 1.e-12, 1.e-12, r1, co); // a @@ -65,7 +62,7 @@ void pAlpha2SAPRampParams(const T& eos, // b b = r0 * (Pc - Pe) / (r1 - rmid); // c - c = (Pc * rmid - Pe * r1) / (r0 * (Pc - Pe) ); + c = (Pc * rmid - Pe * r1) / (r0 * (Pc - Pe)); return; } @@ -97,9 +94,7 @@ class SAPRampEOS : public EosBase> { PORTABLE_INLINE_FUNCTION Real get_ramp_dpdrho(Real rho) const { - const Real dpdr{rho < r0_ ? 0.0 - : rho < rmid_ ? a_ / r0_ - : b_ / r0_}; + const Real dpdr{rho < r0_ ? 0.0 : rho < rmid_ ? a_ / r0_ : b_ / r0_}; return dpdr; } @@ -133,10 +128,9 @@ class SAPRampEOS : public EosBase> { Real *lambda = nullptr) const { const Real p_ramp{get_ramp_pressure(rho)}; const Real p_eos{t_.PressureFromDensityInternalEnergy(rho, sie, lambda)}; - - return - p_eos < p_ramp ? rho*get_ramp_dpdrho(rho) - : t_.BulkModulusFromDensityInternalEnergy(rho, sie, lambda); + + return p_eos < p_ramp ? rho * get_ramp_dpdrho(rho) + : t_.BulkModulusFromDensityInternalEnergy(rho, sie, lambda); } PORTABLE_FUNCTION Real GruneisenParamFromDensityInternalEnergy(const Real rho, const Real sie, @@ -160,9 +154,9 @@ class SAPRampEOS : public EosBase> { Real *lambda = nullptr) const { const Real p_ramp{get_ramp_pressure(rho)}; const Real p_eos{t_.PressureFromDensityTemperature(rho, temperature, lambda)}; - return - p_eos < p_ramp ? rho*get_ramp_dpdrho(rho) - : t_.BulkModulusFromDensityTemperature(rho, temperature, lambda); + return p_eos < p_ramp + ? rho * get_ramp_dpdrho(rho) + : t_.BulkModulusFromDensityTemperature(rho, temperature, lambda); } PORTABLE_FUNCTION Real GruneisenParamFromDensityTemperature(const Real rho, const Real temperature, diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index e053babdcf..cecac983a0 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -35,17 +35,15 @@ int init_sg_eos(const int nmat, EOS *&eos) { return 0; } -// apply everything but ramp in order to possibly calculate the +// apply everything but ramp in order to possibly calculate the // SAP ramp parameters from p-alhpa ramp parameters -#define SGAPPLYMODSIMPLE(A) \ - EOSBuilder::applyShiftAndScale(A, enabled[0] == 1, enabled[1] == 1, \ - vals[0], vals[1]) +#define SGAPPLYMODSIMPLE(A) \ + EOSBuilder::applyShiftAndScale(A, enabled[0] == 1, enabled[1] == 1, vals[0], vals[1]) -#define SGAPPLYMOD(A) \ - EOSBuilder::applyShiftAndScaleAndSAPRamp(A, enabled[0] == 1, enabled[1] == 1, \ - enabled[2] == 1 || enabled[3] == 1, \ - vals[0], vals[1], vals[2], vals[3], \ - vals[4], vals[5]) +#define SGAPPLYMOD(A) \ + EOSBuilder::applyShiftAndScaleAndSAPRamp(A, enabled[0] == 1, enabled[1] == 1, \ + enabled[2] == 1 || enabled[3] == 1, vals[0], \ + vals[1], vals[2], vals[3], vals[4], vals[5]) int def_en[4] = {0, 0, 0, 0}; double def_v[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; @@ -55,8 +53,8 @@ int init_sg_IdealGas(const int matindex, EOS *eos, const double gm1, const doubl assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(IdealGas(gm1, Cv)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], - vals[3], vals[4], vals[5]); + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(IdealGas(gm1, Cv)); eos[matindex] = eos_.GetOnDevice(); @@ -70,13 +68,12 @@ int init_sg_IdealGas(const int matindex, EOS *eos, const double gm1, const doubl int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const double s1, const double s2, const double s3, const double G0, const double b, const double rho0, const double T0, const double P0, - const double Cv, int const *const enabled, - double *const vals) { + const double Cv, int const *const enabled, double *const vals) { assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(Gruneisen(C0, s1, s2, s3, G0, b, rho0, T0, P0, Cv)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], - vals[3], vals[4], vals[5]); + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(Gruneisen(C0, s1, s2, s3, G0, b, rho0, T0, P0, Cv)); eos[matindex] = eos_.GetOnDevice(); @@ -97,8 +94,8 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(JWL(A, B, R1, R2, w, rho0, Cv)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], - vals[3], vals[4], vals[5]); + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(JWL(A, B, R1, R2, w, rho0, Cv)); eos[matindex] = eos_.GetOnDevice(); @@ -118,8 +115,8 @@ int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const do assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], - vals[3], vals[4], vals[5]); + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); eos[matindex] = eos_.GetOnDevice(); @@ -139,10 +136,11 @@ int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, const double Cv0, int const *const enabled, double *const vals) { assert(matindex >= 0); - EOS eosi = SGAPPLYMODSIMPLE(DavisReactants(rho0, e0, P0, T0, A, B, C, G0, Z, alpha, Cv0)); + EOS eosi = + SGAPPLYMODSIMPLE(DavisReactants(rho0, e0, P0, T0, A, B, C, G0, Z, alpha, Cv0)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], - vals[3], vals[4], vals[5]); + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(DavisReactants(rho0, e0, P0, T0, A, B, C, G0, Z, alpha, Cv0)); eos[matindex] = eos_.GetOnDevice(); @@ -165,8 +163,8 @@ int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(SpinerEOSDependsRhoT(std::string(filename), matid)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], - vals[3], vals[4], vals[5]); + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(SpinerEOSDependsRhoT(std::string(filename), matid)); eos[matindex] = eos_.GetOnDevice(); @@ -184,8 +182,8 @@ int init_sg_SpinerDependsRhoSie(const int matindex, EOS *eos, const char *filena assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(SpinerEOSDependsRhoSie(std::string(filename), matid)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], - vals[3], vals[4], vals[5]); + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(SpinerEOSDependsRhoSie(std::string(filename), matid)); eos[matindex] = eos_.GetOnDevice(); @@ -203,8 +201,8 @@ int init_sg_eospac(const int matindex, EOS *eos, const int id, int const *const assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(EOSPAC(id)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], - vals[3], vals[4], vals[5]); + singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(EOSPAC(id)); eos[matindex] = eos_.GetOnDevice(); diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index 0c8b4ba75e..5ea962c1c7 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/singularity-eos/eos/singularity_eos.hpp @@ -36,8 +36,7 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const double s1, const double s2, const double s3, const double G0, const double b, const double rho0, const double T0, const double P0, - const double Cv, int const *const enabled, - double *const vals); + const double Cv, int const *const enabled, double *const vals); int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, @@ -53,8 +52,7 @@ int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, #ifdef SPINER_USE_HDF int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename, - const int id, int const *const enabled, - double *const vals); + const int id, int const *const enabled, double *const vals); int init_sg_SpinerDependsRhoSie(const int matindex, EOS *eos, const char *filename, const int id, int const *const enabled, diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 14d6606516..265f61d569 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -248,11 +248,11 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { EOS ig = IdealGas(gm1, Cv); EOS igsh = ScaledEOS(IdealGas(gm1, Cv), 1.0); EOS igsc = ShiftedEOS(IdealGas(gm1, Cv), 0.0); - EOS igra;// = SAPRampEOS(IdealGas(gm1, Cv), 1.e9, 1.0, 2.0, 1.0); + EOS igra; // = SAPRampEOS(IdealGas(gm1, Cv), 1.e9, 1.0, 2.0, 1.0); // test out the c interface int enabled[4] = {0, 0, 1, 0}; Real vals[6] = {0.0, 0.0, 1.e9, 1.0, 2.0, 1.0}; - Real rho0 = 1.e6 / (gm1*Cv*293.0); + Real rho0 = 1.e6 / (gm1 * Cv * 293.0); init_sg_IdealGas(0, &igra, gm1, Cv, enabled, vals); EOS igra2; enabled[2] = 0; @@ -263,10 +263,11 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { vals[3] = Pe; vals[4] = Pc; vals[5] = 0.0; - Real r1 = Pc / (gm1*Cv*293.0); - Real rmid = Pe / (gm1*Cv*293.0*alpha0); + Real r1 = Pc / (gm1 * Cv * 293.0); + Real rmid = Pe / (gm1 * Cv * 293.0 * alpha0); init_sg_IdealGas(0, &igra2, gm1, Cv, enabled, vals); - REQUIRE(isClose(Pe, igra2.PressureFromDensityTemperature(alpha0*rmid, 293.0), 1.e-12)); + REQUIRE(isClose(Pe, igra2.PressureFromDensityTemperature(alpha0 * rmid, 293.0), + 1.e-12)); REQUIRE(isClose(Pc, igra2.PressureFromDensityTemperature(r1, 293.0), 1.e-12)); THEN("The modified EOS should produce equivalent results") { compare_two_eoss(igsh, ig); From e043200b382651ca9601a516fa9a209bc04723ef Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Tue, 14 Jun 2022 16:15:15 -0600 Subject: [PATCH 14/34] try a little cleaner format script --- .github/workflows/formatting.yml | 4 ++-- utils/scripts/format.sh | 31 +++++++++++++++++-------------- 2 files changed, 19 insertions(+), 16 deletions(-) diff --git a/.github/workflows/formatting.yml b/.github/workflows/formatting.yml index 86708de546..a1c441815d 100644 --- a/.github/workflows/formatting.yml +++ b/.github/workflows/formatting.yml @@ -18,6 +18,6 @@ jobs: run: export DEBIAN_FRONTEND=noninteractive - name: install dependencies run: | - sudo apt-get install -y --force-yes -qq git clang-format + sudo apt-get install -y --force-yes -qq git clang-format-12 - name: check formatting - run: find . -regex '.*\.\(cpp\|hpp\)' | xargs clang-format -style=file -i && git diff --exit-code --ignore-submodules + run: find . -regex '.*\.\(cpp\|hpp\)' | xargs clang-format-12 -style=file -i && git diff --exit-code --ignore-submodules diff --git a/utils/scripts/format.sh b/utils/scripts/format.sh index 797158cb6e..188e77c13b 100755 --- a/utils/scripts/format.sh +++ b/utils/scripts/format.sh @@ -1,7 +1,7 @@ #!/bin/bash #------------------------------------------------------------------------------ -# © 2021-2022. Triad National Security, LLC. All rights reserved. This +# © 2022. Triad National Security, LLC. All rights reserved. This # program was produced under U.S. Government contract 89233218CNA000001 # for Los Alamos National Laboratory (LANL), which is operated by Triad # National Security, LLC for the U.S. Department of Energy/National @@ -15,27 +15,30 @@ #------------------------------------------------------------------------------ -# Run from anywhere in the repo to format your code. - -# You can set the clang-format binary you're using with, e.g., -# CFM=clang-format-13 format.sh -: ${CFM:=$(command -v clang-format)} +: ${CFM:=clang-format} +: ${VERBOSE:=0} if ! command -v ${CFM} &> /dev/null; then - >&2 echo "Error: No clang format found!" + >&2 echo "Error: No clang format found! Looked for ${CFM}" exit 1 +else + CFM=$(command -v ${CFM}) + echo "Clang format found: ${CFM}" fi # clang format major version TARGET_CF_VRSN=12 -CF_VRSN=$(${CFM} --version | cut -d ' ' -f 4 | cut -d '.' -f 1) - -if [ "${CF_VRSN}" != "${TARGET_CF_VRSN}" ]; then - >&2 echo "Warning! Your clang format version ${CF_VRSN} is not the same as the pinned version ${TARGET_CF_VRSN}." - >&2 echo "Results may be unstable." -fi +CF_VRSN=$(${CFM} --version) +echo "Note we assume clang format version ${TARGET_CF_VRSN}." +echo "You are using ${CF_VRSN}." +echo "If these differ, results may not be stable." +echo "Formatting..." REPO=$(git rev-parse --show-toplevel) -for f in $(git grep --cached -Il res -- :/*.hpp :/*.cpp); do +for f in $(git grep --untracked -ail res -- :/*.hpp :/*.cpp); do + if [ ${VERBOSE} -ge 1 ]; then + echo ${f} + fi ${CFM} -i ${REPO}/${f} done +echo "...Done" From 18b0f2ebcd81d0b3fcbdb6752524c700d092574c Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 14 Jun 2022 18:20:15 -0600 Subject: [PATCH 15/34] Adding modifications to fill eos to address jhps comments. --- singularity-eos/eos/modifiers/ramps_eos.hpp | 22 ++++++++++++++++++--- 1 file changed, 19 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 81bb637119..38e37cc1f8 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -72,7 +72,7 @@ class SAPRampEOS : public EosBase> { // move semantics ensures dynamic memory comes along for the ride SAPRampEOS(T &&t, const Real r0, const Real a, const Real b, const Real c) : t_(std::forward(t)), r0_(r0), a_(a), b_(b), c_(c), - rmid_(r0 * (a - b * c) / (a - b)) { + rmid_(r0 * (a - b * c) / (a - b)), Pmid_(a * (rmid_ / r0 - 1.0)) { // add input parameter checks to ensure validity of the ramp assert(r0 >= 0.0); assert(a > 0.0); @@ -92,6 +92,13 @@ class SAPRampEOS : public EosBase> { return p_ramp; } + PORTABLE_INLINE_FUNCTION + Real get_ramp_density(Real P) const { + const Real rho_ramp{P < Pmid_ ? r0_ * (P / a_ + 1.0) + : r0_ * (P / b_ + c_)}; + return rho_ramp; + } + PORTABLE_INLINE_FUNCTION Real get_ramp_dpdrho(Real rho) const { const Real dpdr{rho < r0_ ? 0.0 : rho < rmid_ ? a_ / r0_ : b_ / r0_}; @@ -166,8 +173,6 @@ class SAPRampEOS : public EosBase> { PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, const unsigned long output, Real *lambda = nullptr) const { - // density must be an input - assert(!(output & thermalqs::density)); // output passed into internal filleos can't include pressure const unsigned long ramp_out = output & ~thermalqs::pressure; // if pressure is output, calculate it first @@ -182,6 +187,16 @@ class SAPRampEOS : public EosBase> { } // call internal filleos t_.FillEos(rho, temp, energy, press, cv, bmod, ramp_out, lambda); + // fill ramp density + Real rho_ramp {rho}; + if (output & thermalqs::density) { + assert(!(output & thermalqs::pressure)); + rho_ramp = get_ramp_density(press); + } + rho = rho_ramp > rho ? rho_ramp : rho; + // bulk modulus + bmod = BulkModulusFromDensityInternalEnergy(rho, energy, lambda); + return; } PORTABLE_INLINE_FUNCTION @@ -217,6 +232,7 @@ class SAPRampEOS : public EosBase> { Real b_; Real c_; Real rmid_; + Real Pmid_; }; } // namespace singularity From 4949ab4fc8b5e60ee0d193c8e6ffcb62a459e125 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 14 Jun 2022 18:26:00 -0600 Subject: [PATCH 16/34] format --- singularity-eos/eos/modifiers/ramps_eos.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 38e37cc1f8..b9446a7b3a 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -94,8 +94,7 @@ class SAPRampEOS : public EosBase> { PORTABLE_INLINE_FUNCTION Real get_ramp_density(Real P) const { - const Real rho_ramp{P < Pmid_ ? r0_ * (P / a_ + 1.0) - : r0_ * (P / b_ + c_)}; + const Real rho_ramp{P < Pmid_ ? r0_ * (P / a_ + 1.0) : r0_ * (P / b_ + c_)}; return rho_ramp; } @@ -188,7 +187,7 @@ class SAPRampEOS : public EosBase> { // call internal filleos t_.FillEos(rho, temp, energy, press, cv, bmod, ramp_out, lambda); // fill ramp density - Real rho_ramp {rho}; + Real rho_ramp{rho}; if (output & thermalqs::density) { assert(!(output & thermalqs::pressure)); rho_ramp = get_ramp_density(press); From 1b88341629390b37bebf2150dd766ff0cf70b0db Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Wed, 15 Jun 2022 12:28:03 -0600 Subject: [PATCH 17/34] Sign change to reflect that ramp model should decrease density rather than increase it. --- singularity-eos/eos/modifiers/ramps_eos.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index b9446a7b3a..8ee70ec41c 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -192,7 +192,7 @@ class SAPRampEOS : public EosBase> { assert(!(output & thermalqs::pressure)); rho_ramp = get_ramp_density(press); } - rho = rho_ramp > rho ? rho_ramp : rho; + rho = rho_ramp < rho ? rho_ramp : rho; // bulk modulus bmod = BulkModulusFromDensityInternalEnergy(rho, energy, lambda); return; From e595bf83c03bc0a4d2d1d3937c589eb9d131e032 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 16 Jun 2022 15:00:25 -0600 Subject: [PATCH 18/34] Update format.sh --- utils/scripts/format.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/scripts/format.sh b/utils/scripts/format.sh index 188e77c13b..4d4e01841d 100755 --- a/utils/scripts/format.sh +++ b/utils/scripts/format.sh @@ -1,7 +1,7 @@ #!/bin/bash #------------------------------------------------------------------------------ -# © 2022. Triad National Security, LLC. All rights reserved. This +# © 2021-2022. Triad National Security, LLC. All rights reserved. This # program was produced under U.S. Government contract 89233218CNA000001 # for Los Alamos National Laboratory (LANL), which is operated by Triad # National Security, LLC for the U.S. Department of Energy/National From ad0888821d305e488d3286627ce2a604096e1c9f Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 16 Jun 2022 15:17:37 -0600 Subject: [PATCH 19/34] add more sphinx docs --- doc/sphinx/index.rst | 1 + doc/sphinx/src/modifiers.rst | 9 +++++++++ 2 files changed, 10 insertions(+) create mode 100644 doc/sphinx/src/modifiers.rst diff --git a/doc/sphinx/index.rst b/doc/sphinx/index.rst index 3583fd9f81..786a06c84c 100644 --- a/doc/sphinx/index.rst +++ b/doc/sphinx/index.rst @@ -19,6 +19,7 @@ Documentation approved for unlimited release. LA-UR-21-31131. src/building src/using-eos src/models + src/modifiers src/using-closures src/contributing src/sphinx-doc diff --git a/doc/sphinx/src/modifiers.rst b/doc/sphinx/src/modifiers.rst new file mode 100644 index 0000000000..c6fc3c067d --- /dev/null +++ b/doc/sphinx/src/modifiers.rst @@ -0,0 +1,9 @@ +.. _modifiers: + +EOS Modifiers +============== + +An equation of state modifier (perhaps intuitively) *modifies* an +equation of state in some pre-prescribed way. For example, one can use +a modifier to shift the assumed zero-point energy of the equation of +state, add a unit system, or account for porosity. From c94e901affdd667a79c7ae3b02a75f02896f2a1f Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Wed, 29 Jun 2022 10:12:23 -0600 Subject: [PATCH 20/34] Better separation of ramp test. --- test/test_eos_unit.cpp | 65 +++++++++++++++++++++++++++++++----------- 1 file changed, 49 insertions(+), 16 deletions(-) diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 265f61d569..bb45aab1a8 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -248,33 +248,66 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { EOS ig = IdealGas(gm1, Cv); EOS igsh = ScaledEOS(IdealGas(gm1, Cv), 1.0); EOS igsc = ShiftedEOS(IdealGas(gm1, Cv), 0.0); - EOS igra; // = SAPRampEOS(IdealGas(gm1, Cv), 1.e9, 1.0, 2.0, 1.0); + EOS igra; // test out the c interface int enabled[4] = {0, 0, 1, 0}; Real vals[6] = {0.0, 0.0, 1.e9, 1.0, 2.0, 1.0}; Real rho0 = 1.e6 / (gm1 * Cv * 293.0); init_sg_IdealGas(0, &igra, gm1, Cv, enabled, vals); - EOS igra2; - enabled[2] = 0; - enabled[3] = 1; - Real Pe = 2.e6, Pc = 1.e7; - Real alpha0 = 0.98; - vals[2] = alpha0; - vals[3] = Pe; - vals[4] = Pc; - vals[5] = 0.0; - Real r1 = Pc / (gm1 * Cv * 293.0); - Real rmid = Pe / (gm1 * Cv * 293.0 * alpha0); - init_sg_IdealGas(0, &igra2, gm1, Cv, enabled, vals); - REQUIRE(isClose(Pe, igra2.PressureFromDensityTemperature(alpha0 * rmid, 293.0), - 1.e-12)); - REQUIRE(isClose(Pc, igra2.PressureFromDensityTemperature(r1, 293.0), 1.e-12)); THEN("The modified EOS should produce equivalent results") { compare_two_eoss(igsh, ig); compare_two_eoss(igsc, ig); compare_two_eoss(igra, ig); } } + WHEN("We construct a ramp from a p-alpha model") { + const Real Pe = 2.e6, Pc = 2.5e6; + const Real alpha0 = 1.1; + const Real T0 = 293.0; + int enabled[4] = {0, 0, 0, 1}; + Real vals[6] = {0.0, 0.0, alpha0, Pe, Pc, 0.0}; + const Real rho0 = 1.e6 / (gm1 * Cv * T0); + EOS igra; + const Real r0 = rho0 / alpha0; + const Real r1 = Pc / (gm1 * Cv * T0); + const Real rmid = Pe / (gm1 * Cv * T0 * alpha0); + init_sg_IdealGas(0, &igra, gm1, Cv, enabled, vals); + // construct ramp params and evaluate directly for test + const Real a = r0 * Pe / (rmid - r0); + const Real b = r0 * (Pc - Pe) / (r1 - rmid); + const Real c = (Pc * rmid - Pe * r1) / (r0 * (Pc - Pe)); + // density in the middle of the first slope + const Real rho_t1 = 0.5 * (r0 + rmid); + // density in the middle of the second slope + const Real rho_t2 = 0.5 * (rmid + r1); + // P (alpha0*rho_t1) note that r0 = rho0 / alpha0 + const Real Prhot1 = a * (rho_t1 / r0 - 1.0); + // P (alpha0*rho_t2) + const Real Prhot2 = b * (rho_t2 / r0 - c); + // P (alpha0*rho_t1) note that r0 = rho0 / alpha0 + const Real Prhot1i = gm1 * Cv * rho_t1 * T0; + // P (alpha0*rho_t2) + const Real Prhot2i = gm1 * Cv * rho_t2 * T0; + printf("Pe=%e Pe_lb=%e\n", Pe, 1.e6 * rho_t2 / (rho0 + (1.0-alpha0)*rho_t2)); + printf("rho1=%e\nrho2=%e\n", rho_t1, rho_t2); + printf("P1=%e P2=%e P1noa=%e P2noa=%e \n", Prhot1, Prhot2, Prhot1i, Prhot2i); + printf("P1eval=%e P2eval=%e P1evalnoa=%e P2evalnoa=%e\n", + igra.PressureFromDensityTemperature(rho_t1, T0), + igra.PressureFromDensityTemperature(rho_t2, T0), + igra.PressureFromDensityTemperature(alpha0*rho_t1, T0), + igra.PressureFromDensityTemperature(alpha0*rho_t2, T0)); + THEN("We obtain Pe and Pc when evaluating at mid and upper density values") { + REQUIRE(isClose(Pe, igra.PressureFromDensityTemperature(alpha0 * rmid, 293.0), + 1.e-12)); + REQUIRE(isClose(Pc, igra.PressureFromDensityTemperature(r1, 293.0), 1.e-12)); + // also check pressures on ramp + REQUIRE(isClose(Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0), 1.e-12)); + REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12)); + } + } + WHEN("We construct a bi-linear ramp") { + + } } } From 011c4133005e791d4327ceb3233c1db6bd233ed6 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Wed, 29 Jun 2022 16:22:09 -0600 Subject: [PATCH 21/34] Add tests in both ramp regions and consistency test. --- test/test_eos_unit.cpp | 26 ++++++++------------------ 1 file changed, 8 insertions(+), 18 deletions(-) diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 245bb92fff..f5c3f37b1b 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -260,8 +260,8 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { } } WHEN("We construct a ramp from a p-alpha model") { - const Real Pe = 2.e6, Pc = 2.5e6; - const Real alpha0 = 1.1; + const Real Pe = 5.e7, Pc = 1.e8; + const Real alpha0 = 1.5; const Real T0 = 293.0; int enabled[4] = {0, 0, 0, 1}; Real vals[6] = {0.0, 0.0, alpha0, Pe, Pc, 0.0}; @@ -270,6 +270,8 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { const Real r0 = rho0 / alpha0; const Real r1 = Pc / (gm1 * Cv * T0); const Real rmid = Pe / (gm1 * Cv * T0 * alpha0); + // P(alpha0 * rmid) + const Real P_armid = alpha0 * gm1 * Cv * rmid * T0; init_sg_IdealGas(0, &igra, gm1, Cv, enabled, vals); // construct ramp params and evaluate directly for test const Real a = r0 * Pe / (rmid - r0); @@ -279,26 +281,14 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { const Real rho_t1 = 0.5 * (r0 + rmid); // density in the middle of the second slope const Real rho_t2 = 0.5 * (rmid + r1); - // P (alpha0*rho_t1) note that r0 = rho0 / alpha0 + // P (rho_t1) note that r0 = rho0 / alpha0 const Real Prhot1 = a * (rho_t1 / r0 - 1.0); - // P (alpha0*rho_t2) + // P (rho_t2) const Real Prhot2 = b * (rho_t2 / r0 - c); - // P (alpha0*rho_t1) note that r0 = rho0 / alpha0 - const Real Prhot1i = gm1 * Cv * rho_t1 * T0; - // P (alpha0*rho_t2) - const Real Prhot2i = gm1 * Cv * rho_t2 * T0; - printf("Pe=%e Pe_lb=%e\n", Pe, 1.e6 * rho_t2 / (rho0 + (1.0-alpha0)*rho_t2)); - printf("rho1=%e\nrho2=%e\n", rho_t1, rho_t2); - printf("P1=%e P2=%e P1noa=%e P2noa=%e \n", Prhot1, Prhot2, Prhot1i, Prhot2i); - printf("P1eval=%e P2eval=%e P1evalnoa=%e P2evalnoa=%e\n", - igra.PressureFromDensityTemperature(rho_t1, T0), - igra.PressureFromDensityTemperature(rho_t2, T0), - igra.PressureFromDensityTemperature(alpha0*rho_t1, T0), - igra.PressureFromDensityTemperature(alpha0*rho_t2, T0)); THEN("We obtain Pe and Pc when evaluating at mid and upper density values") { - REQUIRE(isClose(Pe, igra.PressureFromDensityTemperature(alpha0 * rmid, 293.0), + REQUIRE(isClose(P_armid, igra.PressureFromDensityTemperature(rmid, 293.0), 1.e-12)); - REQUIRE(isClose(Pc, igra.PressureFromDensityTemperature(r1, 293.0), 1.e-12)); + //REQUIRE(isClose(Pc, igra.PressureFromDensityTemperature(r1, 293.0), 1.e-12)); // also check pressures on ramp REQUIRE(isClose(Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0), 1.e-12)); REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12)); From 32ee97c36fccae077e2a0f9ff66f8526d8b630ac Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Wed, 29 Jun 2022 16:55:28 -0600 Subject: [PATCH 22/34] Address some of Jonahs concerns. --- singularity-eos/eos/modifiers/ramps_eos.hpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 8ee70ec41c..2d5e1751ff 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -39,7 +39,8 @@ void pAlpha2SAPRampParams(const T &eos, const Real alpha0, const Real Pe, const Real rho0, T0, sie0, P0, cv0, bmod0, dpde0, dvdt0; Real rmid, r1; eos.ValuesAtReferenceState(rho0, T0, sie0, P0, cv0, bmod0, dpde0, dvdt0); - // calculate r0 + // calculate r0, ensure alpha0 > 1 + assert(alpha0 > 1.0); r0 = rho0 / alpha0; // calculate rmid auto rmid_func = [&](const Real x) { @@ -74,7 +75,7 @@ class SAPRampEOS : public EosBase> { : t_(std::forward(t)), r0_(r0), a_(a), b_(b), c_(c), rmid_(r0 * (a - b * c) / (a - b)), Pmid_(a * (rmid_ / r0 - 1.0)) { // add input parameter checks to ensure validity of the ramp - assert(r0 >= 0.0); + assert(r0 > 0.0); assert(a > 0.0); assert(b >= 0); assert(a != b); @@ -179,9 +180,9 @@ class SAPRampEOS : public EosBase> { // maybe switch case on preferred input and check for not output of the input // for now sie lookup is prioritized if (!(output & thermalqs::specific_internal_energy)) { - press = this->PressureFromDensityInternalEnergy(rho, energy, lambda); + press = PressureFromDensityInternalEnergy(rho, energy, lambda); } else if (!(output & thermalqs::temperature)) { - press = this->PressureFromDensityTemperature(rho, temp, lambda); + press = PressureFromDensityTemperature(rho, temp, lambda); } } // call internal filleos @@ -206,7 +207,7 @@ class SAPRampEOS : public EosBase> { PORTABLE_FUNCTION void PrintParams() const { t_.PrintParams(); - printf("r0=%e\na=%e\nb=%e\nc=%e\n", r0_, a_, b_, c_); + printf("Ramp Params:\nr0=%.14e\na=%.14e\nb=%.14e\nc=%.14e\n", r0_, a_, b_, c_); } PORTABLE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, From 905901679768b714056794d71f5179916937897b Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Wed, 29 Jun 2022 17:40:23 -0600 Subject: [PATCH 23/34] Formatting. --- test/test_eos_unit.cpp | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index f5c3f37b1b..57627514de 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -286,17 +286,15 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { // P (rho_t2) const Real Prhot2 = b * (rho_t2 / r0 - c); THEN("We obtain Pe and Pc when evaluating at mid and upper density values") { - REQUIRE(isClose(P_armid, igra.PressureFromDensityTemperature(rmid, 293.0), - 1.e-12)); - //REQUIRE(isClose(Pc, igra.PressureFromDensityTemperature(r1, 293.0), 1.e-12)); - // also check pressures on ramp - REQUIRE(isClose(Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0), 1.e-12)); - REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12)); + REQUIRE( + isClose(P_armid, igra.PressureFromDensityTemperature(rmid, 293.0), 1.e-12)); + // REQUIRE(isClose(Pc, igra.PressureFromDensityTemperature(r1, 293.0), 1.e-12)); + // also check pressures on ramp + REQUIRE(isClose(Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0), 1.e-12)); + REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12)); } } - WHEN("We construct a bi-linear ramp") { - - } + WHEN("We construct a bi-linear ramp") {} } } From ee7fcab07e2ee7fe3be94872dc8b85c73f96e1dc Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Thu, 30 Jun 2022 09:54:03 -0600 Subject: [PATCH 24/34] Add additional outside ramp domain tests and bulk modulus tests. --- test/test_eos_unit.cpp | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 57627514de..3f463361bc 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -285,6 +285,10 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { const Real Prhot1 = a * (rho_t1 / r0 - 1.0); // P (rho_t2) const Real Prhot2 = b * (rho_t2 / r0 - c); + // bmod (rho_t1) + const Real bmodrt1 = rho_t1 * a / r0; + // bmod (rho_t2) + const Real bmodrt2 = rho_t2 * b / r0; THEN("We obtain Pe and Pc when evaluating at mid and upper density values") { REQUIRE( isClose(P_armid, igra.PressureFromDensityTemperature(rmid, 293.0), 1.e-12)); @@ -292,9 +296,19 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { // also check pressures on ramp REQUIRE(isClose(Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0), 1.e-12)); REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12)); + // check pressure below and beyond ramp matches unmodified ideal gas + REQUIRE(isClose(0.8*r0 * gm1 * Cv * T0, igra.PressureFromDensityTemperature(0.8*r0, T0), 1.e-12)); + REQUIRE(isClose(1.2*r1 * gm1 * Cv * T0, igra.PressureFromDensityTemperature(1.2*r1, T0), 1.e-12)); + // check bulk moduli on both pieces of ramp + REQUIRE(isClose(bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0), 1.e-12)); + REQUIRE(isClose(bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0), 1.e-12)); + // check bulk modulus below and beyond ramp matches unmodified ideal gas + REQUIRE(isClose(0.8*r0 * gm1 * (gm1 + 1.0) * Cv * T0, + igra.BulkModulusFromDensityTemperature(0.8*r0, T0), 1.e-12)); + REQUIRE(isClose(1.2*r1 * gm1 * (gm1 + 1.0) * Cv * T0, + igra.BulkModulusFromDensityTemperature(1.2*r1, T0), 1.e-12)); } } - WHEN("We construct a bi-linear ramp") {} } } From da18f8563061db90428b7ab34978a67ed6884308 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Thu, 30 Jun 2022 10:32:43 -0600 Subject: [PATCH 25/34] formatting --- test/test_eos_unit.cpp | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 3f463361bc..48d626be77 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -296,17 +296,21 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { // also check pressures on ramp REQUIRE(isClose(Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0), 1.e-12)); REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12)); - // check pressure below and beyond ramp matches unmodified ideal gas - REQUIRE(isClose(0.8*r0 * gm1 * Cv * T0, igra.PressureFromDensityTemperature(0.8*r0, T0), 1.e-12)); - REQUIRE(isClose(1.2*r1 * gm1 * Cv * T0, igra.PressureFromDensityTemperature(1.2*r1, T0), 1.e-12)); - // check bulk moduli on both pieces of ramp - REQUIRE(isClose(bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0), 1.e-12)); - REQUIRE(isClose(bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0), 1.e-12)); - // check bulk modulus below and beyond ramp matches unmodified ideal gas - REQUIRE(isClose(0.8*r0 * gm1 * (gm1 + 1.0) * Cv * T0, - igra.BulkModulusFromDensityTemperature(0.8*r0, T0), 1.e-12)); - REQUIRE(isClose(1.2*r1 * gm1 * (gm1 + 1.0) * Cv * T0, - igra.BulkModulusFromDensityTemperature(1.2*r1, T0), 1.e-12)); + // check pressure below and beyond ramp matches unmodified ideal gas + REQUIRE(isClose(0.8 * r0 * gm1 * Cv * T0, + igra.PressureFromDensityTemperature(0.8 * r0, T0), 1.e-12)); + REQUIRE(isClose(1.2 * r1 * gm1 * Cv * T0, + igra.PressureFromDensityTemperature(1.2 * r1, T0), 1.e-12)); + // check bulk moduli on both pieces of ramp + REQUIRE( + isClose(bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0), 1.e-12)); + REQUIRE( + isClose(bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0), 1.e-12)); + // check bulk modulus below and beyond ramp matches unmodified ideal gas + REQUIRE(isClose(0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0, + igra.BulkModulusFromDensityTemperature(0.8 * r0, T0), 1.e-12)); + REQUIRE(isClose(1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0, + igra.BulkModulusFromDensityTemperature(1.2 * r1, T0), 1.e-12)); } } } From 5db2cc3c6c7ebe204b27ab254601ccb5b96dfd3e Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Fri, 1 Jul 2022 13:32:21 -0600 Subject: [PATCH 26/34] Testing improvements as suggested by Jeff. --- test/test_eos_unit.cpp | 113 +++++++++++++++++++++++++++++++---------- 1 file changed, 87 insertions(+), 26 deletions(-) diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 48d626be77..21334660b1 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -73,28 +73,58 @@ inline void compare_two_eoss(E1 &&test_e, E2 &&ref_e) { // modifiers that can be initialized in such a way as to // be equivalent of an unmodified eos. Best used with analytic // eoss. - REQUIRE(isClose(test_e.TemperatureFromDensityInternalEnergy(1, 1), - ref_e.TemperatureFromDensityInternalEnergy(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.InternalEnergyFromDensityTemperature(1, 1), - ref_e.InternalEnergyFromDensityTemperature(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.PressureFromDensityInternalEnergy(1, 1), - ref_e.PressureFromDensityInternalEnergy(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.SpecificHeatFromDensityInternalEnergy(1, 1), - ref_e.SpecificHeatFromDensityInternalEnergy(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.BulkModulusFromDensityInternalEnergy(1, 1), - ref_e.BulkModulusFromDensityInternalEnergy(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.GruneisenParamFromDensityInternalEnergy(1, 1), - ref_e.GruneisenParamFromDensityInternalEnergy(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.PressureFromDensityTemperature(1, 1), - ref_e.PressureFromDensityTemperature(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.SpecificHeatFromDensityTemperature(1, 1), - ref_e.SpecificHeatFromDensityTemperature(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.BulkModulusFromDensityTemperature(1, 1), - ref_e.BulkModulusFromDensityTemperature(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.GruneisenParamFromDensityTemperature(1, 1), - ref_e.GruneisenParamFromDensityTemperature(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.MinimumDensity(), ref_e.MinimumDensity(), 1.e-15)); - REQUIRE(isClose(test_e.MinimumTemperature(), ref_e.MinimumTemperature(), 1.e-15)); + INFO("reference T: " << ref_e.TemperatureFromDensityInternalEnergy(1, 1) << " test T: " + << test_e.TemperatureFromDensityInternalEnergy(1, 1)); + CHECK(isClose(test_e.TemperatureFromDensityInternalEnergy(1, 1), + ref_e.TemperatureFromDensityInternalEnergy(1, 1), 1.e-15)); + INFO("reference sie: " << ref_e.InternalEnergyFromDensityTemperature(1, 1) + << " test sie: " + << test_e.InternalEnergyFromDensityTemperature(1, 1)); + CHECK(isClose(test_e.InternalEnergyFromDensityTemperature(1, 1), + ref_e.InternalEnergyFromDensityTemperature(1, 1), 1.e-15)); + INFO("reference P: " << ref_e.PressureFromDensityInternalEnergy(1, 1) + << " test P: " << test_e.PressureFromDensityInternalEnergy(1, 1)); + CHECK(isClose(test_e.PressureFromDensityInternalEnergy(1, 1), + ref_e.PressureFromDensityInternalEnergy(1, 1), 1.e-15)); + INFO("reference Cv: " << ref_e.SpecificHeatFromDensityInternalEnergy(1, 1) + << " test Cv: " + << test_e.SpecificHeatFromDensityInternalEnergy(1, 1)); + CHECK(isClose(test_e.SpecificHeatFromDensityInternalEnergy(1, 1), + ref_e.SpecificHeatFromDensityInternalEnergy(1, 1), 1.e-15)); + INFO("reference bmod: " << ref_e.BulkModulusFromDensityInternalEnergy(1, 1) + << " test bmod: " + << test_e.BulkModulusFromDensityInternalEnergy(1, 1)); + CHECK(isClose(test_e.BulkModulusFromDensityInternalEnergy(1, 1), + ref_e.BulkModulusFromDensityInternalEnergy(1, 1), 1.e-15)); + INFO("reference Grun. Param.: " + << ref_e.GruneisenParamFromDensityInternalEnergy(1, 1) + << " test Grun. Param.: " << test_e.GruneisenParamFromDensityInternalEnergy(1, 1)); + CHECK(isClose(test_e.GruneisenParamFromDensityInternalEnergy(1, 1), + ref_e.GruneisenParamFromDensityInternalEnergy(1, 1), 1.e-15)); + INFO("reference P: " << ref_e.PressureFromDensityTemperature(1, 1) + << " test P: " << test_e.PressureFromDensityTemperature(1, 1)); + CHECK(isClose(test_e.PressureFromDensityTemperature(1, 1), + ref_e.PressureFromDensityTemperature(1, 1), 1.e-15)); + INFO("reference Cv: " << ref_e.SpecificHeatFromDensityTemperature(1, 1) << " test Cv: " + << test_e.SpecificHeatFromDensityTemperature(1, 1)); + CHECK(isClose(test_e.SpecificHeatFromDensityTemperature(1, 1), + ref_e.SpecificHeatFromDensityTemperature(1, 1), 1.e-15)); + INFO("reference bmod: " << ref_e.BulkModulusFromDensityTemperature(1, 1) + << " test bmod: " + << test_e.BulkModulusFromDensityTemperature(1, 1)); + CHECK(isClose(test_e.BulkModulusFromDensityTemperature(1, 1), + ref_e.BulkModulusFromDensityTemperature(1, 1), 1.e-15)); + INFO("reference Grun. Param.: " << ref_e.GruneisenParamFromDensityTemperature(1, 1) + << " test Grun. Param.: " + << test_e.GruneisenParamFromDensityTemperature(1, 1)); + CHECK(isClose(test_e.GruneisenParamFromDensityTemperature(1, 1), + ref_e.GruneisenParamFromDensityTemperature(1, 1), 1.e-15)); + INFO("reference rho min.: " << ref_e.MinimumDensity() + << " test rho min.: " << test_e.MinimumDensity()); + CHECK(isClose(test_e.MinimumDensity(), ref_e.MinimumDensity(), 1.e-15)); + INFO("reference T min.: " << ref_e.MinimumTemperature() + << " test T min.: " << test_e.MinimumTemperature()); + CHECK(isClose(test_e.MinimumTemperature(), ref_e.MinimumTemperature(), 1.e-15)); return; } @@ -289,26 +319,57 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { const Real bmodrt1 = rho_t1 * a / r0; // bmod (rho_t2) const Real bmodrt2 = rho_t2 * b / r0; - THEN("We obtain Pe and Pc when evaluating at mid and upper density values") { - REQUIRE( - isClose(P_armid, igra.PressureFromDensityTemperature(rmid, 293.0), 1.e-12)); - // REQUIRE(isClose(Pc, igra.PressureFromDensityTemperature(r1, 293.0), 1.e-12)); + THEN("P_eos(alpha_0*rmid, T0) = P_ramp(rmid,T0)") { + INFO("P_eos(alpha_0*rmid, T0): " + << P_armid + << " P_ramp(rmid, T0): " << igra.PressureFromDensityTemperature(rmid, T0)); + REQUIRE(isClose(P_armid, igra.PressureFromDensityTemperature(rmid, T0), 1.e-12)); + } + THEN("We obtain correct ramp behavior in P(rho) for rho r1") { // also check pressures on ramp + INFO("reference P((r0+rmid)/2, T0): " + << Prhot1 << " test P((r0+rmid)/2, T0): " + << igra.PressureFromDensityTemperature(rho_t1, T0)); REQUIRE(isClose(Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0), 1.e-12)); + INFO("reference P((rmid+r1)/2, T0): " + << Prhot2 << " test P((rmid+r1)/2, T0): " + << igra.PressureFromDensityTemperature(rho_t2, T0)); REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12)); // check pressure below and beyond ramp matches unmodified ideal gas + INFO("reference P(0.8*r0, T0): " + << 0.8 * r0 * gm1 * Cv * T0 << " test P(0.8*r0, T0): " + << igra.PressureFromDensityTemperature(0.8 * r0, T0)); REQUIRE(isClose(0.8 * r0 * gm1 * Cv * T0, igra.PressureFromDensityTemperature(0.8 * r0, T0), 1.e-12)); + INFO("reference P(1.2*r1, T0): " + << 1.2 * r1 * gm1 * Cv * T0 << " test P(1.2*r1, T0): " + << igra.PressureFromDensityTemperature(1.2 * r1, T0)); REQUIRE(isClose(1.2 * r1 * gm1 * Cv * T0, igra.PressureFromDensityTemperature(1.2 * r1, T0), 1.e-12)); + } + THEN("We obtain correct ramp behavior in bmod(rho) for rho r1") { // check bulk moduli on both pieces of ramp + INFO("reference bmod((r0+rmid)/2, T0): " + << bmodrt1 << " test bmod((r0+rmid)/2, T0): " + << igra.BulkModulusFromDensityTemperature(rho_t1, T0)); REQUIRE( isClose(bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0), 1.e-12)); + INFO("reference bmod((rmid+r1)/2, T0): " + << bmodrt2 << " test bmod((rmid+r1)/2, T0): " + << igra.BulkModulusFromDensityTemperature(rho_t2, T0)); REQUIRE( isClose(bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0), 1.e-12)); // check bulk modulus below and beyond ramp matches unmodified ideal gas + INFO("reference bmod(0.8*r0, T0): " + << 0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0 << " test bmod(0.8*r0, T0): " + << igra.BulkModulusFromDensityTemperature(0.8 * r0, T0)); REQUIRE(isClose(0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0, igra.BulkModulusFromDensityTemperature(0.8 * r0, T0), 1.e-12)); + INFO("reference bmod(1.2*r1, T0): " + << 1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0 << " test bmod(1.2*r1, T0): " + << igra.BulkModulusFromDensityTemperature(1.2 * r1, T0)); REQUIRE(isClose(1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0, igra.BulkModulusFromDensityTemperature(1.2 * r1, T0), 1.e-12)); } From 7126128211fa776653d36471bce69f2fbc7b9375 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Fri, 1 Jul 2022 14:27:39 -0600 Subject: [PATCH 27/34] Add a robust utils header. Making divisions robust in ramp. --- singularity-eos/CMakeLists.txt | 1 + singularity-eos/base/robust_utils.hpp | 57 +++++++++++++++++++++ singularity-eos/eos/modifiers/ramps_eos.hpp | 24 +++++---- 3 files changed, 72 insertions(+), 10 deletions(-) create mode 100644 singularity-eos/base/robust_utils.hpp diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 6684d5c2e1..f044efcca2 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -16,6 +16,7 @@ set(EOS_HEADERS base/constants.hpp base/eos_error.hpp base/fast-math/logs.hpp + base/robust_utils.hpp base/root-finding-1d/root_finding.hpp eos/eos.hpp eos/eos_builder.hpp diff --git a/singularity-eos/base/robust_utils.hpp b/singularity-eos/base/robust_utils.hpp new file mode 100644 index 0000000000..01ec09bbbd --- /dev/null +++ b/singularity-eos/base/robust_utils.hpp @@ -0,0 +1,57 @@ +//------------------------------------------------------------------------------ +// © 2021-2022. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#ifndef SINGULARITY_EOS_BASE_ROBUST_UTILS_HPP_ +#define SINGULARITY_EOS_BASE_ROBUST_UTILS_HPP_ + +#include +#include + +namespace singularity { +namespace robust { + +template +PORTABLE_FORCEINLINE_FUNCTION constexpr auto SMALL() { + return 10 * std::numeric_limits::min(); +} + +template +PORTABLE_FORCEINLINE_FUNCTION constexpr auto EPS() { + return 10 * std::numeric_limits::epsilon(); +} + +template +PORTABLE_FORCEINLINE_FUNCTION auto make_positive(const T val) { + return std::max(val, EPS()); +} + +PORTABLE_FORCEINLINE_FUNCTION +Real make_bounded(const Real val, const Real vmin, const Real vmax) { + return std::min(std::max(val, vmin + EPS()), vmax * (1.0 - EPS())); +} + +template +PORTABLE_FORCEINLINE_FUNCTION int sgn(const T &val) { + return (T(0) <= val) - (val < T(0)); +} + +template +PORTABLE_FORCEINLINE_FUNCTION auto ratio(const A &a, const B &b) { + return a / (b + sgn(b) * SMALL()); +} + +} // namespace robust +} // namespace singularity + +#endif // SINGULARITY_EOS_BASE_ROBUST_UTILS_HPP_ diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 2d5e1751ff..a29a09acc8 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -26,12 +26,15 @@ #include #include #include +#include #include namespace singularity { using namespace eos_base; +using singularity::robust::ratio; + template void pAlpha2SAPRampParams(const T &eos, const Real alpha0, const Real Pe, const Real Pc, Real &r0, Real &a, Real &b, Real &c) { @@ -41,7 +44,7 @@ void pAlpha2SAPRampParams(const T &eos, const Real alpha0, const Real Pe, const eos.ValuesAtReferenceState(rho0, T0, sie0, P0, cv0, bmod0, dpde0, dvdt0); // calculate r0, ensure alpha0 > 1 assert(alpha0 > 1.0); - r0 = rho0 / alpha0; + r0 = ratio(rho0, alpha0); // calculate rmid auto rmid_func = [&](const Real x) { return eos.PressureFromDensityTemperature(alpha0 * x, T0); @@ -50,7 +53,7 @@ void pAlpha2SAPRampParams(const T &eos, const Real alpha0, const Real Pe, const // get upper bound to density informed by the reference // bulk modulus const Real max_exp_arg = std::log(std::numeric_limits::max() * 0.99); - const Real exp_arg = std::min(max_exp_arg, (2.0 * Pc - P0) / bmod0); + const Real exp_arg = std::min(max_exp_arg, ratio((2.0 * Pc - P0), bmod0)); const Real rho_ub = rho0 * std::exp(exp_arg); // finds where rmid_func = Pe RootFinding1D::findRoot(rmid_func, Pe, rho0, r0, rho_ub, 1.e-12, 1.e-12, rmid, co); @@ -59,11 +62,11 @@ void pAlpha2SAPRampParams(const T &eos, const Real alpha0, const Real Pe, const // finds where r1_func = Pc RootFinding1D::findRoot(r1_func, Pc, rmid, r0, rho_ub, 1.e-12, 1.e-12, r1, co); // a - a = r0 * Pe / (rmid - r0); + a = ratio(r0 * Pe, rmid - r0); // b - b = r0 * (Pc - Pe) / (r1 - rmid); + b = ratio(r0 * (Pc - Pe), r1 - rmid); // c - c = (Pc * rmid - Pe * r1) / (r0 * (Pc - Pe)); + c = ratio(Pc * rmid - Pe * r1, r0 * (Pc - Pe)); return; } @@ -87,21 +90,22 @@ class SAPRampEOS : public EosBase> { PORTABLE_INLINE_FUNCTION Real get_ramp_pressure(Real rho) const { - const Real p_ramp{rho < r0_ ? 0.0 - : rho < rmid_ ? a_ * (rho / r0_ - 1.0) - : b_ * (rho / r0_ - c_)}; + const Real r_r0{ratio(rho, r0_)}; + const Real p_ramp{rho < r0_ ? 0.0 + : rho < rmid_ ? a_ * (r_r0 - 1.0) + : b_ * (r_r0 - c_)}; return p_ramp; } PORTABLE_INLINE_FUNCTION Real get_ramp_density(Real P) const { - const Real rho_ramp{P < Pmid_ ? r0_ * (P / a_ + 1.0) : r0_ * (P / b_ + c_)}; + const Real rho_ramp{P < Pmid_ ? r0_ * (ratio(P, a_) + 1.0) : r0_ * (ratio(P, b_) + c_)}; return rho_ramp; } PORTABLE_INLINE_FUNCTION Real get_ramp_dpdrho(Real rho) const { - const Real dpdr{rho < r0_ ? 0.0 : rho < rmid_ ? a_ / r0_ : b_ / r0_}; + const Real dpdr{rho < r0_ ? 0.0 : rho < rmid_ ? ratio(a_, r0_) : ratio(b_, r0_)}; return dpdr; } From 936e6046b3055ac1848b0afe169e9213501efbe9 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Fri, 1 Jul 2022 14:29:22 -0600 Subject: [PATCH 28/34] formatting. --- singularity-eos/eos/modifiers/ramps_eos.hpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index a29a09acc8..bde83cfc53 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -91,15 +91,16 @@ class SAPRampEOS : public EosBase> { PORTABLE_INLINE_FUNCTION Real get_ramp_pressure(Real rho) const { const Real r_r0{ratio(rho, r0_)}; - const Real p_ramp{rho < r0_ ? 0.0 - : rho < rmid_ ? a_ * (r_r0 - 1.0) - : b_ * (r_r0 - c_)}; + const Real p_ramp{rho < r0_ ? 0.0 + : rho < rmid_ ? a_ * (r_r0 - 1.0) + : b_ * (r_r0 - c_)}; return p_ramp; } PORTABLE_INLINE_FUNCTION Real get_ramp_density(Real P) const { - const Real rho_ramp{P < Pmid_ ? r0_ * (ratio(P, a_) + 1.0) : r0_ * (ratio(P, b_) + c_)}; + const Real rho_ramp{P < Pmid_ ? r0_ * (ratio(P, a_) + 1.0) + : r0_ * (ratio(P, b_) + c_)}; return rho_ramp; } From b7038e34ddee24277bbd169ae1cdd21212417d87 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 5 Jul 2022 11:35:56 -0600 Subject: [PATCH 29/34] Minor cmake fixes and gpu build warning suppression. --- CMakeLists.txt | 2 +- example/CMakeLists.txt | 16 +++++++--------- singularity-eos/eos/modifiers/ramps_eos.hpp | 4 ++-- 3 files changed, 10 insertions(+), 12 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 275f4f7148..1923e57e67 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -211,7 +211,7 @@ if(SINGULARITY_USE_KOKKOS) #### # TODO: Once it's oppurtune, we will # not support Kokkos < 3.3 - if(Kokkos_VERSION_MAJOR VERSION_LESS "3.3") + if(${Kokkos_VERSION_MAJOR}.${Kokkos_VERSION_MINOR} VERSION_LESS "3.3") message(WARNING "`Kokkos` version [${Kokkos_VERSION_MAJOR}.${Kokkos_VERSION_MINOR}] is DEPRECATED, and `singularity-eos` will not support versions < '3.3' in the very near future.") if(SINGULARITY_USE_CUDA) get_filename_component(_compiler_exe ${CMAKE_CXX_COMPILER} NAME) diff --git a/example/CMakeLists.txt b/example/CMakeLists.txt index c716507d56..7598557677 100644 --- a/example/CMakeLists.txt +++ b/example/CMakeLists.txt @@ -15,13 +15,11 @@ add_executable(get_sound_speed_press get_sound_speed_press.cpp) target_link_libraries(get_sound_speed_press PRIVATE - eos - singularity-eos::libs - singularity-eos::flags) + singularity-eos::singularity-eos) -add_executable(get_sesame_state - get_sesame_state.cpp) -target_link_libraries(get_sesame_state PRIVATE - eos - singularity-eos::libs - singularity-eos::flags) +if(SINGULARITY_USE_EOSPAC) + add_executable(get_sesame_state + get_sesame_state.cpp) + target_link_libraries(get_sesame_state PRIVATE + singularity-eos::singularity-eos) +endif() diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index bde83cfc53..f8b363c37a 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -46,7 +46,7 @@ void pAlpha2SAPRampParams(const T &eos, const Real alpha0, const Real Pe, const assert(alpha0 > 1.0); r0 = ratio(rho0, alpha0); // calculate rmid - auto rmid_func = [&](const Real x) { + auto rmid_func = PORTABLE_LAMBDA(const Real x) { return eos.PressureFromDensityTemperature(alpha0 * x, T0); }; RootFinding1D::RootCounts co{}; @@ -58,7 +58,7 @@ void pAlpha2SAPRampParams(const T &eos, const Real alpha0, const Real Pe, const // finds where rmid_func = Pe RootFinding1D::findRoot(rmid_func, Pe, rho0, r0, rho_ub, 1.e-12, 1.e-12, rmid, co); // calculate r1 - auto r1_func = [&](const Real x) { return eos.PressureFromDensityTemperature(x, T0); }; + auto r1_func = PORTABLE_LAMBDA(const Real x) { return eos.PressureFromDensityTemperature(x, T0); }; // finds where r1_func = Pc RootFinding1D::findRoot(r1_func, Pc, rmid, r0, rho_ub, 1.e-12, 1.e-12, r1, co); // a From f3916c4b64595386be91dfa50a57328fba4e377f Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 5 Jul 2022 11:40:39 -0600 Subject: [PATCH 30/34] format --- singularity-eos/eos/modifiers/ramps_eos.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index f8b363c37a..5c74eb369e 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -58,7 +58,9 @@ void pAlpha2SAPRampParams(const T &eos, const Real alpha0, const Real Pe, const // finds where rmid_func = Pe RootFinding1D::findRoot(rmid_func, Pe, rho0, r0, rho_ub, 1.e-12, 1.e-12, rmid, co); // calculate r1 - auto r1_func = PORTABLE_LAMBDA(const Real x) { return eos.PressureFromDensityTemperature(x, T0); }; + auto r1_func = PORTABLE_LAMBDA(const Real x) { + return eos.PressureFromDensityTemperature(x, T0); + }; // finds where r1_func = Pc RootFinding1D::findRoot(r1_func, Pc, rmid, r0, rho_ub, 1.e-12, 1.e-12, r1, co); // a From fb5c9676df53f3b478b21b18c6e21394e1168a06 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 5 Jul 2022 16:12:59 -0600 Subject: [PATCH 31/34] Fix compile time errors that occur in ci and ci-like build. --- test/test_eos_unit.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 21334660b1..6d9b931450 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -263,10 +263,10 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { Real a = 1; Real b = 0; Real c = 0; - ramp_params["r0"] = r0; - ramp_params["a"] = a; - ramp_params["b"] = b; - ramp_params["c"] = c; + ramp_params["r0"].emplace(r0); + ramp_params["a"].emplace(a); + ramp_params["b"].emplace(b); + ramp_params["c"].emplace(c); modifiers[EOSBuilder::EOSModifier::SAPRamp] = ramp_params; THEN("The EOS is constructed correctly") { auto eos_ramped = EOSBuilder::buildEOS(type, base_params, modifiers); From b7a9d779571eff7e379797afdbfa8557b282ebcd Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 5 Jul 2022 16:55:52 -0600 Subject: [PATCH 32/34] Further GPU build warning suppression by adding apppriate device markings for nlambda and other associated functions in spiner eoss. --- singularity-eos/eos/eos.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos.hpp b/singularity-eos/eos/eos.hpp index 8857b234ca..5dc7692f96 100644 --- a/singularity-eos/eos/eos.hpp +++ b/singularity-eos/eos/eos.hpp @@ -529,9 +529,12 @@ class SpinerEOSDependsRhoT : public EosBase { } PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return T_(lTMin_); } + PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return _n_lambda; } - inline RootFinding1D::Status rootStatus() const { return status_; } - inline TableStatus tableStatus() const { return whereAmI_; } + PORTABLE_INLINE_FUNCTION + RootFinding1D::Status rootStatus() const { return status_; } + PORTABLE_INLINE_FUNCTION + TableStatus tableStatus() const { return whereAmI_; } RootFinding1D::RootCounts counts; void Finalize(); static std::string EosType() { return std::string("SpinerEOSDependsRhoT"); } @@ -745,6 +748,7 @@ class SpinerEOSDependsRhoSie : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return TMin(); } + PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return _n_lambda; } PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char s1[]{"SpinerEOS Parameters:"}; @@ -753,6 +757,7 @@ class SpinerEOSDependsRhoSie : public EosBase { printf("%s\n\t%s\n\t%s%i\n", s1, s2, s3, matid_); return; } + PORTABLE_INLINE_FUNCTION RootFinding1D::Status rootStatus() const { return status_; } RootFinding1D::RootCounts counts; static std::string EosType() { return std::string("SpinerEOSDependsRhoSie"); } From 5cbc54b92d0cf1445eb270f5a08d592f6e0b6344 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 5 Jul 2022 18:25:07 -0600 Subject: [PATCH 33/34] Changing name of SAPRamp to more general and descriptive BilinearRamp name --- singularity-eos/eos/eos.hpp | 2 +- singularity-eos/eos/eos_builder.cpp | 18 +++++++++--------- singularity-eos/eos/eos_builder.hpp | 10 +++++----- singularity-eos/eos/modifiers/ramps_eos.hpp | 12 ++++++------ singularity-eos/eos/singularity_eos.cpp | 18 +++++++++--------- test/test_eos_unit.cpp | 4 ++-- 6 files changed, 32 insertions(+), 32 deletions(-) diff --git a/singularity-eos/eos/eos.hpp b/singularity-eos/eos/eos.hpp index 5dc7692f96..9e7d00a4ea 100644 --- a/singularity-eos/eos/eos.hpp +++ b/singularity-eos/eos/eos.hpp @@ -1306,7 +1306,7 @@ static constexpr const auto combined_list_1 = singularity::detail::concat( full_eos_list, shifted, scaled, scaled_of_shifted, unit_or_rel); // make a ramped eos of everything static constexpr const auto ramped_all = - transform_variadic_list(combined_list_1, al{}); + transform_variadic_list(combined_list_1, al{}); // final combined list static constexpr const auto combined_list = singularity::detail::concat( full_eos_list, shifted, scaled, scaled_of_shifted, unit_or_rel, ramped_all); diff --git a/singularity-eos/eos/eos_builder.cpp b/singularity-eos/eos/eos_builder.cpp index f73fe5b28b..64ae241bda 100644 --- a/singularity-eos/eos/eos_builder.cpp +++ b/singularity-eos/eos/eos_builder.cpp @@ -27,7 +27,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par bool shifted = (modifiers.count(EOSModifier::Shifted) > 0); bool relativistic = (modifiers.count(EOSModifier::Relativistic) > 0); bool units = (modifiers.count(EOSModifier::UnitSystem) > 0); - bool ramped = (modifiers.count(EOSModifier::SAPRamp) > 0); + bool ramped = (modifiers.count(EOSModifier::BilinearRamp) > 0); if ((shifted || scaled || relativistic || units || ramped) && !(isModifiable(type))) { EOS_ERROR("Modifiers not supported for this EOS"); } @@ -59,10 +59,10 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par Real b = 0; Real c = 0; if (ramped) { - r0 = mpark::get(modifiers[EOSModifier::SAPRamp]["r0"]); - a = mpark::get(modifiers[EOSModifier::SAPRamp]["a"]); - b = mpark::get(modifiers[EOSModifier::SAPRamp]["b"]); - c = mpark::get(modifiers[EOSModifier::SAPRamp]["c"]); + r0 = mpark::get(modifiers[EOSModifier::BilinearRamp]["r0"]); + a = mpark::get(modifiers[EOSModifier::BilinearRamp]["a"]); + b = mpark::get(modifiers[EOSModifier::BilinearRamp]["b"]); + c = mpark::get(modifiers[EOSModifier::BilinearRamp]["c"]); } Real rho_unit = 1; Real sie_unit = 1; @@ -104,7 +104,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(g), cl); } - return applyShiftAndScaleAndSAPRamp(std::move(g), scaled, shifted, ramped, scale, + return applyShiftAndScaleAndBilinearRamp(std::move(g), scaled, shifted, ramped, scale, shift, r0, a, b, c); } #ifdef SPINER_USE_HDF @@ -132,7 +132,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, + return applyShiftAndScaleAndBilinearRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); } } else { @@ -156,7 +156,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, + return applyShiftAndScaleAndBilinearRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); } } @@ -179,7 +179,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return applyShiftAndScaleAndSAPRamp(std::move(s), scaled, shifted, ramped, scale, + return applyShiftAndScaleAndBilinearRamp(std::move(s), scaled, shifted, ramped, scale, shift, r0, a, b, c); } #endif diff --git a/singularity-eos/eos/eos_builder.hpp b/singularity-eos/eos/eos_builder.hpp index 56b1d009bc..9e52d0be05 100644 --- a/singularity-eos/eos/eos_builder.hpp +++ b/singularity-eos/eos/eos_builder.hpp @@ -43,7 +43,7 @@ enum class EOSType { StellarCollapse #endif }; -enum class EOSModifier { Scaled, Shifted, Relativistic, UnitSystem, SAPRamp }; +enum class EOSModifier { Scaled, Shifted, Relativistic, UnitSystem, BilinearRamp }; // evil type erasure using param_t = mpark::variant; @@ -85,8 +85,8 @@ EOS makeRelativistic(T &&eos, Real cl) { } template -EOS makeSAPRamp(T &&eos, Real r0, Real a, Real b, Real c) { - return SAPRampEOS(std::forward(eos), r0, a, b, c); +EOS makeBilinearRamp(T &&eos, Real r0, Real a, Real b, Real c) { + return BilinearRampEOS(std::forward(eos), r0, a, b, c); } template @@ -126,11 +126,11 @@ EOS applyWrappedShiftAndScale(T &&eos, bool scaled, bool shifted, Real scale, Re } template -EOS applyShiftAndScaleAndSAPRamp(T &&eos, bool scaled, bool shifted, bool ramped, +EOS applyShiftAndScaleAndBilinearRamp(T &&eos, bool scaled, bool shifted, bool ramped, Real scale, Real shift, Real r0, Real a, Real b, Real c) { if (ramped) { - return applyWrappedShiftAndScale(std::forward(eos), scaled, shifted, + return applyWrappedShiftAndScale(std::forward(eos), scaled, shifted, scale, shift, r0, a, b, c); } else { return applyShiftAndScale(std::forward(eos), scaled, shifted, scale, shift); diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 5c74eb369e..e1fc476d2b 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -36,7 +36,7 @@ using namespace eos_base; using singularity::robust::ratio; template -void pAlpha2SAPRampParams(const T &eos, const Real alpha0, const Real Pe, const Real Pc, +void pAlpha2BilinearRampParams(const T &eos, const Real alpha0, const Real Pe, const Real Pc, Real &r0, Real &a, Real &b, Real &c) { // get reference conditions Real rho0, T0, sie0, P0, cv0, bmod0, dpde0, dvdt0; @@ -73,10 +73,10 @@ void pAlpha2SAPRampParams(const T &eos, const Real alpha0, const Real Pe, const } template -class SAPRampEOS : public EosBase> { +class BilinearRampEOS : public EosBase> { public: // move semantics ensures dynamic memory comes along for the ride - SAPRampEOS(T &&t, const Real r0, const Real a, const Real b, const Real c) + BilinearRampEOS(T &&t, const Real r0, const Real a, const Real b, const Real c) : t_(std::forward(t)), r0_(r0), a_(a), b_(b), c_(c), rmid_(r0 * (a - b * c) / (a - b)), Pmid_(a * (rmid_ / r0 - 1.0)) { // add input parameter checks to ensure validity of the ramp @@ -85,9 +85,9 @@ class SAPRampEOS : public EosBase> { assert(b >= 0); assert(a != b); } - SAPRampEOS() = default; + BilinearRampEOS() = default; - auto GetOnDevice() { return SAPRampEOS(t_.GetOnDevice(), r0_, a_, b_, c_); } + auto GetOnDevice() { return BilinearRampEOS(t_.GetOnDevice(), r0_, a_, b_, c_); } inline void Finalize() { t_.Finalize(); } PORTABLE_INLINE_FUNCTION @@ -230,7 +230,7 @@ class SAPRampEOS : public EosBase> { } // Vector functions that overload the scalar versions declared here. - SG_ADD_BASE_CLASS_USINGS(SAPRampEOS) + SG_ADD_BASE_CLASS_USINGS(BilinearRampEOS) private: T t_; diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index a19a42aad3..7f28ee7cef 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -41,7 +41,7 @@ int init_sg_eos(const int nmat, EOS *&eos) { EOSBuilder::applyShiftAndScale(A, enabled[0] == 1, enabled[1] == 1, vals[0], vals[1]) #define SGAPPLYMOD(A) \ - EOSBuilder::applyShiftAndScaleAndSAPRamp(A, enabled[0] == 1, enabled[1] == 1, \ + EOSBuilder::applyShiftAndScaleAndBilinearRamp(A, enabled[0] == 1, enabled[1] == 1, \ enabled[2] == 1 || enabled[3] == 1, vals[0], \ vals[1], vals[2], vals[3], vals[4], vals[5]) @@ -53,7 +53,7 @@ int init_sg_IdealGas(const int matindex, EOS *eos, const double gm1, const doubl assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(IdealGas(gm1, Cv)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(IdealGas(gm1, Cv)); @@ -72,7 +72,7 @@ int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const doubl assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(Gruneisen(C0, s1, s2, s3, G0, b, rho0, T0, P0, Cv)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(Gruneisen(C0, s1, s2, s3, G0, b, rho0, T0, P0, Cv)); @@ -94,7 +94,7 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(JWL(A, B, R1, R2, w, rho0, Cv)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(JWL(A, B, R1, R2, w, rho0, Cv)); @@ -115,7 +115,7 @@ int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const do assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); @@ -139,7 +139,7 @@ int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, EOS eosi = SGAPPLYMODSIMPLE(DavisReactants(rho0, e0, P0, T0, A, B, C, G0, Z, alpha, Cv0)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(DavisReactants(rho0, e0, P0, T0, A, B, C, G0, Z, alpha, Cv0)); @@ -163,7 +163,7 @@ int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(SpinerEOSDependsRhoT(std::string(filename), matid)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(SpinerEOSDependsRhoT(std::string(filename), matid)); @@ -182,7 +182,7 @@ int init_sg_SpinerDependsRhoSie(const int matindex, EOS *eos, const char *filena assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(SpinerEOSDependsRhoSie(std::string(filename), matid)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(SpinerEOSDependsRhoSie(std::string(filename), matid)); @@ -201,7 +201,7 @@ int init_sg_eospac(const int matindex, EOS *eos, const int id, int const *const assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(EOSPAC(id)); if (enabled[3] == 1) { - singularity::pAlpha2SAPRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(EOSPAC(id)); diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 6d9b931450..610d633a55 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -36,7 +36,7 @@ using singularity::EOS; using singularity::Gruneisen; using singularity::IdealGas; -using singularity::SAPRampEOS; +using singularity::BilinearRampEOS; using singularity::ScaledEOS; using singularity::ShiftedEOS; @@ -267,7 +267,7 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder],[Modifiers][IdealGas]") { ramp_params["a"].emplace(a); ramp_params["b"].emplace(b); ramp_params["c"].emplace(c); - modifiers[EOSBuilder::EOSModifier::SAPRamp] = ramp_params; + modifiers[EOSBuilder::EOSModifier::BilinearRamp] = ramp_params; THEN("The EOS is constructed correctly") { auto eos_ramped = EOSBuilder::buildEOS(type, base_params, modifiers); } From 2ae3f16fd5f49ac05c3ac512f312bef669fda189 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 5 Jul 2022 18:29:14 -0600 Subject: [PATCH 34/34] format --- singularity-eos/eos/eos_builder.cpp | 12 +++---- singularity-eos/eos/eos_builder.hpp | 8 ++--- singularity-eos/eos/modifiers/ramps_eos.hpp | 4 +-- singularity-eos/eos/singularity_eos.cpp | 38 ++++++++++----------- test/test_eos_unit.cpp | 2 +- 5 files changed, 32 insertions(+), 32 deletions(-) diff --git a/singularity-eos/eos/eos_builder.cpp b/singularity-eos/eos/eos_builder.cpp index 64ae241bda..7e64464eb9 100644 --- a/singularity-eos/eos/eos_builder.cpp +++ b/singularity-eos/eos/eos_builder.cpp @@ -105,7 +105,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par return makeRelativistic(std::move(g), cl); } return applyShiftAndScaleAndBilinearRamp(std::move(g), scaled, shifted, ramped, scale, - shift, r0, a, b, c); + shift, r0, a, b, c); } #ifdef SPINER_USE_HDF if (type == EOSType::SpinerEOSDependsRhoT || type == EOSType::SpinerEOSDependsRhoSie) { @@ -132,8 +132,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return applyShiftAndScaleAndBilinearRamp(std::move(s), scaled, shifted, ramped, scale, - shift, r0, a, b, c); + return applyShiftAndScaleAndBilinearRamp(std::move(s), scaled, shifted, ramped, + scale, shift, r0, a, b, c); } } else { string materialName = mpark::get(base_params["materialName"]); @@ -156,8 +156,8 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par if (relativistic) { return makeRelativistic(std::move(s), cl); } - return applyShiftAndScaleAndBilinearRamp(std::move(s), scaled, shifted, ramped, scale, - shift, r0, a, b, c); + return applyShiftAndScaleAndBilinearRamp(std::move(s), scaled, shifted, ramped, + scale, shift, r0, a, b, c); } } } @@ -180,7 +180,7 @@ EOS EOSBuilder::buildEOS(EOSBuilder::EOSType type, EOSBuilder::params_t base_par return makeRelativistic(std::move(s), cl); } return applyShiftAndScaleAndBilinearRamp(std::move(s), scaled, shifted, ramped, scale, - shift, r0, a, b, c); + shift, r0, a, b, c); } #endif if (type == EOSType::Gruneisen) { diff --git a/singularity-eos/eos/eos_builder.hpp b/singularity-eos/eos/eos_builder.hpp index 9e52d0be05..25f9a75e20 100644 --- a/singularity-eos/eos/eos_builder.hpp +++ b/singularity-eos/eos/eos_builder.hpp @@ -127,11 +127,11 @@ EOS applyWrappedShiftAndScale(T &&eos, bool scaled, bool shifted, Real scale, Re template EOS applyShiftAndScaleAndBilinearRamp(T &&eos, bool scaled, bool shifted, bool ramped, - Real scale, Real shift, Real r0, Real a, Real b, - Real c) { + Real scale, Real shift, Real r0, Real a, Real b, + Real c) { if (ramped) { - return applyWrappedShiftAndScale(std::forward(eos), scaled, shifted, - scale, shift, r0, a, b, c); + return applyWrappedShiftAndScale( + std::forward(eos), scaled, shifted, scale, shift, r0, a, b, c); } else { return applyShiftAndScale(std::forward(eos), scaled, shifted, scale, shift); } diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index e1fc476d2b..b972b6996c 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -36,8 +36,8 @@ using namespace eos_base; using singularity::robust::ratio; template -void pAlpha2BilinearRampParams(const T &eos, const Real alpha0, const Real Pe, const Real Pc, - Real &r0, Real &a, Real &b, Real &c) { +void pAlpha2BilinearRampParams(const T &eos, const Real alpha0, const Real Pe, + const Real Pc, Real &r0, Real &a, Real &b, Real &c) { // get reference conditions Real rho0, T0, sie0, P0, cv0, bmod0, dpde0, dvdt0; Real rmid, r1; diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index 7f28ee7cef..ea80b283c5 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -41,9 +41,9 @@ int init_sg_eos(const int nmat, EOS *&eos) { EOSBuilder::applyShiftAndScale(A, enabled[0] == 1, enabled[1] == 1, vals[0], vals[1]) #define SGAPPLYMOD(A) \ - EOSBuilder::applyShiftAndScaleAndBilinearRamp(A, enabled[0] == 1, enabled[1] == 1, \ - enabled[2] == 1 || enabled[3] == 1, vals[0], \ - vals[1], vals[2], vals[3], vals[4], vals[5]) + EOSBuilder::applyShiftAndScaleAndBilinearRamp( \ + A, enabled[0] == 1, enabled[1] == 1, enabled[2] == 1 || enabled[3] == 1, vals[0], \ + vals[1], vals[2], vals[3], vals[4], vals[5]) int def_en[4] = {0, 0, 0, 0}; double def_v[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; @@ -53,8 +53,8 @@ int init_sg_IdealGas(const int matindex, EOS *eos, const double gm1, const doubl assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(IdealGas(gm1, Cv)); if (enabled[3] == 1) { - singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], - vals[4], vals[5]); + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(IdealGas(gm1, Cv)); eos[matindex] = eos_.GetOnDevice(); @@ -72,8 +72,8 @@ int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const doubl assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(Gruneisen(C0, s1, s2, s3, G0, b, rho0, T0, P0, Cv)); if (enabled[3] == 1) { - singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], - vals[4], vals[5]); + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(Gruneisen(C0, s1, s2, s3, G0, b, rho0, T0, P0, Cv)); eos[matindex] = eos_.GetOnDevice(); @@ -94,8 +94,8 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(JWL(A, B, R1, R2, w, rho0, Cv)); if (enabled[3] == 1) { - singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], - vals[4], vals[5]); + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(JWL(A, B, R1, R2, w, rho0, Cv)); eos[matindex] = eos_.GetOnDevice(); @@ -115,8 +115,8 @@ int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const do assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); if (enabled[3] == 1) { - singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], - vals[4], vals[5]); + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); eos[matindex] = eos_.GetOnDevice(); @@ -139,8 +139,8 @@ int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, EOS eosi = SGAPPLYMODSIMPLE(DavisReactants(rho0, e0, P0, T0, A, B, C, G0, Z, alpha, Cv0)); if (enabled[3] == 1) { - singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], - vals[4], vals[5]); + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(DavisReactants(rho0, e0, P0, T0, A, B, C, G0, Z, alpha, Cv0)); eos[matindex] = eos_.GetOnDevice(); @@ -163,8 +163,8 @@ int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(SpinerEOSDependsRhoT(std::string(filename), matid)); if (enabled[3] == 1) { - singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], - vals[4], vals[5]); + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(SpinerEOSDependsRhoT(std::string(filename), matid)); eos[matindex] = eos_.GetOnDevice(); @@ -182,8 +182,8 @@ int init_sg_SpinerDependsRhoSie(const int matindex, EOS *eos, const char *filena assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(SpinerEOSDependsRhoSie(std::string(filename), matid)); if (enabled[3] == 1) { - singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], - vals[4], vals[5]); + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(SpinerEOSDependsRhoSie(std::string(filename), matid)); eos[matindex] = eos_.GetOnDevice(); @@ -201,8 +201,8 @@ int init_sg_eospac(const int matindex, EOS *eos, const int id, int const *const assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(EOSPAC(id)); if (enabled[3] == 1) { - singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], - vals[4], vals[5]); + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); } EOS eos_ = SGAPPLYMOD(EOSPAC(id)); eos[matindex] = eos_.GetOnDevice(); diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index 610d633a55..f7c5fec632 100644 --- a/test/test_eos_unit.cpp +++ b/test/test_eos_unit.cpp @@ -33,10 +33,10 @@ #include +using singularity::BilinearRampEOS; using singularity::EOS; using singularity::Gruneisen; using singularity::IdealGas; -using singularity::BilinearRampEOS; using singularity::ScaledEOS; using singularity::ShiftedEOS;