From 0e9147e7691bbe05081be2733a3b11cc7db7b786 Mon Sep 17 00:00:00 2001 From: Benjamin Joel Musick Date: Thu, 10 Aug 2023 14:53:04 -0600 Subject: [PATCH 01/20] 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 02/20] 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 03/20] 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 04/20] 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 05/20] 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 06/20] 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 07/20] 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 08/20] 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 09/20] 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 d59ba1857e97b2bb09d5bbe6873ee85144f436ca Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 30 May 2024 19:07:47 -0600 Subject: [PATCH 10/20] Add minimum energy --- singularity-eos/eos/eos_carnahan_starling.hpp | 4 ++++ singularity-eos/eos/eos_noble_abel.hpp | 3 +-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_carnahan_starling.hpp b/singularity-eos/eos/eos_carnahan_starling.hpp index 4e242b5712..9e3147d80f 100644 --- a/singularity-eos/eos/eos_carnahan_starling.hpp +++ b/singularity-eos/eos/eos_carnahan_starling.hpp @@ -170,6 +170,10 @@ class CarnahanStarling : public EosBase { const Real rho, const Real sie, Real *lambda = nullptr) const { return ZedFromDensity(rho) * _gm1; } + PORTABLE_INLINE_FUNCTION Real + MinInternalEnergyFromDensity(const Real rho, Real *lambda = nullptr) const { + return _qq; + } 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; diff --git a/singularity-eos/eos/eos_noble_abel.hpp b/singularity-eos/eos/eos_noble_abel.hpp index fc2edd49c4..c4550e1612 100644 --- a/singularity-eos/eos/eos_noble_abel.hpp +++ b/singularity-eos/eos/eos_noble_abel.hpp @@ -82,8 +82,7 @@ class NobleAbel : public EosBase { } PORTABLE_INLINE_FUNCTION Real MinInternalEnergyFromDensity(const Real rho, Real *lambda = nullptr) const { - MinInternalEnergyIsNotEnabled("Noble Abel"); - return 0.0; + return _qq; } PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature( From 8c44ab4874eee7f1b43ff3c2a823c26382fbff58 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 30 May 2024 19:26:45 -0600 Subject: [PATCH 11/20] Add carnahan starling EOS to models --- singularity-eos/eos/eos_models.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/singularity-eos/eos/eos_models.hpp b/singularity-eos/eos/eos_models.hpp index 989afebd55..ebdc0de1f7 100644 --- a/singularity-eos/eos/eos_models.hpp +++ b/singularity-eos/eos/eos_models.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include From 24697262c136f63f3d98b3831e73ccbe020f8b16 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 30 May 2024 19:28:24 -0600 Subject: [PATCH 12/20] Clang format --- singularity-eos/eos/eos_models.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_models.hpp b/singularity-eos/eos/eos_models.hpp index ebdc0de1f7..306de6fb28 100644 --- a/singularity-eos/eos/eos_models.hpp +++ b/singularity-eos/eos/eos_models.hpp @@ -16,6 +16,7 @@ #define _SINGULARITY_EOS_EOS_EOS_MODELS_HPP_ // EOS models +#include #include #include #include @@ -24,7 +25,6 @@ #include #include #include -#include #include #include #include From f7efe0d53aadf9a033d98d62e82bcbba8858e308 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 30 May 2024 19:30:32 -0600 Subject: [PATCH 13/20] Upgrade Catch --- test/test_eos_carnahan_starling.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_eos_carnahan_starling.cpp b/test/test_eos_carnahan_starling.cpp index efca42313a..191a423dcc 100644 --- a/test/test_eos_carnahan_starling.cpp +++ b/test/test_eos_carnahan_starling.cpp @@ -17,8 +17,9 @@ #include #include #include -#ifndef CATCH_CONFIG_RUNNER -#include "catch2/catch.hpp" +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include #endif #include From 70d341a4e469dab51a5de0c5b4b07c609d5ba75b Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 30 May 2024 19:35:33 -0600 Subject: [PATCH 14/20] Fix bug in EOS variant specification --- test/test_eos_carnahan_starling.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_eos_carnahan_starling.cpp b/test/test_eos_carnahan_starling.cpp index 191a423dcc..17855ce8aa 100644 --- a/test/test_eos_carnahan_starling.cpp +++ b/test/test_eos_carnahan_starling.cpp @@ -26,8 +26,9 @@ #include #include +using singularity::IdealGas; using singularity::CarnahanStarling; -using EOS = singularity::Variant; +using EOS = singularity::Variant; SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") { GIVEN("Parameters for a CarnahanStarling EOS") { From 4f0ad14ed5b2927804a1779d7e1638b73918ae8a Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 30 May 2024 19:39:37 -0600 Subject: [PATCH 15/20] Copy/paste typo --- test/test_eos_carnahan_starling.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_eos_carnahan_starling.cpp b/test/test_eos_carnahan_starling.cpp index 17855ce8aa..169d9434c2 100644 --- a/test/test_eos_carnahan_starling.cpp +++ b/test/test_eos_carnahan_starling.cpp @@ -28,7 +28,7 @@ using singularity::IdealGas; using singularity::CarnahanStarling; -using EOS = singularity::Variant; +using EOS = singularity::Variant; SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") { GIVEN("Parameters for a CarnahanStarling EOS") { From a93135fce2c8f38238cb50740452efabef82db9c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 30 May 2024 19:43:10 -0600 Subject: [PATCH 16/20] Clang format --- test/test_eos_carnahan_starling.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_eos_carnahan_starling.cpp b/test/test_eos_carnahan_starling.cpp index 169d9434c2..1c9a610031 100644 --- a/test/test_eos_carnahan_starling.cpp +++ b/test/test_eos_carnahan_starling.cpp @@ -26,8 +26,8 @@ #include #include -using singularity::IdealGas; using singularity::CarnahanStarling; +using singularity::IdealGas; using EOS = singularity::Variant; SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") { From a81a8572651efa4a112a8ffa477d2fcde49f216b Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 30 May 2024 19:56:21 -0600 Subject: [PATCH 17/20] Add to Carnahan-Staring documentation --- doc/sphinx/src/models.rst | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index d299737133..205aa421c6 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -556,21 +556,36 @@ the form .. math:: - Z(\rho) = \frac{1+\eta+\eta^2-\eta^3}{(1-\eta)^3} + Z(\rho) = \frac{1+\eta+\eta^2-\eta^3}{(1-\eta)^3}, + +where :math:`\eta` is the packing fraction given by .. math:: - \eta = b\rho + \eta = b\rho. + +The energy is related to the temperature through .. math:: e = C_V T + q, -the inputs for the equation of state are similar to the Noble-Abel model. +where :math:`q` is an energy offset. + 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 Carnahan-Starling EOS is intended to represent a hard sphere fluid, and the +covolume parameter, :math:`b`, can be related to the hard sphere +diameter, :math:`\sigma`, through + +.. math:: + + b = \frac{\pi}{6}\frac{\sigma^3}{M}, + +where :math:`M` is the molar mass of the gas. + The entropy for the Carnahan-Starling EOS is given by .. math:: From 0e9b8a496105edab2a335dc23a9cc8eccfaccde3 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 31 May 2024 08:41:27 -0600 Subject: [PATCH 18/20] Correct copy/paste errors --- test/test_eos_carnahan_starling.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/test/test_eos_carnahan_starling.cpp b/test/test_eos_carnahan_starling.cpp index 1c9a610031..78a63671a3 100644 --- a/test/test_eos_carnahan_starling.cpp +++ b/test/test_eos_carnahan_starling.cpp @@ -535,7 +535,6 @@ SCENARIO("Test C-S Entropy Calls", "[CarnahanStarling][CarnahanStarling5]") { #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}; @@ -631,8 +630,8 @@ SCENARIO("CarnahanStarling6", "[CarnahanStarling][CarnahanStarling6]") { // 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); + Kokkos::deep_copy(pressure, v_pressure); + Kokkos::deep_copy(temperature, v_temperature); #endif // PORTABILITY_STRATEGY_KOKKOS // Gold standard values for a subset of lookups From a56df1fdfe9570765b7097053ff28c9ff910c0b9 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 31 May 2024 10:59:02 -0600 Subject: [PATCH 19/20] Change for to portableFor --- test/test_eos_carnahan_starling.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_eos_carnahan_starling.cpp b/test/test_eos_carnahan_starling.cpp index 78a63671a3..b3f024ca3d 100644 --- a/test/test_eos_carnahan_starling.cpp +++ b/test/test_eos_carnahan_starling.cpp @@ -656,14 +656,14 @@ SCENARIO("CarnahanStarling6", "[CarnahanStarling][CarnahanStarling6]") { #endif // PORTABILITY_STRATEGY_KOKKOS WHEN("A rho(P, T) lookup is performed") { - for (int i = 0; i < num; i++) { + portableFor("rho(P, T) FillEos lookup", 0, num, PORTABLE_LAMBDA(int 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(); From 6af50ffeed81dc1c83515f15abbff741e3451ed0 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 31 May 2024 10:59:39 -0600 Subject: [PATCH 20/20] Clang format --- test/test_eos_carnahan_starling.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/test/test_eos_carnahan_starling.cpp b/test/test_eos_carnahan_starling.cpp index b3f024ca3d..fddba3b14b 100644 --- a/test/test_eos_carnahan_starling.cpp +++ b/test/test_eos_carnahan_starling.cpp @@ -656,14 +656,15 @@ SCENARIO("CarnahanStarling6", "[CarnahanStarling][CarnahanStarling6]") { #endif // PORTABILITY_STRATEGY_KOKKOS WHEN("A rho(P, T) lookup is performed") { - portableFor("rho(P, T) FillEos lookup", 0, num, PORTABLE_LAMBDA(int 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); - }); + portableFor( + "rho(P, T) FillEos lookup", 0, num, PORTABLE_LAMBDA(int 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();