From 0cdf72a01f5362b979590f280dddf3e640e5fd6d Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Fri, 5 Apr 2024 18:16:02 -0600 Subject: [PATCH 01/20] Added two PTE tests based on Sn beta and gamma phase Vinet EOS. --- doc/sphinx/src/using-closures.rst | 61 ++++++- test/CMakeLists.txt | 7 + test/pte_longtest_2phaseVinetSn.hpp | 123 ++++++++++++++ test/pte_test_2phaseVinetSn.hpp | 100 ++++++++++++ test/test_pte_phase.cpp | 245 ++++++++++++++++++++++++++++ 5 files changed, 534 insertions(+), 2 deletions(-) create mode 100644 test/pte_longtest_2phaseVinetSn.hpp create mode 100644 test/pte_test_2phaseVinetSn.hpp create mode 100644 test/test_pte_phase.cpp diff --git a/doc/sphinx/src/using-closures.rst b/doc/sphinx/src/using-closures.rst index fee7a9c1896..8bbfdf5e176 100644 --- a/doc/sphinx/src/using-closures.rst +++ b/doc/sphinx/src/using-closures.rst @@ -11,14 +11,21 @@ to assume pressure-temperature equilibrium (PTE). This implies that there's a single pressure in the cell and a single temperature, which are the same for all materials. +In Lagrangian codes, a PTE solver is needed, for example, when using +multi-phase models. In this case the situation corresponds to an +Eulerian multi-material simulation since each phase has its on set +of materials models, in particular the EOS is different in every phase. + ``singularity-eos`` provides several methods for finding a PTE -soluution. These methods differ in what they treat as independent +solution. These methods differ in what they treat as independent variables, and thus what precise system of equations they solve for. However they all share some common characteristics. Governing Equations and Assumptions ------------------------------------ +Using volume fractions +`````````````````````` Given a set of :math:`N` materials ranging from :math:`i = 0` to :math:`N-1`, each material can, in principle, have a distinct pressure :math:`P_i`, a distinct temperature :math:`T_i`, a distinct material @@ -53,6 +60,56 @@ where :math:`u = E/V` is the energy density by volume for total energy :math:`E` and total volume :math:`V`, and :math:`\varepsilon` is the total specific internal energy within that volume. +Using mass fractions +```````````````````` + +When using a PTE solver for a mixture of phases, it is common to use +*mass fractions* :math:`\lambda_i` instead of volume fractions, since +this is the quantity directly related to mass conservation in the sum of +all phases. We have + +.. math:: + + \sum_{i=0}^{N - 1} \lambda_i = 1 + +The density of the mixture is + +.. math:: + + \frac{1}{\rho} = V = \sum_{i=0}^{N - 1} \lambda_i V_i = \sum_{i=0}^{N - 1} \lambda_i \frac{1}{\rho_i} + +where :math:`V` is the specific volume, and the total specific internal energy is + +.. math:: + + \varepsilon = \sum_{i=0}^{N - 1} \lambda_i \varepsilon_i + +For a PTE solver, :math:`\varepsilon`, :math:`\rho`, and the set of :math:`\lambda_i` are +obtained from the rest of the host code and are used as input. Of course, we also give guesses +for :math:`\rho_i` and :math:`\varepsilon_i` to obtain fast convergance. Written in terms of the +volume fraction set of parameters, we have + +.. math:: + + \bar{\rho}_i = \lambda_i \rho + +Note that this gives a direct translation between the volume fraction and mass fraction equations. +It is only the nomenclature that is different. + +In the current version of the PTE solver :math:`f_i` and :math:`\rho_i` are given as independent +input but these parameters are adjusted during the calculation, while keeping the product, :math:`\bar{\rho_i}`, constant. +In order to use this when having only the set of :math:`\lambda_i` and :math:`\rho` available, +set a guess for :math:`\rho_i` (or :math:`f_i`) and set + +.. math:: + + f_i = \frac{\lambda_i \rho}{\rho_i} \qquad \qquad \left(\mbox{or } \rho_i = \frac{\lambda_i \rho}{f_i} \right) + + + +Pressure equilibrium +```````````````````` + The assumption of pressure equilibrium implies that .. math:: @@ -85,7 +142,7 @@ The Density-Energy Formulation One choice is to treat volume fractions and material energies as independent quantities. The material energies provide :math:`N` more unknowns. The equality of pressures provides :math:`N-1` additional -constraints. Additionally, the euqality of material temperatures, evaluated as +constraints. Additionally, the equality of material temperatures, evaluated as .. math:: diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 7169015d351..6ba3d884ddd 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -99,6 +99,13 @@ if(SINGULARITY_BUILD_CLOSURE) target_include_directories(test_pte PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) add_test(pte test_pte) endif() + if(SINGULARITY_USE_SPINER) + add_executable(test_pte_phase test_pte_phase.cpp) + target_link_libraries(test_pte_phase PRIVATE Catch2::Catch2 + singularity-eos::singularity-eos) + target_include_directories(test_pte_phase PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + add_test(pte_phase test_pte_phase) + endif() if(SINGULARITY_USE_KOKKOS AND SINGULARITY_USE_FORTRAN) add_executable(test_get_sg_eos test_get_sg_eos.cpp) target_link_libraries(test_get_sg_eos PRIVATE Catch2::Catch2 diff --git a/test/pte_longtest_2phaseVinetSn.hpp b/test/pte_longtest_2phaseVinetSn.hpp new file mode 100644 index 00000000000..2874db60989 --- /dev/null +++ b/test/pte_longtest_2phaseVinetSn.hpp @@ -0,0 +1,123 @@ +//------------------------------------------------------------------------------ +// © 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. +//------------------------------------------------------------------------------ + +#ifndef _SINGULARITY_EOS_TEST_PTE_TEST_2PHASEVINETSN_ +#define _SINGULARITY_EOS_TEST_PTE_TEST_2PHASEVINETSN_ + +#include + +#include +#include +#include + +constexpr int NMAT = 2; +constexpr int NTRIAL = 20; +constexpr int NPTS = NTRIAL * NMAT; +constexpr int HIST_SIZE = 10; + +constexpr Real in_rho_tot[NTRIAL] = {7.278,7.46221983,7.65110201,7.84476448,8.04332821, + 8.24691719,8.45565857,8.66968272,8.88912327,9.11411728,9.34480523, + 9.58133117,9.8238428,10.0724915,10.3274327,10.5888253,10.8568327, + 11.1316222,11.4133653,11.702238}; +constexpr Real in_sie_tot[NTRIAL] = {8.41323509e8,8.62513154e8,9.28665206e8,1.04377144e9, + 1.21199344e9,1.43767572e9,1.72535639e9,2.07977629e9,2.50588681e9, + 3.00885704e9,3.59408011e9,4.2671792100000005e9,5.03401308e9,5.90068122e9, + 6.87352878e9,7.95915117e9,9.16439846e9,1.04963795e10,1.19624656e10,1.35702945e10}; +constexpr Real in_lambda0 = 0.500480901; +constexpr Real in_lambda1 = 0.499519099; + +constexpr Real trial_vfrac0 = 47.0/100.0; +constexpr Real trial_vfrac1 = 53.0/100.0; + +constexpr Real out_press[NTRIAL] = {-3.29164604e-6,1.284232715e10,2.753234765e10,4.423652945000001e10, + 6.313670939999999e10,8.443021825e10,1.083303075e11,1.3506674800000002e11, + 1.64886557e11,1.9805482549999997e11,2.3485562650000003e11,2.755929975e11, + 3.20591986e11,3.70199756e11,4.247867485e11,4.84747905e11,5.505039405000001e11, + 6.225026735e11,7.012204135e11,7.871634035e11}; +constexpr Real out_temp[NTRIAL] = {298.,317.657834,338.05429449999997,359.18724,381.0524905, + 403.64396550000004,426.953795,450.97240750000003,475.6886215,501.0897275, + 527.1615735,553.8886454999999,581.254151,609.2401025,637.8273995,666.995915, + 696.7245780000001,726.9914615,757.7738645,789.048397}; +constexpr Real out_rho0[NTRIAL] = {7.285,7.44383086,7.61073863,7.78523193,7.9669718,8.15572753, + 8.35135008,8.55375288,8.76289814,8.97878664,9.20145031,9.43094663,9.6673544, + 9.91077059,10.1613078,10.4190924,10.6842631,10.9569698,11.2373726,11.5256413}; +constexpr Real out_rho1[NTRIAL] = {7.271,7.48073556,7.69197476,7.90533181,8.12131372,8.34035068, + 8.56281419,8.78903065,9.01929178,9.25386249,9.49298694,9.73689325,9.98579716, + 10.2399049,10.4994157,10.7645232,11.0354173, 11.3122855,11.5953136,11.8846866}; +constexpr Real out_vfrac0[NTRIAL] = {0.5,0.501717271,0.50313519,0.504308007,0.50527757, + 0.506076807,0.506731916,0.507263967,0.507690077,0.508024281,0.508278194, + 0.508461499,0.508582337,0.508647597,0.508663148,0.50863402, + 0.508564548,0.508458492,0.508319122,0.508149303}; +constexpr Real out_vfrac1[NTRIAL] = {0.5,0.498282729,0.49686481,0.495691993,0.49472243, + 0.493923193,0.493268084,0.492736033,0.492309923,0.491975719,0.491721806, + 0.491538501,0.491417663,0.491352403,0.491336852,0.49136598, + 0.491435452,0.491541508,0.491680878,0.491850697}; + +template +class LinearIndexer { + public: + PORTABLE_FUNCTION LinearIndexer() = default; + LinearIndexer(const T &t) : data_(t) {} + PORTABLE_INLINE_FUNCTION + auto &operator[](const int i) const { return data_(i); } + + private: + T data_; +}; + +template +class Indexer2D { + public: + PORTABLE_FUNCTION Indexer2D() = default; + PORTABLE_FUNCTION Indexer2D(const int j, const T &t) : j_(j), data_(t) {} + Indexer2D(const int j, const T &&t) = delete; // prevents r-value binding + PORTABLE_INLINE_FUNCTION + auto &operator[](const int i) const { return data_(j_, i); } + + private: + const int j_; + const T &data_; +}; + +template +inline void set_eos(T *eos) { + constexpr Real d2to40[39] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; + + singularity::EOS Snbeta = singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, + 0.658e09, 0.4419e07, d2to40); + singularity::EOS Sngamma = singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, 0.2161e07, + 1.025e09, 0.5051e07, d2to40); + eos[0] = Snbeta.GetOnDevice(); + eos[1] = Sngamma.GetOnDevice(); + return; +} + +template +inline void set_state(RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, + RealIndexer &&temp, EOSIndexer &&eos) { + + vfrac[0] = 47.0/100.0; + vfrac[1] = 53.0/100.0; + rho[0] = 7.278 * 0.500480901 / vfrac[0]; + rho[1] = 7.278 * 0.499519099 / vfrac[1]; + sie[0] = 8.41323509e08; + sie[1] = 8.41323509e08; + + + return; +} + +#endif // _SINGULARITY_EOS_TEST_PTE_TEST_2PHASEVINETSN_ diff --git a/test/pte_test_2phaseVinetSn.hpp b/test/pte_test_2phaseVinetSn.hpp new file mode 100644 index 00000000000..4f59d78eea2 --- /dev/null +++ b/test/pte_test_2phaseVinetSn.hpp @@ -0,0 +1,100 @@ +//------------------------------------------------------------------------------ +// © 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. +//------------------------------------------------------------------------------ + +#ifndef _SINGULARITY_EOS_TEST_PTE_TEST_2PHASEVINETSN_ +#define _SINGULARITY_EOS_TEST_PTE_TEST_2PHASEVINETSN_ + +#include + +#include +#include +#include + +constexpr int NMAT = 2; +constexpr int NTRIAL = 5; +constexpr int NPTS = NTRIAL * NMAT; +constexpr int HIST_SIZE = 10; + +constexpr Real in_rho_tot[5] = {7.27800000, 7.27981950, 7.28163945, 7.28345986, 7.28528073}; +constexpr Real in_sie_tot[5] = {8.41323509e08, 8.41325565e08, 8.41331734e08, 8.41342022e08, 8.41356432e08}; +constexpr Real in_lambda0 = 0.500480901; +constexpr Real in_lambda1 = 0.499519099; + +constexpr Real trial_vfrac0 = 47.0/100.0; +constexpr Real trial_vfrac1 = 53.0/100.0; + +constexpr Real out_press[5] = {-3.29164604e-6,1.19722694e8,2.3968114450000003e8,3.598077865e8,4.80104422e8}; +constexpr Real out_temp[5] = {298.,298.192952,298.38594450000005,298.5790125,298.7721535}; +constexpr Real out_rho0[5] = {7.28500000, 7.28653963, 7.28808705, 7.28963516, 7.29118416}; +constexpr Real out_rho1[5] = {7.27100000, 7.27309885, 7.27519088, 7.27728316, 7.27937551}; +constexpr Real out_vfrac0[5] = {0.5, 0.500019325, 0.500038138, 0.500056927, 0.500075679}; +constexpr Real out_vfrac1[5] = {0.5, 0.499980675, 0.499961862, 0.499943073, 0.499924321}; + +template +class LinearIndexer { + public: + PORTABLE_FUNCTION LinearIndexer() = default; + LinearIndexer(const T &t) : data_(t) {} + PORTABLE_INLINE_FUNCTION + auto &operator[](const int i) const { return data_(i); } + + private: + T data_; +}; + +template +class Indexer2D { + public: + PORTABLE_FUNCTION Indexer2D() = default; + PORTABLE_FUNCTION Indexer2D(const int j, const T &t) : j_(j), data_(t) {} + Indexer2D(const int j, const T &&t) = delete; // prevents r-value binding + PORTABLE_INLINE_FUNCTION + auto &operator[](const int i) const { return data_(j_, i); } + + private: + const int j_; + const T &data_; +}; + +template +inline void set_eos(T *eos) { + constexpr Real d2to40[39] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; + + singularity::EOS Snbeta = singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, + 0.658e09, 0.4419e07, d2to40); + singularity::EOS Sngamma = singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, 0.2161e07, + 1.025e09, 0.5051e07, d2to40); + eos[0] = Snbeta.GetOnDevice(); + eos[1] = Sngamma.GetOnDevice(); + return; +} + +template +inline void set_state(RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, + RealIndexer &&temp, EOSIndexer &&eos) { + + vfrac[0] = 47.0/100.0; + vfrac[1] = 53.0/100.0; + rho[0] = 7.278 * 0.500480901 / vfrac[0]; + rho[1] = 7.278 * 0.499519099 / vfrac[1]; + sie[0] = 8.41323509e08; + sie[1] = 8.41323509e08; + + + return; +} + +#endif // _SINGULARITY_EOS_TEST_PTE_TEST_2PHASEVINETSN_ diff --git a/test/test_pte_phase.cpp b/test/test_pte_phase.cpp new file mode 100644 index 00000000000..280ea64f5d0 --- /dev/null +++ b/test/test_pte_phase.cpp @@ -0,0 +1,245 @@ +//------------------------------------------------------------------------------ +// © 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. +//------------------------------------------------------------------------------ + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +// #include +#include +#include +#include + +using namespace singularity; + +using DataBox = Spiner::DataBox; + +int main(int argc, char *argv[]) { + + int nsuccess = 0; +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::initialize(); +#endif + { + // EOS +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View eos_v("eos", NMAT); + auto eos_hv = Kokkos::create_mirror_view(eos_v); +#else + std::vector eos_vec(NMAT); + PortableMDArray eos_hv(eos_vec.data(), NMAT); + PortableMDArray eos_v(eos_vec.data(), NMAT); +#endif + + LinearIndexer eos_h(eos_hv); + set_eos(eos_hv.data()); + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(eos_v, eos_hv); +#endif + + using EOSAccessor = LinearIndexer; + EOSAccessor eos(eos_v); + + // scratch required for PTE solver + auto nscratch_vars = PTESolverRhoTRequiredScratch(NMAT); + + // state vars +#ifdef PORTABILITY_STRATEGY_KOKKOS + using RView = Kokkos::View; + using atomic_view = Kokkos::MemoryTraits; + RView rho_v("rho", NPTS); + RView vfrac_v("vfrac", NPTS); + RView sie_v("sie", NPTS); + RView temp_v("temp", NPTS); + RView press_v("press", NPTS); + RView scratch_v("scratch", NTRIAL * nscratch_vars); + Kokkos::View hist_d("histogram", HIST_SIZE); + auto rho_vh = Kokkos::create_mirror_view(rho_v); + auto vfrac_vh = Kokkos::create_mirror_view(vfrac_v); + auto sie_vh = Kokkos::create_mirror_view(sie_v); + auto temp_vh = Kokkos::create_mirror_view(temp_v); + auto press_vh = Kokkos::create_mirror_view(press_v); + auto scratch_vh = Kokkos::create_mirror_view(scratch_v); + auto hist_vh = Kokkos::create_mirror_view(hist_d); + DataBox rho_d(rho_v.data(), NTRIAL, NMAT); + DataBox vfrac_d(vfrac_v.data(), NTRIAL, NMAT); + DataBox sie_d(sie_v.data(), NTRIAL, NMAT); + DataBox temp_d(temp_v.data(), NTRIAL, NMAT); + DataBox press_d(press_v.data(), NTRIAL, NMAT); + DataBox scratch_d(scratch_v.data(), NTRIAL * nscratch_vars); + DataBox rho_hm(rho_vh.data(), NTRIAL, NMAT); + DataBox vfrac_hm(vfrac_vh.data(), NTRIAL, NMAT); + DataBox sie_hm(sie_vh.data(), NTRIAL, NMAT); + DataBox temp_hm(temp_vh.data(), NTRIAL, NMAT); + DataBox press_hm(press_vh.data(), NTRIAL, NMAT); + DataBox scratch_hm(scratch_vh.data(), NTRIAL * nscratch_vars); +#else + DataBox rho_d(NTRIAL, NMAT); + DataBox vfrac_d(NTRIAL, NMAT); + DataBox sie_d(NTRIAL, NMAT); + DataBox temp_d(NTRIAL, NMAT); + DataBox press_d(NTRIAL, NMAT); + DataBox scratch_d(NTRIAL, nscratch_vars); + DataBox rho_hm = rho_d.slice(2, 0, NTRIAL); + DataBox vfrac_hm = vfrac_d.slice(2, 0, NTRIAL); + DataBox sie_hm = sie_d.slice(2, 0, NTRIAL); + DataBox temp_hm = temp_d.slice(2, 0, NTRIAL); + DataBox press_hm = press_d.slice(2, 0, NTRIAL); + DataBox scratch_hm = scratch_d.slice(2, 0, NTRIAL); + int hist_vh[HIST_SIZE]; + int *hist_d = hist_vh; +#endif + + // setup state + srand(time(NULL)); + for (int n = 0; n < NTRIAL; n++) { +// Indexer2D r(n, rho_hm); +// Indexer2D vf(n, vfrac_hm); +// Indexer2D e(n, sie_hm); +// Indexer2D t(n, temp_hm); +// set_state(r, vf, e, t, eos_h); + + vfrac_hm(n,0) = trial_vfrac0; + vfrac_hm(n,1) = trial_vfrac1; + rho_hm(n,0) = in_lambda0 * in_rho_tot[n] / vfrac_hm(n,0); + rho_hm(n,1) = in_lambda1 * in_rho_tot[n] / vfrac_hm(n,1); + // same sie in both phases gives sie_tot=sie below + sie_hm(n,0) = in_sie_tot[n]; + sie_hm(n,1) = in_sie_tot[n]; + } + for (int i = 0; i < HIST_SIZE; ++i) { + hist_vh[i] = 0; + } +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(rho_v, rho_vh); + Kokkos::deep_copy(vfrac_v, vfrac_vh); + Kokkos::deep_copy(sie_v, sie_vh); + Kokkos::deep_copy(temp_v, temp_vh); + Kokkos::deep_copy(press_v, press_vh); + Kokkos::deep_copy(hist_d, hist_vh); +#endif + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View nsuccess_d("n successes"); +#else + PortableMDArray nsuccess_d(&nsuccess, 1); +#endif + + auto start = std::chrono::high_resolution_clock::now(); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); +#endif + std::cout << "Starting PTE with " << NTRIAL << " trials." << std::endl << std::endl; + std::cout << "The trials have input set to: " << std::endl; + + for (int n = 0; n < NTRIAL; n++) { + std::cout << "Trial number: " << n << std::endl; + std::cout << "Total Specific Internal energy: \t\t" << in_sie_tot[n] << std::endl; + std::cout << "Total density: \t\t\t\t\t" << in_rho_tot[n] << std::endl; + std::cout << "Mass fractions: beta, gamma: \t\t\t" << in_lambda0 << ", " << in_lambda1 << std::endl; + std::cout << "Assuming volume fractions: beta, gamma: \t" << vfrac_hm(n,0) << ", " << vfrac_hm(n,1) << std::endl; + std::cout << "gives starting phase densities: beta, gamma: \t" << rho_hm(n,0) << ", " << rho_hm(n,1) << std::endl << std::endl; + } + + portableFor( + "PTE!", 0, NTRIAL, PORTABLE_LAMBDA(const int &t) { + Real *lambda[NMAT]; + for (int i = 0; i < NMAT; i++) { + lambda[i] = nullptr; + } + + Indexer2D rho(t, rho_d); + Indexer2D vfrac(t, vfrac_d); + Indexer2D sie(t, sie_d); + Indexer2D temp(t, temp_d); + Indexer2D press(t, press_d); + + Real sie_tot = 0.0; + Real rho_tot = 0.0; + for (int i = 0; i < NMAT; i++) { + rho_tot += rho[i] * vfrac[i]; + sie_tot += rho[i] * vfrac[i] * sie[i]; + } + sie_tot /= rho_tot; + + + const Real Tguess = + ApproxTemperatureFromRhoMatU(NMAT, eos, rho_tot * sie_tot, rho, vfrac); + + std::cout << "Tguess " << Tguess << std::endl; + + auto method = + PTESolverRhoT, decltype(lambda)>( + NMAT, eos, 1.0, sie_tot, rho, vfrac, sie, temp, press, lambda, + &scratch_d(t * nscratch_vars), Tguess); + bool success = PTESolver(method); + if (success) { + nsuccess_d() += 1; + } + hist_d[method.Niter()] += 1; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); +#endif + auto stop = std::chrono::high_resolution_clock::now(); + auto sum_time = std::chrono::duration_cast(stop - start); + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(nsuccess, nsuccess_d); + Kokkos::deep_copy(hist_vh, hist_d); +#endif + + Real milliseconds = sum_time.count() / 1e3; + + std::cout << std::endl; + std::cout << "Finished " << NTRIAL << " solves in " << milliseconds << " milliseconds" + << std::endl; + std::cout << "Solves/ms = " << NTRIAL / milliseconds << std::endl << std::endl; + + std::cout << "Results are: " << std::endl; + + for (int n = 0; n < NTRIAL; n++) { + std::cout << "Trial number: " << n << std::endl; + std::cout << "Total Specific Internal energy: \t" << sie_hm(n,0)*in_lambda0+sie_hm(n,1)*in_lambda1 << ", (" << in_sie_tot[n] << ")" << std::endl; + std::cout << "Total density: \t\t\t\t" << 1.0/(1.0/rho_hm(n,0)*in_lambda0+1.0/rho_hm(n,1)*in_lambda1) << ", (" << in_rho_tot[n] << ")" << std::endl; + std::cout << "Volume fractions: beta, gamma: \t\t" << vfrac_hm(n,0) << ", " << vfrac_hm(n,1) << ", (" << out_vfrac0[n] << ", " << out_vfrac1[n] << ")" << std::endl; + std::cout << "Density: beta, gamma: \t\t\t" << rho_hm(n,0) << ", " << rho_hm(n,1) << ", (" << out_rho0[n] << ", " << out_rho1[n] << ")" << std::endl; + std::cout << "Pressure: beta, gamma: \t\t\t" << press_hm(n,0) << ", " << press_hm(n,1) << ", (" << out_press[n] << ")" << std::endl; + std::cout << "Temperature: beta, gamma: \t\t" << temp_hm(n,0) << ", " << temp_hm(n,1) << ", (" << out_temp[n] << ")" << std::endl << std::endl; + } + + std::cout << "Success: " << nsuccess << " Failure: " << NTRIAL - nsuccess + << std::endl; + std::cout << "Histogram:\n" + << "iters\tcount\n" + << "----------------------\n"; + for (int i = 0; i < HIST_SIZE; ++i) { + std::cout << i << "\t" << hist_vh[i] << "\n"; + } + std::cout << std::endl; + } +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::finalize(); +#endif + + // poor-man's ctest integration + return (nsuccess >= 0.5 * NTRIAL) ? 0 : 1; +} From 8657626ddc4c629c9c9bb974667d5ad3201111d5 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Thu, 11 Apr 2024 18:48:00 -0600 Subject: [PATCH 02/20] Two test cases for testing PTESolveRhoT, added (pte_2phase and pte_3phase) for testing that EOSs used in Flag gives the same as when used in sg. --- test/CMakeLists.txt | 15 +- test/pte_test_2phaseVinetSn.hpp | 15 -- test/pte_test_5phaseSesameSn.hpp | 104 ++++++++ ...test_pte_phase.cpp => test_pte_2phase.cpp} | 11 +- test/test_pte_3phase.cpp | 251 ++++++++++++++++++ 5 files changed, 370 insertions(+), 26 deletions(-) create mode 100644 test/pte_test_5phaseSesameSn.hpp rename test/{test_pte_phase.cpp => test_pte_2phase.cpp} (95%) create mode 100644 test/test_pte_3phase.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 6ba3d884ddd..4a9a1442298 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -100,11 +100,18 @@ if(SINGULARITY_BUILD_CLOSURE) add_test(pte test_pte) endif() if(SINGULARITY_USE_SPINER) - add_executable(test_pte_phase test_pte_phase.cpp) - target_link_libraries(test_pte_phase PRIVATE Catch2::Catch2 + add_executable(test_pte_2phase test_pte_2phase.cpp) + target_link_libraries(test_pte_2phase PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) - target_include_directories(test_pte_phase PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) - add_test(pte_phase test_pte_phase) + target_include_directories(test_pte_2phase PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + add_test(pte_2phase test_pte_2phase) + endif() + if(SINGULARITY_USE_SPINER) + add_executable(test_pte_3phase test_pte_3phase.cpp) + target_link_libraries(test_pte_3phase PRIVATE Catch2::Catch2 + singularity-eos::singularity-eos) + target_include_directories(test_pte_3phase PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + add_test(pte_3phase test_pte_3phase) endif() if(SINGULARITY_USE_KOKKOS AND SINGULARITY_USE_FORTRAN) add_executable(test_get_sg_eos test_get_sg_eos.cpp) diff --git a/test/pte_test_2phaseVinetSn.hpp b/test/pte_test_2phaseVinetSn.hpp index 4f59d78eea2..086fffbc439 100644 --- a/test/pte_test_2phaseVinetSn.hpp +++ b/test/pte_test_2phaseVinetSn.hpp @@ -82,19 +82,4 @@ inline void set_eos(T *eos) { return; } -template -inline void set_state(RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, - RealIndexer &&temp, EOSIndexer &&eos) { - - vfrac[0] = 47.0/100.0; - vfrac[1] = 53.0/100.0; - rho[0] = 7.278 * 0.500480901 / vfrac[0]; - rho[1] = 7.278 * 0.499519099 / vfrac[1]; - sie[0] = 8.41323509e08; - sie[1] = 8.41323509e08; - - - return; -} - #endif // _SINGULARITY_EOS_TEST_PTE_TEST_2PHASEVINETSN_ diff --git a/test/pte_test_5phaseSesameSn.hpp b/test/pte_test_5phaseSesameSn.hpp new file mode 100644 index 00000000000..364a79434fd --- /dev/null +++ b/test/pte_test_5phaseSesameSn.hpp @@ -0,0 +1,104 @@ +//------------------------------------------------------------------------------ +// © 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. +//------------------------------------------------------------------------------ + +#ifndef _SINGULARITY_EOS_TEST_PTE_TEST_5PHASESESAMESN_ +#define _SINGULARITY_EOS_TEST_PTE_TEST_5PHASESESAMESN_ + +#include + +#include + +#include +#include +#include + +constexpr int NMAT = 3; +constexpr int NTRIAL = 5; +constexpr int NPTS = NTRIAL * NMAT; +constexpr int HIST_SIZE = 30; + +constexpr Real in_rho_tot[5] = {8.41382588, 8.41550873, 8.41719191, 8.41887543, 8.42055928}; +constexpr Real in_sie_tot[5] = {1.67602914e09, 1.67876139e09, 1.68149578e09, 1.68423198e09, 1.68696932e09}; +constexpr Real in_lambda0[5] = {0.663852654,0.660579238,0.656880511,0.652218675,0.646232106}; +constexpr Real in_lambda1[5] = {0.336093276,0.339262144,0.342703726,0.346800126,0.352135167}; +constexpr Real in_lambda2[5] = {0.000054070,0.000158618,0.000415763,0.000982199,0.001632727}; + +constexpr Real trial_vfrac0 = 650.0/1000.0; +constexpr Real trial_vfrac1 = 349.5/1000.0; +constexpr Real trial_vfrac2 = 0.5/1000.0; + + +constexpr Real out_press[5] = {1.0946046e11,1.0959183e11,1.0971241e11,1.0980857e11,1.0987166e11}; +constexpr Real out_temp[5] = {438.96969,438.98046,438.95507,438.84886,438.64322}; +constexpr Real out_rho0[5] = {8.35319483, 8.35425038, 8.35523024, 8.35603912, 8.35661359}; +constexpr Real out_rho1[5] = {8.53618806, 8.53734120, 8.53841145, 8.53929455, 8.53992118}; +constexpr Real out_rho2[5] = {8.53854993, 8.53970495, 8.54077710, 8.54166211, 8.54229063}; +constexpr Real out_sie0[5] = {1.54379257e09, 1.54520613e09, 1.54644403e09, 1.54728448e09, 1.54760327e09}; +constexpr Real out_sie1[5] = {1.93716882e09, 1.93864981e09, 1.93994981e09, 1.94084038e09, 1.94119298e09}; +constexpr Real out_sie2[5] = {2.01473728e09, 2.01621655e09, 2.01751522e09, 2.01840528e09, 2.01875845e09}; + +template +class LinearIndexer { + public: + PORTABLE_FUNCTION LinearIndexer() = default; + LinearIndexer(const T &t) : data_(t) {} + PORTABLE_INLINE_FUNCTION + auto &operator[](const int i) const { return data_(i); } + + private: + T data_; +}; + +template +class Indexer2D { + public: + PORTABLE_FUNCTION Indexer2D() = default; + PORTABLE_FUNCTION Indexer2D(const int j, const T &t) : j_(j), data_(t) {} + Indexer2D(const int j, const T &&t) = delete; // prevents r-value binding + PORTABLE_INLINE_FUNCTION + auto &operator[](const int i) const { return data_(j_, i); } + + private: + const int j_; + const T &data_; +}; + +template +inline void set_eos(T *eos) { + + symlink("/projects/shavano/dev/aematts/SEOS/PTEsolver_test_cases/sn2162-v01.bin","sesameu"); + + int SnbetaID = 2102; + int SngammaID = 2103; + int SndeltaID = 2104; + int SnhcpID = 2105; + int SnliquidID = 2106; + + bool invert_at_setup = true; + + singularity::EOS Snbeta = singularity::EOSPAC(SnbetaID); + singularity::EOS Sngamma = singularity::EOSPAC(SngammaID); +// singularity::EOS Sndelta = singularity::EOSPAC(SndeltaID); + singularity::EOS Snhcp = singularity::EOSPAC(SnhcpID); +// singularity::EOS Snliquid = singularity::EOSPAC(SnliquidID); + + eos[0] = Snbeta.GetOnDevice(); + eos[1] = Sngamma.GetOnDevice(); +// eos[x] = Sndelta.GetOnDevice(); + eos[2] = Snhcp.GetOnDevice(); +// eos[x] = Snliquid.GetOnDevice(); + return; +} + +#endif // _SINGULARITY_EOS_TEST_PTE_TEST_5PHASESESAMESN_ diff --git a/test/test_pte_phase.cpp b/test/test_pte_2phase.cpp similarity index 95% rename from test/test_pte_phase.cpp rename to test/test_pte_2phase.cpp index 280ea64f5d0..89e79b1c2ac 100644 --- a/test/test_pte_phase.cpp +++ b/test/test_pte_2phase.cpp @@ -49,7 +49,9 @@ int main(int argc, char *argv[]) { #endif LinearIndexer eos_h(eos_hv); + std::cout << "Before set_eos" << std::endl; set_eos(eos_hv.data()); + std::cout << "After set_eos" << std::endl; #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::deep_copy(eos_v, eos_hv); @@ -111,11 +113,6 @@ int main(int argc, char *argv[]) { // setup state srand(time(NULL)); for (int n = 0; n < NTRIAL; n++) { -// Indexer2D r(n, rho_hm); -// Indexer2D vf(n, vfrac_hm); -// Indexer2D e(n, sie_hm); -// Indexer2D t(n, temp_hm); -// set_state(r, vf, e, t, eos_h); vfrac_hm(n,0) = trial_vfrac0; vfrac_hm(n,1) = trial_vfrac1; @@ -220,10 +217,10 @@ int main(int argc, char *argv[]) { std::cout << "Trial number: " << n << std::endl; std::cout << "Total Specific Internal energy: \t" << sie_hm(n,0)*in_lambda0+sie_hm(n,1)*in_lambda1 << ", (" << in_sie_tot[n] << ")" << std::endl; std::cout << "Total density: \t\t\t\t" << 1.0/(1.0/rho_hm(n,0)*in_lambda0+1.0/rho_hm(n,1)*in_lambda1) << ", (" << in_rho_tot[n] << ")" << std::endl; - std::cout << "Volume fractions: beta, gamma: \t\t" << vfrac_hm(n,0) << ", " << vfrac_hm(n,1) << ", (" << out_vfrac0[n] << ", " << out_vfrac1[n] << ")" << std::endl; + std::cout << "Volume fractions: beta, gamma: \t\t" << vfrac_hm(n,0) << ", " << vfrac_hm(n,1) << std::endl; std::cout << "Density: beta, gamma: \t\t\t" << rho_hm(n,0) << ", " << rho_hm(n,1) << ", (" << out_rho0[n] << ", " << out_rho1[n] << ")" << std::endl; std::cout << "Pressure: beta, gamma: \t\t\t" << press_hm(n,0) << ", " << press_hm(n,1) << ", (" << out_press[n] << ")" << std::endl; - std::cout << "Temperature: beta, gamma: \t\t" << temp_hm(n,0) << ", " << temp_hm(n,1) << ", (" << out_temp[n] << ")" << std::endl << std::endl; + std::cout << "Temperature: beta, gamma: \t\t" << temp_hm(n,0) << ", " << temp_hm(n,1) << ", (" << out_temp[n] << ")" << std::endl << std::endl; } std::cout << "Success: " << nsuccess << " Failure: " << NTRIAL - nsuccess diff --git a/test/test_pte_3phase.cpp b/test/test_pte_3phase.cpp new file mode 100644 index 00000000000..54dced50708 --- /dev/null +++ b/test/test_pte_3phase.cpp @@ -0,0 +1,251 @@ +//------------------------------------------------------------------------------ +// © 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. +//------------------------------------------------------------------------------ + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace singularity; + +using DataBox = Spiner::DataBox; + +int main(int argc, char *argv[]) { + + int nsuccess = 0; +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::initialize(); +#endif + { + // EOS +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View eos_v("eos", NMAT); + auto eos_hv = Kokkos::create_mirror_view(eos_v); +#else + std::vector eos_vec(NMAT); + PortableMDArray eos_hv(eos_vec.data(), NMAT); + PortableMDArray eos_v(eos_vec.data(), NMAT); +#endif + + LinearIndexer eos_h(eos_hv); + std::cout << "Before set_eos" << std::endl; + set_eos(eos_hv.data()); + std::cout << "After set_eos" << std::endl; + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(eos_v, eos_hv); +#endif + + using EOSAccessor = LinearIndexer; + EOSAccessor eos(eos_v); + + // scratch required for PTE solver + auto nscratch_vars = PTESolverRhoTRequiredScratch(NMAT); + + // state vars +#ifdef PORTABILITY_STRATEGY_KOKKOS + using RView = Kokkos::View; + using atomic_view = Kokkos::MemoryTraits; + RView rho_v("rho", NPTS); + RView vfrac_v("vfrac", NPTS); + RView sie_v("sie", NPTS); + RView temp_v("temp", NPTS); + RView press_v("press", NPTS); + RView scratch_v("scratch", NTRIAL * nscratch_vars); + Kokkos::View hist_d("histogram", HIST_SIZE); + auto rho_vh = Kokkos::create_mirror_view(rho_v); + auto vfrac_vh = Kokkos::create_mirror_view(vfrac_v); + auto sie_vh = Kokkos::create_mirror_view(sie_v); + auto temp_vh = Kokkos::create_mirror_view(temp_v); + auto press_vh = Kokkos::create_mirror_view(press_v); + auto scratch_vh = Kokkos::create_mirror_view(scratch_v); + auto hist_vh = Kokkos::create_mirror_view(hist_d); + DataBox rho_d(rho_v.data(), NTRIAL, NMAT); + DataBox vfrac_d(vfrac_v.data(), NTRIAL, NMAT); + DataBox sie_d(sie_v.data(), NTRIAL, NMAT); + DataBox temp_d(temp_v.data(), NTRIAL, NMAT); + DataBox press_d(press_v.data(), NTRIAL, NMAT); + DataBox scratch_d(scratch_v.data(), NTRIAL * nscratch_vars); + DataBox rho_hm(rho_vh.data(), NTRIAL, NMAT); + DataBox vfrac_hm(vfrac_vh.data(), NTRIAL, NMAT); + DataBox sie_hm(sie_vh.data(), NTRIAL, NMAT); + DataBox temp_hm(temp_vh.data(), NTRIAL, NMAT); + DataBox press_hm(press_vh.data(), NTRIAL, NMAT); + DataBox scratch_hm(scratch_vh.data(), NTRIAL * nscratch_vars); +#else + DataBox rho_d(NTRIAL, NMAT); + DataBox vfrac_d(NTRIAL, NMAT); + DataBox sie_d(NTRIAL, NMAT); + DataBox temp_d(NTRIAL, NMAT); + DataBox press_d(NTRIAL, NMAT); + DataBox scratch_d(NTRIAL, nscratch_vars); + DataBox rho_hm = rho_d.slice(2, 0, NTRIAL); + DataBox vfrac_hm = vfrac_d.slice(2, 0, NTRIAL); + DataBox sie_hm = sie_d.slice(2, 0, NTRIAL); + DataBox temp_hm = temp_d.slice(2, 0, NTRIAL); + DataBox press_hm = press_d.slice(2, 0, NTRIAL); + DataBox scratch_hm = scratch_d.slice(2, 0, NTRIAL); + int hist_vh[HIST_SIZE]; + int *hist_d = hist_vh; +#endif + + // setup state + srand(time(NULL)); + for (int n = 0; n < NTRIAL; n++) { + + vfrac_hm(n,0) = trial_vfrac0; + vfrac_hm(n,1) = trial_vfrac1; +// vfrac_hm(n,x) = 0.0; + vfrac_hm(n,2) = trial_vfrac2; +// vfrac_hm(n,x) = 0.0; + rho_hm(n,0) = in_lambda0[n] * in_rho_tot[n] / vfrac_hm(n,0); + rho_hm(n,1) = in_lambda1[n] * in_rho_tot[n] / vfrac_hm(n,1); +// rho_hm(n,x) = 7.380; + rho_hm(n,2) = in_lambda2[n] * in_rho_tot[n] / vfrac_hm(n,2); +// rho_hm(n,x) = 7.116; + // same sie in both phases gives sie_tot=sie below + sie_hm(n,0) = in_sie_tot[n]; + sie_hm(n,1) = in_sie_tot[n]; +// sie_hm(n,x) = in_sie_tot[n]; + sie_hm(n,2) = in_sie_tot[n]; +// sie_hm(n,x) = in_sie_tot[n]; + } + for (int i = 0; i < HIST_SIZE; ++i) { + hist_vh[i] = 0; + } +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(rho_v, rho_vh); + Kokkos::deep_copy(vfrac_v, vfrac_vh); + Kokkos::deep_copy(sie_v, sie_vh); + Kokkos::deep_copy(temp_v, temp_vh); + Kokkos::deep_copy(press_v, press_vh); + Kokkos::deep_copy(hist_d, hist_vh); +#endif + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View nsuccess_d("n successes"); +#else + PortableMDArray nsuccess_d(&nsuccess, 1); +#endif + + auto start = std::chrono::high_resolution_clock::now(); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); +#endif + std::cout << "Starting PTE with " << NTRIAL << " trials." << std::endl << std::endl; + std::cout << "The trials have input set to: " << std::endl; + + for (int n = 0; n < NTRIAL; n++) { + std::cout << "Trial number: " << n << std::endl; + std::cout << "Total Specific Internal energy: \t\t\t" << in_sie_tot[n] << std::endl; + std::cout << "Total density: \t\t\t\t\t\t" << in_rho_tot[n] << std::endl; + std::cout << "Mass fractions: beta, gamma, hcp: \t\t\t" << in_lambda0[n] << ", " << in_lambda1[n] << ", " << in_lambda2[n] << std::endl; + std::cout << "Assuming volume fractions: beta, gamma, hcp: \t\t" << vfrac_hm(n,0) << ", " << vfrac_hm(n,1) << ", " << vfrac_hm(n,2) << std::endl; + std::cout << "gives starting phase densities: beta, gamma, hcp: \t" << rho_hm(n,0) << ", " << rho_hm(n,1) << ", " << rho_hm(n,2) << std::endl << std::endl; + } + + portableFor( + "PTE!", 0, NTRIAL, PORTABLE_LAMBDA(const int &t) { + Real *lambda[NMAT]; + for (int i = 0; i < NMAT; i++) { + lambda[i] = nullptr; + } + + Indexer2D rho(t, rho_d); + Indexer2D vfrac(t, vfrac_d); + Indexer2D sie(t, sie_d); + Indexer2D temp(t, temp_d); + Indexer2D press(t, press_d); + + Real sie_tot = 0.0; + Real rho_tot = 0.0; + for (int i = 0; i < NMAT; i++) { + rho_tot += rho[i] * vfrac[i]; + sie_tot += rho[i] * vfrac[i] * sie[i]; + } + sie_tot /= rho_tot; + + + const Real Tguess = + ApproxTemperatureFromRhoMatU(NMAT, eos, rho_tot * sie_tot, rho, vfrac); + + std::cout << "Tguess " << Tguess << std::endl; + + auto method = + PTESolverRhoT, decltype(lambda)>( + NMAT, eos, 1.0, sie_tot, rho, vfrac, sie, temp, press, lambda, + &scratch_d(t * nscratch_vars), Tguess); + bool success = PTESolver(method); + if (success) { + nsuccess_d() += 1; + } + hist_d[method.Niter()] += 1; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); +#endif + auto stop = std::chrono::high_resolution_clock::now(); + auto sum_time = std::chrono::duration_cast(stop - start); + +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::deep_copy(nsuccess, nsuccess_d); + Kokkos::deep_copy(hist_vh, hist_d); +#endif + + Real milliseconds = sum_time.count() / 1e3; + + std::cout << std::endl; + std::cout << "Finished " << NTRIAL << " solves in " << milliseconds << " milliseconds" + << std::endl; + std::cout << "Solves/ms = " << NTRIAL / milliseconds << std::endl << std::endl; + + std::cout << "Results are: " << std::endl; + + for (int n = 0; n < NTRIAL; n++) { + std::cout << "Trial number: " << n << std::endl; + std::cout << "Total Specific Internal energy: \t" << sie_hm(n,0)*in_lambda0[n]+sie_hm(n,1)*in_lambda1[n]+sie_hm(n,2)*in_lambda2[n] << ", (" << in_sie_tot[n] << ")" << std::endl; + std::cout << "Total density: \t\t\t\t" << 1.0/(1.0/rho_hm(n,0)*in_lambda0[n]+1.0/rho_hm(n,1)*in_lambda1[n]+1.0/rho_hm(n,2)*in_lambda2[n]) << ", (" << in_rho_tot[n] << ")" << std::endl; + std::cout << "Volume fractions: beta, gamma, hcp: \t" << vfrac_hm(n,0) << ", " << vfrac_hm(n,1) << ", " << vfrac_hm(n,2) << std::endl; + std::cout << "Density: beta, gamma, hcp: \t\t" << rho_hm(n,0) << ", " << rho_hm(n,1) << ", " << rho_hm(n,2) << ", (" << out_rho0[n] << ", " << out_rho1[n] << ", " << out_rho2[n] << ")" << std::endl; + std::cout << "Pressure: beta, gamma, hcp: \t\t" << press_hm(n,0) << ", " << press_hm(n,1) << ", " << press_hm(n,2) << ", (" << out_press[n] << ")" << std::endl; + std::cout << "Temperature: beta, gamma, hcp: \t\t" << temp_hm(n,0) << ", " << temp_hm(n,1) << ", " << temp_hm(n,2) << ", (" << out_temp[n] << ")" << std::endl; + std::cout << "Internal energy: beta, gamma, hcp: \t" << sie_hm(n,0) << ", " << sie_hm(n,1) << ", " << sie_hm(n,2) << ", (" << out_sie0[n] << ", " << out_sie1[n] << ", " << out_sie2[n] << ")" << std::endl << std::endl; + } + + std::cout << "Success: " << nsuccess << " Failure: " << NTRIAL - nsuccess + << std::endl; + std::cout << "Histogram:\n" + << "iters\tcount\n" + << "----------------------\n"; + for (int i = 0; i < HIST_SIZE; ++i) { + std::cout << i << "\t" << hist_vh[i] << "\n"; + } + std::cout << std::endl; + } +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::finalize(); +#endif + + // poor-man's ctest integration + return (nsuccess >= 0.5 * NTRIAL) ? 0 : 1; +} From 70983beae3f5edb0cc8725102a0584b78c8b3951 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Tue, 16 Apr 2024 19:22:35 -0600 Subject: [PATCH 03/20] Correctly formated... --- test/pte_longtest_2phaseVinetSn.hpp | 100 ++++++++++++++++------------ test/pte_test_2phaseVinetSn.hpp | 30 +++++---- test/pte_test_5phaseSesameSn.hpp | 51 ++++++++------ test/test_pte_2phase.cpp | 47 ++++++++----- test/test_pte_3phase.cpp | 74 ++++++++++++-------- 5 files changed, 181 insertions(+), 121 deletions(-) diff --git a/test/pte_longtest_2phaseVinetSn.hpp b/test/pte_longtest_2phaseVinetSn.hpp index 2874db60989..78282bdbf05 100644 --- a/test/pte_longtest_2phaseVinetSn.hpp +++ b/test/pte_longtest_2phaseVinetSn.hpp @@ -26,43 +26,54 @@ constexpr int NTRIAL = 20; constexpr int NPTS = NTRIAL * NMAT; constexpr int HIST_SIZE = 10; -constexpr Real in_rho_tot[NTRIAL] = {7.278,7.46221983,7.65110201,7.84476448,8.04332821, - 8.24691719,8.45565857,8.66968272,8.88912327,9.11411728,9.34480523, - 9.58133117,9.8238428,10.0724915,10.3274327,10.5888253,10.8568327, - 11.1316222,11.4133653,11.702238}; -constexpr Real in_sie_tot[NTRIAL] = {8.41323509e8,8.62513154e8,9.28665206e8,1.04377144e9, - 1.21199344e9,1.43767572e9,1.72535639e9,2.07977629e9,2.50588681e9, - 3.00885704e9,3.59408011e9,4.2671792100000005e9,5.03401308e9,5.90068122e9, - 6.87352878e9,7.95915117e9,9.16439846e9,1.04963795e10,1.19624656e10,1.35702945e10}; +constexpr Real in_rho_tot[NTRIAL] = { + 7.278, 7.46221983, 7.65110201, 7.84476448, 8.04332821, 8.24691719, 8.45565857, + 8.66968272, 8.88912327, 9.11411728, 9.34480523, 9.58133117, 9.8238428, 10.0724915, + 10.3274327, 10.5888253, 10.8568327, 11.1316222, 11.4133653, 11.702238}; +constexpr Real in_sie_tot[NTRIAL] = { + 8.41323509e8, 8.62513154e8, 9.28665206e8, 1.04377144e9, 1.21199344e9, + 1.43767572e9, 1.72535639e9, 2.07977629e9, 2.50588681e9, 3.00885704e9, + 3.59408011e9, 4.2671792100000005e9, 5.03401308e9, 5.90068122e9, 6.87352878e9, + 7.95915117e9, 9.16439846e9, 1.04963795e10, 1.19624656e10, 1.35702945e10}; constexpr Real in_lambda0 = 0.500480901; constexpr Real in_lambda1 = 0.499519099; -constexpr Real trial_vfrac0 = 47.0/100.0; -constexpr Real trial_vfrac1 = 53.0/100.0; - -constexpr Real out_press[NTRIAL] = {-3.29164604e-6,1.284232715e10,2.753234765e10,4.423652945000001e10, - 6.313670939999999e10,8.443021825e10,1.083303075e11,1.3506674800000002e11, - 1.64886557e11,1.9805482549999997e11,2.3485562650000003e11,2.755929975e11, - 3.20591986e11,3.70199756e11,4.247867485e11,4.84747905e11,5.505039405000001e11, - 6.225026735e11,7.012204135e11,7.871634035e11}; -constexpr Real out_temp[NTRIAL] = {298.,317.657834,338.05429449999997,359.18724,381.0524905, - 403.64396550000004,426.953795,450.97240750000003,475.6886215,501.0897275, - 527.1615735,553.8886454999999,581.254151,609.2401025,637.8273995,666.995915, - 696.7245780000001,726.9914615,757.7738645,789.048397}; -constexpr Real out_rho0[NTRIAL] = {7.285,7.44383086,7.61073863,7.78523193,7.9669718,8.15572753, - 8.35135008,8.55375288,8.76289814,8.97878664,9.20145031,9.43094663,9.6673544, - 9.91077059,10.1613078,10.4190924,10.6842631,10.9569698,11.2373726,11.5256413}; -constexpr Real out_rho1[NTRIAL] = {7.271,7.48073556,7.69197476,7.90533181,8.12131372,8.34035068, - 8.56281419,8.78903065,9.01929178,9.25386249,9.49298694,9.73689325,9.98579716, - 10.2399049,10.4994157,10.7645232,11.0354173, 11.3122855,11.5953136,11.8846866}; -constexpr Real out_vfrac0[NTRIAL] = {0.5,0.501717271,0.50313519,0.504308007,0.50527757, - 0.506076807,0.506731916,0.507263967,0.507690077,0.508024281,0.508278194, - 0.508461499,0.508582337,0.508647597,0.508663148,0.50863402, - 0.508564548,0.508458492,0.508319122,0.508149303}; -constexpr Real out_vfrac1[NTRIAL] = {0.5,0.498282729,0.49686481,0.495691993,0.49472243, - 0.493923193,0.493268084,0.492736033,0.492309923,0.491975719,0.491721806, - 0.491538501,0.491417663,0.491352403,0.491336852,0.49136598, - 0.491435452,0.491541508,0.491680878,0.491850697}; +constexpr Real trial_vfrac0 = 47.0 / 100.0; +constexpr Real trial_vfrac1 = 53.0 / 100.0; + +constexpr Real out_press[NTRIAL] = { + -3.29164604e-6, 1.284232715e10, 2.753234765e10, + 4.423652945000001e10, 6.313670939999999e10, 8.443021825e10, + 1.083303075e11, 1.3506674800000002e11, 1.64886557e11, + 1.9805482549999997e11, 2.3485562650000003e11, 2.755929975e11, + 3.20591986e11, 3.70199756e11, 4.247867485e11, + 4.84747905e11, 5.505039405000001e11, 6.225026735e11, + 7.012204135e11, 7.871634035e11}; +constexpr Real out_temp[NTRIAL] = {298., 317.657834, 338.05429449999997, + 359.18724, 381.0524905, 403.64396550000004, + 426.953795, 450.97240750000003, 475.6886215, + 501.0897275, 527.1615735, 553.8886454999999, + 581.254151, 609.2401025, 637.8273995, + 666.995915, 696.7245780000001, 726.9914615, + 757.7738645, 789.048397}; +constexpr Real out_rho0[NTRIAL] = { + 7.285, 7.44383086, 7.61073863, 7.78523193, 7.9669718, 8.15572753, 8.35135008, + 8.55375288, 8.76289814, 8.97878664, 9.20145031, 9.43094663, 9.6673544, 9.91077059, + 10.1613078, 10.4190924, 10.6842631, 10.9569698, 11.2373726, 11.5256413}; +constexpr Real out_rho1[NTRIAL] = { + 7.271, 7.48073556, 7.69197476, 7.90533181, 8.12131372, 8.34035068, 8.56281419, + 8.78903065, 9.01929178, 9.25386249, 9.49298694, 9.73689325, 9.98579716, 10.2399049, + 10.4994157, 10.7645232, 11.0354173, 11.3122855, 11.5953136, 11.8846866}; +constexpr Real out_vfrac0[NTRIAL] = {0.5, 0.501717271, 0.50313519, 0.504308007, + 0.50527757, 0.506076807, 0.506731916, 0.507263967, + 0.507690077, 0.508024281, 0.508278194, 0.508461499, + 0.508582337, 0.508647597, 0.508663148, 0.50863402, + 0.508564548, 0.508458492, 0.508319122, 0.508149303}; +constexpr Real out_vfrac1[NTRIAL] = {0.5, 0.498282729, 0.49686481, 0.495691993, + 0.49472243, 0.493923193, 0.493268084, 0.492736033, + 0.492309923, 0.491975719, 0.491721806, 0.491538501, + 0.491417663, 0.491352403, 0.491336852, 0.49136598, + 0.491435452, 0.491541508, 0.491680878, 0.491850697}; template class LinearIndexer { @@ -93,13 +104,15 @@ class Indexer2D { template inline void set_eos(T *eos) { constexpr Real d2to40[39] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; - - singularity::EOS Snbeta = singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, - 0.658e09, 0.4419e07, d2to40); - singularity::EOS Sngamma = singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, 0.2161e07, - 1.025e09, 0.5051e07, d2to40); + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; + + singularity::EOS Snbeta = + singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, 0.658e09, + 0.4419e07, d2to40); + singularity::EOS Sngamma = + singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, 0.2161e07, + 1.025e09, 0.5051e07, d2to40); eos[0] = Snbeta.GetOnDevice(); eos[1] = Sngamma.GetOnDevice(); return; @@ -109,14 +122,13 @@ template inline void set_state(RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, RealIndexer &&temp, EOSIndexer &&eos) { - vfrac[0] = 47.0/100.0; - vfrac[1] = 53.0/100.0; + vfrac[0] = 47.0 / 100.0; + vfrac[1] = 53.0 / 100.0; rho[0] = 7.278 * 0.500480901 / vfrac[0]; rho[1] = 7.278 * 0.499519099 / vfrac[1]; sie[0] = 8.41323509e08; sie[1] = 8.41323509e08; - return; } diff --git a/test/pte_test_2phaseVinetSn.hpp b/test/pte_test_2phaseVinetSn.hpp index 086fffbc439..316dbbb635b 100644 --- a/test/pte_test_2phaseVinetSn.hpp +++ b/test/pte_test_2phaseVinetSn.hpp @@ -26,16 +26,20 @@ constexpr int NTRIAL = 5; constexpr int NPTS = NTRIAL * NMAT; constexpr int HIST_SIZE = 10; -constexpr Real in_rho_tot[5] = {7.27800000, 7.27981950, 7.28163945, 7.28345986, 7.28528073}; -constexpr Real in_sie_tot[5] = {8.41323509e08, 8.41325565e08, 8.41331734e08, 8.41342022e08, 8.41356432e08}; +constexpr Real in_rho_tot[5] = {7.27800000, 7.27981950, 7.28163945, 7.28345986, + 7.28528073}; +constexpr Real in_sie_tot[5] = {8.41323509e08, 8.41325565e08, 8.41331734e08, + 8.41342022e08, 8.41356432e08}; constexpr Real in_lambda0 = 0.500480901; constexpr Real in_lambda1 = 0.499519099; -constexpr Real trial_vfrac0 = 47.0/100.0; -constexpr Real trial_vfrac1 = 53.0/100.0; +constexpr Real trial_vfrac0 = 47.0 / 100.0; +constexpr Real trial_vfrac1 = 53.0 / 100.0; -constexpr Real out_press[5] = {-3.29164604e-6,1.19722694e8,2.3968114450000003e8,3.598077865e8,4.80104422e8}; -constexpr Real out_temp[5] = {298.,298.192952,298.38594450000005,298.5790125,298.7721535}; +constexpr Real out_press[5] = {-3.29164604e-6, 1.19722694e8, 2.3968114450000003e8, + 3.598077865e8, 4.80104422e8}; +constexpr Real out_temp[5] = {298., 298.192952, 298.38594450000005, 298.5790125, + 298.7721535}; constexpr Real out_rho0[5] = {7.28500000, 7.28653963, 7.28808705, 7.28963516, 7.29118416}; constexpr Real out_rho1[5] = {7.27100000, 7.27309885, 7.27519088, 7.27728316, 7.27937551}; constexpr Real out_vfrac0[5] = {0.5, 0.500019325, 0.500038138, 0.500056927, 0.500075679}; @@ -70,13 +74,15 @@ class Indexer2D { template inline void set_eos(T *eos) { constexpr Real d2to40[39] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., - 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; - singularity::EOS Snbeta = singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, - 0.658e09, 0.4419e07, d2to40); - singularity::EOS Sngamma = singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, 0.2161e07, - 1.025e09, 0.5051e07, d2to40); + singularity::EOS Snbeta = + singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, 0.658e09, + 0.4419e07, d2to40); + singularity::EOS Sngamma = + singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, 0.2161e07, + 1.025e09, 0.5051e07, d2to40); eos[0] = Snbeta.GetOnDevice(); eos[1] = Sngamma.GetOnDevice(); return; diff --git a/test/pte_test_5phaseSesameSn.hpp b/test/pte_test_5phaseSesameSn.hpp index 364a79434fd..615b09703d7 100644 --- a/test/pte_test_5phaseSesameSn.hpp +++ b/test/pte_test_5phaseSesameSn.hpp @@ -28,25 +28,33 @@ constexpr int NTRIAL = 5; constexpr int NPTS = NTRIAL * NMAT; constexpr int HIST_SIZE = 30; -constexpr Real in_rho_tot[5] = {8.41382588, 8.41550873, 8.41719191, 8.41887543, 8.42055928}; -constexpr Real in_sie_tot[5] = {1.67602914e09, 1.67876139e09, 1.68149578e09, 1.68423198e09, 1.68696932e09}; -constexpr Real in_lambda0[5] = {0.663852654,0.660579238,0.656880511,0.652218675,0.646232106}; -constexpr Real in_lambda1[5] = {0.336093276,0.339262144,0.342703726,0.346800126,0.352135167}; -constexpr Real in_lambda2[5] = {0.000054070,0.000158618,0.000415763,0.000982199,0.001632727}; - -constexpr Real trial_vfrac0 = 650.0/1000.0; -constexpr Real trial_vfrac1 = 349.5/1000.0; -constexpr Real trial_vfrac2 = 0.5/1000.0; - - -constexpr Real out_press[5] = {1.0946046e11,1.0959183e11,1.0971241e11,1.0980857e11,1.0987166e11}; -constexpr Real out_temp[5] = {438.96969,438.98046,438.95507,438.84886,438.64322}; +constexpr Real in_rho_tot[5] = {8.41382588, 8.41550873, 8.41719191, 8.41887543, + 8.42055928}; +constexpr Real in_sie_tot[5] = {1.67602914e09, 1.67876139e09, 1.68149578e09, + 1.68423198e09, 1.68696932e09}; +constexpr Real in_lambda0[5] = {0.663852654, 0.660579238, 0.656880511, 0.652218675, + 0.646232106}; +constexpr Real in_lambda1[5] = {0.336093276, 0.339262144, 0.342703726, 0.346800126, + 0.352135167}; +constexpr Real in_lambda2[5] = {0.000054070, 0.000158618, 0.000415763, 0.000982199, + 0.001632727}; + +constexpr Real trial_vfrac0 = 650.0 / 1000.0; +constexpr Real trial_vfrac1 = 349.5 / 1000.0; +constexpr Real trial_vfrac2 = 0.5 / 1000.0; + +constexpr Real out_press[5] = {1.0946046e11, 1.0959183e11, 1.0971241e11, 1.0980857e11, + 1.0987166e11}; +constexpr Real out_temp[5] = {438.96969, 438.98046, 438.95507, 438.84886, 438.64322}; constexpr Real out_rho0[5] = {8.35319483, 8.35425038, 8.35523024, 8.35603912, 8.35661359}; constexpr Real out_rho1[5] = {8.53618806, 8.53734120, 8.53841145, 8.53929455, 8.53992118}; constexpr Real out_rho2[5] = {8.53854993, 8.53970495, 8.54077710, 8.54166211, 8.54229063}; -constexpr Real out_sie0[5] = {1.54379257e09, 1.54520613e09, 1.54644403e09, 1.54728448e09, 1.54760327e09}; -constexpr Real out_sie1[5] = {1.93716882e09, 1.93864981e09, 1.93994981e09, 1.94084038e09, 1.94119298e09}; -constexpr Real out_sie2[5] = {2.01473728e09, 2.01621655e09, 2.01751522e09, 2.01840528e09, 2.01875845e09}; +constexpr Real out_sie0[5] = {1.54379257e09, 1.54520613e09, 1.54644403e09, 1.54728448e09, + 1.54760327e09}; +constexpr Real out_sie1[5] = {1.93716882e09, 1.93864981e09, 1.93994981e09, 1.94084038e09, + 1.94119298e09}; +constexpr Real out_sie2[5] = {2.01473728e09, 2.01621655e09, 2.01751522e09, 2.01840528e09, + 2.01875845e09}; template class LinearIndexer { @@ -77,7 +85,8 @@ class Indexer2D { template inline void set_eos(T *eos) { - symlink("/projects/shavano/dev/aematts/SEOS/PTEsolver_test_cases/sn2162-v01.bin","sesameu"); + symlink("/projects/shavano/dev/aematts/SEOS/PTEsolver_test_cases/sn2162-v01.bin", + "sesameu"); int SnbetaID = 2102; int SngammaID = 2103; @@ -89,15 +98,15 @@ inline void set_eos(T *eos) { singularity::EOS Snbeta = singularity::EOSPAC(SnbetaID); singularity::EOS Sngamma = singularity::EOSPAC(SngammaID); -// singularity::EOS Sndelta = singularity::EOSPAC(SndeltaID); + // singularity::EOS Sndelta = singularity::EOSPAC(SndeltaID); singularity::EOS Snhcp = singularity::EOSPAC(SnhcpID); -// singularity::EOS Snliquid = singularity::EOSPAC(SnliquidID); + // singularity::EOS Snliquid = singularity::EOSPAC(SnliquidID); eos[0] = Snbeta.GetOnDevice(); eos[1] = Sngamma.GetOnDevice(); -// eos[x] = Sndelta.GetOnDevice(); + // eos[x] = Sndelta.GetOnDevice(); eos[2] = Snhcp.GetOnDevice(); -// eos[x] = Snliquid.GetOnDevice(); + // eos[x] = Snliquid.GetOnDevice(); return; } diff --git a/test/test_pte_2phase.cpp b/test/test_pte_2phase.cpp index 89e79b1c2ac..d7be8bb5a95 100644 --- a/test/test_pte_2phase.cpp +++ b/test/test_pte_2phase.cpp @@ -114,13 +114,13 @@ int main(int argc, char *argv[]) { srand(time(NULL)); for (int n = 0; n < NTRIAL; n++) { - vfrac_hm(n,0) = trial_vfrac0; - vfrac_hm(n,1) = trial_vfrac1; - rho_hm(n,0) = in_lambda0 * in_rho_tot[n] / vfrac_hm(n,0); - rho_hm(n,1) = in_lambda1 * in_rho_tot[n] / vfrac_hm(n,1); + vfrac_hm(n, 0) = trial_vfrac0; + vfrac_hm(n, 1) = trial_vfrac1; + rho_hm(n, 0) = in_lambda0 * in_rho_tot[n] / vfrac_hm(n, 0); + rho_hm(n, 1) = in_lambda1 * in_rho_tot[n] / vfrac_hm(n, 1); // same sie in both phases gives sie_tot=sie below - sie_hm(n,0) = in_sie_tot[n]; - sie_hm(n,1) = in_sie_tot[n]; + sie_hm(n, 0) = in_sie_tot[n]; + sie_hm(n, 1) = in_sie_tot[n]; } for (int i = 0; i < HIST_SIZE; ++i) { hist_vh[i] = 0; @@ -151,9 +151,13 @@ int main(int argc, char *argv[]) { std::cout << "Trial number: " << n << std::endl; std::cout << "Total Specific Internal energy: \t\t" << in_sie_tot[n] << std::endl; std::cout << "Total density: \t\t\t\t\t" << in_rho_tot[n] << std::endl; - std::cout << "Mass fractions: beta, gamma: \t\t\t" << in_lambda0 << ", " << in_lambda1 << std::endl; - std::cout << "Assuming volume fractions: beta, gamma: \t" << vfrac_hm(n,0) << ", " << vfrac_hm(n,1) << std::endl; - std::cout << "gives starting phase densities: beta, gamma: \t" << rho_hm(n,0) << ", " << rho_hm(n,1) << std::endl << std::endl; + std::cout << "Mass fractions: beta, gamma: \t\t\t" << in_lambda0 << ", " + << in_lambda1 << std::endl; + std::cout << "Assuming volume fractions: beta, gamma: \t" << vfrac_hm(n, 0) << ", " + << vfrac_hm(n, 1) << std::endl; + std::cout << "gives starting phase densities: beta, gamma: \t" << rho_hm(n, 0) + << ", " << rho_hm(n, 1) << std::endl + << std::endl; } portableFor( @@ -177,11 +181,10 @@ int main(int argc, char *argv[]) { } sie_tot /= rho_tot; - const Real Tguess = ApproxTemperatureFromRhoMatU(NMAT, eos, rho_tot * sie_tot, rho, vfrac); - std::cout << "Tguess " << Tguess << std::endl; + std::cout << "Tguess " << Tguess << std::endl; auto method = PTESolverRhoT, decltype(lambda)>( @@ -215,12 +218,22 @@ int main(int argc, char *argv[]) { for (int n = 0; n < NTRIAL; n++) { std::cout << "Trial number: " << n << std::endl; - std::cout << "Total Specific Internal energy: \t" << sie_hm(n,0)*in_lambda0+sie_hm(n,1)*in_lambda1 << ", (" << in_sie_tot[n] << ")" << std::endl; - std::cout << "Total density: \t\t\t\t" << 1.0/(1.0/rho_hm(n,0)*in_lambda0+1.0/rho_hm(n,1)*in_lambda1) << ", (" << in_rho_tot[n] << ")" << std::endl; - std::cout << "Volume fractions: beta, gamma: \t\t" << vfrac_hm(n,0) << ", " << vfrac_hm(n,1) << std::endl; - std::cout << "Density: beta, gamma: \t\t\t" << rho_hm(n,0) << ", " << rho_hm(n,1) << ", (" << out_rho0[n] << ", " << out_rho1[n] << ")" << std::endl; - std::cout << "Pressure: beta, gamma: \t\t\t" << press_hm(n,0) << ", " << press_hm(n,1) << ", (" << out_press[n] << ")" << std::endl; - std::cout << "Temperature: beta, gamma: \t\t" << temp_hm(n,0) << ", " << temp_hm(n,1) << ", (" << out_temp[n] << ")" << std::endl << std::endl; + std::cout << "Total Specific Internal energy: \t" + << sie_hm(n, 0) * in_lambda0 + sie_hm(n, 1) * in_lambda1 << ", (" + << in_sie_tot[n] << ")" << std::endl; + std::cout << "Total density: \t\t\t\t" + << 1.0 / + (1.0 / rho_hm(n, 0) * in_lambda0 + 1.0 / rho_hm(n, 1) * in_lambda1) + << ", (" << in_rho_tot[n] << ")" << std::endl; + std::cout << "Volume fractions: beta, gamma: \t\t" << vfrac_hm(n, 0) << ", " + << vfrac_hm(n, 1) << std::endl; + std::cout << "Density: beta, gamma: \t\t\t" << rho_hm(n, 0) << ", " << rho_hm(n, 1) + << ", (" << out_rho0[n] << ", " << out_rho1[n] << ")" << std::endl; + std::cout << "Pressure: beta, gamma: \t\t\t" << press_hm(n, 0) << ", " + << press_hm(n, 1) << ", (" << out_press[n] << ")" << std::endl; + std::cout << "Temperature: beta, gamma: \t\t" << temp_hm(n, 0) << ", " + << temp_hm(n, 1) << ", (" << out_temp[n] << ")" << std::endl + << std::endl; } std::cout << "Success: " << nsuccess << " Failure: " << NTRIAL - nsuccess diff --git a/test/test_pte_3phase.cpp b/test/test_pte_3phase.cpp index 54dced50708..a44a8dc90ba 100644 --- a/test/test_pte_3phase.cpp +++ b/test/test_pte_3phase.cpp @@ -113,22 +113,22 @@ int main(int argc, char *argv[]) { srand(time(NULL)); for (int n = 0; n < NTRIAL; n++) { - vfrac_hm(n,0) = trial_vfrac0; - vfrac_hm(n,1) = trial_vfrac1; -// vfrac_hm(n,x) = 0.0; - vfrac_hm(n,2) = trial_vfrac2; -// vfrac_hm(n,x) = 0.0; - rho_hm(n,0) = in_lambda0[n] * in_rho_tot[n] / vfrac_hm(n,0); - rho_hm(n,1) = in_lambda1[n] * in_rho_tot[n] / vfrac_hm(n,1); -// rho_hm(n,x) = 7.380; - rho_hm(n,2) = in_lambda2[n] * in_rho_tot[n] / vfrac_hm(n,2); -// rho_hm(n,x) = 7.116; + vfrac_hm(n, 0) = trial_vfrac0; + vfrac_hm(n, 1) = trial_vfrac1; + // vfrac_hm(n,x) = 0.0; + vfrac_hm(n, 2) = trial_vfrac2; + // vfrac_hm(n,x) = 0.0; + rho_hm(n, 0) = in_lambda0[n] * in_rho_tot[n] / vfrac_hm(n, 0); + rho_hm(n, 1) = in_lambda1[n] * in_rho_tot[n] / vfrac_hm(n, 1); + // rho_hm(n,x) = 7.380; + rho_hm(n, 2) = in_lambda2[n] * in_rho_tot[n] / vfrac_hm(n, 2); + // rho_hm(n,x) = 7.116; // same sie in both phases gives sie_tot=sie below - sie_hm(n,0) = in_sie_tot[n]; - sie_hm(n,1) = in_sie_tot[n]; -// sie_hm(n,x) = in_sie_tot[n]; - sie_hm(n,2) = in_sie_tot[n]; -// sie_hm(n,x) = in_sie_tot[n]; + sie_hm(n, 0) = in_sie_tot[n]; + sie_hm(n, 1) = in_sie_tot[n]; + // sie_hm(n,x) = in_sie_tot[n]; + sie_hm(n, 2) = in_sie_tot[n]; + // sie_hm(n,x) = in_sie_tot[n]; } for (int i = 0; i < HIST_SIZE; ++i) { hist_vh[i] = 0; @@ -159,9 +159,13 @@ int main(int argc, char *argv[]) { std::cout << "Trial number: " << n << std::endl; std::cout << "Total Specific Internal energy: \t\t\t" << in_sie_tot[n] << std::endl; std::cout << "Total density: \t\t\t\t\t\t" << in_rho_tot[n] << std::endl; - std::cout << "Mass fractions: beta, gamma, hcp: \t\t\t" << in_lambda0[n] << ", " << in_lambda1[n] << ", " << in_lambda2[n] << std::endl; - std::cout << "Assuming volume fractions: beta, gamma, hcp: \t\t" << vfrac_hm(n,0) << ", " << vfrac_hm(n,1) << ", " << vfrac_hm(n,2) << std::endl; - std::cout << "gives starting phase densities: beta, gamma, hcp: \t" << rho_hm(n,0) << ", " << rho_hm(n,1) << ", " << rho_hm(n,2) << std::endl << std::endl; + std::cout << "Mass fractions: beta, gamma, hcp: \t\t\t" << in_lambda0[n] << ", " + << in_lambda1[n] << ", " << in_lambda2[n] << std::endl; + std::cout << "Assuming volume fractions: beta, gamma, hcp: \t\t" << vfrac_hm(n, 0) + << ", " << vfrac_hm(n, 1) << ", " << vfrac_hm(n, 2) << std::endl; + std::cout << "gives starting phase densities: beta, gamma, hcp: \t" << rho_hm(n, 0) + << ", " << rho_hm(n, 1) << ", " << rho_hm(n, 2) << std::endl + << std::endl; } portableFor( @@ -185,11 +189,10 @@ int main(int argc, char *argv[]) { } sie_tot /= rho_tot; - const Real Tguess = ApproxTemperatureFromRhoMatU(NMAT, eos, rho_tot * sie_tot, rho, vfrac); - std::cout << "Tguess " << Tguess << std::endl; + std::cout << "Tguess " << Tguess << std::endl; auto method = PTESolverRhoT, decltype(lambda)>( @@ -223,13 +226,30 @@ int main(int argc, char *argv[]) { for (int n = 0; n < NTRIAL; n++) { std::cout << "Trial number: " << n << std::endl; - std::cout << "Total Specific Internal energy: \t" << sie_hm(n,0)*in_lambda0[n]+sie_hm(n,1)*in_lambda1[n]+sie_hm(n,2)*in_lambda2[n] << ", (" << in_sie_tot[n] << ")" << std::endl; - std::cout << "Total density: \t\t\t\t" << 1.0/(1.0/rho_hm(n,0)*in_lambda0[n]+1.0/rho_hm(n,1)*in_lambda1[n]+1.0/rho_hm(n,2)*in_lambda2[n]) << ", (" << in_rho_tot[n] << ")" << std::endl; - std::cout << "Volume fractions: beta, gamma, hcp: \t" << vfrac_hm(n,0) << ", " << vfrac_hm(n,1) << ", " << vfrac_hm(n,2) << std::endl; - std::cout << "Density: beta, gamma, hcp: \t\t" << rho_hm(n,0) << ", " << rho_hm(n,1) << ", " << rho_hm(n,2) << ", (" << out_rho0[n] << ", " << out_rho1[n] << ", " << out_rho2[n] << ")" << std::endl; - std::cout << "Pressure: beta, gamma, hcp: \t\t" << press_hm(n,0) << ", " << press_hm(n,1) << ", " << press_hm(n,2) << ", (" << out_press[n] << ")" << std::endl; - std::cout << "Temperature: beta, gamma, hcp: \t\t" << temp_hm(n,0) << ", " << temp_hm(n,1) << ", " << temp_hm(n,2) << ", (" << out_temp[n] << ")" << std::endl; - std::cout << "Internal energy: beta, gamma, hcp: \t" << sie_hm(n,0) << ", " << sie_hm(n,1) << ", " << sie_hm(n,2) << ", (" << out_sie0[n] << ", " << out_sie1[n] << ", " << out_sie2[n] << ")" << std::endl << std::endl; + std::cout << "Total Specific Internal energy: \t" + << sie_hm(n, 0) * in_lambda0[n] + sie_hm(n, 1) * in_lambda1[n] + + sie_hm(n, 2) * in_lambda2[n] + << ", (" << in_sie_tot[n] << ")" << std::endl; + std::cout << "Total density: \t\t\t\t" + << 1.0 / (1.0 / rho_hm(n, 0) * in_lambda0[n] + + 1.0 / rho_hm(n, 1) * in_lambda1[n] + + 1.0 / rho_hm(n, 2) * in_lambda2[n]) + << ", (" << in_rho_tot[n] << ")" << std::endl; + std::cout << "Volume fractions: beta, gamma, hcp: \t" << vfrac_hm(n, 0) << ", " + << vfrac_hm(n, 1) << ", " << vfrac_hm(n, 2) << std::endl; + std::cout << "Density: beta, gamma, hcp: \t\t" << rho_hm(n, 0) << ", " + << rho_hm(n, 1) << ", " << rho_hm(n, 2) << ", (" << out_rho0[n] << ", " + << out_rho1[n] << ", " << out_rho2[n] << ")" << std::endl; + std::cout << "Pressure: beta, gamma, hcp: \t\t" << press_hm(n, 0) << ", " + << press_hm(n, 1) << ", " << press_hm(n, 2) << ", (" << out_press[n] + << ")" << std::endl; + std::cout << "Temperature: beta, gamma, hcp: \t\t" << temp_hm(n, 0) << ", " + << temp_hm(n, 1) << ", " << temp_hm(n, 2) << ", (" << out_temp[n] << ")" + << std::endl; + std::cout << "Internal energy: beta, gamma, hcp: \t" << sie_hm(n, 0) << ", " + << sie_hm(n, 1) << ", " << sie_hm(n, 2) << ", (" << out_sie0[n] << ", " + << out_sie1[n] << ", " << out_sie2[n] << ")" << std::endl + << std::endl; } std::cout << "Success: " << nsuccess << " Failure: " << NTRIAL - nsuccess From c158ff355f019afa0885687dbd0563c0437b3e6a Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Thu, 25 Apr 2024 11:45:01 -0600 Subject: [PATCH 04/20] First clean up for better structure of tests. --- test/CMakeLists.txt | 6 +- test/pte_longtest_2phaseVinetSn.hpp | 64 ++++++----------- test/pte_test_2phaseVinetSn.hpp | 72 +++++++++----------- test/pte_test_5phaseSesameSn.hpp | 102 +++++++++++++++------------- test/pte_test_first.hpp | 65 ++++++++++++++++++ test/pte_test_utils.hpp | 44 +----------- test/test_pte.cpp | 1 + test/test_pte_2phase.cpp | 27 ++++---- test/test_pte_3phase.cpp | 40 ++++------- 9 files changed, 202 insertions(+), 219 deletions(-) create mode 100644 test/pte_test_first.hpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4a9a1442298..5a507fa45d5 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -98,15 +98,13 @@ if(SINGULARITY_BUILD_CLOSURE) singularity-eos::singularity-eos) target_include_directories(test_pte PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) add_test(pte test_pte) - endif() - if(SINGULARITY_USE_SPINER) +##### add_executable(test_pte_2phase test_pte_2phase.cpp) target_link_libraries(test_pte_2phase PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) target_include_directories(test_pte_2phase PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) add_test(pte_2phase test_pte_2phase) - endif() - if(SINGULARITY_USE_SPINER) +##### add_executable(test_pte_3phase test_pte_3phase.cpp) target_link_libraries(test_pte_3phase PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) diff --git a/test/pte_longtest_2phaseVinetSn.hpp b/test/pte_longtest_2phaseVinetSn.hpp index 78282bdbf05..67aa63b29b0 100644 --- a/test/pte_longtest_2phaseVinetSn.hpp +++ b/test/pte_longtest_2phaseVinetSn.hpp @@ -12,8 +12,8 @@ // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ -#ifndef _SINGULARITY_EOS_TEST_PTE_TEST_2PHASEVINETSN_ -#define _SINGULARITY_EOS_TEST_PTE_TEST_2PHASEVINETSN_ +#ifndef _SINGULARITY_EOS_TEST_PTE_LONGTEST_2PHASEVINETSN_ +#define _SINGULARITY_EOS_TEST_PTE_LONGTEST_2PHASEVINETSN_ #include @@ -21,6 +21,7 @@ #include #include +namespace pte_longtest_2phaseVinetSn { constexpr int NMAT = 2; constexpr int NTRIAL = 20; constexpr int NPTS = NTRIAL * NMAT; @@ -35,12 +36,8 @@ constexpr Real in_sie_tot[NTRIAL] = { 1.43767572e9, 1.72535639e9, 2.07977629e9, 2.50588681e9, 3.00885704e9, 3.59408011e9, 4.2671792100000005e9, 5.03401308e9, 5.90068122e9, 6.87352878e9, 7.95915117e9, 9.16439846e9, 1.04963795e10, 1.19624656e10, 1.35702945e10}; -constexpr Real in_lambda0 = 0.500480901; -constexpr Real in_lambda1 = 0.499519099; - -constexpr Real trial_vfrac0 = 47.0 / 100.0; -constexpr Real trial_vfrac1 = 53.0 / 100.0; - +constexpr Real in_lambda[NMAT] = {0.500480901,0.499519099}; +constexpr Real trial_vfrac[NMAT] = {47.0 / 100.0, 53.0/100.0}; constexpr Real out_press[NTRIAL] = { -3.29164604e-6, 1.284232715e10, 2.753234765e10, 4.423652945000001e10, 6.313670939999999e10, 8.443021825e10, @@ -75,32 +72,6 @@ constexpr Real out_vfrac1[NTRIAL] = {0.5, 0.498282729, 0.49686481, 0.49 0.491417663, 0.491352403, 0.491336852, 0.49136598, 0.491435452, 0.491541508, 0.491680878, 0.491850697}; -template -class LinearIndexer { - public: - PORTABLE_FUNCTION LinearIndexer() = default; - LinearIndexer(const T &t) : data_(t) {} - PORTABLE_INLINE_FUNCTION - auto &operator[](const int i) const { return data_(i); } - - private: - T data_; -}; - -template -class Indexer2D { - public: - PORTABLE_FUNCTION Indexer2D() = default; - PORTABLE_FUNCTION Indexer2D(const int j, const T &t) : j_(j), data_(t) {} - Indexer2D(const int j, const T &&t) = delete; // prevents r-value binding - PORTABLE_INLINE_FUNCTION - auto &operator[](const int i) const { return data_(j_, i); } - - private: - const int j_; - const T &data_; -}; - template inline void set_eos(T *eos) { constexpr Real d2to40[39] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., @@ -119,17 +90,24 @@ inline void set_eos(T *eos) { } template -inline void set_state(RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, - RealIndexer &&temp, EOSIndexer &&eos) { +inline void set_trial_state(int n, RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, + EOSIndexer &&eos) { - vfrac[0] = 47.0 / 100.0; - vfrac[1] = 53.0 / 100.0; - rho[0] = 7.278 * 0.500480901 / vfrac[0]; - rho[1] = 7.278 * 0.499519099 / vfrac[1]; - sie[0] = 8.41323509e08; - sie[1] = 8.41323509e08; + Real vsum = 0.; + for (int i = 0; i < NMAT; i++) { + vfrac[i] = trial_vfrac[i]; + rho[i] = in_lambda[i]*in_rho_tot[n]/vfrac[i]; + // same sie in both phases gives sie_tot=sie + sie[i] = in_sie_tot[n]; + vsum += vfrac[i]; + } + + for (int i = 0; i < NMAT; i++) + vfrac[i] *= 1.0 / vsum; return; } -#endif // _SINGULARITY_EOS_TEST_PTE_TEST_2PHASEVINETSN_ +} // namespace pte_longtest_2phaseVinetSn + +#endif // _SINGULARITY_EOS_TEST_PTE_LONGTEST_2PHASEVINETSN_ diff --git a/test/pte_test_2phaseVinetSn.hpp b/test/pte_test_2phaseVinetSn.hpp index 316dbbb635b..580b8dda7ed 100644 --- a/test/pte_test_2phaseVinetSn.hpp +++ b/test/pte_test_2phaseVinetSn.hpp @@ -21,55 +21,28 @@ #include #include +namespace pte_test_2phaseVinetSn { + constexpr int NMAT = 2; constexpr int NTRIAL = 5; constexpr int NPTS = NTRIAL * NMAT; constexpr int HIST_SIZE = 10; -constexpr Real in_rho_tot[5] = {7.27800000, 7.27981950, 7.28163945, 7.28345986, +constexpr Real in_rho_tot[NTRIAL] = {7.27800000, 7.27981950, 7.28163945, 7.28345986, 7.28528073}; -constexpr Real in_sie_tot[5] = {8.41323509e08, 8.41325565e08, 8.41331734e08, +constexpr Real in_sie_tot[NTRIAL] = {8.41323509e08, 8.41325565e08, 8.41331734e08, 8.41342022e08, 8.41356432e08}; -constexpr Real in_lambda0 = 0.500480901; -constexpr Real in_lambda1 = 0.499519099; - -constexpr Real trial_vfrac0 = 47.0 / 100.0; -constexpr Real trial_vfrac1 = 53.0 / 100.0; +constexpr Real in_lambda[NMAT] = {0.500480901,0.499519099}; +constexpr Real trial_vfrac[NMAT] = {47.0 / 100.0, 53.0/100.0}; -constexpr Real out_press[5] = {-3.29164604e-6, 1.19722694e8, 2.3968114450000003e8, +constexpr Real out_press[NTRIAL] = {-3.29164604e-6, 1.19722694e8, 2.3968114450000003e8, 3.598077865e8, 4.80104422e8}; -constexpr Real out_temp[5] = {298., 298.192952, 298.38594450000005, 298.5790125, +constexpr Real out_temp[NTRIAL] = {298., 298.192952, 298.38594450000005, 298.5790125, 298.7721535}; -constexpr Real out_rho0[5] = {7.28500000, 7.28653963, 7.28808705, 7.28963516, 7.29118416}; -constexpr Real out_rho1[5] = {7.27100000, 7.27309885, 7.27519088, 7.27728316, 7.27937551}; -constexpr Real out_vfrac0[5] = {0.5, 0.500019325, 0.500038138, 0.500056927, 0.500075679}; -constexpr Real out_vfrac1[5] = {0.5, 0.499980675, 0.499961862, 0.499943073, 0.499924321}; - -template -class LinearIndexer { - public: - PORTABLE_FUNCTION LinearIndexer() = default; - LinearIndexer(const T &t) : data_(t) {} - PORTABLE_INLINE_FUNCTION - auto &operator[](const int i) const { return data_(i); } - - private: - T data_; -}; - -template -class Indexer2D { - public: - PORTABLE_FUNCTION Indexer2D() = default; - PORTABLE_FUNCTION Indexer2D(const int j, const T &t) : j_(j), data_(t) {} - Indexer2D(const int j, const T &&t) = delete; // prevents r-value binding - PORTABLE_INLINE_FUNCTION - auto &operator[](const int i) const { return data_(j_, i); } - - private: - const int j_; - const T &data_; -}; +constexpr Real out_rho0[NTRIAL] = {7.28500000, 7.28653963, 7.28808705, 7.28963516, 7.29118416}; +constexpr Real out_rho1[NTRIAL] = {7.27100000, 7.27309885, 7.27519088, 7.27728316, 7.27937551}; +constexpr Real out_vfrac0[NTRIAL] = {0.5, 0.500019325, 0.500038138, 0.500056927, 0.500075679}; +constexpr Real out_vfrac1[NTRIAL] = {0.5, 0.499980675, 0.499961862, 0.499943073, 0.499924321}; template inline void set_eos(T *eos) { @@ -88,4 +61,25 @@ inline void set_eos(T *eos) { return; } +template +inline void set_trial_state(int n, RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, + EOSIndexer &&eos) { + + Real vsum = 0.; + for (int i = 0; i < NMAT; i++) { + vfrac[i] = trial_vfrac[i]; + rho[i] = in_lambda[i]*in_rho_tot[n]/vfrac[i]; + // same sie in both phases gives sie_tot=sie + sie[i] = in_sie_tot[n]; + vsum += vfrac[i]; + } + + for (int i = 0; i < NMAT; i++) + vfrac[i] *= 1.0 / vsum; + + return; +} + +} // namespace pte_test_2phaseVinetSn + #endif // _SINGULARITY_EOS_TEST_PTE_TEST_2PHASEVINETSN_ diff --git a/test/pte_test_5phaseSesameSn.hpp b/test/pte_test_5phaseSesameSn.hpp index 615b09703d7..483686c895e 100644 --- a/test/pte_test_5phaseSesameSn.hpp +++ b/test/pte_test_5phaseSesameSn.hpp @@ -28,64 +28,34 @@ constexpr int NTRIAL = 5; constexpr int NPTS = NTRIAL * NMAT; constexpr int HIST_SIZE = 30; -constexpr Real in_rho_tot[5] = {8.41382588, 8.41550873, 8.41719191, 8.41887543, +constexpr Real in_rho_tot[NTRIAL] = {8.41382588, 8.41550873, 8.41719191, 8.41887543, 8.42055928}; -constexpr Real in_sie_tot[5] = {1.67602914e09, 1.67876139e09, 1.68149578e09, +constexpr Real in_sie_tot[NTRIAL] = {1.67602914e09, 1.67876139e09, 1.68149578e09, 1.68423198e09, 1.68696932e09}; -constexpr Real in_lambda0[5] = {0.663852654, 0.660579238, 0.656880511, 0.652218675, - 0.646232106}; -constexpr Real in_lambda1[5] = {0.336093276, 0.339262144, 0.342703726, 0.346800126, - 0.352135167}; -constexpr Real in_lambda2[5] = {0.000054070, 0.000158618, 0.000415763, 0.000982199, - 0.001632727}; - -constexpr Real trial_vfrac0 = 650.0 / 1000.0; -constexpr Real trial_vfrac1 = 349.5 / 1000.0; -constexpr Real trial_vfrac2 = 0.5 / 1000.0; - -constexpr Real out_press[5] = {1.0946046e11, 1.0959183e11, 1.0971241e11, 1.0980857e11, +constexpr Real in_lambda[NMAT][NTRIAL] = { {0.663852654, 0.660579238, 0.656880511, 0.652218675, + 0.646232106}, {0.336093276, 0.339262144, 0.342703726, 0.346800126, + 0.352135167}, {0.000054070, 0.000158618, 0.000415763, 0.000982199, + 0.001632727}}; + +constexpr Real trial_vfrac[NMAT] = {650.0 / 1000.0, 349.5 / 1000.0, 0.5 / 1000.0}; + +constexpr Real out_press[NTRIAL] = {1.0946046e11, 1.0959183e11, 1.0971241e11, 1.0980857e11, 1.0987166e11}; -constexpr Real out_temp[5] = {438.96969, 438.98046, 438.95507, 438.84886, 438.64322}; -constexpr Real out_rho0[5] = {8.35319483, 8.35425038, 8.35523024, 8.35603912, 8.35661359}; -constexpr Real out_rho1[5] = {8.53618806, 8.53734120, 8.53841145, 8.53929455, 8.53992118}; -constexpr Real out_rho2[5] = {8.53854993, 8.53970495, 8.54077710, 8.54166211, 8.54229063}; -constexpr Real out_sie0[5] = {1.54379257e09, 1.54520613e09, 1.54644403e09, 1.54728448e09, +constexpr Real out_temp[NTRIAL] = {438.96969, 438.98046, 438.95507, 438.84886, 438.64322}; +constexpr Real out_rho0[NTRIAL] = {8.35319483, 8.35425038, 8.35523024, 8.35603912, 8.35661359}; +constexpr Real out_rho1[NTRIAL] = {8.53618806, 8.53734120, 8.53841145, 8.53929455, 8.53992118}; +constexpr Real out_rho2[NTRIAL] = {8.53854993, 8.53970495, 8.54077710, 8.54166211, 8.54229063}; +constexpr Real out_sie0[NTRIAL] = {1.54379257e09, 1.54520613e09, 1.54644403e09, 1.54728448e09, 1.54760327e09}; -constexpr Real out_sie1[5] = {1.93716882e09, 1.93864981e09, 1.93994981e09, 1.94084038e09, +constexpr Real out_sie1[NTRIAL] = {1.93716882e09, 1.93864981e09, 1.93994981e09, 1.94084038e09, 1.94119298e09}; -constexpr Real out_sie2[5] = {2.01473728e09, 2.01621655e09, 2.01751522e09, 2.01840528e09, +constexpr Real out_sie2[NTRIAL] = {2.01473728e09, 2.01621655e09, 2.01751522e09, 2.01840528e09, 2.01875845e09}; -template -class LinearIndexer { - public: - PORTABLE_FUNCTION LinearIndexer() = default; - LinearIndexer(const T &t) : data_(t) {} - PORTABLE_INLINE_FUNCTION - auto &operator[](const int i) const { return data_(i); } - - private: - T data_; -}; - -template -class Indexer2D { - public: - PORTABLE_FUNCTION Indexer2D() = default; - PORTABLE_FUNCTION Indexer2D(const int j, const T &t) : j_(j), data_(t) {} - Indexer2D(const int j, const T &&t) = delete; // prevents r-value binding - PORTABLE_INLINE_FUNCTION - auto &operator[](const int i) const { return data_(j_, i); } - - private: - const int j_; - const T &data_; -}; - template inline void set_eos(T *eos) { - symlink("/projects/shavano/dev/aematts/SEOS/PTEsolver_test_cases/sn2162-v01.bin", + int sl = symlink("/projects/shavano/dev/aematts/SEOS/PTEsolver_test_cases/sn2162-v01.bin", "sesameu"); int SnbetaID = 2102; @@ -110,4 +80,40 @@ inline void set_eos(T *eos) { return; } + //vfrac_hm(n, 0) = trial_vfrac0; + //vfrac_hm(n, 1) = trial_vfrac1; + // vfrac_hm(n,x) = 0.0; + //vfrac_hm(n, 2) = trial_vfrac2; + // vfrac_hm(n,x) = 0.0; + //rho_hm(n, 0) = in_lambda0[n] * in_rho_tot[n] / vfrac_hm(n, 0); + //rho_hm(n, 1) = in_lambda1[n] * in_rho_tot[n] / vfrac_hm(n, 1); + // rho_hm(n,x) = 7.380; + //rho_hm(n, 2) = in_lambda2[n] * in_rho_tot[n] / vfrac_hm(n, 2); + // rho_hm(n,x) = 7.116; + // same sie in both phases gives sie_tot=sie below + //sie_hm(n, 0) = in_sie_tot[n]; + //sie_hm(n, 1) = in_sie_tot[n]; + // sie_hm(n,x) = in_sie_tot[n]; + //sie_hm(n, 2) = in_sie_tot[n]; + // sie_hm(n,x) = in_sie_tot[n]; + +template +inline void set_trial_3state(int n, RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, + EOSIndexer &&eos) { + + Real vsum = 0.; + for (int i = 0; i < NMAT; i++) { + vfrac[i] = trial_vfrac[i]; + rho[i] = in_lambda[i][n]*in_rho_tot[n]/vfrac[i]; + // same sie in both phases gives sie_tot=sie + sie[i] = in_sie_tot[n]; + vsum += vfrac[i]; + } + + for (int i = 0; i < NMAT; i++) + vfrac[i] *= 1.0 / vsum; + + return; +} + #endif // _SINGULARITY_EOS_TEST_PTE_TEST_5PHASESESAMESN_ diff --git a/test/pte_test_first.hpp b/test/pte_test_first.hpp new file mode 100644 index 00000000000..d81795aeb11 --- /dev/null +++ b/test/pte_test_first.hpp @@ -0,0 +1,65 @@ +//------------------------------------------------------------------------------ +// © 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. +//------------------------------------------------------------------------------ + +#ifndef _SINGULARITY_EOS_TEST_PTE_TEST_FIRST_ +#define _SINGULARITY_EOS_TEST_PTE_TEST_FIRST_ + +#include + +#include +#include +#include + +constexpr int NMAT = 3; +constexpr int NTRIAL = 100; +constexpr int NPTS = NTRIAL * NMAT; +constexpr int HIST_SIZE = 10; + +template +inline void set_eos(T *eos) { + singularity::EOS gr = singularity::Gruneisen(394000.0, 1.489, 0.0, 0.0, 2.02, 0.47, + 8.93, 297.0, 1.0e6, 0.383e7); + singularity::EOS dr = singularity::DavisReactants( + 1.890, 4.115e10, 1.0e6, 297.0, 1.8e5, 4.6, 0.34, 0.56, 0.0, 0.4265, 0.001074e10); + singularity::EOS dp = singularity::DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419, + 3.2e10, 0.001072e10, 0.0); + eos[0] = gr.GetOnDevice(); + eos[1] = dr.GetOnDevice(); + eos[2] = dp.GetOnDevice(); + return; +} + +template +inline void set_state(RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, + RealIndexer &&temp, EOSIndexer &&eos) { + + rho[0] = 8.93; + rho[1] = 1.89; + rho[2] = 2.5; + + Real vsum = 0.; + for (int i = 0; i < NMAT; i++) { + temp[i] = 600.0; + sie[i] = eos[i].InternalEnergyFromDensityTemperature(rho[i], temp[i]); + vfrac[i] = rand() / (1.0 * RAND_MAX); + vsum += vfrac[i]; + } + + for (int i = 0; i < NMAT; i++) + vfrac[i] *= 1.0 / vsum; + + return; +} + +#endif // _SINGULARITY_EOS_TEST_PTE_TEST_FIRST_ diff --git a/test/pte_test_utils.hpp b/test/pte_test_utils.hpp index d61b09747e3..d997ceb410b 100644 --- a/test/pte_test_utils.hpp +++ b/test/pte_test_utils.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2024. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National @@ -19,12 +19,6 @@ #include #include -#include - -constexpr int NMAT = 3; -constexpr int NTRIAL = 100; -constexpr int NPTS = NTRIAL * NMAT; -constexpr int HIST_SIZE = 10; template class LinearIndexer { @@ -52,40 +46,4 @@ class Indexer2D { const T &data_; }; -template -inline void set_eos(T *eos) { - singularity::EOS gr = singularity::Gruneisen(394000.0, 1.489, 0.0, 0.0, 2.02, 0.47, - 8.93, 297.0, 1.0e6, 0.383e7); - singularity::EOS dr = singularity::DavisReactants( - 1.890, 4.115e10, 1.0e6, 297.0, 1.8e5, 4.6, 0.34, 0.56, 0.0, 0.4265, 0.001074e10); - singularity::EOS dp = singularity::DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419, - 3.2e10, 0.001072e10, 0.0); - eos[0] = gr.GetOnDevice(); - eos[1] = dr.GetOnDevice(); - eos[2] = dp.GetOnDevice(); - return; -} - -template -inline void set_state(RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, - RealIndexer &&temp, EOSIndexer &&eos) { - - rho[0] = 8.93; - rho[1] = 1.89; - rho[2] = 2.5; - - Real vsum = 0.; - for (int i = 0; i < NMAT; i++) { - temp[i] = 600.0; - sie[i] = eos[i].InternalEnergyFromDensityTemperature(rho[i], temp[i]); - vfrac[i] = rand() / (1.0 * RAND_MAX); - vsum += vfrac[i]; - } - - for (int i = 0; i < NMAT; i++) - vfrac[i] *= 1.0 / vsum; - - return; -} - #endif // _SINGULARITY_EOS_TEST_PTE_TEST_UTILS_ diff --git a/test/test_pte.cpp b/test/test_pte.cpp index bfed57224f3..c69dad20b3c 100644 --- a/test/test_pte.cpp +++ b/test/test_pte.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include diff --git a/test/test_pte_2phase.cpp b/test/test_pte_2phase.cpp index d7be8bb5a95..ed6f5f6cb43 100644 --- a/test/test_pte_2phase.cpp +++ b/test/test_pte_2phase.cpp @@ -22,13 +22,16 @@ #include #include #include -// #include +#include +#include #include #include #include using namespace singularity; +using namespace pte_longtest_2phaseVinetSn; + using DataBox = Spiner::DataBox; int main(int argc, char *argv[]) { @@ -49,9 +52,7 @@ int main(int argc, char *argv[]) { #endif LinearIndexer eos_h(eos_hv); - std::cout << "Before set_eos" << std::endl; set_eos(eos_hv.data()); - std::cout << "After set_eos" << std::endl; #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::deep_copy(eos_v, eos_hv); @@ -113,14 +114,10 @@ int main(int argc, char *argv[]) { // setup state srand(time(NULL)); for (int n = 0; n < NTRIAL; n++) { - - vfrac_hm(n, 0) = trial_vfrac0; - vfrac_hm(n, 1) = trial_vfrac1; - rho_hm(n, 0) = in_lambda0 * in_rho_tot[n] / vfrac_hm(n, 0); - rho_hm(n, 1) = in_lambda1 * in_rho_tot[n] / vfrac_hm(n, 1); - // same sie in both phases gives sie_tot=sie below - sie_hm(n, 0) = in_sie_tot[n]; - sie_hm(n, 1) = in_sie_tot[n]; + Indexer2D r(n, rho_hm); + Indexer2D vf(n, vfrac_hm); + Indexer2D e(n, sie_hm); + set_trial_state(n, r, vf, e, eos_h); } for (int i = 0; i < HIST_SIZE; ++i) { hist_vh[i] = 0; @@ -151,8 +148,8 @@ int main(int argc, char *argv[]) { std::cout << "Trial number: " << n << std::endl; std::cout << "Total Specific Internal energy: \t\t" << in_sie_tot[n] << std::endl; std::cout << "Total density: \t\t\t\t\t" << in_rho_tot[n] << std::endl; - std::cout << "Mass fractions: beta, gamma: \t\t\t" << in_lambda0 << ", " - << in_lambda1 << std::endl; + std::cout << "Mass fractions: beta, gamma: \t\t\t" << in_lambda[0] << ", " + << in_lambda[1] << std::endl; std::cout << "Assuming volume fractions: beta, gamma: \t" << vfrac_hm(n, 0) << ", " << vfrac_hm(n, 1) << std::endl; std::cout << "gives starting phase densities: beta, gamma: \t" << rho_hm(n, 0) @@ -219,11 +216,11 @@ int main(int argc, char *argv[]) { for (int n = 0; n < NTRIAL; n++) { std::cout << "Trial number: " << n << std::endl; std::cout << "Total Specific Internal energy: \t" - << sie_hm(n, 0) * in_lambda0 + sie_hm(n, 1) * in_lambda1 << ", (" + << sie_hm(n, 0) * in_lambda[0] + sie_hm(n, 1) * in_lambda[1] << ", (" << in_sie_tot[n] << ")" << std::endl; std::cout << "Total density: \t\t\t\t" << 1.0 / - (1.0 / rho_hm(n, 0) * in_lambda0 + 1.0 / rho_hm(n, 1) * in_lambda1) + (1.0 / rho_hm(n, 0) * in_lambda[0] + 1.0 / rho_hm(n, 1) * in_lambda[1]) << ", (" << in_rho_tot[n] << ")" << std::endl; std::cout << "Volume fractions: beta, gamma: \t\t" << vfrac_hm(n, 0) << ", " << vfrac_hm(n, 1) << std::endl; diff --git a/test/test_pte_3phase.cpp b/test/test_pte_3phase.cpp index a44a8dc90ba..583b02f5d31 100644 --- a/test/test_pte_3phase.cpp +++ b/test/test_pte_3phase.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -48,9 +49,7 @@ int main(int argc, char *argv[]) { #endif LinearIndexer eos_h(eos_hv); - std::cout << "Before set_eos" << std::endl; set_eos(eos_hv.data()); - std::cout << "After set_eos" << std::endl; #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::deep_copy(eos_v, eos_hv); @@ -109,26 +108,13 @@ int main(int argc, char *argv[]) { int *hist_d = hist_vh; #endif - // setup state + // setup state srand(time(NULL)); for (int n = 0; n < NTRIAL; n++) { - - vfrac_hm(n, 0) = trial_vfrac0; - vfrac_hm(n, 1) = trial_vfrac1; - // vfrac_hm(n,x) = 0.0; - vfrac_hm(n, 2) = trial_vfrac2; - // vfrac_hm(n,x) = 0.0; - rho_hm(n, 0) = in_lambda0[n] * in_rho_tot[n] / vfrac_hm(n, 0); - rho_hm(n, 1) = in_lambda1[n] * in_rho_tot[n] / vfrac_hm(n, 1); - // rho_hm(n,x) = 7.380; - rho_hm(n, 2) = in_lambda2[n] * in_rho_tot[n] / vfrac_hm(n, 2); - // rho_hm(n,x) = 7.116; - // same sie in both phases gives sie_tot=sie below - sie_hm(n, 0) = in_sie_tot[n]; - sie_hm(n, 1) = in_sie_tot[n]; - // sie_hm(n,x) = in_sie_tot[n]; - sie_hm(n, 2) = in_sie_tot[n]; - // sie_hm(n,x) = in_sie_tot[n]; + Indexer2D r(n, rho_hm); + Indexer2D vf(n, vfrac_hm); + Indexer2D e(n, sie_hm); + set_trial_3state(n, r, vf, e, eos_h); } for (int i = 0; i < HIST_SIZE; ++i) { hist_vh[i] = 0; @@ -159,8 +145,8 @@ int main(int argc, char *argv[]) { std::cout << "Trial number: " << n << std::endl; std::cout << "Total Specific Internal energy: \t\t\t" << in_sie_tot[n] << std::endl; std::cout << "Total density: \t\t\t\t\t\t" << in_rho_tot[n] << std::endl; - std::cout << "Mass fractions: beta, gamma, hcp: \t\t\t" << in_lambda0[n] << ", " - << in_lambda1[n] << ", " << in_lambda2[n] << std::endl; + std::cout << "Mass fractions: beta, gamma, hcp: \t\t\t" << in_lambda[0][n] << ", " + << in_lambda[1][n] << ", " << in_lambda[2][n] << std::endl; std::cout << "Assuming volume fractions: beta, gamma, hcp: \t\t" << vfrac_hm(n, 0) << ", " << vfrac_hm(n, 1) << ", " << vfrac_hm(n, 2) << std::endl; std::cout << "gives starting phase densities: beta, gamma, hcp: \t" << rho_hm(n, 0) @@ -227,13 +213,13 @@ int main(int argc, char *argv[]) { for (int n = 0; n < NTRIAL; n++) { std::cout << "Trial number: " << n << std::endl; std::cout << "Total Specific Internal energy: \t" - << sie_hm(n, 0) * in_lambda0[n] + sie_hm(n, 1) * in_lambda1[n] + - sie_hm(n, 2) * in_lambda2[n] + << sie_hm(n, 0) * in_lambda[0][n] + sie_hm(n, 1) * in_lambda[1][n] + + sie_hm(n, 2) * in_lambda[2][n] << ", (" << in_sie_tot[n] << ")" << std::endl; std::cout << "Total density: \t\t\t\t" - << 1.0 / (1.0 / rho_hm(n, 0) * in_lambda0[n] + - 1.0 / rho_hm(n, 1) * in_lambda1[n] + - 1.0 / rho_hm(n, 2) * in_lambda2[n]) + << 1.0 / (1.0 / rho_hm(n, 0) * in_lambda[0][n] + + 1.0 / rho_hm(n, 1) * in_lambda[1][n] + + 1.0 / rho_hm(n, 2) * in_lambda[2][n]) << ", (" << in_rho_tot[n] << ")" << std::endl; std::cout << "Volume fractions: beta, gamma, hcp: \t" << vfrac_hm(n, 0) << ", " << vfrac_hm(n, 1) << ", " << vfrac_hm(n, 2) << std::endl; From c285ad1a4278183247b6657000ece464548a3dbe Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Thu, 25 Apr 2024 13:08:32 -0600 Subject: [PATCH 05/20] Cleanup of the first cleanup..... --- test/CMakeLists.txt | 14 +++--- test/pte_longtest_2phaseVinetSn.hpp | 14 ++++-- test/pte_test_2phaseVinetSn.hpp | 36 +++++++++------ test/pte_test_5phaseSesameSn.hpp | 71 +++++++++++++---------------- test/test_pte.cpp | 2 +- test/test_pte_2phase.cpp | 4 +- test/test_pte_3phase.cpp | 2 +- 7 files changed, 75 insertions(+), 68 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5a507fa45d5..3290fc4a941 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -104,12 +104,14 @@ if(SINGULARITY_BUILD_CLOSURE) singularity-eos::singularity-eos) target_include_directories(test_pte_2phase PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) add_test(pte_2phase test_pte_2phase) -##### - add_executable(test_pte_3phase test_pte_3phase.cpp) - target_link_libraries(test_pte_3phase PRIVATE Catch2::Catch2 - singularity-eos::singularity-eos) - target_include_directories(test_pte_3phase PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) - add_test(pte_3phase test_pte_3phase) + + if(SINGULARITY_USE_EOSPAC) + add_executable(test_pte_3phase test_pte_3phase.cpp) + target_link_libraries(test_pte_3phase PRIVATE Catch2::Catch2 + singularity-eos::singularity-eos) + target_include_directories(test_pte_3phase PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) + add_test(pte_3phase test_pte_3phase) + endif() endif() if(SINGULARITY_USE_KOKKOS AND SINGULARITY_USE_FORTRAN) add_executable(test_get_sg_eos test_get_sg_eos.cpp) diff --git a/test/pte_longtest_2phaseVinetSn.hpp b/test/pte_longtest_2phaseVinetSn.hpp index 67aa63b29b0..6aa78c7371f 100644 --- a/test/pte_longtest_2phaseVinetSn.hpp +++ b/test/pte_longtest_2phaseVinetSn.hpp @@ -21,6 +21,10 @@ #include #include +// This data is taken from the .ult file for a test case in Flag, MPMat_KPT_PTsolv.ult.std +// Note that at this point in time (April 2024) the global and phase times are not alinged +// in .ult but they will be in the near future. + namespace pte_longtest_2phaseVinetSn { constexpr int NMAT = 2; constexpr int NTRIAL = 20; @@ -36,8 +40,8 @@ constexpr Real in_sie_tot[NTRIAL] = { 1.43767572e9, 1.72535639e9, 2.07977629e9, 2.50588681e9, 3.00885704e9, 3.59408011e9, 4.2671792100000005e9, 5.03401308e9, 5.90068122e9, 6.87352878e9, 7.95915117e9, 9.16439846e9, 1.04963795e10, 1.19624656e10, 1.35702945e10}; -constexpr Real in_lambda[NMAT] = {0.500480901,0.499519099}; -constexpr Real trial_vfrac[NMAT] = {47.0 / 100.0, 53.0/100.0}; +constexpr Real in_lambda[NMAT] = {0.500480901, 0.499519099}; +constexpr Real trial_vfrac[NMAT] = {47.0 / 100.0, 53.0 / 100.0}; constexpr Real out_press[NTRIAL] = { -3.29164604e-6, 1.284232715e10, 2.753234765e10, 4.423652945000001e10, 6.313670939999999e10, 8.443021825e10, @@ -90,13 +94,13 @@ inline void set_eos(T *eos) { } template -inline void set_trial_state(int n, RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, - EOSIndexer &&eos) { +inline void set_trial_state(int n, RealIndexer &&rho, RealIndexer &&vfrac, + RealIndexer &&sie, EOSIndexer &&eos) { Real vsum = 0.; for (int i = 0; i < NMAT; i++) { vfrac[i] = trial_vfrac[i]; - rho[i] = in_lambda[i]*in_rho_tot[n]/vfrac[i]; + rho[i] = in_lambda[i] * in_rho_tot[n] / vfrac[i]; // same sie in both phases gives sie_tot=sie sie[i] = in_sie_tot[n]; vsum += vfrac[i]; diff --git a/test/pte_test_2phaseVinetSn.hpp b/test/pte_test_2phaseVinetSn.hpp index 580b8dda7ed..595620979ea 100644 --- a/test/pte_test_2phaseVinetSn.hpp +++ b/test/pte_test_2phaseVinetSn.hpp @@ -21,6 +21,10 @@ #include #include +// This data is taken from the .ult file for a test case in Flag, MPMat_KPT_PTsolv.ult.std +// Note that at this point in time (April 2024) the global and phase times are not alinged +// in .ult but they will be in the near future. + namespace pte_test_2phaseVinetSn { constexpr int NMAT = 2; @@ -29,20 +33,24 @@ constexpr int NPTS = NTRIAL * NMAT; constexpr int HIST_SIZE = 10; constexpr Real in_rho_tot[NTRIAL] = {7.27800000, 7.27981950, 7.28163945, 7.28345986, - 7.28528073}; + 7.28528073}; constexpr Real in_sie_tot[NTRIAL] = {8.41323509e08, 8.41325565e08, 8.41331734e08, - 8.41342022e08, 8.41356432e08}; -constexpr Real in_lambda[NMAT] = {0.500480901,0.499519099}; -constexpr Real trial_vfrac[NMAT] = {47.0 / 100.0, 53.0/100.0}; + 8.41342022e08, 8.41356432e08}; +constexpr Real in_lambda[NMAT] = {0.500480901, 0.499519099}; +constexpr Real trial_vfrac[NMAT] = {47.0 / 100.0, 53.0 / 100.0}; constexpr Real out_press[NTRIAL] = {-3.29164604e-6, 1.19722694e8, 2.3968114450000003e8, - 3.598077865e8, 4.80104422e8}; + 3.598077865e8, 4.80104422e8}; constexpr Real out_temp[NTRIAL] = {298., 298.192952, 298.38594450000005, 298.5790125, - 298.7721535}; -constexpr Real out_rho0[NTRIAL] = {7.28500000, 7.28653963, 7.28808705, 7.28963516, 7.29118416}; -constexpr Real out_rho1[NTRIAL] = {7.27100000, 7.27309885, 7.27519088, 7.27728316, 7.27937551}; -constexpr Real out_vfrac0[NTRIAL] = {0.5, 0.500019325, 0.500038138, 0.500056927, 0.500075679}; -constexpr Real out_vfrac1[NTRIAL] = {0.5, 0.499980675, 0.499961862, 0.499943073, 0.499924321}; + 298.7721535}; +constexpr Real out_rho0[NTRIAL] = {7.28500000, 7.28653963, 7.28808705, 7.28963516, + 7.29118416}; +constexpr Real out_rho1[NTRIAL] = {7.27100000, 7.27309885, 7.27519088, 7.27728316, + 7.27937551}; +constexpr Real out_vfrac0[NTRIAL] = {0.5, 0.500019325, 0.500038138, 0.500056927, + 0.500075679}; +constexpr Real out_vfrac1[NTRIAL] = {0.5, 0.499980675, 0.499961862, 0.499943073, + 0.499924321}; template inline void set_eos(T *eos) { @@ -62,14 +70,14 @@ inline void set_eos(T *eos) { } template -inline void set_trial_state(int n, RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, - EOSIndexer &&eos) { +inline void set_trial_state(int n, RealIndexer &&rho, RealIndexer &&vfrac, + RealIndexer &&sie, EOSIndexer &&eos) { Real vsum = 0.; for (int i = 0; i < NMAT; i++) { vfrac[i] = trial_vfrac[i]; - rho[i] = in_lambda[i]*in_rho_tot[n]/vfrac[i]; - // same sie in both phases gives sie_tot=sie + rho[i] = in_lambda[i] * in_rho_tot[n] / vfrac[i]; + // same sie in both phases gives sie_tot=sie sie[i] = in_sie_tot[n]; vsum += vfrac[i]; } diff --git a/test/pte_test_5phaseSesameSn.hpp b/test/pte_test_5phaseSesameSn.hpp index 483686c895e..75664066b03 100644 --- a/test/pte_test_5phaseSesameSn.hpp +++ b/test/pte_test_5phaseSesameSn.hpp @@ -23,40 +23,50 @@ #include #include +// This data is taken from the .ult file for a test case in Flag, +// MPMat_KPT_3phases.ult.std in MPMat_KPT_3phases_NW. Note that 3 phases appear by mistake +// (no windows) in sn2162 where only 2 phases should be present (with windows). Note that +// at this point in time (April 2024) the global and phase times are not alinged in .ult +// but they will be in the near future. + constexpr int NMAT = 3; constexpr int NTRIAL = 5; constexpr int NPTS = NTRIAL * NMAT; constexpr int HIST_SIZE = 30; constexpr Real in_rho_tot[NTRIAL] = {8.41382588, 8.41550873, 8.41719191, 8.41887543, - 8.42055928}; + 8.42055928}; constexpr Real in_sie_tot[NTRIAL] = {1.67602914e09, 1.67876139e09, 1.68149578e09, - 1.68423198e09, 1.68696932e09}; -constexpr Real in_lambda[NMAT][NTRIAL] = { {0.663852654, 0.660579238, 0.656880511, 0.652218675, - 0.646232106}, {0.336093276, 0.339262144, 0.342703726, 0.346800126, - 0.352135167}, {0.000054070, 0.000158618, 0.000415763, 0.000982199, - 0.001632727}}; + 1.68423198e09, 1.68696932e09}; +constexpr Real in_lambda[NMAT][NTRIAL] = { + {0.663852654, 0.660579238, 0.656880511, 0.652218675, 0.646232106}, + {0.336093276, 0.339262144, 0.342703726, 0.346800126, 0.352135167}, + {0.000054070, 0.000158618, 0.000415763, 0.000982199, 0.001632727}}; constexpr Real trial_vfrac[NMAT] = {650.0 / 1000.0, 349.5 / 1000.0, 0.5 / 1000.0}; -constexpr Real out_press[NTRIAL] = {1.0946046e11, 1.0959183e11, 1.0971241e11, 1.0980857e11, - 1.0987166e11}; +constexpr Real out_press[NTRIAL] = {1.0946046e11, 1.0959183e11, 1.0971241e11, + 1.0980857e11, 1.0987166e11}; constexpr Real out_temp[NTRIAL] = {438.96969, 438.98046, 438.95507, 438.84886, 438.64322}; -constexpr Real out_rho0[NTRIAL] = {8.35319483, 8.35425038, 8.35523024, 8.35603912, 8.35661359}; -constexpr Real out_rho1[NTRIAL] = {8.53618806, 8.53734120, 8.53841145, 8.53929455, 8.53992118}; -constexpr Real out_rho2[NTRIAL] = {8.53854993, 8.53970495, 8.54077710, 8.54166211, 8.54229063}; -constexpr Real out_sie0[NTRIAL] = {1.54379257e09, 1.54520613e09, 1.54644403e09, 1.54728448e09, - 1.54760327e09}; -constexpr Real out_sie1[NTRIAL] = {1.93716882e09, 1.93864981e09, 1.93994981e09, 1.94084038e09, - 1.94119298e09}; -constexpr Real out_sie2[NTRIAL] = {2.01473728e09, 2.01621655e09, 2.01751522e09, 2.01840528e09, - 2.01875845e09}; +constexpr Real out_rho0[NTRIAL] = {8.35319483, 8.35425038, 8.35523024, 8.35603912, + 8.35661359}; +constexpr Real out_rho1[NTRIAL] = {8.53618806, 8.53734120, 8.53841145, 8.53929455, + 8.53992118}; +constexpr Real out_rho2[NTRIAL] = {8.53854993, 8.53970495, 8.54077710, 8.54166211, + 8.54229063}; +constexpr Real out_sie0[NTRIAL] = {1.54379257e09, 1.54520613e09, 1.54644403e09, + 1.54728448e09, 1.54760327e09}; +constexpr Real out_sie1[NTRIAL] = {1.93716882e09, 1.93864981e09, 1.93994981e09, + 1.94084038e09, 1.94119298e09}; +constexpr Real out_sie2[NTRIAL] = {2.01473728e09, 2.01621655e09, 2.01751522e09, + 2.01840528e09, 2.01875845e09}; template inline void set_eos(T *eos) { - int sl = symlink("/projects/shavano/dev/aematts/SEOS/PTEsolver_test_cases/sn2162-v01.bin", - "sesameu"); + int sl = + symlink("/projects/shavano/dev/aematts/SEOS/PTEsolver_test_cases/sn2162-v01.bin", + "sesameu"); int SnbetaID = 2102; int SngammaID = 2103; @@ -80,31 +90,14 @@ inline void set_eos(T *eos) { return; } - //vfrac_hm(n, 0) = trial_vfrac0; - //vfrac_hm(n, 1) = trial_vfrac1; - // vfrac_hm(n,x) = 0.0; - //vfrac_hm(n, 2) = trial_vfrac2; - // vfrac_hm(n,x) = 0.0; - //rho_hm(n, 0) = in_lambda0[n] * in_rho_tot[n] / vfrac_hm(n, 0); - //rho_hm(n, 1) = in_lambda1[n] * in_rho_tot[n] / vfrac_hm(n, 1); - // rho_hm(n,x) = 7.380; - //rho_hm(n, 2) = in_lambda2[n] * in_rho_tot[n] / vfrac_hm(n, 2); - // rho_hm(n,x) = 7.116; - // same sie in both phases gives sie_tot=sie below - //sie_hm(n, 0) = in_sie_tot[n]; - //sie_hm(n, 1) = in_sie_tot[n]; - // sie_hm(n,x) = in_sie_tot[n]; - //sie_hm(n, 2) = in_sie_tot[n]; - // sie_hm(n,x) = in_sie_tot[n]; - template -inline void set_trial_3state(int n, RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, - EOSIndexer &&eos) { +inline void set_trial_3state(int n, RealIndexer &&rho, RealIndexer &&vfrac, + RealIndexer &&sie, EOSIndexer &&eos) { Real vsum = 0.; for (int i = 0; i < NMAT; i++) { vfrac[i] = trial_vfrac[i]; - rho[i] = in_lambda[i][n]*in_rho_tot[n]/vfrac[i]; + rho[i] = in_lambda[i][n] * in_rho_tot[n] / vfrac[i]; // same sie in both phases gives sie_tot=sie sie[i] = in_sie_tot[n]; vsum += vfrac[i]; diff --git a/test/test_pte.cpp b/test/test_pte.cpp index c69dad20b3c..44eb859241d 100644 --- a/test/test_pte.cpp +++ b/test/test_pte.cpp @@ -21,8 +21,8 @@ #include #include -#include #include +#include #include #include #include diff --git a/test/test_pte_2phase.cpp b/test/test_pte_2phase.cpp index ed6f5f6cb43..9fa55a6007b 100644 --- a/test/test_pte_2phase.cpp +++ b/test/test_pte_2phase.cpp @@ -219,8 +219,8 @@ int main(int argc, char *argv[]) { << sie_hm(n, 0) * in_lambda[0] + sie_hm(n, 1) * in_lambda[1] << ", (" << in_sie_tot[n] << ")" << std::endl; std::cout << "Total density: \t\t\t\t" - << 1.0 / - (1.0 / rho_hm(n, 0) * in_lambda[0] + 1.0 / rho_hm(n, 1) * in_lambda[1]) + << 1.0 / (1.0 / rho_hm(n, 0) * in_lambda[0] + + 1.0 / rho_hm(n, 1) * in_lambda[1]) << ", (" << in_rho_tot[n] << ")" << std::endl; std::cout << "Volume fractions: beta, gamma: \t\t" << vfrac_hm(n, 0) << ", " << vfrac_hm(n, 1) << std::endl; diff --git a/test/test_pte_3phase.cpp b/test/test_pte_3phase.cpp index 583b02f5d31..a605985d10c 100644 --- a/test/test_pte_3phase.cpp +++ b/test/test_pte_3phase.cpp @@ -108,7 +108,7 @@ int main(int argc, char *argv[]) { int *hist_d = hist_vh; #endif - // setup state + // setup state srand(time(NULL)); for (int n = 0; n < NTRIAL; n++) { Indexer2D r(n, rho_hm); From 9825064eb38c6127a3fa76a0fc098dd83c77a9ca Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Thu, 1 Aug 2024 17:28:38 -0600 Subject: [PATCH 06/20] Fixed an shortened argumentlist. Implemented Carl Greeff's rate equation for mass fraction update and added a test case. This is the first ingredient of the Kinetic PhaseTransiton (KPT) scheme. --- singularity-eos/CMakeLists.txt | 1 + .../kinetic_phasetransition_models.hpp | 62 +++++++++++++++++++ test/CMakeLists.txt | 1 + test/pte_test_first.hpp | 2 +- 4 files changed, 65 insertions(+), 1 deletion(-) create mode 100644 singularity-eos/closure/kinetic_phasetransition_models.hpp diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 399702f989d..34c576534f5 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -62,6 +62,7 @@ register_headers( eos/singularity_eos_init_utils.hpp eos/variant_utils.hpp eos/eos_carnahan_starling.hpp + closure/kinetic_phasetransition_models.hpp ) if (SINGULARITY_BUILD_CLOSURE) diff --git a/singularity-eos/closure/kinetic_phasetransition_models.hpp b/singularity-eos/closure/kinetic_phasetransition_models.hpp new file mode 100644 index 00000000000..161deebd140 --- /dev/null +++ b/singularity-eos/closure/kinetic_phasetransition_models.hpp @@ -0,0 +1,62 @@ +//------------------------------------------------------------------------------ +// © 2024. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#ifndef _SINGULARITY_EOS_CLOSURE_KINETIC_PHSETRANSITION_MODELS_ +#define _SINGULARITY_EOS_CLOSURE_KINETIC_PHSETRANSITION_MODELS_ + +#include +#include +#include +#include +#include + +#include + +#ifdef SINGULARITY_USE_KOKKOSKERNELS +#include +#include +#include +#else +#include +#endif // SINGULARITY_USE_KOKKOSKERNELS + +namespace singularity { + + +PORTABLE_INLINE_FUNCTION void LogRatesCGModel(const Real *w, const Real *b, + const int num_phases, const Real *gibbs, const int *gibbsorder, Real *logRjk, int *fromto) { + +// begin the calculation of logRik + int jk = 0; + for (int j = 0; j < num_phases-1; j++){ + // from high Gibb phase to low + for (int k = num_phases-1; k > j; k--){ + // from low Gibb phase to high, k>j always, gibbs_j > gibbs_k always + // always mass transport from j->k + Real hgj=gibbs[gibbsorder[j]]; + Real hgk=gibbs[gibbsorder[k]]; +// calc phase dx from hgi to hgk + Real dgdb=(hgj-hgk)/(b[gibbsorder[j]*num_phases + gibbsorder[k]]); + // First is highest level to levels below, largest gibbs energy diff first. Then follows 2nd highest to levels below. + logRjk[jk] = log(w[gibbsorder[j]*num_phases + gibbsorder[k]]*dgdb)+dgdb*dgdb; + fromto[jk] = gibbsorder[j]*10+gibbsorder[k]; + jk++; + } + } + return; +} + +} // namespace singularity + +#endif // _SINGULARITY_EOS_CLOSURE_KINETIC_PHSETRANSITION_MODELS_ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3a7c3ef3242..6d70f552c8d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -49,6 +49,7 @@ add_executable( closure_unit_tests catch2_define.cpp test_closure_pte.cpp + test_kpt_models.cpp ) get_property(plugin_tests GLOBAL PROPERTY PLUGIN_TESTS) diff --git a/test/pte_test_first.hpp b/test/pte_test_first.hpp index d81795aeb11..77d6dbb6901 100644 --- a/test/pte_test_first.hpp +++ b/test/pte_test_first.hpp @@ -33,7 +33,7 @@ inline void set_eos(T *eos) { singularity::EOS dr = singularity::DavisReactants( 1.890, 4.115e10, 1.0e6, 297.0, 1.8e5, 4.6, 0.34, 0.56, 0.0, 0.4265, 0.001074e10); singularity::EOS dp = singularity::DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419, - 3.2e10, 0.001072e10, 0.0); + 3.2e10, 0.001072e10); eos[0] = gr.GetOnDevice(); eos[1] = dr.GetOnDevice(); eos[2] = dp.GetOnDevice(); From c9a18214f650ae3bbe0281f77323f00b8a76c04e Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Thu, 1 Aug 2024 17:43:34 -0600 Subject: [PATCH 07/20] Added forgotten file and formatted. --- .../kinetic_phasetransition_models.hpp | 43 +++--- test/test_kpt_models.cpp | 141 ++++++++++++++++++ 2 files changed, 164 insertions(+), 20 deletions(-) create mode 100644 test/test_kpt_models.cpp diff --git a/singularity-eos/closure/kinetic_phasetransition_models.hpp b/singularity-eos/closure/kinetic_phasetransition_models.hpp index 161deebd140..14573db345e 100644 --- a/singularity-eos/closure/kinetic_phasetransition_models.hpp +++ b/singularity-eos/closure/kinetic_phasetransition_models.hpp @@ -33,28 +33,31 @@ namespace singularity { - PORTABLE_INLINE_FUNCTION void LogRatesCGModel(const Real *w, const Real *b, - const int num_phases, const Real *gibbs, const int *gibbsorder, Real *logRjk, int *fromto) { - -// begin the calculation of logRik - int jk = 0; - for (int j = 0; j < num_phases-1; j++){ - // from high Gibb phase to low - for (int k = num_phases-1; k > j; k--){ - // from low Gibb phase to high, k>j always, gibbs_j > gibbs_k always - // always mass transport from j->k - Real hgj=gibbs[gibbsorder[j]]; - Real hgk=gibbs[gibbsorder[k]]; -// calc phase dx from hgi to hgk - Real dgdb=(hgj-hgk)/(b[gibbsorder[j]*num_phases + gibbsorder[k]]); - // First is highest level to levels below, largest gibbs energy diff first. Then follows 2nd highest to levels below. - logRjk[jk] = log(w[gibbsorder[j]*num_phases + gibbsorder[k]]*dgdb)+dgdb*dgdb; - fromto[jk] = gibbsorder[j]*10+gibbsorder[k]; - jk++; - } + const int num_phases, const Real *gibbs, + const int *gibbsorder, Real *logRjk, + int *fromto) { + + // begin the calculation of logRik + int jk = 0; + for (int j = 0; j < num_phases - 1; j++) { + // from high Gibb phase to low + for (int k = num_phases - 1; k > j; k--) { + // from low Gibb phase to high, k>j always, gibbs_j > gibbs_k always + // always mass transport from j->k + Real hgj = gibbs[gibbsorder[j]]; + Real hgk = gibbs[gibbsorder[k]]; + // calc phase dx from hgi to hgk + Real dgdb = (hgj - hgk) / (b[gibbsorder[j] * num_phases + gibbsorder[k]]); + // First is highest level to levels below, largest gibbs energy diff first. Then + // follows 2nd highest to levels below. + logRjk[jk] = + log(w[gibbsorder[j] * num_phases + gibbsorder[k]] * dgdb) + dgdb * dgdb; + fromto[jk] = gibbsorder[j] * 10 + gibbsorder[k]; + jk++; } - return; + } + return; } } // namespace singularity diff --git a/test/test_kpt_models.cpp b/test/test_kpt_models.cpp new file mode 100644 index 00000000000..f8a813af7fc --- /dev/null +++ b/test/test_kpt_models.cpp @@ -0,0 +1,141 @@ +//------------------------------------------------------------------------------ +// © 2024. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#include +#include +#include +#include +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +#include +#include +#include +#include +#include +#include + +using singularity::LogRatesCGModel; + +SCENARIO("First log rate test") { + GIVEN("Parameters for a 5 phase system") { + constexpr Real w[25] = {0., 0.8510E+01, 20., 20., 20., 0.8510E+01, 0., 20., 20., + 20., 20., 20., 0., 20., 20., 20., 20., 20., + 0., 20., 20., 20., 20., 20., 0.}; + + constexpr Real b[25] = {0., 0.3060E-04, 0.1000E-05, 0.1000E-05, 0.1000E-05, + 0.3060E-04, 0., 0.1000E-05, 0.1000E-05, 0.1000E-05, + 0.1000E-05, 0.1000E-05, 0., 0.1000E-05, 0.1000E-05, + 0.1000E-05, 0.1000E-05, 0.1000E-05, 0., 0.1000E-05, + 0.1000E-05, 0.1000E-05, 0.1000E-05, 0.1000E-05, 0.}; + + GIVEN("Gibbs free energy and order") { + constexpr int num = 5; + constexpr int mnum = num * (num - 1) / 2; + constexpr int knum = 5; +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_gibbs("gibbs"); + // order phases from high to low gibbs + Kokkos::View v_order("order"); +#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 gibbs; + std::array order; + auto v_gibbs = gibbs.data(); + auto v_order = order.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize gibbs and order", 0, 1, PORTABLE_LAMBDA(int i) { + v_gibbs[0] = 0.0292971; + v_gibbs[1] = 0.0286937; + v_gibbs[2] = 0.0287409; + v_gibbs[3] = 1.0e99; + v_gibbs[4] = 0.02940076; + v_order[0] = 3; + v_order[1] = 4; + v_order[2] = 0; + v_order[3] = 2; + v_order[4] = 1; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create host-side mirrors of the inputs and copy the inputs. These are + // just used for the comparisons + auto gibbs = Kokkos::create_mirror_view(v_gibbs); + auto order = Kokkos::create_mirror_view(v_order); + Kokkos::deep_copy(gibbs, v_gibbs); + Kokkos::deep_copy(order, v_order); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Gold standard values for a subset of lookups + constexpr std::array logrates_true{ + 1.00000000000000007e+210, 1.00000000000000007e+210, 1.00000000000000007e+210, + 1.00000000000000007e+210, 4.99943400447804946e+05, 4.35424707359967171e+05, + 1.07530324485867186e+04, 3.93959978909504514e+02, 3.09367756860215100e+05, + 2.23469012616621285e+03}; + // create deltagibbs and fromphasetophase +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::View v_dgibbs("dgibbs"); + Kokkos::View v_fromto("fromto"); +#else + std::array dgibbs; + std::array fromto; + auto v_dgibbs = dgibbs.data(); + auto v_fromto = fromto.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + int ik = 0; + for (int i = 0; i < num - 1; i++) { + for (int k = num - 1; k > i; k--) { + v_dgibbs[ik] = v_gibbs[v_order[i]] - v_gibbs[v_order[k]]; + ik++; + } + } +#ifdef PORTABILITY_STRATEGY_KOKKOS + auto dgibbs = Kokkos::create_mirror_view(v_dgibbs); + auto fromto = Kokkos::create_mirror_view(v_fromto); + Kokkos::deep_copy(dgibbs, v_dgibbs); + Kokkos::deep_copy(fromto, v_fromto); +#endif // PORTABILITY_STRATEGY_KOKKOS + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create device views for outputs and mirror those views on the host + Kokkos::View v_logrates("LogRates"); + auto h_logrates = Kokkos::create_mirror_view(v_logrates); +#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_logrates; + // Just alias the existing pointers + auto v_logrates = h_logrates.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + WHEN("A LogRij(w,b,num,gibbs,order,fromto) lookup is performed") { + LogRatesCGModel(w, b, num, v_gibbs, v_order, v_logrates, v_fromto); +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_logrates, v_logrates); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned LogRij(w,b,gibbs,order) should be equal to the true value") { + array_compare(mnum, dgibbs, fromto, h_logrates, logrates_true, "DeltaGibbs", + "FromTo"); + } + } + } + } +} From 69a79490d3fada9691cb7d57a4e2b5b20561515a Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Wed, 7 Aug 2024 18:18:11 -0600 Subject: [PATCH 08/20] Added sorting routine for Gibbs free energy and included it in test. --- singularity-eos/CMakeLists.txt | 1 + .../kinetic_phasetransition_methods.hpp | 59 +++++++++++++++++++ test/test_kpt_models.cpp | 35 +++++++++-- 3 files changed, 89 insertions(+), 6 deletions(-) create mode 100644 singularity-eos/closure/kinetic_phasetransition_methods.hpp diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 34c576534f5..44ab1ec09d2 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -63,6 +63,7 @@ register_headers( eos/variant_utils.hpp eos/eos_carnahan_starling.hpp closure/kinetic_phasetransition_models.hpp + closure/kinetic_phasetransition_methods.hpp ) if (SINGULARITY_BUILD_CLOSURE) diff --git a/singularity-eos/closure/kinetic_phasetransition_methods.hpp b/singularity-eos/closure/kinetic_phasetransition_methods.hpp new file mode 100644 index 00000000000..cceb172281e --- /dev/null +++ b/singularity-eos/closure/kinetic_phasetransition_methods.hpp @@ -0,0 +1,59 @@ +//------------------------------------------------------------------------------ +// © 2024. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#ifndef _SINGULARITY_EOS_CLOSURE_KINETIC_PHSETRANSITION_METHODS_ +#define _SINGULARITY_EOS_CLOSURE_KINETIC_PHSETRANSITION_METHODS_ + +#include +#include +#include +#include +#include + +#include + +#ifdef SINGULARITY_USE_KOKKOSKERNELS +#include +#include +#include +#else +#include +#endif // SINGULARITY_USE_KOKKOSKERNELS + +namespace singularity { + +PORTABLE_INLINE_FUNCTION void SortGibbs(const int num_phases, const Real *gibbs, + int *gibbsorder) { + int itmp = 0; + // initiate gibbsorder + for (int j = 0; j < num_phases; j++) { + gibbsorder[j] = j; + } + // sort from high Gibb phase to low + for (int j = 0; j < num_phases - 1; j++) { + for (int k = num_phases - 1; k > j; k--) { + // from low Gibb phase to high, k>j always, gibbs_j > gibbs_k always + if (gibbs[gibbsorder[k]] > gibbs[gibbsorder[j]]) { + itmp = gibbsorder[j]; + gibbsorder[j] = gibbsorder[k]; + gibbsorder[k] = itmp; + } + } + } + return; +} + +} // namespace singularity + +#endif // _SINGULARITY_EOS_CLOSURE_KINETIC_PHSETRANSITION_METHODS_ diff --git a/test/test_kpt_models.cpp b/test/test_kpt_models.cpp index f8a813af7fc..f9c1fccce8c 100644 --- a/test/test_kpt_models.cpp +++ b/test/test_kpt_models.cpp @@ -25,10 +25,12 @@ #include #include #include +#include #include #include using singularity::LogRatesCGModel; +using singularity::SortGibbs; SCENARIO("First log rate test") { GIVEN("Parameters for a 5 phase system") { @@ -45,7 +47,6 @@ SCENARIO("First log rate test") { GIVEN("Gibbs free energy and order") { constexpr int num = 5; constexpr int mnum = num * (num - 1) / 2; - constexpr int knum = 5; #ifdef PORTABILITY_STRATEGY_KOKKOS // Create Kokkos views on device for the input arrays Kokkos::View v_gibbs("gibbs"); @@ -68,11 +69,6 @@ SCENARIO("First log rate test") { v_gibbs[2] = 0.0287409; v_gibbs[3] = 1.0e99; v_gibbs[4] = 0.02940076; - v_order[0] = 3; - v_order[1] = 4; - v_order[2] = 0; - v_order[3] = 2; - v_order[4] = 1; }); #ifdef PORTABILITY_STRATEGY_KOKKOS // Create host-side mirrors of the inputs and copy the inputs. These are @@ -89,6 +85,27 @@ SCENARIO("First log rate test") { 1.00000000000000007e+210, 4.99943400447804946e+05, 4.35424707359967171e+05, 1.07530324485867186e+04, 3.93959978909504514e+02, 3.09367756860215100e+05, 2.23469012616621285e+03}; + constexpr std::array phaseorder{0, 1, 2, 3, 4}; + constexpr std::array order_true{3, 4, 0, 2, 1}; + + WHEN("Using SortGibbs(num,gibbs,order) to give phases in order of largest to " + "smallest gibbs") { + SortGibbs(num, v_gibbs, v_order); + + std::cout << "Order obtained with SortGibbs, largest to smallest Gibbs: " + << std::endl; + for (int l = 0; l < num; l++) { + std::cout << "Phase " << v_order[l] << " Gibbs: " << v_gibbs[v_order[l]] + << std::endl; + } +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(order, v_order); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned order should be equal to the true order") { + array_compare(num, gibbs, phaseorder, order, order_true, "Gibbs", "phase"); + } + } // create deltagibbs and fromphasetophase #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::View v_dgibbs("dgibbs"); @@ -127,6 +144,12 @@ SCENARIO("First log rate test") { WHEN("A LogRij(w,b,num,gibbs,order,fromto) lookup is performed") { LogRatesCGModel(w, b, num, v_gibbs, v_order, v_logrates, v_fromto); + + std::cout << "LogRates obtained with LogRatesCGModel: " << std::endl; + for (int l = 0; l < mnum; l++) { + std::cout << "From phase i to phase k, ik: " << v_fromto[l] + << " LogRik: " << v_logrates[l] << std::endl; + } #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::fence(); Kokkos::deep_copy(h_logrates, v_logrates); From f5cab057ad800d9c07dff1fee0a32b75b9ac5b9b Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Tue, 13 Aug 2024 14:41:04 -0600 Subject: [PATCH 09/20] Added suggestion for max time step. --- .../kinetic_phasetransition_models.hpp | 22 +++++++ test/test_kpt_models.cpp | 57 ++++++++++++++++++- 2 files changed, 76 insertions(+), 3 deletions(-) diff --git a/singularity-eos/closure/kinetic_phasetransition_models.hpp b/singularity-eos/closure/kinetic_phasetransition_models.hpp index 14573db345e..95197aa7c35 100644 --- a/singularity-eos/closure/kinetic_phasetransition_models.hpp +++ b/singularity-eos/closure/kinetic_phasetransition_models.hpp @@ -60,6 +60,28 @@ PORTABLE_INLINE_FUNCTION void LogRatesCGModel(const Real *w, const Real *b, return; } +PORTABLE_INLINE_FUNCTION Real LogMaxTimeStep(const int num_phases, const Real *mfs, + const int *gibbsorder, const Real *logRjk) { + + Real minmassfraction = 1.e-10; + Real logtimestep = 1.; + int jk = 0; + for (int j = 0; j < num_phases - 1; j++) { + // from high Gibb phase to low + for (int k = num_phases - 1; k > j; k--) { + // First is highest level to levels below, largest gibbs energy diff first. Then + // follows 2nd highest to levels below. + // But largest deltaGibbs does not mean max rate. + if (mfs[gibbsorder[j]] > minmassfraction) { + logtimestep = std::min(logtimestep, -logRjk[jk++]); + } else { + jk++; + } + } + } + return logtimestep; +} + } // namespace singularity #endif // _SINGULARITY_EOS_CLOSURE_KINETIC_PHSETRANSITION_MODELS_ diff --git a/test/test_kpt_models.cpp b/test/test_kpt_models.cpp index f9c1fccce8c..f3405421ef0 100644 --- a/test/test_kpt_models.cpp +++ b/test/test_kpt_models.cpp @@ -29,6 +29,7 @@ #include #include +using singularity::LogMaxTimeStep; using singularity::LogRatesCGModel; using singularity::SortGibbs; @@ -44,7 +45,7 @@ SCENARIO("First log rate test") { 0.1000E-05, 0.1000E-05, 0.1000E-05, 0., 0.1000E-05, 0.1000E-05, 0.1000E-05, 0.1000E-05, 0.1000E-05, 0.}; - GIVEN("Gibbs free energy and order") { + GIVEN("Gibbs free energy") { constexpr int num = 5; constexpr int mnum = num * (num - 1) / 2; #ifdef PORTABILITY_STRATEGY_KOKKOS @@ -63,7 +64,7 @@ SCENARIO("First log rate test") { // Populate the input arrays portableFor( - "Initialize gibbs and order", 0, 1, PORTABLE_LAMBDA(int i) { + "Initialize gibbs", 0, 1, PORTABLE_LAMBDA(int i) { v_gibbs[0] = 0.0292971; v_gibbs[1] = 0.0286937; v_gibbs[2] = 0.0287409; @@ -147,7 +148,7 @@ SCENARIO("First log rate test") { std::cout << "LogRates obtained with LogRatesCGModel: " << std::endl; for (int l = 0; l < mnum; l++) { - std::cout << "From phase i to phase k, ik: " << v_fromto[l] + std::cout << "From phase i to phase k, ik ( x means 0x): " << v_fromto[l] << " LogRik: " << v_logrates[l] << std::endl; } #ifdef PORTABILITY_STRATEGY_KOKKOS @@ -159,6 +160,56 @@ SCENARIO("First log rate test") { "FromTo"); } } + + constexpr Real logmts_true = -2234.69; + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + Kokkos::View v_massfractions("massfractions"); +#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 massfractions; + auto v_massfractions = massfractions.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + // Populate the input arrays + portableFor( + "Initialize mass fractions", 0, 1, PORTABLE_LAMBDA(int i) { + v_massfractions[0] = 0.0; + v_massfractions[1] = 0.5; + v_massfractions[2] = 0.5; + v_massfractions[3] = 0.0; + v_massfractions[4] = 0.0; + }); +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create host-side mirrors of the inputs and copy the inputs. These are + // just used for the comparisons + auto massfractions = Kokkos::create_mirror_view(v_massfractions); + Kokkos::deep_copy(massfractions, v_massfractions); +#endif // PORTABILITY_STRATEGY_KOKKOS + + WHEN("A LogMaxTimeStep(num,order,logrates) lookup is performed") { + Real logmts = LogMaxTimeStep(num, v_massfractions, v_order, v_logrates); + for (int l = 0; l < num; l++) { + std::cout << "Phase: " << l + << " initial mass fractions: " << v_massfractions[l] << std::endl; + } + std::cout << "Log(MaxTimeStep) from Rates obtained with CGModel: " << logmts + << std::endl; + for (int l = 0; l < mnum; l++) { + std::cout << "From phase i to phase k, ik ( x means 0x): " << v_fromto[l] + << " Max mass fraction transfer: " + << std::exp(logmts + v_logrates[l]) << std::endl; + } +#ifdef PORTABILITY_STRATEGY_KOKKOS + Kokkos::fence(); + Kokkos::deep_copy(h_logrates, v_logrates); +#endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned Log(MaxTimeStep) should be equal to the true value") { + CHECK(isClose(logmts, logmts_true)); + } + } } } } From d187609fdcba90f6ba67da7045f50090dc4c2806 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Thu, 22 Aug 2024 13:41:07 -0600 Subject: [PATCH 10/20] All ingredients for a strawman kinetic phasetransition framework implemented. --- singularity-eos/CMakeLists.txt | 1 + .../kinetic_phasetransition_methods.hpp | 67 +++++++++++--- .../kinetic_phasetransition_models.hpp | 4 + .../closure/kinetic_phasetransition_utils.hpp | 59 ++++++++++++ test/test_kpt_models.cpp | 90 +++++++++++++++---- 5 files changed, 194 insertions(+), 27 deletions(-) create mode 100644 singularity-eos/closure/kinetic_phasetransition_utils.hpp diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 44ab1ec09d2..45d29e28418 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -63,6 +63,7 @@ register_headers( eos/variant_utils.hpp eos/eos_carnahan_starling.hpp closure/kinetic_phasetransition_models.hpp + closure/kinetic_phasetransition_utils.hpp closure/kinetic_phasetransition_methods.hpp ) diff --git a/singularity-eos/closure/kinetic_phasetransition_methods.hpp b/singularity-eos/closure/kinetic_phasetransition_methods.hpp index cceb172281e..542018834c4 100644 --- a/singularity-eos/closure/kinetic_phasetransition_methods.hpp +++ b/singularity-eos/closure/kinetic_phasetransition_methods.hpp @@ -33,24 +33,67 @@ namespace singularity { -PORTABLE_INLINE_FUNCTION void SortGibbs(const int num_phases, const Real *gibbs, - int *gibbsorder) { - int itmp = 0; - // initiate gibbsorder - for (int j = 0; j < num_phases; j++) { - gibbsorder[j] = j; - } - // sort from high Gibb phase to low +PORTABLE_INLINE_FUNCTION void SmallStepMFUpdate(const Real logdt, const int num_phases, + const Real *massfractions, + const int *gibbsorder, const Real *logRjk, + Real *dmfs, Real *newmfs) { + + Real minmassfraction = 1.e-10; + + // In logRjk: First is highest level to levels below, largest gibbs energy diff first. + // Then follows 2nd highest to levels below. logRjk[jk], with jk = (j+1)(num_phases-1) - + // (j-1)j/2 - k, contains rate for phase j+1 to k+1 where j=0 is highest gibbs free + // energy state, and j=num_phases-1 is lowest gibbs free energy state. jk=0 is rate from + // highest (j=0) to lowest (k=num_phases-1) gibbs energy states. dmfs[jk] is mass + // transport from phase j+1 to k+1, in the same way. Fill lowest level first and stop + // when origin level is exhausted. This will allow for jumping phases if needed. I also + // think it is most physical. + int jk = 0; + Real dx = 0.; + Real test = 0.; + Real newtest = 0.; + int offset = 0; for (int j = 0; j < num_phases - 1; j++) { + // from high Gibb phase to low + test = massfractions[gibbsorder[j]]; + // all phases except the higest gibbs free energy one (j=0) can get contributions from + // phases with higher gibbs free energy. + for (int k = 0; k < j; k++) { + offset = (k + 1) * (num_phases - 1) - (k - 1) * k / 2; + test = test + dmfs[offset - j]; + } for (int k = num_phases - 1; k > j; k--) { // from low Gibb phase to high, k>j always, gibbs_j > gibbs_k always - if (gibbs[gibbsorder[k]] > gibbs[gibbsorder[j]]) { - itmp = gibbsorder[j]; - gibbsorder[j] = gibbsorder[k]; - gibbsorder[k] = itmp; + // always mass transport from j->k + if (test < minmassfraction) { + dx = 0.; + dmfs[jk++] = dx; + } else { + dx = std::exp(logdt + logRjk[jk] + std::log(massfractions[gibbsorder[j]])); + newtest = test - dx; + if (newtest < minmassfraction) { + dmfs[jk++] = test; + } else { + dmfs[jk++] = dx; + } + test = newtest; } } + if (test < minmassfraction) { + newmfs[gibbsorder[j]] = 0.; + } else { + newmfs[gibbsorder[j]] = test; + } + } + // lowest gibbs free energy phase update + // lowest gibbs free energy phase (j=num_phases-1) can only amass what is coming from + // larger gibbs free energy levels + test = massfractions[gibbsorder[num_phases - 1]]; + for (int k = 0; k < num_phases - 1; k++) { + offset = (k + 1) * (num_phases - 1) - (k - 1) * k / 2; + test = test + dmfs[offset - num_phases + 1]; } + newmfs[gibbsorder[num_phases - 1]] = test; return; } diff --git a/singularity-eos/closure/kinetic_phasetransition_models.hpp b/singularity-eos/closure/kinetic_phasetransition_models.hpp index 95197aa7c35..ef30f36543c 100644 --- a/singularity-eos/closure/kinetic_phasetransition_models.hpp +++ b/singularity-eos/closure/kinetic_phasetransition_models.hpp @@ -57,6 +57,10 @@ PORTABLE_INLINE_FUNCTION void LogRatesCGModel(const Real *w, const Real *b, jk++; } } + // logRjk[jk], with jk = (j+1)(num_phases-1) - (j-1)j/2 - k, contains rate for phase j+1 + // to k+1 where j=0 is highest gibbs free energy state, and j=num_phases-1 is lowest + // gibbs free energy state. jk=0 is rate from highest (j=0) to lowest (k=num_phases-1) + // gibbs energy states. return; } diff --git a/singularity-eos/closure/kinetic_phasetransition_utils.hpp b/singularity-eos/closure/kinetic_phasetransition_utils.hpp new file mode 100644 index 00000000000..c0e1eb49cd2 --- /dev/null +++ b/singularity-eos/closure/kinetic_phasetransition_utils.hpp @@ -0,0 +1,59 @@ +//------------------------------------------------------------------------------ +// © 2024. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#ifndef _SINGULARITY_EOS_CLOSURE_KINETIC_PHSETRANSITION_UTILS_ +#define _SINGULARITY_EOS_CLOSURE_KINETIC_PHSETRANSITION_UTILS_ + +#include +#include +#include +#include +#include + +#include + +#ifdef SINGULARITY_USE_KOKKOSKERNELS +#include +#include +#include +#else +#include +#endif // SINGULARITY_USE_KOKKOSKERNELS + +namespace singularity { + +PORTABLE_INLINE_FUNCTION void SortGibbs(const int num_phases, const Real *gibbs, + int *gibbsorder) { + int itmp = 0; + // initiate gibbsorder + for (int j = 0; j < num_phases; j++) { + gibbsorder[j] = j; + } + // sort from high Gibb phase to low + for (int j = 0; j < num_phases - 1; j++) { + for (int k = num_phases - 1; k > j; k--) { + // from low Gibb phase to high, k>j always, gibbs_j > gibbs_k always + if (gibbs[gibbsorder[k]] > gibbs[gibbsorder[j]]) { + itmp = gibbsorder[j]; + gibbsorder[j] = gibbsorder[k]; + gibbsorder[k] = itmp; + } + } + } + return; +} + +} // namespace singularity + +#endif // _SINGULARITY_EOS_CLOSURE_KINETIC_PHSETRANSITION_UTILS_ diff --git a/test/test_kpt_models.cpp b/test/test_kpt_models.cpp index f3405421ef0..502de270865 100644 --- a/test/test_kpt_models.cpp +++ b/test/test_kpt_models.cpp @@ -27,10 +27,12 @@ #include #include #include +#include #include using singularity::LogMaxTimeStep; using singularity::LogRatesCGModel; +using singularity::SmallStepMFUpdate; using singularity::SortGibbs; SCENARIO("First log rate test") { @@ -50,9 +52,12 @@ SCENARIO("First log rate test") { constexpr int mnum = num * (num - 1) / 2; #ifdef PORTABILITY_STRATEGY_KOKKOS // Create Kokkos views on device for the input arrays + // Create host-side mirrors of the inputs Kokkos::View v_gibbs("gibbs"); + auto gibbs = Kokkos::create_mirror_view(v_gibbs); // order phases from high to low gibbs Kokkos::View v_order("order"); + auto order = Kokkos::create_mirror_view(v_order); #else // Otherwise just create arrays to contain values and create pointers to // be passed to the functions in place of the Kokkos views @@ -72,12 +77,8 @@ SCENARIO("First log rate test") { v_gibbs[4] = 0.02940076; }); #ifdef PORTABILITY_STRATEGY_KOKKOS - // Create host-side mirrors of the inputs and copy the inputs. These are - // just used for the comparisons - auto gibbs = Kokkos::create_mirror_view(v_gibbs); - auto order = Kokkos::create_mirror_view(v_order); + // copy the inputs Kokkos::deep_copy(gibbs, v_gibbs); - Kokkos::deep_copy(order, v_order); #endif // PORTABILITY_STRATEGY_KOKKOS // Gold standard values for a subset of lookups @@ -103,6 +104,7 @@ SCENARIO("First log rate test") { Kokkos::fence(); Kokkos::deep_copy(order, v_order); #endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned order should be equal to the true order") { array_compare(num, gibbs, phaseorder, order, order_true, "Gibbs", "phase"); } @@ -110,7 +112,9 @@ SCENARIO("First log rate test") { // create deltagibbs and fromphasetophase #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::View v_dgibbs("dgibbs"); + auto dgibbs = Kokkos::create_mirror_view(v_dgibbs); Kokkos::View v_fromto("fromto"); + auto fromto = Kokkos::create_mirror_view(v_fromto); #else std::array dgibbs; std::array fromto; @@ -125,10 +129,7 @@ SCENARIO("First log rate test") { } } #ifdef PORTABILITY_STRATEGY_KOKKOS - auto dgibbs = Kokkos::create_mirror_view(v_dgibbs); - auto fromto = Kokkos::create_mirror_view(v_fromto); Kokkos::deep_copy(dgibbs, v_dgibbs); - Kokkos::deep_copy(fromto, v_fromto); #endif // PORTABILITY_STRATEGY_KOKKOS #ifdef PORTABILITY_STRATEGY_KOKKOS @@ -153,8 +154,10 @@ SCENARIO("First log rate test") { } #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::fence(); + Kokkos::deep_copy(fromto, v_fromto); Kokkos::deep_copy(h_logrates, v_logrates); #endif // PORTABILITY_STRATEGY_KOKKOS + THEN("The returned LogRij(w,b,gibbs,order) should be equal to the true value") { array_compare(mnum, dgibbs, fromto, h_logrates, logrates_true, "DeltaGibbs", "FromTo"); @@ -165,7 +168,9 @@ SCENARIO("First log rate test") { #ifdef PORTABILITY_STRATEGY_KOKKOS // Create Kokkos views on device for the input arrays + // Create host-side mirrors of the inputs Kokkos::View v_massfractions("massfractions"); + auto massfractions = Kokkos::create_mirror_view(v_massfractions); #else // Otherwise just create arrays to contain values and create pointers to // be passed to the functions in place of the Kokkos views @@ -183,14 +188,14 @@ SCENARIO("First log rate test") { v_massfractions[4] = 0.0; }); #ifdef PORTABILITY_STRATEGY_KOKKOS - // Create host-side mirrors of the inputs and copy the inputs. These are - // just used for the comparisons - auto massfractions = Kokkos::create_mirror_view(v_massfractions); + // copy the inputs Kokkos::deep_copy(massfractions, v_massfractions); #endif // PORTABILITY_STRATEGY_KOKKOS + Real logmts; + WHEN("A LogMaxTimeStep(num,order,logrates) lookup is performed") { - Real logmts = LogMaxTimeStep(num, v_massfractions, v_order, v_logrates); + logmts = LogMaxTimeStep(num, v_massfractions, v_order, v_logrates); for (int l = 0; l < num; l++) { std::cout << "Phase: " << l << " initial mass fractions: " << v_massfractions[l] << std::endl; @@ -202,12 +207,67 @@ SCENARIO("First log rate test") { << " Max mass fraction transfer: " << std::exp(logmts + v_logrates[l]) << std::endl; } + + THEN("The returned Log(MaxTimeStep) should be equal to the true value") { + CHECK(isClose(logmts, logmts_true)); + } + } + + // reset the input massfractions arrays + portableFor( + "Initialize mass fractions", 0, 1, PORTABLE_LAMBDA(int i) { + v_massfractions[0] = 0.2; + v_massfractions[1] = 0.0; + v_massfractions[2] = 0.1; + v_massfractions[3] = 0.3; + v_massfractions[4] = 0.4; + }); + +#ifdef PORTABILITY_STRATEGY_KOKKOS + // Create Kokkos views on device for the input arrays + // Create host-side mirrors + Kokkos::View v_newmfs("newmassfractions"); + auto newmassfractions = Kokkos::create_mirror_view(v_newmfs); + Kokkos::View v_deltamfs("massfractiontransfers"); + auto h_deltamfs = Kokkos::create_mirror_view(v_deltamfs); +#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 newmassfractions; + auto v_newmfs = newmassfractions.data(); + std::array h_deltamfs; + auto v_deltamfs = h_deltamfs.data(); +#endif // PORTABILITY_STRATEGY_KOKKOS + + constexpr std::array newmassfractions_true{0., 0.800012617417208613, + 0.1999873825827914, 0., 0.}; + constexpr std::array deltamfs_true{ + 0.3, 0., 0., 0., 0.4, 0., 0., 0., 0.2, 0.100012617417208613}; + + WHEN("A massfraction update with SmallStepMFUpdate is performed") { + SmallStepMFUpdate(logmts_true, num, v_massfractions, v_order, v_logrates, + v_deltamfs, v_newmfs); + for (int l = 0; l < num; l++) { + std::cout << "Phase: " << v_order[l] + << " initial mass fractions: " << v_massfractions[v_order[l]] + << std::endl; + std::cout << " " + << " final mass fractions: " << v_newmfs[v_order[l]] << std::endl; + } + for (int l = 0; l < mnum; l++) { + std::cout << "From phase i to phase k, ik ( x means 0x): " << v_fromto[l] + << " Mass fraction transfer: " << v_deltamfs[l] << std::endl; + } #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::fence(); - Kokkos::deep_copy(h_logrates, v_logrates); + Kokkos::deep_copy(newmassfractions, v_newmfs); + Kokkos::deep_copy(h_deltamfs, v_deltamfs); #endif // PORTABILITY_STRATEGY_KOKKOS - THEN("The returned Log(MaxTimeStep) should be equal to the true value") { - CHECK(isClose(logmts, logmts_true)); + THEN("The new mass fractions should be equal to the true values") { + array_compare(num, phaseorder, massfractions, newmassfractions, + newmassfractions_true, "Phase", "old massfractions"); + array_compare(mnum, dgibbs, fromto, h_deltamfs, deltamfs_true, "DeltaGibbs", + "FromTo"); } } } From 09d8077150078997b37590e581bbe717d6c67771 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Tue, 27 Aug 2024 20:50:15 -0600 Subject: [PATCH 11/20] Documentation added. --- doc/sphinx/GibbsOrder.pdf | Bin 0 -> 22734 bytes doc/sphinx/index.rst | 1 + doc/sphinx/src/using-closures.rst | 3 + doc/sphinx/src/using-kpt.rst | 244 ++++++++++++++++++++++++++++++ 4 files changed, 248 insertions(+) create mode 100644 doc/sphinx/GibbsOrder.pdf create mode 100644 doc/sphinx/src/using-kpt.rst diff --git a/doc/sphinx/GibbsOrder.pdf b/doc/sphinx/GibbsOrder.pdf new file mode 100644 index 0000000000000000000000000000000000000000..e77baaad9059574e80af6ba29fd9e562a67bda30 GIT binary patch literal 22734 zcmZsC1C-~?)9zS1p4l;eW81cE+qS)9+qP})7(2FY+q(Py@tyD9b90jJ^ponYs!r$R zsZNJPURacdftCe^r1`Su9flb|53n<|gyH4}&`F!vnmL;TSpGtaFaQ96PSnEM*~Ia0 zZf)RfB5Y!0XKVt)%M0V=>}X?Of6g{jfBzmHg>Wy|*f=z}V|F}dJbgdXUdL7+MEmOW zCBv%gYJH~MMWLPKnaJ!?x|m)lK|%3IisWMA{9XCR{$l(q@@${993S6MBy*@azZ*t% z;0eSS=_h&?0k$$&1~D-s&iK=KC_-c8%$NmyeVX8%Lenr)mCTa;?aOGeF}s;~mPAek zahfy|vK@K21i#t7VNF4Q4P~FT*w1Y9)sq=K^eTrfitI55RsykZB*82q3g(eM_VZ?u zuy~rPs$tP3V=%Wz>)vPfEm`bqyPnm-C6r9_$Ox-MpiyE@EY)dw3fm*?C@RN@oG~z2 zm4shKg`ok^W!hB~?UGPyH|dc~@o~y1!ce{i5r;G?#rhHBD#(&G37}rw1iieJ$zzF3 z8nBA*E1rh9)=#|HEehBNfdtm{Q8``5p2m}nz|?+y*~aNezyr%SYPa zXn;#l)x(Vr0L|if(1eTOf9mg$`b4QrxyBCFNw(-iDaWLQM~Oa=HS`@&#YTVbe-e)6 zE05*#9LQhuQe=l~m+_uOn{zb9h^k9XSw!ml_6xuk+9U{Nb{N6*Bu-p5xqOQouX65Z zPpy*X$DNnsM^>@(F9|?~9tWciX`>So?-igDMJCs!3I;XzNP;+m{xY_%py)^2m1~b6 zq|+lnLShkUZ4!sSG_X2fn@NxJBrUTk8}u{=jg+Kc5;2LnX?16oCM|L*DDN!pPZT5d}*hn1|#o4^B{2thAgD5Dkz}p|G3dI(xWwwIYevcE8eE+GoE2-D8QrXuD5RR&;QTGx0Iv=G85zE4zf2t3d`X-TPO>lNO2+Ie+0m4! zu2PzzDiPj5!BA!>b#mU*yQdHj9nDiuds&See%4i~QdEyBaoZD1$8rS&8{N$)+zRQl z`CeorS{lc!xsid|R&gP;?1IYVsgl5bqHvM$*n2V}KiULy4= z25-)|7}#}-5zB4h!FVQ{%(aG+Dy({37zMSr6=S@G|5TP8JKKoFs8TA-%MMghc|whJ zs#7SntQV~3RB07GRJ5!|z*T0!f<2FB6eA@nJPs;P)|ZO!9yjYM`(L5ZHBY-DB;H~ zjM#+JJCPbM(ekcuX?Tz8y~AlFqiC_o4A-t(w>p7!mcwI;1)#vq=A&igc71^_uJgd+ z!1K9|!!U*?@FIq&wcH)xHa0)%!tAD;c*xwq7qf9R9I!z$&7ZQnNFpf47br~@x1Pol zm-)9I=b1+UZ>Trn^qM~7zZ&Ai^-wn948{4|jx6uL5t&hjP^Zybc1JrO+;*uP(VELt zeljyAYFshfO~ZewYBnMGJq2>S-PdWyx_2Df`W)+)d(FG4X2`;!GHiOvELiRKd|mtP6QVM9z3xo8nxFhl77YV`#EX8I z5CczEOCtkNZ`i^ES-*3KOoC+FxDfa0`m(KToBSgs|Cs;Del6774Wx<3w83P-V9Jus z%(O6HQnhX(i+x_HuVV%a&G_;BNNRp=#QOYhS3~xxuPv|39g%{l&8%6>!X_SH51QYTTj!0(^V9Ao zZrS5*rL~Eh-qkfPpRMb!v*N)ro#52@ z+jg5b;?>}p1})9Cb_sXgm#NJ1V(HdxJMtE?7fI|l_!a0V-MpV$l{N*8iLLQ}Uhsd< z{-Fl{@PU8wzl?y9o}Tp|fq#lj|HB;q#T7){oyC-#|MCO?H}^ll-;K_I{%-+5CnNx1 z0MHp3{7wI|hrbv8ozjWg**gC(X;T_{8hQZhf7!=B4Vcm}{Ds*53;nJ4cmI!m#{ba| zpi}m+Hv!Pe8IBXX9no} zOI+x$Eq~>m0RIr8f3qJkM>`k$|7G95mHZ3+uY*c|RXZBkI@$lD%*f-vumr%_(Z%Gy zC857k!X~a3Mkb14g8zX3abMBI$m$ixZoPgg4bFW1=rmj548#_>O9{r^bWKYjIY z_Qn8UWME;Y|G#M3WmmWns>m9>uFB>-&-25EnW-mZy2L0L2nSSvIF>jN9HK1>F*0p? zFe1NlOrACbsxcz9vziU6`ko%WjH9Jpr$ft2l1R9{j&nO*AX(pg<4nwMvJ!wk~@^1*spN;ZyPA+ zeLD?it3It|n{hGmS9HM&JG05(pDCMCcJGK=H^_n{T4vK(oco2XQ@DcWC#CYsZmi1c zFpJv&<_U;rjYf284d3mbpB8P`kk3}fbXmMV9T1j#KMmTp+w{-wL@InMn4sWl{k=S9 z_tGQ#@ng6_o|zEmKJkz^HGXQ-Ya!+5ahSskg&@WH#e)>XY}9wk<|yzxua7+ecPfEH zmekReXIETnPrGtVm=x$nr?Ritf4xyBM)`)Ic6#6Gom~Jw zzW+%aPyl})FTY|k3Pl(?h7~l$T1nlnQ$mR|$eYa+W@9?~k)e-h;t)r`35XFbW-O=9 z_!`#S(xc=*tRp0TuA`(>MRjTd&Vs9Aq<>7`B#K^JmL*K}o#P zlB(T?1M`>B`~tGID3}dcQBpZP<3_f;?b5^I>1dGUr0S%sjP!gBM`3#zQaov%Mv=wb zg}6%gA~ExtKLio=UC=z#24xgja<3IQLzUnMye}F4{`zZU*UXgO0rZ~97wX2Xx7OZ> zsstel(2IrxRy5Eps!mr*%%62<&VUvl&}DhJ4QgobsSQsL817=k*|iH_|<@c1$Hh=;@rIuzp7qCMuG; z9ONLf50(}=XI&0_SRy-0f)1{;`GsMWMWK*>tDm|=LT15cJjezZ%W=_EWioD3<4$fm zVY!?g!1W{y>lZ<UI?0xXf?F#LW4&lhku$&t+UDy?Pe zHwl*!Q5Bxf#H40^HD5L9x%6hZb)F@#|In7zU3tgrT+sz!3RCuPt{FQU%NSc5>lk|* z3u!-V&uCw3?`VH(4_OkoYN6mtNX}j?v|A#!skRk07M)6%NvJBR+BU~uW=3V#o<6zbB<@6=O>3VE(Tr-r8P6ff+`NEVvd_ldC{(aks>Q^9Pp&M zjwf=}U8WkC7uA(#-1I7WdrW3WORK;Lh z%Tnk!`f@2+^&Jr~C;P8So|`qsvJ!tT^N0^qP7`g@c@Iq;S1$9muWgCgn!_~OZ%%mi zPRq5kq{Q;&fb4T5+#qJifL=`be8`svLSRgf(H7C0(BtlYm4GM$45Atlu~8(VJ`N8k z6QN=AQ1bluN)qE5rI9}=Q8IHRN;ajHsmGk96XW5{=jNBw7^eZ+F%<&#ZSSZg;oUU` z#t}PQ010_&Jd#Bg;dFh=u$YUsKk&Gb0ay9d*XDVe!H$rJKUAGh(ajXQ$ya$k36g2( z3{$8RRQ7_T1204}aAqS_cnI#-cUa9~Jrqnw^VFVaJRM|PM5%;5=wQS#b)ytDgH6&- zlTR{OtmN6YiW=&tH21_AHC%~1(N8jZNPGgoar0~llP)iM1MPvO*1!moFe&U9Hl2mX$8hG^rZgCqcSoOVc`ck8BQFIyKIByEoPJE2IEIb{(B{Ov? z?D)fuu1=Sd2}10Q8a)A0S>#F48LQ0-J1Wn1LV?@dgm-II;NCi&+u|`%5^_ zCyu`|ex1$7QF10=lb+Iqtq$9iQNDWK|Lp7xs|pzg?FOwgnqEeiFuxzfn3k}v!JKm| z?7b&2v}Xlt+wZ{mW2dIoepoD5|N8Q~?OX0PcGRH@i0{esBVBHK(Cy0WhuoXbN9c=$ zc2Kuh{T16Soq^=J{%%ZwHB(`0GVrS^^)34#`mNB5!w1a=vhGaUB3ARTK8qD>n356h zA`8fzrDID)J*QHOA=MajD1%{@?n>(kVHe8S2`cB0CrwR0C4L5_<~T(%5R;!}pTQ1D z1=@qk0+u!lQJbO~A-LVD9D#5yvK5_8FF58WckoPQ$BL#h$%LnuaKp)B7beJU8$5EVP zV4T`y8A#*~vd5Hz-NcD9UxjAoAM<<5=4v=!+!8R?Sfjr^%g6T&|i zd?9rLrY1x?T1+LYoE0WfAp=kQGg8$XV#3IVm<$CEKn>f~X&ei3dvu){ET5S6OwTlp z38fzBC0*YmMcw8n0Zp784-DNjiJZMBJYr+5Q6-fb?fLMiP7EfWP#6Wu%-pF*!*T!*PCv|gsm+2gtD-!Dg0a+kf@_h|xyFZc3>9?viMG1iU;+>3&f(%EaW)t*nQK;%XgNxOv}Er&*n z(=Fq8J+Z~VJ7ewQOp46n?4nE>ysGRfOaja*?Aq%#eABgEB}M~eAWEH@3b*P<`#8Ky zs^?_G#gVTcEkvB@%c19=5bA(++$wnRn%>;Al^B{xcuns5<4|!xz}j38iYxxTFIj++ObQ#VUE9N!8bD^;@t9p&R45Awx+6OVNw zObsK{;4{8!3_GD9>)o5Ga5ZL)*-`-PV%jYR2s;;2QZ#+n`m=0EH3?@+< zl9GbpV(2LQdOI?S5{s*xshF%iZM5<9;ot8K4Rj@phv z-t`2Jc5IVP5keCGAXeI|VvI&-jf?JxE@_RXj4pGEmpq4hhr1%!x?LJ`%_FD=<-YL`(j=^SVhSRv&l?!M%d;YHQGQM($QoSEZmVa6E_KW9))g@Z^Jx zbWkyp4`xVccaoJyv6zE^1m4)*>+z9-8B1OgxQ_OFv3r$z;exDz>SDZjL4C6^fx0F1 z&cy$nUO*D9!Pr)%I$W2kIcW;vU1?5Uiv7WwoSRUV=z8;RQMqjvgvN)G^Tp~YK8sR4 zPXS5HMBwNj0=~s3aF0TsyXL51I*Ut_Hz?Hsv7kjQZK|!vOv)q?>!doXQDyvVTlD;- zQ)_22Nu^93Xh#o5Yv`L9tJsVa4>uaFEO0xSL9w`pCrzTVmZ0+%(DM14@;)9nteXQ`# zjYDb~8_V=k%zF2*Wy!?aYOu`v1|CKe3XRg;GW&{G#ais*9aaTgDA50j+U=pokuHdq zuoy%JDS-yFL*tb&=nEajqe-d*VhKr>4O8c|FlTaf1z}WM=1}a6O2#6a0s5v;YwuuV zvw|1pF9GwNuqHZT`b5DZCh{^8`8x^G?4?|c_F~QVSSNN5bP50aW9x8ty8G+E_e-paB?XT-|yaW@<(Xy7JDn29|GCEPijck!#}JtH045`$ibRpzd`W!q-U){J%M$a|DP zWdesN2WNNRXe7q)5Gn(VcQ3((AWs}HdZ1v#eNEusnKjbJWb3qIm=cRRX4!c3Kf&)r-tgE#R+wjSa1K`aCM#_}Uxn zgC+~CthM;urcAZ0ZZg$U@N?b9qdHrLr~O_+!Ftx4eS|J*e#))jk&Ua#{WilTNI9nFKNhF6 zFOaO<7cPPf;hVZ;joIe_CG-TPa||6NC=vy`R&EMA&i}z#bGSZAQ(K{|NkSxFBqJqc zqsYMkupl*Gw+}Z&OrDUBUqPitYt4WnAY-jWvmVG?)%=d6 z^reJN|0z%EL&Bgzu0z?Wox+;M_7HK`%SxIEjpjXikO4Xx_tG__2xEBcBci7M$K;9n z@FfZsO}z@-nIC_sM9sPQtz@>?lSHv%`axi|F-c-$EssMzs{9j|(aBUA?@U_??|rh~ zsQT8B`w-Lin~see#SaAIy>H5c-&6S;Qf2yAR(9G!HYjzGzJ5GWuEKfix(^y9tutxP z^E<3)lrTC=cG(|&pnWFFW%Ar2km_u)05)=~hdSaH@Pm$ z=tYx>NBfMeimi@7z&QI3!dOMKR@KtJXyT-KGo}oLKJt(xSpkL{R`{&B;srZPxta{O zZ!9CCIIBp>Jd<2Jnq*}$4Kf{YW*bIk8@;0h?Ve;~xK*V!M^t1yG#`}kV(n1d_1`dT zN7J?a550@YguPU{jj1>QW6DPUIidn6G7OrJ^tpEr!rIvkX9vw{2JgV%C~lKfQOmpE z>gn?S+S144K04DW=Xb34#;fJy+WUIF=U8{@*4dzp;(bY>9PgZ>&K2O>m!iko@#8Lc zyX*I%^_(6SbX#NK#yEd_8z;#t2x;9@iGTO;LqEhf2C{rocI5uBc%jas=bp zbrLyzw)7Q8I_|MFrViN_AHgR>^mqsSsX6YD$+U;9HUZ@mj9W1G7^7~^$Rwg!*|DbS zsMLa#&-uc+4}QhATbf9^#*;*iemxzcCo0SzirNors}`Juz@|QFGmbTW{T$toEz}4} zMDL7;EQT^xXPojw>ST}}Z=vYdMu@YP_!>Cj082IW{Xps~BB^lu`>{5)DAg=DvKizg%T%0m% zE?f?E`0V!bpyx;P1{yWF?6{hrB464^t~kE$%(6)qpQAND9unQQZbA5_*jt}Rw?Kb? zI(z;l*}s|Hs`^#j!F|6FI(5Gf8#~Y~*YUPnyx8W=lmw1Al&1+mHH$D}vFH)L!!nlr zxeK#aGWB%}*_a)1@*@=|m085NNnHw$a?7yl_=4>Z<*5fN`@)>Dm69e5wM~fVT%rqQnDQ?}YUbrn|Bo&vcof47HzzaCB^4Ab0F z+2nRb{qD2lm2?Xo28QbhfXF5f@R08td-GD@fPTVApQ%hGIYC$E@Qv)NomFxy z-0Jl`k`b0+_8-`DvH*<$shU05Gbj2vq1}PLEYho;q9fP2zKVSr<$tNx9b%<}Z0X@L zV-u8Qc1M{R7MBwimo>_9l%F=|$j!oAy>DK`ov1ozaGK~!I(Km)RqLzgNssa4xRf+b zHYrpx0=F%6-*ez7`E?F`Kl9m~e*!qRZ_Gt-9Tt0O`m)ym*T=8+Ox9*`KXQw*kaS8+cce{t}@yY zJ4T48MbRnqrusCQs%^EPKu~ z&Xf24i)0VAR;?5FmHTSn<28d2U{(1!#6L}2KqJoSV{w>5dG_rj?xlNHB7Z}Y42eBV zDjwdr;Z<|7gC$n-OqH^1Qm4sSbz-oK<;~PBGciF?%x{cyM6{~*1dT12quJ^Dk)guD=A#zYjmfMwhEjCPoiva zr;hYy`t*pK5TV_0o2J>^Z^)ac%}(Ux&yYH%O|C!DBg05bQFJoMe_(y=D4eZh!u((g zjrcK1GQn${B~CmY)?jgezf`l?X_w2AvF~RX%4~m>X#8eMJr%p%*?!MPbkChsYb&00 z?DN5&6euY!MUT7*qsS??!~?Ll6~wE z$U{c4o~~ifi}~QO^eH&y`<$Ouc$LWESMQfnSYxlx_oH-fx>?=4#xS8J?N2*2DhxYl zjk&I^0)<|N8xJ41T|oz{XQZdBbGk1mI~X6PZ$ySvL{wVpU1LNRxIm_rWNa{xID%eJ z{YKkR+dvw9g)4>|HXV>zJw^&cJKdzeR77@7XZ#pYA-@`)x1wc>O_L_gZ;XSRTc5gWTKX32$|ci1g=!DcYpXhI@V}@>$PS|>@OAx^ z9$i`)w5&ie81#T;xNaneD)zvrR!N#A{EmBCbPs5?z;hZ5(x>E|JA+y!O0~^&pwa0n zZA3E%_LRca?J=4l3b)b@70J(ZtV8)zg$%Nz38los`zBbwFd0r5H4q|bTZHl&6^i;k z^znu%UQLl3fdPEqOMbI@4)=vM$>IOy}P;kBSQ2;)s;2jv6U9%p|>+TG_gk>e6+;6w-C7f zw(HXsjiSE9gtVDFS+s?15{tTEsT!XYR!t7oa{;i4j`f?6d;a@ly% zcv8K5O7R%uTr*EIFP(s~)wp)TPH?usWJiB(8)2GY)>D0M} zn(ha;2PQ$7PM=D=4Q#k##gn0D>D}Lh(ZaOEl?ijEM!m6Lq+7(kjbw%_Yvni!-OJ+h zG3g$d)s&z+UeuOG8jRrw$SX$3sw(_m<6^CO!w(M{> z9P`z>9+T@1IGVwT#%nFK?}Gb5zUzKpYxn1Dx9^V;w-_c}W?;In%vb=lU&fv3C219k zdjT+m=?r0Yc(wrjq1fq`QG0p>=5;0M4_1>7DR9pM!=_H>gj0eqf=AbVq7qk?&KcSs zqq#IxvGUMH%L|BD<*4`brJZ-@Gw<0 zZ$h8)W_87Svur$#3DQ(xG{w9oKA;z8u zw@tuJ?G(F1vsUt6Y_O~yk2}3Oe?AjFUq8fR<*=3Ugy0BGY+fC0I+Jqb@GFF@^`6S` z(zGO&P5yBLnY*ep(CV+anWB$nzXLM$I`(>hpfcQ}JNtoZBeGImhe8(-5Fv++Qp#an zNSp}0(6t~s=EjekbZbqK|0d(f_enK8r##W14(|EqkT(jU0MHTDJo+u??UE#ZdIwk=m?m(Q+7weq@8y)(TR zNUK`+u1DcMfj(YI$8d{qaH>xj_&D1}G3-3--C{nvN?WHF_kDIk=na#yq55`2UYbtc zsz~z#vI$Y1g`Q1NEcveTtop|B$&vwoj+rc|n0!2|2XrxL7eaJulw_$newuYRj*T7n zB_u+7NY8PuigqKAfc)P^6 z?e!94X6D;)|GaJO=Xw~SL$_V;trl%2-KIZWl(0go{j|y&#b#&qW^xrvmil}2@T>iE zw9^j9EAY=MVCP*ct`YA;NOH>5+=f(kwaoiFbr?gmc8Vvu?t^^A&y__sJ#1^xkyb^qi-~W=*GKV2-$OPOFKv zjp0#-{GAV}`j4R{(P63%$4Em0C2BG1!MPB~ zYqni!<3P!Ym*RvGN8$9-3*V(5q_Sy zS~7akdX5rK%gRk8-rV5%sL#7-M_k>FXDDSfV(BLG(iWd`(iR&KSv6f0f%x!okrW%S zGMzYJap4^H96>R`Iey6Pw^b!(|Dm9xk_pxlsPb~JgOP?H46_udRtyZ7IPV_KWA#E} ziFI;X&g1kvpmVV$!_IPxgOx?&JC`!3@g9O?<5QtI)pR-18i{JR<3BH~9za({XG?(O zeBN3&TS;@!Vogx=6vn!_Cl{pj6x9q=yQ)7)D`;7ZT%hd^98#I45`Za(9vLfGBaN^m zbj;EmKq)|4|93$ZD?K*)lqd zJHyn!?_?H~<1k&@)XtD!GWjb@vDo7K0Ftz7b&Q3E$>h8V)UuO~*+p2tSSiygd*YGY zh}YcW6rx8Bl{&NN?;_jh6y3$l$OT4Y{dSPwGobGa846TFQoHx-&j>UsRNxYIdCFj! z3I68m*AR68=(#mS(4RVG%FJs3jnmKKPgNkM1bqm9SD3H1k(?Gk`v!edCz=;d414tc zBeCs+(LOf6@Je`Hje*)P5Atp2V(8*eMFJ2`j-c(L_kc0=Vouy9FwT!kM&yk-dcPRH z(%1gi;;J`(7Y$O`W}LaUE=Kieg5Pp=O3C25mQcbmPs$i}aM!Sjb@QhTA*~tp>1Ka> zIgX7<*2rwdni_dbR!RP+Iq3y4_C4VQwYEXYmB|~o)|obXjnYHSjp~i$%6XSg?R~?t z$EgAwlTUUk!E$ABZQz=Ff=_ppPxmQ%SZdF7fVcs9lpI0V&r9K;4w}U zf(w5b7vHq?eQ_k&R;FlH)1dlXtfCJBH`O*cx@~NnB^^64C9@5#X&0`iEt(r(w(LSm zGga@!It+z<*cW$V9z=at*Ke9|1-nL3YtKe>E!d}|$WBt1G%6oLc@aZX?Cwu+H+Y!AL>}>UB$mg1m$_4-hd__xhJV{t$c;h3b6|27C_pX|H{v!!I#2;<{r0kTJ__w=G!R@p zD?KcGn(QKdc7*{Fg^ZdoW$`L^6aKulbfaaO>7uiRYBZ~h_4S=^$(<<2$}pZvc}ygH z!)X0Hk`pmbSK!LmJLbVRa*X0PCC=dqS`6VB=;)K)HzzsCec`4$Q@JA&*R3UOHVvMeeCFSkPn+`6?Xy{fdTvivI*UCQ$rVe*_oCZcC7dg-o1iXod}|!R-1WOt z4$svW++f=ef*(%aY@7j2<_o+up@ZnbtlXO6`ZZCARYNA?QtWL4f_(LUEE&Jh3Kjuf z%I2k4;Wjt8sb}hGIS^y!vIys{#MJ@)dx9NiDVwc!V_ztRcGTpTEY#lpz6a=sOJ9&a$bE5YvdhmgwC#3VzfzA6l*?fa z3Efz56(hf8DY|cpnwF6FI4B``h2U&0Mi%Ux8Fz>G#KaY--H*IkuyqD%E@;(Hrs|#2 z4)FakC(WOjcZA1`M*DN&8?6(injD`8*xd`neKg!wxN=6|ptoj_&p%<=;kr-h+~1IK z7Ql?kg4M>OM2v zI~^3s8{&56@s^`{L$N*H%eLA(bk0ZOIZ%V4`7wJvn|Xo?zsnaFOxaz*UtD%2-i7e6 z0V(6Z$_pvy8|(+D+|iuG-47*&07a-4dY)et^%K_9N0@@7VUK?R{21HcFx>Iyc5gw* z2srBXvw#vVqt!?mj!d_5$$r}N^m0aKXO?>F$>yS;mT`xE3s25gT+AOTt}jGDAHemg;I zOrlLT-Y-6voh)@sq}{M4s`3!6zlzmfYjB?@y5=AGWt}MuZo{$|vdqkU9orOI*{x zKchWI_R88f{2NFF^Bdx92VtP~6MCylduCtd;8&g)4Mw~o#(k`|9q*s+Kv^fVx>3FW z-reSvaIR6O4fDKT%P#>HSf~ehBDMB#?j&nslt^>B3U&~jc)@x=?euu)JNWHjd;H(c z1;|&|C$`M7A3t4y@a=>V_Ax)Go$)CLsb-c>GtOofDP|rrT5XpUB3jHpIUQEU9zgaO#k&#o)qOPFz*szR|gV zIJgtS+?xq>RPVsN0Fo*NDx@LOTkzZ70iyu|8{%aSz=u&l9ce(wdq{090a-f%zm9-n zx1kg-CJL%)L^iq=7X-WxvZ4|<$)BDCL^4@eb0DnroC!?bOM0+^R z-mWy$Wo{LdbtTUF`3Lq7uwTgZxkG1`zMvM3>6kI|iHl~mCaN;G82z70Sh(s`qxd%r zd}uiQt;DU+62YMb+ntv+Mjkb5qJU>r7#%nz($8sdqgfYvPRBbZesEbXYsv8+t%=Tg6Io zk+B$+)|n74XU~8ZfKi=7&@TI?0@8p7HCFu)1S~pO>|}uS*4{hkWp6Y?{MB&%oDrFYP&t0GJ{GM#U7ai znM!19!X1>XPBD(iqqLBmwj8Eh9?^7dF%=v0y(&|!{+btmYU6Yd!O`L*E%Uqtn6$Fm z#l4}1MyR&KgAuX|5vjr=8!$K4v5|qT4J0smBYCJy*k)r=masZI{=*<+Zg1{-zv0Or zUCI$MyoG0d>z)J}-Z-j6!4_K!4Q_O}C4ad9vyF#siRbmM+e+A3Yjc|3D#l2L?ixFr z1o9Aa)jLtc3Q>QBT7k9R+>ka!?#zQIES*C;AZ;mlBlBs_+K6wNvxU&nr$-Y~YRCjc z2*Y7MU+27C=mjZf)E|joUOUeho8cSQ3SQMP?Be1k79+Tl-QLz%WU^U z@vg9$K~YX?ZJy={x~Z6AbMDdI2a3Es8FT+BjdSjCnwX<1x$%9ZZr{%$C*8Uh3%vr1 zsa&+ix<;(J^@v@3&#q)AZueWeQV*w$Q$DeI;za9Dcmwo zQZ?`}pP}5>ow0TiBJevEp_znsu2kLOv3`Ns=(BKn3!&}kE93MulOv73>`%pQVn$~3 zauaYl31kVI@Hi@|(0ZS@(+Ll2bz%Z}2m9XdvGM5p*Dv4w# zf?RW^-V}!4VN!%KK{#1Wal-@2I}zz zU2s+o=}aX^{_7p;G2j_h(2Jodwe95leL$5gtOvC6f#Jpr+I{iWoN&UFH5d2`oWp8*grDDHU(k!@31mJm2M2EGxKTUAo1kbjE{1<^ z7Ig~Y5&58!fJI;p?i8v>DGn?MP+HAc6ddQu&nczV20OJ(2w5bsnuVR?GJqqKg`(r1 zti=|R>jZG{2hKVVuloWdb6q@FcUrz$?Q^5GA{^4}La1skg zgBa4BxHz~X1PHU@Diqkqx2nYZVu)~J>BZP#uhamuQ@dz^?fP*k#G~NG0>F1Pr{XPf zG2Nn0eK7)TaliBb<>Bb|&f72V0M*1s=Sc#(&?eA1BylnV9e?7uM*Hl5bEi8F`Aai$ zd_|5Tiho=NIRzn(_l9T}<}F_6cZ5k5B0NwUCAow4ULuRj9JTP5O`NB1e%6$f`~go7 z$jAi&*+t!XC;^Hr%!~dgKKvUHDIoPS%P$bnqQD6UYnw{7$>l{mH789%H^u+HPA&|{ z4)jcXR)9qgNbf=*`WD=ZV&BDe-k0RdeG>j<479V6t4J{MZX zMhre@E^g^q0UlyOM5aX;2pWj9GujoNRkuers#WEkieq%hbr4(SjV;NSG%{9?<_82f zC&2Z~jF_B+{61etRj+V}V4#*H8l>rVjjX%2u5?B(&^JAgnl2_zG7!Q;_U{x<(*H5+$ROx`b z-60zY_T4+(@Lf*G{6X~;hI|>1ZO$2%3BO-b-*zgo!HmLI-7jNhuY_safypM~kchzc z>X7&=6MomX4pm0oqWh>t**!TgW>j2I@0Q$~RTnW%D)B!{REaAVEPJc^7%_fk8P>wT)Dqi`|K&m^{ z_7|S%UqE$h3_mO-)seX##kQa27mfoc!)JYUXy)4`BN^QVGO%H&oU~b~ z`-U8d<->?GjSJZjYsrJj7JcN_z;vsr4+zI^DjcNR(<}4 zhZY(^-UlM&KQz5jb{UkUi0SMRu}_uWLHS;%z`u5_h=6jQabdiPTW{M)Zh#lM+W2mI zu5)_VUu`VQIB+Zm|Zs=?Ow91XCC>g6C-xst||Dl>m$&4{Vtx{(Bv7`RmTqQH@ z|Fm)?;83;gKPf^lDZ4O2WSzx~8QIrl-X;4IGHO)m&_&bF((%h*o!H&2Cy<+)vJ32eo1+Q)$lWsUho$=nu(r%rW* z2p3#BQMVRyNeKebjoSZ;!zl`8uQXhW6kNXY!FM!PonbVC{1xn{&v)s?01-+EcoXz%z^mE2C8{G zu>&2Fi&>KM1h3N2Yc(IT{Vz5sgcN=P&Q{^2fMpp^4{72&V`h@VZP`1scN=ft_>4&) zHG9{_r#DO|e=0uqP})e*b;-p_FNrhzkLH{(T_)2P7ie@=nVMs8GnWl&d8JG77n_qC z8ebAgRE?J#Tic0mo>ptH-?D~v%G+RAudOmkaqFHoQ>*7Po92txcA5dpn$huO2VaNl zaMsnSMAT_^jT@#`759(j;ibs9K-;mWRd1>oabzWpa_T>GO?!jRJQ!Dz^;r3Kw+ccn z_@>>{;FWR_-FYH{6Bo?iF=Pb?jmn4l+=~~?M1<0$jz-_RrX7`y*HfQ9E%3K;NAx{H zRg7=Z*UA;;Q&^qbR^)vpuBPHX)}}d!%Zc(S`-*43#n`R1dCA%cEg%~jZ4G+vnBlF> zHLc7w#-!s->esHW-OGwjy5W2zbk43eGWM@JE6I_ET>45Z5&ETg+A31Y=ezlA3VF<$ z=P*k=S-1?wb4tpUs3elOt}^h-Gg>r8|GB`mPePxTQo~ulSmg^%seD_bDtE8e_?4_n z;fa1H>q4ZF5a=db7DC~??%%C-g3J|NBaqDZbG+ug{;w3u_SnkwfhP-yw{W$CJPz%mYx1K~1MP})@db&j+mCIS zUC@pIIT`w+jmidnS^|NP21&r6(iA>8z@!F);on)*S_HHU4x{Ym=z;}-*SjfuVkk6h z6d1l z?Tl#(YkV`iH+b5Vwj0BF*%8+1;lRsdETwMUs)F{!uJg}->kk~>lc4H^#&6!d4*Gi; z0^OMQ&mIOrC^6VzjjOu>GC4RmN05{`&Q00P6SrwtzpOdVp5Wy9oq|sJN$xzLG8DEu z<$%Kx|IxW~gsrm7kFhpvGEgXJr){hC_HUcl-kJqvUYk}ZFwRD6Bnl+xz4l)81wF?Za9jT5D;>pc6%5hTcQG+gmw^Qi;@4kr_^@pyEFoU{F(V~ zq3yv^&pDCo?28notJ$O&K?tZf0u|v-f40Aelvntp^UVgqjRLt9q1X#40=|+SgZD#>S*9 z)B-8I_hoW&ZERwVyfUA|;TQqWS7$?~AkE_{=VBj?a10976$RDC2-hT^zQlqR$J6CR z-Ch*%6xPnrQ&|npCrt=wo{nFv*JPJe_2CU*cgSG4cdasEl^w$~#?RE8f~l_SRCuAU zG?A2aJioQzYUN&EUi>ndow(P`L0{EYP;7*rVdkRhS;lF7{_w@Ywi$Wf2RzW`01_@Q zMxAXCmRAxe)ElXZYNIiFD;se=jf@u8=*|e|%FbJOS^tp@1B)nx-LiZj?9DvBP*As^ zL=5V}db_mwiPkTUb>*5maHseFy&!7spn)Y-nh)2hPV|>BdbajYn%b88-MaJO<*Ho8 zu&zZ)K4$zC*(NIg2@h8fm4XAGbV`{5C-LvV)uKmzE8pIgML0x)i3C`&Z;TYTMz^1O3l{RJSRUsS}3pTeP-*B&rKG@-0Ep*3wY1B z>@iF;q-o}6!WbofTS)MD|A-v_QOUf{20G3@pVnJ%y-GV~9Mu*tixt~{=%OQ?sZH^< zklEiL*KinC7OW#%kwkhI{rrxWcJo0YeyL$>(PaN?v)-X0B8u}>X-I_N(=63oG?!lw zEA) z-Agpb%x^|rqCTm4Q=OAx4|Tp^Vj^Ls+%|67x(z1BuQ^}Z`f13@3+QNBP^Be-uFO7B z%O78IH7vnf**VlvP|A*qb4lb{G^o!P)V>xrwg(rt=Uqpk%K9?v zE67h!?7APs<@O#ryznL5x<}l=)}%a64{}Q-Gzu;j#Z!3Wk64RGbsjz8hB|96i#vki zC5BsHKS|#g{U9Nv%RLPn&gx{PN22d3aw%3~Q9&G!Z1dRFHfPVEWqiX5bPvm&jx*XigJEIC)j4!wV&CV1P2nx413xaozI#9Q|`4zp^%%nW`cefwZl z%(-sG>%9ra#TvV*>tuD+nxb6h4}P>BzTaj?)rMcT>2Y{-M7Rq*14pC#*c3@hdzWxAChe)N>l z+LeIQ!2tnPWl>Vi z&Um{&(_~=OvdA1p(#f#r%9v6hv+^pk5eGXY z-YSQuI3H!zMaX7-!q>6br9~`yztds9?IFpj&J?GzWFoC0{2r^T7F>}ttByF2V12#E zxge{1gyVXi=a}!ethg)lAtMVy#u|w#bk5DkMxQCL7M!YQIm6aih6=Hdw?OY-y=$0# zrsBA`U5F5jC>khOlUHc>sF}#wkeUYF-M5%$a;bz3Us=dT8`x`c8qPSfPo*e3H=6Dq zAFLi{*w}kQ{^F1h{JBeI0<*wle8ZA+2B^&WSeP*X7aRYw(1xp*cx?0}*A>G4E~i_ba=|sehUd=T z^4uNd)8V*6if#{!;ZI}kW8^2ZgY|Ds8ST3|m2}|XhqC_b?H`wuqghr=_FY}Y11G(P z8C^KP;PZLa6Ty6!mrwXt@QG@X16Wj2D!Eir);I*RgYN|2*p~rRcjXM%Wtv^EEH;Kh zBm|x=UU?&*=#aCdowTy5|H|Dz*L7$zNBFK}r=r4a0#oMG)vChs=(V@_*j`_HJIo5Z zL|v_;=-Y0Lji1okE2c66yzR3v7uLn_ZB}vU`axDlENjJ-q zD@jzc-LD_`O$mHyWa+h79D9tdelE|KM{@!AkJn4Po0{p8_Zj+F-wpPMK#5nm1G&G1 z=MV+NxxY=f#bokjyayRI()bJ7=!ipkdE{+s$`B$8y2BEWj-K~i+gmYxKO(ud`HszX z;6+%%71SDyTcU6LZ?M9qt|PEKfn^`!gEQ@?n zgm0>%g2mXdV}>Dbm5}@!hi)PaibG-2Km18qt^W!<{tXzFO|T8QbQ2~Jp7Gy#V} z8w2uD>(F2w+-hjzNKSv;+hWlL0{yFeC^n4MY4J9+N?#)-Sq+$DkYFYYUH2 zhFNd(pWw0eh5PoE{O|p;G2Xubvj5LV{|L-bupds`ay|$`na2MFmu*K!zyA(LQozf9 zdu-r7%HLm6(uN@Gmhu^8+vew-qbQ$kWRW(`o&qn|@#yy?P<8hL5I1C-z8GY?z?cEv z-QEiW$h?YR&^Ww12qFnVNPqp9+Xh%2@n{bx9L7`99q)K-BOO%%5m;b% z0VHDX?&5$)yMneQ?Y0&;cQ-YPDx3%u43!2Wp)wF8pihRh+GwKdZvPLp+irJrPmA)* z+5q>eXkeqk-ElqS>-YX$=WVOgc|6vE9)O`>dSHv<`v-)8!(nic18Bnrg`nVoYbhU) z+jkon3Y7*fx5EaY{+%{381idhFbZ(&PkkXs2n=`>|73$g5MUrLf3g9-LHyEJ1}3xf z8W0&6YUlhQGDtX(c0Y{=*kpio{Mm+tBX{8+I$@Wt_z=qm69uxut?;H=XL4WBB zMERG#fTOqhfPhB>TR8Z2#bINdAC?j$3WzsycL%bR66=j0td^Stkc>YzdOQheJYhYm R!1jZ*484$$y1oYe{{eyZR9gT5 literal 0 HcmV?d00001 diff --git a/doc/sphinx/index.rst b/doc/sphinx/index.rst index ea170e3546f..40b2402f504 100644 --- a/doc/sphinx/index.rst +++ b/doc/sphinx/index.rst @@ -24,6 +24,7 @@ Documentation approved for unlimited release. LA-UR-21-31131. src/modifiers src/customization src/using-closures + src/using-kpt src/python src/contributing src/sphinx-doc diff --git a/doc/sphinx/src/using-closures.rst b/doc/sphinx/src/using-closures.rst index 9675f3a518b..913926b7f88 100644 --- a/doc/sphinx/src/using-closures.rst +++ b/doc/sphinx/src/using-closures.rst @@ -115,6 +115,8 @@ closure state. \rho = \sum_{i=0}^{N - 1} \overline{\rho}_i = \sum_{i=0}^{N-1} \rho_i f_i + +.. _massandvolumefractions: .. note:: Since mass fraction information is encoded in the specification of the @@ -276,6 +278,7 @@ choice of the second independent variable is discussed below and has implications for both the number of additional unknowns and the stability of the method. +.. _density-energy-formalism: The Density-Energy Formulation '''''''''''''''''''''''''''''' diff --git a/doc/sphinx/src/using-kpt.rst b/doc/sphinx/src/using-kpt.rst new file mode 100644 index 00000000000..b822afa25a8 --- /dev/null +++ b/doc/sphinx/src/using-kpt.rst @@ -0,0 +1,244 @@ +.. _using-kpt: + +.. _WillsKPT: https://www.osti.gov/biblio/2202604 + +.. _WillsParameterstudy: https://www.osti.gov/biblio/1900445 + +Kinetic Phase Transition Framework +================================== + +A step in adapting to more realistic calculations is to +get away from that phase transitions are instantanious and the +system is in, at least, local thermal equilibrium (LTE) all the time. +Allowing for superheated (the material is hotter than the equilibrium state) +and supercooled (the material is cooler than the equilibrium state) materials is +a large step towards more realistic modelling. Gas bubbles "exploding" when a +spoon is put into a microwave heated cup of coffee is an example of superheating and +raindrops turning directly to ice when hitting the asphalt is an example +of supercooling. + +Phase transitions are governed by the Gibbs free energy, :math:`G(P,T)`. The equilibrium phase is the phase with lowest +Gibbs free energy. If the material is heated or cooled from this equilibrium state and another +phase emerges to have a lower Gibbs free energy, a phase transition occurs and the system eventually transform +to this new phase. However, all phase transitions need some type of *nucleation event* to get started. +For the liquid raindrop cooling down while falling to the ground through very cold air, the hit against the +asphalt is the nucleation event, and the supercooled liquid in the drop instantaniously transform to the equilibrium ice phase. +For the very hot coffee, the spoon insertion is the "nucleation" event when the very hot liquid coffee +instantaneously transforms into its equilibrium phase, gas. + +In a solid state phase transition, say between an face-centered cubic (fcc) phase and a body-centered +cubic (bcc) phase, the mere rearrangement of atoms can take some time, also resulting in a phase transition +that doesn't go directly from the previously lowest Gibbs free energy phase to the presently +lowest Gibbs free energy phase. + +``singularity-eos`` provides the ingredients to build a kinetic phase transition framework in a host code. + +A comprehensive guide for implementing KPT into a hydro code is available `here `_. + +Mass and volume fractions +------------------------- + +The first ingredient needed in a KPT framework is a quantity describing what phase a material currently +is in. *Mass fractions*, :math:`\mu_i`, gives a measure of how much of the total mass of the material in a cell, :math:`M`, +is in a specific phase + +.. math:: + + \mu_i = \frac{M_i}{M} \qquad \qquad \sum_i \mu_i = 1 + +For completeness we also introduce *volume fractions*, :math:`f_i`, that give how much of the volume of +the material in a cell, :math:`V`, is in a specific phase + +.. math:: + + f_i = \frac{V_i}{V} \qquad \qquad \sum_i f_i = 1 + +With these definitions it is clear that both :math:`0 \leq \mu_i \leq 1` and :math:`0 \leq f_i \leq 1`. +Mass and volume fractions are discussed more closely in the :ref:`mixed cell closure ` section. +Note, however, that mass and volume fractions in the KPT framework give how much of each *phase* of *one* material is present, +while in mixed cells they give how much of one *material* is present. Even though we can use some of the same routines +for mixed material cells and KPT, implementing KPT into a mixed cell environment need to use separate mass and volume +fractions for the mixed cell and within each material in the cell. + +Homogenization +-------------- + +A cell in a simulation can contain several different materials in several different phases. In order to have a well defined state +in this heterogeneous cell we need to have some rules, preferably based in physics, on how to *homogenize* the cell into one average +or mixture "material" with a well defined state. The natural and most common homogenization method for phase transitions is +*pressure-temperature equilibrium* (PTE), which also is a common method for mixed materials cells. + +PTE +''' + +It is natural to use PTE for phase transitions since the Gibbs free energy, :math:`G(P,T)`, is uniquely defined if :math:`P` and :math:`T` +are given. PTE means that every phase/material present in a cell has the same pressure and temperature and +given the specific internal energy, :math:`E`, specific volume, :math:`V`, the set of mass fractions :math:`\{\mu_i\}`, and one Equations of state (EOS) for each phase :math:`i`, we have: + +.. math:: + + V &=& \sum_i \mu_i V_i \\ + E &=& \sum_i \mu_i E_i \\ + P &=& P_i(V_i,E_i) \\ + T &=& T_i(V_i,E_i) + +where :math:`i` denotes the different phases of one material. For :math:`N` phases this gives enough of equations (:math:`2N`) to determine all +:math:`V_i` and :math:`E_i`, and thus :math:`P` and :math:`T` and Gibbs Free energy :math:`G_i(P,T)` for each phase: + +.. math:: + + G_i &=& E_i - P V_i - T S_i \\ + G &=& \sum_i \mu_i G_i + +Note that we need information about the specific entropy, :math:`S`, for each phase. This means that each phase EOS needs to be a :ref:`complete EOS ` with a full description of the states. + +We recommend using the :ref:`density-energy formalism ` (``PTESolverRhoU``) described in the :ref:`mixed cell closure ` section, as PTE solver in the KPT framework. + +Mass fraction update +-------------------- + +In a computational cell at a time :math:`t_0`, the set of current mass fractions, :math:`\{\mu_i\}_{t_0}`, the current specific internal energy, :math:`E_{t_0}`, and specific volume, :math:`V_{t_0}`, +give all the information needed, including the state of the phases, about the current state of the cell via a PTE solve. This state is the ground for advancing the state in the cell to a +state at time :math:`t = t_0 + dt` via mass, momentum, and energy conservation in a hydro code. If a system is considered to be in equilibrium also through a phase transition, this state also gives +the new mass fractions. However, in order to allow kinetics to influence the phase transition a new material specific model is needed, a mass fraction update model. + +Equilibrium phase transitions +''''''''''''''''''''''''''''' + +Tables of equilibrium mass fractions as a function of :math:`P` and :math:`T`, can be constructed at the time of constructing complete EOSs for the different phases of a system. +Phases participation in an equilibrium phase transition all have the same Gibbs free energy even though their specific internal energy, specific volume, and specific entropy, are different +(see equation for Gibbs free energy above). Equilibrium phase transition mass fraction tables will be made available through ``singularity-eos`` in the near future. + +Kinetic phase transition models +''''''''''''''''''''''''''''''' + +We write the new mass fractions at time :math:`t` as + +.. math:: + + \mu_i^{t} = \mu_i^{t_0} + \frac{d\mu_i}{dt} dt + +which with discretized time becomes + +.. math:: + + \mu_i^{t} = \mu_i^{t_0} + \dot{\mu_i} \Delta t + +Using mass conservation (:math:`\sum_i \mu_i = 1`), we see that + +.. math:: + + 0 = \sum_i \dot{\mu_i} = \sum_i \sum_j ( \mu_j R_{ji} - \mu_i R_{ij} ) + +where :math:`R_{ij}` is the mass transportation rate from phase :math:`i` to phase :math:`j`, and we can derive the master equation + +.. math:: + + \dot{\mu_i} = \sum_j ( \mu_j R_{ji} - \mu_i R_{ij} ) + +This equation simply states the fact that all mass transforming from one phase ends up in another phase, it is +just mass conservation. + +A KPT model is a model for the :math:`R_{ij}`. + +Carl Greeff's KPT model +''''''''''''''''''''''' + +Carl Greeff formulated an empirical mass fraction update model as + +.. math:: + + R_{ij} = \nu_{ij} \theta(G_i-G_j) \frac{G_i-G_j}{B_{ij}} \exp\left[ \left( \frac{G_i-G_j}{B_{ij}}\right)^2 \right] + +where :math:`\nu_{ij}` and :math:`B_{ij}` are material dependent fitting constants. +Mattsson-Wills performed a `parameter study `_ of this model, +which can be used as a guide for how to choose these parameters for a specific material and a specific phase transition. + +This model is included in ``singularity-eos`` and its signature is + +.. code-block:: cpp + + LogRatesCGModel(const Real *w, const Real *b, const int num_phases, const Real *gibbs, + const int *gibbsorder, Real *logRij, int *fromto) + +where ``w`` is :math:`\nu_{ij}`, ``b`` is :math:`B_{ij}`, ``num_phases`` is +:math:`N`, the number of phases, ``gibbs`` is :math:`G_i`, ``gibbsorder`` is an array of length :math:`N`, where the ``gibbs`` phase indices are ordered +from the largest Gibbs free energy phase to the lowest Gibbs free energy phase (see figure below), +``logRij`` is :math:`\log(R_{ij})`, and ``fromto`` is a map between the phase indices :math:`i` and :math:`j` and the :math:`ij` phase transition indices in :math:`R_{ij}`. +The ``gibbs`` and ``gibbsorder`` arrays are of length :math:`N` while ``w`` and ``b`` are arrays of length :math:`N^2` representing the :math:`N \times N` matrices, row by row. +Note that it is *NOT* assumed that the phase transition parameters are the same when going from :math:`i \rightarrow j` and :math:`j \rightarrow i`, that is, :math:`\{\nu,B\}_{ij} \neq \{\nu,B\}_{ji}`. +Also note that :math:`R_{ii} = 0`. +``logRij`` is an array of length :math:`N (N-1)/2` giving the logarithm of the non-zero mass transportation rates between phases. The ``fromto`` array +gives to which two phases in the ``gibbs`` array, each rate in ``logRij`` is associated: ``logRij[k]`` is the logarithm of the mass +transportation rate from/to phases ``fromto[k]``, with ``k`` a phase transition index according to the figure below. The integer +in ``fromto``, "ij", is composed from the ``gibbs`` index of the "from" phase, :math:`i`, and the ``gibbs`` index of the "to" phase, +:math:`j`, as :math:`i*10+j`, and with a single digit integer, "x", interpreted as "0x". + +.. image:: ../GibbsOrder.pdf + :width: 500 + :alt: Figure: How the phase transition index used in several arrays relate to the phase index in the gibbsorder array. + +``gibbsorder`` can be obtained with any sorting algorithm. In ``singularity-eos``, ``SortGibbs`` can be used + +.. code-block:: cpp + + SortGibbs(const int num_phases, const Real *gibbs, int *gibbsorder) + +where ``num_phases`` is :math:`N`, ``gibbs`` is an array with length :math:`N`, with Gibbs free energy for each phase. +``gibbsorder`` gives the indices of ``gibbs``, the phase indices, in order from highest Gibbs free energy to lowest Gibbs free energy (see figure above). This means +that ``gibbs[gibbsorder[0]]`` is the highest Gibbs free energy and ``gibbs[gibbsorder[N]]`` is the lowest, that is, the Gibbs free energy of the equilibrium phase. + + +The time step +''''''''''''' + +If a timestep would be truely infinitesimal, :math:`R_{ij} dt \leq 1` would always hold, since however big the +rate :math:`R_{ij}` is, :math:`dt < \frac{1}{R_{ij}}`. This means that the new +mass fractions would always obey :math:`0 \leq \mu_i \leq 1`. However, with a discretized time step, :math:`R_{ij} \Delta t` can become larger than :math:`1`, and it can be that even +if the master equation holds, it results in some phase mass fractions becomming negative and some being above :math:`1`, which is unphysical. + +One way of dealing with this is to use a time step, :math:`\Delta t`, that is smaller than the inverse of the largest rate from an active phase. A routine +suggesting a maximum timestep is available in ``singularity-eos``: + +.. code-block:: cpp + + Real LogMaxTimeStep(const int num_phases, const Real *mfs, + const int *gibbsorder, const Real *logRij) + +where ``num_phases`` is :math:`N`, ``mfs`` is the array containing each phase's (old) mass fraction, :math:`\mu_i`, ``gibbsorder`` contains the ``gibbs`` indices of the phases, +ordered from the phase with the largest to the phase with the smallest Gibbs free energy (see figure above), and ``logRij`` contains :math:`\log(R_{ij})`. The function +gives out :math:`\log(\Delta t_{max})` since it will be used together with :math:`\log(R_{ij})` and both can be very large numbers but with opposite signs so that the difference is +small and can be safely evaluated inside an exponential. + +The update method +''''''''''''''''' + +Because of the numerical sensitivity to the size of the time step, several different methods have been developed +for how to perform the update. The first method made available in ``singularity-eos`` is suitable for simulations where a small timestep can be used: + +.. code-block:: cpp + + SmallStepMFUpdate(const Real logdt, const int num_phases, const Real *massfractions, + const int *gibbsorder, const Real *logRij, Real *dmfs, Real *newmfs) + +where ``logdt`` is :math:`\log(\Delta t)`, ``num_phases`` is :math:`N`, ``massfractions`` is the array containing each phase's (old) mass fraction, +``gibbsorder`` contains the indices of the phases, ordered from the phase with the largest to the phase with +the smallest Gibbs free energy (see figure above), ``logRij`` is :math:`\log(R_{ij})`, ``dmfs`` is the mass transformed from one phase to another, +:math:`\mu_i R_{ij} \Delta t`, for each phase transition in the order described in the figure above, and ``newmfs`` is containing the new, updated, +mass fractions. + +The advantage of using Gibbs ordered phases in the internal calculations is shown in the figure above. +All phase transitions will always go from a higher Gibbs free energy phase to a smaller Gibbs free energy phase, +and by using the indexing scheme in the figure all mass transformed will always go from a phase with a lower index +to a phase with a larger index. In addition, the rates are usually larger when the Gibbs free energy difference is larger +(even though the material and phase transition fitting constants could reverse the order of the rates), and dealing with the phase transitions +in the order shown in the figure facilitates the calculations. Using the Gibbs order indices, the connection between these indices, :math:`j` and :math:`k`, and +the phase transition indices :math:`jk` is + +.. math:: + + jk = (j+1)(N-1) - (j-1)j/2 - k + +as can be verified by hand in the figure above. + + From f567449bc67a8c72cdd49c901f927da80b89d7cd Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Wed, 28 Aug 2024 13:04:07 -0600 Subject: [PATCH 12/20] Addressed requested changes. --- doc/sphinx/src/using-kpt.rst | 6 ++- singularity-eos/CMakeLists.txt | 10 +++-- .../kinetic_phasetransition_methods.hpp | 45 +++++++++---------- .../kinetic_phasetransition_models.hpp | 9 +++- .../closure/kinetic_phasetransition_utils.hpp | 2 + test/CMakeLists.txt | 2 +- 6 files changed, 42 insertions(+), 32 deletions(-) diff --git a/doc/sphinx/src/using-kpt.rst b/doc/sphinx/src/using-kpt.rst index b822afa25a8..be0de529f63 100644 --- a/doc/sphinx/src/using-kpt.rst +++ b/doc/sphinx/src/using-kpt.rst @@ -186,7 +186,7 @@ in ``fromto``, "ij", is composed from the ``gibbs`` index of the "from" phase, : where ``num_phases`` is :math:`N`, ``gibbs`` is an array with length :math:`N`, with Gibbs free energy for each phase. ``gibbsorder`` gives the indices of ``gibbs``, the phase indices, in order from highest Gibbs free energy to lowest Gibbs free energy (see figure above). This means -that ``gibbs[gibbsorder[0]]`` is the highest Gibbs free energy and ``gibbs[gibbsorder[N]]`` is the lowest, that is, the Gibbs free energy of the equilibrium phase. +that ``gibbs[gibbsorder[0]]`` is the highest Gibbs free energy and ``gibbs[gibbsorder[N-1]]`` is the lowest, that is, the Gibbs free energy of the equilibrium phase. The time step @@ -241,4 +241,8 @@ the phase transition indices :math:`jk` is as can be verified by hand in the figure above. +Note that this method depleats phases in order of mass transfer to the lowest Gibbs free energy first, then the next lowest, and so on (see figure above), +but stops once the originating phase is depleated. If this is reflecting the physical reality, is up to the user to decide. The size of the time step +problem is taken care of by never transfering more that the existing mass in a phase to any other phase. + diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt index 4e7578c7e31..5f6d68f9b7a 100644 --- a/singularity-eos/CMakeLists.txt +++ b/singularity-eos/CMakeLists.txt @@ -63,13 +63,15 @@ register_headers( eos/singularity_eos_init_utils.hpp eos/variant_utils.hpp eos/eos_carnahan_starling.hpp - closure/kinetic_phasetransition_models.hpp - closure/kinetic_phasetransition_utils.hpp - closure/kinetic_phasetransition_methods.hpp ) if (SINGULARITY_BUILD_CLOSURE) - register_headers(closure/mixed_cell_models.hpp) + register_headers( + closure/mixed_cell_models.hpp + closure/kinetic_phasetransition_models.hpp + closure/kinetic_phasetransition_utils.hpp + closure/kinetic_phasetransition_methods.hpp + ) if (SINGULARITY_USE_FORTRAN OR SINGULARITY_BUILD_TESTS) # while these are C++ files, they # are only needed for the fortran backend or unit testing diff --git a/singularity-eos/closure/kinetic_phasetransition_methods.hpp b/singularity-eos/closure/kinetic_phasetransition_methods.hpp index 542018834c4..bf29793a23c 100644 --- a/singularity-eos/closure/kinetic_phasetransition_methods.hpp +++ b/singularity-eos/closure/kinetic_phasetransition_methods.hpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -38,8 +39,6 @@ PORTABLE_INLINE_FUNCTION void SmallStepMFUpdate(const Real logdt, const int num_ const int *gibbsorder, const Real *logRjk, Real *dmfs, Real *newmfs) { - Real minmassfraction = 1.e-10; - // In logRjk: First is highest level to levels below, largest gibbs energy diff first. // Then follows 2nd highest to levels below. logRjk[jk], with jk = (j+1)(num_phases-1) - // (j-1)j/2 - k, contains rate for phase j+1 to k+1 where j=0 is highest gibbs free @@ -49,51 +48,49 @@ PORTABLE_INLINE_FUNCTION void SmallStepMFUpdate(const Real logdt, const int num_ // when origin level is exhausted. This will allow for jumping phases if needed. I also // think it is most physical. int jk = 0; - Real dx = 0.; - Real test = 0.; - Real newtest = 0.; - int offset = 0; for (int j = 0; j < num_phases - 1; j++) { // from high Gibb phase to low - test = massfractions[gibbsorder[j]]; + Real new_massfrac = massfractions[gibbsorder[j]]; // all phases except the higest gibbs free energy one (j=0) can get contributions from - // phases with higher gibbs free energy. + // phases with higher gibbs free energy (lower gibbsorder index) for (int k = 0; k < j; k++) { - offset = (k + 1) * (num_phases - 1) - (k - 1) * k / 2; - test = test + dmfs[offset - j]; + int offset = (k + 1) * (num_phases - 1) - (k - 1) * k / 2; + new_massfrac = new_massfrac + dmfs[offset - j]; } for (int k = num_phases - 1; k > j; k--) { // from low Gibb phase to high, k>j always, gibbs_j > gibbs_k always // always mass transport from j->k - if (test < minmassfraction) { - dx = 0.; + if (new_massfrac < KPT_MIN_MASS_FRACTION) { + // Phase depleated, nothing more to transfer + Real dx = 0.; dmfs[jk++] = dx; } else { - dx = std::exp(logdt + logRjk[jk] + std::log(massfractions[gibbsorder[j]])); - newtest = test - dx; - if (newtest < minmassfraction) { - dmfs[jk++] = test; + Real dx = std::exp(logdt + logRjk[jk] + std::log(massfractions[gibbsorder[j]])); + Real new_trial_massfrac = new_massfrac - dx; + if (new_trial_massfrac < KPT_MIN_MASS_FRACTION) { + // Phase will be depleated, only take away what is left + dmfs[jk++] = new_massfrac; } else { dmfs[jk++] = dx; } - test = newtest; + new_massfrac = new_trial_massfrac; } } - if (test < minmassfraction) { + if (new_massfrac < KPT_MIN_MASS_FRACTION) { newmfs[gibbsorder[j]] = 0.; } else { - newmfs[gibbsorder[j]] = test; + newmfs[gibbsorder[j]] = new_massfrac; } } // lowest gibbs free energy phase update // lowest gibbs free energy phase (j=num_phases-1) can only amass what is coming from - // larger gibbs free energy levels - test = massfractions[gibbsorder[num_phases - 1]]; + // larger gibbs free energy levels (lower gibbsorder indices) + Real new_massfrac = massfractions[gibbsorder[num_phases - 1]]; for (int k = 0; k < num_phases - 1; k++) { - offset = (k + 1) * (num_phases - 1) - (k - 1) * k / 2; - test = test + dmfs[offset - num_phases + 1]; + int offset = (k + 1) * (num_phases - 1) - (k - 1) * k / 2; + new_massfrac = new_massfrac + dmfs[offset - num_phases + 1]; } - newmfs[gibbsorder[num_phases - 1]] = test; + newmfs[gibbsorder[num_phases - 1]] = new_massfrac; return; } diff --git a/singularity-eos/closure/kinetic_phasetransition_models.hpp b/singularity-eos/closure/kinetic_phasetransition_models.hpp index ef30f36543c..f3c9f9b4d04 100644 --- a/singularity-eos/closure/kinetic_phasetransition_models.hpp +++ b/singularity-eos/closure/kinetic_phasetransition_models.hpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include @@ -53,6 +54,11 @@ PORTABLE_INLINE_FUNCTION void LogRatesCGModel(const Real *w, const Real *b, // follows 2nd highest to levels below. logRjk[jk] = log(w[gibbsorder[j] * num_phases + gibbsorder[k]] * dgdb) + dgdb * dgdb; + // The fromto is just a way of labeling what gibbs (user/host) indexed phases + // (assuming max 10 phases) the phase transition index (jk), used in logRjk, are + // refering to. Internally the gibbsorder indexing (0 for highest Gibbs phase and + // N-1 for lowest Gibbs phase) is always used. + // x should be interpreted as 0x fromto[jk] = gibbsorder[j] * 10 + gibbsorder[k]; jk++; } @@ -67,7 +73,6 @@ PORTABLE_INLINE_FUNCTION void LogRatesCGModel(const Real *w, const Real *b, PORTABLE_INLINE_FUNCTION Real LogMaxTimeStep(const int num_phases, const Real *mfs, const int *gibbsorder, const Real *logRjk) { - Real minmassfraction = 1.e-10; Real logtimestep = 1.; int jk = 0; for (int j = 0; j < num_phases - 1; j++) { @@ -76,7 +81,7 @@ PORTABLE_INLINE_FUNCTION Real LogMaxTimeStep(const int num_phases, const Real *m // First is highest level to levels below, largest gibbs energy diff first. Then // follows 2nd highest to levels below. // But largest deltaGibbs does not mean max rate. - if (mfs[gibbsorder[j]] > minmassfraction) { + if (mfs[gibbsorder[j]] > KPT_MIN_MASS_FRACTION) { logtimestep = std::min(logtimestep, -logRjk[jk++]); } else { jk++; diff --git a/singularity-eos/closure/kinetic_phasetransition_utils.hpp b/singularity-eos/closure/kinetic_phasetransition_utils.hpp index c0e1eb49cd2..6ab10a1a349 100644 --- a/singularity-eos/closure/kinetic_phasetransition_utils.hpp +++ b/singularity-eos/closure/kinetic_phasetransition_utils.hpp @@ -33,6 +33,8 @@ namespace singularity { +constexpr Real KPT_MIN_MASS_FRACTION = 1.0e-10; + PORTABLE_INLINE_FUNCTION void SortGibbs(const int num_phases, const Real *gibbs, int *gibbsorder) { int itmp = 0; diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b67f8c19729..1eb6563d126 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -111,7 +111,7 @@ if(SINGULARITY_BUILD_CLOSURE) singularity-eos::singularity-eos) target_include_directories(test_pte PRIVATE ${CMAKE_CURRENT_SOURCE_DIR}) add_test(pte test_pte) -##### + add_executable(test_pte_2phase test_pte_2phase.cpp) target_link_libraries(test_pte_2phase PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) From 5ca3c2267aea76e5a3ba0392478b80f6c963a3ec Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Thu, 29 Aug 2024 10:56:00 -0600 Subject: [PATCH 13/20] Minor changes, mostly comments and documentation. --- CHANGELOG.md | 3 ++- doc/sphinx/src/using-kpt.rst | 17 +++++++++++++---- .../closure/kinetic_phasetransition_methods.hpp | 1 + .../closure/kinetic_phasetransition_models.hpp | 1 + test/test_kpt_models.cpp | 4 ++++ 5 files changed, 21 insertions(+), 5 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a920d422aab..cdd2e1bdcc7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR361]](https://github.com/lanl/singularity-eos/pull/361) Added tests for PTEsolver and added a strawman kinetic phase transition framework - [[PR330]](https://github.com/lanl/singularity-eos/pull/330) Piecewise grids for Spiner EOS. ### Fixed (Repair bugs, etc) @@ -183,7 +184,7 @@ Date: 07/07/2022 This is the start of changelog -© 2021-2023. Triad National Security, LLC. All rights reserved. This +© 2021-2024. Triad National Security, LLC. All rights reserved. This program was produced under U.S. Government contract 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC for the U.S. Department of Energy/National diff --git a/doc/sphinx/src/using-kpt.rst b/doc/sphinx/src/using-kpt.rst index be0de529f63..c4265603e2a 100644 --- a/doc/sphinx/src/using-kpt.rst +++ b/doc/sphinx/src/using-kpt.rst @@ -150,7 +150,8 @@ Carl Greeff formulated an empirical mass fraction update model as R_{ij} = \nu_{ij} \theta(G_i-G_j) \frac{G_i-G_j}{B_{ij}} \exp\left[ \left( \frac{G_i-G_j}{B_{ij}}\right)^2 \right] -where :math:`\nu_{ij}` and :math:`B_{ij}` are material dependent fitting constants. +where :math:`\nu_{ij}` and :math:`B_{ij}` are material dependent fitting constants, and :math:`\theta` is the Heaviside step function that is +:math:`0` for a negative argument and :math:`1` for a positive argument. Mattsson-Wills performed a `parameter study `_ of this model, which can be used as a guide for how to choose these parameters for a specific material and a specific phase transition. @@ -162,14 +163,22 @@ This model is included in ``singularity-eos`` and its signature is const int *gibbsorder, Real *logRij, int *fromto) where ``w`` is :math:`\nu_{ij}`, ``b`` is :math:`B_{ij}`, ``num_phases`` is -:math:`N`, the number of phases, ``gibbs`` is :math:`G_i`, ``gibbsorder`` is an array of length :math:`N`, where the ``gibbs`` phase indices are ordered +:math:`N`, the number of phases, and ``gibbs`` is :math:`G_i`. +``gibbsorder`` is an array of length :math:`N`, where the ``gibbs`` phase indices are ordered from the largest Gibbs free energy phase to the lowest Gibbs free energy phase (see figure below), ``logRij`` is :math:`\log(R_{ij})`, and ``fromto`` is a map between the phase indices :math:`i` and :math:`j` and the :math:`ij` phase transition indices in :math:`R_{ij}`. + The ``gibbs`` and ``gibbsorder`` arrays are of length :math:`N` while ``w`` and ``b`` are arrays of length :math:`N^2` representing the :math:`N \times N` matrices, row by row. Note that it is *NOT* assumed that the phase transition parameters are the same when going from :math:`i \rightarrow j` and :math:`j \rightarrow i`, that is, :math:`\{\nu,B\}_{ij} \neq \{\nu,B\}_{ji}`. Also note that :math:`R_{ii} = 0`. -``logRij`` is an array of length :math:`N (N-1)/2` giving the logarithm of the non-zero mass transportation rates between phases. The ``fromto`` array -gives to which two phases in the ``gibbs`` array, each rate in ``logRij`` is associated: ``logRij[k]`` is the logarithm of the mass + +``logRij`` is an array of length :math:`N (N-1)/2` giving the logarithm of the non-zero mass transportation rates between phases. +Using ``gibbsorder`` indices, :math:`j` and :math:`k`, we see that all :math:`R_{jk}` with :math:`j \geq k` are :math:`0` because of +:math:`\theta(G_j - G_k) = 0` (since :math:`G_k > G_j`). Writing :math:`R_{jk}` on matrix form would give that only the upper triangular part is non-zero, +giving :math:`N (N-1)/2` non-zero elements. + +The ``fromto`` array +gives which two phases in the ``gibbs`` array, each rate in ``logRij`` is associated with: ``logRij[k]`` is the logarithm of the mass transportation rate from/to phases ``fromto[k]``, with ``k`` a phase transition index according to the figure below. The integer in ``fromto``, "ij", is composed from the ``gibbs`` index of the "from" phase, :math:`i`, and the ``gibbs`` index of the "to" phase, :math:`j`, as :math:`i*10+j`, and with a single digit integer, "x", interpreted as "0x". diff --git a/singularity-eos/closure/kinetic_phasetransition_methods.hpp b/singularity-eos/closure/kinetic_phasetransition_methods.hpp index bf29793a23c..a8ccb68b944 100644 --- a/singularity-eos/closure/kinetic_phasetransition_methods.hpp +++ b/singularity-eos/closure/kinetic_phasetransition_methods.hpp @@ -65,6 +65,7 @@ PORTABLE_INLINE_FUNCTION void SmallStepMFUpdate(const Real logdt, const int num_ Real dx = 0.; dmfs[jk++] = dx; } else { + // we do this multiplication in log space to avoid overflow Real dx = std::exp(logdt + logRjk[jk] + std::log(massfractions[gibbsorder[j]])); Real new_trial_massfrac = new_massfrac - dx; if (new_trial_massfrac < KPT_MIN_MASS_FRACTION) { diff --git a/singularity-eos/closure/kinetic_phasetransition_models.hpp b/singularity-eos/closure/kinetic_phasetransition_models.hpp index f3c9f9b4d04..786b7a9cac7 100644 --- a/singularity-eos/closure/kinetic_phasetransition_models.hpp +++ b/singularity-eos/closure/kinetic_phasetransition_models.hpp @@ -59,6 +59,7 @@ PORTABLE_INLINE_FUNCTION void LogRatesCGModel(const Real *w, const Real *b, // refering to. Internally the gibbsorder indexing (0 for highest Gibbs phase and // N-1 for lowest Gibbs phase) is always used. // x should be interpreted as 0x + // TODO: Change to struct or some more generalizable form fromto[jk] = gibbsorder[j] * 10 + gibbsorder[k]; jk++; } diff --git a/test/test_kpt_models.cpp b/test/test_kpt_models.cpp index 502de270865..25356827911 100644 --- a/test/test_kpt_models.cpp +++ b/test/test_kpt_models.cpp @@ -12,6 +12,9 @@ // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ +// TODO: Need to make this working on GPUs +#ifdef PORTABILITY_STRATEGY_NONE + #include #include #include @@ -273,3 +276,4 @@ SCENARIO("First log rate test") { } } } +#endif // PORTABILITY_STRATEGY_NONE From f3e89b61a579d01d1d66d2d58a765c33f3c39a64 Mon Sep 17 00:00:00 2001 From: Ann Elisabet Wills - 298385 Date: Thu, 29 Aug 2024 14:56:28 -0600 Subject: [PATCH 14/20] Fixing so that kpt test case runs. --- test/test_kpt_models.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/test/test_kpt_models.cpp b/test/test_kpt_models.cpp index 25356827911..911443bc2be 100644 --- a/test/test_kpt_models.cpp +++ b/test/test_kpt_models.cpp @@ -13,7 +13,8 @@ //------------------------------------------------------------------------------ // TODO: Need to make this working on GPUs -#ifdef PORTABILITY_STRATEGY_NONE +//#ifdef PORTABILITY_STRATEGY_NONE +#ifndef PORTABILITY_STRATEGY_KOKKOS #include #include @@ -38,7 +39,7 @@ using singularity::LogRatesCGModel; using singularity::SmallStepMFUpdate; using singularity::SortGibbs; -SCENARIO("First log rate test") { +SCENARIO("Testing the KPT framework") { GIVEN("Parameters for a 5 phase system") { constexpr Real w[25] = {0., 0.8510E+01, 20., 20., 20., 0.8510E+01, 0., 20., 20., 20., 20., 20., 0., 20., 20., 20., 20., 20., @@ -276,4 +277,5 @@ SCENARIO("First log rate test") { } } } -#endif // PORTABILITY_STRATEGY_NONE +//#endif // PORTABILITY_STRATEGY_NONE +#endif // PORTABILITY_STRATEGY_KOKKOS From 25e41c7836cb13357a7e90b920d93c3eccf3c855 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 4 Sep 2024 16:02:16 -0600 Subject: [PATCH 15/20] update kokkos version to support newer cuda versions --- utils/kokkos | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utils/kokkos b/utils/kokkos index 5ad609661e5..08ceff92bcf 160000 --- a/utils/kokkos +++ b/utils/kokkos @@ -1 +1 @@ -Subproject commit 5ad609661e570ba6aa7716a26a91cb67d559f8a2 +Subproject commit 08ceff92bcf3a828844480bc1e6137eb74028517 From 209953e45eeaa9b0fbb3ae63c1e3e80de392d7ca Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 4 Sep 2024 16:08:11 -0600 Subject: [PATCH 16/20] update kokkos kernels --- cmake/singularity-eos/kokkos.cmake | 3 +++ utils/kokkos-kernels | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/cmake/singularity-eos/kokkos.cmake b/cmake/singularity-eos/kokkos.cmake index 6d56a8be321..832d9d190e2 100644 --- a/cmake/singularity-eos/kokkos.cmake +++ b/cmake/singularity-eos/kokkos.cmake @@ -49,6 +49,9 @@ macro(singularity_import_kokkoskernels) set(KokkosKernels_ENABLE_TPL_CUSPARSE OFF CACHE BOOL "" FORCE) + set(KokkosKernels_ENABLE_TPL_CUSOLVER + OFF + CACHE BOOL "" FORCE) set(KokkosKernels_ENABLE_TPL_MAGMA OFF CACHE BOOL "" FORCE) diff --git a/utils/kokkos-kernels b/utils/kokkos-kernels index ddf1b3df43f..0608a337389 160000 --- a/utils/kokkos-kernels +++ b/utils/kokkos-kernels @@ -1 +1 @@ -Subproject commit ddf1b3df43fed1674e166516e4f7169ccce3d224 +Subproject commit 0608a337389ceb7e2695570db830230c2703bfe0 From 1a9f49fddc1a72b6194524d0baa9cd28c950b22f Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 6 Sep 2024 14:02:17 -0600 Subject: [PATCH 17/20] aematts tests now build/run on device. Small issues with cout vs printf and issue that eospac does not run on device easily. Comment out eospac for device builds. Fix printf etc. Speed up builds by limiting what the variant is allowed to be for the tests. --- test/pte_longtest_2phaseVinetSn.hpp | 5 +++-- test/pte_test_2phaseVinetSn.hpp | 9 +++++---- test/pte_test_5phaseSesameSn.hpp | 23 ++++++++++++++--------- test/pte_test_first.hpp | 15 +++++++++++---- test/test_pte.cpp | 6 ++++-- test/test_pte_2phase.cpp | 18 ++++++++++++------ test/test_pte_3phase.cpp | 26 +++++++++++++++++++------- 7 files changed, 68 insertions(+), 34 deletions(-) diff --git a/test/pte_longtest_2phaseVinetSn.hpp b/test/pte_longtest_2phaseVinetSn.hpp index 6aa78c7371f..eac2ce665a8 100644 --- a/test/pte_longtest_2phaseVinetSn.hpp +++ b/test/pte_longtest_2phaseVinetSn.hpp @@ -26,6 +26,7 @@ // in .ult but they will be in the near future. namespace pte_longtest_2phaseVinetSn { + constexpr int NMAT = 2; constexpr int NTRIAL = 20; constexpr int NPTS = NTRIAL * NMAT; @@ -82,10 +83,10 @@ inline void set_eos(T *eos) { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; - singularity::EOS Snbeta = + EOS Snbeta = singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, 0.658e09, 0.4419e07, d2to40); - singularity::EOS Sngamma = + EOS Sngamma = singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, 0.2161e07, 1.025e09, 0.5051e07, d2to40); eos[0] = Snbeta.GetOnDevice(); diff --git a/test/pte_test_2phaseVinetSn.hpp b/test/pte_test_2phaseVinetSn.hpp index 595620979ea..efde8919f37 100644 --- a/test/pte_test_2phaseVinetSn.hpp +++ b/test/pte_test_2phaseVinetSn.hpp @@ -19,12 +19,13 @@ #include #include -#include + +#include +#include // This data is taken from the .ult file for a test case in Flag, MPMat_KPT_PTsolv.ult.std // Note that at this point in time (April 2024) the global and phase times are not alinged // in .ult but they will be in the near future. - namespace pte_test_2phaseVinetSn { constexpr int NMAT = 2; @@ -58,10 +59,10 @@ inline void set_eos(T *eos) { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; - singularity::EOS Snbeta = + EOS Snbeta = singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, 0.658e09, 0.4419e07, d2to40); - singularity::EOS Sngamma = + EOS Sngamma = singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, 0.2161e07, 1.025e09, 0.5051e07, d2to40); eos[0] = Snbeta.GetOnDevice(); diff --git a/test/pte_test_5phaseSesameSn.hpp b/test/pte_test_5phaseSesameSn.hpp index 75664066b03..723296a7ba0 100644 --- a/test/pte_test_5phaseSesameSn.hpp +++ b/test/pte_test_5phaseSesameSn.hpp @@ -21,7 +21,8 @@ #include #include -#include +#include +#include // This data is taken from the .ult file for a test case in Flag, // MPMat_KPT_3phases.ult.std in MPMat_KPT_3phases_NW. Note that 3 phases appear by mistake @@ -61,6 +62,10 @@ constexpr Real out_sie1[NTRIAL] = {1.93716882e09, 1.93864981e09, 1.93994981e09, constexpr Real out_sie2[NTRIAL] = {2.01473728e09, 2.01621655e09, 2.01751522e09, 2.01840528e09, 2.01875845e09}; +using singularity::EOSPAC; +using singularity::Variant; +using EOS = Variant; + template inline void set_eos(T *eos) { @@ -70,17 +75,17 @@ inline void set_eos(T *eos) { int SnbetaID = 2102; int SngammaID = 2103; - int SndeltaID = 2104; + // int SndeltaID = 2104; int SnhcpID = 2105; - int SnliquidID = 2106; + // int SnliquidID = 2106; - bool invert_at_setup = true; + // bool invert_at_setup = true; - singularity::EOS Snbeta = singularity::EOSPAC(SnbetaID); - singularity::EOS Sngamma = singularity::EOSPAC(SngammaID); - // singularity::EOS Sndelta = singularity::EOSPAC(SndeltaID); - singularity::EOS Snhcp = singularity::EOSPAC(SnhcpID); - // singularity::EOS Snliquid = singularity::EOSPAC(SnliquidID); + EOS Snbeta = singularity::EOSPAC(SnbetaID); + EOS Sngamma = singularity::EOSPAC(SngammaID); + // EOS Sndelta = singularity::EOSPAC(SndeltaID); + EOS Snhcp = singularity::EOSPAC(SnhcpID); + // EOS Snliquid = singularity::EOSPAC(SnliquidID); eos[0] = Snbeta.GetOnDevice(); eos[1] = Sngamma.GetOnDevice(); diff --git a/test/pte_test_first.hpp b/test/pte_test_first.hpp index 77d6dbb6901..f2815e1e2de 100644 --- a/test/pte_test_first.hpp +++ b/test/pte_test_first.hpp @@ -19,20 +19,27 @@ #include #include -#include +#include +#include constexpr int NMAT = 3; constexpr int NTRIAL = 100; constexpr int NPTS = NTRIAL * NMAT; constexpr int HIST_SIZE = 10; +using singularity::Variant; +using singularity::Gruneisen; +using singularity::DavisReactants; +using singularity::DavisProducts; +using EOS = Variant; + template inline void set_eos(T *eos) { - singularity::EOS gr = singularity::Gruneisen(394000.0, 1.489, 0.0, 0.0, 2.02, 0.47, + EOS gr = singularity::Gruneisen(394000.0, 1.489, 0.0, 0.0, 2.02, 0.47, 8.93, 297.0, 1.0e6, 0.383e7); - singularity::EOS dr = singularity::DavisReactants( + EOS dr = singularity::DavisReactants( 1.890, 4.115e10, 1.0e6, 297.0, 1.8e5, 4.6, 0.34, 0.56, 0.0, 0.4265, 0.001074e10); - singularity::EOS dp = singularity::DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419, + EOS dp = singularity::DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419, 3.2e10, 0.001072e10); eos[0] = gr.GetOnDevice(); eos[1] = dr.GetOnDevice(); diff --git a/test/test_pte.cpp b/test/test_pte.cpp index 44eb859241d..23fa5955589 100644 --- a/test/test_pte.cpp +++ b/test/test_pte.cpp @@ -24,12 +24,14 @@ #include #include #include -#include #include -using namespace singularity; +#include +#include using DataBox = Spiner::DataBox; +using singularity::PTESolverRhoTRequiredScratch; +using singularity::PTESolverRhoT; int main(int argc, char *argv[]) { diff --git a/test/test_pte_2phase.cpp b/test/test_pte_2phase.cpp index 9fa55a6007b..197c90292aa 100644 --- a/test/test_pte_2phase.cpp +++ b/test/test_pte_2phase.cpp @@ -21,16 +21,21 @@ #include #include -#include -#include #include #include -#include #include -using namespace singularity; +#include +#include + +// these two headers share a variant +using EOS = singularity::Variant; +#include +#include using namespace pte_longtest_2phaseVinetSn; +using singularity::PTESolverRhoTRequiredScratch; +using singularity::PTESolverRhoT; using DataBox = Spiner::DataBox; @@ -180,8 +185,9 @@ int main(int argc, char *argv[]) { const Real Tguess = ApproxTemperatureFromRhoMatU(NMAT, eos, rho_tot * sie_tot, rho, vfrac); - - std::cout << "Tguess " << Tguess << std::endl; + if (t == 0) { + printf("Tguess %.14e\n", Tguess); + } auto method = PTESolverRhoT, decltype(lambda)>( diff --git a/test/test_pte_3phase.cpp b/test/test_pte_3phase.cpp index a605985d10c..09374e846d4 100644 --- a/test/test_pte_3phase.cpp +++ b/test/test_pte_3phase.cpp @@ -21,15 +21,19 @@ #include #include -#include #include #include -#include #include -using namespace singularity; +#include +#include + +using EOS = singularity::Variant; +#include using DataBox = Spiner::DataBox; +using singularity::PTESolverRhoTRequiredScratch; +using singularity::PTESolverRhoT; int main(int argc, char *argv[]) { @@ -37,7 +41,8 @@ int main(int argc, char *argv[]) { #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::initialize(); #endif - { +// JMM: EOSPAC tests do not work on device. +if constexpr (PortsOfCall::EXECUTION_IS_HOST) { // EOS #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::View eos_v("eos", NMAT); @@ -54,7 +59,6 @@ int main(int argc, char *argv[]) { #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::deep_copy(eos_v, eos_hv); #endif - using EOSAccessor = LinearIndexer; EOSAccessor eos(eos_v); @@ -109,6 +113,8 @@ int main(int argc, char *argv[]) { #endif // setup state + printf("pte_3phase setup state\n"); + srand(time(NULL)); for (int n = 0; n < NTRIAL; n++) { Indexer2D r(n, rho_hm); @@ -127,6 +133,7 @@ int main(int argc, char *argv[]) { Kokkos::deep_copy(press_v, press_vh); Kokkos::deep_copy(hist_d, hist_vh); #endif + printf("pte_3phase state set\n"); #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::View nsuccess_d("n successes"); @@ -177,8 +184,9 @@ int main(int argc, char *argv[]) { const Real Tguess = ApproxTemperatureFromRhoMatU(NMAT, eos, rho_tot * sie_tot, rho, vfrac); - - std::cout << "Tguess " << Tguess << std::endl; + if (t == 0) { + printf("Tguess %.14e\n", Tguess); + } auto method = PTESolverRhoT, decltype(lambda)>( @@ -253,5 +261,9 @@ int main(int argc, char *argv[]) { #endif // poor-man's ctest integration + if constexpr (PortsOfCall::EXECUTION_IS_HOST) { return (nsuccess >= 0.5 * NTRIAL) ? 0 : 1; + } else { + return 0; + } } From de3dc040eacb0f2dc6d4cee25cae7557eb543214 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 6 Sep 2024 14:03:04 -0600 Subject: [PATCH 18/20] formatting --- test/pte_longtest_2phaseVinetSn.hpp | 10 ++++------ test/pte_test_2phaseVinetSn.hpp | 10 ++++------ test/pte_test_first.hpp | 18 +++++++++--------- test/test_pte.cpp | 2 +- test/test_pte_2phase.cpp | 8 ++++---- test/test_pte_3phase.cpp | 20 ++++++++++---------- 6 files changed, 32 insertions(+), 36 deletions(-) diff --git a/test/pte_longtest_2phaseVinetSn.hpp b/test/pte_longtest_2phaseVinetSn.hpp index eac2ce665a8..51e66c74f34 100644 --- a/test/pte_longtest_2phaseVinetSn.hpp +++ b/test/pte_longtest_2phaseVinetSn.hpp @@ -83,12 +83,10 @@ inline void set_eos(T *eos) { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; - EOS Snbeta = - singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, 0.658e09, - 0.4419e07, d2to40); - EOS Sngamma = - singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, 0.2161e07, - 1.025e09, 0.5051e07, d2to40); + EOS Snbeta = singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, + 0.658e09, 0.4419e07, d2to40); + EOS Sngamma = singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, + 0.2161e07, 1.025e09, 0.5051e07, d2to40); eos[0] = Snbeta.GetOnDevice(); eos[1] = Sngamma.GetOnDevice(); return; diff --git a/test/pte_test_2phaseVinetSn.hpp b/test/pte_test_2phaseVinetSn.hpp index efde8919f37..881c0c5ee04 100644 --- a/test/pte_test_2phaseVinetSn.hpp +++ b/test/pte_test_2phaseVinetSn.hpp @@ -59,12 +59,10 @@ inline void set_eos(T *eos) { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; - EOS Snbeta = - singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, 0.658e09, - 0.4419e07, d2to40); - EOS Sngamma = - singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, 0.2161e07, - 1.025e09, 0.5051e07, d2to40); + EOS Snbeta = singularity::Vinet(7.285, 298.0, 0.529e12, 5.3345, 0.000072977, 0.2149e07, + 0.658e09, 0.4419e07, d2to40); + EOS Sngamma = singularity::Vinet(7.271, 298.0, 0.3878e12, 6.0532, 0.0001085405, + 0.2161e07, 1.025e09, 0.5051e07, d2to40); eos[0] = Snbeta.GetOnDevice(); eos[1] = Sngamma.GetOnDevice(); return; diff --git a/test/pte_test_first.hpp b/test/pte_test_first.hpp index f2815e1e2de..4f975596306 100644 --- a/test/pte_test_first.hpp +++ b/test/pte_test_first.hpp @@ -27,20 +27,20 @@ constexpr int NTRIAL = 100; constexpr int NPTS = NTRIAL * NMAT; constexpr int HIST_SIZE = 10; -using singularity::Variant; -using singularity::Gruneisen; -using singularity::DavisReactants; using singularity::DavisProducts; +using singularity::DavisReactants; +using singularity::Gruneisen; +using singularity::Variant; using EOS = Variant; template inline void set_eos(T *eos) { - EOS gr = singularity::Gruneisen(394000.0, 1.489, 0.0, 0.0, 2.02, 0.47, - 8.93, 297.0, 1.0e6, 0.383e7); - EOS dr = singularity::DavisReactants( - 1.890, 4.115e10, 1.0e6, 297.0, 1.8e5, 4.6, 0.34, 0.56, 0.0, 0.4265, 0.001074e10); - EOS dp = singularity::DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419, - 3.2e10, 0.001072e10); + EOS gr = singularity::Gruneisen(394000.0, 1.489, 0.0, 0.0, 2.02, 0.47, 8.93, 297.0, + 1.0e6, 0.383e7); + EOS dr = singularity::DavisReactants(1.890, 4.115e10, 1.0e6, 297.0, 1.8e5, 4.6, 0.34, + 0.56, 0.0, 0.4265, 0.001074e10); + EOS dp = singularity::DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419, 3.2e10, + 0.001072e10); eos[0] = gr.GetOnDevice(); eos[1] = dr.GetOnDevice(); eos[2] = dp.GetOnDevice(); diff --git a/test/test_pte.cpp b/test/test_pte.cpp index 23fa5955589..367a739784d 100644 --- a/test/test_pte.cpp +++ b/test/test_pte.cpp @@ -30,8 +30,8 @@ #include using DataBox = Spiner::DataBox; -using singularity::PTESolverRhoTRequiredScratch; using singularity::PTESolverRhoT; +using singularity::PTESolverRhoTRequiredScratch; int main(int argc, char *argv[]) { diff --git a/test/test_pte_2phase.cpp b/test/test_pte_2phase.cpp index 197c90292aa..0a759fc27f4 100644 --- a/test/test_pte_2phase.cpp +++ b/test/test_pte_2phase.cpp @@ -34,8 +34,8 @@ using EOS = singularity::Variant; #include using namespace pte_longtest_2phaseVinetSn; -using singularity::PTESolverRhoTRequiredScratch; using singularity::PTESolverRhoT; +using singularity::PTESolverRhoTRequiredScratch; using DataBox = Spiner::DataBox; @@ -185,9 +185,9 @@ int main(int argc, char *argv[]) { const Real Tguess = ApproxTemperatureFromRhoMatU(NMAT, eos, rho_tot * sie_tot, rho, vfrac); - if (t == 0) { - printf("Tguess %.14e\n", Tguess); - } + if (t == 0) { + printf("Tguess %.14e\n", Tguess); + } auto method = PTESolverRhoT, decltype(lambda)>( diff --git a/test/test_pte_3phase.cpp b/test/test_pte_3phase.cpp index 09374e846d4..c4a89808b97 100644 --- a/test/test_pte_3phase.cpp +++ b/test/test_pte_3phase.cpp @@ -32,8 +32,8 @@ using EOS = singularity::Variant; #include using DataBox = Spiner::DataBox; -using singularity::PTESolverRhoTRequiredScratch; using singularity::PTESolverRhoT; +using singularity::PTESolverRhoTRequiredScratch; int main(int argc, char *argv[]) { @@ -41,8 +41,8 @@ int main(int argc, char *argv[]) { #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::initialize(); #endif -// JMM: EOSPAC tests do not work on device. -if constexpr (PortsOfCall::EXECUTION_IS_HOST) { + // JMM: EOSPAC tests do not work on device. + if constexpr (PortsOfCall::EXECUTION_IS_HOST) { // EOS #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::View eos_v("eos", NMAT); @@ -113,7 +113,7 @@ if constexpr (PortsOfCall::EXECUTION_IS_HOST) { #endif // setup state - printf("pte_3phase setup state\n"); + printf("pte_3phase setup state\n"); srand(time(NULL)); for (int n = 0; n < NTRIAL; n++) { @@ -133,7 +133,7 @@ if constexpr (PortsOfCall::EXECUTION_IS_HOST) { Kokkos::deep_copy(press_v, press_vh); Kokkos::deep_copy(hist_d, hist_vh); #endif - printf("pte_3phase state set\n"); + printf("pte_3phase state set\n"); #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::View nsuccess_d("n successes"); @@ -184,9 +184,9 @@ if constexpr (PortsOfCall::EXECUTION_IS_HOST) { const Real Tguess = ApproxTemperatureFromRhoMatU(NMAT, eos, rho_tot * sie_tot, rho, vfrac); - if (t == 0) { - printf("Tguess %.14e\n", Tguess); - } + if (t == 0) { + printf("Tguess %.14e\n", Tguess); + } auto method = PTESolverRhoT, decltype(lambda)>( @@ -262,8 +262,8 @@ if constexpr (PortsOfCall::EXECUTION_IS_HOST) { // poor-man's ctest integration if constexpr (PortsOfCall::EXECUTION_IS_HOST) { - return (nsuccess >= 0.5 * NTRIAL) ? 0 : 1; + return (nsuccess >= 0.5 * NTRIAL) ? 0 : 1; } else { - return 0; + return 0; } } From dcb5f752676cc990cffdeffed4570adacad20eff Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 9 Sep 2024 16:07:39 -0600 Subject: [PATCH 19/20] try using eos-developmental --- test/pte_test_5phaseSesameSn.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/test/pte_test_5phaseSesameSn.hpp b/test/pte_test_5phaseSesameSn.hpp index 723296a7ba0..5d23afe6203 100644 --- a/test/pte_test_5phaseSesameSn.hpp +++ b/test/pte_test_5phaseSesameSn.hpp @@ -69,9 +69,8 @@ using EOS = Variant; template inline void set_eos(T *eos) { - int sl = - symlink("/projects/shavano/dev/aematts/SEOS/PTEsolver_test_cases/sn2162-v01.bin", - "sesameu"); + int sl = symlink("/usr/projects/data/eos/eos-developmental/Sn2162/v01/sn2162-v01.bin", + "sesameu"); int SnbetaID = 2102; int SngammaID = 2103; From c3517ac5b698f68d587450b67dbbffb2d5011b75 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 9 Sep 2024 16:30:02 -0600 Subject: [PATCH 20/20] upload artifact version --- .github/workflows/draft-pdf.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml index ceeed6452a7..d39204208cd 100644 --- a/.github/workflows/draft-pdf.yml +++ b/.github/workflows/draft-pdf.yml @@ -23,7 +23,7 @@ jobs: # This should be the path to the paper within your repo. paper-path: joss-paper/paper.md - name: Upload - uses: actions/upload-artifact@v1 + uses: actions/upload-artifact@v4 with: name: paper # This is the output path where Pandoc will write the compiled