diff --git a/CHANGELOG.md b/CHANGELOG.md index 8f4dfaf436..cd69599b40 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,7 @@ - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface - [[PR373]](https://github.com/lanl/singularity-eos/pull/373) Initialize cache in `get_sg_eos*` functions - [[PR374]](https://github.com/lanl/singularity-eos/pull/374) Make the Davis EOS more numerically robust +- [[PR383]](https://github.com/lanl/singularity-eos/pull/383) Fix bug in step scaling for PTE solver ### Changed (changing behavior/API/variables/...) - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls diff --git a/CMakeLists.txt b/CMakeLists.txt index f0daaab319..124b2d780e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -259,7 +259,6 @@ if(SINGULARITY_USE_HIGH_RISK_MATH) target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_USE_HIGH_RISK_MATH) endif() - if(SINGULARITY_TEST_SESAME) target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_TEST_SESAME) endif() @@ -269,6 +268,9 @@ endif() if(SINGULARITY_USE_HELMHOLTZ) target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_USE_HELMHOLTZ) endif() +if(SINGULARITY_USE_SPINER_WITH_HDF5) + target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_USE_SPINER_WITH_HDF5) +endif() # ------------------------------------------------------------------------------# # Handle dependencies diff --git a/sesame2spiner/examples/unit_tests/copper.dat b/sesame2spiner/examples/unit_tests/copper.dat new file mode 100644 index 0000000000..16ba63d428 --- /dev/null +++ b/sesame2spiner/examples/unit_tests/copper.dat @@ -0,0 +1,20 @@ +# ====================================================================== +# Example input deck for sesame2spiner, +# a tool for converting eospac to spiner +# Author: Jonah Miller (jonahm@lanl.gov) +# © 2021-2024. Triad National Security, LLC. All rights reserved. This +# program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National +# 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. +# ====================================================================== + +matid=3337 +name=copper +rhomin=0. diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index d1ca6f89db..f8c0335e38 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -170,9 +170,9 @@ class PTESolverBase { PORTABLE_INLINE_FUNCTION virtual void Fixup() const { Real vsum = 0; - for (int m = 0; m < nmat; ++m) + for (int m = 0; m < nmat; ++m) { vsum += vfrac[m]; - PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive"); + } for (int m = 0; m < nmat; ++m) vfrac[m] *= robust::ratio(vfrac_total, vsum); } @@ -182,7 +182,6 @@ class PTESolverBase { void Finalize() { for (int m = 0; m < nmat; m++) { temp[m] *= Tnorm; - PORTABLE_REQUIRE(temp[m] >= 0., "Non-positive temperature returned"); u[m] *= uscale; press[m] *= uscale; } @@ -219,8 +218,9 @@ class PTESolverBase { // material m averaged over the full PTE volume rho_total = 0.0; for (int m = 0; m < nmat; ++m) { - PORTABLE_REQUIRE(vfrac[m] > 0., "Non-positive volume fraction"); - PORTABLE_REQUIRE(rho[m] > 0., "Non-positive density"); + PORTABLE_REQUIRE(vfrac[m] > 0., + "Non-positive volume fraction provided to PTE solver"); + PORTABLE_REQUIRE(rho[m] > 0., "Non-positive density provided to PTE solver"); rhobar[m] = rho[m] * vfrac[m]; rho_total += rhobar[m]; } @@ -353,12 +353,9 @@ class PTESolverBase { Asum += vfrac[m] * robust::ratio(press[m], temp[m]); rhoBsum += rho[m] * vfrac[m] * robust::ratio(sie[m], temp[m]); } - PORTABLE_REQUIRE(Tnorm > 0., "Non-positive temperature guess"); Asum *= uscale / Tnorm; rhoBsum /= Tnorm; - PORTABLE_REQUIRE(rhoBsum > 0., "Non-positive energy density"); Tideal = uscale / rhoBsum / Tnorm; - PORTABLE_REQUIRE(vfrac_total > 0., "Non-positive volume fraction sum"); Pideal = robust::ratio(Tnorm * Tideal * Asum, uscale * vfrac_total); } @@ -391,7 +388,8 @@ class PTESolverBase { const Real alpha = robust::ratio(Pideal, Tideal); for (int m = 0; m < nmat; ++m) { vfrac[m] *= robust::ratio(press[m], (temp[m] * alpha)); - if (robust::ratio(rhobar[m], vfrac[m]) < eos[m].RhoPmin(Tnorm * Tideal)) { + if (Tnorm * Tideal < 0 || + robust::ratio(rhobar[m], vfrac[m]) < eos[m].RhoPmin(Tnorm * Tideal)) { // abort because this is putting this material into a bad state for (int n = m; n >= 0; n--) vfrac[n] = vtemp[n]; @@ -671,6 +669,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { PORTABLE_INLINE_FUNCTION Real ScaleDx() const { using namespace mix_params; + // Each check reduces the scale further if necessary Real scale = 1.0; // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { @@ -684,8 +683,15 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { const Real rho_min = std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); const Real alpha_max = robust::ratio(rhobar[m], rho_min); + if (alpha_max < vfrac[m]) { + // Despite our best efforts, we're already in the unstable regime (i.e. + // dPdV_T > 0) so we would actually want to *increase* the step instead + // of decreasing it. As a result, this code doesn't work as intended and + // should be skipped. + continue; + } if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = robust::ratio(0.5 * alpha_max - vfrac[m], dx[m]); + scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]); } } // control how big of a step toward T = 0 is allowed @@ -701,8 +707,10 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale) { - if (scale == 1.0) { + Real TestUpdate(const Real scale, bool const cache_state = false) { + if (cache_state) { + // Store the current state in temp variables for first iteration of line + // search Ttemp = Tequil; for (int m = 0; m < nmat; ++m) vtemp[m] = vfrac[m]; @@ -879,6 +887,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase PORTABLE_INLINE_FUNCTION Real ScaleDx() const { using namespace mix_params; + // Each check reduces the scale further if necessary Real scale = 1.0; // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { @@ -890,6 +899,13 @@ class PTESolverFixedT : public mix_impl::PTESolverBase for (int m = 0; m < nmat; m++) { const Real rho_min = eos[m].RhoPmin(Tequil); const Real alpha_max = robust::ratio(rhobar[m], rho_min); + if (alpha_max < vfrac[m]) { + // Despite our best efforts, we're already in the unstable regime (i.e. + // dPdV_T > 0) so we would actually want to *increase* the step instead + // of decreasing it. As a result, this code doesn't work as intended and + // should be skipped. + continue; + } if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]); } @@ -903,8 +919,10 @@ class PTESolverFixedT : public mix_impl::PTESolverBase // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale) { - if (scale == 1.0) { + Real TestUpdate(const Real scale, bool const cache_state = false) { + if (cache_state) { + // Store the current state in temp variables for first iteration of line + // search for (int m = 0; m < nmat; ++m) vtemp[m] = vfrac[m]; } @@ -1093,6 +1111,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase PORTABLE_INLINE_FUNCTION Real ScaleDx() const { using namespace mix_params; + // Each check reduces the scale further if necessary Real scale = 1.0; // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { @@ -1106,6 +1125,13 @@ class PTESolverFixedP : public mix_impl::PTESolverBase const Real rho_min = std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); const Real alpha_max = robust::ratio(rhobar[m], rho_min); + if (alpha_max < vfrac[m]) { + // Despite our best efforts, we're already in the unstable regime (i.e. + // dPdV_T > 0) so we would actually want to *increase* the step instead + // of decreasing it. As a result, this code doesn't work as intended and + // should be skipped. + continue; + } if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); } @@ -1123,8 +1149,10 @@ class PTESolverFixedP : public mix_impl::PTESolverBase // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale) { - if (scale == 1.0) { + Real TestUpdate(const Real scale, bool const cache_state = false) { + if (cache_state) { + // Store the current state in temp variables for first iteration of line + // search Ttemp = Tequil; for (int m = 0; m < nmat; ++m) vtemp[m] = vfrac[m]; @@ -1332,6 +1360,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { PORTABLE_INLINE_FUNCTION Real ScaleDx() const { using namespace mix_params; + // Each check reduces the scale further if necessary Real scale = 1.0; for (int m = 0; m < nmat; ++m) { // control how big of a step toward vfrac = 0 is allowed @@ -1348,6 +1377,13 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { std::max(eos[m].RhoPmin(Tnorm * temp[m]), eos[m].RhoPmin(Tnorm * tt)); const Real alpha_max = robust::ratio(rhobar[m], rho_min); // control how big of a step toward rho = rho(Pmin) is allowed + if (alpha_max < vfrac[m]) { + // Despite our best efforts, we're already in the unstable regime (i.e. + // dPdV_T > 0) so we would actually want to *increase* the step instead + // of decreasing it. As a result, this code doesn't work as intended and + // should be skipped. + continue; + } if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); } @@ -1361,8 +1397,10 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale) const { - if (scale == 1.0) { + Real TestUpdate(const Real scale, bool const cache_state = false) const { + if (cache_state) { + // Store the current state in temp variables for first iteration of line + // search for (int m = 0; m < nmat; ++m) { vtemp[m] = vfrac[m]; utemp[m] = u[m]; @@ -1417,25 +1455,28 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { // possibly scale the update to stay within reasonable bounds Real scale = s.ScaleDx(); - // const Real scale_save = scale; + PORTABLE_REQUIRE(scale <= 1.0, "PTE Solver is attempting to increase the step size"); // Line search Real gradfdx = -2.0 * scale * err; - scale = 1.0; + scale = 1.0; // New scale for line search Real err_old = err; - err = s.TestUpdate(scale); + // Test the update and reset the cache the current state + err = s.TestUpdate(scale, true /* cache_state */); if (err > err_old + line_search_alpha * gradfdx) { - // backtrack - Real err_mid = s.TestUpdate(0.5); + // backtrack to middle of step + scale = 0.5; + Real err_mid = s.TestUpdate(scale); if (err_mid < err && err_mid < err_old) { + // We know the half step is better than both the full step and the + // prior result, so try a pseudo parabolic fit to the error and find its + // minimum. The `scale` value is bound between 0.75 and 0.25. scale = 0.75 + 0.5 * robust::ratio(err_mid - err, err - 2.0 * err_mid + err_old); - } else { - scale = line_search_fac; } - for (int line_iter = 0; line_iter < line_search_max_iter; line_iter++) { err = s.TestUpdate(scale); if (err < err_old + line_search_alpha * scale * gradfdx) break; + // shrink the step if the error isn't reduced enough scale *= line_search_fac; } } diff --git a/singularity-eos/eos/get_sg_eos_functors.hpp b/singularity-eos/eos/get_sg_eos_functors.hpp index 7550d13894..aee80f1933 100644 --- a/singularity-eos/eos/get_sg_eos_functors.hpp +++ b/singularity-eos/eos/get_sg_eos_functors.hpp @@ -63,7 +63,8 @@ struct init_functor { PORTABLE_INLINE_FUNCTION void operator()(const int i, const int tid, double &mass_sum, int &npte, - const Real t_mult, const Real s_mult, const Real p_mult) const { + double &vfrac_sum, const Real t_mult, const Real s_mult, + const Real p_mult) const { /* normalize mass fractions */ /* first find the mass sum */ /* also set idxs as the decrement of the eos offsets */ @@ -71,13 +72,19 @@ struct init_functor { for (int m = 0; m < nmat; ++m) { mass_sum += frac_mass_v(i, m); pte_idxs(tid, m) = eos_offsets_v(m) - 1; - frac_vol_v(i, m) = 0.0; + // Assume non-participating materials are all at cell density. PTE volumes + // will be overwritten + frac_vol_v(i, m) = frac_mass_v(i, m); + // Initialize all materials to the cell sie (if provided). PTE sie values + // will be overwritten. This also means that the sie for the PTE + // materials is the same as the cell PTE after removing small materials + frac_ie_v(i, m) = sie_v(i) * frac_mass_v(i, m) * s_mult; } for (int m = 0; m < nmat; ++m) { frac_mass_v(i, m) /= mass_sum; } check_all_vals(i); - // count the number of participating materials and zero the inputs + // count the number of participating materials npte = 0; for (int m = 0; m < nmat; ++m) { if (frac_mass_v(i, m) > mass_frac_cutoff) { @@ -85,10 +92,8 @@ struct init_functor { pte_idxs(tid, npte) = eos_offsets_v(m) - 1; pte_mats(tid, npte) = m; npte += 1; - } else { - frac_ie_v(i, m) = 0.0; } - // zero the inputs + // zero the PTE solver material inputs vfrac_pte(tid, m) = 0.0; sie_pte(tid, m) = 0.0; temp_pte(tid, m) = 0.0; @@ -98,10 +103,12 @@ struct init_functor { // NOTE: the volume fractions and densities need to be consistent with the // total specific volume since they are used to calculate internal // quantities for the PTE solver + vfrac_sum = 0.; for (int mp = 0; mp < npte; ++mp) { const int m = pte_mats(tid, mp); // Need to guess volume fractions vfrac_pte(tid, mp) = frac_mass_v(i, m); + vfrac_sum += vfrac_pte(tid, mp); // Calculate densities to be consistent with these volume fractions rho_pte(tid, mp) = frac_mass_v(i, m) / spvol_v(i) / vfrac_pte(tid, mp); temp_pte(tid, mp) = temp_v(i) * ev2k * t_mult; diff --git a/singularity-eos/eos/get_sg_eos_rho_e.cpp b/singularity-eos/eos/get_sg_eos_rho_e.cpp index 84db423ad3..9518b3176b 100644 --- a/singularity-eos/eos/get_sg_eos_rho_e.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_e.cpp @@ -37,8 +37,9 @@ void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v, const int32_t tid{small_loop ? iloop : token}; double mass_sum{0.0}; int npte{0}; + double vfrac_sum{0.0}; // initialize values for solver / lookup - i_func(i, tid, mass_sum, npte, 0.0, 1.0, 0.0); + i_func(i, tid, mass_sum, npte, vfrac_sum, 0.0, 1.0, 0.0); // need to initialize the scratch before it's used to avoid undefined behavior for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { solver_scratch(tid, idx) = 0.0; @@ -54,7 +55,7 @@ void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v, singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0)); // reset inputs PTESolverRhoT method( - npte, eos_inx, 1.0, sie_v(i), &rho_pte(tid, 0), &vfrac_pte(tid, 0), + npte, eos_inx, vfrac_sum, sie_v(i), &rho_pte(tid, 0), &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache, &solver_scratch(tid, 0)); pte_converged = PTESolver(method); diff --git a/singularity-eos/eos/get_sg_eos_rho_p.cpp b/singularity-eos/eos/get_sg_eos_rho_p.cpp index ea376b62b1..f735ba83ca 100644 --- a/singularity-eos/eos/get_sg_eos_rho_p.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_p.cpp @@ -38,8 +38,9 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v, const int32_t tid{small_loop ? iloop : token}; double mass_sum{0.0}; int npte{0}; + double vfrac_sum{0.0}; // initialize values for solver / lookup - i_func(i, tid, mass_sum, npte, 0.0, 0.0, 1.0); + i_func(i, tid, mass_sum, npte, vfrac_sum, 0.0, 0.0, 1.0); Real sie_tot_true{0.0}; // need to initialize the scratch before it's used to avoid undefined behavior for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { @@ -54,9 +55,9 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v, // eos accessor singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0)); PTESolverFixedP method( - npte, eos_inx, 1.0, press_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0), - &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache[0], - &solver_scratch(tid, 0)); + npte, eos_inx, vfrac_sum, press_pte(tid, 0), &rho_pte(tid, 0), + &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), + cache[0], &solver_scratch(tid, 0)); pte_converged = PTESolver(method); // calculate total sie for (int mp = 0; mp < npte; ++mp) { diff --git a/singularity-eos/eos/get_sg_eos_rho_t.cpp b/singularity-eos/eos/get_sg_eos_rho_t.cpp index df3d14c50b..cd292c811a 100644 --- a/singularity-eos/eos/get_sg_eos_rho_t.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_t.cpp @@ -38,8 +38,9 @@ void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v, const int32_t tid{small_loop ? iloop : token}; double mass_sum{0.0}; int npte{0}; + double vfrac_sum{0.0}; // initialize values for solver / lookup - i_func(i, tid, mass_sum, npte, 1.0, 0.0, 0.0); + i_func(i, tid, mass_sum, npte, vfrac_sum, 1.0, 0.0, 0.0); // calculate pte condition (lookup for 1 mat cell) Real sie_tot_true{0.0}; // need to initialize the scratch before it's used to avoid undefined behavior @@ -55,9 +56,9 @@ void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v, // eos accessor singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0)); PTESolverFixedT method( - npte, eos_inx, 1.0, temp_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0), - &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache, - &solver_scratch(tid, 0)); + npte, eos_inx, vfrac_sum, temp_pte(tid, 0), &rho_pte(tid, 0), + &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), + cache, &solver_scratch(tid, 0)); pte_converged = PTESolver(method); // calculate total internal energy for (int mp = 0; mp < npte; ++mp) { diff --git a/spack-repo/packages/singularity-eos/package.py b/spack-repo/packages/singularity-eos/package.py index ef59aa9b76..916294a1d1 100644 --- a/spack-repo/packages/singularity-eos/package.py +++ b/spack-repo/packages/singularity-eos/package.py @@ -112,6 +112,7 @@ class SingularityEos(CMakePackage, CudaPackage): depends_on("ports-of-call@1.4.2,1.5.2:", when="@:1.7.0") depends_on("ports-of-call@1.5.2:", when="@1.7.1:") + depends_on("ports-of-call@1.6.0:", when="@1.8.1:") # request HEAD of main branch depends_on("ports-of-call@main", when="@main") diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3adfb58278..4c654975ec 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -45,6 +45,12 @@ add_executable( test_eos_stellar_collapse.cpp ) +add_executable( + closure_unit_tests + catch2_define.cpp + test_closure_pte.cpp + ) + get_property(plugin_tests GLOBAL PROPERTY PLUGIN_TESTS) if (plugin_tests) add_executable( @@ -62,6 +68,7 @@ endif() if(SINGULARITY_TEST_SESAME) target_compile_definitions(eos_tabulated_unit_tests PRIVATE SINGULARITY_TEST_SESAME) + target_compile_definitions(closure_unit_tests PRIVATE SINGULARITY_TEST_SESAME) endif() if(SINGULARITY_TEST_STELLAR_COLLAPSE) target_compile_definitions(eos_tabulated_unit_tests @@ -74,6 +81,8 @@ 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) +target_link_libraries(closure_unit_tests PRIVATE Catch2::Catch2 + singularity-eos::singularity-eos) if (plugin_tests) target_link_libraries(eos_plugin_tests PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) @@ -82,6 +91,7 @@ include(Catch) catch_discover_tests(eos_analytic_unit_tests PROPERTIES TIMEOUT 60) catch_discover_tests(eos_infrastructure_tests PROPERTIES TIMEOUT 60) catch_discover_tests(eos_tabulated_unit_tests PROPERTIES TIMEOUT 60) +catch_discover_tests(closure_unit_tests PROPERTIES TIMEOUT 120) if (plugin_tests) catch_discover_tests(eos_plugin_tests PROPERTIES TIMEOUT 60) endif() diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp new file mode 100644 index 0000000000..ca01906c0c --- /dev/null +++ b/test/test_closure_pte.cpp @@ -0,0 +1,230 @@ +//------------------------------------------------------------------------------ +// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// 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. + +#ifdef SINGULARITY_BUILD_CLOSURE + +#include +#include +#include +#include + +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +#include +#include +#include +#include + +constexpr Real GPa = 1.0e10; +constexpr Real MJ_per_kg = 1.0e10; + +#ifdef SINGULARITY_TEST_SESAME +#ifdef SINGULARITY_USE_SPINER_WITH_HDF5 + +using singularity::DavisProducts; +using singularity::DavisReactants; +using singularity::IdealGas; +using singularity::MAX_NUM_LAMBDAS; +using singularity::PTESolverRhoT; +using singularity::PTESolverRhoTRequiredScratch; +using singularity::ShiftedEOS; +using singularity::SpinerEOSDependsRhoT; +using singularity::mix_impl::CacheAccessor; + +using EOS = singularity::Variant, DavisProducts, + DavisReactants, SpinerEOSDependsRhoT>; + +template +EOS *copy_eos_arr_to_device(const int num_eos, EOSArrT eos_arr) { + // Assumes that GetOnDevice() has already been called for each EOS in eos_arr + const size_t EOS_bytes = num_eos * sizeof(EOS); + EOS *v_EOS = (EOS *)PORTABLE_MALLOC(EOS_bytes); + const size_t bytes = num_eos * sizeof(EOS); + portableCopyToDevice(v_EOS, eos_arr.data(), bytes); + return v_EOS; +} + +template +void finalize_eos_arr(EOSArrT eos_arr) { + // Call Finalize on each EOS on the host + for (auto eos : eos_arr) { + eos.Finalize(); + } +} + +template +bool run_PTE_from_state(const int num_pte, EOS *v_EOS, const Real spvol_bulk, + const Real sie_bulk, ArrT mass_frac) { + // Calculate material densities (and corresponding energies) and the total + // volume fraction + std::vector vol_frac(num_pte); + std::vector densities(num_pte); + std::vector sies(num_pte); + std::vector temperatures(num_pte); + std::vector pressures(num_pte); + Real vfrac_sum = 0; + for (auto i = 0; i < num_pte; ++i) { + // Initial guess: all materials at the cell density + vol_frac[i] = mass_frac[i]; + densities[i] = mass_frac[i] / spvol_bulk / vol_frac[i]; + sies[i] = sie_bulk * mass_frac[i]; + vfrac_sum += vol_frac[i]; + // Initialize pressure and temperature arrays to zero + temperatures[i] = 0; + pressures[i] = 0; + } + + // Copy values to device (when available) + const size_t bytes = num_pte * sizeof(Real); + Real *v_densities = (Real *)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_densities, densities.data(), bytes); + Real *v_vol_frac = (Real *)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_vol_frac, vol_frac.data(), bytes); + Real *v_sies = (Real *)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_sies, sies.data(), bytes); + Real *v_temperatures = (Real *)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_temperatures, temperatures.data(), bytes); + Real *v_pressures = (Real *)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_pressures, pressures.data(), bytes); + + // Allocate scratch space for the PTE solver + const int pte_solver_scratch_size = PTESolverRhoTRequiredScratch(num_pte); + const size_t scratch_bytes = pte_solver_scratch_size * sizeof(Real); + Real *scratch = (double *)PORTABLE_MALLOC(scratch_bytes); + + // Allocate lambdas for all EOS and use an accessor to index into it + const size_t lambda_bytes = num_pte * MAX_NUM_LAMBDAS * sizeof(Real *); + Real *lambda_memory = (Real *)PORTABLE_MALLOC(lambda_bytes); + CacheAccessor lambdas = CacheAccessor(lambda_memory); + + // Solve the PTE system on device using a one-teration portableFor + bool pte_converged; + constexpr size_t bool_bytes = 1 * sizeof(bool); + bool *pte_converged_d = (bool *)PORTABLE_MALLOC(bool_bytes); + portableFor( + "Device execution of PTE Test", 0, 1, PORTABLE_LAMBDA(int i) { + PTESolverRhoT method( + num_pte, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, + v_temperatures, v_pressures, lambdas, scratch); + pte_converged_d[0] = PTESolver(method); + }); + portableCopyToHost(&pte_converged, pte_converged_d, bool_bytes); + + // Free temp memory + PORTABLE_FREE(lambda_memory); // Free entire lambda array + PORTABLE_FREE(scratch); + + // Free PTE values + PORTABLE_FREE(v_densities); + PORTABLE_FREE(v_vol_frac); + PORTABLE_FREE(v_sies); + PORTABLE_FREE(v_temperatures); + PORTABLE_FREE(v_pressures); + + return pte_converged; +} + +SCENARIO("Density-Temperature PTE Solver", "[PTE]") { + + GIVEN("Four equations of state") { + // Reference state + constexpr Real P0 = 1.0e6; // 1 bar + constexpr Real T0 = 296; // K + // Ideal gas air + constexpr Real rho0_air = 1e-03; // g/cc + constexpr Real Gruneisen_air = 0.4; + constexpr Real Cv_air = P0 / rho0_air / (Gruneisen_air * T0); + EOS air_eos = IdealGas(Gruneisen_air, Cv_air); + // Spiner copper EOS + constexpr int Cu_matid = 3337; + const std::string eos_file = "../materials.sp5"; + EOS copper_eos = SpinerEOSDependsRhoT(eos_file, Cu_matid); + // Davis Reactants EOS + constexpr Real rho0_DP = 1.890; + constexpr Real e0_DP = 0.; // erg / g + constexpr Real P0_DP = 0.; // microbar + constexpr Real T0_DP = 297; // K + constexpr Real A = 1.8 * 1.0e5; // 1.8 * sqrt(GPa) -> sqrt(microbar) + constexpr Real B = 4.6; + constexpr Real C = 0.34; + constexpr Real G0 = 0.56; + constexpr Real Z = 0.0; + constexpr Real alpha = 0.4265; + constexpr Real Cv_DP = 0.001074 * MJ_per_kg; // erg / g / K + EOS davis_r_eos = + DavisReactants(rho0_DP, e0_DP, P0_DP, T0_DP, A, B, C, G0, Z, alpha, Cv_DP); + // Davis Products EOS + constexpr Real a = 0.798311; + constexpr Real b = 0.58; + constexpr Real k = 1.35; + constexpr Real n = 2.66182; + constexpr Real vc = 0.75419; // cm^3 / g + constexpr Real pc = 3.2 * GPa; // microbar + constexpr Real Cv = 0.001072 * MJ_per_kg; // erg / g / K + constexpr Real Q = 4.115 * MJ_per_kg; + EOS davis_p_eos = DavisProducts(a, b, k, n, vc, pc, Cv).Modify(-Q); + + GIVEN("A difficult three-material state") { + // Define the state + constexpr Real spvol_bulk = 6.256037280402093e-01; + constexpr Real sie_bulk = -1.441692060406520e+10; + constexpr int num_pte = 3; + const std::array mass_frac = { + 1.10382442033331e-10, 0.124935312146569, 0.875064687743048}; + std::array eos_arr = {air_eos.GetOnDevice(), copper_eos.GetOnDevice(), + davis_p_eos.GetOnDevice()}; + + THEN("The PTE solver should converge") { + EOS *v_EOS = copy_eos_arr_to_device(num_pte, eos_arr); + const bool pte_converged = + run_PTE_from_state(num_pte, v_EOS, spvol_bulk, sie_bulk, mass_frac); + CHECK(pte_converged); + // Free EOS copies on device + PORTABLE_FREE(v_EOS); + } + finalize_eos_arr(eos_arr); + } + GIVEN("A difficult two-material state") { + // TODO: make this test converge + // Define the state + constexpr Real spvol_bulk = 4.010467628234189e-01; + constexpr Real sie_bulk = 3.290180957185173e+07; + constexpr int num_pte = 2; + const std::array mass_frac = {0.000312273191678158, + 0.999687726808322}; + std::array eos_arr = {air_eos.GetOnDevice(), + copper_eos.GetOnDevice()}; + THEN("The PTE solver should converge") { + EOS *v_EOS = copy_eos_arr_to_device(num_pte, eos_arr); + const bool pte_converged = + run_PTE_from_state(num_pte, v_EOS, spvol_bulk, sie_bulk, mass_frac); + // CHECK(pte_converged); + // Free EOS copies on device + PORTABLE_FREE(v_EOS); + } + // Deallocates device memory for each EOS (if applicable) + finalize_eos_arr(eos_arr); + } + // Clean up original EOS objects before they go out of scope + air_eos.Finalize(); // irrelevant because no data allocated + copper_eos.Finalize(); // actually does something + davis_r_eos.Finalize(); // irrelevant because no data allocated + davis_p_eos.Finalize(); // irrelevant because no data allocated + } +} +#endif // SINGULARITY_USE_SPINER_WITH_HDF5 +#endif // SINGULARITY_TEST_SESAME +#endif // SINGULARITY_BUILD_CLOSURE diff --git a/test/test_get_sg_eos.cpp b/test/test_get_sg_eos.cpp index 85c664ddf6..7a9ea0c6b6 100644 --- a/test/test_get_sg_eos.cpp +++ b/test/test_get_sg_eos.cpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2024. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -58,7 +58,7 @@ int run_sg_get_eos_tests() { Real vfrac_true[NMAT], ie_true[NMAT]; get_sg_eos(NMAT, 1, 1, -1, eos_offset, eoss, &cell_offset, &P_true, &pmax, &v_true, &spvol, &sie_tot_true, &T_true_ev, &bmod, &dpde, &cv, mfrac, vfrac_true, - ie_true, nullptr, nullptr, nullptr); + ie_true, nullptr, nullptr, nullptr, 1.0e-12); Real sie_tot_check = 0.0; for (int m = 0; m < NMAT; ++m) { const Real r_m = mfrac[m] / vfrac_true[m]; @@ -85,7 +85,7 @@ int run_sg_get_eos_tests() { Real p_check, vfrac_check[NMAT], ie_check[NMAT]; get_sg_eos(NMAT, 1, 1, -3, eos_offset, eoss, &cell_offset, &p_check, &pmax, &v_true, &spvol, &sie_tot_check, &T_true_ev, &bmod, &dpde, &cv, mfrac, vfrac_check, - ie_check, nullptr, nullptr, nullptr); + ie_check, nullptr, nullptr, nullptr, 1.0e-12); // check output pressure and sie, indicate failure if relative err is too large if (std::abs(P_true - p_check) / std::abs(P_true) > 1.e-5 || std::abs(sie_tot_true - sie_tot_check) / std::abs(sie_tot_true) > 1.e-5) { @@ -109,7 +109,7 @@ int run_sg_get_eos_tests() { Real t_check; get_sg_eos(NMAT, 1, 1, -2, eos_offset, eoss, &cell_offset, &P_true, &pmax, &v_true, &spvol, &sie_tot_check, &t_check, &bmod, &dpde, &cv, mfrac, vfrac_check, - ie_check, nullptr, nullptr, nullptr); + ie_check, nullptr, nullptr, nullptr, 1.0e-12); // check output temperature and sie, indicate failure if relative err is too large if (std::abs(T_true_ev - t_check) / std::abs(T_true_ev) > 1.e-5 || std::abs(sie_tot_true - sie_tot_check) / std::abs(sie_tot_true) > 1.e-5) {