From ade7a101dcfc329e4db2e77755b8cd0e252e8b46 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Wed, 11 Oct 2023 19:55:06 -0600 Subject: [PATCH 01/20] Enable fortran preprocessor for fortran interface and then add same logic in fortran interface as is in the c/c++ interface. --- singularity-eos/CMakeLists.txt | 2 ++ singularity-eos/eos/singularity_eos.f90 | 26 +++++++++++++++++-------- 2 files changed, 20 insertions(+), 8 deletions(-) diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 12bfcff8193..d798b83c3ee 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -69,6 +69,8 @@ endif() if (SINGULARITY_USE_FORTRAN) list(APPEND EOS_SRCS eos/singularity_eos.f90) + set_source_files_properties(eos/singularity_eos.f90 + PROPERTIES Fortran_PREPROCESS ON) list(APPEND EOS_SRCS eos/singularity_eos.cpp) list(APPEND EOS_HEADERS eos/singularity_eos.hpp) # would rather handle this more robustly, being sloppy for now diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index 7f908f8880c..5a482db56bb 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -40,10 +40,16 @@ module singularity_eos init_sg_NobleAbel_f,& init_sg_SAP_Polynomial_f,& init_sg_StiffGas_f,& +#ifdef SINGULARITY_USE_SPINER_WITH_HDF5 +#ifdef SINGULARITY_USE_HELMHOLTZ init_sg_Helmholtz_f,& +#endif ! SINGULARITY_USE_HELMHOLTZ init_sg_SpinerDependsRhoT_f,& init_sg_SpinerDependsRhoSie_f,& +#endif ! SINGULARITY_USE_SPINER_WITH_HDF5 +#ifdef SINGULARITY_USE_EOSPAC init_sg_eospac_f,& +#endif ! SINGULARITY_USE_EOSPAC get_sg_eos_f,& finalize_sg_eos_f @@ -164,7 +170,8 @@ end function init_sg_StiffGas type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values end function init_sg_SAP_Polynomial end interface - +#ifdef SINGULARITY_USE_SPINER_WITH_HDF5 +#ifdef SINGULARITY_USE_HELMHOLTZ interface integer(kind=c_int) function & init_sg_Helmholtz(matindex, eos, filename, rad, gas, coul, ion, ele, & @@ -179,7 +186,7 @@ end function init_sg_SAP_Polynomial type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values end function init_sg_Helmholtz end interface - +#endif ! SINGULARITY_USE_HELMHOLTZ interface integer(kind=c_int) function & init_sg_SpinerDependsRhoT(matindex, eos, filename, id, sg_mods_enabled, & @@ -205,7 +212,8 @@ end function init_sg_SpinerDependsRhoT type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values end function init_sg_SpinerDependsRhoSie end interface - +#endif ! SINGULARITY_USE_SPINER_WITH_HDF5 +#ifdef SINGULARITY_USE_EOSPAC interface integer(kind=c_int) function & init_sg_eospac(matindex, eos, id, sg_mods_enabled, sg_mods_values) & @@ -216,7 +224,7 @@ end function init_sg_SpinerDependsRhoSie type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values end function init_sg_eospac end interface - +#endif ! SINGULARITY_USE_EOSPAC interface integer(kind=c_int) function & get_sg_eos(nmat, ncell, cell_dim,& @@ -443,7 +451,8 @@ integer function init_sg_NobleAbel_f(matindex, eos, gm1, Cv, & err = init_sg_NobleAbel(matindex-1, eos%ptr, gm1, Cv, bb, qq, & c_loc(sg_mods_enabled), c_loc(sg_mods_values)) end function init_sg_NobleAbel_f - +#ifdef SINGULARITY_USE_SPINER_WITH_HDF5 +#ifdef SINGULARITY_USE_HELMHOLTZ integer function init_sg_Helmholtz_f(matindex, eos, filename, rad, gas, coul, ion, ele, & verbose, sg_mods_enabled, sg_mods_values) & result(err) @@ -458,7 +467,7 @@ integer function init_sg_Helmholtz_f(matindex, eos, filename, rad, gas, coul, io rad, gas, coul, ion, ele, verbose, & c_loc(sg_mods_enabled), c_loc(sg_mods_values)) end function init_sg_Helmholtz_f - +#endif ! SINGULARITY_USE_HELMHOLTZ integer function init_sg_SpinerDependsRhoT_f(matindex, eos, filename, id, & sg_mods_enabled, & sg_mods_values) & @@ -489,7 +498,8 @@ integer function init_sg_SpinerDependsRhoSie_f(matindex, eos, filename, id, & c_loc(sg_mods_enabled), & c_loc(sg_mods_values)) end function init_sg_SpinerDependsRhoSie_f - +#endif ! SINGULARITY_USE_SPINER_WITH_HDF5 +#ifdef SINGULARITY_USE_EOSPAC integer function init_sg_eospac_f(matindex, eos, id, sg_mods_enabled, & sg_mods_values) & result(err) @@ -500,7 +510,7 @@ integer function init_sg_eospac_f(matindex, eos, id, sg_mods_enabled, & err = init_sg_eospac(matindex-1, eos%ptr, id, c_loc(sg_mods_enabled), & c_loc(sg_mods_values)) end function init_sg_eospac_f - +#endif ! SINGULARITY_USE_EOSPAC integer function finalize_sg_eos_f(nmat, eos) & result(err) integer(c_int), value, intent(in) :: nmat From 7c9edccd3765bdeba45f091e72c1f67f0af67554 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Thu, 12 Oct 2023 09:16:39 -0600 Subject: [PATCH 02/20] Actually enabled fortran preprocessor. --- CMakeLists.txt | 1 + singularity-eos/CMakeLists.txt | 2 -- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 59ad2c2c3d1..f294330d614 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -390,6 +390,7 @@ if(SINGULARITY_USE_FORTRAN) target_include_directories( singularity-eos INTERFACE $ $) + set_target_properties(singularity-eos PROPERTIES Fortran_PREPROCESS ON) endif() # SINGULARITY_USE_FORTRAN target_include_directories( diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index d798b83c3ee..12bfcff8193 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -69,8 +69,6 @@ endif() if (SINGULARITY_USE_FORTRAN) list(APPEND EOS_SRCS eos/singularity_eos.f90) - set_source_files_properties(eos/singularity_eos.f90 - PROPERTIES Fortran_PREPROCESS ON) list(APPEND EOS_SRCS eos/singularity_eos.cpp) list(APPEND EOS_HEADERS eos/singularity_eos.hpp) # would rather handle this more robustly, being sloppy for now From 1f0e6a01cc14cb326f0c71de8109515d45f8036e Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Wed, 18 Oct 2023 17:43:39 -0600 Subject: [PATCH 03/20] Fix of intel compiler error. Thanks Richard for the fixgit add singularity-eos/eos/modifiers/ramps_eos.hpp --- singularity-eos/eos/modifiers/ramps_eos.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index a0111b8ea5c..935c32a6dff 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -75,6 +75,9 @@ void pAlpha2BilinearRampParams(const T &eos, const Real alpha0, const Real Pe, template class BilinearRampEOS : public EosBase> { public: + // Vector functions that overload the scalar versions declared here. + SG_ADD_BASE_CLASS_USINGS(BilinearRampEOS) + // move semantics ensures dynamic memory comes along for the ride BilinearRampEOS(T &&t, const Real r0, const Real a, const Real b, const Real c) : t_(std::forward(t)), r0_(r0), a_(a), b_(b), c_(c), @@ -425,9 +428,6 @@ class BilinearRampEOS : public EosBase> { t_.ValuesAtReferenceState(rho, temp, sie, press, cv, bmod, dpde, dvdt, lambda); } - // Vector functions that overload the scalar versions declared here. - SG_ADD_BASE_CLASS_USINGS(BilinearRampEOS) - inline constexpr bool IsModified() const { return true; } inline constexpr T UnmodifyOnce() { return t_; } From c0e11fdf38c7c00cb04c3a160012cfb1b932e3b2 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Wed, 18 Oct 2023 17:44:12 -0600 Subject: [PATCH 04/20] Add more precision to catch info. --- test/eos_unit_test_helpers.hpp | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index edd3f5279bc..ba0b8f82121 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -22,6 +22,10 @@ #include #include +#include +#include +#include + // typename demangler #ifdef __GNUG__ #include @@ -64,9 +68,16 @@ PORTABLE_INLINE_FUNCTION Real myAtan(Real x, Real shift, Real scale, Real offset template inline void array_compare(int num, X &&x, Y &&y, Z &&z, ZT &&ztrue, XN xname, YN yname, Real tol = 1e-12) { + assert(num > 0); + using underlying_t = typename std::remove_cv::type>::type; for (int i = 0; i < num; i++) { - INFO("i: " << i << ", " << xname << ": " << x[i] << ", " << yname << ": " << y[i] - << ", Value: " << z[i] << ", True Value: " << ztrue[i]); + auto s = std::ostringstream{}; + s << std::setprecision(std::numeric_limits::max_digits10) + << std::scientific + << "i: " << i << ", " << xname << ": " << x[i] << ", " + << yname << ": " << y[i] << ", Value: " << z[i] + << ", True Value: " << ztrue[i]; + INFO(s.str()); CHECK(isClose(z[i], ztrue[i], 1e-12)); } } From f50b4224e258f4a2e9e2e79e76eb60e2db0d03ed Mon Sep 17 00:00:00 2001 From: jdolence Date: Tue, 10 Oct 2023 08:00:37 -0600 Subject: [PATCH 05/20] fix PTE when downstream code defines its own variant --- singularity-eos/closure/mixed_cell_models.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index e9c2b8d2104..5fd7c10172e 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -272,8 +272,9 @@ class PTESolverBase { return Tguess; } + template PORTABLE_FORCEINLINE_FUNCTION - static Real GetPressureFromPreferred(const EOS &eos, const Real rho, const Real T, + static Real GetPressureFromPreferred(const EOS_t &eos, const Real rho, const Real T, Real sie, Real *lambda, const bool do_e_lookup) { Real P{}; if (eos.PreferredInput() == From 47790a504c9df822b70b3688e8f56c378e86395c Mon Sep 17 00:00:00 2001 From: jdolence Date: Tue, 10 Oct 2023 11:34:22 -0600 Subject: [PATCH 06/20] formatting --- singularity-eos/closure/mixed_cell_models.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 5fd7c10172e..365e1dd4a6c 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -273,9 +273,9 @@ class PTESolverBase { } template - PORTABLE_FORCEINLINE_FUNCTION - static Real GetPressureFromPreferred(const EOS_t &eos, const Real rho, const Real T, - Real sie, Real *lambda, const bool do_e_lookup) { + PORTABLE_FORCEINLINE_FUNCTION static Real + GetPressureFromPreferred(const EOS_t &eos, const Real rho, const Real T, Real sie, + Real *lambda, const bool do_e_lookup) { Real P{}; if (eos.PreferredInput() == (thermalqs::density | thermalqs::specific_internal_energy)) { From ff608ee5332f637f30b3aef4eb364aa922f3fc74 Mon Sep 17 00:00:00 2001 From: jdolence Date: Tue, 10 Oct 2023 12:53:20 -0600 Subject: [PATCH 07/20] make get a host device function --- singularity-eos/eos/eos_variant.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index 3e72d87662f..09b2e4b9178 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -62,7 +62,7 @@ class Variant { typename std::enable_if< !std::is_same::type>::value, bool>::type = true> - EOSChoice get() { + PORTABLE_INLINE_FUNCTION EOSChoice get() { return mpark::get(eos_); } From fbf8cd82ad4d131b15dc09a64ddcfc10a5481011 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 11 Oct 2023 19:15:22 -0600 Subject: [PATCH 08/20] 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 db952eda522..5f71e9a5a48 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 ba0b8f82121..d7a74dc98f0 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -53,6 +53,7 @@ inline std::string demangle(const char *name) { return name; } #endif #include +#include using singularity::EOS; @@ -82,34 +83,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 00000000000..117c61e9142 --- /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 00000000000..21400018306 --- /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 00000000000..8979554fabe --- /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 00000000000..bce04d0d954 --- /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 04f5ff009c9..00000000000 --- 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 00000000000..b82eec4760e --- /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 00000000000..1e7c15c7bba --- /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 d1ba3cb3dba8ee2982237e25f6bbcd36b56502f6 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 11 Oct 2023 20:16:34 -0600 Subject: [PATCH 09/20] 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 5f71e9a5a48..d30b5b0b8e0 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 65437e8f74f..c06cf7a5c0d 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 d7a74dc98f0..90747589bda 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 c91ea9bd37e..04936f106b8 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 b1048ee5c87..319463cc61e 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 117c61e9142..37a18187724 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 21400018306..74674f668d5 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 5d6b4264bd4..e2a00e84142 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 d88ec89b10b..b916a138e9f 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 8979554fabe..859031d5b59 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 548af96a939..2c8bb319d99 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 bce04d0d954..dfaea1a05d5 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 b82eec4760e..14106e940bb 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 926159b3c4e..ee2e23ed43c 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 1e7c15c7bba..301c34da7f2 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 9e47c88ce94d5859a66ca2a249d80beb4e5c80e8 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Wed, 11 Oct 2023 20:19:58 -0600 Subject: [PATCH 10/20] 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 14106e940bb..7dcb1cff655 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 9f0dabfde8bf2df632c0b80d922e1aa7fec46f1e Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 12 Oct 2023 09:58:21 -0600 Subject: [PATCH 11/20] 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 90747589bda..9d55b4c30ff 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -56,8 +56,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; } @@ -85,7 +83,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 1a731c6c1ea..9eae5fa62b8 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 04936f106b8..fcefd1fe3ae 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 37a18187724..72f769f8946 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 e2a00e84142..b97dc833d73 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 b916a138e9f..a531e2817f0 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 2c8bb319d99..4e184aa29f4 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 7dcb1cff655..7f839bf643f 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 ee2e23ed43c..08db7e4669e 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 301c34da7f2..211cdb332fe 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 b77437d69c6b7dc5c0ee08718a02a2bc643a47b0 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 12 Oct 2023 10:18:12 -0600 Subject: [PATCH 12/20] 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 8aef7c8d8d2..a5f6a5dccd0 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 9d55b4c30ff..e7f36813363 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -82,7 +82,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 9eae5fa62b8..dfdff7d6d10 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 74674f668d5..2c6f9f3d296 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 b97dc833d73..c59fee244b9 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 4e184aa29f4..6871f0d4069 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 aea1e951bdb92e4381aead368d6af0ea523dc3fa Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 12 Oct 2023 10:24:16 -0600 Subject: [PATCH 13/20] 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 a3536910e27..138a3819354 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 d0df5ab1e8931ca6e7396724d757e0a10f7684aa Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 12 Oct 2023 10:42:06 -0600 Subject: [PATCH 14/20] 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 138a3819354..9a53e738881 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 From 429cc0f7cff8aad6ff8dbc890d91e6c9c37ee9f2 Mon Sep 17 00:00:00 2001 From: Peter Brady Date: Mon, 16 Oct 2023 16:53:26 -0500 Subject: [PATCH 15/20] Add eos_sap_polynomial to list of headers for install --- singularity-eos/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 12bfcff8193..6d6daa7d284 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -33,6 +33,7 @@ set(EOS_HEADERS eos/eos_builder.hpp eos/eos_jwl.hpp eos/eos_helmholtz.hpp + eos/eos_sap_polynomial.hpp eos/modifiers/relativistic_eos.hpp eos/modifiers/scaled_eos.hpp eos/modifiers/ramps_eos.hpp From 479095827b27255954c24b1ea77b4dca8a08767e Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Mon, 23 Oct 2023 16:46:35 -0600 Subject: [PATCH 16/20] The golds were not calculated to enough precision. They have been updated. The array compare now showns all digits on printout. Change == to isclose, intel tests were failing otherwise. Also format. --- test/eos_unit_test_helpers.hpp | 13 ++++++------- test/test_eos_gruneisen.cpp | 2 +- test/test_eos_stellar_collapse.cpp | 2 +- test/test_eos_stiff.cpp | 12 ++++++------ 4 files changed, 14 insertions(+), 15 deletions(-) diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index e7f36813363..b23f5cf93e3 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -23,9 +23,9 @@ #include #include +#include #include #include -#include // typename demangler #ifdef __GNUG__ @@ -69,14 +69,13 @@ template 0); - using underlying_t = typename std::remove_cv::type>::type; + using underlying_t = + typename std::remove_cv::type>::type; for (int i = 0; i < num; i++) { auto s = std::ostringstream{}; - s << std::setprecision(std::numeric_limits::max_digits10) - << std::scientific - << "i: " << i << ", " << xname << ": " << x[i] << ", " - << yname << ": " << y[i] << ", Value: " << z[i] - << ", True Value: " << ztrue[i]; + s << std::setprecision(std::numeric_limits::max_digits10) + << std::scientific << "i: " << i << ", " << xname << ": " << x[i] << ", " << yname + << ": " << y[i] << ", Value: " << z[i] << ", True Value: " << ztrue[i]; INFO(s.str()); CHECK(isClose(z[i], ztrue[i], 1e-12)); } diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index fcefd1fe3ae..e0c0b9bbfc3 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -479,7 +479,7 @@ SCENARIO("Gruneisen EOS density limit") { INFO("FillEos bmod: " << bmod << ", Lookup bmod: " << bmod_true); CHECK(bmod == bmod_true); INFO("FillEos pressure: " << P << ", Lookup pressure: " << pres_true); - CHECK(P == pres_true); + CHECK(isClose(P, pres_true, 1.e-14)); } } } diff --git a/test/test_eos_stellar_collapse.cpp b/test/test_eos_stellar_collapse.cpp index 859031d5b59..dc067f43582 100644 --- a/test/test_eos_stellar_collapse.cpp +++ b/test/test_eos_stellar_collapse.cpp @@ -44,7 +44,7 @@ SCENARIO("Test 3D reinterpolation to fast log grid", "[StellarCollapse]") { constexpr int N2 = 100; constexpr int N1 = 101; constexpr int N0 = 102; - StellarCollapse::Grid_t g2(0, 1, N2); + StellarCollapse::Grid_t g2(1.0/N2, 1, N2); StellarCollapse::Grid_t g1(1, 3, N1); StellarCollapse::Grid_t g0(2, 4, N0); StellarCollapse::DataBox db(N2, N1, N0); diff --git a/test/test_eos_stiff.cpp b/test/test_eos_stiff.cpp index 6871f0d4069..0ab248a2367 100644 --- a/test/test_eos_stiff.cpp +++ b/test/test_eos_stiff.cpp @@ -79,14 +79,14 @@ SCENARIO("StiffGas1", "[StiffGas][StiffGas1]") { // Gold standard values for a subset of lookups constexpr std::array pressure_true{ - 1.0132500000019073e+06, 3.4450500000000000e+07, 6.7887750000001907e+07, - 1.0132500000000000e+08}; + 1.01324999999964016e+06, 3.44504999999983162e+07, 6.78877500000010729e+07, + 1.01324999999999583e+08}; constexpr std::array bulkmodulus_true{ - 2.3502381137500000e+10, 2.3580958675000000e+10, 2.3659536212500004e+10, - 2.3738113750000004e+10}; + 2.35023811375000000e+10, 2.35809586749999962e+10, 2.36595362125000038e+10, + 2.37381137500000000e+10}; constexpr std::array temperature_true{ - 2.9814999999999998e+02, 1.5320999999999999e+03, 2.7660500000000002e+03, - 4.0000000000000000e+03}; + 2.98149999999999920e+02, 1.53209999999999968e+03, 2.76605000000000064e+03, + 3.99999999999999955e+03}; constexpr std::array gruneisen_true{1.35, 1.35, 1.35, 1.35}; #ifdef PORTABILITY_STRATEGY_KOKKOS From 653cc1423425999db3e915140559e082b5e9f69d Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Mon, 23 Oct 2023 17:12:39 -0600 Subject: [PATCH 17/20] format. --- test/test_eos_stellar_collapse.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_eos_stellar_collapse.cpp b/test/test_eos_stellar_collapse.cpp index dc067f43582..66366fe186e 100644 --- a/test/test_eos_stellar_collapse.cpp +++ b/test/test_eos_stellar_collapse.cpp @@ -44,7 +44,7 @@ SCENARIO("Test 3D reinterpolation to fast log grid", "[StellarCollapse]") { constexpr int N2 = 100; constexpr int N1 = 101; constexpr int N0 = 102; - StellarCollapse::Grid_t g2(1.0/N2, 1, N2); + StellarCollapse::Grid_t g2(1.0 / N2, 1, N2); StellarCollapse::Grid_t g1(1, 3, N1); StellarCollapse::Grid_t g0(2, 4, N0); StellarCollapse::DataBox db(N2, N1, N0); From 67327b7f7cde44daffd273917f6c1516eec609f2 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 31 Oct 2023 09:52:35 -0600 Subject: [PATCH 18/20] Modify a currently unbuilt fortran interface test file to test the fortran interface and introduce it to cmake. Make the fortran interface have optional modifier parameters. --- singularity-eos/eos/singularity_eos.f90 | 192 ++++++++++++++++++------ test/CMakeLists.txt | 7 + test/test_f_iface.f90 | 85 +++-------- 3 files changed, 174 insertions(+), 110 deletions(-) diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index 5a482db56bb..96a605bdda2 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -349,10 +349,19 @@ integer function init_sg_IdealGas_f(matindex, eos, gm1, Cv, & integer(c_int), value, intent(in) :: matindex type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: gm1, Cv - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + err = init_sg_IdealGas(matindex-1, eos%ptr, gm1, Cv, & - c_loc(sg_mods_enabled), c_loc(sg_mods_values)) + c_loc(sg_mods_enabled_use), c_loc(sg_mods_values_use)) end function init_sg_IdealGas_f integer function init_sg_Gruneisen_f(matindex, eos, C0, s1, s2, s3, G0, b,& @@ -363,11 +372,20 @@ integer function init_sg_Gruneisen_f(matindex, eos, C0, s1, s2, s3, G0, b,& type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: C0, s1, s2, s3, G0, b, rho0 real(kind=8), value, intent(in) :: T0, P0, Cv - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + err = init_sg_Gruneisen(matindex-1, eos%ptr, C0, s1, s2, s3, G0, b, rho0,& - T0, P0, Cv, c_loc(sg_mods_enabled), & - c_loc(sg_mods_values)) + T0, P0, Cv, c_loc(sg_mods_enabled_use), & + c_loc(sg_mods_values_use)) end function init_sg_Gruneisen_f integer function init_sg_JWL_f(matindex, eos, A, B, R1, R2, w, rho0, Cv, & @@ -376,10 +394,19 @@ integer function init_sg_JWL_f(matindex, eos, A, B, R1, R2, w, rho0, Cv, & integer(c_int), value, intent(in) :: matindex type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: A, B, R1, R2, w, rho0, Cv - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + err = init_sg_JWL(matindex-1, eos%ptr, A, B, R1, R2, w, rho0, Cv, & - c_loc(sg_mods_enabled), c_loc(sg_mods_values)) + c_loc(sg_mods_enabled_use), c_loc(sg_mods_values_use)) end function init_sg_JWL_f integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, & @@ -389,11 +416,20 @@ integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, & integer(c_int), value, intent(in) :: matindex type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: a, b, k, n, vc, pc, Cv, E0 - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + err = init_sg_DavisProducts(matindex-1, eos%ptr, a, b, k, n, vc, pc, Cv, & - E0, c_loc(sg_mods_enabled), & - c_loc(sg_mods_values)) + E0, c_loc(sg_mods_enabled_use), & + c_loc(sg_mods_values_use)) end function init_sg_DavisProducts_f integer function init_sg_DavisReactants_f(matindex, eos, rho0, e0, P0, T0, & @@ -404,11 +440,20 @@ integer function init_sg_DavisReactants_f(matindex, eos, rho0, e0, P0, T0, & type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: rho0, e0, P0, T0, A, B, C, G0, Z real(kind=8), value, intent(in) :: alpha, Cv0 - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + err = init_sg_DavisReactants(matindex-1, eos%ptr, rho0, e0, P0, T0, A, B, & - C, G0, Z, alpha, Cv0, c_loc(sg_mods_enabled), & - c_loc(sg_mods_values)) + C, G0, Z, alpha, Cv0, c_loc(sg_mods_enabled_use), & + c_loc(sg_mods_values_use)) end function init_sg_DavisReactants_f integer function init_sg_SAP_Polynomial_f(matindex, eos, rho0, a0, a1, a2c, & @@ -419,11 +464,20 @@ integer function init_sg_SAP_Polynomial_f(matindex, eos, rho0, a0, a1, a2c, & type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: rho0, a0, a1, a2c, a2e, a3,& b0, b1, b2c, b2e, b3 - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + err = init_sg_SAP_Polynomial(matindex-1, eos%ptr, rho0, a0, a1, a2c, a2e, & a3, b0, b1, b2c, b2e, b3, & - c_loc(sg_mods_enabled), c_loc(sg_mods_values)) + c_loc(sg_mods_enabled_use), c_loc(sg_mods_values_use)) end function init_sg_SAP_Polynomial_f integer function init_sg_StiffGas_f(matindex, eos, gm1, Cv, & @@ -433,10 +487,19 @@ integer function init_sg_StiffGas_f(matindex, eos, gm1, Cv, & integer(c_int), value, intent(in) :: matindex type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: gm1, Cv, Pinf, qq - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + err = init_sg_StiffGas(matindex-1, eos%ptr, gm1, Cv, Pinf, qq, & - c_loc(sg_mods_enabled), c_loc(sg_mods_values)) + c_loc(sg_mods_enabled_use), c_loc(sg_mods_values_use)) end function init_sg_StiffGas_f integer function init_sg_NobleAbel_f(matindex, eos, gm1, Cv, & @@ -446,10 +509,19 @@ integer function init_sg_NobleAbel_f(matindex, eos, gm1, Cv, & integer(c_int), value, intent(in) :: matindex type(sg_eos_ary_t), intent(in) :: eos real(kind=8), value, intent(in) :: gm1, Cv, bb, qq - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + err = init_sg_NobleAbel(matindex-1, eos%ptr, gm1, Cv, bb, qq, & - c_loc(sg_mods_enabled), c_loc(sg_mods_values)) + c_loc(sg_mods_enabled_use), c_loc(sg_mods_values_use)) end function init_sg_NobleAbel_f #ifdef SINGULARITY_USE_SPINER_WITH_HDF5 #ifdef SINGULARITY_USE_HELMHOLTZ @@ -461,11 +533,20 @@ integer function init_sg_Helmholtz_f(matindex, eos, filename, rad, gas, coul, io character(len=*, kind=c_char), intent(in) :: filename logical(c_bool), value, intent(in) :: rad, gas, coul, ion, ele, & verbose - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + err = init_sg_Helmholtz(matindex-1, eos%ptr, trim(filename)//C_NULL_CHAR, & rad, gas, coul, ion, ele, verbose, & - c_loc(sg_mods_enabled), c_loc(sg_mods_values)) + c_loc(sg_mods_enabled_use), c_loc(sg_mods_values_use)) end function init_sg_Helmholtz_f #endif ! SINGULARITY_USE_HELMHOLTZ integer function init_sg_SpinerDependsRhoT_f(matindex, eos, filename, id, & @@ -476,12 +557,21 @@ integer function init_sg_SpinerDependsRhoT_f(matindex, eos, filename, id, & type(sg_eos_ary_t), intent(in) :: eos character(len=*, kind=c_char), intent(in) :: filename integer(c_int), intent(inout) :: id - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + err = init_sg_SpinerDependsRhoT(matindex-1, eos%ptr,& trim(filename)//C_NULL_CHAR, id, & - c_loc(sg_mods_enabled), & - c_loc(sg_mods_values)) + c_loc(sg_mods_enabled_use), & + c_loc(sg_mods_values_use)) end function init_sg_SpinerDependsRhoT_f integer function init_sg_SpinerDependsRhoSie_f(matindex, eos, filename, id, & @@ -491,12 +581,21 @@ integer function init_sg_SpinerDependsRhoSie_f(matindex, eos, filename, id, & integer(c_int), value, intent(in) :: matindex, id type(sg_eos_ary_t), intent(in) :: eos character(len=*, kind=c_char), intent(in) :: filename - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + err = init_sg_SpinerDependsRhoSie(matindex-1, eos%ptr,& trim(filename)//C_NULL_CHAR, id, & - c_loc(sg_mods_enabled), & - c_loc(sg_mods_values)) + c_loc(sg_mods_enabled_use), & + c_loc(sg_mods_values_use)) end function init_sg_SpinerDependsRhoSie_f #endif ! SINGULARITY_USE_SPINER_WITH_HDF5 #ifdef SINGULARITY_USE_EOSPAC @@ -505,10 +604,19 @@ integer function init_sg_eospac_f(matindex, eos, id, sg_mods_enabled, & result(err) integer(c_int), value, intent(in) :: matindex, id type(sg_eos_ary_t), intent(in) :: eos - integer(kind=c_int), dimension(:), target, intent(inout) :: sg_mods_enabled - real(kind=8), dimension(:), target, intent(inout) :: sg_mods_values - err = init_sg_eospac(matindex-1, eos%ptr, id, c_loc(sg_mods_enabled), & - c_loc(sg_mods_values)) + integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled + real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values + ! local vars + integer(kind=c_int), dimension(4) :: sg_mods_enabled_use + real(kind=8), dimension(6) :: sg_mods_values_use + + sg_mods_enabled_use = 0 + sg_mods_values_use = 0.d0 + if(present(sg_mods_enabled)) sg_mods_enabled_use = sg_mods_enabled + if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values + + err = init_sg_eospac(matindex-1, eos%ptr, id, c_loc(sg_mods_enabled_use), & + c_loc(sg_mods_values_use)) end function init_sg_eospac_f #endif ! SINGULARITY_USE_EOSPAC integer function finalize_sg_eos_f(nmat, eos) & diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d30b5b0b8e0..80b10973295 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -131,6 +131,13 @@ if(SINGULARITY_USE_SPINER) target_link_libraries(profile_eos singularity-eos::singularity-eos) endif() +if(SINGULARITY_USE_FORTRAN) + add_executable(ftn_interface test_f_iface.f90) + target_link_libraries(ftn_interface singularity-eos::singularity-eos) + set_target_properties(ftn_interface PROPERTIES LINKER_LANGUAGE Fortran) + add_test(test_fortran_interface ftn_interface) +endif() + if(SINGULARITY_USE_EOSPAC AND SINGULARITY_TEST_SESAME AND NOT SINGULARITY_USE_CUDA) diff --git a/test/test_f_iface.f90 b/test/test_f_iface.f90 index 61ca03a6a83..b69bd7d33f5 100644 --- a/test/test_f_iface.f90 +++ b/test/test_f_iface.f90 @@ -12,81 +12,30 @@ ! publicly and display publicly, and to permit others to do so. !------------------------------------------------------------------------------ -program test -use bedroom_door_eos -integer :: res, nmat, ncell, i, j, matid -type(bd_eos_ary_t) :: eos -integer, allocatable, dimension(:) :: offsets -real(kind=8), allocatable, dimension(:) :: spvol, sie_tot, vsum, press, pmax -real(kind=8), allocatable, dimension(:) :: temp, bmod, dpde, cv -real(kind=8), allocatable, dimension(:,:) :: frac_mass, frac_vol, frac_sie -real(kind=8) :: vsumint -character(len=64) :: filename +program test_sg_fortran_interface +! modules +use singularity_eos_types +use singularity_eos +! no implicit vars +implicit none +! variable declaration +integer :: nmat, res +type(sg_eos_ary_t) :: eos -ncell = 1 +! set test parameters nmat = 3 -filename = "../data/materials.sp5" -! allocate arrays -allocate(spvol(ncell), sie_tot(ncell), temp(ncell), press(ncell), vsum(ncell)) -allocate(bmod(ncell), dpde(ncell), cv(ncell), pmax(ncell)) -allocate(offsets(ncell)) -allocate(frac_mass(ncell, nmat)) -allocate(frac_vol(ncell, nmat), frac_sie(ncell, nmat)) -! initialize eos's -res = init_bd_eos_f(nmat, eos) +! allocate and initialize eos's +res = init_sg_eos_f(nmat, eos) -res = init_Gruneisen_f(1, eos, 394000.d0, 1.489d0, 0.d0, 0.d0, 2.02d0, 0.47d0,& +res = init_sg_Gruneisen_f(1, eos, 394000.d0, 1.489d0, 0.d0, 0.d0, 2.02d0, 0.47d0,& 8.93d0, 297.0d0, 1.0d6, 0.383d7) -res = init_DavisReactants_f(2, eos, 1.890d0, 4.115d10, 1.0d6, 297.0d0, 1.8d0,& +res = init_sg_DavisReactants_f(2, eos, 1.890d0, 4.115d10, 1.0d6, 297.0d0, 1.8d0,& 4.6d0, 0.34d0, 0.56d0,0.d0, 0.4265d0, 0.001074d10) -res = init_DavisProducts_f(3, eos, 0.798311d0, 0.58d0, 1.35d0, 2.66182d0,& +res = init_sg_DavisProducts_f(3, eos, 0.798311d0, 0.58d0, 1.35d0, 2.66182d0,& 0.75419d0, 3.2d10, 0.001072d10, 0.d0) -!matid = 4272 -!res = init_SpinerDependsRhoT_f(2, eos, filename, matid) -!matid = 5030 -!res = init_SpinerDependsRhoT_f(3, eos, filename, matid) -! initialize state variables -call srand(10) -press = 0.d0 -do i=1,ncell - vsum(i) = 1.0d-2 - frac_mass(i,1) = 8.93d0*vsum(i)*0.33 - frac_mass(i,2) = 1.89d0*vsum(i)*.033 - frac_mass(i,3) = 2.5d0*vsum(i)*0.34 - sie_tot(i) = 0.d0 - spvol(i) = 0.d0 - !vsum(i) = 0.d0 - offsets(i) = i - frac_sie(i,1) = 1.16049000000000000d+09 - frac_sie(i,2) = 4.50111706015744858d+10 - frac_sie(i,3) = 3.29477034927098885d+10 - temp(i) = 0.d0 - ! frac_vols will be normalized in eos call - do j=1,nmat - frac_vol(i,j) = 0.d0 - enddo - do j=1,nmat - spvol(i) = spvol(i) + frac_mass(i,j) - sie_tot(i) = sie_tot(i) + frac_mass(i,j)*frac_sie(i,j) - frac_sie(i,j) = frac_sie(i,j)*frac_mass(i,j) ! needs to be energy - enddo - sie_tot(i) = sie_tot(i)/spvol(i) - spvol(i) = vsum(i)/spvol(i) -enddo -!write(*,*) 'spvol', spvol -write(*,*) frac_mass, frac_vol, frac_sie -! calculate PTE eos -res = get_bd_eos_f(nmat, ncell, ncell, 0, eos, offsets, press, pmax, vsum,& - spvol, sie_tot, temp, bmod, dpde, cv, frac_mass, frac_vol,& - frac_sie) - -write(*,*) press, cv, bmod, temp -write(*,*) frac_mass, frac_vol, frac_sie ! cleanup -res = finalize_bd_eos_f(nmat, eos) +res = finalize_sg_eos_f(nmat, eos) -deallocate(vsum, sie_tot, spvol, frac_mass, temp, press, frac_vol, frac_sie) -deallocate(bmod, dpde, cv, offsets, pmax) -end program test +end program test_sg_fortran_interface From 67a2ad6fbda04db006bed1d0a4fd9349c038d612 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 31 Oct 2023 10:02:28 -0600 Subject: [PATCH 19/20] fix gnu issues that did not appear with intel. --- singularity-eos/eos/singularity_eos.f90 | 66 ++++++++++++++----------- 1 file changed, 36 insertions(+), 30 deletions(-) diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index 96a605bdda2..d04b98d68f9 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -43,13 +43,16 @@ module singularity_eos #ifdef SINGULARITY_USE_SPINER_WITH_HDF5 #ifdef SINGULARITY_USE_HELMHOLTZ init_sg_Helmholtz_f,& -#endif ! SINGULARITY_USE_HELMHOLTZ +#endif +! SINGULARITY_USE_HELMHOLTZ init_sg_SpinerDependsRhoT_f,& init_sg_SpinerDependsRhoSie_f,& -#endif ! SINGULARITY_USE_SPINER_WITH_HDF5 +#endif +! SINGULARITY_USE_SPINER_WITH_HDF5 #ifdef SINGULARITY_USE_EOSPAC init_sg_eospac_f,& -#endif ! SINGULARITY_USE_EOSPAC +#endif +! SINGULARITY_USE_EOSPAC get_sg_eos_f,& finalize_sg_eos_f @@ -186,7 +189,8 @@ end function init_sg_SAP_Polynomial type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values end function init_sg_Helmholtz end interface -#endif ! SINGULARITY_USE_HELMHOLTZ +#endif +! SINGULARITY_USE_HELMHOLTZ interface integer(kind=c_int) function & init_sg_SpinerDependsRhoT(matindex, eos, filename, id, sg_mods_enabled, & @@ -212,7 +216,8 @@ end function init_sg_SpinerDependsRhoT type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values end function init_sg_SpinerDependsRhoSie end interface -#endif ! SINGULARITY_USE_SPINER_WITH_HDF5 +#endif +! SINGULARITY_USE_SPINER_WITH_HDF5 #ifdef SINGULARITY_USE_EOSPAC interface integer(kind=c_int) function & @@ -224,7 +229,8 @@ end function init_sg_SpinerDependsRhoSie type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values end function init_sg_eospac end interface -#endif ! SINGULARITY_USE_EOSPAC +#endif +! SINGULARITY_USE_EOSPAC interface integer(kind=c_int) function & get_sg_eos(nmat, ncell, cell_dim,& @@ -352,8 +358,8 @@ integer function init_sg_IdealGas_f(matindex, eos, gm1, Cv, & integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 @@ -375,8 +381,8 @@ integer function init_sg_Gruneisen_f(matindex, eos, C0, s1, s2, s3, G0, b,& integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 @@ -397,8 +403,8 @@ integer function init_sg_JWL_f(matindex, eos, A, B, R1, R2, w, rho0, Cv, & integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 @@ -419,8 +425,8 @@ integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, & integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 @@ -443,8 +449,8 @@ integer function init_sg_DavisReactants_f(matindex, eos, rho0, e0, P0, T0, & integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 @@ -467,8 +473,8 @@ integer function init_sg_SAP_Polynomial_f(matindex, eos, rho0, a0, a1, a2c, & integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 @@ -490,8 +496,8 @@ integer function init_sg_StiffGas_f(matindex, eos, gm1, Cv, & integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 @@ -512,8 +518,8 @@ integer function init_sg_NobleAbel_f(matindex, eos, gm1, Cv, & integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 @@ -536,8 +542,8 @@ integer function init_sg_Helmholtz_f(matindex, eos, filename, rad, gas, coul, io integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 @@ -560,8 +566,8 @@ integer function init_sg_SpinerDependsRhoT_f(matindex, eos, filename, id, & integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 @@ -584,8 +590,8 @@ integer function init_sg_SpinerDependsRhoSie_f(matindex, eos, filename, id, & integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 @@ -607,8 +613,8 @@ integer function init_sg_eospac_f(matindex, eos, id, sg_mods_enabled, & integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars - integer(kind=c_int), dimension(4) :: sg_mods_enabled_use - real(kind=8), dimension(6) :: sg_mods_values_use + integer(kind=c_int), target, dimension(4) :: sg_mods_enabled_use + real(kind=8), target, dimension(6) :: sg_mods_values_use sg_mods_enabled_use = 0 sg_mods_values_use = 0.d0 From cc23a92d4faac852c82995b2c29b66cb69165084 Mon Sep 17 00:00:00 2001 From: Daniel Holladay Date: Tue, 31 Oct 2023 15:13:57 -0600 Subject: [PATCH 20/20] Modified changelog. --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9bba35e26b0..f57aad0db4b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ - [[PR228]](https://github.com/lanl/singularity-eos/pull/228) added untracked header files in cmake - [[PR215]](https://github.com/lanl/singularity-eos/pull/215) and [[PR216]](https://github.com/lanl/singularity-eos/pull/216) fix duplicate definition of EPS and fix CI - [[PR232]](https://github.com/lanl/singularity-eos/pull/228) Fixed uninitialized cmake path variables +- [[PR308]](https://github.com/lanl/singularity-eos/pull/308) spack builds +fortran now compile via correct blocking out of interfaces via preprocessor ifdef ### Added (new features/APIs/variables/...) - [[PR306]](https://github.com/lanl/singularity-eos/pull/306) Added generic Evaluate method @@ -48,10 +49,12 @@ - [[PR234]](https://github.com/lanl/singularity-eos/pull/234) update ports-of-call to correct for undefined behavior in error handling - [[PR219]](https://github.com/lanl/singularity-eos/pull/219) Removed static analysis from re-git pipeline - [[PR233]](https://github.com/lanl/singularity-eos/pull/233) Exposed entropy for the EOS type (now required for future EOS) +- [[PR308]](https://github.com/lanl/singularity-eos/pull/308) Fortran initialization interface functions no longer require modifier arrays, they are optional parameters. ### Infrastructure (changes irrelevant to downstream codes) - [[PR190]](https://github.com/lanl/singularity-eos/pull/190) update CI on re-git - [[PR245]](https://github.com/lanl/singularity-eos/pull/245) Separating get_sg_eos to other files. Build/compilation improvements, warning fixes/suppression. +- [[PR308]](https://github.com/lanl/singularity-eos/pull/308) Added a fortran test. ### Removed (removing behavior/API/varaibles/...) - [[PR293]](https://github.com/lanl/singularity-eos/pull/293) Removing PTofRE function. This will no longer be callable downstream.