From ef8815425c24679b77c3dfbfb78a078aa2a8a375 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 27 Jun 2024 11:41:09 -0600 Subject: [PATCH 01/61] Simplify line search and add comments --- singularity-eos/closure/mixed_cell_models.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index d1ca6f89db..4da20aa442 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -1426,16 +1426,16 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { err = s.TestUpdate(scale); if (err > err_old + line_search_alpha * gradfdx) { // backtrack - Real err_mid = s.TestUpdate(0.5); + scale = 0.5; + Real err_mid = s.TestUpdate(scale); if (err_mid < err && err_mid < err_old) { + // try a larger step since a half step reduces the error scale = 0.75 + 0.5 * robust::ratio(err_mid - err, err - 2.0 * err_mid + err_old); - } else { - scale = line_search_fac; } - for (int line_iter = 0; line_iter < line_search_max_iter; line_iter++) { err = s.TestUpdate(scale); if (err < err_old + line_search_alpha * scale * gradfdx) break; + // shrink the step if the error isn't reduced scale *= line_search_fac; } } From 8b99c75bebb6f7108f64e23a1785890cd9e724e0 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 27 Jun 2024 12:37:14 -0600 Subject: [PATCH 02/61] Make sure the scale that comes out is the minimum over all checks --- singularity-eos/closure/mixed_cell_models.hpp | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 4da20aa442..d4e9d902eb 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -675,7 +675,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); + scale = std::min(scale, -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m])); } } const Real Tnew = Tequil + scale * dx[nmat]; @@ -685,12 +685,12 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = robust::ratio(0.5 * alpha_max - vfrac[m], dx[m]); + scale = std::min(scale, robust::ratio(0.5 * alpha_max - vfrac[m], dx[m])); } } // control how big of a step toward T = 0 is allowed if (scale * dx[nmat] < -0.95 * Tequil) { - scale = robust::ratio(-0.95 * Tequil, dx[nmat]); + scale = std::min(scale, robust::ratio(-0.95 * Tequil, dx[nmat])); } // Now apply the overall scaling for (int i = 0; i < neq; ++i) @@ -883,7 +883,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = robust::ratio(-vfrac_safety_fac * vfrac[m], dx[m]); + scale = std::min(scale, robust::ratio(-vfrac_safety_fac * vfrac[m], dx[m])); } } // control how big of a step toward rho = rho(Pmin) is allowed @@ -891,7 +891,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase const Real rho_min = eos[m].RhoPmin(Tequil); const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]); + scale = std::min(scale, robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m])); } } // Now apply the overall scaling @@ -1097,7 +1097,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); + scale = std::min(scale, -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m])); } } const Real Tnew = Tequil + scale * dx[nmat]; @@ -1107,12 +1107,12 @@ class PTESolverFixedP : public mix_impl::PTESolverBase std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); + scale = std::min(scale, 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m])); } } // control how big of a step toward T = 0 is allowed if (scale * dx[nmat] < -0.95 * Tequil) { - scale = -0.95 * robust::ratio(Tequil, dx[nmat]); + scale = std::min(scale, -0.95 * robust::ratio(Tequil, dx[nmat])); } // Now apply the overall scaling for (int i = 0; i < neq; ++i) @@ -1336,12 +1336,12 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { for (int m = 0; m < nmat; ++m) { // control how big of a step toward vfrac = 0 is allowed if (scale * dx[m] < -0.1 * vfrac[m]) { - scale = -0.1 * robust::ratio(vfrac[m], dx[m]); + scale = std::min(scale, -0.1 * robust::ratio(vfrac[m], dx[m])); } // try to control steps toward T = 0 const Real dt = (dtdv[m] * dx[m] + dtde[m] * dx[m + nmat]); if (scale * dt < -0.1 * temp[m]) { - scale = -0.1 * robust::ratio(temp[m], dt); + scale = std::min(scale, -0.1 * robust::ratio(temp[m], dt)); } const Real tt = temp[m] + scale * dt; const Real rho_min = @@ -1349,7 +1349,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { const Real alpha_max = robust::ratio(rhobar[m], rho_min); // control how big of a step toward rho = rho(Pmin) is allowed if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); + scale = std::min(scale, 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m])); } } // Now apply the overall scaling From f60e4b4836823734cf132275d99377cab6e2fb4c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 27 Jun 2024 13:09:58 -0600 Subject: [PATCH 03/61] Remove commented variable and add clarifying comment --- singularity-eos/closure/mixed_cell_models.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index d4e9d902eb..d21822db52 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -1417,11 +1417,10 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { // possibly scale the update to stay within reasonable bounds Real scale = s.ScaleDx(); - // const Real scale_save = scale; // Line search Real gradfdx = -2.0 * scale * err; - scale = 1.0; + scale = 1.0; // New scale for line search Real err_old = err; err = s.TestUpdate(scale); if (err > err_old + line_search_alpha * gradfdx) { From f501b2aae6e32a52d527e1039c77b2fa5d869521 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 27 Jun 2024 13:13:55 -0600 Subject: [PATCH 04/61] Clang format --- singularity-eos/closure/mixed_cell_models.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index d21822db52..f079b04fb4 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -1420,7 +1420,7 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { // Line search Real gradfdx = -2.0 * scale * err; - scale = 1.0; // New scale for line search + scale = 1.0; // New scale for line search Real err_old = err; err = s.TestUpdate(scale); if (err > err_old + line_search_alpha * gradfdx) { From ca3d16a769def4a651ac9af60eb7b95da1d567a1 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 27 Jun 2024 18:51:58 -0600 Subject: [PATCH 05/61] Change behavior to assume small mass fraction materials are at the cell sie and density --- singularity-eos/eos/get_sg_eos_functors.hpp | 18 ++++++++++++------ singularity-eos/eos/get_sg_eos_rho_e.cpp | 5 +++-- singularity-eos/eos/get_sg_eos_rho_p.cpp | 5 +++-- singularity-eos/eos/get_sg_eos_rho_t.cpp | 5 +++-- 4 files changed, 21 insertions(+), 12 deletions(-) diff --git a/singularity-eos/eos/get_sg_eos_functors.hpp b/singularity-eos/eos/get_sg_eos_functors.hpp index 8bfa21138c..95cf22b745 100644 --- a/singularity-eos/eos/get_sg_eos_functors.hpp +++ b/singularity-eos/eos/get_sg_eos_functors.hpp @@ -62,7 +62,7 @@ struct init_functor { mass_frac_cutoff_} {} PORTABLE_INLINE_FUNCTION - void operator()(const int i, const int tid, double &mass_sum, int &npte, + void operator()(const int i, const int tid, double &mass_sum, int &npte, double &vfrac_sum, const Real t_mult, const Real s_mult, const Real p_mult) const { /* normalize mass fractions */ /* first find the mass sum */ @@ -71,13 +71,19 @@ struct init_functor { for (int m = 0; m < nmat; ++m) { mass_sum += frac_mass_v(i, m); pte_idxs(tid, m) = eos_offsets_v(m) - 1; - frac_vol_v(i, m) = 0.0; + // Assume non-participating materials are all at cell density. PTE volumes + // will be overwritten + frac_vol_v(i, m) = frac_mass_v(i, m); + // Initialize all materials to the cell sie (if provided). PTE sie values + // will be overwritten. This also means that the sie for the PTE + // materials is the same as the cell PTE after removing small materials + frac_ie_v(i, m) = sie_v(i) * frac_mass_v(i, m) * s_mult; } for (int m = 0; m < nmat; ++m) { frac_mass_v(i, m) /= mass_sum; } check_all_vals(i); - // count the number of participating materials and zero the inputs + // count the number of participating materials npte = 0; for (int m = 0; m < nmat; ++m) { if (frac_mass_v(i, m) > mass_frac_cutoff) { @@ -85,10 +91,8 @@ struct init_functor { pte_idxs(tid, npte) = eos_offsets_v(m) - 1; pte_mats(tid, npte) = m; npte += 1; - } else { - frac_ie_v(i, m) = 0.0; } - // zero the inputs + // zero the PTE solver material inputs vfrac_pte(tid, m) = 0.0; sie_pte(tid, m) = 0.0; temp_pte(tid, m) = 0.0; @@ -98,10 +102,12 @@ struct init_functor { // NOTE: the volume fractions and densities need to be consistent with the // total specific volume since they are used to calculate internal // quantities for the PTE solver + vfrac_sum = 0.; for (int mp = 0; mp < npte; ++mp) { const int m = pte_mats(tid, mp); // Need to guess volume fractions vfrac_pte(tid, mp) = frac_mass_v(i, m); + vfrac_sum += vfrac_pte(tid, mp); // Calculate densities to be consistent with these volume fractions rho_pte(tid, mp) = frac_mass_v(i, m) / spvol_v(i) / vfrac_pte(tid, mp); temp_pte(tid, mp) = temp_v(i) * ev2k * t_mult; diff --git a/singularity-eos/eos/get_sg_eos_rho_e.cpp b/singularity-eos/eos/get_sg_eos_rho_e.cpp index 84db423ad3..9518b3176b 100644 --- a/singularity-eos/eos/get_sg_eos_rho_e.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_e.cpp @@ -37,8 +37,9 @@ void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v, const int32_t tid{small_loop ? iloop : token}; double mass_sum{0.0}; int npte{0}; + double vfrac_sum{0.0}; // initialize values for solver / lookup - i_func(i, tid, mass_sum, npte, 0.0, 1.0, 0.0); + i_func(i, tid, mass_sum, npte, vfrac_sum, 0.0, 1.0, 0.0); // need to initialize the scratch before it's used to avoid undefined behavior for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { solver_scratch(tid, idx) = 0.0; @@ -54,7 +55,7 @@ void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v, singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0)); // reset inputs PTESolverRhoT method( - npte, eos_inx, 1.0, sie_v(i), &rho_pte(tid, 0), &vfrac_pte(tid, 0), + npte, eos_inx, vfrac_sum, sie_v(i), &rho_pte(tid, 0), &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache, &solver_scratch(tid, 0)); pte_converged = PTESolver(method); diff --git a/singularity-eos/eos/get_sg_eos_rho_p.cpp b/singularity-eos/eos/get_sg_eos_rho_p.cpp index ea376b62b1..73a349b145 100644 --- a/singularity-eos/eos/get_sg_eos_rho_p.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_p.cpp @@ -38,8 +38,9 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v, const int32_t tid{small_loop ? iloop : token}; double mass_sum{0.0}; int npte{0}; + double vfrac_sum{0.0}; // initialize values for solver / lookup - i_func(i, tid, mass_sum, npte, 0.0, 0.0, 1.0); + i_func(i, tid, mass_sum, npte, vfrac_sum, 0.0, 0.0, 1.0); Real sie_tot_true{0.0}; // need to initialize the scratch before it's used to avoid undefined behavior for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { @@ -54,7 +55,7 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v, // eos accessor singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0)); PTESolverFixedP method( - npte, eos_inx, 1.0, press_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0), + npte, eos_inx, vfrac_sum, press_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache[0], &solver_scratch(tid, 0)); pte_converged = PTESolver(method); diff --git a/singularity-eos/eos/get_sg_eos_rho_t.cpp b/singularity-eos/eos/get_sg_eos_rho_t.cpp index df3d14c50b..313c388cd3 100644 --- a/singularity-eos/eos/get_sg_eos_rho_t.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_t.cpp @@ -38,8 +38,9 @@ void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v, const int32_t tid{small_loop ? iloop : token}; double mass_sum{0.0}; int npte{0}; + double vfrac_sum{0.0}; // initialize values for solver / lookup - i_func(i, tid, mass_sum, npte, 1.0, 0.0, 0.0); + i_func(i, tid, mass_sum, npte, vfrac_sum, 1.0, 0.0, 0.0); // calculate pte condition (lookup for 1 mat cell) Real sie_tot_true{0.0}; // need to initialize the scratch before it's used to avoid undefined behavior @@ -55,7 +56,7 @@ void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v, // eos accessor singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0)); PTESolverFixedT method( - npte, eos_inx, 1.0, temp_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0), + npte, eos_inx, vfrac_sum, temp_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache, &solver_scratch(tid, 0)); pte_converged = PTESolver(method); From a853cbc1aaa60932cafb043a4bd0bc907499abf7 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 27 Jun 2024 18:55:26 -0600 Subject: [PATCH 06/61] Clang format --- singularity-eos/eos/get_sg_eos_functors.hpp | 5 +++-- singularity-eos/eos/get_sg_eos_rho_p.cpp | 6 +++--- singularity-eos/eos/get_sg_eos_rho_t.cpp | 6 +++--- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/singularity-eos/eos/get_sg_eos_functors.hpp b/singularity-eos/eos/get_sg_eos_functors.hpp index 95cf22b745..d8633b2046 100644 --- a/singularity-eos/eos/get_sg_eos_functors.hpp +++ b/singularity-eos/eos/get_sg_eos_functors.hpp @@ -62,8 +62,9 @@ struct init_functor { mass_frac_cutoff_} {} PORTABLE_INLINE_FUNCTION - void operator()(const int i, const int tid, double &mass_sum, int &npte, double &vfrac_sum, - const Real t_mult, const Real s_mult, const Real p_mult) const { + void operator()(const int i, const int tid, double &mass_sum, int &npte, + double &vfrac_sum, const Real t_mult, const Real s_mult, + const Real p_mult) const { /* normalize mass fractions */ /* first find the mass sum */ /* also set idxs as the decrement of the eos offsets */ diff --git a/singularity-eos/eos/get_sg_eos_rho_p.cpp b/singularity-eos/eos/get_sg_eos_rho_p.cpp index 73a349b145..f735ba83ca 100644 --- a/singularity-eos/eos/get_sg_eos_rho_p.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_p.cpp @@ -55,9 +55,9 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v, // eos accessor singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0)); PTESolverFixedP method( - npte, eos_inx, vfrac_sum, press_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0), - &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache[0], - &solver_scratch(tid, 0)); + npte, eos_inx, vfrac_sum, press_pte(tid, 0), &rho_pte(tid, 0), + &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), + cache[0], &solver_scratch(tid, 0)); pte_converged = PTESolver(method); // calculate total sie for (int mp = 0; mp < npte; ++mp) { diff --git a/singularity-eos/eos/get_sg_eos_rho_t.cpp b/singularity-eos/eos/get_sg_eos_rho_t.cpp index 313c388cd3..cd292c811a 100644 --- a/singularity-eos/eos/get_sg_eos_rho_t.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_t.cpp @@ -56,9 +56,9 @@ void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v, // eos accessor singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0)); PTESolverFixedT method( - npte, eos_inx, vfrac_sum, temp_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0), - &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache, - &solver_scratch(tid, 0)); + npte, eos_inx, vfrac_sum, temp_pte(tid, 0), &rho_pte(tid, 0), + &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), + cache, &solver_scratch(tid, 0)); pte_converged = PTESolver(method); // calculate total internal energy for (int mp = 0; mp < npte; ++mp) { From e8851136b203c549914edb6135eb7dc48d7c808f Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 1 Jul 2024 18:31:42 -0600 Subject: [PATCH 07/61] min not needed since scales are checked sequentially _after_ being modified --- singularity-eos/closure/mixed_cell_models.hpp | 24 +++++++++---------- 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index f079b04fb4..13f6ecb63a 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -675,7 +675,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = std::min(scale, -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m])); + scale = -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); } } const Real Tnew = Tequil + scale * dx[nmat]; @@ -685,12 +685,12 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = std::min(scale, robust::ratio(0.5 * alpha_max - vfrac[m], dx[m])); + scale = robust::ratio(0.5 * alpha_max - vfrac[m], dx[m]); } } // control how big of a step toward T = 0 is allowed if (scale * dx[nmat] < -0.95 * Tequil) { - scale = std::min(scale, robust::ratio(-0.95 * Tequil, dx[nmat])); + scale = robust::ratio(-0.95 * Tequil, dx[nmat]); } // Now apply the overall scaling for (int i = 0; i < neq; ++i) @@ -883,7 +883,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = std::min(scale, robust::ratio(-vfrac_safety_fac * vfrac[m], dx[m])); + scale = robust::ratio(-vfrac_safety_fac * vfrac[m], dx[m]); } } // control how big of a step toward rho = rho(Pmin) is allowed @@ -891,7 +891,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase const Real rho_min = eos[m].RhoPmin(Tequil); const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = std::min(scale, robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m])); + scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]); } } // Now apply the overall scaling @@ -1097,7 +1097,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = std::min(scale, -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m])); + scale = -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); } } const Real Tnew = Tequil + scale * dx[nmat]; @@ -1107,12 +1107,12 @@ class PTESolverFixedP : public mix_impl::PTESolverBase std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); const Real alpha_max = robust::ratio(rhobar[m], rho_min); if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = std::min(scale, 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m])); + scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); } } // control how big of a step toward T = 0 is allowed if (scale * dx[nmat] < -0.95 * Tequil) { - scale = std::min(scale, -0.95 * robust::ratio(Tequil, dx[nmat])); + scale = -0.95 * robust::ratio(Tequil, dx[nmat]); } // Now apply the overall scaling for (int i = 0; i < neq; ++i) @@ -1336,12 +1336,12 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { for (int m = 0; m < nmat; ++m) { // control how big of a step toward vfrac = 0 is allowed if (scale * dx[m] < -0.1 * vfrac[m]) { - scale = std::min(scale, -0.1 * robust::ratio(vfrac[m], dx[m])); + scale = -0.1 * robust::ratio(vfrac[m], dx[m]); } // try to control steps toward T = 0 const Real dt = (dtdv[m] * dx[m] + dtde[m] * dx[m + nmat]); if (scale * dt < -0.1 * temp[m]) { - scale = std::min(scale, -0.1 * robust::ratio(temp[m], dt)); + scale = -0.1 * robust::ratio(temp[m], dt); } const Real tt = temp[m] + scale * dt; const Real rho_min = @@ -1349,7 +1349,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { const Real alpha_max = robust::ratio(rhobar[m], rho_min); // control how big of a step toward rho = rho(Pmin) is allowed if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = std::min(scale, 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m])); + scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); } } // Now apply the overall scaling @@ -1424,7 +1424,7 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { Real err_old = err; err = s.TestUpdate(scale); if (err > err_old + line_search_alpha * gradfdx) { - // backtrack + // backtrack to middle of step scale = 0.5; Real err_mid = s.TestUpdate(scale); if (err_mid < err && err_mid < err_old) { From 71250a5419e94a1d3c3717723780327d2ae32896 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 1 Jul 2024 18:56:41 -0600 Subject: [PATCH 08/61] Add sanity check that initial volume fractions add up to 1 or less --- singularity-eos/eos/get_sg_eos_functors.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/singularity-eos/eos/get_sg_eos_functors.hpp b/singularity-eos/eos/get_sg_eos_functors.hpp index d8633b2046..0561617500 100644 --- a/singularity-eos/eos/get_sg_eos_functors.hpp +++ b/singularity-eos/eos/get_sg_eos_functors.hpp @@ -115,6 +115,7 @@ struct init_functor { press_pte(tid, mp) = press_v(i) * p_mult; sie_pte(tid, mp) = sie_v(i) * frac_mass_v(i, m) * s_mult; } + PORTABLE_REQUIRE(vfrac_sum <= 1.0, "Volume fraction sum is greater than 1"); return; } From a4e1f0c1c5ab48ee08e107ea6f243ffd8baa3ec2 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 1 Jul 2024 19:00:02 -0600 Subject: [PATCH 09/61] Add more asserts and label what is going on if an assert is tripped --- singularity-eos/closure/mixed_cell_models.hpp | 23 ++++++++++++------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 13f6ecb63a..4557e0109f 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -170,9 +170,11 @@ class PTESolverBase { PORTABLE_INLINE_FUNCTION virtual void Fixup() const { Real vsum = 0; - for (int m = 0; m < nmat; ++m) + for (int m = 0; m < nmat; ++m) { + PORTABLE_REQUIRE(vfrac[m] > 0, "Negative volume fraction in Fixup"); vsum += vfrac[m]; - PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive"); + } + PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive in Fixup"); for (int m = 0; m < nmat; ++m) vfrac[m] *= robust::ratio(vfrac_total, vsum); } @@ -182,7 +184,7 @@ class PTESolverBase { void Finalize() { for (int m = 0; m < nmat; m++) { temp[m] *= Tnorm; - PORTABLE_REQUIRE(temp[m] >= 0., "Non-positive temperature returned"); + PORTABLE_REQUIRE(temp[m] >= 0., "Non-positive temperature returned in Finalize"); u[m] *= uscale; press[m] *= uscale; } @@ -219,8 +221,9 @@ class PTESolverBase { // material m averaged over the full PTE volume rho_total = 0.0; for (int m = 0; m < nmat; ++m) { - PORTABLE_REQUIRE(vfrac[m] > 0., "Non-positive volume fraction"); - PORTABLE_REQUIRE(rho[m] > 0., "Non-positive density"); + PORTABLE_REQUIRE(vfrac[m] > 0., + "Non-positive volume fraction provided to PTE solver"); + PORTABLE_REQUIRE(rho[m] > 0., "Non-positive density provided to PTE solver"); rhobar[m] = rho[m] * vfrac[m]; rho_total += rhobar[m]; } @@ -353,12 +356,12 @@ class PTESolverBase { Asum += vfrac[m] * robust::ratio(press[m], temp[m]); rhoBsum += rho[m] * vfrac[m] * robust::ratio(sie[m], temp[m]); } - PORTABLE_REQUIRE(Tnorm > 0., "Non-positive temperature guess"); + PORTABLE_REQUIRE(Tnorm > 0., "Non-positive Ideal PTE temperature guess"); Asum *= uscale / Tnorm; rhoBsum /= Tnorm; - PORTABLE_REQUIRE(rhoBsum > 0., "Non-positive energy density"); + PORTABLE_REQUIRE(rhoBsum > 0., "Non-positive Ideal PTE energy density guess"); Tideal = uscale / rhoBsum / Tnorm; - PORTABLE_REQUIRE(vfrac_total > 0., "Non-positive volume fraction sum"); + PORTABLE_REQUIRE(vfrac_total > 0., "Non-positive Ideal PTE volume fraction sum"); Pideal = robust::ratio(Tnorm * Tideal * Asum, uscale * vfrac_total); } @@ -629,6 +632,8 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { const Real vf_pert = vfrac[m] + dv; const Real rho_pert = robust::ratio(rhobar[m], vf_pert); + PORTABLE_REQUIRE(vfrac[m] > 0, "Negative volume fraction in rho-T PTE iteration"); + Real p_pert{}; Real e_pert = eos[m].InternalEnergyFromDensityTemperature(rho_pert, Tnorm * Tequil, Cache[m]); @@ -708,6 +713,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { vtemp[m] = vfrac[m]; } Tequil = Ttemp + scale * dx[nmat]; + PORTABLE_REQUIRE(Tequil >= 0., "Negative temperature update in rho-T PTE solver iteration"); for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; rho[m] = robust::ratio(rhobar[m], vfrac[m]); @@ -1130,6 +1136,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase vtemp[m] = vfrac[m]; } Tequil = Ttemp + scale * dx[nmat]; + PORTABLE_REQUIRE(Tequil >= 0., "Negative temperature in Fixed P PTE solver iteration"); for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; rho[m] = robust::ratio(rhobar[m], vfrac[m]); From 5e947cd6ba691ce50d7c5432279318f2ec458331 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 1 Jul 2024 20:08:56 -0600 Subject: [PATCH 10/61] Remove volume fraction sum check since it's really not relevant and produces read herrings --- singularity-eos/eos/get_sg_eos_functors.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/singularity-eos/eos/get_sg_eos_functors.hpp b/singularity-eos/eos/get_sg_eos_functors.hpp index 0561617500..d8633b2046 100644 --- a/singularity-eos/eos/get_sg_eos_functors.hpp +++ b/singularity-eos/eos/get_sg_eos_functors.hpp @@ -115,7 +115,6 @@ struct init_functor { press_pte(tid, mp) = press_v(i) * p_mult; sie_pte(tid, mp) = sie_v(i) * frac_mass_v(i, m) * s_mult; } - PORTABLE_REQUIRE(vfrac_sum <= 1.0, "Volume fraction sum is greater than 1"); return; } From 9b3f9003766a3a00b865e5355280fd52fd2bfece Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 1 Jul 2024 21:00:02 -0600 Subject: [PATCH 11/61] Clang format --- singularity-eos/closure/mixed_cell_models.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 4557e0109f..51cb3aa927 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -713,7 +713,8 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { vtemp[m] = vfrac[m]; } Tequil = Ttemp + scale * dx[nmat]; - PORTABLE_REQUIRE(Tequil >= 0., "Negative temperature update in rho-T PTE solver iteration"); + PORTABLE_REQUIRE(Tequil >= 0., + "Negative temperature update in rho-T PTE solver iteration"); for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; rho[m] = robust::ratio(rhobar[m], vfrac[m]); @@ -1136,7 +1137,8 @@ class PTESolverFixedP : public mix_impl::PTESolverBase vtemp[m] = vfrac[m]; } Tequil = Ttemp + scale * dx[nmat]; - PORTABLE_REQUIRE(Tequil >= 0., "Negative temperature in Fixed P PTE solver iteration"); + PORTABLE_REQUIRE(Tequil >= 0., + "Negative temperature in Fixed P PTE solver iteration"); for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; rho[m] = robust::ratio(rhobar[m], vfrac[m]); From 3538048e4022a7d08088e295153c83bb87a50146 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 1 Jul 2024 21:04:09 -0600 Subject: [PATCH 12/61] Add comment about ScaleDx --- singularity-eos/closure/mixed_cell_models.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 51cb3aa927..4cae2f1c52 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -676,6 +676,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { PORTABLE_INLINE_FUNCTION Real ScaleDx() const { using namespace mix_params; + // Each check reduces the scale further if necessary Real scale = 1.0; // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { @@ -886,6 +887,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase PORTABLE_INLINE_FUNCTION Real ScaleDx() const { using namespace mix_params; + // Each check reduces the scale further if necessary Real scale = 1.0; // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { @@ -1100,6 +1102,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase PORTABLE_INLINE_FUNCTION Real ScaleDx() const { using namespace mix_params; + // Each check reduces the scale further if necessary Real scale = 1.0; // control how big of a step toward vfrac = 0 is allowed for (int m = 0; m < nmat; ++m) { @@ -1341,6 +1344,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { PORTABLE_INLINE_FUNCTION Real ScaleDx() const { using namespace mix_params; + // Each check reduces the scale further if necessary Real scale = 1.0; for (int m = 0; m < nmat; ++m) { // control how big of a step toward vfrac = 0 is allowed From 19162671c7a37b69ab382d0928abfada5bbede59 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 3 Jul 2024 18:05:17 -0600 Subject: [PATCH 13/61] Add mass fraction cutoff argument --- test/test_get_sg_eos.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_get_sg_eos.cpp b/test/test_get_sg_eos.cpp index 85c664ddf6..082c035011 100644 --- a/test/test_get_sg_eos.cpp +++ b/test/test_get_sg_eos.cpp @@ -58,7 +58,7 @@ int run_sg_get_eos_tests() { Real vfrac_true[NMAT], ie_true[NMAT]; get_sg_eos(NMAT, 1, 1, -1, eos_offset, eoss, &cell_offset, &P_true, &pmax, &v_true, &spvol, &sie_tot_true, &T_true_ev, &bmod, &dpde, &cv, mfrac, vfrac_true, - ie_true, nullptr, nullptr, nullptr); + ie_true, nullptr, nullptr, nullptr, 1.0e-12); Real sie_tot_check = 0.0; for (int m = 0; m < NMAT; ++m) { const Real r_m = mfrac[m] / vfrac_true[m]; @@ -85,7 +85,7 @@ int run_sg_get_eos_tests() { Real p_check, vfrac_check[NMAT], ie_check[NMAT]; get_sg_eos(NMAT, 1, 1, -3, eos_offset, eoss, &cell_offset, &p_check, &pmax, &v_true, &spvol, &sie_tot_check, &T_true_ev, &bmod, &dpde, &cv, mfrac, vfrac_check, - ie_check, nullptr, nullptr, nullptr); + ie_check, nullptr, nullptr, nullptr, 1.0e-12); // check output pressure and sie, indicate failure if relative err is too large if (std::abs(P_true - p_check) / std::abs(P_true) > 1.e-5 || std::abs(sie_tot_true - sie_tot_check) / std::abs(sie_tot_true) > 1.e-5) { @@ -109,7 +109,7 @@ int run_sg_get_eos_tests() { Real t_check; get_sg_eos(NMAT, 1, 1, -2, eos_offset, eoss, &cell_offset, &P_true, &pmax, &v_true, &spvol, &sie_tot_check, &t_check, &bmod, &dpde, &cv, mfrac, vfrac_check, - ie_check, nullptr, nullptr, nullptr); + ie_check, nullptr, nullptr, nullptr, 1.0e-12); // check output temperature and sie, indicate failure if relative err is too large if (std::abs(T_true_ev - t_check) / std::abs(T_true_ev) > 1.e-5 || std::abs(sie_tot_true - sie_tot_check) / std::abs(sie_tot_true) > 1.e-5) { From bb61c27bec8eef29b5082bce547303b528ddf9db Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 3 Jul 2024 18:05:43 -0600 Subject: [PATCH 14/61] Add closure unit test --- test/CMakeLists.txt | 17 ++++++ test/test_closure_pte.cpp | 125 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 142 insertions(+) create mode 100644 test/test_closure_pte.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3adfb58278..90e1e2b15d 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -45,6 +45,11 @@ add_executable( test_eos_stellar_collapse.cpp ) +add_executable( + closure_unit_tests + test_closure_pte.cpp + ) + get_property(plugin_tests GLOBAL PROPERTY PLUGIN_TESTS) if (plugin_tests) add_executable( @@ -62,11 +67,20 @@ endif() if(SINGULARITY_TEST_SESAME) target_compile_definitions(eos_tabulated_unit_tests PRIVATE SINGULARITY_TEST_SESAME) + target_compile_definitions(closure_unit_tests PRIVATE SINGULARITY_TEST_SESAME) endif() if(SINGULARITY_TEST_STELLAR_COLLAPSE) target_compile_definitions(eos_tabulated_unit_tests PRIVATE SINGULARITY_TEST_STELLAR_COLLAPSE) endif() +if(SINGULARITY_USE_SPINER) + target_compile_definitions(closure_unit_tests + PRIVATE SINGULARITY_USE_SPINER) +endif() +if(SINGULARITY_BUILD_CLOSURE) + target_compile_definitions(closure_unit_tests + PRIVATE SINGULARITY_BUILD_CLOSURE) +endif() target_link_libraries(eos_analytic_unit_tests PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) @@ -74,6 +88,8 @@ target_link_libraries(eos_infrastructure_tests PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) target_link_libraries(eos_tabulated_unit_tests PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) +target_link_libraries(closure_unit_tests PRIVATE Catch2::Catch2 + singularity-eos::singularity-eos) if (plugin_tests) target_link_libraries(eos_plugin_tests PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) @@ -82,6 +98,7 @@ include(Catch) catch_discover_tests(eos_analytic_unit_tests PROPERTIES TIMEOUT 60) catch_discover_tests(eos_infrastructure_tests PROPERTIES TIMEOUT 60) catch_discover_tests(eos_tabulated_unit_tests PROPERTIES TIMEOUT 60) +catch_discover_tests(closure_unit_tests PROPERTIES TIMEOUT 120) if (plugin_tests) catch_discover_tests(eos_plugin_tests PROPERTIES TIMEOUT 60) endif() diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp new file mode 100644 index 0000000000..f97646bd87 --- /dev/null +++ b/test/test_closure_pte.cpp @@ -0,0 +1,125 @@ +//------------------------------------------------------------------------------ +// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. + +#include +#include + +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +#include +#include +#include + +using namespace singularity; + +constexpr Real GPa = 1.0e10; +constexpr Real MJ_per_kg = 1.0e10; + +#ifdef SINGULARITY_BUILD_CLOSURE +#ifdef SINGULARITY_TEST_SESAME +#ifdef SINGULARITY_USE_SPINER +#ifdef SPINER_USE_HDF + +SCENARIO("Density-Temperature PTE Solver", "[PTE]") { + + GIVEN("Equations of state") { + // Set up the three EOS + constexpr int num_eos = 3; + std::array eos_arr; + // Reference state + constexpr Real P0 = 1.0e6; // 1 bar + constexpr Real T0 = 296; // K + // Ideal gas air + constexpr Real rho0_air = 1e-03; // g/cc + constexpr Real Gruneisen_air = 0.4; + constexpr Real CV_air = P0 / rho0_air / (Gruneisen_air * T0); + eos_arr[0] = IdealGas(Gruneisen_air, CV_air); + // Spiner copper EOS + constexpr int Cu_matid = 3337; + const std::string eos_file = "../materials.sp5"; + eos_arr[1] = SpinerEOSDependsRhoT(eos_file, Cu_matid); + // Davis Products EOS + constexpr Real a = 0.798311; + constexpr Real b = 0.58; + constexpr Real k = 1.35; + constexpr Real n = 2.66182; + constexpr Real vc = 0.75419; // cm^3 / g + constexpr Real pc = 3.2 * GPa; // microbar + constexpr Real Cv = 0.001072 * MJ_per_kg; // erg / g / K + constexpr Real Q = 4.115 * MJ_per_kg; + eos_arr[2] = DavisProducts(a, b, k, n, vc, pc, Cv).Modify(-Q); + + // Move EOS array from host to device + for (auto eos : eos_arr) { + eos.GetOnDevice(); + } + constexpr size_t EOS_bytes = num_eos * sizeof(EOS); + EOS* v_EOS = (EOS*)PORTABLE_MALLOC(EOS_bytes); + portableCopyToDevice(v_EOS, eos_arr.data(), EOS_bytes); + + GIVEN("A difficult state for the density-temperature PTE solver") { + // Define the state + constexpr Real spvol_bulk = 6.256037280402093e-01; + constexpr Real sie_bulk = -1.441692060406520e+10; + const std::array mass_frac = + {1.10382442033331e-10, 0.124935312146569, 0.875064687743048}; + // Initial guess: all materials at the cell density + const std::array vol_frac = mass_frac; + // Calculate material densities (and corresponding energies) and the total + // volume fraction + std::array densities; + std::array sies; + Real vfrac_sum = 0; + for (auto i = 0; i < num_eos; ++i) { + densities[i] = mass_frac[i] / spvol_bulk / vol_frac[i]; + sies[i] = sie_bulk * mass_frac[i]; + vfrac_sum += vol_frac[i]; + } + // Initialize pressure and temperature arrays to zero + std::array temperatures = {0., 0., 0.}; + std::array pressures = {0., 0., 0.}; + // Copy values to device (when available) + constexpr size_t bytes = num_eos * sizeof(Real); + Real* v_densities = (Real*)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_densities, densities.data(), bytes); + Real* v_vol_frac = (Real*)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_vol_frac, vol_frac.data(), bytes); + Real* v_sies = (Real*)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_sies, sies.data(), bytes); + Real* v_temperatures = (Real*)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_temperatures, temperatures.data(), bytes); + Real* v_pressures = (Real*)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_pressures, pressures.data(), bytes); + + THEN("The PTE solver will converge") { + // Allocate scratch space for the PTE solver + int pte_solver_scratch_size = PTESolverRhoTRequiredScratch(num_eos); + double* scratch = (double*)PORTABLE_MALLOC(pte_solver_scratch_size); + double** lambdas = (double**)PORTABLE_MALLOC(MAX_NUM_LAMBDAS * num_eos); + PTESolverRhoT method( + num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, + v_sies, v_temperatures, v_pressures, lambdas, + scratch); + bool pte_converged = PTESolver(method); + REQUIRE(pte_converged); + } + } + } +} +#endif // SPINER_USE_HDF +#endif // SINGULARITY_USE_SPINER +#endif // SINGULARITY_TEST_SESAME +#endif // SINGULARITY_BUILD_CLOSURE From 1e2952baa4d95bfa9285f860b10fd2d254130ce6 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 3 Jul 2024 18:08:44 -0600 Subject: [PATCH 15/61] Change TestUpdate() to take an argument to reset the cached variables and temporarily comment out vfrac checks --- singularity-eos/closure/mixed_cell_models.hpp | 31 ++++++++++++------- 1 file changed, 20 insertions(+), 11 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 4cae2f1c52..d6a40440f0 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -171,7 +171,7 @@ class PTESolverBase { virtual void Fixup() const { Real vsum = 0; for (int m = 0; m < nmat; ++m) { - PORTABLE_REQUIRE(vfrac[m] > 0, "Negative volume fraction in Fixup"); + // PORTABLE_REQUIRE(vfrac[m] > 0, "Negative volume fraction in Fixup"); vsum += vfrac[m]; } PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive in Fixup"); @@ -632,7 +632,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { const Real vf_pert = vfrac[m] + dv; const Real rho_pert = robust::ratio(rhobar[m], vf_pert); - PORTABLE_REQUIRE(vfrac[m] > 0, "Negative volume fraction in rho-T PTE iteration"); + // PORTABLE_REQUIRE(vfrac[m] > 0, "Negative volume fraction in rho-T PTE iteration"); Real p_pert{}; Real e_pert = @@ -707,8 +707,10 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale) { - if (scale == 1.0) { + Real TestUpdate(const Real scale, bool const cache_state=false) { + if (cache_state) { + // Store the current state in temp variables for first iteration of line + // search Ttemp = Tequil; for (int m = 0; m < nmat; ++m) vtemp[m] = vfrac[m]; @@ -912,8 +914,10 @@ class PTESolverFixedT : public mix_impl::PTESolverBase // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale) { - if (scale == 1.0) { + Real TestUpdate(const Real scale, bool const cache_state=false) { + if (cache_state) { + // Store the current state in temp variables for first iteration of line + // search for (int m = 0; m < nmat; ++m) vtemp[m] = vfrac[m]; } @@ -1133,8 +1137,10 @@ class PTESolverFixedP : public mix_impl::PTESolverBase // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale) { - if (scale == 1.0) { + Real TestUpdate(const Real scale, bool const cache_state=false) { + if (cache_state) { + // Store the current state in temp variables for first iteration of line + // search Ttemp = Tequil; for (int m = 0; m < nmat; ++m) vtemp[m] = vfrac[m]; @@ -1374,8 +1380,10 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale) const { - if (scale == 1.0) { + Real TestUpdate(const Real scale, bool const cache_state=false) const { + if (cache_state) { + // Store the current state in temp variables for first iteration of line + // search for (int m = 0; m < nmat; ++m) { vtemp[m] = vfrac[m]; utemp[m] = u[m]; @@ -1435,7 +1443,8 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { Real gradfdx = -2.0 * scale * err; scale = 1.0; // New scale for line search Real err_old = err; - err = s.TestUpdate(scale); + // Test the update and reset the cache the current state + err = s.TestUpdate(scale, true /* cache_state */); if (err > err_old + line_search_alpha * gradfdx) { // backtrack to middle of step scale = 0.5; From 52de0c47d60420c60a33ba768728556543a829df Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 3 Jul 2024 18:09:25 -0600 Subject: [PATCH 16/61] Move #endif to be in proper place --- singularity-eos/eos/get_sg_eos_functors.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/get_sg_eos_functors.hpp b/singularity-eos/eos/get_sg_eos_functors.hpp index d8633b2046..a6af7a3183 100644 --- a/singularity-eos/eos/get_sg_eos_functors.hpp +++ b/singularity-eos/eos/get_sg_eos_functors.hpp @@ -383,8 +383,8 @@ struct final_functor { PORTABLE_ALWAYS_ABORT( "Bad value RETURNED from singularity-eos. See output for details"); } - } #endif // #ifndef NDEBUG + } }; } // namespace singularity From ac21c030a16e5e3c2354159cded834a9e72266af Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 3 Jul 2024 18:23:06 -0600 Subject: [PATCH 17/61] Fix another misplaced #endif --- singularity-eos/eos/get_sg_eos_functors.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/get_sg_eos_functors.hpp b/singularity-eos/eos/get_sg_eos_functors.hpp index a6af7a3183..aee80f1933 100644 --- a/singularity-eos/eos/get_sg_eos_functors.hpp +++ b/singularity-eos/eos/get_sg_eos_functors.hpp @@ -151,8 +151,8 @@ struct init_functor { PORTABLE_ALWAYS_ABORT( "Bad values INPUT to singularity-eos interface. See output for details"); } - } #endif // #ifndef NDEBUG + } }; struct final_functor { From 6dd5411704271726bf4711bc6ae9aa0e92085f57 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 3 Jul 2024 18:25:55 -0600 Subject: [PATCH 18/61] Clang format --- singularity-eos/closure/mixed_cell_models.hpp | 11 +++--- test/test_closure_pte.cpp | 37 +++++++++---------- 2 files changed, 24 insertions(+), 24 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index d6a40440f0..57bffc7562 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -632,7 +632,8 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { const Real vf_pert = vfrac[m] + dv; const Real rho_pert = robust::ratio(rhobar[m], vf_pert); - // PORTABLE_REQUIRE(vfrac[m] > 0, "Negative volume fraction in rho-T PTE iteration"); + // PORTABLE_REQUIRE(vfrac[m] > 0, "Negative volume fraction in rho-T PTE + // iteration"); Real p_pert{}; Real e_pert = @@ -707,7 +708,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale, bool const cache_state=false) { + Real TestUpdate(const Real scale, bool const cache_state = false) { if (cache_state) { // Store the current state in temp variables for first iteration of line // search @@ -914,7 +915,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale, bool const cache_state=false) { + Real TestUpdate(const Real scale, bool const cache_state = false) { if (cache_state) { // Store the current state in temp variables for first iteration of line // search @@ -1137,7 +1138,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale, bool const cache_state=false) { + Real TestUpdate(const Real scale, bool const cache_state = false) { if (cache_state) { // Store the current state in temp variables for first iteration of line // search @@ -1380,7 +1381,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { // Update the solution and return new residual. Possibly called repeatedly with // different scale factors as part of a line search PORTABLE_INLINE_FUNCTION - Real TestUpdate(const Real scale, bool const cache_state=false) const { + Real TestUpdate(const Real scale, bool const cache_state = false) const { if (cache_state) { // Store the current state in temp variables for first iteration of line // search diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index f97646bd87..6aa941a7c8 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -20,8 +20,8 @@ #endif #include -#include #include +#include using namespace singularity; @@ -40,10 +40,10 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { constexpr int num_eos = 3; std::array eos_arr; // Reference state - constexpr Real P0 = 1.0e6; // 1 bar - constexpr Real T0 = 296; // K + constexpr Real P0 = 1.0e6; // 1 bar + constexpr Real T0 = 296; // K // Ideal gas air - constexpr Real rho0_air = 1e-03; // g/cc + constexpr Real rho0_air = 1e-03; // g/cc constexpr Real Gruneisen_air = 0.4; constexpr Real CV_air = P0 / rho0_air / (Gruneisen_air * T0); eos_arr[0] = IdealGas(Gruneisen_air, CV_air); @@ -56,8 +56,8 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { constexpr Real b = 0.58; constexpr Real k = 1.35; constexpr Real n = 2.66182; - constexpr Real vc = 0.75419; // cm^3 / g - constexpr Real pc = 3.2 * GPa; // microbar + constexpr Real vc = 0.75419; // cm^3 / g + constexpr Real pc = 3.2 * GPa; // microbar constexpr Real Cv = 0.001072 * MJ_per_kg; // erg / g / K constexpr Real Q = 4.115 * MJ_per_kg; eos_arr[2] = DavisProducts(a, b, k, n, vc, pc, Cv).Modify(-Q); @@ -67,15 +67,15 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { eos.GetOnDevice(); } constexpr size_t EOS_bytes = num_eos * sizeof(EOS); - EOS* v_EOS = (EOS*)PORTABLE_MALLOC(EOS_bytes); + EOS *v_EOS = (EOS *)PORTABLE_MALLOC(EOS_bytes); portableCopyToDevice(v_EOS, eos_arr.data(), EOS_bytes); GIVEN("A difficult state for the density-temperature PTE solver") { // Define the state constexpr Real spvol_bulk = 6.256037280402093e-01; constexpr Real sie_bulk = -1.441692060406520e+10; - const std::array mass_frac = - {1.10382442033331e-10, 0.124935312146569, 0.875064687743048}; + const std::array mass_frac = { + 1.10382442033331e-10, 0.124935312146569, 0.875064687743048}; // Initial guess: all materials at the cell density const std::array vol_frac = mass_frac; // Calculate material densities (and corresponding energies) and the total @@ -93,26 +93,25 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { std::array pressures = {0., 0., 0.}; // Copy values to device (when available) constexpr size_t bytes = num_eos * sizeof(Real); - Real* v_densities = (Real*)PORTABLE_MALLOC(bytes); + Real *v_densities = (Real *)PORTABLE_MALLOC(bytes); portableCopyToDevice(v_densities, densities.data(), bytes); - Real* v_vol_frac = (Real*)PORTABLE_MALLOC(bytes); + Real *v_vol_frac = (Real *)PORTABLE_MALLOC(bytes); portableCopyToDevice(v_vol_frac, vol_frac.data(), bytes); - Real* v_sies = (Real*)PORTABLE_MALLOC(bytes); + Real *v_sies = (Real *)PORTABLE_MALLOC(bytes); portableCopyToDevice(v_sies, sies.data(), bytes); - Real* v_temperatures = (Real*)PORTABLE_MALLOC(bytes); + Real *v_temperatures = (Real *)PORTABLE_MALLOC(bytes); portableCopyToDevice(v_temperatures, temperatures.data(), bytes); - Real* v_pressures = (Real*)PORTABLE_MALLOC(bytes); + Real *v_pressures = (Real *)PORTABLE_MALLOC(bytes); portableCopyToDevice(v_pressures, pressures.data(), bytes); THEN("The PTE solver will converge") { // Allocate scratch space for the PTE solver int pte_solver_scratch_size = PTESolverRhoTRequiredScratch(num_eos); - double* scratch = (double*)PORTABLE_MALLOC(pte_solver_scratch_size); - double** lambdas = (double**)PORTABLE_MALLOC(MAX_NUM_LAMBDAS * num_eos); + double *scratch = (double *)PORTABLE_MALLOC(pte_solver_scratch_size); + double **lambdas = (double **)PORTABLE_MALLOC(MAX_NUM_LAMBDAS * num_eos); PTESolverRhoT method( - num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, - v_sies, v_temperatures, v_pressures, lambdas, - scratch); + num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, + v_temperatures, v_pressures, lambdas, scratch); bool pte_converged = PTESolver(method); REQUIRE(pte_converged); } From bf284978710fc7c43ba82e7964080d8a8cb6d982 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 8 Jul 2024 18:36:22 -0600 Subject: [PATCH 19/61] Expliclty depend on ports-of-call1.6.0 for new printf --- spack-repo/packages/singularity-eos/package.py | 1 + 1 file changed, 1 insertion(+) diff --git a/spack-repo/packages/singularity-eos/package.py b/spack-repo/packages/singularity-eos/package.py index ef59aa9b76..8148ac4ba8 100644 --- a/spack-repo/packages/singularity-eos/package.py +++ b/spack-repo/packages/singularity-eos/package.py @@ -112,6 +112,7 @@ class SingularityEos(CMakePackage, CudaPackage): depends_on("ports-of-call@1.4.2,1.5.2:", when="@:1.7.0") depends_on("ports-of-call@1.5.2:", when="@1.7.1:") + depends_on("ports-of-call@1.6.0:", when="@1.7.2:") # request HEAD of main branch depends_on("ports-of-call@main", when="@main") From f017752371dba494d549c68845117d50f190bc65 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Jul 2024 18:03:24 -0600 Subject: [PATCH 20/61] Fix typo in scale update and skip update if appropriate --- singularity-eos/closure/mixed_cell_models.hpp | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 57bffc7562..6a761283b9 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -691,8 +691,15 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { const Real rho_min = std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); const Real alpha_max = robust::ratio(rhobar[m], rho_min); + if (alpha_max < vfrac[m]) { + // Despite our best efforts, we're already in the unstable regime (i.e. + // dPdV_T > 0) so we would actually want to *increase* the step instead + // of decreasing it. As a result, this code doesn't work as intended and + // should be skipped. + continue + } if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { - scale = robust::ratio(0.5 * alpha_max - vfrac[m], dx[m]); + scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]); } } // control how big of a step toward T = 0 is allowed From f913f97578d63877e695aa6ea2c785a9089cc84c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Jul 2024 18:03:53 -0600 Subject: [PATCH 21/61] Link Catch2 to closure --- test/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 90e1e2b15d..f67d60e0b3 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -47,6 +47,7 @@ add_executable( add_executable( closure_unit_tests + catch2_define.cpp test_closure_pte.cpp ) From d834c84e03ba13c132fa7f8a308f61731eb72d1a Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Jul 2024 18:04:39 -0600 Subject: [PATCH 22/61] Add PORTABLE_FREE and change check condition temporarily --- test/test_closure_pte.cpp | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 6aa941a7c8..47ff09e1fd 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -22,6 +22,7 @@ #include #include #include +#include using namespace singularity; @@ -113,9 +114,14 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, v_temperatures, v_pressures, lambdas, scratch); bool pte_converged = PTESolver(method); - REQUIRE(pte_converged); + CHECK(!pte_converged); + + // Free scratch and lambda memory + PORTABLE_FREE(scratch); + PORTABLE_FREE(lambdas); } } + PORTABLE_FREE(v_EOS); } } #endif // SPINER_USE_HDF From 724639e3517598e418aec81d84cc2ec43e97cd6b Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Jul 2024 18:08:50 -0600 Subject: [PATCH 23/61] Missing semicolon --- singularity-eos/closure/mixed_cell_models.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 6a761283b9..1beddf8614 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -696,7 +696,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { // dPdV_T > 0) so we would actually want to *increase* the step instead // of decreasing it. As a result, this code doesn't work as intended and // should be skipped. - continue + continue; } if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]); From 33269a40c2710c67c82450adc6225b96505112b7 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Jul 2024 18:55:48 -0600 Subject: [PATCH 24/61] Add back in a couple free statements --- test/test_closure_pte.cpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 47ff09e1fd..0306520302 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -77,8 +77,10 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { constexpr Real sie_bulk = -1.441692060406520e+10; const std::array mass_frac = { 1.10382442033331e-10, 0.124935312146569, 0.875064687743048}; + // Initial guess: all materials at the cell density const std::array vol_frac = mass_frac; + // Calculate material densities (and corresponding energies) and the total // volume fraction std::array densities; @@ -89,9 +91,11 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { sies[i] = sie_bulk * mass_frac[i]; vfrac_sum += vol_frac[i]; } + // Initialize pressure and temperature arrays to zero std::array temperatures = {0., 0., 0.}; std::array pressures = {0., 0., 0.}; + // Copy values to device (when available) constexpr size_t bytes = num_eos * sizeof(Real); Real *v_densities = (Real *)PORTABLE_MALLOC(bytes); @@ -105,7 +109,7 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { Real *v_pressures = (Real *)PORTABLE_MALLOC(bytes); portableCopyToDevice(v_pressures, pressures.data(), bytes); - THEN("The PTE solver will converge") { + THEN("The PTE solver should converge") { // Allocate scratch space for the PTE solver int pte_solver_scratch_size = PTESolverRhoTRequiredScratch(num_eos); double *scratch = (double *)PORTABLE_MALLOC(pte_solver_scratch_size); @@ -113,14 +117,14 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { PTESolverRhoT method( num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, v_temperatures, v_pressures, lambdas, scratch); - bool pte_converged = PTESolver(method); - CHECK(!pte_converged); - // Free scratch and lambda memory - PORTABLE_FREE(scratch); + // Solve the PTE system and ensure it converged + bool pte_converged = PTESolver(method); + CHECK(pte_converged); PORTABLE_FREE(lambdas); } } + // Free EOS memory PORTABLE_FREE(v_EOS); } } From ae40f6a4237221b441a33c03a07e5dc77afa4781 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Jul 2024 19:04:11 -0600 Subject: [PATCH 25/61] Fix bug in other PTE solvers --- singularity-eos/closure/mixed_cell_models.hpp | 21 +++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 1beddf8614..55d0ce0f38 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -909,6 +909,13 @@ class PTESolverFixedT : public mix_impl::PTESolverBase for (int m = 0; m < nmat; m++) { const Real rho_min = eos[m].RhoPmin(Tequil); const Real alpha_max = robust::ratio(rhobar[m], rho_min); + if (alpha_max < vfrac[m]) { + // Despite our best efforts, we're already in the unstable regime (i.e. + // dPdV_T > 0) so we would actually want to *increase* the step instead + // of decreasing it. As a result, this code doesn't work as intended and + // should be skipped. + continue; + } if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { scale = robust::ratio(0.5 * (alpha_max - vfrac[m]), dx[m]); } @@ -1128,6 +1135,13 @@ class PTESolverFixedP : public mix_impl::PTESolverBase const Real rho_min = std::max(eos[m].RhoPmin(Tnorm * Tequil), eos[m].RhoPmin(Tnorm * Tnew)); const Real alpha_max = robust::ratio(rhobar[m], rho_min); + if (alpha_max < vfrac[m]) { + // Despite our best efforts, we're already in the unstable regime (i.e. + // dPdV_T > 0) so we would actually want to *increase* the step instead + // of decreasing it. As a result, this code doesn't work as intended and + // should be skipped. + continue; + } if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); } @@ -1375,6 +1389,13 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { std::max(eos[m].RhoPmin(Tnorm * temp[m]), eos[m].RhoPmin(Tnorm * tt)); const Real alpha_max = robust::ratio(rhobar[m], rho_min); // control how big of a step toward rho = rho(Pmin) is allowed + if (alpha_max < vfrac[m]) { + // Despite our best efforts, we're already in the unstable regime (i.e. + // dPdV_T > 0) so we would actually want to *increase* the step instead + // of decreasing it. As a result, this code doesn't work as intended and + // should be skipped. + continue; + } if (scale * dx[m] > 0.5 * (alpha_max - vfrac[m])) { scale = 0.5 * robust::ratio(alpha_max - vfrac[m], dx[m]); } From b994fdf5a59674d24a01e542ac57600c97d26a8f Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Jul 2024 19:06:41 -0600 Subject: [PATCH 26/61] Remove requirements on volume fractions within iteration since there's better information afterwards --- singularity-eos/closure/mixed_cell_models.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 55d0ce0f38..29707edeab 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -171,7 +171,6 @@ class PTESolverBase { virtual void Fixup() const { Real vsum = 0; for (int m = 0; m < nmat; ++m) { - // PORTABLE_REQUIRE(vfrac[m] > 0, "Negative volume fraction in Fixup"); vsum += vfrac[m]; } PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive in Fixup"); @@ -632,9 +631,6 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { const Real vf_pert = vfrac[m] + dv; const Real rho_pert = robust::ratio(rhobar[m], vf_pert); - // PORTABLE_REQUIRE(vfrac[m] > 0, "Negative volume fraction in rho-T PTE - // iteration"); - Real p_pert{}; Real e_pert = eos[m].InternalEnergyFromDensityTemperature(rho_pert, Tnorm * Tequil, Cache[m]); From af9e22a0bf5b8c4d4206dd8e4a33eee2eb31b633 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Jul 2024 19:13:42 -0600 Subject: [PATCH 27/61] Add check to ensure step size is always less than or equal to Newton step --- singularity-eos/closure/mixed_cell_models.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 29707edeab..649f8cbbdc 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -1463,6 +1463,7 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { // possibly scale the update to stay within reasonable bounds Real scale = s.ScaleDx(); + PORTABLE_REQUIRE(scale <= 1.0, "PTE Solver is attempting to increase the step size"); // Line search Real gradfdx = -2.0 * scale * err; From ddb1969a3a6f754444430ff06d9ac7c160560442 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Jul 2024 19:17:47 -0600 Subject: [PATCH 28/61] Add PTE solver bug --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8f4dfaf436..cd69599b40 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -26,6 +26,7 @@ - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface - [[PR373]](https://github.com/lanl/singularity-eos/pull/373) Initialize cache in `get_sg_eos*` functions - [[PR374]](https://github.com/lanl/singularity-eos/pull/374) Make the Davis EOS more numerically robust +- [[PR383]](https://github.com/lanl/singularity-eos/pull/383) Fix bug in step scaling for PTE solver ### Changed (changing behavior/API/variables/...) - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls From 1c5d5b72678b9a5da6357e4bba1a8cc8d230a69b Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 10 Jul 2024 19:24:40 -0600 Subject: [PATCH 29/61] Change singularity version for updated POC --- spack-repo/packages/singularity-eos/package.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/spack-repo/packages/singularity-eos/package.py b/spack-repo/packages/singularity-eos/package.py index 8148ac4ba8..916294a1d1 100644 --- a/spack-repo/packages/singularity-eos/package.py +++ b/spack-repo/packages/singularity-eos/package.py @@ -112,7 +112,7 @@ class SingularityEos(CMakePackage, CudaPackage): depends_on("ports-of-call@1.4.2,1.5.2:", when="@:1.7.0") depends_on("ports-of-call@1.5.2:", when="@1.7.1:") - depends_on("ports-of-call@1.6.0:", when="@1.7.2:") + depends_on("ports-of-call@1.6.0:", when="@1.8.1:") # request HEAD of main branch depends_on("ports-of-call@main", when="@main") From 65788e5933eb93508839ef6a7fd7a249cef164db Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 11 Jul 2024 13:46:11 -0600 Subject: [PATCH 30/61] Remove debug checks from inside the PTE solver iteration --- singularity-eos/closure/mixed_cell_models.hpp | 6 ------ 1 file changed, 6 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 649f8cbbdc..88af8d5b55 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -173,7 +173,6 @@ class PTESolverBase { for (int m = 0; m < nmat; ++m) { vsum += vfrac[m]; } - PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive in Fixup"); for (int m = 0; m < nmat; ++m) vfrac[m] *= robust::ratio(vfrac_total, vsum); } @@ -183,7 +182,6 @@ class PTESolverBase { void Finalize() { for (int m = 0; m < nmat; m++) { temp[m] *= Tnorm; - PORTABLE_REQUIRE(temp[m] >= 0., "Non-positive temperature returned in Finalize"); u[m] *= uscale; press[m] *= uscale; } @@ -720,8 +718,6 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { vtemp[m] = vfrac[m]; } Tequil = Ttemp + scale * dx[nmat]; - PORTABLE_REQUIRE(Tequil >= 0., - "Negative temperature update in rho-T PTE solver iteration"); for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; rho[m] = robust::ratio(rhobar[m], vfrac[m]); @@ -1164,8 +1160,6 @@ class PTESolverFixedP : public mix_impl::PTESolverBase vtemp[m] = vfrac[m]; } Tequil = Ttemp + scale * dx[nmat]; - PORTABLE_REQUIRE(Tequil >= 0., - "Negative temperature in Fixed P PTE solver iteration"); for (int m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; rho[m] = robust::ratio(rhobar[m], vfrac[m]); From 49fbd27d14ae837899206300920a87312e8ad863 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 11 Jul 2024 14:53:27 -0600 Subject: [PATCH 31/61] Add PORTABLE_FREE and finalize... this should work but doesn't --- test/test_closure_pte.cpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 0306520302..ad3ca66585 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -121,9 +121,16 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { // Solve the PTE system and ensure it converged bool pte_converged = PTESolver(method); CHECK(pte_converged); + + // Free temp memory PORTABLE_FREE(lambdas); + PORTABLE_FREE(scratch); } } + // Call Finalize on each EOS + for (auto eos : eos_arr) { + eos.Finalize(); + } // Free EOS memory PORTABLE_FREE(v_EOS); } From 0974b1142814bbce39ff4b63b70dab5d792b70da Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 11 Jul 2024 17:06:05 -0600 Subject: [PATCH 32/61] Fix memory allocation --- test/test_closure_pte.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index ad3ca66585..8681061f3d 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -111,9 +111,11 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { THEN("The PTE solver should converge") { // Allocate scratch space for the PTE solver - int pte_solver_scratch_size = PTESolverRhoTRequiredScratch(num_eos); - double *scratch = (double *)PORTABLE_MALLOC(pte_solver_scratch_size); - double **lambdas = (double **)PORTABLE_MALLOC(MAX_NUM_LAMBDAS * num_eos); + const int pte_solver_scratch_size = PTESolverRhoTRequiredScratch(num_eos); + const size_t scratch_bytes = pte_solver_scratch_size * sizeof(Real); + double *scratch = (double *)PORTABLE_MALLOC(scratch_bytes); + constexpr size_t lambda_bytes = MAX_NUM_LAMBDAS * num_eos * sizeof(Real); + double **lambdas = (double **)PORTABLE_MALLOC(lambda_bytes); PTESolverRhoT method( num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, v_temperatures, v_pressures, lambdas, scratch); From 1a048ad8dc3dcbd461b50d8f36debf7c46ef0952 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 15 Jul 2024 18:34:53 -0600 Subject: [PATCH 33/61] Break out parts of tests into functions --- test/test_closure_pte.cpp | 199 ++++++++++++++++++++++++++------------ 1 file changed, 136 insertions(+), 63 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 8681061f3d..d300f1210d 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -11,8 +11,10 @@ // prepare derivative works, distribute copies to the public, perform // publicly and display publicly, and to permit others to do so. -#include +#ifdef SINGULARITY_BUILD_CLOSURE + #include +#include #ifndef CATCH_CONFIG_FAST_COMPILE #define CATCH_CONFIG_FAST_COMPILE @@ -29,17 +31,111 @@ using namespace singularity; constexpr Real GPa = 1.0e10; constexpr Real MJ_per_kg = 1.0e10; -#ifdef SINGULARITY_BUILD_CLOSURE #ifdef SINGULARITY_TEST_SESAME #ifdef SINGULARITY_USE_SPINER #ifdef SPINER_USE_HDF +using myEOS = singularity::Variant, DavisProducts, + SpinerEOSDependsRhoT>; + +template +myEOS *copy_eos_arr_to_device(const int num_eos, EOSArrT eos_arr) { + // Move EOS array from host to device + for (auto i = 0; i < num_eos; i++) { + eos_arr[i] = eos_arr[i].GetOnDevice(); + } + const size_t EOS_bytes = num_eos * sizeof(myEOS); + myEOS *v_EOS = (myEOS *)PORTABLE_MALLOC(EOS_bytes); + portableCopyToDevice(v_EOS, eos_arr.data(), EOS_bytes); + + return v_EOS; +} + +void finalize_eos_arr(const int num_eos, myEOS *v_EOS) { + // Call Finalize on each EOS + for (auto i = 0; i < num_eos; i++) { + v_EOS[i].Finalize(); + } + // Free EOS memory + PORTABLE_FREE(v_EOS); +} + +template +bool run_PTE_from_state(const int num_eos, myEOS *v_EOS, const Real spvol_bulk, + const Real sie_bulk, ArrT mass_frac) { + // Calculate material densities (and corresponding energies) and the total + // volume fraction + std::vector vol_frac(num_eos); + std::vector densities(num_eos); + std::vector sies(num_eos); + std::vector temperatures(num_eos); + std::vector pressures(num_eos); + Real vfrac_sum = 0; + for (auto i = 0; i < num_eos; ++i) { + // Initial guess: all materials at the cell density + vol_frac[i] = mass_frac[i]; + densities[i] = mass_frac[i] / spvol_bulk / vol_frac[i]; + sies[i] = sie_bulk * mass_frac[i]; + vfrac_sum += vol_frac[i]; + // Initialize pressure and temperature arrays to zero + temperatures[i] = 0; + pressures[i] = 0; + } + + // Copy values to device (when available) + const size_t bytes = num_eos * sizeof(Real); + Real *v_densities = (Real *)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_densities, densities.data(), bytes); + Real *v_vol_frac = (Real *)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_vol_frac, vol_frac.data(), bytes); + Real *v_sies = (Real *)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_sies, sies.data(), bytes); + Real *v_temperatures = (Real *)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_temperatures, temperatures.data(), bytes); + Real *v_pressures = (Real *)PORTABLE_MALLOC(bytes); + portableCopyToDevice(v_pressures, pressures.data(), bytes); + + // Allocate scratch space for the PTE solver + const int pte_solver_scratch_size = PTESolverRhoTRequiredScratch(num_eos); + const size_t scratch_bytes = pte_solver_scratch_size * sizeof(Real); + Real *scratch = (double *)PORTABLE_MALLOC(scratch_bytes); + + // Allocate lambdas for each EOS + constexpr size_t lambda_bytes = MAX_NUM_LAMBDAS * sizeof(Real); // EOS lambda size + const size_t lambda_size = num_eos * sizeof(Real *); // Size of lambda array + Real **lambdas = (Real **)PORTABLE_MALLOC(lambda_size); + for (auto i = 0; i < num_eos; i++) { + lambdas[i] = (Real *)PORTABLE_MALLOC(lambda_bytes); + } + + // Set up the PTE solver + PTESolverRhoT method( + num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, + v_temperatures, v_pressures, lambdas, scratch); + + // Solve the PTE system + bool pte_converged = PTESolver(method); + + // Free temp memory + PORTABLE_FREE(lambdas); + PORTABLE_FREE(scratch); + + // Free PTE values + PORTABLE_FREE(v_densities); + PORTABLE_FREE(v_vol_frac); + PORTABLE_FREE(v_sies); + PORTABLE_FREE(v_temperatures); + PORTABLE_FREE(v_pressures); + + return pte_converged; +} + SCENARIO("Density-Temperature PTE Solver", "[PTE]") { - GIVEN("Equations of state") { + GIVEN("Three equations of state") { // Set up the three EOS constexpr int num_eos = 3; - std::array eos_arr; + std::array eos_arr; // Reference state constexpr Real P0 = 1.0e6; // 1 bar constexpr Real T0 = 296; // K @@ -63,14 +159,6 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { constexpr Real Q = 4.115 * MJ_per_kg; eos_arr[2] = DavisProducts(a, b, k, n, vc, pc, Cv).Modify(-Q); - // Move EOS array from host to device - for (auto eos : eos_arr) { - eos.GetOnDevice(); - } - constexpr size_t EOS_bytes = num_eos * sizeof(EOS); - EOS *v_EOS = (EOS *)PORTABLE_MALLOC(EOS_bytes); - portableCopyToDevice(v_EOS, eos_arr.data(), EOS_bytes); - GIVEN("A difficult state for the density-temperature PTE solver") { // Define the state constexpr Real spvol_bulk = 6.256037280402093e-01; @@ -78,63 +166,48 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { const std::array mass_frac = { 1.10382442033331e-10, 0.124935312146569, 0.875064687743048}; - // Initial guess: all materials at the cell density - const std::array vol_frac = mass_frac; - - // Calculate material densities (and corresponding energies) and the total - // volume fraction - std::array densities; - std::array sies; - Real vfrac_sum = 0; - for (auto i = 0; i < num_eos; ++i) { - densities[i] = mass_frac[i] / spvol_bulk / vol_frac[i]; - sies[i] = sie_bulk * mass_frac[i]; - vfrac_sum += vol_frac[i]; + THEN("The PTE solver should converge") { + myEOS *v_EOS = copy_eos_arr_to_device(num_eos, eos_arr); + const bool pte_converged = + run_PTE_from_state(num_eos, v_EOS, spvol_bulk, sie_bulk, mass_frac); + CHECK(pte_converged); + finalize_eos_arr(num_eos, v_EOS); } + } + } + GIVEN("Two equations of state") { + // Set up the three EOS + constexpr int num_eos = 2; + std::array eos_arr; + // Reference state + constexpr Real P0 = 1.0e6; // 1 bar + constexpr Real T0 = 296; // K + // Ideal gas air + constexpr Real rho0_air = 1e-03; // g/cc + constexpr Real Gruneisen_air = 0.4; + constexpr Real CV_air = P0 / rho0_air / (Gruneisen_air * T0); + eos_arr[0] = IdealGas(Gruneisen_air, CV_air); + // Spiner copper EOS + constexpr int Cu_matid = 3337; + const std::string eos_file = "../materials.sp5"; + eos_arr[1] = SpinerEOSDependsRhoT(eos_file, Cu_matid); - // Initialize pressure and temperature arrays to zero - std::array temperatures = {0., 0., 0.}; - std::array pressures = {0., 0., 0.}; - - // Copy values to device (when available) - constexpr size_t bytes = num_eos * sizeof(Real); - Real *v_densities = (Real *)PORTABLE_MALLOC(bytes); - portableCopyToDevice(v_densities, densities.data(), bytes); - Real *v_vol_frac = (Real *)PORTABLE_MALLOC(bytes); - portableCopyToDevice(v_vol_frac, vol_frac.data(), bytes); - Real *v_sies = (Real *)PORTABLE_MALLOC(bytes); - portableCopyToDevice(v_sies, sies.data(), bytes); - Real *v_temperatures = (Real *)PORTABLE_MALLOC(bytes); - portableCopyToDevice(v_temperatures, temperatures.data(), bytes); - Real *v_pressures = (Real *)PORTABLE_MALLOC(bytes); - portableCopyToDevice(v_pressures, pressures.data(), bytes); + GIVEN("A state that would cause a negative temperature in the PTE solver") { + // Define the state + constexpr Real spvol_bulk = 4.010467628234189e-01; + constexpr Real sie_bulk = 3.290180957185173e+07; + const std::array mass_frac = {0.000312273191678158, + 0.999687726808322}; THEN("The PTE solver should converge") { - // Allocate scratch space for the PTE solver - const int pte_solver_scratch_size = PTESolverRhoTRequiredScratch(num_eos); - const size_t scratch_bytes = pte_solver_scratch_size * sizeof(Real); - double *scratch = (double *)PORTABLE_MALLOC(scratch_bytes); - constexpr size_t lambda_bytes = MAX_NUM_LAMBDAS * num_eos * sizeof(Real); - double **lambdas = (double **)PORTABLE_MALLOC(lambda_bytes); - PTESolverRhoT method( - num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, - v_temperatures, v_pressures, lambdas, scratch); - - // Solve the PTE system and ensure it converged - bool pte_converged = PTESolver(method); - CHECK(pte_converged); - - // Free temp memory - PORTABLE_FREE(lambdas); - PORTABLE_FREE(scratch); + myEOS *v_EOS = copy_eos_arr_to_device(num_eos, eos_arr); + const bool pte_converged = + run_PTE_from_state(num_eos, v_EOS, spvol_bulk, sie_bulk, mass_frac); + // TODO: make this test converge + // CHECK(pte_converged); + finalize_eos_arr(num_eos, v_EOS); } } - // Call Finalize on each EOS - for (auto eos : eos_arr) { - eos.Finalize(); - } - // Free EOS memory - PORTABLE_FREE(v_EOS); } } #endif // SPINER_USE_HDF From b03cb4414779cce6dc1f025b8b53ff8bb3a82490 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 15 Jul 2024 18:40:02 -0600 Subject: [PATCH 34/61] Remove portable copy and instead just output GetOnDevice() to device array --- test/test_closure_pte.cpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index d300f1210d..c5da71aaf1 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -41,12 +41,11 @@ using myEOS = singularity::Variant, DavisPro template myEOS *copy_eos_arr_to_device(const int num_eos, EOSArrT eos_arr) { // Move EOS array from host to device - for (auto i = 0; i < num_eos; i++) { - eos_arr[i] = eos_arr[i].GetOnDevice(); - } const size_t EOS_bytes = num_eos * sizeof(myEOS); myEOS *v_EOS = (myEOS *)PORTABLE_MALLOC(EOS_bytes); - portableCopyToDevice(v_EOS, eos_arr.data(), EOS_bytes); + for (auto i = 0; i < num_eos; i++) { + v_EOS[i] = eos_arr[i].GetOnDevice(); + } return v_EOS; } From 843057c5bc1456306c36abf944ef1842691c9540 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 15 Jul 2024 18:45:00 -0600 Subject: [PATCH 35/61] add missing header --- test/test_closure_pte.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index c5da71aaf1..d9afb7327c 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -13,6 +13,7 @@ #ifdef SINGULARITY_BUILD_CLOSURE +#include #include #include From ae688b8f34f85fb98e45bc6615354a41743bf651 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 15 Jul 2024 18:47:45 -0600 Subject: [PATCH 36/61] Consolidate ifdef statements --- test/test_closure_pte.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index d9afb7327c..4524086039 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -33,8 +33,7 @@ constexpr Real GPa = 1.0e10; constexpr Real MJ_per_kg = 1.0e10; #ifdef SINGULARITY_TEST_SESAME -#ifdef SINGULARITY_USE_SPINER -#ifdef SPINER_USE_HDF +#ifdef SINGULARITY_USE_SPINER_WITH_HDF5 using myEOS = singularity::Variant, DavisProducts, SpinerEOSDependsRhoT>; @@ -210,7 +209,6 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { } } } -#endif // SPINER_USE_HDF -#endif // SINGULARITY_USE_SPINER +#endif // SINGULARITY_USE_SPINER_WITH_HDF5 #endif // SINGULARITY_TEST_SESAME #endif // SINGULARITY_BUILD_CLOSURE From 50edd2290fc20b2b0f4a6d15b97ee8d8506ef3eb Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 16 Jul 2024 12:10:34 -0600 Subject: [PATCH 37/61] Get on device and _then_ copy EOS array to device. Also free each lambda separately --- test/test_closure_pte.cpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 4524086039..b1c5034bca 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -42,11 +42,12 @@ template myEOS *copy_eos_arr_to_device(const int num_eos, EOSArrT eos_arr) { // Move EOS array from host to device const size_t EOS_bytes = num_eos * sizeof(myEOS); - myEOS *v_EOS = (myEOS *)PORTABLE_MALLOC(EOS_bytes); for (auto i = 0; i < num_eos; i++) { - v_EOS[i] = eos_arr[i].GetOnDevice(); + eos_arr[i] = eos_arr[i].GetOnDevice(); } - + myEOS *v_EOS = (myEOS *)PORTABLE_MALLOC(EOS_bytes); + const size_t bytes = num_eos * sizeof(myEOS); + portableCopyToDevice(v_EOS, eos_arr.data(), bytes); return v_EOS; } @@ -116,7 +117,11 @@ bool run_PTE_from_state(const int num_eos, myEOS *v_EOS, const Real spvol_bulk, bool pte_converged = PTESolver(method); // Free temp memory - PORTABLE_FREE(lambdas); + for (auto i = 0; i < num_eos; i++) { + // Free each lambda separately first + PORTABLE_FREE(lambdas[i]); + } + PORTABLE_FREE(lambdas); // Free entire lambda array PORTABLE_FREE(scratch); // Free PTE values From 10fb97a4ce0c854988afe492a26f9bfb37a60cb5 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 16 Jul 2024 12:23:19 -0600 Subject: [PATCH 38/61] Remove using namespace singularity and instead explicitly bring things in one by one --- test/test_closure_pte.cpp | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index b1c5034bca..b387cd15cb 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -27,31 +27,37 @@ #include #include -using namespace singularity; - constexpr Real GPa = 1.0e10; constexpr Real MJ_per_kg = 1.0e10; #ifdef SINGULARITY_TEST_SESAME #ifdef SINGULARITY_USE_SPINER_WITH_HDF5 -using myEOS = singularity::Variant, DavisProducts, +using singularity::IdealGas; +using singularity::ShiftedEOS; +using singularity::DavisProducts; +using singularity::SpinerEOSDependsRhoT; +using singularity::PTESolverRhoT; +using singularity::PTESolverRhoTRequiredScratch; +using singularity::MAX_NUM_LAMBDAS; + +using EOS = singularity::Variant, DavisProducts, SpinerEOSDependsRhoT>; template -myEOS *copy_eos_arr_to_device(const int num_eos, EOSArrT eos_arr) { +EOS *copy_eos_arr_to_device(const int num_eos, EOSArrT eos_arr) { // Move EOS array from host to device - const size_t EOS_bytes = num_eos * sizeof(myEOS); + const size_t EOS_bytes = num_eos * sizeof(EOS); for (auto i = 0; i < num_eos; i++) { eos_arr[i] = eos_arr[i].GetOnDevice(); } - myEOS *v_EOS = (myEOS *)PORTABLE_MALLOC(EOS_bytes); - const size_t bytes = num_eos * sizeof(myEOS); + EOS *v_EOS = (EOS *)PORTABLE_MALLOC(EOS_bytes); + const size_t bytes = num_eos * sizeof(EOS); portableCopyToDevice(v_EOS, eos_arr.data(), bytes); return v_EOS; } -void finalize_eos_arr(const int num_eos, myEOS *v_EOS) { +void finalize_eos_arr(const int num_eos, EOS *v_EOS) { // Call Finalize on each EOS for (auto i = 0; i < num_eos; i++) { v_EOS[i].Finalize(); @@ -61,7 +67,7 @@ void finalize_eos_arr(const int num_eos, myEOS *v_EOS) { } template -bool run_PTE_from_state(const int num_eos, myEOS *v_EOS, const Real spvol_bulk, +bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, const Real sie_bulk, ArrT mass_frac) { // Calculate material densities (and corresponding energies) and the total // volume fraction @@ -139,7 +145,7 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { GIVEN("Three equations of state") { // Set up the three EOS constexpr int num_eos = 3; - std::array eos_arr; + std::array eos_arr; // Reference state constexpr Real P0 = 1.0e6; // 1 bar constexpr Real T0 = 296; // K @@ -171,7 +177,7 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { 1.10382442033331e-10, 0.124935312146569, 0.875064687743048}; THEN("The PTE solver should converge") { - myEOS *v_EOS = copy_eos_arr_to_device(num_eos, eos_arr); + EOS *v_EOS = copy_eos_arr_to_device(num_eos, eos_arr); const bool pte_converged = run_PTE_from_state(num_eos, v_EOS, spvol_bulk, sie_bulk, mass_frac); CHECK(pte_converged); @@ -182,7 +188,7 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { GIVEN("Two equations of state") { // Set up the three EOS constexpr int num_eos = 2; - std::array eos_arr; + std::array eos_arr; // Reference state constexpr Real P0 = 1.0e6; // 1 bar constexpr Real T0 = 296; // K @@ -204,7 +210,7 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { 0.999687726808322}; THEN("The PTE solver should converge") { - myEOS *v_EOS = copy_eos_arr_to_device(num_eos, eos_arr); + EOS *v_EOS = copy_eos_arr_to_device(num_eos, eos_arr); const bool pte_converged = run_PTE_from_state(num_eos, v_EOS, spvol_bulk, sie_bulk, mass_frac); // TODO: make this test converge From a9edf92349e885ed6087f96688d7c68f07b47fb8 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 16 Jul 2024 12:23:49 -0600 Subject: [PATCH 39/61] Clang format --- test/test_closure_pte.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index b387cd15cb..b3f81b4493 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -33,16 +33,16 @@ constexpr Real MJ_per_kg = 1.0e10; #ifdef SINGULARITY_TEST_SESAME #ifdef SINGULARITY_USE_SPINER_WITH_HDF5 -using singularity::IdealGas; -using singularity::ShiftedEOS; using singularity::DavisProducts; -using singularity::SpinerEOSDependsRhoT; +using singularity::IdealGas; +using singularity::MAX_NUM_LAMBDAS; using singularity::PTESolverRhoT; using singularity::PTESolverRhoTRequiredScratch; -using singularity::MAX_NUM_LAMBDAS; +using singularity::ShiftedEOS; +using singularity::SpinerEOSDependsRhoT; using EOS = singularity::Variant, DavisProducts, - SpinerEOSDependsRhoT>; + SpinerEOSDependsRhoT>; template EOS *copy_eos_arr_to_device(const int num_eos, EOSArrT eos_arr) { From 8e89d457766c0f26ae4fb1d892f4c28f2ede1248 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 16 Jul 2024 12:53:11 -0600 Subject: [PATCH 40/61] Make PTE call on device --- test/test_closure_pte.cpp | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index b3f81b4493..f99ed43d56 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -114,13 +114,17 @@ bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, lambdas[i] = (Real *)PORTABLE_MALLOC(lambda_bytes); } - // Set up the PTE solver - PTESolverRhoT method( - num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, - v_temperatures, v_pressures, lambdas, scratch); - - // Solve the PTE system - bool pte_converged = PTESolver(method); + // Solve the PTE system on device using a one-teration portableFor + bool pte_converged; + constexpr size_t bool_bytes = 1 * sizeof(bool); + bool *pte_converged_d = (bool *)PORTABLE_MALLOC(bool_bytes); + portableFor("Device execution of PTE Test", 0, 1, PORTABLE_LAMBDA(auto i) { + PTESolverRhoT method( + num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, + v_temperatures, v_pressures, lambdas, scratch); + pte_converged_d[0] = PTESolver(method); + }); + portableCopyToHost(&pte_converged, pte_converged_d, bool_bytes); // Free temp memory for (auto i = 0; i < num_eos; i++) { From a27eecb8a328e2bf5d7169ace3dd0869a38ad45d Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 16 Jul 2024 12:53:36 -0600 Subject: [PATCH 41/61] Clang format --- test/test_closure_pte.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index f99ed43d56..864abd0683 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -118,12 +118,13 @@ bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, bool pte_converged; constexpr size_t bool_bytes = 1 * sizeof(bool); bool *pte_converged_d = (bool *)PORTABLE_MALLOC(bool_bytes); - portableFor("Device execution of PTE Test", 0, 1, PORTABLE_LAMBDA(auto i) { - PTESolverRhoT method( - num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, - v_temperatures, v_pressures, lambdas, scratch); - pte_converged_d[0] = PTESolver(method); - }); + portableFor( + "Device execution of PTE Test", 0, 1, PORTABLE_LAMBDA(auto i) { + PTESolverRhoT method( + num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, + v_temperatures, v_pressures, lambdas, scratch); + pte_converged_d[0] = PTESolver(method); + }); portableCopyToHost(&pte_converged, pte_converged_d, bool_bytes); // Free temp memory From a825f1e73317f996c8196547ba832cf2011fd04b Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 16 Jul 2024 17:45:29 -0600 Subject: [PATCH 42/61] Change type of iteration variable --- test/test_closure_pte.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 864abd0683..bfd8b3fbe0 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -119,7 +119,7 @@ bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, constexpr size_t bool_bytes = 1 * sizeof(bool); bool *pte_converged_d = (bool *)PORTABLE_MALLOC(bool_bytes); portableFor( - "Device execution of PTE Test", 0, 1, PORTABLE_LAMBDA(auto i) { + "Device execution of PTE Test", 0, 1, PORTABLE_LAMBDA(int i) { PTESolverRhoT method( num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, v_temperatures, v_pressures, lambdas, scratch); From 5e89dd3aee2ffcd9ad0a3f7f2a363bb1d0f00d6a Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Tue, 16 Jul 2024 19:15:57 -0600 Subject: [PATCH 43/61] Add copper material for closure unit test --- sesame2spiner/examples/unit_tests/copper.dat | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) create mode 100644 sesame2spiner/examples/unit_tests/copper.dat diff --git a/sesame2spiner/examples/unit_tests/copper.dat b/sesame2spiner/examples/unit_tests/copper.dat new file mode 100644 index 0000000000..3ccdcdc544 --- /dev/null +++ b/sesame2spiner/examples/unit_tests/copper.dat @@ -0,0 +1,20 @@ +# ====================================================================== +# Example input deck for sesame2spiner, +# a tool for converting eospac to spiner +# Author: Jonah Miller (jonahm@lanl.gov) +# © 2021-2023. Triad National Security, LLC. All rights reserved. This +# program was produced under U.S. Government contract 89233218CNA000001 +# for Los Alamos National Laboratory (LANL), which is operated by Triad +# National Security, LLC for the U.S. Department of Energy/National +# Nuclear Security Administration. All rights in the program are +# reserved by Triad National Security, LLC, and the U.S. Department of +# Energy/National Nuclear Security Administration. The Government is +# granted for itself and others acting on its behalf a nonexclusive, +# paid-up, irrevocable worldwide license in this material to reproduce, +# prepare derivative works, distribute copies to the public, perform +# publicly and display publicly, and to permit others to do so. +# ====================================================================== + +matid=3337 +name=copper +rhomin=0. From 14a5303b27bad65039d5f36c4a070b51c51bb6b7 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 17 Jul 2024 13:31:54 -0600 Subject: [PATCH 44/61] Allocate lambda memory all at once and then use accessor to index into it --- test/test_closure_pte.cpp | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index bfd8b3fbe0..caf56a2e33 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -33,6 +33,7 @@ constexpr Real MJ_per_kg = 1.0e10; #ifdef SINGULARITY_TEST_SESAME #ifdef SINGULARITY_USE_SPINER_WITH_HDF5 +using singularity::miximpl::CacheAccessor; using singularity::DavisProducts; using singularity::IdealGas; using singularity::MAX_NUM_LAMBDAS; @@ -106,13 +107,10 @@ bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, const size_t scratch_bytes = pte_solver_scratch_size * sizeof(Real); Real *scratch = (double *)PORTABLE_MALLOC(scratch_bytes); - // Allocate lambdas for each EOS - constexpr size_t lambda_bytes = MAX_NUM_LAMBDAS * sizeof(Real); // EOS lambda size - const size_t lambda_size = num_eos * sizeof(Real *); // Size of lambda array - Real **lambdas = (Real **)PORTABLE_MALLOC(lambda_size); - for (auto i = 0; i < num_eos; i++) { - lambdas[i] = (Real *)PORTABLE_MALLOC(lambda_bytes); - } + // Allocate lambdas for all EOS and use an accessor to index into it + const size_t lambda_bytes = num_eos * MAX_NUM_LAMBDAS * sizeof(Real *); + Real *lambda_memory = (Real *)PORTABLE_MALLOC(lambda_bytes); + CacheAccessor lambdas = CacheAccessor(lambda_memory); // Solve the PTE system on device using a one-teration portableFor bool pte_converged; @@ -128,10 +126,6 @@ bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, portableCopyToHost(&pte_converged, pte_converged_d, bool_bytes); // Free temp memory - for (auto i = 0; i < num_eos; i++) { - // Free each lambda separately first - PORTABLE_FREE(lambdas[i]); - } PORTABLE_FREE(lambdas); // Free entire lambda array PORTABLE_FREE(scratch); From 73b6ea466aead368cf21af058bbda6fbfc9459f8 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 17 Jul 2024 13:32:28 -0600 Subject: [PATCH 45/61] Clang format --- test/test_closure_pte.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index caf56a2e33..80c3ea9e3e 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -33,7 +33,6 @@ constexpr Real MJ_per_kg = 1.0e10; #ifdef SINGULARITY_TEST_SESAME #ifdef SINGULARITY_USE_SPINER_WITH_HDF5 -using singularity::miximpl::CacheAccessor; using singularity::DavisProducts; using singularity::IdealGas; using singularity::MAX_NUM_LAMBDAS; @@ -41,6 +40,7 @@ using singularity::PTESolverRhoT; using singularity::PTESolverRhoTRequiredScratch; using singularity::ShiftedEOS; using singularity::SpinerEOSDependsRhoT; +using singularity::miximpl::CacheAccessor; using EOS = singularity::Variant, DavisProducts, SpinerEOSDependsRhoT>; From 200531ec460a11d05d63d314fdf6cf7c71c11c57 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 17 Jul 2024 13:35:04 -0600 Subject: [PATCH 46/61] fix typo --- test/test_closure_pte.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 80c3ea9e3e..9736530772 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -40,7 +40,7 @@ using singularity::PTESolverRhoT; using singularity::PTESolverRhoTRequiredScratch; using singularity::ShiftedEOS; using singularity::SpinerEOSDependsRhoT; -using singularity::miximpl::CacheAccessor; +using singularity::mix_impl::CacheAccessor; using EOS = singularity::Variant, DavisProducts, SpinerEOSDependsRhoT>; From 592c2d3693919c4814eed12078a08bd6795f24b6 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 17 Jul 2024 13:36:44 -0600 Subject: [PATCH 47/61] Free lambda memory, not the struct itself --- test/test_closure_pte.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 9736530772..4d1182c646 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -126,7 +126,7 @@ bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, portableCopyToHost(&pte_converged, pte_converged_d, bool_bytes); // Free temp memory - PORTABLE_FREE(lambdas); // Free entire lambda array + PORTABLE_FREE(lambda_memory); // Free entire lambda array PORTABLE_FREE(scratch); // Free PTE values From b68f524b716396662d98b40eb324494ae7de4984 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 17 Jul 2024 17:15:31 -0600 Subject: [PATCH 48/61] Ensure Finalize() is called from the host and combine tests --- test/test_closure_pte.cpp | 89 +++++++++++++++++++++------------------ 1 file changed, 47 insertions(+), 42 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 4d1182c646..2a0b96405c 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -43,7 +43,7 @@ using singularity::SpinerEOSDependsRhoT; using singularity::mix_impl::CacheAccessor; using EOS = singularity::Variant, DavisProducts, - SpinerEOSDependsRhoT>; + DavisReactants, SpinerEOSDependsRhoT>; template EOS *copy_eos_arr_to_device(const int num_eos, EOSArrT eos_arr) { @@ -58,13 +58,12 @@ EOS *copy_eos_arr_to_device(const int num_eos, EOSArrT eos_arr) { return v_EOS; } -void finalize_eos_arr(const int num_eos, EOS *v_EOS) { - // Call Finalize on each EOS - for (auto i = 0; i < num_eos; i++) { - v_EOS[i].Finalize(); +template +void finalize_eos_arr(EOSArrT eos_arr) { + // Call Finalize on each EOS on the host + for (auto eos : eos_arr) { + eos.Finalize(); } - // Free EOS memory - PORTABLE_FREE(v_EOS); } template @@ -141,22 +140,35 @@ bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, SCENARIO("Density-Temperature PTE Solver", "[PTE]") { - GIVEN("Three equations of state") { - // Set up the three EOS - constexpr int num_eos = 3; - std::array eos_arr; + GIVEN("Four equations of state") { + // A series of tests using some subset of these four EOS + constexpr int num_eos = 4; + // Reference state constexpr Real P0 = 1.0e6; // 1 bar constexpr Real T0 = 296; // K // Ideal gas air constexpr Real rho0_air = 1e-03; // g/cc constexpr Real Gruneisen_air = 0.4; - constexpr Real CV_air = P0 / rho0_air / (Gruneisen_air * T0); - eos_arr[0] = IdealGas(Gruneisen_air, CV_air); + constexpr Real Cv_air = P0 / rho0_air / (Gruneisen_air * T0); + air_eos = IdealGas(Gruneisen_air, Cv_air); // Spiner copper EOS constexpr int Cu_matid = 3337; const std::string eos_file = "../materials.sp5"; - eos_arr[1] = SpinerEOSDependsRhoT(eos_file, Cu_matid); + copper_eos = SpinerEOSDependsRhoT(eos_file, Cu_matid); + // Davis Reactants EOS + constexpr rho0_DP = 1.890; + constexpr e0_DP = 0.; // erg / g + constexpr P0_DP = 0.; // microbar + constexpr T0_DP = 297; // K + constexpr A = 1.8 * std::sqrt(GPa); // sqrt(microbar) + constexpr B = 4.6; + constexpr C = 0.34; + constexpr G0 = 0.56; + constexpr Z = 0.0; + constexpr alpha = 0.4265; + constexpr Cv_DP = 0.001074 * MJ_per_kg; // erg / g / K + davis_r_eos = DavisReactants(rho0_DP, e0_DP, P0_DP, T0_DP, A, B, C, G0, Z, alpha, Cv_DP); // Davis Products EOS constexpr Real a = 0.798311; constexpr Real b = 0.58; @@ -166,57 +178,50 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { constexpr Real pc = 3.2 * GPa; // microbar constexpr Real Cv = 0.001072 * MJ_per_kg; // erg / g / K constexpr Real Q = 4.115 * MJ_per_kg; - eos_arr[2] = DavisProducts(a, b, k, n, vc, pc, Cv).Modify(-Q); + davis_p_eos = DavisProducts(a, b, k, n, vc, pc, Cv).Modify(-Q); - GIVEN("A difficult state for the density-temperature PTE solver") { + GIVEN("A difficult three-material state") { // Define the state constexpr Real spvol_bulk = 6.256037280402093e-01; constexpr Real sie_bulk = -1.441692060406520e+10; - const std::array mass_frac = { + constexpr int num_pte = 3; + const std::array mass_frac = { 1.10382442033331e-10, 0.124935312146569, 0.875064687743048}; + std::array eos_arr = {air_eos.GetOnDevice(), copper_eos.GetOnDevice(), davis_p_eos.GetOnDevice()}; THEN("The PTE solver should converge") { - EOS *v_EOS = copy_eos_arr_to_device(num_eos, eos_arr); + EOS *v_EOS = copy_eos_arr_to_device(num_pte, eos_arr); const bool pte_converged = - run_PTE_from_state(num_eos, v_EOS, spvol_bulk, sie_bulk, mass_frac); + run_PTE_from_state(num_pte, v_EOS, spvol_bulk, sie_bulk, mass_frac); CHECK(pte_converged); - finalize_eos_arr(num_eos, v_EOS); + // Free EOS copies on device + PORTABLE_FREE(v_EOS); } + finalize_eos(eos_arr); } - } - GIVEN("Two equations of state") { - // Set up the three EOS - constexpr int num_eos = 2; - std::array eos_arr; - // Reference state - constexpr Real P0 = 1.0e6; // 1 bar - constexpr Real T0 = 296; // K - // Ideal gas air - constexpr Real rho0_air = 1e-03; // g/cc - constexpr Real Gruneisen_air = 0.4; - constexpr Real CV_air = P0 / rho0_air / (Gruneisen_air * T0); - eos_arr[0] = IdealGas(Gruneisen_air, CV_air); - // Spiner copper EOS - constexpr int Cu_matid = 3337; - const std::string eos_file = "../materials.sp5"; - eos_arr[1] = SpinerEOSDependsRhoT(eos_file, Cu_matid); - - GIVEN("A state that would cause a negative temperature in the PTE solver") { + GIVEN("A difficult two-material state") { + // TODO: make this test converge // Define the state constexpr Real spvol_bulk = 4.010467628234189e-01; constexpr Real sie_bulk = 3.290180957185173e+07; const std::array mass_frac = {0.000312273191678158, 0.999687726808322}; - + std::array eos_arr = {air_eos.GetOnDevice(), copper_eos.GetOnDevice()}; THEN("The PTE solver should converge") { EOS *v_EOS = copy_eos_arr_to_device(num_eos, eos_arr); const bool pte_converged = run_PTE_from_state(num_eos, v_EOS, spvol_bulk, sie_bulk, mass_frac); - // TODO: make this test converge // CHECK(pte_converged); - finalize_eos_arr(num_eos, v_EOS); + // Free EOS copies on device + PORTABLE_FREE(v_EOS); } + finalize_eos(eos_arr); } + // Clean up original EOS objects before they go out of scope + air_eos.Finalize(); // irrelevant because no data allocated + copper_eos.Finalize(); // actually does something + davis_r_eos.Finalize(); // irrelevant because no data allocated + davis_p_eos.Finalize(); // irrelevant because no data allocated } } #endif // SINGULARITY_USE_SPINER_WITH_HDF5 From e59140f9a4b97fd2f3862b0eef82be2c78650d8c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 17 Jul 2024 17:16:06 -0600 Subject: [PATCH 49/61] Clang format --- test/test_closure_pte.cpp | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 2a0b96405c..5fcd89e7f5 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -158,9 +158,9 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { copper_eos = SpinerEOSDependsRhoT(eos_file, Cu_matid); // Davis Reactants EOS constexpr rho0_DP = 1.890; - constexpr e0_DP = 0.; // erg / g - constexpr P0_DP = 0.; // microbar - constexpr T0_DP = 297; // K + constexpr e0_DP = 0.; // erg / g + constexpr P0_DP = 0.; // microbar + constexpr T0_DP = 297; // K constexpr A = 1.8 * std::sqrt(GPa); // sqrt(microbar) constexpr B = 4.6; constexpr C = 0.34; @@ -168,7 +168,8 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { constexpr Z = 0.0; constexpr alpha = 0.4265; constexpr Cv_DP = 0.001074 * MJ_per_kg; // erg / g / K - davis_r_eos = DavisReactants(rho0_DP, e0_DP, P0_DP, T0_DP, A, B, C, G0, Z, alpha, Cv_DP); + davis_r_eos = + DavisReactants(rho0_DP, e0_DP, P0_DP, T0_DP, A, B, C, G0, Z, alpha, Cv_DP); // Davis Products EOS constexpr Real a = 0.798311; constexpr Real b = 0.58; @@ -187,7 +188,8 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { constexpr int num_pte = 3; const std::array mass_frac = { 1.10382442033331e-10, 0.124935312146569, 0.875064687743048}; - std::array eos_arr = {air_eos.GetOnDevice(), copper_eos.GetOnDevice(), davis_p_eos.GetOnDevice()}; + std::array eos_arr = {air_eos.GetOnDevice(), copper_eos.GetOnDevice(), + davis_p_eos.GetOnDevice()}; THEN("The PTE solver should converge") { EOS *v_EOS = copy_eos_arr_to_device(num_pte, eos_arr); @@ -206,7 +208,8 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { constexpr Real sie_bulk = 3.290180957185173e+07; const std::array mass_frac = {0.000312273191678158, 0.999687726808322}; - std::array eos_arr = {air_eos.GetOnDevice(), copper_eos.GetOnDevice()}; + std::array eos_arr = {air_eos.GetOnDevice(), + copper_eos.GetOnDevice()}; THEN("The PTE solver should converge") { EOS *v_EOS = copy_eos_arr_to_device(num_eos, eos_arr); const bool pte_converged = @@ -218,8 +221,8 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { finalize_eos(eos_arr); } // Clean up original EOS objects before they go out of scope - air_eos.Finalize(); // irrelevant because no data allocated - copper_eos.Finalize(); // actually does something + air_eos.Finalize(); // irrelevant because no data allocated + copper_eos.Finalize(); // actually does something davis_r_eos.Finalize(); // irrelevant because no data allocated davis_p_eos.Finalize(); // irrelevant because no data allocated } From 53f419851b1de0fadc0b04c0b1e0c6c7e5e1a1ad Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 17 Jul 2024 17:20:32 -0600 Subject: [PATCH 50/61] Fix silly bugs --- test/test_closure_pte.cpp | 38 ++++++++++++++++++++------------------ 1 file changed, 20 insertions(+), 18 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 5fcd89e7f5..fc97935a40 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -34,6 +34,7 @@ constexpr Real MJ_per_kg = 1.0e10; #ifdef SINGULARITY_USE_SPINER_WITH_HDF5 using singularity::DavisProducts; +using singularity::DavisReactants; using singularity::IdealGas; using singularity::MAX_NUM_LAMBDAS; using singularity::PTESolverRhoT; @@ -151,24 +152,24 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { constexpr Real rho0_air = 1e-03; // g/cc constexpr Real Gruneisen_air = 0.4; constexpr Real Cv_air = P0 / rho0_air / (Gruneisen_air * T0); - air_eos = IdealGas(Gruneisen_air, Cv_air); + EOS air_eos = IdealGas(Gruneisen_air, Cv_air); // Spiner copper EOS constexpr int Cu_matid = 3337; const std::string eos_file = "../materials.sp5"; - copper_eos = SpinerEOSDependsRhoT(eos_file, Cu_matid); + EOS copper_eos = SpinerEOSDependsRhoT(eos_file, Cu_matid); // Davis Reactants EOS - constexpr rho0_DP = 1.890; - constexpr e0_DP = 0.; // erg / g - constexpr P0_DP = 0.; // microbar - constexpr T0_DP = 297; // K - constexpr A = 1.8 * std::sqrt(GPa); // sqrt(microbar) - constexpr B = 4.6; - constexpr C = 0.34; - constexpr G0 = 0.56; - constexpr Z = 0.0; - constexpr alpha = 0.4265; - constexpr Cv_DP = 0.001074 * MJ_per_kg; // erg / g / K - davis_r_eos = + constexpr Real rho0_DP = 1.890; + constexpr Real e0_DP = 0.; // erg / g + constexpr Real P0_DP = 0.; // microbar + constexpr Real T0_DP = 297; // K + constexpr Real A = 1.8 * std::sqrt(GPa); // sqrt(microbar) + constexpr Real B = 4.6; + constexpr Real C = 0.34; + constexpr Real G0 = 0.56; + constexpr Real Z = 0.0; + constexpr Real alpha = 0.4265; + constexpr Real Cv_DP = 0.001074 * MJ_per_kg; // erg / g / K + EOS davis_r_eos = DavisReactants(rho0_DP, e0_DP, P0_DP, T0_DP, A, B, C, G0, Z, alpha, Cv_DP); // Davis Products EOS constexpr Real a = 0.798311; @@ -179,7 +180,7 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { constexpr Real pc = 3.2 * GPa; // microbar constexpr Real Cv = 0.001072 * MJ_per_kg; // erg / g / K constexpr Real Q = 4.115 * MJ_per_kg; - davis_p_eos = DavisProducts(a, b, k, n, vc, pc, Cv).Modify(-Q); + EOS davis_p_eos = DavisProducts(a, b, k, n, vc, pc, Cv).Modify(-Q); GIVEN("A difficult three-material state") { // Define the state @@ -199,14 +200,15 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { // Free EOS copies on device PORTABLE_FREE(v_EOS); } - finalize_eos(eos_arr); + finalize_eos_arr(eos_arr); } GIVEN("A difficult two-material state") { // TODO: make this test converge // Define the state constexpr Real spvol_bulk = 4.010467628234189e-01; constexpr Real sie_bulk = 3.290180957185173e+07; - const std::array mass_frac = {0.000312273191678158, + constexpr int num_pte = 2; + const std::array mass_frac = {0.000312273191678158, 0.999687726808322}; std::array eos_arr = {air_eos.GetOnDevice(), copper_eos.GetOnDevice()}; @@ -218,7 +220,7 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { // Free EOS copies on device PORTABLE_FREE(v_EOS); } - finalize_eos(eos_arr); + finalize_eos_arr(eos_arr); } // Clean up original EOS objects before they go out of scope air_eos.Finalize(); // irrelevant because no data allocated From e9db000cfd29e4d8209b7871dec3eb5f36c58835 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 17 Jul 2024 19:42:50 -0600 Subject: [PATCH 51/61] Switch to using SINGULARITY_USE_SPINER_WITH_HDF5 for closure test --- test/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index f67d60e0b3..daf0cc8046 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -74,9 +74,9 @@ if(SINGULARITY_TEST_STELLAR_COLLAPSE) target_compile_definitions(eos_tabulated_unit_tests PRIVATE SINGULARITY_TEST_STELLAR_COLLAPSE) endif() -if(SINGULARITY_USE_SPINER) +if(SINGULARITY_USE_SPINER_WITH_HDF5) target_compile_definitions(closure_unit_tests - PRIVATE SINGULARITY_USE_SPINER) + PRIVATE SINGULARITY_USE_SPINER_WITH_HDF5) endif() if(SINGULARITY_BUILD_CLOSURE) target_compile_definitions(closure_unit_tests From e9dbf6e8ab8f33c42dca42ca28a713d6e4984bbb Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 17 Jul 2024 19:43:13 -0600 Subject: [PATCH 52/61] Make parameter constexpr --- test/test_closure_pte.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index fc97935a40..36a38a35ac 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -16,6 +16,7 @@ #include #include #include +#include #ifndef CATCH_CONFIG_FAST_COMPILE #define CATCH_CONFIG_FAST_COMPILE @@ -162,7 +163,7 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { constexpr Real e0_DP = 0.; // erg / g constexpr Real P0_DP = 0.; // microbar constexpr Real T0_DP = 297; // K - constexpr Real A = 1.8 * std::sqrt(GPa); // sqrt(microbar) + constexpr Real A = 1.8 * 1.0e5; // 1.8 * sqrt(GPa) -> sqrt(microbar) constexpr Real B = 4.6; constexpr Real C = 0.34; constexpr Real G0 = 0.56; From 3adea33ed4e964f5f7e82ddc1068046524ddeff6 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 17 Jul 2024 19:43:43 -0600 Subject: [PATCH 53/61] Clang format --- test/test_closure_pte.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 36a38a35ac..1e7b5a8ca2 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -14,9 +14,9 @@ #ifdef SINGULARITY_BUILD_CLOSURE #include +#include #include #include -#include #ifndef CATCH_CONFIG_FAST_COMPILE #define CATCH_CONFIG_FAST_COMPILE @@ -160,9 +160,9 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { EOS copper_eos = SpinerEOSDependsRhoT(eos_file, Cu_matid); // Davis Reactants EOS constexpr Real rho0_DP = 1.890; - constexpr Real e0_DP = 0.; // erg / g - constexpr Real P0_DP = 0.; // microbar - constexpr Real T0_DP = 297; // K + constexpr Real e0_DP = 0.; // erg / g + constexpr Real P0_DP = 0.; // microbar + constexpr Real T0_DP = 297; // K constexpr Real A = 1.8 * 1.0e5; // 1.8 * sqrt(GPa) -> sqrt(microbar) constexpr Real B = 4.6; constexpr Real C = 0.34; From f32ca9044848c6c57d91b02a51ca3925d2cb6763 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 18 Jul 2024 15:14:23 -0600 Subject: [PATCH 54/61] Enforce convention that '*_TEST_*' variables are defined in the test CMakeLists and other '*_USE_*' or '*_BUILD_*' variables are in the main project --- CMakeLists.txt | 4 +++- test/CMakeLists.txt | 8 -------- 2 files changed, 3 insertions(+), 9 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 1a7ec465ba..8eea80ed30 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -258,7 +258,6 @@ if(SINGULARITY_USE_HIGH_RISK_MATH) target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_USE_HIGH_RISK_MATH) endif() - if(SINGULARITY_TEST_SESAME) target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_TEST_SESAME) endif() @@ -268,6 +267,9 @@ endif() if(SINGULARITY_USE_HELMHOLTZ) target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_USE_HELMHOLTZ) endif() +if(SINGULARITY_USE_SPINER_WITH_HDF5) + target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_USE_SPINER_WITH_HDF5) +endif() # ------------------------------------------------------------------------------# # Handle dependencies diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index daf0cc8046..4c654975ec 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -74,14 +74,6 @@ if(SINGULARITY_TEST_STELLAR_COLLAPSE) target_compile_definitions(eos_tabulated_unit_tests PRIVATE SINGULARITY_TEST_STELLAR_COLLAPSE) endif() -if(SINGULARITY_USE_SPINER_WITH_HDF5) - target_compile_definitions(closure_unit_tests - PRIVATE SINGULARITY_USE_SPINER_WITH_HDF5) -endif() -if(SINGULARITY_BUILD_CLOSURE) - target_compile_definitions(closure_unit_tests - PRIVATE SINGULARITY_BUILD_CLOSURE) -endif() target_link_libraries(eos_analytic_unit_tests PRIVATE Catch2::Catch2 singularity-eos::singularity-eos) From 4588d76cc425da143de769b896c1c176cae14e8c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 19 Jul 2024 16:13:28 -0600 Subject: [PATCH 55/61] Make comments more descriptive --- singularity-eos/closure/mixed_cell_models.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 88af8d5b55..202a072060 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -1470,13 +1470,15 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { scale = 0.5; Real err_mid = s.TestUpdate(scale); if (err_mid < err && err_mid < err_old) { - // try a larger step since a half step reduces the error + // We know the half step is better than both the full step and the + // prior result, so try a pseudo parabolic fit to the error and find its + // minimum. The `scale` value is bound between 0.75 and 0.25. scale = 0.75 + 0.5 * robust::ratio(err_mid - err, err - 2.0 * err_mid + err_old); } for (int line_iter = 0; line_iter < line_search_max_iter; line_iter++) { err = s.TestUpdate(scale); if (err < err_old + line_search_alpha * scale * gradfdx) break; - // shrink the step if the error isn't reduced + // shrink the step if the error isn't reduced enough scale *= line_search_fac; } } From 7d35c432bf0877e4e15449584f8e36cf33a04254 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 19 Jul 2024 17:59:13 -0600 Subject: [PATCH 56/61] Don't call GetOnDevice() twice --- test/test_closure_pte.cpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 1e7b5a8ca2..2c6e2be681 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -49,11 +49,8 @@ using EOS = singularity::Variant, DavisProdu template EOS *copy_eos_arr_to_device(const int num_eos, EOSArrT eos_arr) { - // Move EOS array from host to device + // Assumes that GetOnDevice() has already been called for each EOS in eos_arr const size_t EOS_bytes = num_eos * sizeof(EOS); - for (auto i = 0; i < num_eos; i++) { - eos_arr[i] = eos_arr[i].GetOnDevice(); - } EOS *v_EOS = (EOS *)PORTABLE_MALLOC(EOS_bytes); const size_t bytes = num_eos * sizeof(EOS); portableCopyToDevice(v_EOS, eos_arr.data(), bytes); @@ -221,6 +218,7 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { // Free EOS copies on device PORTABLE_FREE(v_EOS); } + // Deallocates device memory for each EOS (if applicable) finalize_eos_arr(eos_arr); } // Clean up original EOS objects before they go out of scope From fd07bdd7e3d5c3da8d4188ac036f821eec8fc6c2 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 19 Jul 2024 18:43:08 -0600 Subject: [PATCH 57/61] Fix array out of bounds issues when confusing num_pte and num_eos --- test/test_closure_pte.cpp | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 2c6e2be681..e2589d04ce 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -66,17 +66,17 @@ void finalize_eos_arr(EOSArrT eos_arr) { } template -bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, +bool run_PTE_from_state(const int num_pte, EOS *v_EOS, const Real spvol_bulk, const Real sie_bulk, ArrT mass_frac) { // Calculate material densities (and corresponding energies) and the total // volume fraction - std::vector vol_frac(num_eos); - std::vector densities(num_eos); - std::vector sies(num_eos); - std::vector temperatures(num_eos); - std::vector pressures(num_eos); + std::vector vol_frac(num_pte); + std::vector densities(num_pte); + std::vector sies(num_pte); + std::vector temperatures(num_pte); + std::vector pressures(num_pte); Real vfrac_sum = 0; - for (auto i = 0; i < num_eos; ++i) { + for (auto i = 0; i < num_pte; ++i) { // Initial guess: all materials at the cell density vol_frac[i] = mass_frac[i]; densities[i] = mass_frac[i] / spvol_bulk / vol_frac[i]; @@ -88,7 +88,7 @@ bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, } // Copy values to device (when available) - const size_t bytes = num_eos * sizeof(Real); + const size_t bytes = num_pte * sizeof(Real); Real *v_densities = (Real *)PORTABLE_MALLOC(bytes); portableCopyToDevice(v_densities, densities.data(), bytes); Real *v_vol_frac = (Real *)PORTABLE_MALLOC(bytes); @@ -101,12 +101,12 @@ bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, portableCopyToDevice(v_pressures, pressures.data(), bytes); // Allocate scratch space for the PTE solver - const int pte_solver_scratch_size = PTESolverRhoTRequiredScratch(num_eos); + const int pte_solver_scratch_size = PTESolverRhoTRequiredScratch(num_pte); const size_t scratch_bytes = pte_solver_scratch_size * sizeof(Real); Real *scratch = (double *)PORTABLE_MALLOC(scratch_bytes); // Allocate lambdas for all EOS and use an accessor to index into it - const size_t lambda_bytes = num_eos * MAX_NUM_LAMBDAS * sizeof(Real *); + const size_t lambda_bytes = num_pte * MAX_NUM_LAMBDAS * sizeof(Real *); Real *lambda_memory = (Real *)PORTABLE_MALLOC(lambda_bytes); CacheAccessor lambdas = CacheAccessor(lambda_memory); @@ -117,7 +117,7 @@ bool run_PTE_from_state(const int num_eos, EOS *v_EOS, const Real spvol_bulk, portableFor( "Device execution of PTE Test", 0, 1, PORTABLE_LAMBDA(int i) { PTESolverRhoT method( - num_eos, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, + num_pte, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, v_temperatures, v_pressures, lambdas, scratch); pte_converged_d[0] = PTESolver(method); }); @@ -211,9 +211,9 @@ SCENARIO("Density-Temperature PTE Solver", "[PTE]") { std::array eos_arr = {air_eos.GetOnDevice(), copper_eos.GetOnDevice()}; THEN("The PTE solver should converge") { - EOS *v_EOS = copy_eos_arr_to_device(num_eos, eos_arr); + EOS *v_EOS = copy_eos_arr_to_device(num_pte, eos_arr); const bool pte_converged = - run_PTE_from_state(num_eos, v_EOS, spvol_bulk, sie_bulk, mass_frac); + run_PTE_from_state(num_pte, v_EOS, spvol_bulk, sie_bulk, mass_frac); // CHECK(pte_converged); // Free EOS copies on device PORTABLE_FREE(v_EOS); From 6de9259a86dff2a57d22897a1feee377bd4c541e Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 19 Jul 2024 18:45:24 -0600 Subject: [PATCH 58/61] Remove unused num_eos --- test/test_closure_pte.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index e2589d04ce..7a1e4e8af6 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -140,9 +140,6 @@ bool run_PTE_from_state(const int num_pte, EOS *v_EOS, const Real spvol_bulk, SCENARIO("Density-Temperature PTE Solver", "[PTE]") { GIVEN("Four equations of state") { - // A series of tests using some subset of these four EOS - constexpr int num_eos = 4; - // Reference state constexpr Real P0 = 1.0e6; // 1 bar constexpr Real T0 = 296; // K From 8a8ac4e2760a26a600a1c8ee7c34bee3647b730f Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 22 Jul 2024 12:52:54 -0600 Subject: [PATCH 59/61] Update copyright --- sesame2spiner/examples/unit_tests/copper.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/sesame2spiner/examples/unit_tests/copper.dat b/sesame2spiner/examples/unit_tests/copper.dat index 3ccdcdc544..16ba63d428 100644 --- a/sesame2spiner/examples/unit_tests/copper.dat +++ b/sesame2spiner/examples/unit_tests/copper.dat @@ -2,7 +2,7 @@ # Example input deck for sesame2spiner, # a tool for converting eospac to spiner # Author: Jonah Miller (jonahm@lanl.gov) -# © 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 From a3c4d6e7eaa33e6a0033eb6232eeaeea2c93b15b Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 22 Jul 2024 12:59:58 -0600 Subject: [PATCH 60/61] Update copyright --- sesame2spiner/examples/unit_tests/copper.dat | 2 +- test/test_closure_pte.cpp | 2 +- test/test_get_sg_eos.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/sesame2spiner/examples/unit_tests/copper.dat b/sesame2spiner/examples/unit_tests/copper.dat index 3ccdcdc544..16ba63d428 100644 --- a/sesame2spiner/examples/unit_tests/copper.dat +++ b/sesame2spiner/examples/unit_tests/copper.dat @@ -2,7 +2,7 @@ # Example input deck for sesame2spiner, # a tool for converting eospac to spiner # Author: Jonah Miller (jonahm@lanl.gov) -# © 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/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 7a1e4e8af6..ca01906c0c 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2024. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National diff --git a/test/test_get_sg_eos.cpp b/test/test_get_sg_eos.cpp index 082c035011..7a9ea0c6b6 100644 --- a/test/test_get_sg_eos.cpp +++ b/test/test_get_sg_eos.cpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2024. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National From 454039a1aef6fd38cb3892d1b7981b0d4f93eb40 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 22 Jul 2024 16:22:01 -0600 Subject: [PATCH 61/61] Remove asserts in Ideal calculation in favor of a positive temperature check to make sure the guess is valid --- singularity-eos/closure/mixed_cell_models.hpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index 202a072060..f8c0335e38 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -353,12 +353,9 @@ class PTESolverBase { Asum += vfrac[m] * robust::ratio(press[m], temp[m]); rhoBsum += rho[m] * vfrac[m] * robust::ratio(sie[m], temp[m]); } - PORTABLE_REQUIRE(Tnorm > 0., "Non-positive Ideal PTE temperature guess"); Asum *= uscale / Tnorm; rhoBsum /= Tnorm; - PORTABLE_REQUIRE(rhoBsum > 0., "Non-positive Ideal PTE energy density guess"); Tideal = uscale / rhoBsum / Tnorm; - PORTABLE_REQUIRE(vfrac_total > 0., "Non-positive Ideal PTE volume fraction sum"); Pideal = robust::ratio(Tnorm * Tideal * Asum, uscale * vfrac_total); } @@ -391,7 +388,8 @@ class PTESolverBase { const Real alpha = robust::ratio(Pideal, Tideal); for (int m = 0; m < nmat; ++m) { vfrac[m] *= robust::ratio(press[m], (temp[m] * alpha)); - if (robust::ratio(rhobar[m], vfrac[m]) < eos[m].RhoPmin(Tnorm * Tideal)) { + if (Tnorm * Tideal < 0 || + robust::ratio(rhobar[m], vfrac[m]) < eos[m].RhoPmin(Tnorm * Tideal)) { // abort because this is putting this material into a bad state for (int n = m; n >= 0; n--) vfrac[n] = vtemp[n];