From 0e9147e7691bbe05081be2733a3b11cc7db7b786 Mon Sep 17 00:00:00 2001 From: Benjamin Joel Musick Date: Thu, 10 Aug 2023 14:53:04 -0600 Subject: [PATCH 001/405] initial carnahan starling commit --- singularity-eos/CMakeLists.txt | 1 + singularity-eos/eos/eos.hpp | 3 +- singularity-eos/eos/eos_carnahan_starling.hpp | 262 ++++++++ singularity-eos/eos/singularity_eos.cpp | 19 + singularity-eos/eos/singularity_eos.f90 | 27 + singularity-eos/eos/singularity_eos.hpp | 7 + test/CMakeLists.txt | 1 + test/test_eos_carnahan_starling.cpp | 589 ++++++++++++++++++ 8 files changed, 908 insertions(+), 1 deletion(-) create mode 100644 singularity-eos/eos/eos_carnahan_starling.hpp create mode 100644 test/test_eos_carnahan_starling.cpp diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 285545c10e..ba55737763 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -42,6 +42,7 @@ set(EOS_HEADERS eos/eos_eospac.hpp eos/eos_noble_abel.hpp eos/eos_stiff.hpp + eos/eos_carnahan_starling.hpp ) set(EOS_SRCS diff --git a/singularity-eos/eos/eos.hpp b/singularity-eos/eos/eos.hpp index c9893d1873..8f0540ab5c 100644 --- a/singularity-eos/eos/eos.hpp +++ b/singularity-eos/eos/eos.hpp @@ -33,6 +33,7 @@ #include // EOS models +#include #include #include #include @@ -67,7 +68,7 @@ using singularity::detail::transform_variadic_list; // all eos's static constexpr const auto full_eos_list = tl +#include +#include + +// Ports-of-call +#include + +// Base stuff +#include +#include +#include +#include + +namespace singularity { + +using namespace eos_base; + +class CarnahanStarling : public EosBase { + public: + CarnahanStarling() = default; + PORTABLE_INLINE_FUNCTION CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq) + : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(ROOM_TEMPERATURE), + _P0(ATMOSPHERIC_PRESSURE), _qp(0.0), + _rho0(DensityFromPressureTemperature(_P0, _T0)), _vol0(robust::ratio(1.0, _rho0)), + _sie0(Cv * _T0 + qq), + _bmod0(_rho0 * Cv * _T0 * + (PartialRhoZedFromDensity(_rho0) + + ZedFromDensity(_rho0) * ZedFromDensity(_rho0) * gm1)), + _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1) { + checkParams(); + } + PORTABLE_INLINE_FUNCTION CarnahanStarling(Real gm1, Real Cv, Real bb, Real qq, Real qp, + Real T0, Real P0) + : _Cv(Cv), _gm1(gm1), _bb(bb), _qq(qq), _T0(T0), _P0(P0), _qp(qp), + _rho0(DensityFromPressureTemperature(P0, T0)), _vol0(robust::ratio(1.0, _rho0)), + _sie0(Cv * T0 + qq), + _bmod0(_rho0 * Cv * T0 * + (PartialRhoZedFromDensity(_rho0) + + ZedFromDensity(_rho0) * ZedFromDensity(_rho0) * gm1)), + _dpde0(_rho0 * ZedFromDensity(_rho0) * gm1) { + checkParams(); + } + CarnahanStarling GetOnDevice() { return *this; } + PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return std::max(robust::SMALL(), (sie - _qq) / _Cv); + } + PORTABLE_INLINE_FUNCTION Real ZedFromDensity(const Real rho, + Real *lambda = nullptr) const { + const Real eta = _bb * rho; + const Real zed = robust::ratio(1.0 + eta + eta * eta - eta * eta * eta, + (1.0 - eta) * (1.0 - eta) * (1.0 - eta)); + return zed; + } + PORTABLE_INLINE_FUNCTION Real PartialRhoZedFromDensity(const Real rho, + Real *lambda = nullptr) const { + const Real eta = _bb * rho; + return robust::ratio(eta * eta * eta * eta - 4.0 * eta * eta * eta + 4.0 * eta * eta + + 4.0 * eta + 1.0, + (1.0 - eta) * (1.0 - eta) * (1.0 - eta) * (1.0 - eta)); + } + PORTABLE_INLINE_FUNCTION Real DensityFromPressureTemperature( + const Real press, const Real temperature, Real *lambda = nullptr) const { + Real real_root; + // iteration related variables + int newton_counter = 0; + const int newton_max = 50; + const Real rel_tol = 1.e-12; + real_root = 0.0; // initial guess + for (int i = 0; i < newton_max; i++) { + Real real_root_old = real_root; + Real newton_f = + _Cv * temperature * _gm1 * real_root_old * ZedFromDensity(real_root_old) - + press; + Real newton_f_prime = + _Cv * temperature * _gm1 * PartialRhoZedFromDensity(real_root_old); + real_root = real_root_old - robust::ratio(newton_f, newton_f_prime); + Real rel_err = std::abs(robust::ratio(real_root - real_root_old, real_root)); + newton_counter++; + if (rel_err < rel_tol) { + return std::max(robust::SMALL(), real_root); + } + } + if (newton_counter == newton_max) { + printf("*** (Warning) DensityFromPressureTemperature :: Convergence not met in " + "Carnahan-Starling EoS ***"); + } + return std::max(robust::SMALL(), real_root); + } + PORTABLE_INLINE_FUNCTION void checkParams() const { + PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive"); + PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive"); + PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be positive"); + } + PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + return std::max(_qq, _Cv * temperature + _qq); + } + PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + const Real zed = ZedFromDensity(rho); + return std::max(robust::SMALL(), zed * rho * temperature * _Cv * _gm1); + } + PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + const Real zed = ZedFromDensity(rho); + return std::max(robust::SMALL(), zed * rho * (sie - _qq) * _gm1); + } + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + const Real vol = robust::ratio(1.0, rho); + return _Cv * std::log(robust::ratio(temperature, _T0) + robust::SMALL()) + + _gm1 * _Cv * std::log(robust::ratio(vol, _vol0) + robust::SMALL()) - + _gm1 * _Cv * _bb * + (4.0 * (robust::ratio(1.0, vol - _bb) - robust::ratio(1.0, _vol0 - _bb)) + + _bb * (robust::ratio(1.0, (vol - _bb) * (vol - _bb)) - + robust::ratio(1.0, (_vol0 - _bb) * (_vol0 - _bb)))) + + _qp; + } + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + const Real vol = robust::ratio(1.0, rho); + return _Cv * std::log(robust::ratio(sie - _qq, _sie0 - _qq) + robust::SMALL()) + + _gm1 * _Cv * std::log(robust::ratio(vol, _vol0) + robust::SMALL()) - + _gm1 * _Cv * _bb * + (4.0 * (robust::ratio(1.0, vol - _bb) - robust::ratio(1.0, _vol0 - _bb)) + + _bb * (robust::ratio(1.0, (vol - _bb) * (vol - _bb)) - + robust::ratio(1.0, (_vol0 - _bb) * (_vol0 - _bb)))) + + _qp; + } + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + return _Cv; + } + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return _Cv; + } + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + return std::max(robust::SMALL(), + rho * _Cv * temperature * _gm1 * + (PartialRhoZedFromDensity(rho) + + ZedFromDensity(rho) * ZedFromDensity(rho) * _gm1)); + } + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return std::max(robust::SMALL(), + rho * (sie - _qq) * _gm1 * + (PartialRhoZedFromDensity(rho) + + ZedFromDensity(rho) * ZedFromDensity(rho) * _gm1)); + } + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + return ZedFromDensity(rho) * _gm1; + } + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return ZedFromDensity(rho) * _gm1; + } + PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, + Real &cv, Real &bmod, const unsigned long output, + Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION + void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, + Real &bmod, Real &dpde, Real &dvdt, + Real *lambda = nullptr) const { + // use STP: 1 atmosphere, room temperature + rho = _rho0; + temp = _T0; + sie = _sie0; + press = _P0; + cv = _Cv; + bmod = _bmod0; + dpde = _dpde0; + } + // Generic functions provided by the base class. These contain e.g. the vector + // overloads that use the scalar versions declared here + SG_ADD_BASE_CLASS_USINGS(CarnahanStarling) + PORTABLE_INLINE_FUNCTION + int nlambda() const noexcept { return 0; } + static constexpr unsigned long PreferredInput() { return _preferred_input; } + static inline unsigned long scratch_size(std::string method, unsigned int nelements) { + return 0; + } + static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } + PORTABLE_INLINE_FUNCTION void PrintParams() const { + printf("Carnahan-Starling Parameters:\nGamma = %g\nCv = %g\nb = %g\nq = " + "%g\n", + _gm1 + 1.0, _Cv, _bb, _qq); + } + PORTABLE_INLINE_FUNCTION void + DensityEnergyFromPressureTemperature(const Real press, const Real temp, Real *lambda, + Real &rho, Real &sie) const { + sie = std::max(_qq, _Cv * temp + _qq); + rho = DensityFromPressureTemperature(press, temp); + } + inline void Finalize() {} + static std::string EosType() { return std::string("CarnahanStarling"); } + static std::string EosPyType() { return EosType(); } + + private: + Real _Cv, _gm1, _bb, _qq; + // optional reference state variables + Real _T0, _P0, _qp; + // reference values + Real _rho0, _vol0, _sie0, _bmod0, _dpde0; + // static constexpr const Real _T0 = ROOM_TEMPERATURE; + // static constexpr const Real _P0 = ATMOSPHERIC_PRESSURE; + static constexpr const unsigned long _preferred_input = + thermalqs::density | thermalqs::specific_internal_energy; +}; + +PORTABLE_INLINE_FUNCTION +void CarnahanStarling::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, + Real &bmod, const unsigned long output, + Real *lambda) const { + if (output & thermalqs::density && output & thermalqs::specific_internal_energy) { + if (output & thermalqs::pressure || output & thermalqs::temperature) { + UNDEFINED_ERROR; + } + DensityEnergyFromPressureTemperature(press, temp, lambda, rho, sie); + } + if (output & thermalqs::pressure && output & thermalqs::specific_internal_energy) { + if (output & thermalqs::density || output & thermalqs::temperature) { + UNDEFINED_ERROR; + } + sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); + } + if (output & thermalqs::temperature && output & thermalqs::specific_internal_energy) { + sie = robust::ratio(press, ZedFromDensity(rho) * rho * _gm1) + _qq; + } + if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::temperature) + temp = TemperatureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::bulk_modulus) + bmod = BulkModulusFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::specific_heat) + cv = SpecificHeatFromDensityInternalEnergy(rho, sie); +} + +} // namespace singularity + +#endif // _SINGULARITY_EOS_EOS_EOS_CARNAHAN_STARLING_HPP_ diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index 54867736f0..540b79dee9 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -212,6 +212,25 @@ int init_sg_NobleAbel(const int matindex, EOS *eos, const double gm1, const doub return init_sg_NobleAbel(matindex, eos, gm1, Cv, bb, qq, def_en, def_v); } +int init_sg_CarnahanStarling(const int matindex, EOS *eos, const double gm1, + const double Cv, const double bb, const double qq, + int const *const enabled, double *const vals) { + assert(matindex >= 0); + EOS eosi = SGAPPLYMODSIMPLE(CarnahanStarling(gm1, Cv, bb, qq)); + if (enabled[3] == 1) { + singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], + vals[3], vals[4], vals[5]); + } + EOS eos_ = SGAPPLYMOD(CarnahanStarling(gm1, Cv, bb, qq)); + eos[matindex] = eos_.GetOnDevice(); + return 0; +} + +int init_sg_CarnahanStarling(const int matindex, EOS *eos, const double gm1, + const double Cv, const double bb, const double qq) { + return init_sg_CarnahanStarling(matindex, eos, gm1, Cv, bb, qq, def_en, def_v); +} + #ifdef SINGULARITY_USE_SPINER_WITH_HDF5 int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename, diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index 826fe8a61e..e0995ae19f 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -40,6 +40,7 @@ module singularity_eos init_sg_NobleAbel_f,& init_sg_SAP_Polynomial_f,& init_sg_StiffGas_f,& + init_sg_CarnahanStarling_f,& init_sg_SpinerDependsRhoT_f,& init_sg_SpinerDependsRhoSie_f,& init_sg_eospac_f,& @@ -136,6 +137,19 @@ end function init_sg_DavisReactants type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values end function init_sg_NobleAbel end interface + + interface + integer(kind=c_int) function & + init_sg_CarnahanStarling(matindex, eos, gm1, Cv, bb, qq, sg_mods_enabled, & + sg_mods_values) & + bind(C, name='init_sg_CarnahanStarling') + import + integer(c_int), value, intent(in) :: matindex + type(c_ptr), value, intent(in) :: eos + real(kind=c_double), value, intent(in) :: gm1, Cv, bb, qq + type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values + end function init_sg_CarnahanStarling + end interface interface integer(kind=c_int) function & @@ -428,6 +442,19 @@ integer function init_sg_NobleAbel_f(matindex, eos, gm1, Cv, & c_loc(sg_mods_enabled), c_loc(sg_mods_values)) end function init_sg_NobleAbel_f + integer function init_sg_CarnahanStarling_f(matindex, eos, gm1, Cv, & + bb, qq, & + sg_mods_enabled, sg_mods_values) & + result(err) + integer(c_int), value, intent(in) :: matindex + type(sg_eos_ary_t), intent(in) :: eos + real(kind=8), value, intent(in) :: gm1, Cv, bb, qq + integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + err = init_sg_CarnahanStarling(matindex-1, eos%ptr, gm1, Cv, bb, qq, & + c_loc(sg_mods_enabled), c_loc(sg_mods_values)) + end function init_sg_CarnahanStarling_f + integer function init_sg_SpinerDependsRhoT_f(matindex, eos, filename, id, & sg_mods_enabled, & sg_mods_values) & diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index 2382909666..ed46da4b62 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/singularity-eos/eos/singularity_eos.hpp @@ -64,6 +64,10 @@ int init_sg_NobleAbel(const int matindex, EOS *eos, const double gm1, const doub const double bb, const double qq, int const *const enabled, double *const vals); +int init_sg_CarnahanStarling(const int matindex, EOS *eos, const double gm1, + const double Cv, const double bb, const double qq, + int const *const enabled, double *const vals); + #ifdef SINGULARITY_USE_SPINER_WITH_HDF5 int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename, @@ -132,6 +136,9 @@ int init_sg_StiffGas(const int matindex, EOS *eos, const double gm1, const doubl int init_sg_NobleAbel(const int matindex, EOS *eos, const double gm1, const double Cv, const double bb, const double qq); +int init_sg_CarnahanStarling(const int matindex, EOS *eos, const double gm1, + const double Cv, const double bb, const double qq); + #ifdef SINGULARITY_USE_SPINER_WITH_HDF5 int init_sg_SpinerDependsRhoT(const int matindex, EOS *eos, const char *filename, diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c0fa96f8dd..48bfce94af 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -21,6 +21,7 @@ add_executable( test_eos_vinet.cpp test_eos_noble_abel.cpp test_eos_stiff.cpp + test_eos_carnahan_starling.cpp ) if(SINGULARITY_TEST_SESAME) diff --git a/test/test_eos_carnahan_starling.cpp b/test/test_eos_carnahan_starling.cpp new file mode 100644 index 0000000000..b0e096b30f --- /dev/null +++ b/test/test_eos_carnahan_starling.cpp @@ -0,0 +1,589 @@ +//------------------------------------------------------------------------------ +// © 2021-2023. 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. +//------------------------------------------------------------------------------ + +#include +#include +#include +#include +#include +#ifndef CATCH_CONFIG_RUNNER +#include "catch2/catch.hpp" +#endif + +#include +#include +#include + +using singularity::CarnahanStarling; +using singularity::EOS; + +SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") { + GIVEN("Parameters for a CarnahanStarling EOS") { + constexpr Real gm1 = 0.4; + constexpr Real Cv = 7180000.0; + constexpr Real bb = 1.e-3; + constexpr Real qq = 0.0; + // Create the EOS + EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq); + EOS eos = host_eos.GetOnDevice(); + GIVEN("Densities and energies") { + constexpr int num = 4; +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_density("density"); + Kokkos::View v_energy("density"); +#else + // Otherwise just create arrays to contain values and create pointers to + // be passed to the functions in place of the Kokkos views + std::array density; + std::array energy; + auto v_density = density.data(); + auto v_energy = energy.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize density and energy", 0, 1, PORTABLE_LAMBDA(int i) { + v_density[0] = 1.1833012071291069e-03; + v_density[1] = 7.8290736890381501e-03; + v_density[2] = 8.5453943327882340e-03; + v_density[3] = 8.8197619601121831e-03; + v_energy[0] = 2.1407169999999998e+09; + v_energy[1] = 1.1000478000000000e+10; + v_energy[2] = 1.9860239000000000e+10; + v_energy[3] = 2.8720000000000000e+10; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create host-side mirrors of the inputs and copy the inputs. These are + // just used for the comparisons + auto density = Kokkos::create_mirror_view(v_density); + auto energy = Kokkos::create_mirror_view(v_energy); + Kokkos::deep_copy(density, v_density); + Kokkos::deep_copy(energy, v_energy); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Gold standard values for a subset of lookups + constexpr std::array pressure_true{ + 1.0132499999999999e+06, 3.4450500000000000e+07, 6.7887750000000000e+07, + 1.0132500000000000e+08}; + constexpr std::array bulkmodulus_true{ + 1.4185567142990597e+06, 4.8232210423710555e+07, 9.5046098754186615e+07, + 1.4186000457238689e+08}; + constexpr std::array temperature_true{ + 2.9814999999999998e+02, 1.5320999999999999e+03, 2.7660500000000002e+03, + 4.0000000000000000e+03}; + constexpr std::array gruneisen_true{ + 4.0000189328753222e-01, 4.0001252676308346e-01, 4.0001367292303214e-01, + 4.0001411193029396e-01}; + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create device views for outputs and mirror those views on the host + Kokkos::View v_temperature("Temperature"); + Kokkos::View v_pressure("Pressure"); + Kokkos::View v_bulkmodulus("bmod"); + Kokkos::View v_gruneisen("Gamma"); + auto h_temperature = Kokkos::create_mirror_view(v_temperature); + auto h_pressure = Kokkos::create_mirror_view(v_pressure); + auto h_bulkmodulus = Kokkos::create_mirror_view(v_bulkmodulus); + auto h_gruneisen = Kokkos::create_mirror_view(v_gruneisen); +#else + // Create arrays for the outputs and then pointers to those arrays that + // will be passed to the functions in place of the Kokkos views + std::array h_temperature; + std::array h_pressure; + std::array h_bulkmodulus; + std::array h_gruneisen; + // Just alias the existing pointers + auto v_temperature = h_temperature.data(); + auto v_pressure = h_pressure.data(); + auto v_bulkmodulus = h_bulkmodulus.data(); + auto v_gruneisen = h_gruneisen.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + WHEN("A P(rho, e) lookup is performed") { + eos.PressureFromDensityInternalEnergy(v_density, v_energy, v_pressure, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_pressure, v_pressure); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned P(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_pressure, pressure_true, "Density", + "Energy"); + } + } + + WHEN("A B_S(rho, e) lookup is performed") { + eos.BulkModulusFromDensityInternalEnergy(v_density, v_energy, v_bulkmodulus, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_bulkmodulus, v_bulkmodulus); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned B_S(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_bulkmodulus, bulkmodulus_true, "Density", + "Energy"); + } + } + + WHEN("A T(rho, e) lookup is performed") { + eos.TemperatureFromDensityInternalEnergy(v_density, v_energy, v_temperature, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_temperature, v_temperature); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned B_S(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_temperature, temperature_true, "Density", + "Energy"); + } + } + + WHEN("A Gamma(rho, e) lookup is performed") { + eos.GruneisenParamFromDensityInternalEnergy(v_density, v_energy, v_gruneisen, + num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_gruneisen, v_gruneisen); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned Gamma(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_gruneisen, gruneisen_true, "Density", + "Energy"); + } + } + } + } +} + +SCENARIO("CarnahanStarling2", "[CarnahanStarling][CarnahanStarling2]") { + GIVEN("Parameters for a CarnahanStarling EOS") { + constexpr Real gm1 = 0.4; + constexpr Real Cv = 7180000.0; + constexpr Real bb = 1.e-3; + constexpr Real qq = 42.00e+09; + constexpr Real qp = -23.0e+7; + constexpr Real T0 = 200.0; + constexpr Real P0 = 1000000.0; + // Create the EOS + EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq, qp, T0, P0); + EOS eos = host_eos.GetOnDevice(); + GIVEN("Densities and energies") { + constexpr int num = 4; +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_density("density"); + Kokkos::View v_energy("density"); +#else + // Otherwise just create arrays to contain values and create pointers to + // be passed to the functions in place of the Kokkos views + std::array density; + std::array energy; + auto v_density = density.data(); + auto v_energy = energy.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize density and energy", 0, 1, PORTABLE_LAMBDA(int i) { + v_density[0] = 1.1833012071291069e-03; + v_density[1] = 7.8290736890381501e-03; + v_density[2] = 8.5453943327882340e-03; + v_density[3] = 8.8197619601121831e-03; + v_energy[0] = 4.4140717000000000e+10; + v_energy[1] = 5.3000478000000000e+10; + v_energy[2] = 6.1860239000000000e+10; + v_energy[3] = 7.0720000000000000e+10; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create host-side mirrors of the inputs and copy the inputs. These are + // just used for the comparisons + auto density = Kokkos::create_mirror_view(v_density); + auto energy = Kokkos::create_mirror_view(v_energy); + Kokkos::deep_copy(density, v_density); + Kokkos::deep_copy(energy, v_energy); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Gold standard values for a subset of lookups + constexpr std::array pressure_true{ + 1.0132499999999999e+06, 3.4450500000000000e+07, 6.7887750000000000e+07, + 1.0132500000000000e+08}; + constexpr std::array bulkmodulus_true{ + 1.4185567142990597e+06, 4.8232210423710555e+07, 9.5046098754186615e+07, + 1.4186000457238689e+08}; + constexpr std::array temperature_true{ + 2.9814999999999998e+02, 1.5320999999999999e+03, 2.7660500000000002e+03, + 4.0000000000000000e+03}; + constexpr std::array gruneisen_true{ + 4.0000189328753222e-01, 4.0001252676308346e-01, 4.0001367292303214e-01, + 4.0001411193029396e-01}; + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create device views for outputs and mirror those views on the host + Kokkos::View v_temperature("Temperature"); + Kokkos::View v_pressure("Pressure"); + Kokkos::View v_bulkmodulus("bmod"); + Kokkos::View v_gruneisen("Gamma"); + auto h_temperature = Kokkos::create_mirror_view(v_temperature); + auto h_pressure = Kokkos::create_mirror_view(v_pressure); + auto h_bulkmodulus = Kokkos::create_mirror_view(v_bulkmodulus); + auto h_gruneisen = Kokkos::create_mirror_view(v_gruneisen); +#else + // Create arrays for the outputs and then pointers to those arrays that + // will be passed to the functions in place of the Kokkos views + std::array h_temperature; + std::array h_pressure; + std::array h_bulkmodulus; + std::array h_gruneisen; + // Just alias the existing pointers + auto v_temperature = h_temperature.data(); + auto v_pressure = h_pressure.data(); + auto v_bulkmodulus = h_bulkmodulus.data(); + auto v_gruneisen = h_gruneisen.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + WHEN("A P(rho, e) lookup is performed") { + eos.PressureFromDensityInternalEnergy(v_density, v_energy, v_pressure, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_pressure, v_pressure); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned P(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_pressure, pressure_true, "Density", + "Energy"); + } + } + + WHEN("A B_S(rho, e) lookup is performed") { + eos.BulkModulusFromDensityInternalEnergy(v_density, v_energy, v_bulkmodulus, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_bulkmodulus, v_bulkmodulus); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned B_S(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_bulkmodulus, bulkmodulus_true, "Density", + "Energy"); + } + } + + WHEN("A T(rho, e) lookup is performed") { + eos.TemperatureFromDensityInternalEnergy(v_density, v_energy, v_temperature, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_temperature, v_temperature); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned B_S(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_temperature, temperature_true, "Density", + "Energy"); + } + } + + WHEN("A Gamma(rho, e) lookup is performed") { + eos.GruneisenParamFromDensityInternalEnergy(v_density, v_energy, v_gruneisen, + num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_gruneisen, v_gruneisen); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned Gamma(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_gruneisen, gruneisen_true, "Density", + "Energy"); + } + } + } + } +} + +SCENARIO("(C-S EoS) Isentropic Bulk Modulus Analytic vs. FD", + "[CarnahanStarling][CarnahanStarling3]") { + GIVEN("Parameters for a C-S Gas EOS") { + constexpr Real gm1 = 0.4; + constexpr Real Cv = 7180000.0; + constexpr Real bb = 1.e-3; + constexpr Real qq = 0.0; + // Create the EOS + EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq); + EOS eos = host_eos.GetOnDevice(); + GIVEN("Density and energy") { + constexpr Real density = 1.1833012071291069e-03; + constexpr Real energy = 2.1407169999999998e+09; + constexpr Real true_sound_speed = 3.4623877142797290e+04; + WHEN("A B_S(rho, e) lookup is performed") { + const Real bulk_modulus = + eos.BulkModulusFromDensityInternalEnergy(density, energy); + THEN("The correct sound speed should be computed") { + const Real sound_speed = std::sqrt(bulk_modulus / density); + INFO("Density: " << density << " Energy: " << energy + << " Sound speed: " << sound_speed + << " cm/s True sound speed: " << true_sound_speed << " cm/s"); + REQUIRE(isClose(sound_speed, true_sound_speed, 1e-12)); + } + } + WHEN("A finite difference approximation is used for the bulk modulus") { + // Bulk modulus approximation: + // B_S = rho * dPdr_e + P / rho * dPde_r + constexpr Real drho = 1e-06 * density; + constexpr Real de = 1e-06 * energy; + const Real P1 = eos.PressureFromDensityInternalEnergy(density, energy); + Real P2 = eos.PressureFromDensityInternalEnergy(density + drho, energy); + const Real dPdr_e = (P2 - P1) / drho; + P2 = eos.PressureFromDensityInternalEnergy(density, energy + de); + const Real dPde_r = (P2 - P1) / de; + const Real bmod_approx = density * dPdr_e + P1 / density * dPde_r; + THEN("The finite difference solution should approximate the exact solution") { + const Real bulk_modulus = + eos.BulkModulusFromDensityInternalEnergy(density, energy); + const Real ss_approx = std::sqrt(bmod_approx / density); + const Real sound_speed = std::sqrt(bulk_modulus / density); + INFO("Density: " << density << " Energy: " << energy + << " Sound speed: " << sound_speed + << " cm/s Approximate sound speed: " << ss_approx << " cm/s"); + REQUIRE(isClose(sound_speed, ss_approx, 1e-5)); + } + } + } + } +} + +SCENARIO("Recover Ideal Gas from C-S", "[CarnahanStarling][CarnahanStarling4]") { + GIVEN("Parameters for a CarnahanStarling EOS") { + constexpr Real gm1 = 0.4; + constexpr Real Cv = 7180000.0; + constexpr Real bb = 0; + constexpr Real qq = 0; + // Create the EOS + EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq); + EOS eos = host_eos.GetOnDevice(); + EOS ideal_eos = singularity::IdealGas(gm1, Cv); + GIVEN("Densities and energies") { + constexpr int num = 1; +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_density("density"); + Kokkos::View v_energy("density"); +#else + // Otherwise just create arrays to contain values and create pointers to + // be passed to the functions in place of the Kokkos views + std::array density; + std::array energy; + auto v_density = density.data(); + auto v_energy = energy.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize density and energy", 0, 1, PORTABLE_LAMBDA(int i) { + v_density[0] = 1.1833068079526625e-03; + v_energy[0] = 2.1407169999999998e+09; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create host-side mirrors of the inputs and copy the inputs. These are + // just used for the comparisons + auto density = Kokkos::create_mirror_view(v_density); + auto energy = Kokkos::create_mirror_view(v_energy); + Kokkos::deep_copy(density, v_density); + Kokkos::deep_copy(energy, v_energy); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // values for a subset of lookups + std::array pressure_true; + std::array bulkmodulus_true; + std::array temperature_true; + std::array gruneisen_true; + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create device views for outputs and mirror those views on the host + Kokkos::View v_temperature("Temperature"); + Kokkos::View v_pressure("Pressure"); + Kokkos::View v_bulkmodulus("bmod"); + Kokkos::View v_gruneisen("Gamma"); + auto h_temperature = Kokkos::create_mirror_view(v_temperature); + auto h_pressure = Kokkos::create_mirror_view(v_pressure); + auto h_bulkmodulus = Kokkos::create_mirror_view(v_bulkmodulus); + auto h_gruneisen = Kokkos::create_mirror_view(v_gruneisen); +#else + // Create arrays for the outputs and then pointers to those arrays that + // will be passed to the functions in place of the Kokkos views + std::array h_temperature; + std::array h_pressure; + std::array h_bulkmodulus; + std::array h_gruneisen; + // Just alias the existing pointers + auto v_temperature = h_temperature.data(); + auto v_pressure = h_pressure.data(); + auto v_bulkmodulus = h_bulkmodulus.data(); + auto v_gruneisen = h_gruneisen.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + WHEN("A P(rho, e) lookup is performed") { + eos.PressureFromDensityInternalEnergy(v_density, v_energy, v_pressure, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_pressure, v_pressure); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned P(rho, e) should be equal to the true value") { + for (int i = 0; i < num; i++) { + pressure_true[i] = + ideal_eos.PressureFromDensityInternalEnergy(density[i], energy[i]); + } + array_compare(num, density, energy, h_pressure, pressure_true, "Density", + "Energy"); + } + } + + WHEN("A B_S(rho, e) lookup is performed") { + eos.BulkModulusFromDensityInternalEnergy(v_density, v_energy, v_bulkmodulus, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_bulkmodulus, v_bulkmodulus); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned B_S(rho, e) should be equal to the true value") { + for (int i = 0; i < num; i++) { + bulkmodulus_true[i] = + ideal_eos.BulkModulusFromDensityInternalEnergy(density[i], energy[i]); + } + array_compare(num, density, energy, h_bulkmodulus, bulkmodulus_true, "Density", + "Energy"); + } + } + + WHEN("A T(rho, e) lookup is performed") { + eos.TemperatureFromDensityInternalEnergy(v_density, v_energy, v_temperature, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_temperature, v_temperature); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned T(rho, e) should be equal to the true value") { + for (int i = 0; i < num; i++) { + temperature_true[i] = + ideal_eos.TemperatureFromDensityInternalEnergy(density[i], energy[i]); + } + array_compare(num, density, energy, h_temperature, temperature_true, "Density", + "Energy"); + } + } + + WHEN("A Gamma(rho, e) lookup is performed") { + eos.GruneisenParamFromDensityInternalEnergy(v_density, v_energy, v_gruneisen, + num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_gruneisen, v_gruneisen); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned Gamma(rho, e) should be equal to the true value") { + for (int i = 0; i < num; i++) { + gruneisen_true[i] = + ideal_eos.GruneisenParamFromDensityInternalEnergy(density[i], energy[i]); + } + array_compare(num, density, energy, h_gruneisen, gruneisen_true, "Density", + "Energy"); + } + } + } + } +} + +SCENARIO("Test C-S Entropy Calls", "[CarnahanStarling][CarnahanStarling5]") { + GIVEN("Parameters for a CarnahanStarling EOS") { + constexpr Real gm1 = 0.4; + constexpr Real Cv = 7180000.0; + constexpr Real bb = 1.e-3; + constexpr Real qq = 42.0e+9; + constexpr Real qp = -23.0e+9; + constexpr Real T0 = 200.0; + constexpr Real P0 = 1000000.0; + // Create the EOS + EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq, qp, T0, P0); + EOS eos = host_eos.GetOnDevice(); + GIVEN("Densities and energies") { + constexpr int num = 1; +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_density("density"); + Kokkos::View v_energy("density"); +#else + // Otherwise just create arrays to contain values and create pointers to + // be passed to the functions in place of the Kokkos views + std::array density; + std::array energy; + auto v_density = density.data(); + auto v_energy = energy.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize density and energy", 0, 1, PORTABLE_LAMBDA(int i) { + v_density[0] = 8.8197619601121831e-03; + v_energy[0] = 7.0720000000000000e+10; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create host-side mirrors of the inputs and copy the inputs. These are + // just used for the comparisons + auto density = Kokkos::create_mirror_view(v_density); + auto energy = Kokkos::create_mirror_view(v_energy); + Kokkos::deep_copy(density, v_density); + Kokkos::deep_copy(energy, v_energy); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Gold standard values for a subset of lookups + constexpr std::array pressure_true{1.0132500000000000e+08}; + constexpr std::array temperature_true{4.0000000000000000e+03}; + constexpr std::array entropy_true{-2.2983150752058342e+10}; + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create device views for outputs and mirror those views on the host + Kokkos::View v_temperature("Temperature"); + Kokkos::View v_pressure("Pressure"); + Kokkos::View v_entropy("entr"); + Kokkos::View v_local_temp("temp"); + auto h_temperature = Kokkos::create_mirror_view(v_temperature); + auto h_pressure = Kokkos::create_mirror_view(v_pressure); + auto h_entropy = Kokkos::create_mirror_view(v_entropy); + auto h_local_temp = Kokkos::create_mirror_view(v_local_temp); +#else + // Create arrays for the outputs and then pointers to those arrays that + // will be passed to the functions in place of the Kokkos views + std::array h_temperature; + std::array h_pressure; + std::array h_entropy; + std::array h_local_temp; + // Just alias the existing pointers + auto v_temperature = h_temperature.data(); + auto v_pressure = h_pressure.data(); + auto v_entropy = h_entropy.data(); + auto v_local_temp = h_local_temp.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + WHEN("A S(rho, e) lookup is performed") { + eos.EntropyFromDensityInternalEnergy(v_density, v_energy, v_entropy, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_entropy, v_entropy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned S(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_entropy, entropy_true, "Density", + "Energy"); + } + } + WHEN("A S(rho, T(rho,e)) lookup is performed") { + eos.TemperatureFromDensityInternalEnergy(v_density, v_energy, v_local_temp, num); + eos.EntropyFromDensityTemperature(v_density, v_local_temp, v_entropy, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_entropy, v_entropy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned S(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_entropy, entropy_true, "Density", + "Energy"); + } + } + } + } +} From d0af2babd4ebc5af1cb1e483db2481ddcff4749c Mon Sep 17 00:00:00 2001 From: Benjamin Joel Musick Date: Thu, 10 Aug 2023 15:40:14 -0600 Subject: [PATCH 002/405] added doc --- doc/sphinx/src/models.rst | 71 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index a6fc12c49d..c73f501185 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -542,6 +542,77 @@ these values are not set, they will be the same as those returned by the conditions are given, the return values of the :code:`ValuesAtReferenceState()` function will not be the same. +Carnahan-Starling +````````````````` + +The (quasi-exact) Carnahan-Starling model in ``singularity-eos`` takes +the form + +.. math:: + + P = Z(\rho) \rho (e-q) (\gamma-1) + +.. math:: + + Z(\rho) = \frac{1+\eta+\eta^2-\eta^3}{(1-\eta)^3} + +.. math:: + + \eta = b\rho + +.. math:: + + e = C_V T + q, + +the inputs for the equation of state are similar to the Noble-Abel model. +As with the Noble-Abel EOS, it should be noted that covolume is physically +significant as it represents the maximum compressibility of the gas, +and as a result it should be non-negative. + +The entropy for the Carnahan-Starling EOS is given by + +.. math:: + + S = C_V \ln\left(\frac{T}{T_0}\right) + C_V (\gamma-1) \left\{ \ln\left(\frac{v} + {v_0}\right) - S^{CS} \right\} + q', + +.. math:: + S^{CS} = b\left(4\left(\frac{1}{v-b} - \frac{1}{v_0-b}\right)+ + b\left(\frac{1}{(v-b)^2} - \frac{1}{(v_0-b)^2}\right)\right) + +where :math:`S(\rho_0,T_0)=q'`. By default, :math:`T_0 = 298` K and the +reference density is given by + +.. math:: + + P_0 = \rho_0 Z(\rho_0) C_V T_0(\gamma-1), + +where :math:`P_0` is by default 1 bar. Denisty is obtained via Newton iteration. + +The settable parameters for this EOS are :math:`\gamma-1`, specific +heat capacity (:math:`C_V`), covolume (:math:`b`) and offset internal energy (:math:`q`). Optionally, the reference state for the entropy calculation can +be provided by setting the reference temperature, pressure, and entropy offset. + +The ``CarnahanStarling`` EOS constructor has four arguments: ``gm1``, which is :math:`\gamma-1`; ``Cv``, the +specific heat :math:`C_V`; :math:`b`, the covolume; and :math:`q`, the internal energy offset. + +.. code-block:: cpp + + CarnahanStarling(Real gm1, Real Cv, Real b, Real q) + +Optionally, the reference state for the entropy calculation, +can be provided in the constructor via ``qp``, ``T0`` and ``P0``: + +.. code-block:: cpp + + CarnahanStarling(Real gm1, Real Cv, Real b, Real q, Real qp, Real T0, Real P0) + +Note that these parameters are provided solely for the entropy calculation. When +these values are not set, they will be the same as those returned by the +:code:`ValuesAtReferenceState()` function. However, if the entropy reference +conditions are given, the return values of the :code:`ValuesAtReferenceState()` +function will not be the same. + Gruneisen EOS ````````````` From 9703866a2c6ae8365144f6bcc201c0c5871824b0 Mon Sep 17 00:00:00 2001 From: Benjamin Joel Musick Date: Thu, 10 Aug 2023 15:45:26 -0600 Subject: [PATCH 003/405] updated changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index e70c1f8242..8586e8f6bb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ - [[PR269]](https://github.com/lanl/singularity-eos/pull/269) Add SAP Polynomial EoS ### Fixed (Repair bugs, etc) +- [[PR292]](https://github.com/lanl/singularity-eos/pull/292) Added Carnahan-Starling EoS - [[PR290]](https://github.com/lanl/singularity-eos/pull/290) Added target guards on export config - [[PR288]](https://github.com/lanl/singularity-eos/pull/288) Don't build tests that depend on spiner when spiner is disabled - [[PR287]](https://github.com/lanl/singularity-eos/pull/287) Fix testing logic with new HDF5 options From 32f461a0604256ebf88f857c1215aea5bf586fbd Mon Sep 17 00:00:00 2001 From: Benjamin Joel Musick Date: Mon, 14 Aug 2023 17:59:50 -0600 Subject: [PATCH 004/405] revised eos header --- singularity-eos/eos/eos_carnahan_starling.hpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/singularity-eos/eos/eos_carnahan_starling.hpp b/singularity-eos/eos/eos_carnahan_starling.hpp index dc14c0f89e..64fd81e9f9 100644 --- a/singularity-eos/eos/eos_carnahan_starling.hpp +++ b/singularity-eos/eos/eos_carnahan_starling.hpp @@ -66,25 +66,27 @@ class CarnahanStarling : public EosBase { PORTABLE_INLINE_FUNCTION Real ZedFromDensity(const Real rho, Real *lambda = nullptr) const { const Real eta = _bb * rho; - const Real zed = robust::ratio(1.0 + eta + eta * eta - eta * eta * eta, + const Real eta2 = eta * eta; + const Real zed = robust::ratio(1.0 + eta + eta2 - eta2 * eta, (1.0 - eta) * (1.0 - eta) * (1.0 - eta)); return zed; } PORTABLE_INLINE_FUNCTION Real PartialRhoZedFromDensity(const Real rho, Real *lambda = nullptr) const { const Real eta = _bb * rho; - return robust::ratio(eta * eta * eta * eta - 4.0 * eta * eta * eta + 4.0 * eta * eta + - 4.0 * eta + 1.0, + const Real eta2 = eta * eta; + return robust::ratio(eta2 * eta2 - 4.0 * eta2 * eta + 4.0 * eta2 + 4.0 * eta + 1.0, (1.0 - eta) * (1.0 - eta) * (1.0 - eta) * (1.0 - eta)); } PORTABLE_INLINE_FUNCTION Real DensityFromPressureTemperature( - const Real press, const Real temperature, Real *lambda = nullptr) const { + const Real press, const Real temperature, const Real guess = robust::SMALL(), + Real *lambda = nullptr) const { Real real_root; // iteration related variables int newton_counter = 0; const int newton_max = 50; const Real rel_tol = 1.e-12; - real_root = 0.0; // initial guess + real_root = guess; for (int i = 0; i < newton_max; i++) { Real real_root_old = real_root; Real newton_f = From 12ff7cf5fed03e3a88ec9261b308a0ef6c53a0d7 Mon Sep 17 00:00:00 2001 From: Benjamin Joel Musick Date: Wed, 16 Aug 2023 10:02:24 -0600 Subject: [PATCH 005/405] uses root finding utils, adds test to check root finding --- singularity-eos/eos/eos_carnahan_starling.hpp | 42 ++++++--- test/test_eos_carnahan_starling.cpp | 90 +++++++++++++++++++ 2 files changed, 122 insertions(+), 10 deletions(-) diff --git a/singularity-eos/eos/eos_carnahan_starling.hpp b/singularity-eos/eos/eos_carnahan_starling.hpp index 64fd81e9f9..ae8c7107b8 100644 --- a/singularity-eos/eos/eos_carnahan_starling.hpp +++ b/singularity-eos/eos/eos_carnahan_starling.hpp @@ -27,6 +27,7 @@ #include #include #include +#include #include namespace singularity { @@ -63,6 +64,11 @@ class CarnahanStarling : public EosBase { const Real rho, const Real sie, Real *lambda = nullptr) const { return std::max(robust::SMALL(), (sie - _qq) / _Cv); } + PORTABLE_INLINE_FUNCTION void checkParams() const { + PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive"); + PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive"); + PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be positive"); + } PORTABLE_INLINE_FUNCTION Real ZedFromDensity(const Real rho, Real *lambda = nullptr) const { const Real eta = _bb * rho; @@ -82,10 +88,29 @@ class CarnahanStarling : public EosBase { const Real press, const Real temperature, const Real guess = robust::SMALL(), Real *lambda = nullptr) const { Real real_root; - // iteration related variables +#ifdef _SINGULARITY_EOS_UTILS_ROOT_FINDING_HPP_ + // use singularity utilities if available + auto poly = [=](Real dens) { + return _Cv * temperature * _gm1 * dens * ZedFromDensity(dens); + }; + using RootFinding1D::findRoot; + using RootFinding1D::Status; + static constexpr Real xtol = 1.0e-12; + static constexpr Real ytol = 1.0e-12; + static constexpr Real lo_bound = robust::SMALL(); + const Real hi_bound = robust::ratio(1.0, _bb); + auto status = findRoot(poly, press, guess, lo_bound, hi_bound, xtol, ytol, real_root); + if (status != Status::SUCCESS) { + // Root finder failed even though the solution was bracketed... this is an error + EOS_ERROR("*** (Warning) DensityFromPressureTemperature :: Convergence not met in " + "Carnahan-Starling EoS (regula_falsi) ***\n"); + real_root = -1.0; + } +#else + // newton-raphson iteration if utilities not available int newton_counter = 0; - const int newton_max = 50; - const Real rel_tol = 1.e-12; + static constexpr int newton_max = 50; + static constexpr Real rel_tol = 1.e-12; real_root = guess; for (int i = 0; i < newton_max; i++) { Real real_root_old = real_root; @@ -102,16 +127,13 @@ class CarnahanStarling : public EosBase { } } if (newton_counter == newton_max) { - printf("*** (Warning) DensityFromPressureTemperature :: Convergence not met in " - "Carnahan-Starling EoS ***"); + EOS_ERROR("*** (Warning) DensityFromPressureTemperature :: Convergence not met in " + "Carnahan-Starling EoS (newton-raphson) ***\n"); + real_root = -1.0; } +#endif return std::max(robust::SMALL(), real_root); } - PORTABLE_INLINE_FUNCTION void checkParams() const { - PORTABLE_ALWAYS_REQUIRE(_Cv >= 0, "Heat capacity must be positive"); - PORTABLE_ALWAYS_REQUIRE(_gm1 >= 0, "Gruneisen parameter must be positive"); - PORTABLE_ALWAYS_REQUIRE(_bb >= 0, "Covolume must be positive"); - } PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { return std::max(_qq, _Cv * temperature + _qq); diff --git a/test/test_eos_carnahan_starling.cpp b/test/test_eos_carnahan_starling.cpp index b0e096b30f..4f0c6da5c8 100644 --- a/test/test_eos_carnahan_starling.cpp +++ b/test/test_eos_carnahan_starling.cpp @@ -587,3 +587,93 @@ SCENARIO("Test C-S Entropy Calls", "[CarnahanStarling][CarnahanStarling5]") { } } } + +SCENARIO("CarnahanStarling6", "[CarnahanStarling][CarnahanStarling6]") { + GIVEN("Parameters for a CarnahanStarling EOS") { + constexpr Real gm1 = 0.4; + constexpr Real Cv = 7180000.0; + constexpr Real bb = 1.e-3; + constexpr Real qq = 0.0; + // Create the EOS + EOS host_eos = CarnahanStarling(gm1, Cv, bb, qq); + EOS eos = host_eos.GetOnDevice(); + GIVEN("Pressure and temperature") { + constexpr int num = 4; +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_pressure("Pressure"); + Kokkos::View v_temperature("Temperature"); +#else + // Otherwise just create arrays to contain values and create pointers to + // be passed to the functions in place of the Kokkos views + std::array pressure; + std::array temperature; + auto v_pressure = pressure.data(); + auto v_temperature = temperature.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize pressure and temperature", 0, 1, PORTABLE_LAMBDA(int i) { + v_pressure[0] = 1.0132499999999999e+06; + v_pressure[1] = 3.4450500000000000e+07; + v_pressure[2] = 6.7887750000000000e+07; + v_pressure[3] = 1.0132500000000000e+08; + v_temperature[0] = 2.9814999999999998e+02; + v_temperature[1] = 1.5320999999999999e+03; + v_temperature[2] = 2.7660500000000002e+03; + v_temperature[3] = 4.0000000000000000e+03; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create host-side mirrors of the inputs and copy the inputs. These are + // just used for the comparisons + auto pressure = Kokkos::create_mirror_view(v_pressure); + auto temperature = Kokkos::create_mirror_view(v_temperature); + Kokkos::deep_copy(density, v_pressure); + Kokkos::deep_copy(energy, v_temperature); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Gold standard values for a subset of lookups + constexpr std::array density_true{ + 1.1833012071291069e-03, 7.8290736890381501e-03, 8.5453943327882340e-03, + 8.8197619601121831e-03}; + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create device views for outputs and mirror those views on the host + Kokkos::View v_density("Density"); + Kokkos::View v_energy("Energy"); + auto h_density = Kokkos::create_mirror_view(v_density); + auto h_energy = Kokkos::create_mirror_view(v_energy); +#else + // Create arrays for the outputs and then pointers to those arrays that + // will be passed to the functions in place of the Kokkos views + std::array h_density; + std::array h_energy; + // Just alias the existing pointers + auto v_density = h_density.data(); + auto v_energy = h_energy.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + WHEN("A rho(P, T) lookup is performed") { + for (int i = 0; i < num; i++) { + Real cv, bmod; + static constexpr const unsigned long _output = + singularity::thermalqs::density | + singularity::thermalqs::specific_internal_energy; + eos.FillEos(v_density[i], v_temperature[i], v_energy[i], v_pressure[i], cv, + bmod, _output); + } + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_density, v_density); +#endif // PORTABILITY_STRATEGY_KOKKOS + + THEN("The returned rho(P, T) should be equal to the true value") { + array_compare(num, pressure, temperature, h_density, density_true, "Pressure", + "Temperature"); + } + } + } + } +} From 0648d59d0222d33e02264d8ce1819c82fca6044f Mon Sep 17 00:00:00 2001 From: Benjamin Joel Musick Date: Wed, 16 Aug 2023 10:10:46 -0600 Subject: [PATCH 006/405] typo in error message --- singularity-eos/eos/eos_carnahan_starling.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_carnahan_starling.hpp b/singularity-eos/eos/eos_carnahan_starling.hpp index ae8c7107b8..1ee06603ea 100644 --- a/singularity-eos/eos/eos_carnahan_starling.hpp +++ b/singularity-eos/eos/eos_carnahan_starling.hpp @@ -103,7 +103,7 @@ class CarnahanStarling : public EosBase { if (status != Status::SUCCESS) { // Root finder failed even though the solution was bracketed... this is an error EOS_ERROR("*** (Warning) DensityFromPressureTemperature :: Convergence not met in " - "Carnahan-Starling EoS (regula_falsi) ***\n"); + "Carnahan-Starling EoS (root finder util) ***\n"); real_root = -1.0; } #else From f47d8d5eb8a3a25c8092690a51e00265649f9a05 Mon Sep 17 00:00:00 2001 From: Benjamin Joel Musick Date: Wed, 16 Aug 2023 12:16:53 -0600 Subject: [PATCH 007/405] removed newton solve, updated docs --- doc/sphinx/src/models.rst | 2 +- singularity-eos/eos/eos_carnahan_starling.hpp | 28 ------------------- 2 files changed, 1 insertion(+), 29 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index c73f501185..133aac7849 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -587,7 +587,7 @@ reference density is given by P_0 = \rho_0 Z(\rho_0) C_V T_0(\gamma-1), -where :math:`P_0` is by default 1 bar. Denisty is obtained via Newton iteration. +where :math:`P_0` is by default 1 bar. Denisty is obtained through root finding methods. The settable parameters for this EOS are :math:`\gamma-1`, specific heat capacity (:math:`C_V`), covolume (:math:`b`) and offset internal energy (:math:`q`). Optionally, the reference state for the entropy calculation can diff --git a/singularity-eos/eos/eos_carnahan_starling.hpp b/singularity-eos/eos/eos_carnahan_starling.hpp index 1ee06603ea..21b38c636c 100644 --- a/singularity-eos/eos/eos_carnahan_starling.hpp +++ b/singularity-eos/eos/eos_carnahan_starling.hpp @@ -88,8 +88,6 @@ class CarnahanStarling : public EosBase { const Real press, const Real temperature, const Real guess = robust::SMALL(), Real *lambda = nullptr) const { Real real_root; -#ifdef _SINGULARITY_EOS_UTILS_ROOT_FINDING_HPP_ - // use singularity utilities if available auto poly = [=](Real dens) { return _Cv * temperature * _gm1 * dens * ZedFromDensity(dens); }; @@ -106,32 +104,6 @@ class CarnahanStarling : public EosBase { "Carnahan-Starling EoS (root finder util) ***\n"); real_root = -1.0; } -#else - // newton-raphson iteration if utilities not available - int newton_counter = 0; - static constexpr int newton_max = 50; - static constexpr Real rel_tol = 1.e-12; - real_root = guess; - for (int i = 0; i < newton_max; i++) { - Real real_root_old = real_root; - Real newton_f = - _Cv * temperature * _gm1 * real_root_old * ZedFromDensity(real_root_old) - - press; - Real newton_f_prime = - _Cv * temperature * _gm1 * PartialRhoZedFromDensity(real_root_old); - real_root = real_root_old - robust::ratio(newton_f, newton_f_prime); - Real rel_err = std::abs(robust::ratio(real_root - real_root_old, real_root)); - newton_counter++; - if (rel_err < rel_tol) { - return std::max(robust::SMALL(), real_root); - } - } - if (newton_counter == newton_max) { - EOS_ERROR("*** (Warning) DensityFromPressureTemperature :: Convergence not met in " - "Carnahan-Starling EoS (newton-raphson) ***\n"); - real_root = -1.0; - } -#endif return std::max(robust::SMALL(), real_root); } PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( From ce22f21dd6abe2ced1669592762f1021b15f2b2c Mon Sep 17 00:00:00 2001 From: Benjamin Joel Musick Date: Thu, 17 Aug 2023 16:21:22 -0600 Subject: [PATCH 008/405] dholladay optimizations --- singularity-eos/eos/eos_carnahan_starling.hpp | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/singularity-eos/eos/eos_carnahan_starling.hpp b/singularity-eos/eos/eos_carnahan_starling.hpp index 21b38c636c..88bd8d129b 100644 --- a/singularity-eos/eos/eos_carnahan_starling.hpp +++ b/singularity-eos/eos/eos_carnahan_starling.hpp @@ -26,6 +26,7 @@ // Base stuff #include #include +#include #include #include #include @@ -72,17 +73,14 @@ class CarnahanStarling : public EosBase { PORTABLE_INLINE_FUNCTION Real ZedFromDensity(const Real rho, Real *lambda = nullptr) const { const Real eta = _bb * rho; - const Real eta2 = eta * eta; - const Real zed = robust::ratio(1.0 + eta + eta2 - eta2 * eta, - (1.0 - eta) * (1.0 - eta) * (1.0 - eta)); + const Real zed = + 1.0 + robust::ratio(eta * (4.0 - 2.0 * eta), math_utils::pow<3>(1.0 - eta)); return zed; } PORTABLE_INLINE_FUNCTION Real PartialRhoZedFromDensity(const Real rho, Real *lambda = nullptr) const { const Real eta = _bb * rho; - const Real eta2 = eta * eta; - return robust::ratio(eta2 * eta2 - 4.0 * eta2 * eta + 4.0 * eta2 + 4.0 * eta + 1.0, - (1.0 - eta) * (1.0 - eta) * (1.0 - eta) * (1.0 - eta)); + return 1.0 + robust::ratio(eta * (8.0 - 2.0 * eta), math_utils::pow<4>(1.0 - eta)); } PORTABLE_INLINE_FUNCTION Real DensityFromPressureTemperature( const Real press, const Real temperature, const Real guess = robust::SMALL(), From 86aa660d37e1f89655538f21661ac528559c3328 Mon Sep 17 00:00:00 2001 From: Benjamin Joel Musick Date: Thu, 17 Aug 2023 16:47:10 -0600 Subject: [PATCH 009/405] dholladay optimizations --- singularity-eos/eos/eos_carnahan_starling.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/singularity-eos/eos/eos_carnahan_starling.hpp b/singularity-eos/eos/eos_carnahan_starling.hpp index 88bd8d129b..4e242b5712 100644 --- a/singularity-eos/eos/eos_carnahan_starling.hpp +++ b/singularity-eos/eos/eos_carnahan_starling.hpp @@ -121,23 +121,23 @@ class CarnahanStarling : public EosBase { PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { const Real vol = robust::ratio(1.0, rho); + const Real one_by_vb = robust::ratio(1.0, vol - _bb); + const Real one_by_v0b = robust::ratio(1.0, _vol0 - _bb); return _Cv * std::log(robust::ratio(temperature, _T0) + robust::SMALL()) + _gm1 * _Cv * std::log(robust::ratio(vol, _vol0) + robust::SMALL()) - - _gm1 * _Cv * _bb * - (4.0 * (robust::ratio(1.0, vol - _bb) - robust::ratio(1.0, _vol0 - _bb)) + - _bb * (robust::ratio(1.0, (vol - _bb) * (vol - _bb)) - - robust::ratio(1.0, (_vol0 - _bb) * (_vol0 - _bb)))) + + _gm1 * _Cv * _bb * (one_by_vb - one_by_v0b) * + (4.0 + _bb * (one_by_vb + one_by_v0b)) + _qp; } PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const { const Real vol = robust::ratio(1.0, rho); + const Real one_by_vb = robust::ratio(1.0, vol - _bb); + const Real one_by_v0b = robust::ratio(1.0, _vol0 - _bb); return _Cv * std::log(robust::ratio(sie - _qq, _sie0 - _qq) + robust::SMALL()) + _gm1 * _Cv * std::log(robust::ratio(vol, _vol0) + robust::SMALL()) - - _gm1 * _Cv * _bb * - (4.0 * (robust::ratio(1.0, vol - _bb) - robust::ratio(1.0, _vol0 - _bb)) + - _bb * (robust::ratio(1.0, (vol - _bb) * (vol - _bb)) - - robust::ratio(1.0, (_vol0 - _bb) * (_vol0 - _bb)))) + + _gm1 * _Cv * _bb * (one_by_vb - one_by_v0b) * + (4.0 + _bb * (one_by_vb + one_by_v0b)) + _qp; } PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( From 351afd308e4c041f0a10f1b6045b80b2dcb174a2 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Tue, 5 Dec 2023 11:52:51 -0700 Subject: [PATCH 010/405] Added files for MGUsup EOS. --- singularity-eos/eos/eos.hpp | 4 +- singularity-eos/eos/eos_mgusup.hpp | 386 +++++++++++++++++++++++ test/CMakeLists.txt | 1 + test/test_eos_mgusup.cpp | 480 +++++++++++++++++++++++++++++ 4 files changed, 869 insertions(+), 2 deletions(-) create mode 100644 singularity-eos/eos/eos_mgusup.hpp create mode 100644 test/test_eos_mgusup.cpp diff --git a/singularity-eos/eos/eos.hpp b/singularity-eos/eos/eos.hpp index 69a839235a..ac929d9d87 100644 --- a/singularity-eos/eos/eos.hpp +++ b/singularity-eos/eos/eos.hpp @@ -45,7 +45,7 @@ #include #include #include - +#include // Modifiers #include #include @@ -67,7 +67,7 @@ using singularity::variadic_utils::transform_variadic_list; // all eos's static constexpr const auto full_eos_list = - tl +#include + +#include + +#include +#include +#include +#include +#include + +namespace singularity { + +using namespace eos_base; +using namespace robust; + +// COMMENT: This is a complete, thermodynamically consistent Mie-Gruneisen +// EOS based on a linear Us-up Hugoniot, Us=Cs+s*up. +class MGUsup : public EosBase { + public: + MGUsup() = default; + // Constructor + MGUsup(const Real rho0, const Real T0, const Real B0, const Real BP0, const Real A0, + const Real Cv0, const Real E0, const Real S0, const Real *expconsts) + : _rho0(rho0), _T0(T0), _B0(B0), _BP0(BP0), _A0(A0), _Cv0(Cv0), _E0(E0), _S0(S0) { + CheckMGUsup(); + InitializeMGUsup(expconsts); + } + + MGUsup GetOnDevice() { return *this; } + PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real + MinInternalEnergyFromDensity(const Real rho, Real *lambda = nullptr) const; + // Entropy added AEM Dec. 2022 + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const { + return _Cv0; + } + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return _Cv0; + } + // Thermal Bulk Modulus added AEM Dec 2022 + PORTABLE_INLINE_FUNCTION Real TBulkModulusFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const; + // Thermal expansion coefficient added AEM 2022 + PORTABLE_INLINE_FUNCTION Real TExpansionCoeffFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( + const Real rho, const Real temp, Real *lambda = nullptr) const { + return robust::ratio(_A0 * _B0, _Cv0 * rho); + } + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return robust::ratio(_A0 * _B0, _Cv0 * rho); + } + PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, + Real &cv, Real &bmod, const unsigned long output, + Real *lambda = nullptr) const; + PORTABLE_INLINE_FUNCTION + void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, + Real &bmod, Real &dpde, Real &dvdt, + Real *lambda = nullptr) const; + // Generic functions provided by the base class. These contain e.g. the vector + // overloads that use the scalar versions declared here + SG_ADD_BASE_CLASS_USINGS(MGUsup) + PORTABLE_INLINE_FUNCTION + int nlambda() const noexcept { return 0; } + static constexpr unsigned long PreferredInput() { return _preferred_input; } + static inline unsigned long scratch_size(std::string method, unsigned int nelements) { + return 0; + } + static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } + PORTABLE_INLINE_FUNCTION void PrintParams() const { + static constexpr char st[]{"MGUsup Params: "}; + printf("%s rho0:%e T0:%e B0:%e BP0:%e\n A0:%e Cv0:%e E0:%e S0:%e\n" + "non-zero elements in d2tod40 array:\n", + st, _rho0, _T0, _B0, _BP0, _A0, _Cv0, _E0, _S0); + for (int i = 0; i < 39; i++) { + if (_d2tod40[i] > 0.0) printf("d%i:%e\t", i + 2, _d2tod40[i]); + } + printf("\n\n"); + } + // Density/Energy from P/T not unique, if used will give error + PORTABLE_INLINE_FUNCTION void + DensityEnergyFromPressureTemperature(const Real press, const Real temp, Real *lambda, + Real &rho, Real &sie) const; + inline void Finalize() {} + static std::string EosType() { return std::string("MGUsup"); } + static std::string EosPyType() { return EosType(); } + + private: + static constexpr const unsigned long _preferred_input = + thermalqs::density | thermalqs::temperature; + Real _rho0, _T0, _B0, _BP0, _A0, _Cv0, _E0, _S0; + static constexpr const int PressureCoeffsd2tod40Size = 39; + static constexpr const int MGUsupInternalParametersSize = PressureCoeffsd2tod40Size + 4; + Real _VIP[MGUsupInternalParametersSize], _d2tod40[PressureCoeffsd2tod40Size]; + void CheckMGUsup(); + void InitializeMGUsup(const Real *expcoeffs); + PORTABLE_INLINE_FUNCTION void MGUsup_F_DT_func(const Real rho, const Real T, + Real *output) const; +}; + +inline void MGUsup::CheckMGUsup() { + + if (_rho0 < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter rho0 < 0"); + } + if (_T0 < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter T0 < 0"); + } + if (_B0 < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter B0 < 0"); + } + if (_BP0 < 1.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter BP0 < 1"); + } + if (_A0 < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter A0 < 0"); + } + if (_Cv0 < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter Cv0 < 0"); + } +} + +inline void MGUsup::InitializeMGUsup(const Real *d2tod40input) { + + // The PressureCoeffsd2tod40Size (=39) allowed d2 to d40 coefficients + // for the pressure reference curve vs rho + // are seldom all used so did not want to crowd the argument list with them. + // Instead I ask the host code to send me a pointer to this array so that I can + // copy it here. Not used coeffs should be set to 0.0 (of course). + for (int ind = 0; ind < PressureCoeffsd2tod40Size; ind++) { + _d2tod40[ind] = d2tod40input[ind]; + } + // Put a couple of much used model parameter combinations and + // the energy coefficients f0 to f40 in an internal parameters array + // of size 2+41=MGUsupInternalParametersSize. + _VIP[0] = robust::ratio(_B0, _rho0) + + robust::ratio(_A0 * _A0 * _B0 * _B0 * _T0, + _rho0 * _rho0 * _Cv0); // sound speed squared + _VIP[1] = 3.0 / 2.0 * (_BP0 - 1.0); // exponent eta0 + // initializing with pressure coeffs to get energy coeffs into VIP + _VIP[2] = 1.0; // prefactor d0 + _VIP[3] = 0.0; // prefactor d1 + for (int ind = 4; ind < MGUsupInternalParametersSize; ind++) { //_d2tod40[0]=d2 + _VIP[ind] = _d2tod40[ind - 4]; // dn is in VIP[n+2] + } + for (int ind = MGUsupInternalParametersSize - 2; ind >= 2; + ind--) { // _VIP[42]=d40=f40 given,first calculated is _VIP[41]=f39 + _VIP[ind] = _VIP[ind] - (ind) / _VIP[1] * _VIP[ind + 1]; // prefactors f40 to f0 + } // _VIP[n+2]=fn, ind=n+2 +} + +PORTABLE_INLINE_FUNCTION void MGUsup::MGUsup_F_DT_func(const Real rho, const Real T, + Real *output) const { + constexpr int pref0vp = -2, pref0vip = 2, maxind = MGUsupInternalParametersSize - 3; + Real sumP = 0.0, sumB = 0.0, sumE = 0.0; + + Real x = std::cbrt(robust::ratio(_rho0, rho)); /*rho-dependent*/ + Real x2inv = robust::ratio(1.0, x * x); + Real onemx = 1.0 - x; + Real etatimes1mx = _VIP[1] * onemx; + Real expetatimes1mx = exp(etatimes1mx); + +#pragma unroll + for (int ind = maxind; ind >= 2; ind--) { //_d2tod40[0]=d2 + sumP = _d2tod40[pref0vp + ind] + onemx * sumP; //_d2tod40[38]=d40 + sumB = _d2tod40[pref0vp + ind] * ind + onemx * sumB; //_d2tod40[-2+40]=d40 + sumE = _VIP[pref0vip + ind] + onemx * sumE; //_VIP[42]=f40 + } //_VIP[2]=f0 +#pragma unroll + for (int ind = 1; ind >= 0; ind--) { + sumE = _VIP[pref0vip + ind] + onemx * sumE; + } + sumP = 1.0 + sumP * (onemx * onemx); + sumB = sumB * onemx; + sumE = _VIP[1] * onemx * sumE; + + /* T0 isotherm */ + Real energy = 9.0 * robust::ratio(_B0, _VIP[1] * _VIP[1] * _rho0) * + (_VIP[pref0vip] - expetatimes1mx * (_VIP[pref0vip] - sumE)) - + (_A0 * _B0) * (robust::ratio(_T0, _rho0) - robust::ratio(_T0, rho)); + Real pressure = 3.0 * _B0 * x2inv * onemx * expetatimes1mx * sumP; + Real temp = (1.0 + onemx * (_VIP[1] * x + 1.0)) * sumP + x * onemx * sumB; + temp = robust::ratio(_B0, _rho0) * x * expetatimes1mx * temp; + + /* Go to required temperature */ + energy = energy + _Cv0 * (T - _T0) + _E0; + pressure = pressure + (_A0 * _B0) * (T - _T0); + Real dpdrho = temp; + Real dpdt = _A0 * _B0; + Real dedt = _Cv0; + Real dedrho = robust::ratio(pressure - T * (_A0 * _B0), rho * rho); + Real entropy; + if (T < 0.0) { +#ifndef NDEBUG + PORTABLE_WARN("Negative temperature input"); +#endif // NDEBUG + entropy = 0.0; + } else { + entropy = (_A0 * _B0) * (robust::ratio(1.0, rho) - robust::ratio(1, _rho0)) + + _Cv0 * std::log(robust::ratio(T, _T0)) + _S0; + } + Real soundspeed = + std::max(0.0, dpdrho + T * robust::ratio(dpdt * dpdt, rho * rho * _Cv0)); + soundspeed = std::sqrt(soundspeed); + + output[0] = energy; + output[1] = pressure; + output[2] = dpdrho; + output[3] = dpdt; + output[4] = dedt; + output[5] = dedrho; + output[6] = entropy; + output[7] = soundspeed; + + return; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::InternalEnergyFromDensityTemperature( + const Real rho, const Real temp, Real *lambda) const { + Real output[8]; + MGUsup_F_DT_func(rho, temp, output); + return output[0]; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::PressureFromDensityTemperature(const Real rho, + const Real temp, + Real *lambda) const { + Real output[8]; + MGUsup_F_DT_func(rho, temp, output); + return output[1]; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::EntropyFromDensityTemperature(const Real rho, + const Real temp, + Real *lambda) const { + Real output[8]; + MGUsup_F_DT_func(rho, temp, output); + return output[6]; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::TExpansionCoeffFromDensityTemperature( + const Real rho, const Real temp, Real *lambda) const { + Real output[8]; + MGUsup_F_DT_func(rho, temp, output); + return robust::ratio(output[3], output[2] * rho); +} +PORTABLE_INLINE_FUNCTION Real MGUsup::TBulkModulusFromDensityTemperature( + const Real rho, const Real temp, Real *lambda) const { + Real output[8]; + MGUsup_F_DT_func(rho, temp, output); + return output[2] * rho; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::BulkModulusFromDensityTemperature( + const Real rho, const Real temp, Real *lambda) const { + Real output[8]; + MGUsup_F_DT_func(rho, temp, output); + return output[7] * output[7] * rho; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::TemperatureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda) const { + Real Tref; + Real output[8]; + Tref = _T0; + MGUsup_F_DT_func(rho, Tref, output); + return robust::ratio(sie - output[0], _Cv0) + Tref; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::PressureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda) const { + Real temp; + Real output[8]; + temp = TemperatureFromDensityInternalEnergy(rho, sie); + MGUsup_F_DT_func(rho, temp, output); + return output[1]; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::MinInternalEnergyFromDensity(const Real rho, + Real *lambda) const { + MinInternalEnergyIsNotEnabled("MGUsup"); + return 0.0; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda) const { + Real temp; + Real output[8]; + temp = TemperatureFromDensityInternalEnergy(rho, sie); + MGUsup_F_DT_func(rho, temp, output); + return output[6]; +} +PORTABLE_INLINE_FUNCTION Real MGUsup::BulkModulusFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda) const { + Real temp; + Real output[8]; + temp = TemperatureFromDensityInternalEnergy(rho, sie); + MGUsup_F_DT_func(rho, temp, output); + return output[7] * output[7] * rho; +} +// AEM: Give error since function is not well defined +PORTABLE_INLINE_FUNCTION void +MGUsup::DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Real *lambda, Real &rho, Real &sie) const { + EOS_ERROR("MGUsup::DensityEnergyFromPressureTemperature: " + "Not implemented.\n"); +} +// AEM: We should add entropy and Gruneissen parameters here so that it is complete +// If we add also alpha and BT, those should also be in here. +PORTABLE_INLINE_FUNCTION void MGUsup::FillEos(Real &rho, Real &temp, Real &sie, + Real &press, Real &cv, Real &bmod, + const unsigned long output, + Real *lambda) const { + const unsigned long input = ~output; // everything that is not output is input + if (thermalqs::density & output) { + EOS_ERROR("MGUsup FillEos: Density is required input.\n"); + } else if (thermalqs::temperature & output && + thermalqs::specific_internal_energy & output) { + EOS_ERROR("MGUsup FillEos: Density and Internal Energy or Density and Temperature " + "are required input parameters.\n"); + } + if (thermalqs::specific_internal_energy & input) { + temp = TemperatureFromDensityInternalEnergy(rho, sie); + } + Real Vout[8]; + MGUsup_F_DT_func(rho, temp, Vout); + if (output & thermalqs::temperature) temp = temp; + if (output & thermalqs::specific_internal_energy) sie = Vout[0]; + if (output & thermalqs::pressure) press = Vout[1]; + if (output & thermalqs::specific_heat) cv = Vout[4]; + if (output & thermalqs::bulk_modulus) bmod = Vout[7] * Vout[7] * rho; +} + +// TODO(JMM): pre-cache these rather than recomputing them each time +PORTABLE_INLINE_FUNCTION +void MGUsup::ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, + Real &cv, Real &bmod, Real &dpde, Real &dvdt, + Real *lambda) const { + // AEM: Added all variables I think should be output eventually + Real tbmod; + // Real entropy, alpha, Gamma; + + rho = _rho0; + temp = _T0; + sie = _E0; + press = PressureFromDensityTemperature(rho, temp, lambda); + // entropy = _S0; + cv = _Cv0; + tbmod = _B0; + // alpha = _A0; + bmod = BulkModulusFromDensityTemperature(rho, temp, lambda); + // Gamma = robust::ratio(_A0 * _B0, _Cv0 * _rho0); + // AEM: I suggest taking the two following away. + dpde = robust::ratio(_A0 * tbmod, _Cv0); + dvdt = robust::ratio(_A0, _rho0); +} +} // namespace singularity + +#endif // _SINGULARITY_EOS_EOS_MGUSUP_HPP_ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d63bbbb51b..6add5dbd39 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -21,6 +21,7 @@ add_executable( test_eos_noble_abel.cpp test_eos_stiff.cpp test_eos_vinet.cpp + test_eos_mgusup.cpp ) add_executable( diff --git a/test/test_eos_mgusup.cpp b/test/test_eos_mgusup.cpp new file mode 100644 index 0000000000..b43194f3d7 --- /dev/null +++ b/test/test_eos_mgusup.cpp @@ -0,0 +1,480 @@ +//------------------------------------------------------------------------------ +// © 2021-2023. 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. +//------------------------------------------------------------------------------ + +#include +#include +#include +#include +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +#include +#include +#include +#include +#include +#include + +using singularity::MGUsup; +using EOS = singularity::Variant; + +SCENARIO("MGUsup EOS rho sie", "[VectorEOS][MGUsupEOS]") { + GIVEN("Parameters for a MGUsup EOS") { + // Unit conversions + constexpr Real Mbcc_per_g = 1e12; + // MGUsup parameters for copper + constexpr Real rho0 = 8.93; + constexpr Real T0 = 298.0; + constexpr Real B0 = 1.3448466 * Mbcc_per_g; + constexpr Real BP0 = 4.956; + constexpr Real A0 = 5.19245e-05; + constexpr Real Cv0 = 0.383e-05 * Mbcc_per_g; + constexpr Real E0 = 0.0; + constexpr Real S0 = 5.05e-04 * Mbcc_per_g; + constexpr Real d2to40[39] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; + // Create the EOS + EOS host_eos = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40); + EOS eos = host_eos.GetOnDevice(); + MGUsup host_eos2 = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40); + MGUsup eos2 = host_eos2.GetOnDevice(); + + eos.PrintParams(); + + GIVEN("Densities and energies") { + constexpr int num = 4; +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_density("density"); + // AEM: I think it should not be "density" again below + Kokkos::View v_energy("density"); +#else + // Otherwise just create arrays to contain values and create pointers to + // be passed to the functions in place of the Kokkos views + std::array density; + std::array energy; + auto v_density = density.data(); + auto v_energy = energy.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize density and energy", 0, 1, PORTABLE_LAMBDA(int i) { + v_density[0] = 8.0; + v_density[1] = 8.5; + v_density[2] = 9.0; + v_density[3] = 8.93; + v_energy[0] = 2.e8; + v_energy[1] = 4.6e9; + v_energy[2] = 4.7e9; + v_energy[3] = 0.0; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create host-side mirrors of the inputs and copy the inputs. These are + // just used for the comparisons + auto density = Kokkos::create_mirror_view(v_density); + auto energy = Kokkos::create_mirror_view(v_energy); + Kokkos::deep_copy(density, v_density); + Kokkos::deep_copy(energy, v_energy); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Gold standard values for a subset of lookups + constexpr std::array temperature_true{ + 6.638294350641033e+1, 1.42266691572769e+03, 1.52867838544081e+03, T0}; + constexpr std::array pressure_true{ + -1.283459624233147e+11, 1.98538351697004e+10, 9.66446220551959e+10, 0.}; + constexpr std::array entropy_true{ + 5.0015771521255749e+08, 5.1138262492866594e+08, 5.1120147992457777e+08, S0}; + constexpr std::array cv_true{Cv0, Cv0, Cv0, Cv0}; + constexpr std::array bulkmodulus_true{ + 7.44411028807209e+11, 1.254880120063067e+12, 1.613822819346433e+12, + (B0 + B0 * B0 * A0 * A0 / Cv0 / rho0 * T0)}; + constexpr std::array gruneisen_true{ + B0 * A0 / Cv0 / 8.0, B0 * A0 / Cv0 / 8.5, B0 * A0 / Cv0 / 9.0, + B0 * A0 / Cv0 / rho0}; + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create device views for outputs and mirror those views on the host + Kokkos::View v_temperature("Temperature"); + Kokkos::View v_pressure("Pressure"); + Kokkos::View v_entropy("Entropy"); + Kokkos::View v_cv("Cv"); + Kokkos::View v_bulkmodulus("bmod"); + Kokkos::View v_gruneisen("Gamma"); + auto h_temperature = Kokkos::create_mirror_view(v_temperature); + auto h_pressure = Kokkos::create_mirror_view(v_pressure); + auto h_entropy = Kokkos::create_mirror_view(v_entropy); + auto h_cv = Kokkos::create_mirror_view(v_cv); + auto h_bulkmodulus = Kokkos::create_mirror_view(v_bulkmodulus); + auto h_gruneisen = Kokkos::create_mirror_view(v_gruneisen); +#else + // Create arrays for the outputs and then pointers to those arrays that + // will be passed to the functions in place of the Kokkos views + std::array h_temperature; + std::array h_pressure; + std::array h_entropy; + std::array h_cv; + std::array h_bulkmodulus; + std::array h_gruneisen; + // Just alias the existing pointers + auto v_temperature = h_temperature.data(); + auto v_pressure = h_pressure.data(); + auto v_entropy = h_entropy.data(); + auto v_cv = h_cv.data(); + auto v_bulkmodulus = h_bulkmodulus.data(); + auto v_gruneisen = h_gruneisen.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + WHEN("A T(rho, e) lookup is performed") { + eos.TemperatureFromDensityInternalEnergy(v_density, v_energy, v_temperature, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_temperature, v_temperature); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned T(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_temperature, temperature_true, "Density", + "Energy"); + } + } + + WHEN("A P(rho, e) lookup is performed") { + eos.PressureFromDensityInternalEnergy(v_density, v_energy, v_pressure, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_pressure, v_pressure); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned P(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_pressure, pressure_true, "Density", + "Energy"); + } + } + + WHEN("A S(rho, e) lookup is performed") { + portableFor( + "Test entropy", 0, num, PORTABLE_LAMBDA(const int i) { + v_entropy[i] = + eos2.EntropyFromDensityInternalEnergy(v_density[i], v_energy[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_entropy, v_entropy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned S(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_entropy, entropy_true, "Density", + "Energy"); + } + } + + WHEN("A B_S(rho, e) lookup is performed") { + eos.BulkModulusFromDensityInternalEnergy(v_density, v_energy, v_bulkmodulus, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_bulkmodulus, v_bulkmodulus); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned B_S(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_bulkmodulus, bulkmodulus_true, "Density", + "Energy"); + } + } + + WHEN("A Cv(rho, e) lookup is performed") { + eos.SpecificHeatFromDensityInternalEnergy(v_density, v_energy, v_cv, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_cv, v_cv); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned T(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_cv, cv_true, "Density", "Energy"); + } + } + + WHEN("A Gamma(rho, e) lookup is performed") { + eos.GruneisenParamFromDensityInternalEnergy(v_density, v_energy, v_gruneisen, + num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_gruneisen, v_gruneisen); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned Gamma(rho, e) should be equal to the true value") { + array_compare(num, density, energy, h_gruneisen, gruneisen_true, "Density", + "Energy"); + } + } + } + } +} +SCENARIO("MGUsup EOS rho T", "[VectorEOS][MGUsupEOS]") { + GIVEN("Parameters for a MGUsup EOS") { + // Unit conversions + constexpr Real Mbcc_per_g = 1e12; + // MGUsup parameters for copper + constexpr Real rho0 = 8.93; + constexpr Real T0 = 298.0; + constexpr Real B0 = 1.3448466 * Mbcc_per_g; + constexpr Real BP0 = 4.956; + constexpr Real A0 = 5.19245e-05; + constexpr Real Cv0 = 0.383e-05 * Mbcc_per_g; + constexpr Real E0 = 0.0; + constexpr Real S0 = 5.05e-04 * Mbcc_per_g; + constexpr Real d2to40[39] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; + // Create the EOS + EOS host_eos = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40); + EOS eos = host_eos.GetOnDevice(); + MGUsup host_eos2 = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40); + MGUsup eos2 = host_eos2.GetOnDevice(); + + eos.PrintParams(); + + GIVEN("Densities and temperatures") { + constexpr int num = 4; +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_density("density"); + Kokkos::View v_temperature("temperature"); +#else + // Otherwise just create arrays to contain values and create pointers to + // be passed to the functions in place of the Kokkos views + std::array density; + std::array temperature; + auto v_density = density.data(); + auto v_temperature = temperature.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize density and temperature", 0, 1, PORTABLE_LAMBDA(int i) { + v_density[0] = 8.0; + v_density[1] = 8.5; + v_density[2] = 9.0; + v_density[3] = 8.93; + v_temperature[0] = 66.383; + v_temperature[1] = 1422.7; + v_temperature[2] = 1528.7; + v_temperature[3] = 298.0; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create host-side mirrors of the inputs and copy the inputs. These are + // just used for the comparisons + auto density = Kokkos::create_mirror_view(v_density); + auto temperature = Kokkos::create_mirror_view(v_temperature); + Kokkos::deep_copy(density, v_density); + Kokkos::deep_copy(temperature, v_temperature); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Gold standard values for a subset of lookups + constexpr std::array sie_true{ + 2.0000021637044835e+8, 4.6001267127629471e+09, 4.7000827837617149e+09, 0.0}; + constexpr std::array pressure_true{ + -1.283459584783398e+11, 1.98561454605571e+10, 9.66461314103968e+10, 0.}; + constexpr std::array entropy_true{ + 5.0015771847198445e+08, 5.1138271399469268e+08, 5.1120153407800680e+08, S0}; + constexpr std::array cv_true{Cv0, Cv0, Cv0, Cv0}; + constexpr std::array tbulkmodulus_true{ + 7.3384631127398950e+11, 1.0417839336794381e+12, 1.3975684023296028e+12, B0}; + constexpr std::array alpha_true{ + 9.5156828083623124e-05, 6.7029721830195901e-05, 4.9965702691402983e-05, A0}; + // B0 * A0 / tbulkmodulus_true[0], B0 * A0 / tbulkmodulus_true[1], B0 * A0 / + // tbulkmodulus_true[2] + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create device views for outputs and mirror those views on the host + Kokkos::View v_energy("Energy"); + Kokkos::View v_pressure("Pressure"); + Kokkos::View v_entropy("Entropy"); + Kokkos::View v_cv("Cv"); + Kokkos::View v_tbulkmodulus("tbmod"); + Kokkos::View v_alpha("alpha"); + auto h_energy = Kokkos::create_mirror_view(v_energy); + auto h_pressure = Kokkos::create_mirror_view(v_pressure); + auto h_entropy = Kokkos::create_mirror_view(v_entropy); + auto h_cv = Kokkos::create_mirror_view(v_cv); + auto h_tbulkmodulus = Kokkos::create_mirror_view(v_tbulkmodulus); + auto h_alpha = Kokkos::create_mirror_view(v_alpha); +#else + // Create arrays for the outputs and then pointers to those arrays that + // will be passed to the functions in place of the Kokkos views + std::array h_energy; + std::array h_pressure; + std::array h_entropy; + std::array h_cv; + std::array h_tbulkmodulus; + std::array h_alpha; + // Just alias the existing pointers + auto v_energy = h_energy.data(); + auto v_pressure = h_pressure.data(); + auto v_entropy = h_entropy.data(); + auto v_cv = h_cv.data(); + auto v_tbulkmodulus = h_tbulkmodulus.data(); + auto v_alpha = h_alpha.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + WHEN("A E(rho, T) lookup is performed") { + eos.InternalEnergyFromDensityTemperature(v_density, v_temperature, v_energy, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_energy, v_energy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned E(rho, T) should be equal to the true value") { + array_compare(num, density, temperature, h_energy, sie_true, "Density", + "Temperature"); + } + } + + WHEN("A P(rho, T) lookup is performed") { + eos.PressureFromDensityTemperature(v_density, v_temperature, v_pressure, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_pressure, v_pressure); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned P(rho, T) should be equal to the true value") { + array_compare(num, density, temperature, h_pressure, pressure_true, "Density", + "Temperature"); + } + } + + portableFor( + "entropy", 0, num, PORTABLE_LAMBDA(const int i) { + v_entropy[i] = + eos2.EntropyFromDensityInternalEnergy(v_density[i], v_energy[i]); + }); + WHEN("A S(rho, T) lookup is performed") { + portableFor( + "entropy", 0, num, PORTABLE_LAMBDA(const int i) { + v_entropy[i] = + eos2.EntropyFromDensityTemperature(v_density[i], v_temperature[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_entropy, v_entropy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned S(rho, T) should be equal to the true value") { + array_compare(num, density, temperature, h_entropy, entropy_true, "Density", + "Temperature"); + } + } + + WHEN("A B_T(rho, T) lookup is performed") { + portableFor( + "bulk modulus", 0, num, PORTABLE_LAMBDA(const int i) { + v_tbulkmodulus[i] = + eos2.TBulkModulusFromDensityTemperature(v_density[i], v_temperature[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_tbulkmodulus, v_tbulkmodulus); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned B_T(rho, T) should be equal to the true value") { + array_compare(num, density, temperature, h_tbulkmodulus, tbulkmodulus_true, + "Density", "Temperature"); + } + } + + WHEN("A Cv(rho, T) lookup is performed") { + eos.SpecificHeatFromDensityTemperature(v_density, v_temperature, v_cv, num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_cv, v_cv); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned T(rho, e) should be equal to the true value") { + array_compare(num, density, temperature, h_cv, cv_true, "Density", + "Temperature"); + } + } + + WHEN("A alpha(rho, T) lookup is performed") { + portableFor( + "alpha", 0, num, PORTABLE_LAMBDA(const int i) { + v_alpha[i] = eos2.TExpansionCoeffFromDensityTemperature(v_density[i], + v_temperature[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_alpha, v_alpha); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned alpha(rho, T) should be equal to the true value") { + array_compare(num, density, temperature, h_alpha, alpha_true, "Density", + "Temperature"); + } + } + } + } +} +SCENARIO("MGUsup EOS SetUp", "[VectorEOS][MGUsupEOS]") { + GIVEN("Parameters for a MGUsup EOS") { + // Unit conversions + constexpr Real Mbcc_per_g = 1e12; + // MGUsup parameters for copper + Real rho0 = 8.93; + Real T0 = 298.0; + Real B0 = 1.3448466 * Mbcc_per_g; + Real BP0 = 4.956; + Real A0 = 5.19245e-05; + Real Cv0 = 0.383e-05 * Mbcc_per_g; + constexpr Real E0 = 0.0; + constexpr Real S0 = 5.05e-04 * Mbcc_per_g; + constexpr Real d2to40[39] = {0., 0., 0.0000000000001, + 0., 0., 0., + 0., 0., 0., + 0., 0., 0., + 0., 0., 0., + 0., 0., 0., + 0., 0., 0., + 0., 0., 0., + 0., 0., 0., + 0., 0., 0., + 0., 0., 0., + 0., 1.0e-26, 0., + 0., 0., 0.}; + // Create the EOS + EOS host_eos = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40); + EOS eos = host_eos.GetOnDevice(); + eos.Finalize(); + + WHEN("Faulty/not set parameter rho0 is given") { + rho0 = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + THEN("An error message should be written out") {} + } + WHEN("Faulty/not set parameter T0 is given") { + T0 = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + THEN("An error message should be written out") {} + } + WHEN("Faulty/not set parameter B0 is given") { + B0 = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + THEN("An error message should be written out") {} + } + WHEN("Faulty/not set parameter BP0 is given") { + BP0 = 0.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + THEN("An error message should be written out") {} + } + WHEN("Faulty/not set parameter Cv0 is given") { + Cv0 = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + THEN("An error message should be written out") {} + } + WHEN("Faulty/not set parameter A0 is given") { + A0 = -10000000.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + THEN("An error message should be written out") {} + } + } +} From be51813ec97353e46df19bdbb90ed5a956807a73 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Tue, 2 Jan 2024 18:55:09 -0700 Subject: [PATCH 011/405] Included MGUsup EOS. --- doc/sphinx/src/models.rst | 101 +++++++++ singularity-eos/eos/eos.hpp | 2 +- singularity-eos/eos/eos_mgusup.hpp | 335 ++++++++++++++--------------- test/test_eos_mgusup.cpp | 197 ++++++++++------- 4 files changed, 381 insertions(+), 254 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 26846724af..3e5819aca3 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -751,6 +751,107 @@ is :math:`S_0`, and ``expconsts`` is a pointer to the constant array of length :math:`d_2` to :math:`d_{40}`. Expansion coefficients not used should be set to 0.0. +Mie-Gruneisen linear :math:`U_s`-:math:`u_p` EOS +```````````````````````````````````````````````` + +One of the most commonly-used EOS is the linear :math:`U_s`-:math:`u_p` version of the Mie-Gruneisen EOS. This EOS +uses the Hugoniot as the reference curve and is extensively used in shock physics. +This version implements the exact thermodynamic tempearture on the Hugoniot and also adds an entropy. + +The pressure follows the traditional Mie-Gruneisen form, + +.. math:: + + P(\rho, e) = P_H(\rho) + \rho\Gamma(\rho) \left(e - e_H(\rho) \right), + +Here the subscript :math:`H` is a reminder that the reference curve is a +Hugoniot. :math:`\Gamma` is the Gruneisen parameter and the first approximation +is that :math:`\rho\Gamma(\rho)=\rho_0\Gamma(\rho_0)` +which is the same assumption as in the Gruneisen EOS when :math:`b=0`. + +The above is an incomplete equation of state because it only relates the +pressure to the density and energy, the minimum required in a solution to the +Euler equations. To complete the EOS and determine the temperature and entropy, a constant +heat capacity is assumed so that + +.. math:: + + T(\rho, e) = \frac{\left(e-e_H(\rho)\right)}{C_V} + T_H(\rho) + +Note the difference from the Gruneisen EOS described above. We still use a constant :math:`C_V`, +and it is usually taken at the reference temperature, but +we now extrapolate from the temperature on the Hugoniot, :math: `T_H(\rho)`, and not +from the reference temperature, :math:`T_0`. + +With this consistent temperature we can derive an entropy in a similar way as for the Vinet EOS, +.. math:: + + S(\rho,T) = S_0 - \Gamma(\rho_0)C_V \eta + {C_V} \ln \frac{T}{T_{ref}} , + + +where :math:`\eta` is a measure of compression given by +.. math:: + + \eta = 1 - \frac{\rho_0}{\rho}. + +This is convenient because :math:`eta = 0` when :math:`\rho = \rho_0`, +:math:`\eta = 1` at the infinite density limit, and :math:`\eta = -\infty` at +the zero density limit. + +The pressure, energy, and temperature, on the Hugoniot are derived from the +shock jump conditions assuming a linear :math:`U_s`-:math:`u_p` relation, +.. math:: + + U_s = C_s + s u_p . + +Here :math:`U_s` is the shock velocity and :math:`u_p` is the particle +velocity. As is pointed out in the description of the Gruneisen EOS, +for many materials, the :math:`U_s`-:math:`u_p` relationship is roughly linear +so only this :math:`s` parameter is needed. The units for :math:`C_s` is velocity while +:math:`s` is unitless. Note that the parameter :math:`s` is related to the +fundamental derivative of shock physics as shown by `Wills `_. + +The reference pressure along the Hugoniot is determined by + +.. math:: + + P_H(\rho) = C_s^2 \rho_0 \frac{\eta}{\left(1 - s \eta \right)^2} . + +The energy along the Hugoniot is given by + +.. math:: + + E_H(\rho) = \frac{P_H \eta }{2 \rho_0} + E_0 . + +The temperature on the Hugoniot is hard to derive but with the help of Mathematica +it is +.. math:: + + T_H(\rho) = T_0 e^{\Gammma(\rho_0) \eta} + \frac{e^{\Gammma(\rho_0) \eta}}{2 C_V \rho_0} + \int_0^\eta e^{-\gamma(\rho_0) z} z^2 \frac{d}{dz} \left( \frac{P_H}{z}\right) dz + + = T_0 e^{\Gammma(\rho_0) \eta} + \frac{C_s^2}}{2 C_V s^2} + \left[\frac{- s \eta}{(1 - s \eta)^2} + \left( \frac{\Gamma(\rho_0)}{s} - 3 \right) + \left( e^{\Gammma(\rho_0) \eta} - \frac{1}{(1-s \eta)} + + e^{-\frac{\Gamma(\rho_0)}{s} (1-s \eta)} + \left( Ei(\frac{\Gamma(\rho_0)}{s}(1-s \eta))-Ei(\frac{\Gamma(\rho_0)}{s}) \right) + \left((\frac{\Gamma(\rho_0)}{s})^2 - 4 \frac{\Gamma(\rho_0)}{s} + 2 \right) \right] + + +The constructor for the ``MGUsup`` EOS has the signature + +.. code-block:: cpp + + MGUsup(const Real rho0, const Real T0, const Real Cs, const Real s, const Real G0, + const Real Cv0, const Real E0, const Real S0) + +where +``rho0`` is :math:`\rho_0`, ``T0`` is :math:`T_0`, +``Cs`` is :math:`C_s`, ``s`` is :math:`s`, +``G0`` is :math:`\Gamma(\rho_0)`, ``Cv0`` is :math:`C_v`, +``E0`` is :math:`E_0`, and ``S0`` is :math:`S_0`. + + JWL EOS `````````` diff --git a/singularity-eos/eos/eos.hpp b/singularity-eos/eos/eos.hpp index ac929d9d87..64d9f14dff 100644 --- a/singularity-eos/eos/eos.hpp +++ b/singularity-eos/eos/eos.hpp @@ -39,13 +39,13 @@ #include #include #include +#include #include #include #include #include #include #include -#include // Modifiers #include #include diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index b1ba0bb0e1..d6d056f99e 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -37,11 +37,10 @@ class MGUsup : public EosBase { public: MGUsup() = default; // Constructor - MGUsup(const Real rho0, const Real T0, const Real B0, const Real BP0, const Real A0, - const Real Cv0, const Real E0, const Real S0, const Real *expconsts) - : _rho0(rho0), _T0(T0), _B0(B0), _BP0(BP0), _A0(A0), _Cv0(Cv0), _E0(E0), _S0(S0) { + MGUsup(const Real rho0, const Real T0, const Real Cs, const Real s, const Real G0, + const Real Cv0, const Real E0, const Real S0) + : _rho0(rho0), _T0(T0), _Cs(Cs), _s(s), _G0(G0), _Cv0(Cv0), _E0(E0), _S0(S0) { CheckMGUsup(); - InitializeMGUsup(expconsts); } MGUsup GetOnDevice() { return *this; } @@ -55,7 +54,6 @@ class MGUsup : public EosBase { const Real rho, const Real sie, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real MinInternalEnergyFromDensity(const Real rho, Real *lambda = nullptr) const; - // Entropy added AEM Dec. 2022 PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( const Real rho, const Real temp, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( @@ -68,6 +66,10 @@ class MGUsup : public EosBase { const Real rho, const Real sie, Real *lambda = nullptr) const { return _Cv0; } + // added for testing AEM Dec 2023 + PORTABLE_INLINE_FUNCTION Real HugPressureFromDensity(const Real rho) const; + PORTABLE_INLINE_FUNCTION Real HugInternalEnergyFromDensity(const Real rho) const; + PORTABLE_INLINE_FUNCTION Real HugTemperatureFromDensity(const Real rho) const; // Thermal Bulk Modulus added AEM Dec 2022 PORTABLE_INLINE_FUNCTION Real TBulkModulusFromDensityTemperature( const Real rho, const Real temp, Real *lambda = nullptr) const; @@ -80,11 +82,11 @@ class MGUsup : public EosBase { const Real rho, const Real temp, Real *lambda = nullptr) const; PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( const Real rho, const Real temp, Real *lambda = nullptr) const { - return robust::ratio(_A0 * _B0, _Cv0 * rho); + return robust::ratio(_G0 * _rho0, rho); } PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const { - return robust::ratio(_A0 * _B0, _Cv0 * rho); + return robust::ratio(_G0 * _rho0, rho); } PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, const unsigned long output, @@ -105,12 +107,8 @@ class MGUsup : public EosBase { static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char st[]{"MGUsup Params: "}; - printf("%s rho0:%e T0:%e B0:%e BP0:%e\n A0:%e Cv0:%e E0:%e S0:%e\n" - "non-zero elements in d2tod40 array:\n", - st, _rho0, _T0, _B0, _BP0, _A0, _Cv0, _E0, _S0); - for (int i = 0; i < 39; i++) { - if (_d2tod40[i] > 0.0) printf("d%i:%e\t", i + 2, _d2tod40[i]); - } + printf("%s rho0:%e T0:%e Cs:%e s:%e\n G0:%e Cv0:%e E0:%e S0:%e\n", st, _rho0, _T0, + _Cs, _s, _G0, _Cv0, _E0, _S0); printf("\n\n"); } // Density/Energy from P/T not unique, if used will give error @@ -123,15 +121,9 @@ class MGUsup : public EosBase { private: static constexpr const unsigned long _preferred_input = - thermalqs::density | thermalqs::temperature; - Real _rho0, _T0, _B0, _BP0, _A0, _Cv0, _E0, _S0; - static constexpr const int PressureCoeffsd2tod40Size = 39; - static constexpr const int MGUsupInternalParametersSize = PressureCoeffsd2tod40Size + 4; - Real _VIP[MGUsupInternalParametersSize], _d2tod40[PressureCoeffsd2tod40Size]; + thermalqs::density | thermalqs::specific_internal_energy; + Real _rho0, _T0, _Cs, _s, _G0, _Cv0, _E0, _S0; void CheckMGUsup(); - void InitializeMGUsup(const Real *expcoeffs); - PORTABLE_INLINE_FUNCTION void MGUsup_F_DT_func(const Real rho, const Real T, - Real *output) const; }; inline void MGUsup::CheckMGUsup() { @@ -142,202 +134,188 @@ inline void MGUsup::CheckMGUsup() { if (_T0 < 0.0) { PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter T0 < 0"); } - if (_B0 < 0.0) { - PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter B0 < 0"); + if (_Cs < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter Cs < 0"); } - if (_BP0 < 1.0) { - PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter BP0 < 1"); + if (_s < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter s < 0"); } - if (_A0 < 0.0) { - PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter A0 < 0"); + if (_G0 < 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter G0 < 0"); } if (_Cv0 < 0.0) { PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter Cv0 < 0"); } } -inline void MGUsup::InitializeMGUsup(const Real *d2tod40input) { - - // The PressureCoeffsd2tod40Size (=39) allowed d2 to d40 coefficients - // for the pressure reference curve vs rho - // are seldom all used so did not want to crowd the argument list with them. - // Instead I ask the host code to send me a pointer to this array so that I can - // copy it here. Not used coeffs should be set to 0.0 (of course). - for (int ind = 0; ind < PressureCoeffsd2tod40Size; ind++) { - _d2tod40[ind] = d2tod40input[ind]; - } - // Put a couple of much used model parameter combinations and - // the energy coefficients f0 to f40 in an internal parameters array - // of size 2+41=MGUsupInternalParametersSize. - _VIP[0] = robust::ratio(_B0, _rho0) + - robust::ratio(_A0 * _A0 * _B0 * _B0 * _T0, - _rho0 * _rho0 * _Cv0); // sound speed squared - _VIP[1] = 3.0 / 2.0 * (_BP0 - 1.0); // exponent eta0 - // initializing with pressure coeffs to get energy coeffs into VIP - _VIP[2] = 1.0; // prefactor d0 - _VIP[3] = 0.0; // prefactor d1 - for (int ind = 4; ind < MGUsupInternalParametersSize; ind++) { //_d2tod40[0]=d2 - _VIP[ind] = _d2tod40[ind - 4]; // dn is in VIP[n+2] - } - for (int ind = MGUsupInternalParametersSize - 2; ind >= 2; - ind--) { // _VIP[42]=d40=f40 given,first calculated is _VIP[41]=f39 - _VIP[ind] = _VIP[ind] - (ind) / _VIP[1] * _VIP[ind + 1]; // prefactors f40 to f0 - } // _VIP[n+2]=fn, ind=n+2 +PORTABLE_INLINE_FUNCTION Real MGUsup::HugPressureFromDensity(Real rho) const { + Real eta = 1.0 - robust::ratio(_rho0, rho); + return _Cs * _Cs * _rho0 * robust::ratio(eta, (1.0 - _s * eta) * (1.0 - _s * eta)); } - -PORTABLE_INLINE_FUNCTION void MGUsup::MGUsup_F_DT_func(const Real rho, const Real T, - Real *output) const { - constexpr int pref0vp = -2, pref0vip = 2, maxind = MGUsupInternalParametersSize - 3; - Real sumP = 0.0, sumB = 0.0, sumE = 0.0; - - Real x = std::cbrt(robust::ratio(_rho0, rho)); /*rho-dependent*/ - Real x2inv = robust::ratio(1.0, x * x); - Real onemx = 1.0 - x; - Real etatimes1mx = _VIP[1] * onemx; - Real expetatimes1mx = exp(etatimes1mx); - -#pragma unroll - for (int ind = maxind; ind >= 2; ind--) { //_d2tod40[0]=d2 - sumP = _d2tod40[pref0vp + ind] + onemx * sumP; //_d2tod40[38]=d40 - sumB = _d2tod40[pref0vp + ind] * ind + onemx * sumB; //_d2tod40[-2+40]=d40 - sumE = _VIP[pref0vip + ind] + onemx * sumE; //_VIP[42]=f40 - } //_VIP[2]=f0 -#pragma unroll - for (int ind = 1; ind >= 0; ind--) { - sumE = _VIP[pref0vip + ind] + onemx * sumE; +PORTABLE_INLINE_FUNCTION Real MGUsup::HugInternalEnergyFromDensity(Real rho) const { + Real eta = 1.0 - robust::ratio(_rho0, rho); + return _E0 + robust::ratio(eta, 2.0 * _rho0) * HugPressureFromDensity(rho); +} +PORTABLE_INLINE_FUNCTION Real MGUsup::HugTemperatureFromDensity(Real rho) const { + int sumkmax = 20; + Real cutoff = 1.e-16; + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real f1 = 1.0 - _s * eta; + if (f1 <= 0.0) { + PORTABLE_ALWAYS_THROW_OR_ABORT("MGUsup model parameters s and rho0 together with rho " + "give a negative argument for a logarithm."); } - sumP = 1.0 + sumP * (onemx * onemx); - sumB = sumB * onemx; - sumE = _VIP[1] * onemx * sumE; - - /* T0 isotherm */ - Real energy = 9.0 * robust::ratio(_B0, _VIP[1] * _VIP[1] * _rho0) * - (_VIP[pref0vip] - expetatimes1mx * (_VIP[pref0vip] - sumE)) - - (_A0 * _B0) * (robust::ratio(_T0, _rho0) - robust::ratio(_T0, rho)); - Real pressure = 3.0 * _B0 * x2inv * onemx * expetatimes1mx * sumP; - Real temp = (1.0 + onemx * (_VIP[1] * x + 1.0)) * sumP + x * onemx * sumB; - temp = robust::ratio(_B0, _rho0) * x * expetatimes1mx * temp; - - /* Go to required temperature */ - energy = energy + _Cv0 * (T - _T0) + _E0; - pressure = pressure + (_A0 * _B0) * (T - _T0); - Real dpdrho = temp; - Real dpdt = _A0 * _B0; - Real dedt = _Cv0; - Real dedrho = robust::ratio(pressure - T * (_A0 * _B0), rho * rho); - Real entropy; - if (T < 0.0) { + Real G0os = robust::ratio(_G0, _s); + // sk, lk, mk, sum, and enough values may change in the loop + Real sk = -_G0 * eta; + Real lk = _G0; + Real mk = 0.0; + Real sum = sk; + int enough = 0; + Real temp; + Real pf = _Cs * _Cs / (2.0 * _Cv0 * _s * _s); + if (eta * eta < 1.e-8) { + temp = _T0 * exp(_G0 * eta) + + pf * (2.0 * std::log(f1) - _s * eta / (f1 * f1) * (3.0 * _s * eta - 2.0)); + } else { + for (int i = 0; ((i < sumkmax) && (enough == 0)); i++) { + mk = G0os / (i + 2) / (i + 2); + sk = (sk * f1 - lk * eta) * (i + 1) * mk; + lk = mk * (i + 1) * lk; + sum = sum + sk; + if (robust::ratio((sk * sk), (sum * sum)) < (cutoff * cutoff)) { + enough = 1; + } + } + if (enough == 0) { #ifndef NDEBUG - PORTABLE_WARN("Negative temperature input"); + PORTABLE_WARN("Hugoniot Temperature not converged"); #endif // NDEBUG - entropy = 0.0; - } else { - entropy = (_A0 * _B0) * (robust::ratio(1.0, rho) - robust::ratio(1, _rho0)) + - _Cv0 * std::log(robust::ratio(T, _T0)) + _S0; + } + // printf("sum=%e\n",sum); + temp = _T0 - pf * ((G0os - 3.0) + exp(-G0os) * (G0os * G0os - 4.0 * G0os + 2.0) * + (std::log(f1) + sum)); + temp = temp * exp(_G0 * eta); + temp = temp - pf / f1 * ((4.0 - G0os) - 1.0 / f1); } - Real soundspeed = - std::max(0.0, dpdrho + T * robust::ratio(dpdt * dpdt, rho * rho * _Cv0)); - soundspeed = std::sqrt(soundspeed); - - output[0] = energy; - output[1] = pressure; - output[2] = dpdrho; - output[3] = dpdt; - output[4] = dedt; - output[5] = dedrho; - output[6] = entropy; - output[7] = soundspeed; - - return; + return temp; } + PORTABLE_INLINE_FUNCTION Real MGUsup::InternalEnergyFromDensityTemperature( const Real rho, const Real temp, Real *lambda) const { - Real output[8]; - MGUsup_F_DT_func(rho, temp, output); - return output[0]; + Real value = + HugInternalEnergyFromDensity(rho) + _Cv0 * (temp - HugTemperatureFromDensity(rho)); + return value; } PORTABLE_INLINE_FUNCTION Real MGUsup::PressureFromDensityTemperature(const Real rho, - const Real temp, - Real *lambda) const { - Real output[8]; - MGUsup_F_DT_func(rho, temp, output); - return output[1]; + const Real temp, + Real *lambda) const { + Real value = HugPressureFromDensity(rho) + + _G0 * _rho0 * _Cv0 * (temp - HugTemperatureFromDensity(rho)); + return value; } PORTABLE_INLINE_FUNCTION Real MGUsup::EntropyFromDensityTemperature(const Real rho, - const Real temp, - Real *lambda) const { - Real output[8]; - MGUsup_F_DT_func(rho, temp, output); - return output[6]; + const Real temp, + Real *lambda) const { + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real value = _S0 - _G0 * _Cv0 * eta + _Cv0 * std::log(temp / _T0); + return value; } PORTABLE_INLINE_FUNCTION Real MGUsup::TExpansionCoeffFromDensityTemperature( const Real rho, const Real temp, Real *lambda) const { - Real output[8]; - MGUsup_F_DT_func(rho, temp, output); - return robust::ratio(output[3], output[2] * rho); + Real value = + robust::ratio(_Cv0 * _rho0 * _G0, TBulkModulusFromDensityTemperature(rho, temp)); + return value; } PORTABLE_INLINE_FUNCTION Real MGUsup::TBulkModulusFromDensityTemperature( const Real rho, const Real temp, Real *lambda) const { - Real output[8]; - MGUsup_F_DT_func(rho, temp, output); - return output[2] * rho; + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real value = robust::ratio((1.0 + _s * eta - _G0 * _s * eta * eta), (1.0 - _s * eta)); + if (eta == 0.0) { + value = _Cs * _Cs * _rho0; + } else { + value = value * robust::ratio(HugPressureFromDensity(rho), eta); + } + value = value - _G0 * _G0 * _Cv0 * _rho0 * HugTemperatureFromDensity(rho); + value = robust::ratio(_rho0, rho) * value; + return value; } PORTABLE_INLINE_FUNCTION Real MGUsup::BulkModulusFromDensityTemperature( const Real rho, const Real temp, Real *lambda) const { - Real output[8]; - MGUsup_F_DT_func(rho, temp, output); - return output[7] * output[7] * rho; + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real value = + robust::ratio((1.0 + _s * eta - _G0 * _s * eta * eta), eta * (1.0 - _s * eta)); + if (eta == 0.0) { + value = _Cs * _Cs * _rho0; + } else { + value = value * robust::ratio(HugPressureFromDensity(rho), eta); + } + value = value - _G0 * _G0 * _Cv0 * _rho0 * (temp - HugTemperatureFromDensity(rho)); + value = robust::ratio(_rho0, rho) * value; + return value; } PORTABLE_INLINE_FUNCTION Real MGUsup::TemperatureFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda) const { - Real Tref; - Real output[8]; - Tref = _T0; - MGUsup_F_DT_func(rho, Tref, output); - return robust::ratio(sie - output[0], _Cv0) + Tref; + Real value = + (sie - HugInternalEnergyFromDensity(rho)) / _Cv0 + HugTemperatureFromDensity(rho); + if (value < 0.0) { +#ifndef NDEBUG + PORTABLE_WARN("Negative temperature"); +#endif // NDEBUG + value = 1.e-12; + } + return value; } PORTABLE_INLINE_FUNCTION Real MGUsup::PressureFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda) const { - Real temp; - Real output[8]; - temp = TemperatureFromDensityInternalEnergy(rho, sie); - MGUsup_F_DT_func(rho, temp, output); - return output[1]; + Real value = HugPressureFromDensity(rho) + + _rho0 * _G0 * (sie - HugInternalEnergyFromDensity(rho)); + return value; } PORTABLE_INLINE_FUNCTION Real MGUsup::MinInternalEnergyFromDensity(const Real rho, - Real *lambda) const { + Real *lambda) const { MinInternalEnergyIsNotEnabled("MGUsup"); return 0.0; } PORTABLE_INLINE_FUNCTION Real MGUsup::EntropyFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda) const { - Real temp; - Real output[8]; - temp = TemperatureFromDensityInternalEnergy(rho, sie); - MGUsup_F_DT_func(rho, temp, output); - return output[6]; + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real value = std::log(TemperatureFromDensityInternalEnergy(rho, sie) / _T0); + value = _S0 - _G0 * _Cv0 * eta + _Cv0 * value; + if (value < 0.0) { +#ifndef NDEBUG + PORTABLE_WARN("Negative entropy"); +#endif // NDEBUG + value = 1.e-12; + } + return value; } PORTABLE_INLINE_FUNCTION Real MGUsup::BulkModulusFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda) const { - Real temp; - Real output[8]; - temp = TemperatureFromDensityInternalEnergy(rho, sie); - MGUsup_F_DT_func(rho, temp, output); - return output[7] * output[7] * rho; + Real eta = 1.0 - robust::ratio(_rho0, rho); + Real value = robust::ratio((1.0 + _s * eta - _G0 * _s * eta * eta), (1.0 - _s * eta)); + if (eta == 0.0) { + value = _Cs * _Cs * _rho0; + } else { + value = value * robust::ratio(HugPressureFromDensity(rho), eta); + } + value = value + _G0 * _G0 * _rho0 * (sie - HugInternalEnergyFromDensity(rho)); + value = robust::ratio(_rho0, rho) * value; + return value; } // AEM: Give error since function is not well defined PORTABLE_INLINE_FUNCTION void MGUsup::DensityEnergyFromPressureTemperature(const Real press, const Real temp, - Real *lambda, Real &rho, Real &sie) const { + Real *lambda, Real &rho, Real &sie) const { EOS_ERROR("MGUsup::DensityEnergyFromPressureTemperature: " "Not implemented.\n"); } // AEM: We should add entropy and Gruneissen parameters here so that it is complete // If we add also alpha and BT, those should also be in here. PORTABLE_INLINE_FUNCTION void MGUsup::FillEos(Real &rho, Real &temp, Real &sie, - Real &press, Real &cv, Real &bmod, - const unsigned long output, - Real *lambda) const { + Real &press, Real &cv, Real &bmod, + const unsigned long output, + Real *lambda) const { const unsigned long input = ~output; // everything that is not output is input if (thermalqs::density & output) { EOS_ERROR("MGUsup FillEos: Density is required input.\n"); @@ -346,23 +324,24 @@ PORTABLE_INLINE_FUNCTION void MGUsup::FillEos(Real &rho, Real &temp, Real &sie, EOS_ERROR("MGUsup FillEos: Density and Internal Energy or Density and Temperature " "are required input parameters.\n"); } - if (thermalqs::specific_internal_energy & input) { - temp = TemperatureFromDensityInternalEnergy(rho, sie); + if (thermalqs::temperature & input) { + sie = InternalEnergyFromDensityTemperature(rho, temp); } - Real Vout[8]; - MGUsup_F_DT_func(rho, temp, Vout); - if (output & thermalqs::temperature) temp = temp; - if (output & thermalqs::specific_internal_energy) sie = Vout[0]; - if (output & thermalqs::pressure) press = Vout[1]; - if (output & thermalqs::specific_heat) cv = Vout[4]; - if (output & thermalqs::bulk_modulus) bmod = Vout[7] * Vout[7] * rho; + if (output & thermalqs::temperature) + temp = TemperatureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::specific_internal_energy) sie = sie; + if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::specific_heat) + cv = SpecificHeatFromDensityInternalEnergy(rho, sie); + if (output & thermalqs::bulk_modulus) + bmod = BulkModulusFromDensityInternalEnergy(rho, sie); } // TODO(JMM): pre-cache these rather than recomputing them each time PORTABLE_INLINE_FUNCTION void MGUsup::ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, - Real &cv, Real &bmod, Real &dpde, Real &dvdt, - Real *lambda) const { + Real &cv, Real &bmod, Real &dpde, Real &dvdt, + Real *lambda) const { // AEM: Added all variables I think should be output eventually Real tbmod; // Real entropy, alpha, Gamma; @@ -373,13 +352,13 @@ void MGUsup::ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &pres press = PressureFromDensityTemperature(rho, temp, lambda); // entropy = _S0; cv = _Cv0; - tbmod = _B0; + tbmod = _Cs * _Cs * _rho0 - _G0 * _G0 * _Cv0 * _rho0 * _T0; // alpha = _A0; - bmod = BulkModulusFromDensityTemperature(rho, temp, lambda); - // Gamma = robust::ratio(_A0 * _B0, _Cv0 * _rho0); + bmod = _Cs * _Cs * _rho0; + // Gamma = _G0; // AEM: I suggest taking the two following away. - dpde = robust::ratio(_A0 * tbmod, _Cv0); - dvdt = robust::ratio(_A0, _rho0); + dpde = _G0 * _rho0; + dvdt = robust::ratio(-1.0, _T0 * _G0 * _rho0); } } // namespace singularity diff --git a/test/test_eos_mgusup.cpp b/test/test_eos_mgusup.cpp index b43194f3d7..273197c668 100644 --- a/test/test_eos_mgusup.cpp +++ b/test/test_eos_mgusup.cpp @@ -35,22 +35,19 @@ SCENARIO("MGUsup EOS rho sie", "[VectorEOS][MGUsupEOS]") { GIVEN("Parameters for a MGUsup EOS") { // Unit conversions constexpr Real Mbcc_per_g = 1e12; - // MGUsup parameters for copper - constexpr Real rho0 = 8.93; + // MGUsup parameters for tin beta phase + constexpr Real rho0 = 7.285; constexpr Real T0 = 298.0; - constexpr Real B0 = 1.3448466 * Mbcc_per_g; - constexpr Real BP0 = 4.956; - constexpr Real A0 = 5.19245e-05; - constexpr Real Cv0 = 0.383e-05 * Mbcc_per_g; - constexpr Real E0 = 0.0; - constexpr Real S0 = 5.05e-04 * Mbcc_per_g; - constexpr Real d2to40[39] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; + constexpr Real Cs = 2766.0e2; + constexpr Real s = 1.5344; + constexpr Real G0 = 2.4659; + constexpr Real Cv0 = 0.2149e-05 * Mbcc_per_g; + constexpr Real E0 = 0.658e-03 * Mbcc_per_g; + constexpr Real S0 = 0.4419e-05 * Mbcc_per_g; // Create the EOS - EOS host_eos = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40); + EOS host_eos = MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0); EOS eos = host_eos.GetOnDevice(); - MGUsup host_eos2 = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40); + MGUsup host_eos2 = MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0); MGUsup eos2 = host_eos2.GetOnDevice(); eos.PrintParams(); @@ -74,14 +71,14 @@ SCENARIO("MGUsup EOS rho sie", "[VectorEOS][MGUsupEOS]") { // Populate the input arrays portableFor( "Initialize density and energy", 0, 1, PORTABLE_LAMBDA(int i) { - v_density[0] = 8.0; - v_density[1] = 8.5; + v_density[0] = 7.0; + v_density[1] = 8.0; v_density[2] = 9.0; - v_density[3] = 8.93; - v_energy[0] = 2.e8; - v_energy[1] = 4.6e9; - v_energy[2] = 4.7e9; - v_energy[3] = 0.0; + v_density[3] = 7.285; + v_energy[0] = 4.6e8; + v_energy[1] = 1.146e9; + v_energy[2] = 3.28e9; + v_energy[3] = 0.658e9; }); #ifdef PORTABILITY_STRATEGY_KOKKOS // Create host-side mirrors of the inputs and copy the inputs. These are @@ -94,29 +91,40 @@ SCENARIO("MGUsup EOS rho sie", "[VectorEOS][MGUsupEOS]") { // Gold standard values for a subset of lookups constexpr std::array temperature_true{ - 6.638294350641033e+1, 1.42266691572769e+03, 1.52867838544081e+03, T0}; + 1.502141159804009e+02, 4.26645103516999e+02, 7.08167099601959e+02, T0}; + constexpr std::array hugtemperature_true{ + 2.684894520258222e+02, 3.90542065829047e+02, 7.78960970777667e+02, T0}; constexpr std::array pressure_true{ - -1.283459624233147e+11, 1.98538351697004e+10, 9.66446220551959e+10, 0.}; + -2.466829656715215e+10, 6.82999362849669e+10, 2.09379236091689e+11, 0.}; + constexpr std::array hugpressure_true{ + -2.010229955618467e+10, 6.690618533331688e+10, 2.121122201185450e+11, 0.}; + constexpr std::array hugenergy_true{ + 7.141736971616102e+8, 1.068414572008593e+09, 3.432136029156598e+09, E0}; constexpr std::array entropy_true{ - 5.0015771521255749e+08, 5.1138262492866594e+08, 5.1120147992457777e+08, S0}; + 3.1626206459885668e+06, 4.7165703738323096e+06, 5.2693499546078807e+06, S0}; constexpr std::array cv_true{Cv0, Cv0, Cv0, Cv0}; constexpr std::array bulkmodulus_true{ - 7.44411028807209e+11, 1.254880120063067e+12, 1.613822819346433e+12, - (B0 + B0 * B0 * A0 * A0 / Cv0 / rho0 * T0)}; - constexpr std::array gruneisen_true{ - B0 * A0 / Cv0 / 8.0, B0 * A0 / Cv0 / 8.5, B0 * A0 / Cv0 / 9.0, - B0 * A0 / Cv0 / rho0}; + 4.38665321162636e+11, 8.776332438691490e+11, 1.465222098683647e+12, + Cs * Cs * rho0}; + constexpr std::array gruneisen_true{G0 * rho0 / 7.0, G0 * rho0 / 8.0, + G0 * rho0 / 9.0, G0}; #ifdef PORTABILITY_STRATEGY_KOKKOS // Create device views for outputs and mirror those views on the host Kokkos::View v_temperature("Temperature"); + Kokkos::View v_hugtemperature("HugTemperature"); Kokkos::View v_pressure("Pressure"); + Kokkos::View v_hugpressure("HugPressure"); + Kokkos::View v_hugenergy("HugEnergy"); Kokkos::View v_entropy("Entropy"); Kokkos::View v_cv("Cv"); Kokkos::View v_bulkmodulus("bmod"); Kokkos::View v_gruneisen("Gamma"); auto h_temperature = Kokkos::create_mirror_view(v_temperature); + auto h_hugtemperature = Kokkos::create_mirror_view(v_hugtemperature); auto h_pressure = Kokkos::create_mirror_view(v_pressure); + auto h_hugpressure = Kokkos::create_mirror_view(v_hugpressure); + auto h_hugenergy = Kokkos::create_mirror_view(v_hugenergy); auto h_entropy = Kokkos::create_mirror_view(v_entropy); auto h_cv = Kokkos::create_mirror_view(v_cv); auto h_bulkmodulus = Kokkos::create_mirror_view(v_bulkmodulus); @@ -125,20 +133,74 @@ SCENARIO("MGUsup EOS rho sie", "[VectorEOS][MGUsupEOS]") { // Create arrays for the outputs and then pointers to those arrays that // will be passed to the functions in place of the Kokkos views std::array h_temperature; + std::array h_hugtemperature; std::array h_pressure; + std::array h_hugpressure; + std::array h_hugenergy; std::array h_entropy; std::array h_cv; std::array h_bulkmodulus; std::array h_gruneisen; // Just alias the existing pointers auto v_temperature = h_temperature.data(); + auto v_hugtemperature = h_hugtemperature.data(); auto v_pressure = h_pressure.data(); + auto v_hugpressure = h_hugpressure.data(); + auto v_hugenergy = h_hugenergy.data(); auto v_entropy = h_entropy.data(); auto v_cv = h_cv.data(); auto v_bulkmodulus = h_bulkmodulus.data(); auto v_gruneisen = h_gruneisen.data(); #endif // PORTABILITY_STRATEGY_KOKKOS + // temporary + + WHEN("A PH(rho) lookup is performed") { + portableFor( + "Test Hugoniot pressure", 0, num, PORTABLE_LAMBDA(const int i) { + v_hugpressure[i] = eos2.HugPressureFromDensity(v_density[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_hugpressure, v_hugpressure); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned PH(rho) should be equal to the true value") { + array_compare(num, density, energy, h_hugpressure, hugpressure_true, "Density", + "Energy"); + } + } + + WHEN("A EH(rho) lookup is performed") { + portableFor( + "Test Hugoniot internal energy", 0, num, PORTABLE_LAMBDA(const int i) { + v_hugenergy[i] = eos2.HugInternalEnergyFromDensity(v_density[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_hugenergy, v_hugenergy); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned EH(rho) should be equal to the true value") { + array_compare(num, density, energy, h_hugenergy, hugenergy_true, "Density", + "Energy"); + } + } + + WHEN("A TH(rho) lookup is performed") { + portableFor( + "Test Hugoniot temperature", 0, num, PORTABLE_LAMBDA(const int i) { + v_hugtemperature[i] = eos2.HugTemperatureFromDensity(v_density[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_hugtemperature, v_hugtemperature); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned TH(rho) should be equal to the true value") { + array_compare(num, density, energy, h_hugtemperature, hugtemperature_true, + "Density", "Energy"); + } + } + // end temporary + WHEN("A T(rho, e) lookup is performed") { eos.TemperatureFromDensityInternalEnergy(v_density, v_energy, v_temperature, num); #ifdef PORTABILITY_STRATEGY_KOKKOS @@ -222,21 +284,18 @@ SCENARIO("MGUsup EOS rho T", "[VectorEOS][MGUsupEOS]") { // Unit conversions constexpr Real Mbcc_per_g = 1e12; // MGUsup parameters for copper - constexpr Real rho0 = 8.93; + constexpr Real rho0 = 7.285; constexpr Real T0 = 298.0; - constexpr Real B0 = 1.3448466 * Mbcc_per_g; - constexpr Real BP0 = 4.956; - constexpr Real A0 = 5.19245e-05; - constexpr Real Cv0 = 0.383e-05 * Mbcc_per_g; - constexpr Real E0 = 0.0; - constexpr Real S0 = 5.05e-04 * Mbcc_per_g; - constexpr Real d2to40[39] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; + constexpr Real Cs = 2766.0e2; + constexpr Real s = 1.5344; + constexpr Real G0 = 2.4659; + constexpr Real Cv0 = 0.2149e-05 * Mbcc_per_g; + constexpr Real E0 = 0.658e-03 * Mbcc_per_g; + constexpr Real S0 = 0.4419e-05 * Mbcc_per_g; // Create the EOS - EOS host_eos = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40); + EOS host_eos = MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0); EOS eos = host_eos.GetOnDevice(); - MGUsup host_eos2 = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40); + MGUsup host_eos2 = MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0); MGUsup eos2 = host_eos2.GetOnDevice(); eos.PrintParams(); @@ -259,13 +318,13 @@ SCENARIO("MGUsup EOS rho T", "[VectorEOS][MGUsupEOS]") { // Populate the input arrays portableFor( "Initialize density and temperature", 0, 1, PORTABLE_LAMBDA(int i) { - v_density[0] = 8.0; - v_density[1] = 8.5; + v_density[0] = 7.0; + v_density[1] = 8.0; v_density[2] = 9.0; - v_density[3] = 8.93; - v_temperature[0] = 66.383; - v_temperature[1] = 1422.7; - v_temperature[2] = 1528.7; + v_density[3] = 7.285; + v_temperature[0] = 150.215; + v_temperature[1] = 426.645; + v_temperature[2] = 708.165; v_temperature[3] = 298.0; }); #ifdef PORTABILITY_STRATEGY_KOKKOS @@ -278,19 +337,20 @@ SCENARIO("MGUsup EOS rho T", "[VectorEOS][MGUsupEOS]") { #endif // PORTABILITY_STRATEGY_KOKKOS // Gold standard values for a subset of lookups - constexpr std::array sie_true{ - 2.0000021637044835e+8, 4.6001267127629471e+09, 4.7000827837617149e+09, 0.0}; + constexpr std::array sie_true{4.6000189975811839e+08, + 1.1459997775419691e+09, + 3.2799954879553909e+09, 6.58e+08}; constexpr std::array pressure_true{ - -1.283459584783398e+11, 1.98561454605571e+10, 9.66461314103968e+10, 0.}; + -2.466826243974248e+10, 6.82999322887127e+10, 2.09379155036952e+11, 0.}; constexpr std::array entropy_true{ - 5.0015771847198445e+08, 5.1138271399469268e+08, 5.1120153407800680e+08, S0}; + 3.1626332929526418e+06, 4.7165698524198849e+06, 5.2693435831578374e+06, S0}; constexpr std::array cv_true{Cv0, Cv0, Cv0, Cv0}; constexpr std::array tbulkmodulus_true{ - 7.3384631127398950e+11, 1.0417839336794381e+12, 1.3975684023296028e+12, B0}; + 4.2378339466579810e+11, 8.4064844786200464e+11, 1.4106538914691243e+12, + (Cs * Cs - G0 * G0 * T0 * Cv0) * rho0}; constexpr std::array alpha_true{ - 9.5156828083623124e-05, 6.7029721830195901e-05, 4.9965702691402983e-05, A0}; - // B0 * A0 / tbulkmodulus_true[0], B0 * A0 / tbulkmodulus_true[1], B0 * A0 / - // tbulkmodulus_true[2] + 9.1095620143267597e-05, 4.5922657969193220e-05, 2.7366607342141916e-05, + G0 * Cv0 / (Cs * Cs - G0 * G0 * T0 * Cv0)}; #ifdef PORTABILITY_STRATEGY_KOKKOS // Create device views for outputs and mirror those views on the host @@ -428,52 +488,39 @@ SCENARIO("MGUsup EOS SetUp", "[VectorEOS][MGUsupEOS]") { Real Cv0 = 0.383e-05 * Mbcc_per_g; constexpr Real E0 = 0.0; constexpr Real S0 = 5.05e-04 * Mbcc_per_g; - constexpr Real d2to40[39] = {0., 0., 0.0000000000001, - 0., 0., 0., - 0., 0., 0., - 0., 0., 0., - 0., 0., 0., - 0., 0., 0., - 0., 0., 0., - 0., 0., 0., - 0., 0., 0., - 0., 0., 0., - 0., 0., 0., - 0., 1.0e-26, 0., - 0., 0., 0.}; // Create the EOS - EOS host_eos = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40); + EOS host_eos = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0); EOS eos = host_eos.GetOnDevice(); eos.Finalize(); WHEN("Faulty/not set parameter rho0 is given") { rho0 = -1.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); THEN("An error message should be written out") {} } WHEN("Faulty/not set parameter T0 is given") { T0 = -1.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); THEN("An error message should be written out") {} } WHEN("Faulty/not set parameter B0 is given") { B0 = -1.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); THEN("An error message should be written out") {} } WHEN("Faulty/not set parameter BP0 is given") { - BP0 = 0.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + BP0 = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); THEN("An error message should be written out") {} } WHEN("Faulty/not set parameter Cv0 is given") { Cv0 = -1.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); THEN("An error message should be written out") {} } WHEN("Faulty/not set parameter A0 is given") { A0 = -10000000.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0, d2to40)); + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); THEN("An error message should be written out") {} } } From 80ba8189d7a1f64c3f594f7423143d682835627e Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Wed, 3 Jan 2024 11:24:33 -0700 Subject: [PATCH 012/405] Updates to PR#331. --- CHANGELOG.md | 2 ++ doc/sphinx/src/models.rst | 30 +++++++++++++++++++++++++----- singularity-eos/eos/eos.hpp | 2 +- singularity-eos/eos/eos_mgusup.hpp | 2 +- test/test_eos_mgusup.cpp | 2 +- 5 files changed, 30 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1f91381416..23369ea421 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,8 @@ ## Current develop ### Fixed (Repair bugs, etc) +### Added (new features/APIs/variables/...) +- [[PR331]](https://github.com/lanl/singularity-eos/pull/331) Included code and documentation for a full, temperature consistent, Mie-Gruneisen EOS based on a linear Us-up relation. ### Added (new features/APIs/variables/...) - [[PR326]](https://github.com/lanl/singularity-eos/pull/326) Document how to do a release diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 3e5819aca3..15294b25c5 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -580,6 +580,9 @@ Given the inconsisetency in the temperature, we have made the choice **not** to expose the entropy for this EOS. **Requesting an entropy value will result in an error.** +If a linear :math:`U_s`-:math:`u_p` relation is enough for your problem, we recommend using the MGUsup +EOS described below. It is a complete EOS with consistent temperature. + Given a reference density, :math:`\rho_0`, we first parameterize the EOS using :math:`\eta` as a measure of compression given by @@ -780,26 +783,42 @@ heat capacity is assumed so that Note the difference from the Gruneisen EOS described above. We still use a constant :math:`C_V`, and it is usually taken at the reference temperature, but -we now extrapolate from the temperature on the Hugoniot, :math: `T_H(\rho)`, and not +we now extrapolate from the temperature on the Hugoniot, :math:`T_H(\rho)`, and not from the reference temperature, :math:`T_0`. -With this consistent temperature we can derive an entropy in a similar way as for the Vinet EOS, +With this consistent temperature we can derive an entropy in a similar way as for the Vinet EOS. Using +thermodynamic derivatives we can show that + +.. math:: + + \Gamma \rho = \frac{\alpha B_T}{C_V} , + +and we arrive at + .. math:: S(\rho,T) = S_0 - \Gamma(\rho_0)C_V \eta + {C_V} \ln \frac{T}{T_{ref}} , where :math:`\eta` is a measure of compression given by + .. math:: \eta = 1 - \frac{\rho_0}{\rho}. -This is convenient because :math:`eta = 0` when :math:`\rho = \rho_0`, +This is convenient because :math:`\eta = 0` when :math:`\rho = \rho_0`, :math:`\eta = 1` at the infinite density limit, and :math:`\eta = -\infty` at the zero density limit. The pressure, energy, and temperature, on the Hugoniot are derived from the -shock jump conditions assuming a linear :math:`U_s`-:math:`u_p` relation, +shock jump conditions, + +.. math:: + \rho_0 U_s = \rho (U_s - u_p) + P_H = \rho_0 U_s u_p , + +assuming a linear :math:`U_s`-:math:`u_p` relation, + .. math:: U_s = C_s + s u_p . @@ -825,6 +844,7 @@ The energy along the Hugoniot is given by The temperature on the Hugoniot is hard to derive but with the help of Mathematica it is + .. math:: T_H(\rho) = T_0 e^{\Gammma(\rho_0) \eta} + \frac{e^{\Gammma(\rho_0) \eta}}{2 C_V \rho_0} @@ -848,7 +868,7 @@ The constructor for the ``MGUsup`` EOS has the signature where ``rho0`` is :math:`\rho_0`, ``T0`` is :math:`T_0`, ``Cs`` is :math:`C_s`, ``s`` is :math:`s`, -``G0`` is :math:`\Gamma(\rho_0)`, ``Cv0`` is :math:`C_v`, +``G0`` is :math:`\Gamma(\rho_0)`, ``Cv0`` is :math:`C_V`, ``E0`` is :math:`E_0`, and ``S0`` is :math:`S_0`. diff --git a/singularity-eos/eos/eos.hpp b/singularity-eos/eos/eos.hpp index 64d9f14dff..f5d45382f6 100644 --- a/singularity-eos/eos/eos.hpp +++ b/singularity-eos/eos/eos.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2024. 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 diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index d6d056f99e..b256726ede 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2024. 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 diff --git a/test/test_eos_mgusup.cpp b/test/test_eos_mgusup.cpp index 273197c668..380bc5da30 100644 --- a/test/test_eos_mgusup.cpp +++ b/test/test_eos_mgusup.cpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2024. 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 51707603bb62b1f11b6640b042b01e870afe0bb5 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Wed, 3 Jan 2024 12:10:32 -0700 Subject: [PATCH 013/405] Update 2 to #331. --- doc/sphinx/src/models.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 15294b25c5..c900cd38f2 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -797,7 +797,7 @@ and we arrive at .. math:: - S(\rho,T) = S_0 - \Gamma(\rho_0)C_V \eta + {C_V} \ln \frac{T}{T_{ref}} , + S(\rho,T) = S_0 - \Gamma(\rho_0)C_V \eta + {C_V} \ln \frac{T}{T_0} , where :math:`\eta` is a measure of compression given by @@ -814,7 +814,7 @@ The pressure, energy, and temperature, on the Hugoniot are derived from the shock jump conditions, .. math:: - \rho_0 U_s = \rho (U_s - u_p) + \rho_0 U_s = \rho (U_s - u_p) \\ P_H = \rho_0 U_s u_p , assuming a linear :math:`U_s`-:math:`u_p` relation, @@ -830,7 +830,7 @@ so only this :math:`s` parameter is needed. The units for :math:`C_s` is velocit :math:`s` is unitless. Note that the parameter :math:`s` is related to the fundamental derivative of shock physics as shown by `Wills `_. -The reference pressure along the Hugoniot is determined by +Solving the jump equations above gives that the reference pressure along the Hugoniot is determined by .. math:: @@ -848,7 +848,7 @@ it is .. math:: T_H(\rho) = T_0 e^{\Gammma(\rho_0) \eta} + \frac{e^{\Gammma(\rho_0) \eta}}{2 C_V \rho_0} - \int_0^\eta e^{-\gamma(\rho_0) z} z^2 \frac{d}{dz} \left( \frac{P_H}{z}\right) dz + \int_0^\eta e^{-\gamma(\rho_0) z} z^2 \frac{d}{dz} \left( \frac{P_H}{z}\right) dz \\ = T_0 e^{\Gammma(\rho_0) \eta} + \frac{C_s^2}}{2 C_V s^2} \left[\frac{- s \eta}{(1 - s \eta)^2} + \left( \frac{\Gamma(\rho_0)}{s} - 3 \right) From aee8452eb60ba90cc8b4757699e06e8f32119ed0 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Wed, 3 Jan 2024 16:21:24 -0700 Subject: [PATCH 014/405] Update 3 of #331. --- doc/sphinx/src/models.rst | 5 ++-- singularity-eos/eos/eos_mgusup.hpp | 6 ++-- test/test_eos_mgusup.cpp | 44 +++++++++++++++--------------- 3 files changed, 27 insertions(+), 28 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index c900cd38f2..6458601b78 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -815,7 +815,7 @@ shock jump conditions, .. math:: \rho_0 U_s = \rho (U_s - u_p) \\ - P_H = \rho_0 U_s u_p , + P_H = \rho_0 U_s u_p , assuming a linear :math:`U_s`-:math:`u_p` relation, @@ -849,10 +849,9 @@ it is T_H(\rho) = T_0 e^{\Gammma(\rho_0) \eta} + \frac{e^{\Gammma(\rho_0) \eta}}{2 C_V \rho_0} \int_0^\eta e^{-\gamma(\rho_0) z} z^2 \frac{d}{dz} \left( \frac{P_H}{z}\right) dz \\ - = T_0 e^{\Gammma(\rho_0) \eta} + \frac{C_s^2}}{2 C_V s^2} \left[\frac{- s \eta}{(1 - s \eta)^2} + \left( \frac{\Gamma(\rho_0)}{s} - 3 \right) - \left( e^{\Gammma(\rho_0) \eta} - \frac{1}{(1-s \eta)} + \left( e^{\Gammma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right) + e^{-\frac{\Gamma(\rho_0)}{s} (1-s \eta)} \left( Ei(\frac{\Gamma(\rho_0)}{s}(1-s \eta))-Ei(\frac{\Gamma(\rho_0)}{s}) \right) \left((\frac{\Gamma(\rho_0)}{s})^2 - 4 \frac{\Gamma(\rho_0)}{s} + 2 \right) \right] diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index b256726ede..58b560268f 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -40,7 +40,7 @@ class MGUsup : public EosBase { MGUsup(const Real rho0, const Real T0, const Real Cs, const Real s, const Real G0, const Real Cv0, const Real E0, const Real S0) : _rho0(rho0), _T0(T0), _Cs(Cs), _s(s), _G0(G0), _Cv0(Cv0), _E0(E0), _S0(S0) { - CheckMGUsup(); + _CheckMGUsup(); } MGUsup GetOnDevice() { return *this; } @@ -123,10 +123,10 @@ class MGUsup : public EosBase { static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; Real _rho0, _T0, _Cs, _s, _G0, _Cv0, _E0, _S0; - void CheckMGUsup(); + void _CheckMGUsup(); }; -inline void MGUsup::CheckMGUsup() { +inline void MGUsup::_CheckMGUsup() { if (_rho0 < 0.0) { PORTABLE_ALWAYS_THROW_OR_ABORT("Required MGUsup model parameter rho0 < 0"); diff --git a/test/test_eos_mgusup.cpp b/test/test_eos_mgusup.cpp index 380bc5da30..df35715eac 100644 --- a/test/test_eos_mgusup.cpp +++ b/test/test_eos_mgusup.cpp @@ -480,47 +480,47 @@ SCENARIO("MGUsup EOS SetUp", "[VectorEOS][MGUsupEOS]") { // Unit conversions constexpr Real Mbcc_per_g = 1e12; // MGUsup parameters for copper - Real rho0 = 8.93; + Real rho0 = 7.285; Real T0 = 298.0; - Real B0 = 1.3448466 * Mbcc_per_g; - Real BP0 = 4.956; - Real A0 = 5.19245e-05; - Real Cv0 = 0.383e-05 * Mbcc_per_g; - constexpr Real E0 = 0.0; - constexpr Real S0 = 5.05e-04 * Mbcc_per_g; + Real Cs = 2766.0e2; + Real s = 1.5344; + Real G0 = 2.4659; + Real Cv0 = 0.2149e-05 * Mbcc_per_g; + constexpr Real E0 = 0.658e-03 * Mbcc_per_g; + constexpr Real S0 = 0.4419e-05 * Mbcc_per_g; // Create the EOS - EOS host_eos = MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0); + EOS host_eos = MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0); EOS eos = host_eos.GetOnDevice(); eos.Finalize(); WHEN("Faulty/not set parameter rho0 is given") { rho0 = -1.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); THEN("An error message should be written out") {} } WHEN("Faulty/not set parameter T0 is given") { T0 = -1.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); THEN("An error message should be written out") {} } - WHEN("Faulty/not set parameter B0 is given") { - B0 = -1.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); + WHEN("Faulty/not set parameter Cs is given") { + Cs = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); THEN("An error message should be written out") {} } - WHEN("Faulty/not set parameter BP0 is given") { - BP0 = -1.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); + WHEN("Faulty/not set parameter s is given") { + s = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); THEN("An error message should be written out") {} } - WHEN("Faulty/not set parameter Cv0 is given") { - Cv0 = -1.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); + WHEN("Faulty/not set parameter G0 is given") { + G0 = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); THEN("An error message should be written out") {} } - WHEN("Faulty/not set parameter A0 is given") { - A0 = -10000000.0; - REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, B0, BP0, A0, Cv0, E0, S0)); + WHEN("Faulty/not set parameter Cv0 is given") { + Cv0 = -1.0; + REQUIRE_MAYBE_THROWS(MGUsup(rho0, T0, Cs, s, G0, Cv0, E0, S0)); THEN("An error message should be written out") {} } } From d741fc4bc0646421a1fbf071e7f53a6f51d41d32 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Wed, 3 Jan 2024 16:35:02 -0700 Subject: [PATCH 015/405] Update 4 of #331. Now it is only the documentation that is updated. --- doc/sphinx/src/models.rst | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 6458601b78..60492810cf 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -759,7 +759,7 @@ Mie-Gruneisen linear :math:`U_s`-:math:`u_p` EOS One of the most commonly-used EOS is the linear :math:`U_s`-:math:`u_p` version of the Mie-Gruneisen EOS. This EOS uses the Hugoniot as the reference curve and is extensively used in shock physics. -This version implements the exact thermodynamic tempearture on the Hugoniot and also adds an entropy. +This version implements the exact thermodynamic temperature on the Hugoniot and also adds an entropy. The pressure follows the traditional Mie-Gruneisen form, @@ -848,7 +848,11 @@ it is .. math:: T_H(\rho) = T_0 e^{\Gammma(\rho_0) \eta} + \frac{e^{\Gammma(\rho_0) \eta}}{2 C_V \rho_0} - \int_0^\eta e^{-\gamma(\rho_0) z} z^2 \frac{d}{dz} \left( \frac{P_H}{z}\right) dz \\ + \int_0^\eta e^{-\gamma(\rho_0) z} z^2 \frac{d}{dz} \left( \frac{P_H}{z}\right) dz + + +testing +.. math:: = T_0 e^{\Gammma(\rho_0) \eta} + \frac{C_s^2}}{2 C_V s^2} \left[\frac{- s \eta}{(1 - s \eta)^2} + \left( \frac{\Gamma(\rho_0)}{s} - 3 \right) \left( e^{\Gammma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right) From f91c0f568b12a5b908e89fef06b8c3c798ab43c8 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Wed, 3 Jan 2024 17:03:11 -0700 Subject: [PATCH 016/405] Update 5 for #331. --- doc/sphinx/src/models.rst | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 60492810cf..d8e852af48 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -754,10 +754,10 @@ is :math:`S_0`, and ``expconsts`` is a pointer to the constant array of length :math:`d_2` to :math:`d_{40}`. Expansion coefficients not used should be set to 0.0. -Mie-Gruneisen linear :math:`U_s`-:math:`u_p` EOS -```````````````````````````````````````````````` +Mie-Gruneisen linear :math:`U_s`- :math:`u_p` EOS +````````````````````````````````````````````````` -One of the most commonly-used EOS is the linear :math:`U_s`-:math:`u_p` version of the Mie-Gruneisen EOS. This EOS +One of the most commonly-used EOS is the linear :math:`U_s`- :math:`u_p` version of the Mie-Gruneisen EOS. This EOS uses the Hugoniot as the reference curve and is extensively used in shock physics. This version implements the exact thermodynamic temperature on the Hugoniot and also adds an entropy. @@ -817,7 +817,7 @@ shock jump conditions, \rho_0 U_s = \rho (U_s - u_p) \\ P_H = \rho_0 U_s u_p , -assuming a linear :math:`U_s`-:math:`u_p` relation, +assuming a linear :math:`U_s`- :math:`u_p` relation, .. math:: @@ -825,7 +825,7 @@ assuming a linear :math:`U_s`-:math:`u_p` relation, Here :math:`U_s` is the shock velocity and :math:`u_p` is the particle velocity. As is pointed out in the description of the Gruneisen EOS, -for many materials, the :math:`U_s`-:math:`u_p` relationship is roughly linear +for many materials, the :math:`U_s`- :math:`u_p` relationship is roughly linear so only this :math:`s` parameter is needed. The units for :math:`C_s` is velocity while :math:`s` is unitless. Note that the parameter :math:`s` is related to the fundamental derivative of shock physics as shown by `Wills `_. @@ -836,26 +836,32 @@ Solving the jump equations above gives that the reference pressure along the Hug P_H(\rho) = C_s^2 \rho_0 \frac{\eta}{\left(1 - s \eta \right)^2} . +Note the singularity at :math:`s \eta = 1` which limits this model's validity to compressions +:math:`\eta << 1/s`. If your problem can be expected to have compressions of this order, you should use the PowerMG +EOS that is explicitely constructed for large compressions. +The assumption of linear :math:`U_s`- :math:`u_p` relation is simply not valid at large compressions. + The energy along the Hugoniot is given by .. math:: E_H(\rho) = \frac{P_H \eta }{2 \rho_0} + E_0 . -The temperature on the Hugoniot is hard to derive but with the help of Mathematica -it is +The temperature on the Hugoniot is hard to derive explicitely but with the help of Mathematica +we can solve .. math:: - T_H(\rho) = T_0 e^{\Gammma(\rho_0) \eta} + \frac{e^{\Gammma(\rho_0) \eta}}{2 C_V \rho_0} + T_H(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{e^{\Gamma(\rho_0) \eta}}{2 C_V \rho_0} \int_0^\eta e^{-\gamma(\rho_0) z} z^2 \frac{d}{dz} \left( \frac{P_H}{z}\right) dz -testing +into the explicit formula + .. math:: - = T_0 e^{\Gammma(\rho_0) \eta} + \frac{C_s^2}}{2 C_V s^2} + T_H(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}}{2 C_V s^2} \left[\frac{- s \eta}{(1 - s \eta)^2} + \left( \frac{\Gamma(\rho_0)}{s} - 3 \right) - \left( e^{\Gammma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right) + \left( e^{\Gamma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right) + e^{-\frac{\Gamma(\rho_0)}{s} (1-s \eta)} \left( Ei(\frac{\Gamma(\rho_0)}{s}(1-s \eta))-Ei(\frac{\Gamma(\rho_0)}{s}) \right) \left((\frac{\Gamma(\rho_0)}{s})^2 - 4 \frac{\Gamma(\rho_0)}{s} + 2 \right) \right] From 61ad5db65fa66ad1c7c29a754b81c320eda29092 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Wed, 3 Jan 2024 17:20:33 -0700 Subject: [PATCH 017/405] Update 6 for #331. --- doc/sphinx/src/models.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index d8e852af48..87a6bdcddd 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -853,13 +853,13 @@ we can solve .. math:: T_H(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{e^{\Gamma(\rho_0) \eta}}{2 C_V \rho_0} - \int_0^\eta e^{-\gamma(\rho_0) z} z^2 \frac{d}{dz} \left( \frac{P_H}{z}\right) dz + \int_0^\eta e^{-\Gamma(\rho_0) z} z^2 \frac{d}{dz} \left( \frac{P_H}{z}\right) dz into the explicit formula .. math:: - T_H(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}}{2 C_V s^2} + T_H(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}{2 C_V s^2} \left[\frac{- s \eta}{(1 - s \eta)^2} + \left( \frac{\Gamma(\rho_0)}{s} - 3 \right) \left( e^{\Gamma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right) + e^{-\frac{\Gamma(\rho_0)}{s} (1-s \eta)} From a35744ef32c86f4b978434b9e0085e3d9f17c96e Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Wed, 3 Jan 2024 17:43:53 -0700 Subject: [PATCH 018/405] Update 7 of #331. --- doc/sphinx/src/models.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 87a6bdcddd..97b5ca8b61 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -815,7 +815,7 @@ shock jump conditions, .. math:: \rho_0 U_s = \rho (U_s - u_p) \\ - P_H = \rho_0 U_s u_p , + P_H = \rho_0 U_s u_p \ \ \ \ \ \ , assuming a linear :math:`U_s`- :math:`u_p` relation, @@ -861,7 +861,7 @@ into the explicit formula .. math:: T_H(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}{2 C_V s^2} \left[\frac{- s \eta}{(1 - s \eta)^2} + \left( \frac{\Gamma(\rho_0)}{s} - 3 \right) - \left( e^{\Gamma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right) + \left( e^{\Gamma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right) \\ + e^{-\frac{\Gamma(\rho_0)}{s} (1-s \eta)} \left( Ei(\frac{\Gamma(\rho_0)}{s}(1-s \eta))-Ei(\frac{\Gamma(\rho_0)}{s}) \right) \left((\frac{\Gamma(\rho_0)}{s})^2 - 4 \frac{\Gamma(\rho_0)}{s} + 2 \right) \right] From b54f84092e3035676ffdef8092023c31a4409a3f Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Wed, 3 Jan 2024 18:36:29 -0700 Subject: [PATCH 019/405] Update 8 for #331. --- doc/sphinx/src/models.rst | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 97b5ca8b61..94755ed1ce 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -815,7 +815,7 @@ shock jump conditions, .. math:: \rho_0 U_s = \rho (U_s - u_p) \\ - P_H = \rho_0 U_s u_p \ \ \ \ \ \ , + P_H = \rho_0 U_s u_p \ , \ \ \ \ assuming a linear :math:`U_s`- :math:`u_p` relation, @@ -861,11 +861,18 @@ into the explicit formula .. math:: T_H(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}{2 C_V s^2} \left[\frac{- s \eta}{(1 - s \eta)^2} + \left( \frac{\Gamma(\rho_0)}{s} - 3 \right) - \left( e^{\Gamma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right) \\ - + e^{-\frac{\Gamma(\rho_0)}{s} (1-s \eta)} + \left( e^{\Gamma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right) \\ + \ \ \ \ \ \ \ + e^{-\frac{\Gamma(\rho_0)}{s} (1-s \eta)} \left( Ei(\frac{\Gamma(\rho_0)}{s}(1-s \eta))-Ei(\frac{\Gamma(\rho_0)}{s}) \right) \left((\frac{\Gamma(\rho_0)}{s})^2 - 4 \frac{\Gamma(\rho_0)}{s} + 2 \right) \right] +where :math:`Ei` is the exponential integral function. We replace the :math:`Ei` difference with a sum with cutoff +giving an error less than machine precision. For :math:`s \eta` close to :math:`0`, there are +severe cancellations in this formula and we use the expansion + +.. math:: + {T_H}_{exp}(\rho) = \ . + The constructor for the ``MGUsup`` EOS has the signature From b9aef7c8865ba7a0916ecb31f73a5d441d8cc5d9 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Wed, 3 Jan 2024 19:06:35 -0700 Subject: [PATCH 020/405] Update 9 for #331. --- doc/sphinx/src/models.rst | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 94755ed1ce..defb8948b8 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -793,7 +793,7 @@ thermodynamic derivatives we can show that \Gamma \rho = \frac{\alpha B_T}{C_V} , -and we arrive at +and we arrive at .. math:: @@ -814,8 +814,9 @@ The pressure, energy, and temperature, on the Hugoniot are derived from the shock jump conditions, .. math:: - \rho_0 U_s = \rho (U_s - u_p) \\ - P_H = \rho_0 U_s u_p \ , \ \ \ \ + + \rho_0 U_s &= \rho (U_s - u_p) \\ + P_H &= \rho_0 U_s u_p \ , assuming a linear :math:`U_s`- :math:`u_p` relation, @@ -859,10 +860,11 @@ we can solve into the explicit formula .. math:: - T_H(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}{2 C_V s^2} + + T_H(\rho) &= T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}{2 C_V s^2} \left[\frac{- s \eta}{(1 - s \eta)^2} + \left( \frac{\Gamma(\rho_0)}{s} - 3 \right) - \left( e^{\Gamma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right) \\ - \ \ \ \ \ \ \ + e^{-\frac{\Gamma(\rho_0)}{s} (1-s \eta)} + \left( e^{\Gamma(\rho_0) \eta} - \frac{1}{(1-s \eta)}\right)\right. \\ + & \ \left. + e^{-\frac{\Gamma(\rho_0)}{s} (1-s \eta)} \left( Ei(\frac{\Gamma(\rho_0)}{s}(1-s \eta))-Ei(\frac{\Gamma(\rho_0)}{s}) \right) \left((\frac{\Gamma(\rho_0)}{s})^2 - 4 \frac{\Gamma(\rho_0)}{s} + 2 \right) \right] @@ -871,7 +873,9 @@ giving an error less than machine precision. For :math:`s \eta` close to :math:` severe cancellations in this formula and we use the expansion .. math:: - {T_H}_{exp}(\rho) = \ . + + {T_H}_{exp}(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}{2 C_V s^2} + \left[ -2 \ln ( 1- s \eta) + \frac{s \eta}{(1 - s \eta)^2} ( 3 s \eta - 2) \ . The constructor for the ``MGUsup`` EOS has the signature From 8ea971abde63482b71e9038b69091e21eb7de8bd Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Thu, 4 Jan 2024 10:34:40 -0700 Subject: [PATCH 021/405] Update 10 for #331. --- doc/sphinx/src/models.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index defb8948b8..8b40073e2d 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -874,8 +874,8 @@ severe cancellations in this formula and we use the expansion .. math:: - {T_H}_{exp}(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}{2 C_V s^2} - \left[ -2 \ln ( 1- s \eta) + \frac{s \eta}{(1 - s \eta)^2} ( 3 s \eta - 2) \ . + {T_H}_{exp}(\rho) = T_0 e^{\Gamma(\rho_0) \eta} + \frac{C_s^2}{2 C_V s^2} + \left[ -2 \ln ( 1- s \eta) + \frac{s \eta}{(1 - s \eta)^2} ( 3 s \eta - 2) \right] \ . The constructor for the ``MGUsup`` EOS has the signature From efa4ad29d28587014870625cd55d0aab20643726 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Thu, 4 Jan 2024 11:00:35 -0700 Subject: [PATCH 022/405] Update 11 for #331. Hopefully the last. --- doc/sphinx/src/models.rst | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 8b40073e2d..e6157ddb7b 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -870,7 +870,7 @@ into the explicit formula where :math:`Ei` is the exponential integral function. We replace the :math:`Ei` difference with a sum with cutoff giving an error less than machine precision. For :math:`s \eta` close to :math:`0`, there are -severe cancellations in this formula and we use the expansion +severe cancellations in this formula and we use the expansion .. math:: @@ -878,6 +878,10 @@ severe cancellations in this formula and we use the expansion \left[ -2 \ln ( 1- s \eta) + \frac{s \eta}{(1 - s \eta)^2} ( 3 s \eta - 2) \right] \ . +The first omitted term in the expansion inside the square brackets is :math:`\Gamma(\rho_0) \eta^4 / 6`. This exapnsion is +in fact even better than the common approximation of replacing the full temperature on the Hugoniot with the temperature on the +isentrope, that is, the first term :math:`T_0 e^{\Gamma(\rho_0) \eta}`. + The constructor for the ``MGUsup`` EOS has the signature .. code-block:: cpp From 6fb8d9ea85c9b18f5660aef5c15f514f6ceb00f7 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Thu, 4 Jan 2024 11:17:05 -0700 Subject: [PATCH 023/405] Update 12 for #331. This is the last. --- doc/sphinx/src/models.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index e6157ddb7b..e32f2a5e47 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -878,7 +878,7 @@ severe cancellations in this formula and we use the expansion \left[ -2 \ln ( 1- s \eta) + \frac{s \eta}{(1 - s \eta)^2} ( 3 s \eta - 2) \right] \ . -The first omitted term in the expansion inside the square brackets is :math:`\Gamma(\rho_0) \eta^4 / 6`. This exapnsion is +The first omitted term in the expansion inside the square brackets is :math:`\Gamma(\rho_0) \eta^4 / 6`. This expansion is in fact even better than the common approximation of replacing the full temperature on the Hugoniot with the temperature on the isentrope, that is, the first term :math:`T_0 e^{\Gamma(\rho_0) \eta}`. From cab0ae5d317b52f65064d29761a55221a6f5ebf7 Mon Sep 17 00:00:00 2001 From: Gopinath Subramanian Date: Tue, 9 Jan 2024 11:04:14 -0700 Subject: [PATCH 024/405] Fixed typo in documentation --- doc/sphinx/src/models.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index e32f2a5e47..7a28eabc00 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1298,7 +1298,7 @@ SAP Polynomial EOS `````````````````` This model is specific to the Safety Applications Project (SAP). It is -an imcomplete EOS, and is a simple analytical form used to fit +an incomplete EOS, and is a simple analytical form used to fit experimental data: .. math:: From b50858b278904a66b449c11724f7931042871df7 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 10 Jan 2024 17:15:56 -0700 Subject: [PATCH 025/405] first pass at plugins --- CMakeLists.txt | 21 +- cmake/plugins.cmake | 49 +++++ .../CMakeDirectoryInformation.cmake | 16 ++ example/plugin/CMakeFiles/progress.marks | 1 + example/plugin/CMakeLists.txt | 22 ++ example/plugin/Makefile | 202 ++++++++++++++++++ example/plugin/cmake_install.cmake | 44 ++++ example/plugin/dust.hpp | 144 +++++++++++++ example/plugin/test_dust.cpp | 48 +++++ singularity-eos/CMakeLists.txt | 23 +- test/CMakeLists.txt | 16 ++ 11 files changed, 564 insertions(+), 22 deletions(-) create mode 100644 cmake/plugins.cmake create mode 100644 example/plugin/CMakeFiles/CMakeDirectoryInformation.cmake create mode 100644 example/plugin/CMakeFiles/progress.marks create mode 100644 example/plugin/CMakeLists.txt create mode 100644 example/plugin/Makefile create mode 100644 example/plugin/cmake_install.cmake create mode 100644 example/plugin/dust.hpp create mode 100644 example/plugin/test_dust.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 7353b4aea1..7452b4d4ff 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -118,6 +118,9 @@ option(SINGULARITY_FORCE_SUBMODULE_MODE "Submodule mode" OFF) option(SINGULARITY_PATCH_MPARK_VARIANT "Apply GPU patch to mpark-variant submodule" ON) +# Plugins +set(SINGULARITY_PLUGINS "" CACHE STRING "List of paths to plugin directories") + # ------------------------------------------------------------------------------# # singularity-eos Library # ------------------------------------------------------------------------------# @@ -359,27 +362,27 @@ if(SINGULARITY_BUILD_TESTS) endif() endif() -# ##########################OLD # ------------------------------------------------------------------------------# # singularity-eos library # ------------------------------------------------------------------------------# -# this subdirectory populates `EOS_HEADERS/EOS_SRCS` NOTE: these include path +# Here we populate `EOS_HEADERS/EOS_SRCS` NOTE: these include path # prefixes of subdirectories on files (e.g. eos/eos.hpp) see # singularity-eos/CMakeLists.txt +include(cmake/plugins.cmake) add_subdirectory(singularity-eos) -foreach(_header ${EOS_HEADERS}) - list(APPEND _install_headers ${_header}) - list(APPEND _headers singularity-eos/${_header}) +foreach(_plugin ${SINGULARITY_PLUGINS}) + message(STATUS "Adding plugin ${_plugin}...") + add_subdirectory(${_plugin} ${_plugin}) endforeach() -foreach(_src ${EOS_SRCS}) - list(APPEND _srcs singularity-eos/${_src}) -endforeach() +# TODO(JMM): Kind of nice to have? +message(STATUS "EOS Headers:\n\t${EOS_HEADERS}") +message(STATUS "EOS Sources:\n\t${EOS_SRCS}") -target_sources(singularity-eos PRIVATE ${_srcs} ${_headers}) +target_sources(singularity-eos PRIVATE ${EOS_SRCS} ${EOS_HEADERS}) if(SINGULARITY_USE_FORTRAN) # Turn on preprocessor for fortran files diff --git a/cmake/plugins.cmake b/cmake/plugins.cmake new file mode 100644 index 0000000000..11ad8a6ce9 --- /dev/null +++ b/cmake/plugins.cmake @@ -0,0 +1,49 @@ +#------------------------------------------------------------------------------# +# © 2024. 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. +#------------------------------------------------------------------------------# + +set(EOS_HEADERS "") +set(_install_headers "") +function(register_headers) + foreach(arg IN LISTS ARGN) + file(RELATIVE_PATH relative_path ${CMAKE_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}) + list(APPEND EOS_HEADERS ${relative_path}/${arg}) + list(APPEND _install_headers ${arg}) + endforeach() + set(EOS_HEADERS ${EOS_HEADERS} PARENT_SCOPE) + set(_install_headers ${_install_headers} PARENT_SCOPE) +endfunction() + +set(EOS_SRCS "") +function(register_srcs) + foreach(arg IN LISTS ARGN) + file(RELATIVE_PATH relative_path ${CMAKE_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}) + list(APPEND EOS_SRCS ${relative_path}/${arg}) + endforeach() + set(EOS_SRCS ${EOS_SRCS} PARENT_SCOPE) +endfunction() + +set(PLUGIN_TESTS "") +function(register_tests) + foreach(arg IN LISTS ARGN) + list(APPEND PLUGIN_TESTS ${CMAKE_CURRENT_SOURCE_DIR}/${arg}) + endforeach() + set(PLUGIN_TESTS ${PLUGIN_TESTS} PARENT_SCOPE) +endfunction() + +macro(export_plugin) + set(EOS_HEADERS ${EOS_HEADERS} PARENT_SCOPE) + set(_install_headers ${_install_headers} PARENT_SCOPE) + set(EOS_SRCS ${EOS_SRCS} PARENT_SCOPE) + set(PLUGIN_TESTS ${PLUGIN_TESTS} PARENT_SCOPE) +endmacro() diff --git a/example/plugin/CMakeFiles/CMakeDirectoryInformation.cmake b/example/plugin/CMakeFiles/CMakeDirectoryInformation.cmake new file mode 100644 index 0000000000..d5f4c60f97 --- /dev/null +++ b/example/plugin/CMakeFiles/CMakeDirectoryInformation.cmake @@ -0,0 +1,16 @@ +# CMAKE generated file: DO NOT EDIT! +# Generated by "Unix Makefiles" Generator, CMake Version 3.25 + +# Relative path conversion top directories. +set(CMAKE_RELATIVE_PATH_TOP_SOURCE "/home/jonahm/programming/singularity-eos") +set(CMAKE_RELATIVE_PATH_TOP_BINARY "/home/jonahm/programming/singularity-eos/example/plugin") + +# Force unix paths in dependencies. +set(CMAKE_FORCE_UNIX_PATHS 1) + + +# The C and CXX include file regular expressions for this directory. +set(CMAKE_C_INCLUDE_REGEX_SCAN "^.*$") +set(CMAKE_C_INCLUDE_REGEX_COMPLAIN "^$") +set(CMAKE_CXX_INCLUDE_REGEX_SCAN ${CMAKE_C_INCLUDE_REGEX_SCAN}) +set(CMAKE_CXX_INCLUDE_REGEX_COMPLAIN ${CMAKE_C_INCLUDE_REGEX_COMPLAIN}) diff --git a/example/plugin/CMakeFiles/progress.marks b/example/plugin/CMakeFiles/progress.marks new file mode 100644 index 0000000000..573541ac97 --- /dev/null +++ b/example/plugin/CMakeFiles/progress.marks @@ -0,0 +1 @@ +0 diff --git a/example/plugin/CMakeLists.txt b/example/plugin/CMakeLists.txt new file mode 100644 index 0000000000..16f4616c6b --- /dev/null +++ b/example/plugin/CMakeLists.txt @@ -0,0 +1,22 @@ +#------------------------------------------------------------------------------# +# © 2024. 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. +#------------------------------------------------------------------------------# + +# Add the header file to the list of headers +register_headers(dust.hpp) + +# Register this test to be built in the test suite +register_tests(test_dust.cpp) + +# Export lists to parent scope +export_plugin() diff --git a/example/plugin/Makefile b/example/plugin/Makefile new file mode 100644 index 0000000000..b4c8c93e1e --- /dev/null +++ b/example/plugin/Makefile @@ -0,0 +1,202 @@ +# CMAKE generated file: DO NOT EDIT! +# Generated by "Unix Makefiles" Generator, CMake Version 3.25 + +# Default target executed when no arguments are given to make. +default_target: all +.PHONY : default_target + +# Allow only one "make -f Makefile2" at a time, but pass parallelism. +.NOTPARALLEL: + +#============================================================================= +# Special targets provided by cmake. + +# Disable implicit rules so canonical targets will work. +.SUFFIXES: + +# Disable VCS-based implicit rules. +% : %,v + +# Disable VCS-based implicit rules. +% : RCS/% + +# Disable VCS-based implicit rules. +% : RCS/%,v + +# Disable VCS-based implicit rules. +% : SCCS/s.% + +# Disable VCS-based implicit rules. +% : s.% + +.SUFFIXES: .hpux_make_needs_suffix_list + +# Command-line flag to silence nested $(MAKE). +$(VERBOSE)MAKESILENT = -s + +#Suppress display of executed commands. +$(VERBOSE).SILENT: + +# A target that is always out of date. +cmake_force: +.PHONY : cmake_force + +#============================================================================= +# Set environment variables for the build. + +# The shell in which to execute make rules. +SHELL = /bin/sh + +# The CMake executable. +CMAKE_COMMAND = /usr/bin/cmake + +# The command to remove a file. +RM = /usr/bin/cmake -E rm -f + +# Escaping for special characters. +EQUALS = = + +# The top-level source directory on which CMake was run. +CMAKE_SOURCE_DIR = /home/jonahm/programming/singularity-eos + +# The top-level build directory on which CMake was run. +CMAKE_BINARY_DIR = /home/jonahm/programming/singularity-eos/build + +#============================================================================= +# Targets provided globally by CMake. + +# Special rule for the target test +test: + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running tests..." + /usr/bin/ctest --force-new-ctest-process $(ARGS) +.PHONY : test + +# Special rule for the target test +test/fast: test +.PHONY : test/fast + +# Special rule for the target edit_cache +edit_cache: + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake cache editor..." + /usr/bin/ccmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) +.PHONY : edit_cache + +# Special rule for the target edit_cache +edit_cache/fast: edit_cache +.PHONY : edit_cache/fast + +# Special rule for the target rebuild_cache +rebuild_cache: + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake to regenerate build system..." + /usr/bin/cmake --regenerate-during-build -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) +.PHONY : rebuild_cache + +# Special rule for the target rebuild_cache +rebuild_cache/fast: rebuild_cache +.PHONY : rebuild_cache/fast + +# Special rule for the target list_install_components +list_install_components: + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Available install components are: \"Unspecified\"" +.PHONY : list_install_components + +# Special rule for the target list_install_components +list_install_components/fast: list_install_components +.PHONY : list_install_components/fast + +# Special rule for the target install +install: preinstall + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Install the project..." + /usr/bin/cmake -P cmake_install.cmake +.PHONY : install + +# Special rule for the target install +install/fast: preinstall/fast + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Install the project..." + /usr/bin/cmake -P cmake_install.cmake +.PHONY : install/fast + +# Special rule for the target install/local +install/local: preinstall + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Installing only the local directory..." + /usr/bin/cmake -DCMAKE_INSTALL_LOCAL_ONLY=1 -P cmake_install.cmake +.PHONY : install/local + +# Special rule for the target install/local +install/local/fast: preinstall/fast + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Installing only the local directory..." + /usr/bin/cmake -DCMAKE_INSTALL_LOCAL_ONLY=1 -P cmake_install.cmake +.PHONY : install/local/fast + +# Special rule for the target install/strip +install/strip: preinstall + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Installing the project stripped..." + /usr/bin/cmake -DCMAKE_INSTALL_DO_STRIP=1 -P cmake_install.cmake +.PHONY : install/strip + +# Special rule for the target install/strip +install/strip/fast: preinstall/fast + @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Installing the project stripped..." + /usr/bin/cmake -DCMAKE_INSTALL_DO_STRIP=1 -P cmake_install.cmake +.PHONY : install/strip/fast + +# The main all target +all: cmake_check_build_system + cd /home/jonahm/programming/singularity-eos/build && $(CMAKE_COMMAND) -E cmake_progress_start /home/jonahm/programming/singularity-eos/build/CMakeFiles /home/jonahm/programming/singularity-eos/example/plugin//CMakeFiles/progress.marks + cd /home/jonahm/programming/singularity-eos/build && $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 /home/jonahm/programming/singularity-eos/example/plugin/all + $(CMAKE_COMMAND) -E cmake_progress_start /home/jonahm/programming/singularity-eos/build/CMakeFiles 0 +.PHONY : all + +# The main clean target +clean: + cd /home/jonahm/programming/singularity-eos/build && $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 /home/jonahm/programming/singularity-eos/example/plugin/clean +.PHONY : clean + +# The main clean target +clean/fast: clean +.PHONY : clean/fast + +# Prepare targets for installation. +preinstall: all + cd /home/jonahm/programming/singularity-eos/build && $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 /home/jonahm/programming/singularity-eos/example/plugin/preinstall +.PHONY : preinstall + +# Prepare targets for installation. +preinstall/fast: + cd /home/jonahm/programming/singularity-eos/build && $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 /home/jonahm/programming/singularity-eos/example/plugin/preinstall +.PHONY : preinstall/fast + +# clear depends +depend: + cd /home/jonahm/programming/singularity-eos/build && $(CMAKE_COMMAND) -P /home/jonahm/programming/singularity-eos/build/CMakeFiles/VerifyGlobs.cmake + cd /home/jonahm/programming/singularity-eos/build && $(CMAKE_COMMAND) -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) --check-build-system CMakeFiles/Makefile.cmake 1 +.PHONY : depend + +# Help Target +help: + @echo "The following are some of the valid targets for this Makefile:" + @echo "... all (the default if no target is provided)" + @echo "... clean" + @echo "... depend" + @echo "... edit_cache" + @echo "... install" + @echo "... install/local" + @echo "... install/strip" + @echo "... list_install_components" + @echo "... rebuild_cache" + @echo "... test" +.PHONY : help + + + +#============================================================================= +# Special targets to cleanup operation of make. + +# Special rule to run CMake to check the build system integrity. +# No rule that depends on this can have commands that come from listfiles +# because they might be regenerated. +cmake_check_build_system: + cd /home/jonahm/programming/singularity-eos/build && $(CMAKE_COMMAND) -P /home/jonahm/programming/singularity-eos/build/CMakeFiles/VerifyGlobs.cmake + cd /home/jonahm/programming/singularity-eos/build && $(CMAKE_COMMAND) -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) --check-build-system CMakeFiles/Makefile.cmake 0 +.PHONY : cmake_check_build_system + diff --git a/example/plugin/cmake_install.cmake b/example/plugin/cmake_install.cmake new file mode 100644 index 0000000000..88c1fd1691 --- /dev/null +++ b/example/plugin/cmake_install.cmake @@ -0,0 +1,44 @@ +# Install script for directory: /home/jonahm/programming/singularity-eos/example/plugin + +# Set the install prefix +if(NOT DEFINED CMAKE_INSTALL_PREFIX) + set(CMAKE_INSTALL_PREFIX "/home/jonahm/test") +endif() +string(REGEX REPLACE "/$" "" CMAKE_INSTALL_PREFIX "${CMAKE_INSTALL_PREFIX}") + +# Set the install configuration name. +if(NOT DEFINED CMAKE_INSTALL_CONFIG_NAME) + if(BUILD_TYPE) + string(REGEX REPLACE "^[^A-Za-z0-9_]+" "" + CMAKE_INSTALL_CONFIG_NAME "${BUILD_TYPE}") + else() + set(CMAKE_INSTALL_CONFIG_NAME "") + endif() + message(STATUS "Install configuration: \"${CMAKE_INSTALL_CONFIG_NAME}\"") +endif() + +# Set the component getting installed. +if(NOT CMAKE_INSTALL_COMPONENT) + if(COMPONENT) + message(STATUS "Install component: \"${COMPONENT}\"") + set(CMAKE_INSTALL_COMPONENT "${COMPONENT}") + else() + set(CMAKE_INSTALL_COMPONENT) + endif() +endif() + +# Install shared libraries without execute permission? +if(NOT DEFINED CMAKE_INSTALL_SO_NO_EXE) + set(CMAKE_INSTALL_SO_NO_EXE "1") +endif() + +# Is this installation the result of a crosscompile? +if(NOT DEFINED CMAKE_CROSSCOMPILING) + set(CMAKE_CROSSCOMPILING "FALSE") +endif() + +# Set default install directory permissions. +if(NOT DEFINED CMAKE_OBJDUMP) + set(CMAKE_OBJDUMP "/usr/bin/objdump") +endif() + diff --git a/example/plugin/dust.hpp b/example/plugin/dust.hpp new file mode 100644 index 0000000000..e99828661f --- /dev/null +++ b/example/plugin/dust.hpp @@ -0,0 +1,144 @@ +//------------------------------------------------------------------------------ +// © 2024. 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 _EXAMPLE_PLUGIN_DUST_ +#define _EXAMPLE_PLUGIN_DUST_ + +// Ports-of-call +#include + +// Base stuff +#include +#include +#include +#include + +namespace singularity { + +using namespace eos_base; + +class Dust : public EosBase { + public: + PORTABLE_INLINE_FUNCTION Dust() : _Cv(1) {} + PORTABLE_INLINE_FUNCTION Dust(Real Cv) : _Cv(Cv) {} + + Dust GetOnDevice() { return *this; } + inline void Finalize() {} + + PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return sie / _Cv; + } + PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) { + return _Cv * temperature; + } + PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + return 0; + } + PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return 0; + } + PORTABLE_INLINE_FUNCTION Real + MinInternalEnergyFromDensity(const Real rho, Real *lambda = nullptr) const { + return 0.0; + }; + + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + return 0; + } + PORTABLE_INLINE_FUNCTION Real EntropyFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + const Real temp = TemperatureFromDensityInternalEnergy(rho, sie, lambda); + return 0; + } + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + return _Cv; + } + PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return _Cv; + } + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + return 0; + } + PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return 0; + } + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityTemperature( + const Real rho, const Real temperature, Real *lambda = nullptr) const { + return 0; + } + PORTABLE_INLINE_FUNCTION Real GruneisenParamFromDensityInternalEnergy( + const Real rho, const Real sie, Real *lambda = nullptr) const { + return 0; + } + PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, + Real &cv, Real &bmod, const unsigned long output, + Real *lambda = nullptr) const { + press = 0; + bmod = 0; + cv = _Cv; + if (output & thermalqs::density) { + UNDEFINED_ERROR; + } + if (output & thermalqs::specific_internal_energy) { + if (!(output & thermalqs::temperature)) { + UNDEFINED_ERROR; + } + energy = _Cv*temp; + } + if (output & thermalqs::temperature) { + if (!(output & thermalqs::specific_internal_energy)) { + temp = energy / _Cv; + } + } + }; + PORTABLE_INLINE_FUNCTION + void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, + Real &bmod, Real &dpde, Real &dvdt, + Real *lambda = nullptr) const { + rho = temp = sie = 1; + press = bmod = dpde = 0; + dvdt = 0; + } + SG_ADD_BASE_CLASS_USINGS(Dust) + + PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; } + static constexpr unsigned long PreferredInput() { return _preferred_input; } + static inline unsigned long scratch_size(std::string method, unsigned int nelements) { + return 0; + } + static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } + PORTABLE_INLINE_FUNCTION void PrintParams() const { + printf("Dust Parameters:\nCv = %g\n", _Cv); + } + static std::string EosType() { return std::string("Dust"); } + static std::string EosPyType() { return EosType(); } + +private: + Real _Cv = 1; + static constexpr const unsigned long _preferred_input = + thermalqs::density | thermalqs::specific_internal_energy | thermalqs::temperature; +}; + +} // namespace singularity + +#endif // _EXAMPLE_PLUGIN_DUST_ diff --git a/example/plugin/test_dust.cpp b/example/plugin/test_dust.cpp new file mode 100644 index 0000000000..bef73438a6 --- /dev/null +++ b/example/plugin/test_dust.cpp @@ -0,0 +1,48 @@ +//------------------------------------------------------------------------------ +// © 2024. 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. +//------------------------------------------------------------------------------ + +#include +#include +#include +#include + +#include +#include +#include +#include + +#ifdef SINGULARITY_BUILD_CLOSURE +#include +#endif + +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +#include +#include "dust.hpp" + +using singularity::Dust; +using EOS = singularity::Variant; + +SCENARIO("Dust EOS", "[Dust]") { + GIVEN("A dust EOS") { + constexpr Real Cv = 1; + EOS dust = Dust(Cv); + THEN("The pressure is zero") { + REQUIRE(dust.PressureFromDensityInternalEnergy(1, 1) == 0); + } + } +} diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 46c86316a7..82deec6e2c 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -12,7 +12,7 @@ # publicly and display publicly, and to permit others to do so. #------------------------------------------------------------------------------# -set(EOS_HEADERS +register_headers( base/fast-math/logs.hpp base/robust_utils.hpp base/root-finding-1d/root_finding.hpp @@ -45,33 +45,30 @@ set(EOS_HEADERS eos/eos_stiff.hpp ) -set(EOS_SRCS eos/eos.cpp) +register_srcs(eos/eos.cpp) if (SINGULARITY_BUILD_CLOSURE) - list(APPEND EOS_HEADERS - closure/mixed_cell_models.hpp) + register_headers(closure/mixed_cell_models.hpp) if (SINGULARITY_USE_FORTRAN) - list(APPEND EOS_SRCS - eos/get_sg_eos.cpp) + register_srcs(eos/get_sg_eos.cpp) if (SINGULARITY_USE_KOKKOS) - list(APPEND EOS_SRCS + register_srcs( eos/get_sg_eos_p_t.cpp eos/get_sg_eos_rho_t.cpp eos/get_sg_eos_rho_p.cpp eos/get_sg_eos_rho_e.cpp) endif() - list(APPEND EOS_HEADERS + register_headers( eos/get_sg_eos.hpp eos/get_sg_eos_functors.hpp) endif() endif() if (SINGULARITY_USE_FORTRAN) - list(APPEND EOS_SRCS eos/singularity_eos.f90) - list(APPEND EOS_SRCS eos/singularity_eos.cpp) - list(APPEND EOS_HEADERS eos/singularity_eos.hpp) + register_srcs(eos/singularity_eos.f90) + register_srcs(eos/singularity_eos.cpp) + register_headers(eos/singularity_eos.hpp) endif() # export to parent scope -set(EOS_HEADERS ${EOS_HEADERS} PARENT_SCOPE) -set(EOS_SRCS ${EOS_SRCS} PARENT_SCOPE) +export_plugin() diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6add5dbd39..5830f0b0cd 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -43,6 +43,15 @@ add_executable( test_eos_stellar_collapse.cpp ) +if (PLUGIN_TESTS) + add_executable( + eos_plugin_tests + catch2_define.cpp + eos_unit_test_helpers.hpp + ${PLUGIN_TESTS} + ) +endif() + if(SINGULARITY_TEST_HELMHOLTZ) configure_file(${PROJECT_SOURCE_DIR}/data/helmholtz/helm_table.dat ${CMAKE_BINARY_DIR}/data/helmholtz/helm_table.dat COPYONLY) target_compile_definitions(eos_tabulated_unit_tests PRIVATE SINGULARITY_TEST_HELMHOLTZ SINGULARITY_USE_HELMHOLTZ) @@ -62,10 +71,17 @@ target_link_libraries(eos_infrastructure_tests PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) target_link_libraries(eos_tabulated_unit_tests PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) +if (PLUGIN_TESTS) + target_link_libraries(eos_plugin_tests PRIVATE Catch2::Catch2 + singularity-eos::singularity-eos) +endif() include(Catch) catch_discover_tests(eos_analytic_unit_tests PROPERTIES TIMEOUT 60) catch_discover_tests(eos_infrastructure_tests PROPERTIES TIMEOUT 60) catch_discover_tests(eos_tabulated_unit_tests PROPERTIES TIMEOUT 60) +if (PLUGIN_TESTS) + catch_discover_tests(eos_plugin_tests PROPERTIES TIMEOUT 60) +endif() if(SINGULARITY_USE_EOSPAC AND SINGULARITY_TEST_SESAME) add_executable(compare_to_eospac compare_to_eospac.cpp) From 900c4ffc7252a5ed1977a26e57da91051a8a2b8d Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 10 Jan 2024 17:20:14 -0700 Subject: [PATCH 026/405] remove extraneous cmake files that got added --- .../CMakeDirectoryInformation.cmake | 16 -- example/plugin/CMakeFiles/progress.marks | 1 - example/plugin/Makefile | 202 ------------------ example/plugin/cmake_install.cmake | 44 ---- 4 files changed, 263 deletions(-) delete mode 100644 example/plugin/CMakeFiles/CMakeDirectoryInformation.cmake delete mode 100644 example/plugin/CMakeFiles/progress.marks delete mode 100644 example/plugin/Makefile delete mode 100644 example/plugin/cmake_install.cmake diff --git a/example/plugin/CMakeFiles/CMakeDirectoryInformation.cmake b/example/plugin/CMakeFiles/CMakeDirectoryInformation.cmake deleted file mode 100644 index d5f4c60f97..0000000000 --- a/example/plugin/CMakeFiles/CMakeDirectoryInformation.cmake +++ /dev/null @@ -1,16 +0,0 @@ -# CMAKE generated file: DO NOT EDIT! -# Generated by "Unix Makefiles" Generator, CMake Version 3.25 - -# Relative path conversion top directories. -set(CMAKE_RELATIVE_PATH_TOP_SOURCE "/home/jonahm/programming/singularity-eos") -set(CMAKE_RELATIVE_PATH_TOP_BINARY "/home/jonahm/programming/singularity-eos/example/plugin") - -# Force unix paths in dependencies. -set(CMAKE_FORCE_UNIX_PATHS 1) - - -# The C and CXX include file regular expressions for this directory. -set(CMAKE_C_INCLUDE_REGEX_SCAN "^.*$") -set(CMAKE_C_INCLUDE_REGEX_COMPLAIN "^$") -set(CMAKE_CXX_INCLUDE_REGEX_SCAN ${CMAKE_C_INCLUDE_REGEX_SCAN}) -set(CMAKE_CXX_INCLUDE_REGEX_COMPLAIN ${CMAKE_C_INCLUDE_REGEX_COMPLAIN}) diff --git a/example/plugin/CMakeFiles/progress.marks b/example/plugin/CMakeFiles/progress.marks deleted file mode 100644 index 573541ac97..0000000000 --- a/example/plugin/CMakeFiles/progress.marks +++ /dev/null @@ -1 +0,0 @@ -0 diff --git a/example/plugin/Makefile b/example/plugin/Makefile deleted file mode 100644 index b4c8c93e1e..0000000000 --- a/example/plugin/Makefile +++ /dev/null @@ -1,202 +0,0 @@ -# CMAKE generated file: DO NOT EDIT! -# Generated by "Unix Makefiles" Generator, CMake Version 3.25 - -# Default target executed when no arguments are given to make. -default_target: all -.PHONY : default_target - -# Allow only one "make -f Makefile2" at a time, but pass parallelism. -.NOTPARALLEL: - -#============================================================================= -# Special targets provided by cmake. - -# Disable implicit rules so canonical targets will work. -.SUFFIXES: - -# Disable VCS-based implicit rules. -% : %,v - -# Disable VCS-based implicit rules. -% : RCS/% - -# Disable VCS-based implicit rules. -% : RCS/%,v - -# Disable VCS-based implicit rules. -% : SCCS/s.% - -# Disable VCS-based implicit rules. -% : s.% - -.SUFFIXES: .hpux_make_needs_suffix_list - -# Command-line flag to silence nested $(MAKE). -$(VERBOSE)MAKESILENT = -s - -#Suppress display of executed commands. -$(VERBOSE).SILENT: - -# A target that is always out of date. -cmake_force: -.PHONY : cmake_force - -#============================================================================= -# Set environment variables for the build. - -# The shell in which to execute make rules. -SHELL = /bin/sh - -# The CMake executable. -CMAKE_COMMAND = /usr/bin/cmake - -# The command to remove a file. -RM = /usr/bin/cmake -E rm -f - -# Escaping for special characters. -EQUALS = = - -# The top-level source directory on which CMake was run. -CMAKE_SOURCE_DIR = /home/jonahm/programming/singularity-eos - -# The top-level build directory on which CMake was run. -CMAKE_BINARY_DIR = /home/jonahm/programming/singularity-eos/build - -#============================================================================= -# Targets provided globally by CMake. - -# Special rule for the target test -test: - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running tests..." - /usr/bin/ctest --force-new-ctest-process $(ARGS) -.PHONY : test - -# Special rule for the target test -test/fast: test -.PHONY : test/fast - -# Special rule for the target edit_cache -edit_cache: - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake cache editor..." - /usr/bin/ccmake -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) -.PHONY : edit_cache - -# Special rule for the target edit_cache -edit_cache/fast: edit_cache -.PHONY : edit_cache/fast - -# Special rule for the target rebuild_cache -rebuild_cache: - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake to regenerate build system..." - /usr/bin/cmake --regenerate-during-build -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) -.PHONY : rebuild_cache - -# Special rule for the target rebuild_cache -rebuild_cache/fast: rebuild_cache -.PHONY : rebuild_cache/fast - -# Special rule for the target list_install_components -list_install_components: - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Available install components are: \"Unspecified\"" -.PHONY : list_install_components - -# Special rule for the target list_install_components -list_install_components/fast: list_install_components -.PHONY : list_install_components/fast - -# Special rule for the target install -install: preinstall - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Install the project..." - /usr/bin/cmake -P cmake_install.cmake -.PHONY : install - -# Special rule for the target install -install/fast: preinstall/fast - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Install the project..." - /usr/bin/cmake -P cmake_install.cmake -.PHONY : install/fast - -# Special rule for the target install/local -install/local: preinstall - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Installing only the local directory..." - /usr/bin/cmake -DCMAKE_INSTALL_LOCAL_ONLY=1 -P cmake_install.cmake -.PHONY : install/local - -# Special rule for the target install/local -install/local/fast: preinstall/fast - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Installing only the local directory..." - /usr/bin/cmake -DCMAKE_INSTALL_LOCAL_ONLY=1 -P cmake_install.cmake -.PHONY : install/local/fast - -# Special rule for the target install/strip -install/strip: preinstall - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Installing the project stripped..." - /usr/bin/cmake -DCMAKE_INSTALL_DO_STRIP=1 -P cmake_install.cmake -.PHONY : install/strip - -# Special rule for the target install/strip -install/strip/fast: preinstall/fast - @$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Installing the project stripped..." - /usr/bin/cmake -DCMAKE_INSTALL_DO_STRIP=1 -P cmake_install.cmake -.PHONY : install/strip/fast - -# The main all target -all: cmake_check_build_system - cd /home/jonahm/programming/singularity-eos/build && $(CMAKE_COMMAND) -E cmake_progress_start /home/jonahm/programming/singularity-eos/build/CMakeFiles /home/jonahm/programming/singularity-eos/example/plugin//CMakeFiles/progress.marks - cd /home/jonahm/programming/singularity-eos/build && $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 /home/jonahm/programming/singularity-eos/example/plugin/all - $(CMAKE_COMMAND) -E cmake_progress_start /home/jonahm/programming/singularity-eos/build/CMakeFiles 0 -.PHONY : all - -# The main clean target -clean: - cd /home/jonahm/programming/singularity-eos/build && $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 /home/jonahm/programming/singularity-eos/example/plugin/clean -.PHONY : clean - -# The main clean target -clean/fast: clean -.PHONY : clean/fast - -# Prepare targets for installation. -preinstall: all - cd /home/jonahm/programming/singularity-eos/build && $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 /home/jonahm/programming/singularity-eos/example/plugin/preinstall -.PHONY : preinstall - -# Prepare targets for installation. -preinstall/fast: - cd /home/jonahm/programming/singularity-eos/build && $(MAKE) $(MAKESILENT) -f CMakeFiles/Makefile2 /home/jonahm/programming/singularity-eos/example/plugin/preinstall -.PHONY : preinstall/fast - -# clear depends -depend: - cd /home/jonahm/programming/singularity-eos/build && $(CMAKE_COMMAND) -P /home/jonahm/programming/singularity-eos/build/CMakeFiles/VerifyGlobs.cmake - cd /home/jonahm/programming/singularity-eos/build && $(CMAKE_COMMAND) -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) --check-build-system CMakeFiles/Makefile.cmake 1 -.PHONY : depend - -# Help Target -help: - @echo "The following are some of the valid targets for this Makefile:" - @echo "... all (the default if no target is provided)" - @echo "... clean" - @echo "... depend" - @echo "... edit_cache" - @echo "... install" - @echo "... install/local" - @echo "... install/strip" - @echo "... list_install_components" - @echo "... rebuild_cache" - @echo "... test" -.PHONY : help - - - -#============================================================================= -# Special targets to cleanup operation of make. - -# Special rule to run CMake to check the build system integrity. -# No rule that depends on this can have commands that come from listfiles -# because they might be regenerated. -cmake_check_build_system: - cd /home/jonahm/programming/singularity-eos/build && $(CMAKE_COMMAND) -P /home/jonahm/programming/singularity-eos/build/CMakeFiles/VerifyGlobs.cmake - cd /home/jonahm/programming/singularity-eos/build && $(CMAKE_COMMAND) -S$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR) --check-build-system CMakeFiles/Makefile.cmake 0 -.PHONY : cmake_check_build_system - diff --git a/example/plugin/cmake_install.cmake b/example/plugin/cmake_install.cmake deleted file mode 100644 index 88c1fd1691..0000000000 --- a/example/plugin/cmake_install.cmake +++ /dev/null @@ -1,44 +0,0 @@ -# Install script for directory: /home/jonahm/programming/singularity-eos/example/plugin - -# Set the install prefix -if(NOT DEFINED CMAKE_INSTALL_PREFIX) - set(CMAKE_INSTALL_PREFIX "/home/jonahm/test") -endif() -string(REGEX REPLACE "/$" "" CMAKE_INSTALL_PREFIX "${CMAKE_INSTALL_PREFIX}") - -# Set the install configuration name. -if(NOT DEFINED CMAKE_INSTALL_CONFIG_NAME) - if(BUILD_TYPE) - string(REGEX REPLACE "^[^A-Za-z0-9_]+" "" - CMAKE_INSTALL_CONFIG_NAME "${BUILD_TYPE}") - else() - set(CMAKE_INSTALL_CONFIG_NAME "") - endif() - message(STATUS "Install configuration: \"${CMAKE_INSTALL_CONFIG_NAME}\"") -endif() - -# Set the component getting installed. -if(NOT CMAKE_INSTALL_COMPONENT) - if(COMPONENT) - message(STATUS "Install component: \"${COMPONENT}\"") - set(CMAKE_INSTALL_COMPONENT "${COMPONENT}") - else() - set(CMAKE_INSTALL_COMPONENT) - endif() -endif() - -# Install shared libraries without execute permission? -if(NOT DEFINED CMAKE_INSTALL_SO_NO_EXE) - set(CMAKE_INSTALL_SO_NO_EXE "1") -endif() - -# Is this installation the result of a crosscompile? -if(NOT DEFINED CMAKE_CROSSCOMPILING) - set(CMAKE_CROSSCOMPILING "FALSE") -endif() - -# Set default install directory permissions. -if(NOT DEFINED CMAKE_OBJDUMP) - set(CMAKE_OBJDUMP "/usr/bin/objdump") -endif() - From 6c43d6971f2b236e1143430d52561b13ac6b961c Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 10 Jan 2024 17:25:25 -0700 Subject: [PATCH 027/405] bin dir should be in build directory --- CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7452b4d4ff..86f6dd96d6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -375,7 +375,8 @@ add_subdirectory(singularity-eos) foreach(_plugin ${SINGULARITY_PLUGINS}) message(STATUS "Adding plugin ${_plugin}...") - add_subdirectory(${_plugin} ${_plugin}) + get_filename_component(BINDIR_LOC ${_plugin} NAME) + add_subdirectory(${_plugin} ${CMAKE_CURRENT_BINARY_DIR}/${BINDIR_LOC}) endforeach() # TODO(JMM): Kind of nice to have? From d508698e66ead3e8ce30680944a0a962e5b3bd8a Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 11 Jan 2024 12:07:22 -0700 Subject: [PATCH 028/405] make include directories and install directories work... I think --- CMakeLists.txt | 22 +++++++++++++++++----- cmake/install.cmake | 4 ++-- cmake/plugins.cmake | 12 +++++++++++- example/plugin/CMakeLists.txt | 4 ++-- example/plugin/{ => dust}/dust.hpp | 14 +++++++------- example/plugin/{ => tst}/test_dust.cpp | 2 +- singularity-eos/CMakeLists.txt | 2 +- 7 files changed, 41 insertions(+), 19 deletions(-) rename example/plugin/{ => dust}/dust.hpp (95%) rename example/plugin/{ => tst}/test_dust.cpp (98%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 86f6dd96d6..bafe34bec6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -121,6 +121,9 @@ option(SINGULARITY_PATCH_MPARK_VARIANT # Plugins set(SINGULARITY_PLUGINS "" CACHE STRING "List of paths to plugin directories") +# CMake specific options +option(SINGULARITY_VERBOSE_CONFIG "Have cmake output additional information" OFF) + # ------------------------------------------------------------------------------# # singularity-eos Library # ------------------------------------------------------------------------------# @@ -362,6 +365,12 @@ if(SINGULARITY_BUILD_TESTS) endif() endif() +# ------------------------------------------------------------------------------# +# Plugin infrastructure +# ------------------------------------------------------------------------------# +include(cmake/plugins.cmake) +set(PLUGIN_BUILD_ROOT "${CMAKE_CURRENT_BINARY_DIR}/plugins") + # ------------------------------------------------------------------------------# # singularity-eos library # ------------------------------------------------------------------------------# @@ -369,19 +378,19 @@ endif() # Here we populate `EOS_HEADERS/EOS_SRCS` NOTE: these include path # prefixes of subdirectories on files (e.g. eos/eos.hpp) see # singularity-eos/CMakeLists.txt -include(cmake/plugins.cmake) - add_subdirectory(singularity-eos) foreach(_plugin ${SINGULARITY_PLUGINS}) message(STATUS "Adding plugin ${_plugin}...") get_filename_component(BINDIR_LOC ${_plugin} NAME) - add_subdirectory(${_plugin} ${CMAKE_CURRENT_BINARY_DIR}/${BINDIR_LOC}) + add_subdirectory(${_plugin} ${PLUGIN_BUILD_ROOT}/${BINDIR_LOC}) endforeach() # TODO(JMM): Kind of nice to have? -message(STATUS "EOS Headers:\n\t${EOS_HEADERS}") -message(STATUS "EOS Sources:\n\t${EOS_SRCS}") +if (SINGULARITY_VERBOSE_CONFIG) + message(STATUS "EOS Headers:\n\t${EOS_HEADERS}") + message(STATUS "EOS Sources:\n\t${EOS_SRCS}") +endif() target_sources(singularity-eos PRIVATE ${EOS_SRCS} ${EOS_HEADERS}) @@ -400,6 +409,9 @@ endif() # SINGULARITY_USE_FORTRAN target_include_directories( singularity-eos PUBLIC $ $) +foreach(path ${PLUGIN_INCLUDE_PATHS}) + target_include_directories(singularity-eos PUBLIC $) +endforeach() # plug in collected includes/libs/definitions diff --git a/cmake/install.cmake b/cmake/install.cmake index 3067aa5ed1..679c2f2671 100644 --- a/cmake/install.cmake +++ b/cmake/install.cmake @@ -72,8 +72,8 @@ install( # install singularity-eos headers foreach(file ${_install_headers}) get_filename_component(DIR ${file} DIRECTORY) - install(FILES singularity-eos/${file} - DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/singularity-eos/${DIR}) + install(FILES ${file} + DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${DIR}) endforeach() # file # install the fortran modules NB: cmake doesn't provide a clean way to handle diff --git a/cmake/plugins.cmake b/cmake/plugins.cmake index 11ad8a6ce9..4a26384fdf 100644 --- a/cmake/plugins.cmake +++ b/cmake/plugins.cmake @@ -18,7 +18,7 @@ function(register_headers) foreach(arg IN LISTS ARGN) file(RELATIVE_PATH relative_path ${CMAKE_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}) list(APPEND EOS_HEADERS ${relative_path}/${arg}) - list(APPEND _install_headers ${arg}) + list(APPEND _install_headers ${relative_path}/${arg}) endforeach() set(EOS_HEADERS ${EOS_HEADERS} PARENT_SCOPE) set(_install_headers ${_install_headers} PARENT_SCOPE) @@ -41,9 +41,19 @@ function(register_tests) set(PLUGIN_TESTS ${PLUGIN_TESTS} PARENT_SCOPE) endfunction() +set(PLUGIN_INCLUDE_PATHS "") +macro(export_source) + set(EOS_HEADERS ${EOS_HEADERS} PARENT_SCOPE) + set(_install_headers ${_install_headers} PARENT_SCOPE) + set(EOS_SRCS ${EOS_SRCS} PARENT_SCOPE) + set(PLUGIN_TESTS ${PLUGIN_TESTS} PARENT_SCOPE) + set(PLUGIN_INCLUDE_PATHS ${PLUGIN_INCLUDE_PATHS} PARENT_SCOPE) +endmacro() macro(export_plugin) + list(APPEND PLUGIN_INCLUDE_PATHS ${CMAKE_CURRENT_SOURCE_DIR}) set(EOS_HEADERS ${EOS_HEADERS} PARENT_SCOPE) set(_install_headers ${_install_headers} PARENT_SCOPE) set(EOS_SRCS ${EOS_SRCS} PARENT_SCOPE) set(PLUGIN_TESTS ${PLUGIN_TESTS} PARENT_SCOPE) + set(PLUGIN_INCLUDE_PATHS ${PLUGIN_INCLUDE_PATHS} PARENT_SCOPE) endmacro() diff --git a/example/plugin/CMakeLists.txt b/example/plugin/CMakeLists.txt index 16f4616c6b..03984c8f63 100644 --- a/example/plugin/CMakeLists.txt +++ b/example/plugin/CMakeLists.txt @@ -13,10 +13,10 @@ #------------------------------------------------------------------------------# # Add the header file to the list of headers -register_headers(dust.hpp) +register_headers(dust/dust.hpp) # Register this test to be built in the test suite -register_tests(test_dust.cpp) +register_tests(tst/test_dust.cpp) # Export lists to parent scope export_plugin() diff --git a/example/plugin/dust.hpp b/example/plugin/dust/dust.hpp similarity index 95% rename from example/plugin/dust.hpp rename to example/plugin/dust/dust.hpp index e99828661f..03c737d4ec 100644 --- a/example/plugin/dust.hpp +++ b/example/plugin/dust/dust.hpp @@ -26,7 +26,7 @@ namespace singularity { -using namespace eos_base; +using namespace eos_base; class Dust : public EosBase { public: @@ -44,9 +44,9 @@ class Dust : public EosBase { const Real rho, const Real temperature, Real *lambda = nullptr) { return _Cv * temperature; } - PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( + PORTABLE_INLINE_FUNCTION Real PressureFromDensityTemperature( const Real rho, const Real temperature, Real *lambda = nullptr) const { - return 0; + return 0; } PORTABLE_INLINE_FUNCTION Real PressureFromDensityInternalEnergy( const Real rho, const Real sie, Real *lambda = nullptr) const { @@ -103,7 +103,7 @@ class Dust : public EosBase { if (!(output & thermalqs::temperature)) { UNDEFINED_ERROR; } - energy = _Cv*temp; + energy = _Cv * temp; } if (output & thermalqs::temperature) { if (!(output & thermalqs::specific_internal_energy)) { @@ -120,7 +120,7 @@ class Dust : public EosBase { dvdt = 0; } SG_ADD_BASE_CLASS_USINGS(Dust) - + PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 0; } static constexpr unsigned long PreferredInput() { return _preferred_input; } static inline unsigned long scratch_size(std::string method, unsigned int nelements) { @@ -133,10 +133,10 @@ class Dust : public EosBase { static std::string EosType() { return std::string("Dust"); } static std::string EosPyType() { return EosType(); } -private: + private: Real _Cv = 1; static constexpr const unsigned long _preferred_input = - thermalqs::density | thermalqs::specific_internal_energy | thermalqs::temperature; + thermalqs::density | thermalqs::specific_internal_energy | thermalqs::temperature; }; } // namespace singularity diff --git a/example/plugin/test_dust.cpp b/example/plugin/tst/test_dust.cpp similarity index 98% rename from example/plugin/test_dust.cpp rename to example/plugin/tst/test_dust.cpp index bef73438a6..19bdd8ddc0 100644 --- a/example/plugin/test_dust.cpp +++ b/example/plugin/tst/test_dust.cpp @@ -31,8 +31,8 @@ #include #endif +#include #include -#include "dust.hpp" using singularity::Dust; using EOS = singularity::Variant; diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 82deec6e2c..87f73c708d 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -71,4 +71,4 @@ if (SINGULARITY_USE_FORTRAN) endif() # export to parent scope -export_plugin() +export_source() From 57ada428a9085f112174d86c77bb25d008402bcb Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 11 Jan 2024 12:56:05 -0700 Subject: [PATCH 029/405] fix installs by adding PLUGIN named argument --- cmake/install.cmake | 10 +++++++--- cmake/plugins.cmake | 15 ++++++++++++--- example/plugin/CMakeLists.txt | 2 +- 3 files changed, 20 insertions(+), 7 deletions(-) diff --git a/cmake/install.cmake b/cmake/install.cmake index 679c2f2671..77fabcd491 100644 --- a/cmake/install.cmake +++ b/cmake/install.cmake @@ -70,9 +70,13 @@ install( # ----------------------------------------------------------------------------# # install singularity-eos headers -foreach(file ${_install_headers}) - get_filename_component(DIR ${file} DIRECTORY) - install(FILES ${file} +list(LENGTH _install_headers length) +math(EXPR max_index "${length} - 1") +foreach(index RANGE ${max_index}) + list(GET EOS_HEADERS ${index} src) + list(GET _install_headers ${index} dst) + get_filename_component(DIR ${dst} DIRECTORY) + install(FILES ${src} DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${DIR}) endforeach() # file diff --git a/cmake/plugins.cmake b/cmake/plugins.cmake index 4a26384fdf..5cb836224b 100644 --- a/cmake/plugins.cmake +++ b/cmake/plugins.cmake @@ -15,10 +15,19 @@ set(EOS_HEADERS "") set(_install_headers "") function(register_headers) - foreach(arg IN LISTS ARGN) + set(keyword_args PLUGIN) + cmake_parse_arguments(ARG "" "${keyword_args}" "" ${ARGN}) + set(variadic_args ${ARG_UNPARSED_ARGUMENTS}) + + foreach(arg IN LISTS variadic_args) file(RELATIVE_PATH relative_path ${CMAKE_SOURCE_DIR} ${CMAKE_CURRENT_SOURCE_DIR}) - list(APPEND EOS_HEADERS ${relative_path}/${arg}) - list(APPEND _install_headers ${relative_path}/${arg}) + if (ARG_PLUGIN) + list(APPEND EOS_HEADERS ${relative_path}/${ARG_PLUGIN}/${arg}) + list(APPEND _install_headers ${ARG_PLUGIN}/${arg}) + else() + list(APPEND EOS_HEADERS ${relative_path}/${arg}) + list(APPEND _install_headers ${relative_path}/${arg}) + endif() endforeach() set(EOS_HEADERS ${EOS_HEADERS} PARENT_SCOPE) set(_install_headers ${_install_headers} PARENT_SCOPE) diff --git a/example/plugin/CMakeLists.txt b/example/plugin/CMakeLists.txt index 03984c8f63..8b9efa1f7f 100644 --- a/example/plugin/CMakeLists.txt +++ b/example/plugin/CMakeLists.txt @@ -13,7 +13,7 @@ #------------------------------------------------------------------------------# # Add the header file to the list of headers -register_headers(dust/dust.hpp) +register_headers(PLUGIN dust dust.hpp) # Register this test to be built in the test suite register_tests(tst/test_dust.cpp) From 29a18d7ec83c3442dcba25091d35c5c8f7381f16 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 11 Jan 2024 17:31:18 -0700 Subject: [PATCH 030/405] add ability to specify variant with cmake filegen --- CMakeLists.txt | 6 + cmake/install.cmake | 7 + example/plugin/CMakeLists.txt | 2 +- example/plugin/dust/dust_variant.hpp | 144 ++++++++++++++++++ example/plugin/tst/test_dust.cpp | 1 + singularity-eos/CMakeLists.txt | 14 +- .../eos/{eos.hpp => default_variant.hpp} | 28 +--- singularity-eos/eos/eos.hpp.in | 20 +++ singularity-eos/eos/eos_models.hpp | 39 +++++ 9 files changed, 236 insertions(+), 25 deletions(-) create mode 100644 example/plugin/dust/dust_variant.hpp rename singularity-eos/eos/{eos.hpp => default_variant.hpp} (83%) create mode 100644 singularity-eos/eos/eos.hpp.in create mode 100644 singularity-eos/eos/eos_models.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index bafe34bec6..6e262e8c95 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -120,6 +120,8 @@ option(SINGULARITY_PATCH_MPARK_VARIANT # Plugins set(SINGULARITY_PLUGINS "" CACHE STRING "List of paths to plugin directories") +set(SINGULARITY_VARIANT "singularity-eos/eos/default_variant.hpp" CACHE STRING + "The include path for the file containing the definition of singularity::EOS.") # CMake specific options option(SINGULARITY_VERBOSE_CONFIG "Have cmake output additional information" OFF) @@ -409,6 +411,10 @@ endif() # SINGULARITY_USE_FORTRAN target_include_directories( singularity-eos PUBLIC $ $) +target_include_directories( + singularity-eos PUBLIC $) + + foreach(path ${PLUGIN_INCLUDE_PATHS}) target_include_directories(singularity-eos PUBLIC $) endforeach() diff --git a/cmake/install.cmake b/cmake/install.cmake index 77fabcd491..cac21733b3 100644 --- a/cmake/install.cmake +++ b/cmake/install.cmake @@ -80,6 +80,13 @@ foreach(index RANGE ${max_index}) DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/${DIR}) endforeach() # file +# Special sauce so generated file has proper include path +install(FILES + ${CMAKE_BINARY_DIR}/generated/singularity-eos/eos/eos.hpp + DESTINATION + ${CMAKE_INSTALL_INCLUDEDIR}/singularity-eos/eos + ) + # install the fortran modules NB: cmake doesn't provide a clean way to handle if(SINGULARITY_USE_FORTRAN) install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/fortran/ diff --git a/example/plugin/CMakeLists.txt b/example/plugin/CMakeLists.txt index 8b9efa1f7f..c7f4340673 100644 --- a/example/plugin/CMakeLists.txt +++ b/example/plugin/CMakeLists.txt @@ -13,7 +13,7 @@ #------------------------------------------------------------------------------# # Add the header file to the list of headers -register_headers(PLUGIN dust dust.hpp) +register_headers(PLUGIN dust dust.hpp dust_variant.hpp) # Register this test to be built in the test suite register_tests(tst/test_dust.cpp) diff --git a/example/plugin/dust/dust_variant.hpp b/example/plugin/dust/dust_variant.hpp new file mode 100644 index 0000000000..cba779a537 --- /dev/null +++ b/example/plugin/dust/dust_variant.hpp @@ -0,0 +1,144 @@ +//------------------------------------------------------------------------------ +// © 2024. 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 _EXAMPLE_PLUGIN_DUST_DUST_VARIANT_HPP_ +#define _EXAMPLE_PLUGIN_DUST_DUST_VARIANT_HPP_ + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +// Base stuff +#include +#include +#include + +// EOS models +#include + +#include + +namespace singularity { + +// recreate variadic list +template +using tl = singularity::variadic_utils::type_list; + +template