From 8c947a31103dbab074b256a14e6f5b2e32d9fdd5 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 11 Oct 2023 19:15:22 -0600 Subject: [PATCH 1/7] split test_eos_unit into multiple, easier to read, faster to compile files --- test/CMakeLists.txt | 7 +- test/eos_unit_test_helpers.hpp | 78 +- test/test_eos_ideal.cpp | 144 +++ test/test_eos_modifiers.cpp | 276 ++++++ test/test_eos_stellar_collapse.cpp | 299 ++++++ test/test_eos_tabulated.cpp | 314 ++++++ test/test_eos_unit.cpp | 1441 ---------------------------- test/test_eos_vector.cpp | 387 ++++++++ test/test_math_utils.cpp | 128 +++ 9 files changed, 1609 insertions(+), 1465 deletions(-) create mode 100644 test/test_eos_ideal.cpp create mode 100644 test/test_eos_modifiers.cpp create mode 100644 test/test_eos_stellar_collapse.cpp create mode 100644 test/test_eos_tabulated.cpp delete mode 100644 test/test_eos_unit.cpp create mode 100644 test/test_eos_vector.cpp create mode 100644 test/test_math_utils.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index db952eda52..5f71e9a5a4 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -15,13 +15,18 @@ add_executable( eos_unit_tests catch2_define.cpp eos_unit_test_helpers.hpp - test_eos_unit.cpp + test_eos_ideal.cpp test_eos_gruneisen.cpp test_eos_sap_polynomial.cpp test_eos_vinet.cpp test_eos_noble_abel.cpp test_eos_stiff.cpp test_eos_helmholtz.cpp + test_eos_modifiers.cpp + test_eos_vector.cpp + test_eos_tabulated.cpp + test_eos_stellar_collapse.cpp + test_math_utils.cpp ) if(SINGULARITY_TEST_HELMHOLTZ) diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index edd3f5279b..cda089b013 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -49,6 +49,7 @@ inline std::string demangle(const char *name) { return name; } #endif #include +#include using singularity::EOS; @@ -71,34 +72,65 @@ inline void array_compare(int num, X &&x, Y &&y, Z &&z, ZT &&ztrue, XN xname, YN } } -inline void compare_two_eoss(const EOS &test_e, const EOS &ref_e) { +template +inline void compare_two_eoss(const EOS &&test_e, const EOS &&ref_e) { // compare all individual member functions with 1 as inputs, // this function is meant to catch mis-implementations of // modifiers that can be initialized in such a way as to // be equivalent of an unmodified eos. Best used with analytic // eoss. - REQUIRE(isClose(test_e.TemperatureFromDensityInternalEnergy(1, 1), - ref_e.TemperatureFromDensityInternalEnergy(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.InternalEnergyFromDensityTemperature(1, 1), - ref_e.InternalEnergyFromDensityTemperature(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.PressureFromDensityInternalEnergy(1, 1), - ref_e.PressureFromDensityInternalEnergy(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.SpecificHeatFromDensityInternalEnergy(1, 1), - ref_e.SpecificHeatFromDensityInternalEnergy(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.BulkModulusFromDensityInternalEnergy(1, 1), - ref_e.BulkModulusFromDensityInternalEnergy(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.GruneisenParamFromDensityInternalEnergy(1, 1), - ref_e.GruneisenParamFromDensityInternalEnergy(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.PressureFromDensityTemperature(1, 1), - ref_e.PressureFromDensityTemperature(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.SpecificHeatFromDensityTemperature(1, 1), - ref_e.SpecificHeatFromDensityTemperature(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.BulkModulusFromDensityTemperature(1, 1), - ref_e.BulkModulusFromDensityTemperature(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.GruneisenParamFromDensityTemperature(1, 1), - ref_e.GruneisenParamFromDensityTemperature(1, 1), 1.e-15)); - REQUIRE(isClose(test_e.MinimumDensity(), ref_e.MinimumDensity(), 1.e-15)); - REQUIRE(isClose(test_e.MinimumTemperature(), ref_e.MinimumTemperature(), 1.e-15)); + INFO("reference T: " << ref_e.TemperatureFromDensityInternalEnergy(1, 1) << " test T: " + << test_e.TemperatureFromDensityInternalEnergy(1, 1)); + CHECK(isClose(test_e.TemperatureFromDensityInternalEnergy(1, 1), + ref_e.TemperatureFromDensityInternalEnergy(1, 1), 1.e-15)); + INFO("reference sie: " << ref_e.InternalEnergyFromDensityTemperature(1, 1) + << " test sie: " + << test_e.InternalEnergyFromDensityTemperature(1, 1)); + CHECK(isClose(test_e.InternalEnergyFromDensityTemperature(1, 1), + ref_e.InternalEnergyFromDensityTemperature(1, 1), 1.e-15)); + INFO("reference P: " << ref_e.PressureFromDensityInternalEnergy(1, 1) + << " test P: " << test_e.PressureFromDensityInternalEnergy(1, 1)); + CHECK(isClose(test_e.PressureFromDensityInternalEnergy(1, 1), + ref_e.PressureFromDensityInternalEnergy(1, 1), 1.e-15)); + INFO("reference Cv: " << ref_e.SpecificHeatFromDensityInternalEnergy(1, 1) + << " test Cv: " + << test_e.SpecificHeatFromDensityInternalEnergy(1, 1)); + CHECK(isClose(test_e.SpecificHeatFromDensityInternalEnergy(1, 1), + ref_e.SpecificHeatFromDensityInternalEnergy(1, 1), 1.e-15)); + INFO("reference bmod: " << ref_e.BulkModulusFromDensityInternalEnergy(1, 1) + << " test bmod: " + << test_e.BulkModulusFromDensityInternalEnergy(1, 1)); + CHECK(isClose(test_e.BulkModulusFromDensityInternalEnergy(1, 1), + ref_e.BulkModulusFromDensityInternalEnergy(1, 1), 1.e-15)); + INFO("reference Grun. Param.: " + << ref_e.GruneisenParamFromDensityInternalEnergy(1, 1) + << " test Grun. Param.: " << test_e.GruneisenParamFromDensityInternalEnergy(1, 1)); + CHECK(isClose(test_e.GruneisenParamFromDensityInternalEnergy(1, 1), + ref_e.GruneisenParamFromDensityInternalEnergy(1, 1), 1.e-15)); + INFO("reference P: " << ref_e.PressureFromDensityTemperature(1, 1) + << " test P: " << test_e.PressureFromDensityTemperature(1, 1)); + CHECK(isClose(test_e.PressureFromDensityTemperature(1, 1), + ref_e.PressureFromDensityTemperature(1, 1), 1.e-15)); + INFO("reference Cv: " << ref_e.SpecificHeatFromDensityTemperature(1, 1) << " test Cv: " + << test_e.SpecificHeatFromDensityTemperature(1, 1)); + CHECK(isClose(test_e.SpecificHeatFromDensityTemperature(1, 1), + ref_e.SpecificHeatFromDensityTemperature(1, 1), 1.e-15)); + INFO("reference bmod: " << ref_e.BulkModulusFromDensityTemperature(1, 1) + << " test bmod: " + << test_e.BulkModulusFromDensityTemperature(1, 1)); + CHECK(isClose(test_e.BulkModulusFromDensityTemperature(1, 1), + ref_e.BulkModulusFromDensityTemperature(1, 1), 1.e-15)); + INFO("reference Grun. Param.: " << ref_e.GruneisenParamFromDensityTemperature(1, 1) + << " test Grun. Param.: " + << test_e.GruneisenParamFromDensityTemperature(1, 1)); + CHECK(isClose(test_e.GruneisenParamFromDensityTemperature(1, 1), + ref_e.GruneisenParamFromDensityTemperature(1, 1), 1.e-15)); + INFO("reference rho min.: " << ref_e.MinimumDensity() + << " test rho min.: " << test_e.MinimumDensity()); + CHECK(isClose(test_e.MinimumDensity(), ref_e.MinimumDensity(), 1.e-15)); + INFO("reference T min.: " << ref_e.MinimumTemperature() + << " test T min.: " << test_e.MinimumTemperature()); + CHECK(isClose(test_e.MinimumTemperature(), ref_e.MinimumTemperature(), 1.e-15)); return; } diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp new file mode 100644 index 0000000000..117c61e914 --- /dev/null +++ b/test/test_eos_ideal.cpp @@ -0,0 +1,144 @@ +//------------------------------------------------------------------------------ +// © 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 // debug +#include + +#include +#include +#include +#include + +#ifdef SINGULARITY_BUILD_CLOSURE +#include +#endif + +#ifndef CATCH_CONFIG_RUNNER +#include "catch2/catch.hpp" +#endif + +#include + +using singularity::EOS; +using singularity::IdealGas; + +SCENARIO("Ideal gas entropy", "[IdealGas][Entropy]") { + GIVEN("Parameters for an ideal gas with entropy reference states") { + // Create ideal gas EOS ojbect + constexpr Real Cv = 5.0; + constexpr Real gm1 = 0.4; + constexpr Real EntropyT0 = 100; + constexpr Real EntropyRho0 = 1e-03; + EOS host_eos = IdealGas(gm1, Cv, EntropyT0, EntropyRho0); + THEN("The entropy at the reference state should be zero") { + auto entropy = host_eos.EntropyFromDensityTemperature(EntropyRho0, EntropyT0); + INFO("Entropy should be zero but it is " << entropy); + CHECK(isClose(entropy, 0.0, 1.e-14)); + } + GIVEN("A state at the reference temperature and a density whose cube root is the " + "reference density") { + constexpr Real T = EntropyT0; + constexpr Real rho = 0.1; // rho**3 = EntropyRho0 + THEN("The entropy should be 2. / 3. * gm1 * Cv * log(EntropyRho0)") { + const Real entropy_true = 2. / 3. * gm1 * Cv * log(EntropyRho0); + auto entropy = host_eos.EntropyFromDensityTemperature(rho, T); + INFO("Entropy: " << entropy << " True entropy: " << entropy_true); + CHECK(isClose(entropy, entropy_true, 1e-12)); + } + } + GIVEN("A state at the reference density and a temperature whose square is the " + "reference temperature") { + constexpr Real T = 10; // T**2 = EntropyT0 + constexpr Real rho = EntropyRho0; + THEN("The entropy should be -1. / 2. * Cv * log(EntropyT0)") { + const Real entropy_true = -1. / 2. * Cv * log(EntropyT0); + auto entropy = host_eos.EntropyFromDensityTemperature(rho, T); + INFO("Entropy: " << entropy << " True entropy: " << entropy_true); + CHECK(isClose(entropy, entropy_true, 1e-12)); + } + } + } +} + +// A non-standard pattern where we do a reduction +class CheckPofRE { + public: + CheckPofRE(Real *P, Real *rho, Real *sie, int N) : P_(P), rho_(rho), sie_(sie), N_(N) {} + template + void operator()(const T &eos) { + Real *P = P_; + Real *rho = rho_; + Real *sie = sie_; + portableReduce( + "MyCheckPofRE", 0, N_, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += !(isClose(P[i], + eos.PressureFromDensityInternalEnergy(rho[i], sie[i], nullptr), + 1e-15)); + }, + nwrong); + } + + // must be initialized to zero because portableReduce is a simple for loop on host + int nwrong = 0; + + private: + int N_; + Real *P_; + Real *rho_; + Real *sie_; +}; +SCENARIO("Ideal gas vector Evaluate call", "[IdealGas][Evaluate]") { + GIVEN("An ideal gas, and some device memory") { + constexpr Real Cv = 2.0; + constexpr Real gm1 = 0.5; + constexpr int N = 100; + constexpr Real rhomin = 1; + constexpr Real rhomax = 11; + constexpr Real drho = (rhomax - rhomin) / (N - 1); + constexpr Real siemin = 1; + constexpr Real siemax = 11; + constexpr Real dsie = (siemax - siemin) / (N - 1); + + EOS eos_host = IdealGas(gm1, Cv); + auto eos_device = eos_host.GetOnDevice(); + + Real *P = (Real *)PORTABLE_MALLOC(N * sizeof(Real)); + Real *rho = (Real *)PORTABLE_MALLOC(N * sizeof(Real)); + Real *sie = (Real *)PORTABLE_MALLOC(N * sizeof(Real)); + + WHEN("We set P, rho, sie via the scalar API") { + portableFor( + "Set rho and sie", 0, N, PORTABLE_LAMBDA(const int i) { + rho[i] = rhomin + drho * i; // just some diagonal in the 2d rho/sie space + sie[i] = siemin + dsie * i; + P[i] = eos_device.PressureFromDensityInternalEnergy(rho[i], sie[i]); + }); + THEN("The vector Evaluate API can be used to compare") { + CheckPofRE my_op(P, rho, sie, N); + eos_device.Evaluate(my_op); + REQUIRE(my_op.nwrong == 0); + } + } + + eos_device.Finalize(); + PORTABLE_FREE(P); + PORTABLE_FREE(rho); + PORTABLE_FREE(sie); + } +} diff --git a/test/test_eos_modifiers.cpp b/test/test_eos_modifiers.cpp new file mode 100644 index 0000000000..2140001830 --- /dev/null +++ b/test/test_eos_modifiers.cpp @@ -0,0 +1,276 @@ +//------------------------------------------------------------------------------ +// © 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 +#include +#include + +#ifndef CATCH_CONFIG_RUNNER +#include "catch2/catch.hpp" +#endif + +#include + +using singularity::EOS; +using singularity::IdealGas; +using singularity::ScaledEOS; +using singularity::ShiftedEOS; +using singularity::BilinearRampEOS; + +namespace EOSBuilder = singularity::EOSBuilder; +namespace thermalqs = singularity::thermalqs; + +SCENARIO("EOS Builder and Modifiers", "[EOSBuilder][Modifiers][IdealGas]") { + + GIVEN("Parameters for a shifted and scaled ideal gas") { + constexpr Real Cv = 2.0; + constexpr Real gm1 = 0.5; + constexpr Real scale = 2.0; + constexpr Real shift = 0.1; + constexpr Real rho = 2.0; + constexpr Real sie = 0.5; + WHEN("We construct a shifted, scaled IdealGas by hand") { + IdealGas a = IdealGas(gm1, Cv); + ShiftedEOS b = ShiftedEOS(std::move(a), shift); + EOS eos = ScaledEOS>(std::move(b), scale); + THEN("The shift and scale parameters pass through correctly") { + + REQUIRE(eos.PressureFromDensityInternalEnergy(rho, sie) == 0.3); + } + THEN("We can UnmodifyOnce to get the shifted EOS object") { + EOS shifted = eos.UnmodifyOnce(); + REQUIRE(shifted.IsType>()); + AND_THEN("We can extract the unmodified object") { + EOS unmod = eos.GetUnmodifiedObject(); + REQUIRE(unmod.IsType()); + } + } + } + WHEN("We use the EOSBuilder") { + EOSBuilder::EOSType type = EOSBuilder::EOSType::IdealGas; + EOSBuilder::modifiers_t modifiers; + EOSBuilder::params_t base_params, shifted_params, scaled_params; + base_params["Cv"].emplace(Cv); + base_params["gm1"].emplace(gm1); + shifted_params["shift"].emplace(shift); + scaled_params["scale"].emplace(scale); + modifiers[EOSBuilder::EOSModifier::Shifted] = shifted_params; + modifiers[EOSBuilder::EOSModifier::Scaled] = scaled_params; + EOS eos = EOSBuilder::buildEOS(type, base_params, modifiers); + THEN("The shift and scale parameters pass through correctly") { + REQUIRE(eos.PressureFromDensityInternalEnergy(rho, sie) == 0.3); + } + WHEN("We add a ramp") { + EOSBuilder::params_t ramp_params; + Real r0 = 1; + Real a = 1; + Real b = 0; + Real c = 0; + ramp_params["r0"].emplace(r0); + ramp_params["a"].emplace(a); + ramp_params["b"].emplace(b); + ramp_params["c"].emplace(c); + modifiers[EOSBuilder::EOSModifier::BilinearRamp] = ramp_params; + THEN("The EOS is constructed correctly") { + auto eos_ramped = EOSBuilder::buildEOS(type, base_params, modifiers); + } + } + } +#ifdef SINGULARITY_BUILD_CLOSURE + WHEN("We construct a non-modifying modifier") { + EOS ig = IdealGas(gm1, Cv); + EOS igsh = ScaledEOS(IdealGas(gm1, Cv), 1.0); + EOS igsc = ShiftedEOS(IdealGas(gm1, Cv), 0.0); + EOS igra; + // test out the c interface + int enabled[4] = {0, 0, 1, 0}; + Real vals[6] = {0.0, 0.0, 1.e9, 1.0, 2.0, 1.0}; + Real rho0 = 1.e6 / (gm1 * Cv * 293.0); + init_sg_IdealGas(0, &igra, gm1, Cv, enabled, vals); + THEN("The modified EOS should produce equivalent results") { + compare_two_eoss(igsh, ig); + compare_two_eoss(igsc, ig); + compare_two_eoss(igra, ig); + } + } + WHEN("We construct a ramp from a p-alpha model") { + const Real Pe = 5.e7, Pc = 1.e8; + const Real alpha0 = 1.5; + const Real T0 = 293.0; + int enabled[4] = {0, 0, 0, 1}; + Real vals[6] = {0.0, 0.0, alpha0, Pe, Pc, 0.0}; + const Real rho0 = 1.e6 / (gm1 * Cv * T0); + EOS igra; + const Real r0 = rho0 / alpha0; + const Real r1 = Pc / (gm1 * Cv * T0); + const Real rmid = Pe / (gm1 * Cv * T0 * alpha0); + // P(alpha0 * rmid) + const Real P_armid = alpha0 * gm1 * Cv * rmid * T0; + init_sg_IdealGas(0, &igra, gm1, Cv, enabled, vals); + // construct ramp params and evaluate directly for test + const Real a = r0 * Pe / (rmid - r0); + const Real b = r0 * (Pc - Pe) / (r1 - rmid); + const Real c = (Pc * rmid - Pe * r1) / (r0 * (Pc - Pe)); + // density in the middle of the first slope + const Real rho_t1 = 0.5 * (r0 + rmid); + // density in the middle of the second slope + const Real rho_t2 = 0.5 * (rmid + r1); + // P (rho_t1) note that r0 = rho0 / alpha0 + const Real Prhot1 = a * (rho_t1 / r0 - 1.0); + // P (rho_t2) + const Real Prhot2 = b * (rho_t2 / r0 - c); + // bmod (rho_t1) + const Real bmodrt1 = rho_t1 * a / r0; + // bmod (rho_t2) + const Real bmodrt2 = rho_t2 * b / r0; + THEN("P_eos(alpha_0*rmid, T0) = P_ramp(rmid,T0)") { + INFO("P_eos(alpha_0*rmid, T0): " + << P_armid + << " P_ramp(rmid, T0): " << igra.PressureFromDensityTemperature(rmid, T0)); + REQUIRE(isClose(P_armid, igra.PressureFromDensityTemperature(rmid, T0), 1.e-12)); + } + THEN("We obtain correct ramp behavior in P(rho) for rho r1") { + // also check pressures on ramp + INFO("reference P((r0+rmid)/2, T0): " + << Prhot1 << " test P((r0+rmid)/2, T0): " + << igra.PressureFromDensityTemperature(rho_t1, T0)); + REQUIRE(isClose(Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0), 1.e-12)); + INFO("reference P((rmid+r1)/2, T0): " + << Prhot2 << " test P((rmid+r1)/2, T0): " + << igra.PressureFromDensityTemperature(rho_t2, T0)); + REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12)); + // check pressure below and beyond ramp matches unmodified ideal gas + INFO("reference P(0.8*r0, T0): " + << 0.8 * r0 * gm1 * Cv * T0 << " test P(0.8*r0, T0): " + << igra.PressureFromDensityTemperature(0.8 * r0, T0)); + REQUIRE(isClose(0.8 * r0 * gm1 * Cv * T0, + igra.PressureFromDensityTemperature(0.8 * r0, T0), 1.e-12)); + INFO("reference P(1.2*r1, T0): " + << 1.2 * r1 * gm1 * Cv * T0 << " test P(1.2*r1, T0): " + << igra.PressureFromDensityTemperature(1.2 * r1, T0)); + REQUIRE(isClose(1.2 * r1 * gm1 * Cv * T0, + igra.PressureFromDensityTemperature(1.2 * r1, T0), 1.e-12)); + } + THEN("We obtain correct ramp behavior in bmod(rho) for rho r1") { + // check bulk moduli on both pieces of ramp + INFO("reference bmod((r0+rmid)/2, T0): " + << bmodrt1 << " test bmod((r0+rmid)/2, T0): " + << igra.BulkModulusFromDensityTemperature(rho_t1, T0)); + REQUIRE( + isClose(bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0), 1.e-12)); + INFO("reference bmod((rmid+r1)/2, T0): " + << bmodrt2 << " test bmod((rmid+r1)/2, T0): " + << igra.BulkModulusFromDensityTemperature(rho_t2, T0)); + REQUIRE( + isClose(bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0), 1.e-12)); + // check bulk modulus below and beyond ramp matches unmodified ideal gas + INFO("reference bmod(0.8*r0, T0): " + << 0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0 << " test bmod(0.8*r0, T0): " + << igra.BulkModulusFromDensityTemperature(0.8 * r0, T0)); + REQUIRE(isClose(0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0, + igra.BulkModulusFromDensityTemperature(0.8 * r0, T0), 1.e-12)); + INFO("reference bmod(1.2*r1, T0): " + << 1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0 << " test bmod(1.2*r1, T0): " + << igra.BulkModulusFromDensityTemperature(1.2 * r1, T0)); + REQUIRE(isClose(1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0, + igra.BulkModulusFromDensityTemperature(1.2 * r1, T0), 1.e-12)); + } + } +#endif // SINGULARITY_BUILD_CLOSURE + } +} + +SCENARIO("Relativistic EOS", "[EOSBuilder][RelativisticEOS][IdealGas]") { + GIVEN("Parameters for an ideal gas") { + constexpr Real Cv = 2.0; + constexpr Real gm1 = 0.5; + WHEN("We construct a relativistic IdealGas with EOSBuilder") { + EOSBuilder::EOSType type = EOSBuilder::EOSType::IdealGas; + EOSBuilder::modifiers_t modifiers; + EOSBuilder::params_t base_params, relativity_params; + base_params["Cv"].emplace(Cv); + base_params["gm1"].emplace(gm1); + relativity_params["cl"].emplace(1.0); + modifiers[EOSBuilder::EOSModifier::Relativistic] = relativity_params; + EOS eos = EOSBuilder::buildEOS(type, base_params, modifiers); + THEN("The EOS has finite sound speeds") { + constexpr Real rho = 1e3; + constexpr Real sie = 1e3; + Real bmod = eos.BulkModulusFromDensityInternalEnergy(rho, sie); + Real cs2 = bmod / rho; + REQUIRE(cs2 < 1); + } + } + } +} + +SCENARIO("EOS Unit System", "[EOSBuilder][UnitSystem][IdealGas]") { + GIVEN("Parameters for an ideal gas") { + constexpr Real Cv = 2.0; + constexpr Real gm1 = 0.5; + EOSBuilder::EOSType type = EOSBuilder::EOSType::IdealGas; + EOSBuilder::modifiers_t modifiers; + EOSBuilder::params_t base_params, units_params; + base_params["Cv"].emplace(Cv); + base_params["gm1"].emplace(gm1); + GIVEN("Units with a thermal unit system") { + constexpr Real rho_unit = 1e1; + constexpr Real sie_unit = 1e-1; + constexpr Real temp_unit = 123; + WHEN("We construct an IdealGas with EOSBuilder") { + units_params["rho_unit"].emplace(rho_unit); + units_params["sie_unit"].emplace(sie_unit); + units_params["temp_unit"].emplace(temp_unit); + modifiers[EOSBuilder::EOSModifier::UnitSystem] = units_params; + EOS eos = EOSBuilder::buildEOS(type, base_params, modifiers); + THEN("Units cancel out for an ideal gas") { + Real rho = 1e3; + Real sie = 1e3; + Real P = eos.PressureFromDensityInternalEnergy(rho, sie); + Real Ptrue = gm1 * rho * sie; + REQUIRE(std::abs(P - Ptrue) / Ptrue < 1e-3); + } + } + } + GIVEN("Units with length and time units") { + constexpr Real time_unit = 456; + constexpr Real length_unit = 1e2; + constexpr Real mass_unit = 1e6; + constexpr Real temp_unit = 789; + WHEN("We construct an IdealGas with EOSBuilder") { + units_params["use_length_time"].emplace(true); + units_params["time_unit"].emplace(time_unit); + units_params["length_unit"].emplace(length_unit); + units_params["mass_unit"].emplace(mass_unit); + units_params["temp_unit"].emplace(temp_unit); + modifiers[EOSBuilder::EOSModifier::UnitSystem] = units_params; + EOS eos = EOSBuilder::buildEOS(type, base_params, modifiers); + THEN("Units cancel out for an ideal gas") { + Real rho = 1e3; + Real sie = 1e3; + Real P = eos.PressureFromDensityInternalEnergy(rho, sie); + Real Ptrue = gm1 * rho * sie; + REQUIRE(std::abs(P - Ptrue) / Ptrue < 1e-3); + } + } + } + } +} + diff --git a/test/test_eos_stellar_collapse.cpp b/test/test_eos_stellar_collapse.cpp new file mode 100644 index 0000000000..8979554fab --- /dev/null +++ b/test/test_eos_stellar_collapse.cpp @@ -0,0 +1,299 @@ +//------------------------------------------------------------------------------ +// © 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 // debug +#include + +#include +#include +#include +#include +#include + +#ifdef SINGULARITY_BUILD_CLOSURE +#include +#endif + +#ifndef CATCH_CONFIG_RUNNER +#include "catch2/catch.hpp" +#endif + +#include + +#ifdef SPINER_USE_HDF +#ifdef SINGULARITY_TEST_STELLAR_COLLAPSE +SCENARIO("Test 3D reinterpolation to fast log grid", "[StellarCollapse]") { + WHEN("We generate a 3D DataBox of a 2D power law times a line") { + using singularity::StellarCollapse; + constexpr int N2 = 100; + constexpr int N1 = 101; + constexpr int N0 = 102; + StellarCollapse::Grid_t g2(0, 1, N2); + StellarCollapse::Grid_t g1(1, 3, N1); + StellarCollapse::Grid_t g0(2, 4, N0); + StellarCollapse::DataBox db(N2, N1, N0); + + db.setRange(2, g2); + db.setRange(1, g1); + db.setRange(0, g0); + + for (int i2 = 0; i2 < N2; ++i2) { + Real x2 = g2.x(i2); + for (int i1 = 0; i1 < N1; ++i1) { + Real lx1 = g1.x(i1); + for (int i0 = 0; i0 < N0; ++i0) { + Real lx0 = g0.x(i0); + db(i2, i1, i0) = std::log10(x2) + 2 * lx1 + 2 * lx0; + } + } + } + + THEN("The databox should interpolate values correctly") { + const Real x2 = 0.5; + const Real lx1 = 2; + const Real x1 = std::pow(10., lx1); + const Real lx0 = 3; + const Real x0 = std::pow(10., lx0); + REQUIRE(isClose(db.interpToReal(x2, lx1, lx0), std::log10(x2 * x1 * x1 * x0 * x0), + 1e-5)); + } + + THEN("We can re-interpolate to fast log space") { + StellarCollapse::DataBox scratch(N2, N1, N0); + StellarCollapse::dataBoxToFastLogs(db, scratch, true); + + AND_THEN("The fast-log gridded table contains correct ranges") { + REQUIRE(db.range(2) == g2); + REQUIRE(db.range(1).nPoints() == N1); + REQUIRE(isClose(db.range(1).min(), + singularity::FastMath::log10(std::pow(10, g1.min())), 1e-12)); + REQUIRE(isClose(db.range(1).max(), + singularity::FastMath::log10(std::pow(10, g1.max())), 1e-12)); + REQUIRE(db.range(0).nPoints() == N0); + REQUIRE(isClose(db.range(0).min(), + singularity::FastMath::log10(std::pow(10, g0.min())), 1e-12)); + REQUIRE(isClose(db.range(0).max(), + singularity::FastMath::log10(std::pow(10, g0.max())), 1e-12)); + } + + AND_THEN("The re-interpolated fast log is a sane number") {} + + AND_THEN("The fast-log table approximately interpolates the power law") { + const Real x2 = 0.5; + const Real x1 = 100; + const Real x0 = 1000; + const Real lx1 = singularity::FastMath::log10(x1); + const Real lx0 = singularity::FastMath::log10(x0); + const Real lval_interp = db.interpToReal(x2, lx1, lx0); + const Real val_interp = singularity::FastMath::pow10(lval_interp); + const Real val_true = x2 * x1 * x1 * x0 * x0; + const Real rel_diff = + 0.5 * std::abs(val_interp - val_true) / (val_true + val_interp); + REQUIRE(rel_diff <= 1e-3); + } + scratch.finalize(); + } + db.finalize(); + } +} + +SCENARIO("Stellar Collapse EOS", "[StellarCollapse][EOSBuilder]") { + using singularity::IdealGas; + using singularity::StellarCollapse; + const std::string savename = "stellar_collapse_ideal_2.sp5"; + GIVEN("A stellar collapse EOS") { + const std::string filename = "./goldfiles/stellar_collapse_ideal.h5"; + THEN("We can load the file") { // don't bother filtering bmod here. + StellarCollapse sc(filename, false, false); + AND_THEN("Some properties we expect for ideal gas hold") { + Real lambda[2]; + Real rho, t, sie, p, cv, b, dpde, dvdt; + sc.ValuesAtReferenceState(rho, t, sie, p, cv, b, dpde, dvdt, lambda); + Real yemin = sc.YeMin(); + Real yemax = sc.YeMax(); + int N = 123; + Real dY = (yemax - yemin) / (N + 1); + for (int i = 0; i < N; ++i) { + Real Ye = yemin + i * dY; + lambda[0] = Ye; + REQUIRE(isClose(sie, sc.InternalEnergyFromDensityTemperature(rho, t, lambda))); + Real Xa, Xh, Xn, Xp, Abar, Zbar; + sc.MassFractionsFromDensityTemperature(rho, t, Xa, Xh, Xn, Xp, Abar, Zbar, + lambda); + REQUIRE(isClose(Ye, Xp)); + } + Real rhomin = sc.rhoMin(); + Real rhomax = sc.rhoMax(); + Real drho = (rhomax - rhomin) / (N + 1); + for (int i = 0; i < N; ++i) { + Real rho = rhomin + i * drho; + REQUIRE(isClose(sie, sc.InternalEnergyFromDensityTemperature(rho, t, lambda))); + } + } + GIVEN("An Ideal Gas equation of state") { + constexpr Real gamma = 1.4; + constexpr Real mp = 1.67262171e-24; + constexpr Real kb = 1.3806505e-16; + constexpr Real Cv = kb / (mp * (gamma - 1)); // mean molecular weight = mp + IdealGas ig(gamma - 1, Cv); + auto ig_d = ig.GetOnDevice(); + THEN("The tabulated gamma Stellar Collapse and the gamma agree roughly") { + Real yemin = sc.YeMin(); + Real yemax = sc.YeMax(); + Real tmin = sc.TMin(); + Real tmax = sc.TMax(); + Real ltmin = std::log10(tmin); + Real ltmax = std::log10(tmax); + Real lrhomin = std::log10(sc.rhoMin()); + Real lrhomax = std::log10(sc.rhoMax()); + auto sc_d = sc.GetOnDevice(); + + int nwrong_h = 0; +#ifdef PORTABILITY_STRATEGY_KOKKOS + using atomic_view = Kokkos::MemoryTraits; + Kokkos::View nwrong_d("wrong"); +#else + PortableMDArray nwrong_d(&nwrong_h, 1); +#endif + + const int N = 123; + const Real dY = (yemax - yemin) / (N + 1); + const Real dlT = (ltmax - ltmin) / (N + 1); + const Real dlR = (lrhomax - lrhomin) / (N + 1); + portableFor( + "fill eos", 0, N, 0, N, 0, N, + PORTABLE_LAMBDA(const int &k, const int &j, const int &i) { + Real lambda[2]; + Real Ye = yemin + k * dY; + Real lT = ltmin + j * dlT; + Real lR = lrhomin + i * dlR; + Real T = std::pow(10., lT); + Real R = std::pow(10., lR); + Real e1, e2, p1, p2, cv1, cv2, b1, b2; + unsigned long output = (singularity::thermalqs::pressure | + singularity::thermalqs::specific_internal_energy | + singularity::thermalqs::specific_heat | + singularity::thermalqs::bulk_modulus); + lambda[0] = Ye; + + sc_d.FillEos(R, T, e1, p1, cv1, b1, output, lambda); + ig_d.FillEos(R, T, e2, p2, cv2, b2, output, lambda); + if (!isClose(e1, e2)) { + nwrong_d() += 1; + } + if (!isClose(p1, p2)) { + nwrong_d() += 1; + } + if (!isClose(cv1, cv2)) { + nwrong_d() += 1; + } + if (!isClose(b1, b2)) { + nwrong_d() += 1; + } + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(nwrong_h, nwrong_d); +#endif + REQUIRE(nwrong_h == 0); + sc_d.Finalize(); + } + ig_d.Finalize(); + ig.Finalize(); + } + AND_THEN("We can save the file to SP5") { + sc.Save(savename); + AND_THEN("We can load the sp5 file") { + StellarCollapse sc2(savename, true); + AND_THEN("The two stellar collapse EOS's agree") { + + Real yemin = sc.YeMin(); + Real yemax = sc.YeMax(); + Real tmin = sc.TMin(); + Real tmax = sc.TMax(); + Real ltmin = std::log10(tmin); + Real ltmax = std::log10(tmax); + Real lrhomin = std::log10(sc.rhoMin()); + Real lrhomax = std::log10(sc.rhoMax()); + REQUIRE(yemin == sc2.YeMin()); + REQUIRE(yemax == sc2.YeMax()); + REQUIRE(sc.TMin() == sc2.TMin()); + REQUIRE(sc.TMax() == sc2.TMax()); + REQUIRE(isClose(lrhomin, std::log10(sc2.rhoMin()), 1e-12)); + REQUIRE(isClose(lrhomax, std::log10(sc2.rhoMax()), 1e-12)); + + auto sc1_d = sc.GetOnDevice(); + auto sc2_d = sc2.GetOnDevice(); + + int nwrong_h = 0; +#ifdef PORTABILITY_STRATEGY_KOKKOS + using atomic_view = Kokkos::MemoryTraits; + Kokkos::View nwrong_d("wrong"); +#else + PortableMDArray nwrong_d(&nwrong_h, 1); +#endif + + const int N = 123; + constexpr Real gamma = 1.4; + const Real dY = (yemax - yemin) / (N + 1); + const Real dlT = (ltmax - ltmin) / (N + 1); + const Real dlR = (lrhomax - lrhomin) / (N + 1); + portableFor( + "fill eos", 0, N, 0, N, 0, N, + PORTABLE_LAMBDA(const int &k, const int &j, const int &i) { + Real lambda[2]; + Real Ye = yemin + k * dY; + Real lT = ltmin + j * dlT; + Real lR = lrhomin + i * dlR; + Real T = std::pow(10., lT); + Real R = std::pow(10., lR); + Real e1, e2, p1, p2, cv1, cv2, b1, b2, s1, s2; + unsigned long output = + (singularity::thermalqs::pressure | + singularity::thermalqs::specific_internal_energy | + singularity::thermalqs::specific_heat | + singularity::thermalqs::bulk_modulus); + lambda[0] = Ye; + + sc1_d.FillEos(R, T, e1, p1, cv1, b1, output, lambda); + sc2_d.FillEos(R, T, e2, p2, cv2, b2, output, lambda); + // Fill entropy. Will need to change later. + s1 = sc1_d.EntropyFromDensityTemperature(R, T, lambda); + s2 = p2 * std::pow(R, -gamma); // ideal + if (!isClose(e1, e2)) nwrong_d() += 1; + if (!isClose(p1, p2)) nwrong_d() += 1; + if (!isClose(cv1, cv2)) nwrong_d() += 1; + if (!isClose(b1, b2)) nwrong_d() += 1; + if (!isClose(s1, s2)) nwrong_d() += 1; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(nwrong_h, nwrong_d); +#endif + REQUIRE(nwrong_h == 0); + + sc1_d.Finalize(); + sc2_d.Finalize(); + } + sc2.Finalize(); + } + } + sc.Finalize(); + } + } +} +#endif // SINGULARITY_TEST_STELLAR_COLLAPSE +#endif // USE_HDF5 diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp new file mode 100644 index 0000000000..bce04d0d95 --- /dev/null +++ b/test/test_eos_tabulated.cpp @@ -0,0 +1,314 @@ +//------------------------------------------------------------------------------ +// © 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 // debug +#include + +#include +#include +#include +#include +#include + +#ifdef SINGULARITY_BUILD_CLOSURE +#include +#endif + +#ifndef CATCH_CONFIG_RUNNER +#include "catch2/catch.hpp" +#endif + +#include + +using singularity::EOS; + +#ifdef SPINER_USE_HDF +using singularity::SpinerEOSDependsRhoSie; +using singularity::SpinerEOSDependsRhoT; +#endif + +#ifdef SINGULARITY_USE_EOSPAC +using singularity::EOSPAC; +#endif + +namespace EOSBuilder = singularity::EOSBuilder; +namespace thermalqs = singularity::thermalqs; + +const std::string eosName = "../materials.sp5"; +const std::string airName = "air"; +const std::string steelName = "stainless steel 347"; + +#ifdef SPINER_USE_HDF +#ifdef SINGULARITY_TEST_SESAME +constexpr int steelID = 4272; +constexpr int airID = 5030; +constexpr int DTID = 5267; +constexpr int gID = 2700; +constexpr Real ev2k = 1.160451812e4; +#endif // SINGULARITY_TEST_SESAME +#endif // SPINER_USE_HDF + +#ifdef SPINER_USE_HDF +#ifdef SINGULARITY_TEST_SESAME +#ifdef SINGULARITY_USE_EOSPAC +SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS],[DependsRhoT][EOSPAC]") { + + GIVEN("SpinerEOS and EOSPAC EOS for steel can be initialized with matid") { + EOS steelEOS_host_polymorphic = SpinerEOSDependsRhoT(eosName, steelID); + EOS steelEOS = steelEOS_host_polymorphic.GetOnDevice(); + auto steelEOS_host = steelEOS_host_polymorphic.get(); + + EOS eospac = EOSPAC(steelID); + + THEN("The correct metadata is read in") { + REQUIRE(steelEOS_host.matid() == steelID); + + AND_THEN("We can get a reference density and temperature") { + Real rho, T, sie, P, cv, bmod, dpde, dvdt; + Real rho_pac, T_pac, sie_pac, P_pac, cv_pac, bmod_pac, dpde_pac, dvdt_pac; + steelEOS_host.ValuesAtReferenceState(rho, T, sie, P, cv, bmod, dpde, dvdt); + eospac.ValuesAtReferenceState(rho_pac, T_pac, sie_pac, P_pac, cv_pac, bmod_pac, + dpde_pac, dvdt_pac); + REQUIRE(isClose(rho, rho_pac)); + REQUIRE(isClose(T, T_pac)); + } + + // TODO: this needs to be a much more rigorous test + AND_THEN("Quantities can be read from density and temperature") { + const Real sie_pac = eospac.InternalEnergyFromDensityTemperature(1e0, 1e6); + + int nw_ie{0}; +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + using atomic_view = Kokkos::MemoryTraits; + Kokkos::View n_wrong_ie("wrong_ie"); +#else + PortableMDArray n_wrong_ie(&nw_ie, 1); +#endif + portableFor( + "calc ie's steel", 0, 100, PORTABLE_LAMBDA(const int &i) { + const double ie = steelEOS.InternalEnergyFromDensityTemperature(1e0, 1e6); + if (!isClose(ie, sie_pac)) n_wrong_ie() += 1; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(nw_ie, n_wrong_ie); +#endif + REQUIRE(nw_ie == 0); + } + + AND_THEN("rho(P,T) correct for P=1atm, T=freezing") { + Real T = 273; // Kelvin + Real P = 1e6; // barye + Real rho, sie; // output + Real rho_pac, sie_pac; + std::vector lambda(steelEOS_host_polymorphic.nlambda()); + steelEOS_host_polymorphic.DensityEnergyFromPressureTemperature( + P, T, lambda.data(), rho, sie); + eospac.DensityEnergyFromPressureTemperature(P, T, nullptr, rho_pac, sie_pac); + REQUIRE(isClose(rho, rho_pac)); + } + } + // Failing to call finalize leads to a memory leak, + // but otherwise behaviour is as expected. + steelEOS_host_polymorphic.Finalize(); // host and device must be + // finalized separately. + steelEOS.Finalize(); // cleans up memory on device. + } + + GIVEN("SpinerEOS and EOSPAC for air can be initialized with matid") { + SpinerEOSDependsRhoT airEOS_host = SpinerEOSDependsRhoT(eosName, airID); + EOS airEOS = airEOS_host.GetOnDevice(); + EOS eospac = EOSPAC(airID); + constexpr Real rho = 1e0; + constexpr Real sie = 2.5e-4; + THEN("We can get a reference state") { + Real rho, T, sie, P, cv, bmod, dpde, dvdt; + Real rho_pac, T_pac, sie_pac, P_pac, cv_pac, bmod_pac, dpde_pac, dvdt_pac; + airEOS_host.ValuesAtReferenceState(rho, T, sie, P, cv, bmod, dpde, dvdt); + eospac.ValuesAtReferenceState(rho_pac, T_pac, sie_pac, P_pac, cv_pac, bmod_pac, + dpde_pac, dvdt_pac); + REQUIRE(isClose(rho, rho_pac)); + REQUIRE(isClose(T, T_pac)); + } + THEN("Quantities of rho and sie look sane on both host and device") { + const Real gm1_host = airEOS_host.GruneisenParamFromDensityInternalEnergy(rho, sie); + const Real T_host = airEOS_host.TemperatureFromDensityInternalEnergy(rho, sie); + const Real cv_host = airEOS_host.SpecificHeatFromDensityInternalEnergy(rho, sie); + + int nw_gm1{0}, nw_cv{0}; +#ifdef PORTABILITY_STRATEGY_KOKKOS + using atomic_view = Kokkos::MemoryTraits; + Kokkos::fence(); + Kokkos::View n_wrong_gm1("wrong_gm1"); + Kokkos::View n_wrong_cv("wrong_cv"); +#else + PortableMDArray n_wrong_gm1(&nw_gm1, 1); + PortableMDArray n_wrong_cv(&nw_cv, 1); +#endif + portableFor( + "calc gm1 and cv", 0, 100, PORTABLE_LAMBDA(const int &i) { + const Real gm1 = airEOS.GruneisenParamFromDensityInternalEnergy(rho, sie); + const Real cv = airEOS.SpecificHeatFromDensityInternalEnergy(rho, sie); + if (!isClose(gm1, gm1_host)) n_wrong_gm1() += 1; + if (!isClose(cv, cv_host)) n_wrong_cv() += 1; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(nw_gm1, n_wrong_gm1); + Kokkos::deep_copy(nw_cv, n_wrong_cv); +#endif + REQUIRE(nw_gm1 == 0); + REQUIRE(nw_cv == 0); + } + AND_THEN("P(rho, sie) correct for extrapolation regime") { + Real rho = 1; + Real sie = 2.43e16; + Real P_pac = eospac.PressureFromDensityInternalEnergy(rho, sie); + Real P_spi = airEOS_host.PressureFromDensityInternalEnergy(rho, sie); + REQUIRE(isClose(P_pac, P_spi)); + } + airEOS_host.Finalize(); + airEOS.Finalize(); + } + + GIVEN("EOS initialized with matid") { + SpinerEOSDependsRhoT eos_spiner = SpinerEOSDependsRhoT(eosName, DTID); + EOS eos_eospac = EOSPAC(DTID); + Real P = 1e8; + Real rho = 1.28e-3; + THEN("Inversion for T(rho,P) works on host") { + Real T, sie, cv, bmod; + Real T_pac, sie_pac, cv_pac, bmod_pac; + const unsigned long output = + (thermalqs::temperature | thermalqs::specific_internal_energy | + thermalqs::specific_heat | thermalqs::bulk_modulus); + eos_spiner.FillEos(rho, T, sie, P, cv, bmod, output); + eos_eospac.FillEos(rho, T_pac, sie_pac, P, cv_pac, bmod_pac, output); + REQUIRE(isClose(T, T_pac)); + REQUIRE(isClose(sie, sie_pac)); + REQUIRE(isClose(cv, cv_pac)); + } + eos_spiner.Finalize(); + } + + GIVEN("EOS initialized with matid") { + SpinerEOSDependsRhoT eos_spiner = SpinerEOSDependsRhoT(eosName, gID); + EOS eos_eospac = EOSPAC(gID); + THEN("PT lookup works on the host") { + Real P = 1e6; // cgs + Real T = 0.025 / ev2k; // K + Real rho, sie; + Real rho_pac, sie_pac; + std::vector lambda(eos_spiner.nlambda()); + eos_spiner.DensityEnergyFromPressureTemperature(P, T, lambda.data(), rho, sie); + eos_eospac.DensityEnergyFromPressureTemperature(P, T, lambda.data(), rho_pac, + sie_pac); + REQUIRE(isClose(rho, rho_pac)); + } + eos_spiner.Finalize(); + } +} + +// Disabling these tests for now as the DependsRhoSie code is not well-maintained +SCENARIO("SpinerEOS depends on rho and sie", "[SpinerEOS],[DependsRhoSie]") { + + GIVEN("SpinerEOSes for steel can be initialised with matid") { + SpinerEOSDependsRhoSie steelEOS_host(eosName, steelID); + EOS steelEOS = steelEOS_host.GetOnDevice(); + THEN("The correct metadata is read in") { + REQUIRE(steelEOS_host.matid() == steelID); + + int nw_ie2{0}, nw_te2{0}; +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + using atomic_view = Kokkos::MemoryTraits; + Kokkos::View n_wrong_ie("wrong_ie"); + Kokkos::View n_wrong_te("wrong_te"); +#else + PortableMDArray n_wrong_ie(&nw_ie2, 1); + PortableMDArray n_wrong_te(&nw_te2, 1); +#endif + portableFor( + "calc ie's steel 2", 0, 100, PORTABLE_LAMBDA(const int &i) { + const Real ie{steelEOS.InternalEnergyFromDensityTemperature(1e0, 1e6)}; + const Real te{steelEOS.TemperatureFromDensityInternalEnergy(1e0, 1e12)}; + if (!isClose(ie, 4.96415e13)) n_wrong_ie() += 1; + if (!isClose(te, 75594.6)) n_wrong_te() += 1; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(nw_ie2, n_wrong_ie); + Kokkos::deep_copy(nw_te2, n_wrong_te); +#endif + REQUIRE(nw_ie2 == 0); + REQUIRE(nw_te2 == 0); + } + // this can be removed with with reference counting or other tricks + steelEOS_host.Finalize(); // cleans up host memory + steelEOS.Finalize(); // cleans up device memory + } +} + +SCENARIO("EOS Builder and SpinerEOS", + "[SpinerEOS],[EOSBuilder],[GetOnDevice],[Finalize]") { + GIVEN("Parameters for shift and scale") { + constexpr Real shift = 0.0; + constexpr Real scale = 1.0; + WHEN("We construct a SpinerEOS with EOSBuilder") { + EOSBuilder::EOSType type = EOSBuilder::EOSType::SpinerEOSDependsRhoT; + EOSBuilder::modifiers_t modifiers; + EOSBuilder::params_t base_params, shifted_params, scaled_params; + base_params["filename"].emplace(eosName); + base_params["matid"].emplace(steelID); + shifted_params["shift"].emplace(shift); + scaled_params["scale"].emplace(scale); + modifiers[EOSBuilder::EOSModifier::Shifted] = shifted_params; + modifiers[EOSBuilder::EOSModifier::Scaled] = scaled_params; + EOS eosHost = EOSBuilder::buildEOS(type, base_params, modifiers); + THEN("The EOS is consistent.") { + REQUIRE( + isClose(4.96416e13, eosHost.InternalEnergyFromDensityTemperature(1e0, 1e6))); + } + WHEN("We get an EOS on device") { + EOS eosDevice = eosHost.GetOnDevice(); + THEN("EOS calls match raw access") { + int nw_bm{0}; +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + using atomic_view = Kokkos::MemoryTraits; + Kokkos::View n_wrong_bm("wrong_bm"); +#else + PortableMDArray n_wrong_bm(&nw_bm, 1); +#endif + portableFor( + "calc ie's steel 3", 0, 100, PORTABLE_LAMBDA(const int &i) { + const Real bm{eosDevice.BulkModulusFromDensityTemperature(1e0, 1e6)}; + if (!isClose(bm, 2.55268e13)) n_wrong_bm() += 1; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(nw_bm, n_wrong_bm); +#endif + REQUIRE(nw_bm == 0); + } + eosDevice.Finalize(); + } + eosHost.Finalize(); + } + } +} +#endif // SINGULARITY_USE_EOSPAC +#endif // SINGULARITY_TEST_SESAME +#endif // SPINER_USE_HDF diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp deleted file mode 100644 index 04f5ff009c..0000000000 --- a/test/test_eos_unit.cpp +++ /dev/null @@ -1,1441 +0,0 @@ -//------------------------------------------------------------------------------ -// © 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 // debug -#include - -#include -#include -#include -#include -#include -#include -#include - -#ifdef SINGULARITY_BUILD_CLOSURE -#include -#endif - -#ifndef CATCH_CONFIG_RUNNER -#include "catch2/catch.hpp" -#endif - -#include - -using singularity::BilinearRampEOS; -using singularity::EOS; -using singularity::Gruneisen; -using singularity::IdealGas; -using singularity::ScaledEOS; -using singularity::ShiftedEOS; - -#ifdef SPINER_USE_HDF -using singularity::SpinerEOSDependsRhoSie; -using singularity::SpinerEOSDependsRhoT; -#endif - -#ifdef SINGULARITY_USE_EOSPAC -using singularity::EOSPAC; -#endif - -namespace EOSBuilder = singularity::EOSBuilder; -namespace thermalqs = singularity::thermalqs; - -const std::string eosName = "../materials.sp5"; -const std::string airName = "air"; -const std::string steelName = "stainless steel 347"; - -#ifdef SPINER_USE_HDF -#ifdef SINGULARITY_TEST_SESAME -constexpr int steelID = 4272; -constexpr int airID = 5030; -constexpr int DTID = 5267; -constexpr int gID = 2700; -constexpr Real ev2k = 1.160451812e4; -#endif // SINGULARITY_TEST_SESAME -#endif // SPINER_USE_HDF - -template -inline void compare_two_eoss(E1 &&test_e, E2 &&ref_e) { - // compare all individual member functions with 1 as inputs, - // this function is meant to catch mis-implementations of - // modifiers that can be initialized in such a way as to - // be equivalent of an unmodified eos. Best used with analytic - // eoss. - INFO("reference T: " << ref_e.TemperatureFromDensityInternalEnergy(1, 1) << " test T: " - << test_e.TemperatureFromDensityInternalEnergy(1, 1)); - CHECK(isClose(test_e.TemperatureFromDensityInternalEnergy(1, 1), - ref_e.TemperatureFromDensityInternalEnergy(1, 1), 1.e-15)); - INFO("reference sie: " << ref_e.InternalEnergyFromDensityTemperature(1, 1) - << " test sie: " - << test_e.InternalEnergyFromDensityTemperature(1, 1)); - CHECK(isClose(test_e.InternalEnergyFromDensityTemperature(1, 1), - ref_e.InternalEnergyFromDensityTemperature(1, 1), 1.e-15)); - INFO("reference P: " << ref_e.PressureFromDensityInternalEnergy(1, 1) - << " test P: " << test_e.PressureFromDensityInternalEnergy(1, 1)); - CHECK(isClose(test_e.PressureFromDensityInternalEnergy(1, 1), - ref_e.PressureFromDensityInternalEnergy(1, 1), 1.e-15)); - INFO("reference Cv: " << ref_e.SpecificHeatFromDensityInternalEnergy(1, 1) - << " test Cv: " - << test_e.SpecificHeatFromDensityInternalEnergy(1, 1)); - CHECK(isClose(test_e.SpecificHeatFromDensityInternalEnergy(1, 1), - ref_e.SpecificHeatFromDensityInternalEnergy(1, 1), 1.e-15)); - INFO("reference bmod: " << ref_e.BulkModulusFromDensityInternalEnergy(1, 1) - << " test bmod: " - << test_e.BulkModulusFromDensityInternalEnergy(1, 1)); - CHECK(isClose(test_e.BulkModulusFromDensityInternalEnergy(1, 1), - ref_e.BulkModulusFromDensityInternalEnergy(1, 1), 1.e-15)); - INFO("reference Grun. Param.: " - << ref_e.GruneisenParamFromDensityInternalEnergy(1, 1) - << " test Grun. Param.: " << test_e.GruneisenParamFromDensityInternalEnergy(1, 1)); - CHECK(isClose(test_e.GruneisenParamFromDensityInternalEnergy(1, 1), - ref_e.GruneisenParamFromDensityInternalEnergy(1, 1), 1.e-15)); - INFO("reference P: " << ref_e.PressureFromDensityTemperature(1, 1) - << " test P: " << test_e.PressureFromDensityTemperature(1, 1)); - CHECK(isClose(test_e.PressureFromDensityTemperature(1, 1), - ref_e.PressureFromDensityTemperature(1, 1), 1.e-15)); - INFO("reference Cv: " << ref_e.SpecificHeatFromDensityTemperature(1, 1) << " test Cv: " - << test_e.SpecificHeatFromDensityTemperature(1, 1)); - CHECK(isClose(test_e.SpecificHeatFromDensityTemperature(1, 1), - ref_e.SpecificHeatFromDensityTemperature(1, 1), 1.e-15)); - INFO("reference bmod: " << ref_e.BulkModulusFromDensityTemperature(1, 1) - << " test bmod: " - << test_e.BulkModulusFromDensityTemperature(1, 1)); - CHECK(isClose(test_e.BulkModulusFromDensityTemperature(1, 1), - ref_e.BulkModulusFromDensityTemperature(1, 1), 1.e-15)); - INFO("reference Grun. Param.: " << ref_e.GruneisenParamFromDensityTemperature(1, 1) - << " test Grun. Param.: " - << test_e.GruneisenParamFromDensityTemperature(1, 1)); - CHECK(isClose(test_e.GruneisenParamFromDensityTemperature(1, 1), - ref_e.GruneisenParamFromDensityTemperature(1, 1), 1.e-15)); - INFO("reference rho min.: " << ref_e.MinimumDensity() - << " test rho min.: " << test_e.MinimumDensity()); - CHECK(isClose(test_e.MinimumDensity(), ref_e.MinimumDensity(), 1.e-15)); - INFO("reference T min.: " << ref_e.MinimumTemperature() - << " test T min.: " << test_e.MinimumTemperature()); - CHECK(isClose(test_e.MinimumTemperature(), ref_e.MinimumTemperature(), 1.e-15)); - return; -} - -SCENARIO("Test that we can either throw an error on host or do nothing on device", - "[RequireMaybe]") { - // TODO(JMM): For whatever reason, the preprocessor does not like it if I - // call `PORTABLE_ALWAYS_THROW_OR_ABORT - REQUIRE_MAYBE_THROWS(PORTABLE_ALWAYS_THROW_OR_ABORT("Error message")); -} - -SCENARIO("Test that fast logs are invertible and run on device", "[FastMath]") { - GIVEN("A set of values to invert over a large dynamic range") { - constexpr Real LXMIN = -20; - constexpr Real LXMAX = 32; - constexpr int NX = 1000; - constexpr Real DLX = (LXMAX - LXMIN) / (NX - 1); - Real *x = (Real *)PORTABLE_MALLOC(NX * sizeof(Real)); - portableFor( - "Set x values", 0, NX, PORTABLE_LAMBDA(const int i) { - const Real lx = LXMIN + i * DLX; - x[i] = std::pow(10., lx); - }); - THEN("The fast exp of the fast log returns the original") { - int nw_ie = 0; -#ifdef PORTABILITY_STRATEGY_KOKKOS - using atomic_view = Kokkos::MemoryTraits; - Kokkos::View n_wrong_ie("wrong_ie"); -#else - PortableMDArray n_wrong_ie(&nw_ie, 1); -#endif - portableFor( - "try out the fast math", 0, NX, PORTABLE_LAMBDA(const int i) { - constexpr Real machine_eps = std::numeric_limits::epsilon(); - constexpr Real acceptable_err = 100 * machine_eps; - const Real lx = singularity::FastMath::log10(x[i]); - const Real elx = singularity::FastMath::pow10(lx); - const Real rel_err = 2.0 * std::abs(x[i] - elx) / - (std::abs(x[i]) + std::abs(elx) + machine_eps); - n_wrong_ie() += (rel_err > acceptable_err); - }); -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(nw_ie, n_wrong_ie); -#endif - REQUIRE(nw_ie == 0); - } - PORTABLE_FREE(x); - } -} - -SCENARIO("Rudimentary test of the root finder", "[RootFinding1D]") { - - GIVEN("Root finding") { - using namespace RootFinding1D; - - THEN("A root can be found for shift = 1, scale = 2, offset = 0.5") { - int ntimes = 100; - Real guess = 0; - Real root; - Status status; - Real shift = 1; - Real scale = 2; - Real offset = 0.5; - - auto f = PORTABLE_LAMBDA(const Real x) { return myAtan(x, shift, scale, offset); }; - Status *statusesp = (Status *)PORTABLE_MALLOC(ntimes * sizeof(Status)); - Real *rootsp = (Real *)PORTABLE_MALLOC(ntimes * sizeof(Real)); -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::View> statuses(statusesp, - ntimes); - Kokkos::View> roots(rootsp, ntimes); -#else - PortableMDArray statuses; - PortableMDArray roots; - statuses.NewPortableMDArray(statusesp, ntimes); - roots.NewPortableMDArray(rootsp, ntimes); -#endif - portableFor( - "find roots", 0, ntimes, PORTABLE_LAMBDA(const int i) { - statuses(i) = regula_falsi(f, 0, guess, -1, 3, 1e-10, 1e-10, roots(i)); - }); -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::View s_copy(statuses, 0); - Kokkos::View r_copy(roots, 0); - Kokkos::deep_copy(root, r_copy); - Kokkos::deep_copy(status, s_copy); -#else - status = statuses(ntimes - 1); - root = roots(ntimes - 1); -#endif - - PORTABLE_FREE(statusesp); - PORTABLE_FREE(rootsp); - REQUIRE(status == Status::SUCCESS); - REQUIRE(isClose(root, 0.744658)); - } - } -} - -SCENARIO("EOS Variant Type", "[Variant][EOS]") { - // print out the eos type - std::cout << demangle(typeid(EOS).name()) << std::endl; -} - -SCENARIO("EOS Builder and Modifiers", "[EOSBuilder][Modifiers][IdealGas]") { - - GIVEN("Parameters for a shifted and scaled ideal gas") { - constexpr Real Cv = 2.0; - constexpr Real gm1 = 0.5; - constexpr Real scale = 2.0; - constexpr Real shift = 0.1; - constexpr Real rho = 2.0; - constexpr Real sie = 0.5; - WHEN("We construct a shifted, scaled IdealGas by hand") { - IdealGas a = IdealGas(gm1, Cv); - ShiftedEOS b = ShiftedEOS(std::move(a), shift); - EOS eos = ScaledEOS>(std::move(b), scale); - THEN("The shift and scale parameters pass through correctly") { - - REQUIRE(eos.PressureFromDensityInternalEnergy(rho, sie) == 0.3); - } - THEN("We can UnmodifyOnce to get the shifted EOS object") { - EOS shifted = eos.UnmodifyOnce(); - REQUIRE(shifted.IsType>()); - AND_THEN("We can extract the unmodified object") { - EOS unmod = eos.GetUnmodifiedObject(); - REQUIRE(unmod.IsType()); - } - } - } - WHEN("We use the EOSBuilder") { - EOSBuilder::EOSType type = EOSBuilder::EOSType::IdealGas; - EOSBuilder::modifiers_t modifiers; - EOSBuilder::params_t base_params, shifted_params, scaled_params; - base_params["Cv"].emplace(Cv); - base_params["gm1"].emplace(gm1); - shifted_params["shift"].emplace(shift); - scaled_params["scale"].emplace(scale); - modifiers[EOSBuilder::EOSModifier::Shifted] = shifted_params; - modifiers[EOSBuilder::EOSModifier::Scaled] = scaled_params; - EOS eos = EOSBuilder::buildEOS(type, base_params, modifiers); - THEN("The shift and scale parameters pass through correctly") { - REQUIRE(eos.PressureFromDensityInternalEnergy(rho, sie) == 0.3); - } - WHEN("We add a ramp") { - EOSBuilder::params_t ramp_params; - Real r0 = 1; - Real a = 1; - Real b = 0; - Real c = 0; - ramp_params["r0"].emplace(r0); - ramp_params["a"].emplace(a); - ramp_params["b"].emplace(b); - ramp_params["c"].emplace(c); - modifiers[EOSBuilder::EOSModifier::BilinearRamp] = ramp_params; - THEN("The EOS is constructed correctly") { - auto eos_ramped = EOSBuilder::buildEOS(type, base_params, modifiers); - } - } - } -#ifdef SINGULARITY_BUILD_CLOSURE - WHEN("We construct a non-modifying modifier") { - EOS ig = IdealGas(gm1, Cv); - EOS igsh = ScaledEOS(IdealGas(gm1, Cv), 1.0); - EOS igsc = ShiftedEOS(IdealGas(gm1, Cv), 0.0); - EOS igra; - // test out the c interface - int enabled[4] = {0, 0, 1, 0}; - Real vals[6] = {0.0, 0.0, 1.e9, 1.0, 2.0, 1.0}; - Real rho0 = 1.e6 / (gm1 * Cv * 293.0); - init_sg_IdealGas(0, &igra, gm1, Cv, enabled, vals); - THEN("The modified EOS should produce equivalent results") { - compare_two_eoss(igsh, ig); - compare_two_eoss(igsc, ig); - compare_two_eoss(igra, ig); - } - } - WHEN("We construct a ramp from a p-alpha model") { - const Real Pe = 5.e7, Pc = 1.e8; - const Real alpha0 = 1.5; - const Real T0 = 293.0; - int enabled[4] = {0, 0, 0, 1}; - Real vals[6] = {0.0, 0.0, alpha0, Pe, Pc, 0.0}; - const Real rho0 = 1.e6 / (gm1 * Cv * T0); - EOS igra; - const Real r0 = rho0 / alpha0; - const Real r1 = Pc / (gm1 * Cv * T0); - const Real rmid = Pe / (gm1 * Cv * T0 * alpha0); - // P(alpha0 * rmid) - const Real P_armid = alpha0 * gm1 * Cv * rmid * T0; - init_sg_IdealGas(0, &igra, gm1, Cv, enabled, vals); - // construct ramp params and evaluate directly for test - const Real a = r0 * Pe / (rmid - r0); - const Real b = r0 * (Pc - Pe) / (r1 - rmid); - const Real c = (Pc * rmid - Pe * r1) / (r0 * (Pc - Pe)); - // density in the middle of the first slope - const Real rho_t1 = 0.5 * (r0 + rmid); - // density in the middle of the second slope - const Real rho_t2 = 0.5 * (rmid + r1); - // P (rho_t1) note that r0 = rho0 / alpha0 - const Real Prhot1 = a * (rho_t1 / r0 - 1.0); - // P (rho_t2) - const Real Prhot2 = b * (rho_t2 / r0 - c); - // bmod (rho_t1) - const Real bmodrt1 = rho_t1 * a / r0; - // bmod (rho_t2) - const Real bmodrt2 = rho_t2 * b / r0; - THEN("P_eos(alpha_0*rmid, T0) = P_ramp(rmid,T0)") { - INFO("P_eos(alpha_0*rmid, T0): " - << P_armid - << " P_ramp(rmid, T0): " << igra.PressureFromDensityTemperature(rmid, T0)); - REQUIRE(isClose(P_armid, igra.PressureFromDensityTemperature(rmid, T0), 1.e-12)); - } - THEN("We obtain correct ramp behavior in P(rho) for rho r1") { - // also check pressures on ramp - INFO("reference P((r0+rmid)/2, T0): " - << Prhot1 << " test P((r0+rmid)/2, T0): " - << igra.PressureFromDensityTemperature(rho_t1, T0)); - REQUIRE(isClose(Prhot1, igra.PressureFromDensityTemperature(rho_t1, T0), 1.e-12)); - INFO("reference P((rmid+r1)/2, T0): " - << Prhot2 << " test P((rmid+r1)/2, T0): " - << igra.PressureFromDensityTemperature(rho_t2, T0)); - REQUIRE(isClose(Prhot2, igra.PressureFromDensityTemperature(rho_t2, T0), 1.e-12)); - // check pressure below and beyond ramp matches unmodified ideal gas - INFO("reference P(0.8*r0, T0): " - << 0.8 * r0 * gm1 * Cv * T0 << " test P(0.8*r0, T0): " - << igra.PressureFromDensityTemperature(0.8 * r0, T0)); - REQUIRE(isClose(0.8 * r0 * gm1 * Cv * T0, - igra.PressureFromDensityTemperature(0.8 * r0, T0), 1.e-12)); - INFO("reference P(1.2*r1, T0): " - << 1.2 * r1 * gm1 * Cv * T0 << " test P(1.2*r1, T0): " - << igra.PressureFromDensityTemperature(1.2 * r1, T0)); - REQUIRE(isClose(1.2 * r1 * gm1 * Cv * T0, - igra.PressureFromDensityTemperature(1.2 * r1, T0), 1.e-12)); - } - THEN("We obtain correct ramp behavior in bmod(rho) for rho r1") { - // check bulk moduli on both pieces of ramp - INFO("reference bmod((r0+rmid)/2, T0): " - << bmodrt1 << " test bmod((r0+rmid)/2, T0): " - << igra.BulkModulusFromDensityTemperature(rho_t1, T0)); - REQUIRE( - isClose(bmodrt1, igra.BulkModulusFromDensityTemperature(rho_t1, T0), 1.e-12)); - INFO("reference bmod((rmid+r1)/2, T0): " - << bmodrt2 << " test bmod((rmid+r1)/2, T0): " - << igra.BulkModulusFromDensityTemperature(rho_t2, T0)); - REQUIRE( - isClose(bmodrt2, igra.BulkModulusFromDensityTemperature(rho_t2, T0), 1.e-12)); - // check bulk modulus below and beyond ramp matches unmodified ideal gas - INFO("reference bmod(0.8*r0, T0): " - << 0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0 << " test bmod(0.8*r0, T0): " - << igra.BulkModulusFromDensityTemperature(0.8 * r0, T0)); - REQUIRE(isClose(0.8 * r0 * gm1 * (gm1 + 1.0) * Cv * T0, - igra.BulkModulusFromDensityTemperature(0.8 * r0, T0), 1.e-12)); - INFO("reference bmod(1.2*r1, T0): " - << 1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0 << " test bmod(1.2*r1, T0): " - << igra.BulkModulusFromDensityTemperature(1.2 * r1, T0)); - REQUIRE(isClose(1.2 * r1 * gm1 * (gm1 + 1.0) * Cv * T0, - igra.BulkModulusFromDensityTemperature(1.2 * r1, T0), 1.e-12)); - } - } -#endif // SINGULARITY_BUILD_CLOSURE - } -} - -SCENARIO("Relativistic EOS", "[EOSBuilder][RelativisticEOS][IdealGas]") { - GIVEN("Parameters for an ideal gas") { - constexpr Real Cv = 2.0; - constexpr Real gm1 = 0.5; - WHEN("We construct a relativistic IdealGas with EOSBuilder") { - EOSBuilder::EOSType type = EOSBuilder::EOSType::IdealGas; - EOSBuilder::modifiers_t modifiers; - EOSBuilder::params_t base_params, relativity_params; - base_params["Cv"].emplace(Cv); - base_params["gm1"].emplace(gm1); - relativity_params["cl"].emplace(1.0); - modifiers[EOSBuilder::EOSModifier::Relativistic] = relativity_params; - EOS eos = EOSBuilder::buildEOS(type, base_params, modifiers); - THEN("The EOS has finite sound speeds") { - constexpr Real rho = 1e3; - constexpr Real sie = 1e3; - Real bmod = eos.BulkModulusFromDensityInternalEnergy(rho, sie); - Real cs2 = bmod / rho; - REQUIRE(cs2 < 1); - } - } - } -} - -SCENARIO("EOS Unit System", "[EOSBuilder][UnitSystem][IdealGas]") { - GIVEN("Parameters for an ideal gas") { - constexpr Real Cv = 2.0; - constexpr Real gm1 = 0.5; - EOSBuilder::EOSType type = EOSBuilder::EOSType::IdealGas; - EOSBuilder::modifiers_t modifiers; - EOSBuilder::params_t base_params, units_params; - base_params["Cv"].emplace(Cv); - base_params["gm1"].emplace(gm1); - GIVEN("Units with a thermal unit system") { - constexpr Real rho_unit = 1e1; - constexpr Real sie_unit = 1e-1; - constexpr Real temp_unit = 123; - WHEN("We construct an IdealGas with EOSBuilder") { - units_params["rho_unit"].emplace(rho_unit); - units_params["sie_unit"].emplace(sie_unit); - units_params["temp_unit"].emplace(temp_unit); - modifiers[EOSBuilder::EOSModifier::UnitSystem] = units_params; - EOS eos = EOSBuilder::buildEOS(type, base_params, modifiers); - THEN("Units cancel out for an ideal gas") { - Real rho = 1e3; - Real sie = 1e3; - Real P = eos.PressureFromDensityInternalEnergy(rho, sie); - Real Ptrue = gm1 * rho * sie; - REQUIRE(std::abs(P - Ptrue) / Ptrue < 1e-3); - } - } - } - GIVEN("Units with length and time units") { - constexpr Real time_unit = 456; - constexpr Real length_unit = 1e2; - constexpr Real mass_unit = 1e6; - constexpr Real temp_unit = 789; - WHEN("We construct an IdealGas with EOSBuilder") { - units_params["use_length_time"].emplace(true); - units_params["time_unit"].emplace(time_unit); - units_params["length_unit"].emplace(length_unit); - units_params["mass_unit"].emplace(mass_unit); - units_params["temp_unit"].emplace(temp_unit); - modifiers[EOSBuilder::EOSModifier::UnitSystem] = units_params; - EOS eos = EOSBuilder::buildEOS(type, base_params, modifiers); - THEN("Units cancel out for an ideal gas") { - Real rho = 1e3; - Real sie = 1e3; - Real P = eos.PressureFromDensityInternalEnergy(rho, sie); - Real Ptrue = gm1 * rho * sie; - REQUIRE(std::abs(P - Ptrue) / Ptrue < 1e-3); - } - } - } - } -} - -SCENARIO("Ideal gas entropy", "[IdealGas][Entropy]") { - GIVEN("Parameters for an ideal gas with entropy reference states") { - // Create ideal gas EOS ojbect - constexpr Real Cv = 5.0; - constexpr Real gm1 = 0.4; - constexpr Real EntropyT0 = 100; - constexpr Real EntropyRho0 = 1e-03; - EOS host_eos = IdealGas(gm1, Cv, EntropyT0, EntropyRho0); - THEN("The entropy at the reference state should be zero") { - auto entropy = host_eos.EntropyFromDensityTemperature(EntropyRho0, EntropyT0); - INFO("Entropy should be zero but it is " << entropy); - CHECK(isClose(entropy, 0.0, 1.e-14)); - } - GIVEN("A state at the reference temperature and a density whose cube root is the " - "reference density") { - constexpr Real T = EntropyT0; - constexpr Real rho = 0.1; // rho**3 = EntropyRho0 - THEN("The entropy should be 2. / 3. * gm1 * Cv * log(EntropyRho0)") { - const Real entropy_true = 2. / 3. * gm1 * Cv * log(EntropyRho0); - auto entropy = host_eos.EntropyFromDensityTemperature(rho, T); - INFO("Entropy: " << entropy << " True entropy: " << entropy_true); - CHECK(isClose(entropy, entropy_true, 1e-12)); - } - } - GIVEN("A state at the reference density and a temperature whose square is the " - "reference temperature") { - constexpr Real T = 10; // T**2 = EntropyT0 - constexpr Real rho = EntropyRho0; - THEN("The entropy should be -1. / 2. * Cv * log(EntropyT0)") { - const Real entropy_true = -1. / 2. * Cv * log(EntropyT0); - auto entropy = host_eos.EntropyFromDensityTemperature(rho, T); - INFO("Entropy: " << entropy << " True entropy: " << entropy_true); - CHECK(isClose(entropy, entropy_true, 1e-12)); - } - } - } -} - -// A non-standard pattern where we do a reduction -class CheckPofRE { - public: - CheckPofRE(Real *P, Real *rho, Real *sie, int N) : P_(P), rho_(rho), sie_(sie), N_(N) {} - template - void operator()(const T &eos) { - Real *P = P_; - Real *rho = rho_; - Real *sie = sie_; - portableReduce( - "MyCheckPofRE", 0, N_, - PORTABLE_LAMBDA(const int i, int &nw) { - nw += !(isClose(P[i], - eos.PressureFromDensityInternalEnergy(rho[i], sie[i], nullptr), - 1e-15)); - }, - nwrong); - } - - // must be initialized to zero because portableReduce is a simple for loop on host - int nwrong = 0; - - private: - int N_; - Real *P_; - Real *rho_; - Real *sie_; -}; -SCENARIO("Ideal gas vector Evaluate call", "[IdealGas][Evaluate]") { - GIVEN("An ideal gas, and some device memory") { - constexpr Real Cv = 2.0; - constexpr Real gm1 = 0.5; - constexpr int N = 100; - constexpr Real rhomin = 1; - constexpr Real rhomax = 11; - constexpr Real drho = (rhomax - rhomin) / (N - 1); - constexpr Real siemin = 1; - constexpr Real siemax = 11; - constexpr Real dsie = (siemax - siemin) / (N - 1); - - EOS eos_host = IdealGas(gm1, Cv); - auto eos_device = eos_host.GetOnDevice(); - - Real *P = (Real *)PORTABLE_MALLOC(N * sizeof(Real)); - Real *rho = (Real *)PORTABLE_MALLOC(N * sizeof(Real)); - Real *sie = (Real *)PORTABLE_MALLOC(N * sizeof(Real)); - - WHEN("We set P, rho, sie via the scalar API") { - portableFor( - "Set rho and sie", 0, N, PORTABLE_LAMBDA(const int i) { - rho[i] = rhomin + drho * i; // just some diagonal in the 2d rho/sie space - sie[i] = siemin + dsie * i; - P[i] = eos_device.PressureFromDensityInternalEnergy(rho[i], sie[i]); - }); - THEN("The vector Evaluate API can be used to compare") { - CheckPofRE my_op(P, rho, sie, N); - eos_device.Evaluate(my_op); - REQUIRE(my_op.nwrong == 0); - } - } - - eos_device.Finalize(); - PORTABLE_FREE(P); - PORTABLE_FREE(rho); - PORTABLE_FREE(sie); - } -} - -SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { - - GIVEN("Parameters for an ideal gas") { - // Create ideal gas EOS ojbect - constexpr Real Cv = 5.0; - constexpr Real gm1 = 0.4; - EOS host_eos = IdealGas(gm1, Cv); - EOS eos = host_eos.GetOnDevice(); - - GIVEN("Energies and densities") { - constexpr int num = 3; -#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; - Real *v_density = density.data(); - Real *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.0; - v_density[1] = 2.0; - v_density[2] = 5.0; - v_energy[0] = 5.0; - v_energy[1] = 10.0; - v_energy[2] = 15.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 - constexpr std::array pressure_true{2.0, 8.0, 30.0}; - constexpr std::array temperature_true{1., 2., 3.}; - constexpr std::array bulkmodulus_true{2.8, 11.2, 42.}; - constexpr std::array heatcapacity_true{Cv, Cv, Cv}; - constexpr std::array gruneisen_true{gm1, gm1, gm1}; - - // Gold standard entropy doesn't produce round numbers so we need to - // calculate it from the device views so this requires a bit more work -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::View v_entropy_true("True Entropy"); -#else - std::array entropy_true; - Real *v_entropy_true = entropy_true.data(); -#endif - constexpr Real P0 = 1e6; // microbar - constexpr Real T0 = 293; // K - constexpr Real rho0 = P0 / (gm1 * Cv * T0); // g/cm^3 - portableFor( - "Calculate true entropy", 0, num, PORTABLE_LAMBDA(const int i) { - v_entropy_true[i] = - Cv * log(v_energy[i] / Cv / T0) + gm1 * Cv * log(rho0 / v_density[i]); - }); -#ifdef PORTABILITY_STRATEGY_KOKKOS - auto entropy_true = Kokkos::create_mirror_view(v_entropy_true); - Kokkos::deep_copy(entropy_true, v_entropy_true); -#endif - -#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_heatcapacity("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_heatcapacity = Kokkos::create_mirror_view(v_heatcapacity); - 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_heatcapacity; - 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_heatcapacity = h_heatcapacity.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 " - "temperature") { - 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 pressure") { - array_compare(num, density, energy, h_pressure, pressure_true, "Density", - "Energy"); - } - } - - WHEN("An 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 entropy") { - array_compare(num, density, energy, h_entropy, entropy_true, "Density", - "Energy"); - } - } - - WHEN("A C_v(rho, e) lookup is performed") { - eos.SpecificHeatFromDensityInternalEnergy(v_density, v_energy, v_heatcapacity, - num); -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::fence(); - Kokkos::deep_copy(h_heatcapacity, v_heatcapacity); -#endif // PORTABILITY_STRATEGY_KOKKOS - THEN("The returned C_v(rho, e) should be constant") { - array_compare(num, density, energy, h_heatcapacity, heatcapacity_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 bulk " - "modulus") { - array_compare(num, density, energy, h_bulkmodulus, bulkmodulus_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 constant") { - array_compare(num, density, energy, h_gruneisen, gruneisen_true, "Density", - "Energy"); - } - } - } - GIVEN("Densities and temperatures") { - constexpr int num = 3; -#ifdef PORTABILITY_STRATEGY_KOKKOS - // Create Kokkos views on device for the input arrays - Kokkos::View v_density("density"); - Kokkos::View v_temperature("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 temperature; - Real *v_density = density.data(); - Real *v_temperature = temperature.data(); -#endif // PORTABILITY_STRATEGY_KOKKOS - - // Populate the input arrays - portableFor( - "Initialize density and energy", 0, 1, PORTABLE_LAMBDA(int i) { - v_density[0] = 1.0; - v_density[1] = 2.0; - v_density[2] = 5.0; - v_temperature[0] = 50.0; - v_temperature[1] = 100.0; - v_temperature[2] = 150.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 - constexpr std::array energy_true{250., 500., 750.}; - constexpr std::array pressure_true{100., 400., 1500.}; - constexpr std::array bulkmodulus_true{140., 560., 2100.}; - constexpr std::array heatcapacity_true{Cv, Cv, Cv}; - constexpr std::array gruneisen_true{gm1, gm1, gm1}; - - // Gold standard entropy doesn't produce round numbers so we need to - // calculate it from the device views so this requires a bit more work -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::View v_entropy_true("True Entropy"); -#else - std::array entropy_true; - Real *v_entropy_true = entropy_true.data(); -#endif - constexpr Real P0 = 1e6; // microbar - constexpr Real T0 = 293; // K - constexpr Real rho0 = P0 / (gm1 * Cv * T0); // g/cm^3 - portableFor( - "Calculate true entropy", 0, num, PORTABLE_LAMBDA(const int i) { - v_entropy_true[i] = - Cv * log(v_temperature[i] / T0) + gm1 * Cv * log(rho0 / v_density[i]); - }); -#ifdef PORTABILITY_STRATEGY_KOKKOS - auto entropy_true = Kokkos::create_mirror_view(v_entropy_true); - Kokkos::deep_copy(entropy_true, v_entropy_true); -#endif - -#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_heatcapacity("Cv"); - Kokkos::View v_bulkmodulus("bmod"); - Kokkos::View v_gruneisen("Gamma"); - 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_heatcapacity = Kokkos::create_mirror_view(v_heatcapacity); - 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_energy; - std::array h_pressure; - std::array h_entropy; - std::array h_heatcapacity; - std::array h_bulkmodulus; - std::array h_gruneisen; - // 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_heatcapacity = h_heatcapacity.data(); - auto v_bulkmodulus = h_bulkmodulus.data(); - auto v_gruneisen = h_gruneisen.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 energy") { - array_compare(num, density, temperature, h_energy, energy_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 pressure") { - array_compare(num, density, temperature, h_pressure, pressure_true, "Density", - "Temperature"); - } - } - - WHEN("An S(rho, T) lookup is performed") { - eos.EntropyFromDensityTemperature(v_density, v_temperature, 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, T) should be equal to the true entropy") { - array_compare(num, density, temperature, h_entropy, entropy_true, "Density", - "Temperature"); - } - } - - WHEN("A C_v(rho, T) lookup is performed") { - eos.SpecificHeatFromDensityTemperature(v_density, v_temperature, v_heatcapacity, - num); -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::fence(); - Kokkos::deep_copy(h_heatcapacity, v_heatcapacity); -#endif // PORTABILITY_STRATEGY_KOKKOS - THEN("The returned C_v(rho, T) should be constant") { - array_compare(num, density, temperature, h_heatcapacity, heatcapacity_true, - "Density", "Temperature"); - } - } - - WHEN("A B_S(rho, T) lookup is performed") { - eos.BulkModulusFromDensityTemperature(v_density, v_temperature, 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, T) should be equal to the true bulk " - "modulus") { - array_compare(num, density, temperature, h_bulkmodulus, bulkmodulus_true, - "Density", "Temperature"); - } - } - - WHEN("A Gamma(rho, T) lookup is performed") { - eos.GruneisenParamFromDensityTemperature(v_density, v_temperature, 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, T) should be constant") { - array_compare(num, density, temperature, h_gruneisen, gruneisen_true, "Density", - "Temperature"); - } - } - } - } -} - -#ifdef SPINER_USE_HDF -#ifdef SINGULARITY_TEST_SESAME -#ifdef SINGULARITY_USE_EOSPAC -SCENARIO("SpinerEOS depends on Rho and T", "[SpinerEOS],[DependsRhoT][EOSPAC]") { - - GIVEN("SpinerEOS and EOSPAC EOS for steel can be initialized with matid") { - EOS steelEOS_host_polymorphic = SpinerEOSDependsRhoT(eosName, steelID); - EOS steelEOS = steelEOS_host_polymorphic.GetOnDevice(); - auto steelEOS_host = steelEOS_host_polymorphic.get(); - - EOS eospac = EOSPAC(steelID); - - THEN("The correct metadata is read in") { - REQUIRE(steelEOS_host.matid() == steelID); - - AND_THEN("We can get a reference density and temperature") { - Real rho, T, sie, P, cv, bmod, dpde, dvdt; - Real rho_pac, T_pac, sie_pac, P_pac, cv_pac, bmod_pac, dpde_pac, dvdt_pac; - steelEOS_host.ValuesAtReferenceState(rho, T, sie, P, cv, bmod, dpde, dvdt); - eospac.ValuesAtReferenceState(rho_pac, T_pac, sie_pac, P_pac, cv_pac, bmod_pac, - dpde_pac, dvdt_pac); - REQUIRE(isClose(rho, rho_pac)); - REQUIRE(isClose(T, T_pac)); - } - - // TODO: this needs to be a much more rigorous test - AND_THEN("Quantities can be read from density and temperature") { - const Real sie_pac = eospac.InternalEnergyFromDensityTemperature(1e0, 1e6); - - int nw_ie{0}; -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::fence(); - using atomic_view = Kokkos::MemoryTraits; - Kokkos::View n_wrong_ie("wrong_ie"); -#else - PortableMDArray n_wrong_ie(&nw_ie, 1); -#endif - portableFor( - "calc ie's steel", 0, 100, PORTABLE_LAMBDA(const int &i) { - const double ie = steelEOS.InternalEnergyFromDensityTemperature(1e0, 1e6); - if (!isClose(ie, sie_pac)) n_wrong_ie() += 1; - }); -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(nw_ie, n_wrong_ie); -#endif - REQUIRE(nw_ie == 0); - } - - AND_THEN("rho(P,T) correct for P=1atm, T=freezing") { - Real T = 273; // Kelvin - Real P = 1e6; // barye - Real rho, sie; // output - Real rho_pac, sie_pac; - std::vector lambda(steelEOS_host_polymorphic.nlambda()); - steelEOS_host_polymorphic.DensityEnergyFromPressureTemperature( - P, T, lambda.data(), rho, sie); - eospac.DensityEnergyFromPressureTemperature(P, T, nullptr, rho_pac, sie_pac); - REQUIRE(isClose(rho, rho_pac)); - } - } - // Failing to call finalize leads to a memory leak, - // but otherwise behaviour is as expected. - steelEOS_host_polymorphic.Finalize(); // host and device must be - // finalized separately. - steelEOS.Finalize(); // cleans up memory on device. - } - - GIVEN("SpinerEOS and EOSPAC for air can be initialized with matid") { - SpinerEOSDependsRhoT airEOS_host = SpinerEOSDependsRhoT(eosName, airID); - EOS airEOS = airEOS_host.GetOnDevice(); - EOS eospac = EOSPAC(airID); - constexpr Real rho = 1e0; - constexpr Real sie = 2.5e-4; - THEN("We can get a reference state") { - Real rho, T, sie, P, cv, bmod, dpde, dvdt; - Real rho_pac, T_pac, sie_pac, P_pac, cv_pac, bmod_pac, dpde_pac, dvdt_pac; - airEOS_host.ValuesAtReferenceState(rho, T, sie, P, cv, bmod, dpde, dvdt); - eospac.ValuesAtReferenceState(rho_pac, T_pac, sie_pac, P_pac, cv_pac, bmod_pac, - dpde_pac, dvdt_pac); - REQUIRE(isClose(rho, rho_pac)); - REQUIRE(isClose(T, T_pac)); - } - THEN("Quantities of rho and sie look sane on both host and device") { - const Real gm1_host = airEOS_host.GruneisenParamFromDensityInternalEnergy(rho, sie); - const Real T_host = airEOS_host.TemperatureFromDensityInternalEnergy(rho, sie); - const Real cv_host = airEOS_host.SpecificHeatFromDensityInternalEnergy(rho, sie); - - int nw_gm1{0}, nw_cv{0}; -#ifdef PORTABILITY_STRATEGY_KOKKOS - using atomic_view = Kokkos::MemoryTraits; - Kokkos::fence(); - Kokkos::View n_wrong_gm1("wrong_gm1"); - Kokkos::View n_wrong_cv("wrong_cv"); -#else - PortableMDArray n_wrong_gm1(&nw_gm1, 1); - PortableMDArray n_wrong_cv(&nw_cv, 1); -#endif - portableFor( - "calc gm1 and cv", 0, 100, PORTABLE_LAMBDA(const int &i) { - const Real gm1 = airEOS.GruneisenParamFromDensityInternalEnergy(rho, sie); - const Real cv = airEOS.SpecificHeatFromDensityInternalEnergy(rho, sie); - if (!isClose(gm1, gm1_host)) n_wrong_gm1() += 1; - if (!isClose(cv, cv_host)) n_wrong_cv() += 1; - }); -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(nw_gm1, n_wrong_gm1); - Kokkos::deep_copy(nw_cv, n_wrong_cv); -#endif - REQUIRE(nw_gm1 == 0); - REQUIRE(nw_cv == 0); - } - AND_THEN("P(rho, sie) correct for extrapolation regime") { - Real rho = 1; - Real sie = 2.43e16; - Real P_pac = eospac.PressureFromDensityInternalEnergy(rho, sie); - Real P_spi = airEOS_host.PressureFromDensityInternalEnergy(rho, sie); - REQUIRE(isClose(P_pac, P_spi)); - } - airEOS_host.Finalize(); - airEOS.Finalize(); - } - - GIVEN("EOS initialized with matid") { - SpinerEOSDependsRhoT eos_spiner = SpinerEOSDependsRhoT(eosName, DTID); - EOS eos_eospac = EOSPAC(DTID); - Real P = 1e8; - Real rho = 1.28e-3; - THEN("Inversion for T(rho,P) works on host") { - Real T, sie, cv, bmod; - Real T_pac, sie_pac, cv_pac, bmod_pac; - const unsigned long output = - (thermalqs::temperature | thermalqs::specific_internal_energy | - thermalqs::specific_heat | thermalqs::bulk_modulus); - eos_spiner.FillEos(rho, T, sie, P, cv, bmod, output); - eos_eospac.FillEos(rho, T_pac, sie_pac, P, cv_pac, bmod_pac, output); - REQUIRE(isClose(T, T_pac)); - REQUIRE(isClose(sie, sie_pac)); - REQUIRE(isClose(cv, cv_pac)); - } - eos_spiner.Finalize(); - } - - GIVEN("EOS initialized with matid") { - SpinerEOSDependsRhoT eos_spiner = SpinerEOSDependsRhoT(eosName, gID); - EOS eos_eospac = EOSPAC(gID); - THEN("PT lookup works on the host") { - Real P = 1e6; // cgs - Real T = 0.025 / ev2k; // K - Real rho, sie; - Real rho_pac, sie_pac; - std::vector lambda(eos_spiner.nlambda()); - eos_spiner.DensityEnergyFromPressureTemperature(P, T, lambda.data(), rho, sie); - eos_eospac.DensityEnergyFromPressureTemperature(P, T, lambda.data(), rho_pac, - sie_pac); - REQUIRE(isClose(rho, rho_pac)); - } - eos_spiner.Finalize(); - } -} - -// Disabling these tests for now as the DependsRhoSie code is not well-maintained -SCENARIO("SpinerEOS depends on rho and sie", "[SpinerEOS],[DependsRhoSie]") { - - GIVEN("SpinerEOSes for steel can be initialised with matid") { - SpinerEOSDependsRhoSie steelEOS_host(eosName, steelID); - EOS steelEOS = steelEOS_host.GetOnDevice(); - THEN("The correct metadata is read in") { - REQUIRE(steelEOS_host.matid() == steelID); - - int nw_ie2{0}, nw_te2{0}; -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::fence(); - using atomic_view = Kokkos::MemoryTraits; - Kokkos::View n_wrong_ie("wrong_ie"); - Kokkos::View n_wrong_te("wrong_te"); -#else - PortableMDArray n_wrong_ie(&nw_ie2, 1); - PortableMDArray n_wrong_te(&nw_te2, 1); -#endif - portableFor( - "calc ie's steel 2", 0, 100, PORTABLE_LAMBDA(const int &i) { - const Real ie{steelEOS.InternalEnergyFromDensityTemperature(1e0, 1e6)}; - const Real te{steelEOS.TemperatureFromDensityInternalEnergy(1e0, 1e12)}; - if (!isClose(ie, 4.96415e13)) n_wrong_ie() += 1; - if (!isClose(te, 75594.6)) n_wrong_te() += 1; - }); -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(nw_ie2, n_wrong_ie); - Kokkos::deep_copy(nw_te2, n_wrong_te); -#endif - REQUIRE(nw_ie2 == 0); - REQUIRE(nw_te2 == 0); - } - // this can be removed with with reference counting or other tricks - steelEOS_host.Finalize(); // cleans up host memory - steelEOS.Finalize(); // cleans up device memory - } -} - -SCENARIO("EOS Builder and SpinerEOS", - "[SpinerEOS],[EOSBuilder],[GetOnDevice],[Finalize]") { - GIVEN("Parameters for shift and scale") { - constexpr Real shift = 0.0; - constexpr Real scale = 1.0; - WHEN("We construct a SpinerEOS with EOSBuilder") { - EOSBuilder::EOSType type = EOSBuilder::EOSType::SpinerEOSDependsRhoT; - EOSBuilder::modifiers_t modifiers; - EOSBuilder::params_t base_params, shifted_params, scaled_params; - base_params["filename"].emplace(eosName); - base_params["matid"].emplace(steelID); - shifted_params["shift"].emplace(shift); - scaled_params["scale"].emplace(scale); - modifiers[EOSBuilder::EOSModifier::Shifted] = shifted_params; - modifiers[EOSBuilder::EOSModifier::Scaled] = scaled_params; - EOS eosHost = EOSBuilder::buildEOS(type, base_params, modifiers); - THEN("The EOS is consistent.") { - REQUIRE( - isClose(4.96416e13, eosHost.InternalEnergyFromDensityTemperature(1e0, 1e6))); - } - WHEN("We get an EOS on device") { - EOS eosDevice = eosHost.GetOnDevice(); - THEN("EOS calls match raw access") { - int nw_bm{0}; -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::fence(); - using atomic_view = Kokkos::MemoryTraits; - Kokkos::View n_wrong_bm("wrong_bm"); -#else - PortableMDArray n_wrong_bm(&nw_bm, 1); -#endif - portableFor( - "calc ie's steel 3", 0, 100, PORTABLE_LAMBDA(const int &i) { - const Real bm{eosDevice.BulkModulusFromDensityTemperature(1e0, 1e6)}; - if (!isClose(bm, 2.55268e13)) n_wrong_bm() += 1; - }); -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(nw_bm, n_wrong_bm); -#endif - REQUIRE(nw_bm == 0); - } - eosDevice.Finalize(); - } - eosHost.Finalize(); - } - } -} -#endif // SINGULARITY_USE_EOSPAC -#endif // SINGULARITY_TEST_SESAME - -#ifdef SINGULARITY_TEST_STELLAR_COLLAPSE -SCENARIO("Test 3D reinterpolation to fast log grid", "[StellarCollapse]") { - WHEN("We generate a 3D DataBox of a 2D power law times a line") { - using singularity::StellarCollapse; - constexpr int N2 = 100; - constexpr int N1 = 101; - constexpr int N0 = 102; - StellarCollapse::Grid_t g2(0, 1, N2); - StellarCollapse::Grid_t g1(1, 3, N1); - StellarCollapse::Grid_t g0(2, 4, N0); - StellarCollapse::DataBox db(N2, N1, N0); - - db.setRange(2, g2); - db.setRange(1, g1); - db.setRange(0, g0); - - for (int i2 = 0; i2 < N2; ++i2) { - Real x2 = g2.x(i2); - for (int i1 = 0; i1 < N1; ++i1) { - Real lx1 = g1.x(i1); - for (int i0 = 0; i0 < N0; ++i0) { - Real lx0 = g0.x(i0); - db(i2, i1, i0) = std::log10(x2) + 2 * lx1 + 2 * lx0; - } - } - } - - THEN("The databox should interpolate values correctly") { - const Real x2 = 0.5; - const Real lx1 = 2; - const Real x1 = std::pow(10., lx1); - const Real lx0 = 3; - const Real x0 = std::pow(10., lx0); - REQUIRE(isClose(db.interpToReal(x2, lx1, lx0), std::log10(x2 * x1 * x1 * x0 * x0), - 1e-5)); - } - - THEN("We can re-interpolate to fast log space") { - StellarCollapse::DataBox scratch(N2, N1, N0); - StellarCollapse::dataBoxToFastLogs(db, scratch, true); - - AND_THEN("The fast-log gridded table contains correct ranges") { - REQUIRE(db.range(2) == g2); - REQUIRE(db.range(1).nPoints() == N1); - REQUIRE(isClose(db.range(1).min(), - singularity::FastMath::log10(std::pow(10, g1.min())), 1e-12)); - REQUIRE(isClose(db.range(1).max(), - singularity::FastMath::log10(std::pow(10, g1.max())), 1e-12)); - REQUIRE(db.range(0).nPoints() == N0); - REQUIRE(isClose(db.range(0).min(), - singularity::FastMath::log10(std::pow(10, g0.min())), 1e-12)); - REQUIRE(isClose(db.range(0).max(), - singularity::FastMath::log10(std::pow(10, g0.max())), 1e-12)); - } - - AND_THEN("The re-interpolated fast log is a sane number") {} - - AND_THEN("The fast-log table approximately interpolates the power law") { - const Real x2 = 0.5; - const Real x1 = 100; - const Real x0 = 1000; - const Real lx1 = singularity::FastMath::log10(x1); - const Real lx0 = singularity::FastMath::log10(x0); - const Real lval_interp = db.interpToReal(x2, lx1, lx0); - const Real val_interp = singularity::FastMath::pow10(lval_interp); - const Real val_true = x2 * x1 * x1 * x0 * x0; - const Real rel_diff = - 0.5 * std::abs(val_interp - val_true) / (val_true + val_interp); - REQUIRE(rel_diff <= 1e-3); - } - scratch.finalize(); - } - db.finalize(); - } -} - -SCENARIO("Stellar Collapse EOS", "[StellarCollapse][EOSBuilder]") { - using singularity::IdealGas; - using singularity::StellarCollapse; - const std::string savename = "stellar_collapse_ideal_2.sp5"; - GIVEN("A stellar collapse EOS") { - const std::string filename = "./goldfiles/stellar_collapse_ideal.h5"; - THEN("We can load the file") { // don't bother filtering bmod here. - StellarCollapse sc(filename, false, false); - AND_THEN("Some properties we expect for ideal gas hold") { - Real lambda[2]; - Real rho, t, sie, p, cv, b, dpde, dvdt; - sc.ValuesAtReferenceState(rho, t, sie, p, cv, b, dpde, dvdt, lambda); - Real yemin = sc.YeMin(); - Real yemax = sc.YeMax(); - int N = 123; - Real dY = (yemax - yemin) / (N + 1); - for (int i = 0; i < N; ++i) { - Real Ye = yemin + i * dY; - lambda[0] = Ye; - REQUIRE(isClose(sie, sc.InternalEnergyFromDensityTemperature(rho, t, lambda))); - Real Xa, Xh, Xn, Xp, Abar, Zbar; - sc.MassFractionsFromDensityTemperature(rho, t, Xa, Xh, Xn, Xp, Abar, Zbar, - lambda); - REQUIRE(isClose(Ye, Xp)); - } - Real rhomin = sc.rhoMin(); - Real rhomax = sc.rhoMax(); - Real drho = (rhomax - rhomin) / (N + 1); - for (int i = 0; i < N; ++i) { - Real rho = rhomin + i * drho; - REQUIRE(isClose(sie, sc.InternalEnergyFromDensityTemperature(rho, t, lambda))); - } - } - GIVEN("An Ideal Gas equation of state") { - constexpr Real gamma = 1.4; - constexpr Real mp = 1.67262171e-24; - constexpr Real kb = 1.3806505e-16; - constexpr Real Cv = kb / (mp * (gamma - 1)); // mean molecular weight = mp - IdealGas ig(gamma - 1, Cv); - auto ig_d = ig.GetOnDevice(); - THEN("The tabulated gamma Stellar Collapse and the gamma agree roughly") { - Real yemin = sc.YeMin(); - Real yemax = sc.YeMax(); - Real tmin = sc.TMin(); - Real tmax = sc.TMax(); - Real ltmin = std::log10(tmin); - Real ltmax = std::log10(tmax); - Real lrhomin = std::log10(sc.rhoMin()); - Real lrhomax = std::log10(sc.rhoMax()); - auto sc_d = sc.GetOnDevice(); - - int nwrong_h = 0; -#ifdef PORTABILITY_STRATEGY_KOKKOS - using atomic_view = Kokkos::MemoryTraits; - Kokkos::View nwrong_d("wrong"); -#else - PortableMDArray nwrong_d(&nwrong_h, 1); -#endif - - const int N = 123; - const Real dY = (yemax - yemin) / (N + 1); - const Real dlT = (ltmax - ltmin) / (N + 1); - const Real dlR = (lrhomax - lrhomin) / (N + 1); - portableFor( - "fill eos", 0, N, 0, N, 0, N, - PORTABLE_LAMBDA(const int &k, const int &j, const int &i) { - Real lambda[2]; - Real Ye = yemin + k * dY; - Real lT = ltmin + j * dlT; - Real lR = lrhomin + i * dlR; - Real T = std::pow(10., lT); - Real R = std::pow(10., lR); - Real e1, e2, p1, p2, cv1, cv2, b1, b2; - unsigned long output = (singularity::thermalqs::pressure | - singularity::thermalqs::specific_internal_energy | - singularity::thermalqs::specific_heat | - singularity::thermalqs::bulk_modulus); - lambda[0] = Ye; - - sc_d.FillEos(R, T, e1, p1, cv1, b1, output, lambda); - ig_d.FillEos(R, T, e2, p2, cv2, b2, output, lambda); - if (!isClose(e1, e2)) { - nwrong_d() += 1; - } - if (!isClose(p1, p2)) { - nwrong_d() += 1; - } - if (!isClose(cv1, cv2)) { - nwrong_d() += 1; - } - if (!isClose(b1, b2)) { - nwrong_d() += 1; - } - }); -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(nwrong_h, nwrong_d); -#endif - REQUIRE(nwrong_h == 0); - sc_d.Finalize(); - } - ig_d.Finalize(); - ig.Finalize(); - } - AND_THEN("We can save the file to SP5") { - sc.Save(savename); - AND_THEN("We can load the sp5 file") { - StellarCollapse sc2(savename, true); - AND_THEN("The two stellar collapse EOS's agree") { - - Real yemin = sc.YeMin(); - Real yemax = sc.YeMax(); - Real tmin = sc.TMin(); - Real tmax = sc.TMax(); - Real ltmin = std::log10(tmin); - Real ltmax = std::log10(tmax); - Real lrhomin = std::log10(sc.rhoMin()); - Real lrhomax = std::log10(sc.rhoMax()); - REQUIRE(yemin == sc2.YeMin()); - REQUIRE(yemax == sc2.YeMax()); - REQUIRE(sc.TMin() == sc2.TMin()); - REQUIRE(sc.TMax() == sc2.TMax()); - REQUIRE(isClose(lrhomin, std::log10(sc2.rhoMin()), 1e-12)); - REQUIRE(isClose(lrhomax, std::log10(sc2.rhoMax()), 1e-12)); - - auto sc1_d = sc.GetOnDevice(); - auto sc2_d = sc2.GetOnDevice(); - - int nwrong_h = 0; -#ifdef PORTABILITY_STRATEGY_KOKKOS - using atomic_view = Kokkos::MemoryTraits; - Kokkos::View nwrong_d("wrong"); -#else - PortableMDArray nwrong_d(&nwrong_h, 1); -#endif - - const int N = 123; - constexpr Real gamma = 1.4; - const Real dY = (yemax - yemin) / (N + 1); - const Real dlT = (ltmax - ltmin) / (N + 1); - const Real dlR = (lrhomax - lrhomin) / (N + 1); - portableFor( - "fill eos", 0, N, 0, N, 0, N, - PORTABLE_LAMBDA(const int &k, const int &j, const int &i) { - Real lambda[2]; - Real Ye = yemin + k * dY; - Real lT = ltmin + j * dlT; - Real lR = lrhomin + i * dlR; - Real T = std::pow(10., lT); - Real R = std::pow(10., lR); - Real e1, e2, p1, p2, cv1, cv2, b1, b2, s1, s2; - unsigned long output = - (singularity::thermalqs::pressure | - singularity::thermalqs::specific_internal_energy | - singularity::thermalqs::specific_heat | - singularity::thermalqs::bulk_modulus); - lambda[0] = Ye; - - sc1_d.FillEos(R, T, e1, p1, cv1, b1, output, lambda); - sc2_d.FillEos(R, T, e2, p2, cv2, b2, output, lambda); - // Fill entropy. Will need to change later. - s1 = sc1_d.EntropyFromDensityTemperature(R, T, lambda); - s2 = p2 * std::pow(R, -gamma); // ideal - if (!isClose(e1, e2)) nwrong_d() += 1; - if (!isClose(p1, p2)) nwrong_d() += 1; - if (!isClose(cv1, cv2)) nwrong_d() += 1; - if (!isClose(b1, b2)) nwrong_d() += 1; - if (!isClose(s1, s2)) nwrong_d() += 1; - }); -#ifdef PORTABILITY_STRATEGY_KOKKOS - Kokkos::deep_copy(nwrong_h, nwrong_d); -#endif - REQUIRE(nwrong_h == 0); - - sc1_d.Finalize(); - sc2_d.Finalize(); - } - sc2.Finalize(); - } - } - sc.Finalize(); - } - } -} -#endif // SINGULARITY_TEST_STELLAR_COLLAPSE -#endif // USE_HDF5 diff --git a/test/test_eos_vector.cpp b/test/test_eos_vector.cpp new file mode 100644 index 0000000000..b82eec4760 --- /dev/null +++ b/test/test_eos_vector.cpp @@ -0,0 +1,387 @@ +//------------------------------------------------------------------------------ +// © 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 +#include +#include +#include + +#ifndef CATCH_CONFIG_RUNNER +#include "catch2/catch.hpp" +#endif + +#include + +using singularity::EOS; +using singularity::IdealGas; + +namespace EOSBuilder = singularity::EOSBuilder; +namespace thermalqs = singularity::thermalqs; + +SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { + + GIVEN("Parameters for an ideal gas") { + // Create ideal gas EOS ojbect + constexpr Real Cv = 5.0; + constexpr Real gm1 = 0.4; + EOS host_eos = IdealGas(gm1, Cv); + EOS eos = host_eos.GetOnDevice(); + + GIVEN("Energies and densities") { + constexpr int num = 3; +#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; + Real *v_density = density.data(); + Real *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.0; + v_density[1] = 2.0; + v_density[2] = 5.0; + v_energy[0] = 5.0; + v_energy[1] = 10.0; + v_energy[2] = 15.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 + constexpr std::array pressure_true{2.0, 8.0, 30.0}; + constexpr std::array temperature_true{1., 2., 3.}; + constexpr std::array bulkmodulus_true{2.8, 11.2, 42.}; + constexpr std::array heatcapacity_true{Cv, Cv, Cv}; + constexpr std::array gruneisen_true{gm1, gm1, gm1}; + + // Gold standard entropy doesn't produce round numbers so we need to + // calculate it from the device views so this requires a bit more work +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View v_entropy_true("True Entropy"); +#else + std::array entropy_true; + Real *v_entropy_true = entropy_true.data(); +#endif + constexpr Real P0 = 1e6; // microbar + constexpr Real T0 = 293; // K + constexpr Real rho0 = P0 / (gm1 * Cv * T0); // g/cm^3 + portableFor( + "Calculate true entropy", 0, num, PORTABLE_LAMBDA(const int i) { + v_entropy_true[i] = + Cv * log(v_energy[i] / Cv / T0) + gm1 * Cv * log(rho0 / v_density[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + auto entropy_true = Kokkos::create_mirror_view(v_entropy_true); + Kokkos::deep_copy(entropy_true, v_entropy_true); +#endif + +#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_heatcapacity("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_heatcapacity = Kokkos::create_mirror_view(v_heatcapacity); + 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_heatcapacity; + 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_heatcapacity = h_heatcapacity.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 " + "temperature") { + 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 pressure") { + array_compare(num, density, energy, h_pressure, pressure_true, "Density", + "Energy"); + } + } + + WHEN("An 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 entropy") { + array_compare(num, density, energy, h_entropy, entropy_true, "Density", + "Energy"); + } + } + + WHEN("A C_v(rho, e) lookup is performed") { + eos.SpecificHeatFromDensityInternalEnergy(v_density, v_energy, v_heatcapacity, + num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_heatcapacity, v_heatcapacity); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned C_v(rho, e) should be constant") { + array_compare(num, density, energy, h_heatcapacity, heatcapacity_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 bulk " + "modulus") { + array_compare(num, density, energy, h_bulkmodulus, bulkmodulus_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 constant") { + array_compare(num, density, energy, h_gruneisen, gruneisen_true, "Density", + "Energy"); + } + } + } + GIVEN("Densities and temperatures") { + constexpr int num = 3; +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_density("density"); + Kokkos::View v_temperature("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 temperature; + Real *v_density = density.data(); + Real *v_temperature = temperature.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize density and energy", 0, 1, PORTABLE_LAMBDA(int i) { + v_density[0] = 1.0; + v_density[1] = 2.0; + v_density[2] = 5.0; + v_temperature[0] = 50.0; + v_temperature[1] = 100.0; + v_temperature[2] = 150.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 + constexpr std::array energy_true{250., 500., 750.}; + constexpr std::array pressure_true{100., 400., 1500.}; + constexpr std::array bulkmodulus_true{140., 560., 2100.}; + constexpr std::array heatcapacity_true{Cv, Cv, Cv}; + constexpr std::array gruneisen_true{gm1, gm1, gm1}; + + // Gold standard entropy doesn't produce round numbers so we need to + // calculate it from the device views so this requires a bit more work +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View v_entropy_true("True Entropy"); +#else + std::array entropy_true; + Real *v_entropy_true = entropy_true.data(); +#endif + constexpr Real P0 = 1e6; // microbar + constexpr Real T0 = 293; // K + constexpr Real rho0 = P0 / (gm1 * Cv * T0); // g/cm^3 + portableFor( + "Calculate true entropy", 0, num, PORTABLE_LAMBDA(const int i) { + v_entropy_true[i] = + Cv * log(v_temperature[i] / T0) + gm1 * Cv * log(rho0 / v_density[i]); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + auto entropy_true = Kokkos::create_mirror_view(v_entropy_true); + Kokkos::deep_copy(entropy_true, v_entropy_true); +#endif + +#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_heatcapacity("Cv"); + Kokkos::View v_bulkmodulus("bmod"); + Kokkos::View v_gruneisen("Gamma"); + 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_heatcapacity = Kokkos::create_mirror_view(v_heatcapacity); + 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_energy; + std::array h_pressure; + std::array h_entropy; + std::array h_heatcapacity; + std::array h_bulkmodulus; + std::array h_gruneisen; + // 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_heatcapacity = h_heatcapacity.data(); + auto v_bulkmodulus = h_bulkmodulus.data(); + auto v_gruneisen = h_gruneisen.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 energy") { + array_compare(num, density, temperature, h_energy, energy_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 pressure") { + array_compare(num, density, temperature, h_pressure, pressure_true, "Density", + "Temperature"); + } + } + + WHEN("An S(rho, T) lookup is performed") { + eos.EntropyFromDensityTemperature(v_density, v_temperature, 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, T) should be equal to the true entropy") { + array_compare(num, density, temperature, h_entropy, entropy_true, "Density", + "Temperature"); + } + } + + WHEN("A C_v(rho, T) lookup is performed") { + eos.SpecificHeatFromDensityTemperature(v_density, v_temperature, v_heatcapacity, + num); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_heatcapacity, v_heatcapacity); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned C_v(rho, T) should be constant") { + array_compare(num, density, temperature, h_heatcapacity, heatcapacity_true, + "Density", "Temperature"); + } + } + + WHEN("A B_S(rho, T) lookup is performed") { + eos.BulkModulusFromDensityTemperature(v_density, v_temperature, 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, T) should be equal to the true bulk " + "modulus") { + array_compare(num, density, temperature, h_bulkmodulus, bulkmodulus_true, + "Density", "Temperature"); + } + } + + WHEN("A Gamma(rho, T) lookup is performed") { + eos.GruneisenParamFromDensityTemperature(v_density, v_temperature, 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, T) should be constant") { + array_compare(num, density, temperature, h_gruneisen, gruneisen_true, "Density", + "Temperature"); + } + } + } + } +} diff --git a/test/test_math_utils.cpp b/test/test_math_utils.cpp new file mode 100644 index 0000000000..1e7c15c7bb --- /dev/null +++ b/test/test_math_utils.cpp @@ -0,0 +1,128 @@ +//------------------------------------------------------------------------------ +// © 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 +#include +#include + +#ifndef CATCH_CONFIG_RUNNER +#include "catch2/catch.hpp" +#endif + +#include + +SCENARIO("EOS Variant Type", "[Variant][EOS]") { + // print out the eos type + std::cout << demangle(typeid(EOS).name()) << std::endl; +} + +SCENARIO("Test that we can either throw an error on host or do nothing on device", + "[RequireMaybe]") { + // TODO(JMM): For whatever reason, the preprocessor does not like it if I + // call `PORTABLE_ALWAYS_THROW_OR_ABORT + REQUIRE_MAYBE_THROWS(PORTABLE_ALWAYS_THROW_OR_ABORT("Error message")); +} + +SCENARIO("Test that fast logs are invertible and run on device", "[FastMath]") { + GIVEN("A set of values to invert over a large dynamic range") { + constexpr Real LXMIN = -20; + constexpr Real LXMAX = 32; + constexpr int NX = 1000; + constexpr Real DLX = (LXMAX - LXMIN) / (NX - 1); + Real *x = (Real *)PORTABLE_MALLOC(NX * sizeof(Real)); + portableFor( + "Set x values", 0, NX, PORTABLE_LAMBDA(const int i) { + const Real lx = LXMIN + i * DLX; + x[i] = std::pow(10., lx); + }); + THEN("The fast exp of the fast log returns the original") { + int nw_ie = 0; +#ifdef PORTABILITY_STRATEGY_KOKKOS + using atomic_view = Kokkos::MemoryTraits; + Kokkos::View n_wrong_ie("wrong_ie"); +#else + PortableMDArray n_wrong_ie(&nw_ie, 1); +#endif + portableFor( + "try out the fast math", 0, NX, PORTABLE_LAMBDA(const int i) { + constexpr Real machine_eps = std::numeric_limits::epsilon(); + constexpr Real acceptable_err = 100 * machine_eps; + const Real lx = singularity::FastMath::log10(x[i]); + const Real elx = singularity::FastMath::pow10(lx); + const Real rel_err = 2.0 * std::abs(x[i] - elx) / + (std::abs(x[i]) + std::abs(elx) + machine_eps); + n_wrong_ie() += (rel_err > acceptable_err); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(nw_ie, n_wrong_ie); +#endif + REQUIRE(nw_ie == 0); + } + PORTABLE_FREE(x); + } +} + +SCENARIO("Rudimentary test of the root finder", "[RootFinding1D]") { + + GIVEN("Root finding") { + using namespace RootFinding1D; + + THEN("A root can be found for shift = 1, scale = 2, offset = 0.5") { + int ntimes = 100; + Real guess = 0; + Real root; + Status status; + Real shift = 1; + Real scale = 2; + Real offset = 0.5; + + auto f = PORTABLE_LAMBDA(const Real x) { return myAtan(x, shift, scale, offset); }; + Status *statusesp = (Status *)PORTABLE_MALLOC(ntimes * sizeof(Status)); + Real *rootsp = (Real *)PORTABLE_MALLOC(ntimes * sizeof(Real)); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View> statuses(statusesp, + ntimes); + Kokkos::View> roots(rootsp, ntimes); +#else + PortableMDArray statuses; + PortableMDArray roots; + statuses.NewPortableMDArray(statusesp, ntimes); + roots.NewPortableMDArray(rootsp, ntimes); +#endif + portableFor( + "find roots", 0, ntimes, PORTABLE_LAMBDA(const int i) { + statuses(i) = regula_falsi(f, 0, guess, -1, 3, 1e-10, 1e-10, roots(i)); + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View s_copy(statuses, 0); + Kokkos::View r_copy(roots, 0); + Kokkos::deep_copy(root, r_copy); + Kokkos::deep_copy(status, s_copy); +#else + status = statuses(ntimes - 1); + root = roots(ntimes - 1); +#endif + + PORTABLE_FREE(statusesp); + PORTABLE_FREE(rootsp); + REQUIRE(status == Status::SUCCESS); + REQUIRE(isClose(root, 0.744658)); + } + } +} From cbf860b8994455130596cb9e19d466962875ee70 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 11 Oct 2023 20:16:34 -0600 Subject: [PATCH 2/7] add CATCH_CONFIG_FAST_COMPILE and split executable into 3 --- test/CMakeLists.txt | 38 ++++++++++++++++++++++-------- test/catch2_define.cpp | 1 + test/eos_unit_test_helpers.hpp | 3 ++- test/test_eos_gruneisen.cpp | 3 ++- test/test_eos_helmholtz.cpp | 3 ++- test/test_eos_ideal.cpp | 5 ++-- test/test_eos_modifiers.cpp | 3 ++- test/test_eos_noble_abel.cpp | 3 ++- test/test_eos_sap_polynomial.cpp | 3 ++- test/test_eos_stellar_collapse.cpp | 3 ++- test/test_eos_stiff.cpp | 4 ++-- test/test_eos_tabulated.cpp | 4 ++-- test/test_eos_vector.cpp | 3 ++- test/test_eos_vinet.cpp | 4 ++-- test/test_math_utils.cpp | 3 ++- 15 files changed, 55 insertions(+), 28 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5f71e9a5a4..d30b5b0b8e 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -12,40 +12,58 @@ # ------------------------------------------------------------------------------ add_executable( - eos_unit_tests + eos_analytic_unit_tests catch2_define.cpp eos_unit_test_helpers.hpp test_eos_ideal.cpp test_eos_gruneisen.cpp test_eos_sap_polynomial.cpp - test_eos_vinet.cpp test_eos_noble_abel.cpp test_eos_stiff.cpp - test_eos_helmholtz.cpp + ) + +add_executable( + eos_infrastructure_tests + catch2_define.cpp + eos_unit_test_helpers.hpp test_eos_modifiers.cpp test_eos_vector.cpp + test_math_utils.cpp + ) + +add_executable( + eos_tabulated_unit_tests + catch2_define.cpp + eos_unit_test_helpers.hpp + test_eos_vinet.cpp + test_eos_helmholtz.cpp test_eos_tabulated.cpp test_eos_stellar_collapse.cpp - test_math_utils.cpp ) 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_unit_tests PRIVATE SINGULARITY_TEST_HELMHOLTZ SINGULARITY_USE_HELMHOLTZ) + target_compile_definitions(eos_tabulated_unit_tests PRIVATE SINGULARITY_TEST_HELMHOLTZ SINGULARITY_USE_HELMHOLTZ) endif() if(SINGULARITY_TEST_SESAME) - target_compile_definitions(eos_unit_tests PRIVATE SINGULARITY_TEST_SESAME) + target_compile_definitions(eos_tabulated_unit_tests PRIVATE SINGULARITY_TEST_SESAME) endif() if(SINGULARITY_TEST_STELLAR_COLLAPSE) - target_compile_definitions(eos_unit_tests + target_compile_definitions(eos_tabulated_unit_tests PRIVATE SINGULARITY_TEST_STELLAR_COLLAPSE) endif() -target_link_libraries(eos_unit_tests PRIVATE Catch2::Catch2 - singularity-eos::singularity-eos) +target_link_libraries(eos_analytic_unit_tests PRIVATE Catch2::Catch2 + singularity-eos::singularity-eos) +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) include(Catch) -catch_discover_tests(eos_unit_tests PROPERTIES TIMEOUT 60) +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(SINGULARITY_USE_EOSPAC AND SINGULARITY_TEST_SESAME) add_executable(compare_to_eospac compare_to_eospac.cpp) diff --git a/test/catch2_define.cpp b/test/catch2_define.cpp index 65437e8f74..c06cf7a5c0 100644 --- a/test/catch2_define.cpp +++ b/test/catch2_define.cpp @@ -16,6 +16,7 @@ #ifndef CATCH_CONFIG_RUNNER #define CATCH_CONFIG_RUNNER +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index cda089b013..a7c40d214a 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -15,7 +15,8 @@ #ifndef _SINGULARITY_EOS_TEST_TEST_HELPERS_ #define _SINGULARITY_EOS_TEST_TEST_HELPERS_ -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif #include diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index c91ea9bd37..04936f106b 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -17,7 +17,8 @@ #include #include #include -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/test_eos_helmholtz.cpp b/test/test_eos_helmholtz.cpp index b1048ee5c8..319463cc61 100644 --- a/test/test_eos_helmholtz.cpp +++ b/test/test_eos_helmholtz.cpp @@ -22,7 +22,8 @@ #include #include -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp index 117c61e914..37a1818772 100644 --- a/test/test_eos_ideal.cpp +++ b/test/test_eos_ideal.cpp @@ -12,11 +12,9 @@ // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ -#include #include #include #include -#include // debug #include #include @@ -28,7 +26,8 @@ #include #endif -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/test_eos_modifiers.cpp b/test/test_eos_modifiers.cpp index 2140001830..74674f668d 100644 --- a/test/test_eos_modifiers.cpp +++ b/test/test_eos_modifiers.cpp @@ -20,7 +20,8 @@ #include #include -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/test_eos_noble_abel.cpp b/test/test_eos_noble_abel.cpp index 5d6b4264bd..e2a00e8414 100644 --- a/test/test_eos_noble_abel.cpp +++ b/test/test_eos_noble_abel.cpp @@ -17,7 +17,8 @@ #include #include #include -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/test_eos_sap_polynomial.cpp b/test/test_eos_sap_polynomial.cpp index d88ec89b10..b916a138e9 100644 --- a/test/test_eos_sap_polynomial.cpp +++ b/test/test_eos_sap_polynomial.cpp @@ -17,7 +17,8 @@ #include #include #include -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/test_eos_stellar_collapse.cpp b/test/test_eos_stellar_collapse.cpp index 8979554fab..859031d5b5 100644 --- a/test/test_eos_stellar_collapse.cpp +++ b/test/test_eos_stellar_collapse.cpp @@ -29,7 +29,8 @@ #include #endif -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/test_eos_stiff.cpp b/test/test_eos_stiff.cpp index 548af96a93..2c8bb319d9 100644 --- a/test/test_eos_stiff.cpp +++ b/test/test_eos_stiff.cpp @@ -16,8 +16,8 @@ #include #include #include -#include -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/test_eos_tabulated.cpp b/test/test_eos_tabulated.cpp index bce04d0d95..dfaea1a05d 100644 --- a/test/test_eos_tabulated.cpp +++ b/test/test_eos_tabulated.cpp @@ -17,7 +17,6 @@ #include #include #include // debug -#include #include #include @@ -29,7 +28,8 @@ #include #endif -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/test_eos_vector.cpp b/test/test_eos_vector.cpp index b82eec4760..14106e940b 100644 --- a/test/test_eos_vector.cpp +++ b/test/test_eos_vector.cpp @@ -22,7 +22,8 @@ #include #include -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/test_eos_vinet.cpp b/test/test_eos_vinet.cpp index 926159b3c4..ee2e23ed43 100644 --- a/test/test_eos_vinet.cpp +++ b/test/test_eos_vinet.cpp @@ -16,8 +16,8 @@ #include #include #include -#include -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif diff --git a/test/test_math_utils.cpp b/test/test_math_utils.cpp index 1e7c15c7bb..301c34da7f 100644 --- a/test/test_math_utils.cpp +++ b/test/test_math_utils.cpp @@ -21,7 +21,8 @@ #include #include -#ifndef CATCH_CONFIG_RUNNER +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE #include "catch2/catch.hpp" #endif From fd42d5a330352834473d77c2b4e9995197846f40 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 11 Oct 2023 20:19:58 -0600 Subject: [PATCH 3/7] remove some unneeded headers --- test/test_eos_vector.cpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/test/test_eos_vector.cpp b/test/test_eos_vector.cpp index 14106e940b..7dcb1cff65 100644 --- a/test/test_eos_vector.cpp +++ b/test/test_eos_vector.cpp @@ -17,10 +17,7 @@ #include #include #include -#include -#include #include -#include #ifndef CATCH_CONFIG_FAST_COMPILE #define CATCH_CONFIG_FAST_COMPILE @@ -32,9 +29,6 @@ using singularity::EOS; using singularity::IdealGas; -namespace EOSBuilder = singularity::EOSBuilder; -namespace thermalqs = singularity::thermalqs; - SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { GIVEN("Parameters for an ideal gas") { From 58a4f9834c555656db1e634bcbce29478989386e Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 12 Oct 2023 09:58:21 -0600 Subject: [PATCH 4/7] shrink the variant in all the translation units --- test/eos_unit_test_helpers.hpp | 4 +--- test/profile_eos.cpp | 10 ++++------ test/test_eos_gruneisen.cpp | 2 +- test/test_eos_ideal.cpp | 2 +- test/test_eos_noble_abel.cpp | 5 +++-- test/test_eos_sap_polynomial.cpp | 2 +- test/test_eos_stiff.cpp | 5 +++-- test/test_eos_vector.cpp | 2 +- test/test_eos_vinet.cpp | 2 +- test/test_math_utils.cpp | 1 + 10 files changed, 17 insertions(+), 18 deletions(-) diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index a7c40d214a..a50206321a 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -52,8 +52,6 @@ inline std::string demangle(const char *name) { return name; } #include #include -using singularity::EOS; - PORTABLE_INLINE_FUNCTION bool isClose(Real a, Real b, Real eps = 5e-2) { return fabs(b - a) / (fabs(a + b) + 1e-20) <= eps; } @@ -74,7 +72,7 @@ inline void array_compare(int num, X &&x, Y &&y, Z &&z, ZT &&ztrue, XN xname, YN } template -inline void compare_two_eoss(const EOS &&test_e, const EOS &&ref_e) { +inline void compare_two_eoss(const E1 &&test_e, const E2 &&ref_e) { // compare all individual member functions with 1 as inputs, // this function is meant to catch mis-implementations of // modifiers that can be initialized in such a way as to diff --git a/test/profile_eos.cpp b/test/profile_eos.cpp index 1a731c6c1e..9eae5fa62b 100644 --- a/test/profile_eos.cpp +++ b/test/profile_eos.cpp @@ -33,7 +33,6 @@ #include #include -#include #include using namespace singularity; @@ -67,6 +66,7 @@ inline double get_duration(Function function) { TODO(JMM): Only profiles Pressure call and does not accept lambdas. */ +template inline void get_timing(int ncycles, const ivec &ncells_1d, const double rho_min, const double rho_max, const double T_min, const double T_max, const std::string &name, EOS &eos_h) { @@ -220,11 +220,9 @@ int main(int argc, char *argv[]) { << "Please note that timings are for 1 node and 1 GPU respectively.\n" << "Beginning profiling..." << std::endl; - EOSBuilder::EOSType type = EOSBuilder::EOSType::IdealGas; - EOSBuilder::params_t params; - params["Cv"].emplace(1e7 * 0.716); // specific heat in ergs/(g K) - params["gm1"].emplace(0.4); // gamma - 1 - EOS eos = EOSBuilder::buildEOS(type, params); + constexpr Real Cv = 1e7 * 0.716; // specific heat in ergs/(g K) + constexpr Real gm1 = 0.4; // gamma - 1 + auto eos = singularity::IdealGas(gm1, Cv); get_timing(ncycles, ncells_1d, RHO_MIN, RHO_MAX, T_MIN, T_MAX, "IdealGas", eos); std::cout << "Done." << std::endl; diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index 04936f106b..fcefd1fe3a 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -26,8 +26,8 @@ #include #include -using singularity::EOS; using singularity::Gruneisen; +using EOS = singularity::Variant; PORTABLE_INLINE_FUNCTION Real QuadFormulaMinus(Real a, Real b, Real c) { return (-b - std::sqrt(b * b - 4 * a * c)) / (2 * a); diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp index 37a1818772..72f769f894 100644 --- a/test/test_eos_ideal.cpp +++ b/test/test_eos_ideal.cpp @@ -33,8 +33,8 @@ #include -using singularity::EOS; using singularity::IdealGas; +using EOS = singularity::Variant; SCENARIO("Ideal gas entropy", "[IdealGas][Entropy]") { GIVEN("Parameters for an ideal gas with entropy reference states") { diff --git a/test/test_eos_noble_abel.cpp b/test/test_eos_noble_abel.cpp index e2a00e8414..b97dc833d7 100644 --- a/test/test_eos_noble_abel.cpp +++ b/test/test_eos_noble_abel.cpp @@ -26,8 +26,9 @@ #include #include -using singularity::EOS; using singularity::NobleAbel; +using singularity::IdealGas; +using EOS = singularity::Variant; SCENARIO("NobleAbel1", "[NobleAbel][NobleAbel1]") { GIVEN("Parameters for a NobleAbel EOS") { @@ -361,7 +362,7 @@ SCENARIO("Recover Ideal Gas from NobleAbel", "[NobleAbel][NobleAbel4]") { // Create the EOS EOS host_eos = NobleAbel(gm1, Cv, bb, qq); EOS eos = host_eos.GetOnDevice(); - EOS ideal_eos = singularity::IdealGas(gm1, Cv); + EOS ideal_eos = IdealGas(gm1, Cv); GIVEN("Densities and energies") { constexpr int num = 1; #ifdef PORTABILITY_STRATEGY_KOKKOS diff --git a/test/test_eos_sap_polynomial.cpp b/test/test_eos_sap_polynomial.cpp index b916a138e9..a531e2817f 100644 --- a/test/test_eos_sap_polynomial.cpp +++ b/test/test_eos_sap_polynomial.cpp @@ -26,8 +26,8 @@ #include #include -using singularity::EOS; using singularity::SAP_Polynomial; +using EOS = singularity::Variant; SCENARIO("SAP_Polynomial EOS", "Check if eos returns expected values") { GIVEN("Parameters for a SAP_Polynomial EOS") { diff --git a/test/test_eos_stiff.cpp b/test/test_eos_stiff.cpp index 2c8bb319d9..4e184aa29f 100644 --- a/test/test_eos_stiff.cpp +++ b/test/test_eos_stiff.cpp @@ -25,8 +25,9 @@ #include #include -using singularity::EOS; using singularity::StiffGas; +using singularity::IdealGas; +using EOS = singularity::Variant; SCENARIO("StiffGas1", "[StiffGas][StiffGas1]") { GIVEN("Parameters for a StiffGas EOS") { @@ -362,7 +363,7 @@ SCENARIO("Recover Ideal Gas from Stiff Gas", "[StiffGas][StiffGas4]") { // Create the EOS EOS host_eos = StiffGas(gm1, Cv, Pinf, qq); EOS eos = host_eos.GetOnDevice(); - EOS ideal_eos = singularity::IdealGas(gm1, Cv); + EOS ideal_eos = IdealGas(gm1, Cv); GIVEN("Densities and energies") { constexpr int num = 1; #ifdef PORTABILITY_STRATEGY_KOKKOS diff --git a/test/test_eos_vector.cpp b/test/test_eos_vector.cpp index 7dcb1cff65..7f839bf643 100644 --- a/test/test_eos_vector.cpp +++ b/test/test_eos_vector.cpp @@ -26,8 +26,8 @@ #include -using singularity::EOS; using singularity::IdealGas; +using EOS = singularity::Variant; SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") { diff --git a/test/test_eos_vinet.cpp b/test/test_eos_vinet.cpp index ee2e23ed43..08db7e4669 100644 --- a/test/test_eos_vinet.cpp +++ b/test/test_eos_vinet.cpp @@ -28,8 +28,8 @@ #include #include -using singularity::EOS; using singularity::Vinet; +using EOS = singularity::Variant; SCENARIO("Vinet EOS rho sie", "[VectorEOS][VinetEOS]") { GIVEN("Parameters for a Vinet EOS") { diff --git a/test/test_math_utils.cpp b/test/test_math_utils.cpp index 301c34da7f..211cdb332f 100644 --- a/test/test_math_utils.cpp +++ b/test/test_math_utils.cpp @@ -29,6 +29,7 @@ #include SCENARIO("EOS Variant Type", "[Variant][EOS]") { + using singularity::EOS; // print out the eos type std::cout << demangle(typeid(EOS).name()) << std::endl; } From 42040c30167e30b54bd152bb04d4b647de641061 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 12 Oct 2023 10:18:12 -0600 Subject: [PATCH 5/7] formatting and changelog --- CHANGELOG.md | 1 + test/eos_unit_test_helpers.hpp | 2 +- test/profile_eos.cpp | 4 ++-- test/test_eos_modifiers.cpp | 3 +-- test/test_eos_noble_abel.cpp | 2 +- test/test_eos_stiff.cpp | 2 +- 6 files changed, 7 insertions(+), 7 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8aef7c8d8d..a5f6a5dccd 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -40,6 +40,7 @@ - [[PR177]](https://github.com/lanl/singularity-eos/pull/177) added EOSPAC vector functions ### Changed (changing behavior/API/variables/...) +- [[PR310]](https://github.com/lanl/singularity-eos/pull/310) Speed up and clean up tests - [[PR295]](https://github.com/lanl/singularity-eos/pull/295) Add fast logs to singularity-eos - [[PR246]](https://github.com/lanl/singularity-eos/pull/246) Update CMake with upstream changes - [[PR223]](https://github.com/lanl/singularity-eos/pull/223) Update ports-of-call and add portable error handling diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index a50206321a..a5616f83eb 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -71,7 +71,7 @@ inline void array_compare(int num, X &&x, Y &&y, Z &&z, ZT &&ztrue, XN xname, YN } } -template +template inline void compare_two_eoss(const E1 &&test_e, const E2 &&ref_e) { // compare all individual member functions with 1 as inputs, // this function is meant to catch mis-implementations of diff --git a/test/profile_eos.cpp b/test/profile_eos.cpp index 9eae5fa62b..dfdff7d6d1 100644 --- a/test/profile_eos.cpp +++ b/test/profile_eos.cpp @@ -66,7 +66,7 @@ inline double get_duration(Function function) { TODO(JMM): Only profiles Pressure call and does not accept lambdas. */ -template +template inline void get_timing(int ncycles, const ivec &ncells_1d, const double rho_min, const double rho_max, const double T_min, const double T_max, const std::string &name, EOS &eos_h) { @@ -221,7 +221,7 @@ int main(int argc, char *argv[]) { << "Beginning profiling..." << std::endl; constexpr Real Cv = 1e7 * 0.716; // specific heat in ergs/(g K) - constexpr Real gm1 = 0.4; // gamma - 1 + constexpr Real gm1 = 0.4; // gamma - 1 auto eos = singularity::IdealGas(gm1, Cv); get_timing(ncycles, ncells_1d, RHO_MIN, RHO_MAX, T_MIN, T_MAX, "IdealGas", eos); diff --git a/test/test_eos_modifiers.cpp b/test/test_eos_modifiers.cpp index 74674f668d..2c6f9f3d29 100644 --- a/test/test_eos_modifiers.cpp +++ b/test/test_eos_modifiers.cpp @@ -27,11 +27,11 @@ #include +using singularity::BilinearRampEOS; using singularity::EOS; using singularity::IdealGas; using singularity::ScaledEOS; using singularity::ShiftedEOS; -using singularity::BilinearRampEOS; namespace EOSBuilder = singularity::EOSBuilder; namespace thermalqs = singularity::thermalqs; @@ -274,4 +274,3 @@ SCENARIO("EOS Unit System", "[EOSBuilder][UnitSystem][IdealGas]") { } } } - diff --git a/test/test_eos_noble_abel.cpp b/test/test_eos_noble_abel.cpp index b97dc833d7..c59fee244b 100644 --- a/test/test_eos_noble_abel.cpp +++ b/test/test_eos_noble_abel.cpp @@ -26,8 +26,8 @@ #include #include -using singularity::NobleAbel; using singularity::IdealGas; +using singularity::NobleAbel; using EOS = singularity::Variant; SCENARIO("NobleAbel1", "[NobleAbel][NobleAbel1]") { diff --git a/test/test_eos_stiff.cpp b/test/test_eos_stiff.cpp index 4e184aa29f..6871f0d406 100644 --- a/test/test_eos_stiff.cpp +++ b/test/test_eos_stiff.cpp @@ -25,8 +25,8 @@ #include #include -using singularity::StiffGas; using singularity::IdealGas; +using singularity::StiffGas; using EOS = singularity::Variant; SCENARIO("StiffGas1", "[StiffGas][StiffGas1]") { From 81df6ea57792f9c511917c837039432f59487db7 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 12 Oct 2023 10:24:16 -0600 Subject: [PATCH 6/7] update docs to reflect changes to test infrastructure --- doc/sphinx/src/contributing.rst | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/doc/sphinx/src/contributing.rst b/doc/sphinx/src/contributing.rst index a3536910e2..138a381935 100644 --- a/doc/sphinx/src/contributing.rst +++ b/doc/sphinx/src/contributing.rst @@ -233,10 +233,9 @@ test in the ``CMakeLists.txt`` file, .. code-block:: cmake - add_executable(eos_unit_tests + add_executable(eos_analytic_unit_tests catch2_define.cpp eos_unit_test_helpers.hpp - test_eos_unit.cpp test_eos_gruneisen.cpp test_eos_vinet.cpp test_my_new_eos.cpp @@ -245,15 +244,23 @@ test in the ``CMakeLists.txt`` file, in order for the test to be compiled. If your EOS requires any special dependencies, be sure to block off the test using ``#IFDEF`` blocks. -**Important:** this is a subtlety that highlights the importance of unit tests! -Since our library is header only, the unit tests are often the only place where -a specific EOS may be instantiated when ``singularity-eos`` is compiled. Unit -tests _must_ make use of the ``EOS`` type, i.e. +.. note:: + + Note that there are three executables, ``eos_analytic_unit_tests``, + ``eos_infrastructure_tests`` and ``eos_tabulated_unit_tests``. Pick + the executable that most closely matches what your model is. + +**This is important.** Since our library is header only, the unit +tests are often the only place where a specific EOS may be +instantiated when ``singularity-eos`` is compiled. Therefore to +exercise all code paths, it is best to create an ``EOS`` type +instantiated as .. code-block:: c++ #include - EOS my_eos = my_new_eos(parameter1, parameter2, ...) + using EOS = singularity::Variant;``. + EOS my_eos = MyNewEOS(parameter1, parameter2, ...) in order to properly test the functionality of a new EOS. Simply using the new class as the type such as From bd4bc244f5a02dc1c768fcfb62c0b426ecee8ee4 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 12 Oct 2023 10:42:06 -0600 Subject: [PATCH 7/7] contributing --- doc/sphinx/src/contributing.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/sphinx/src/contributing.rst b/doc/sphinx/src/contributing.rst index 138a381935..9a53e73888 100644 --- a/doc/sphinx/src/contributing.rst +++ b/doc/sphinx/src/contributing.rst @@ -250,7 +250,7 @@ dependencies, be sure to block off the test using ``#IFDEF`` blocks. ``eos_infrastructure_tests`` and ``eos_tabulated_unit_tests``. Pick the executable that most closely matches what your model is. -**This is important.** Since our library is header only, the unit +**Important:** Since our library is header only, the unit tests are often the only place where a specific EOS may be instantiated when ``singularity-eos`` is compiled. Therefore to exercise all code paths, it is best to create an ``EOS`` type