Skip to content

Commit

Permalink
Merge pull request #383 from lanl/jhp/robust_pte
Browse files Browse the repository at this point in the history 
Fix Bug in PTE Solver and Several Code Quality Improvements
  • Loading branch information
@jhp-lanl
jhp-lanl authored Jul 23, 2024
2 parents 03ccf13 + ff35cea commit f5e52ba
Show file tree
Hide file tree
Showing 12 changed files with 362 additions and 47 deletions.
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 3 additions & 1 deletion CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,7 +259,6 @@ if(SINGULARITY_USE_HIGH_RISK_MATH)
target_compile_definitions(singularity-eos_Interface
INTERFACE SINGULARITY_USE_HIGH_RISK_MATH)
endif()

if(SINGULARITY_TEST_SESAME)
target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_TEST_SESAME)
endif()
Expand All @@ -269,6 +268,9 @@ endif()
if(SINGULARITY_USE_HELMHOLTZ)
target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_USE_HELMHOLTZ)
endif()
if(SINGULARITY_USE_SPINER_WITH_HDF5)
target_compile_definitions(singularity-eos_Interface INTERFACE SINGULARITY_USE_SPINER_WITH_HDF5)
endif()

# ------------------------------------------------------------------------------#
# Handle dependencies
Expand Down
20 changes: 20 additions & 0 deletions sesame2spiner/examples/unit_tests/copper.dat
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,20 @@
# ======================================================================
# Example input deck for sesame2spiner,
# a tool for converting eospac to spiner
# Author: Jonah Miller ([email protected])
# © 2021-2024. Triad National Security, LLC. All rights reserved. This
# program was produced under U.S. Government contract 89233218CNA000001
# for Los Alamos National Laboratory (LANL), which is operated by Triad
# National Security, LLC for the U.S. Department of Energy/National
# Nuclear Security Administration. All rights in the program are
# reserved by Triad National Security, LLC, and the U.S. Department of
# Energy/National Nuclear Security Administration. The Government is
# granted for itself and others acting on its behalf a nonexclusive,
# paid-up, irrevocable worldwide license in this material to reproduce,
# prepare derivative works, distribute copies to the public, perform
# publicly and display publicly, and to permit others to do so.
# ======================================================================

matid=3337
name=copper
rhomin=0.
93 changes: 67 additions & 26 deletions singularity-eos/closure/mixed_cell_models.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,9 +170,9 @@ class PTESolverBase {
PORTABLE_INLINE_FUNCTION
virtual void Fixup() const {
Real vsum = 0;
for (int m = 0; m < nmat; ++m)
for (int m = 0; m < nmat; ++m) {
vsum += vfrac[m];
PORTABLE_REQUIRE(vsum > 0., "Volume fraction sum is non-positive");
}
for (int m = 0; m < nmat; ++m)
vfrac[m] *= robust::ratio(vfrac_total, vsum);
}
Expand All @@ -182,7 +182,6 @@ class PTESolverBase {
void Finalize() {
for (int m = 0; m < nmat; m++) {
temp[m] *= Tnorm;
PORTABLE_REQUIRE(temp[m] >= 0., "Non-positive temperature returned");
u[m] *= uscale;
press[m] *= uscale;
}
Expand Down Expand Up @@ -219,8 +218,9 @@ class PTESolverBase {
// material m averaged over the full PTE volume
rho_total = 0.0;
for (int m = 0; m < nmat; ++m) {
PORTABLE_REQUIRE(vfrac[m] > 0., "Non-positive volume fraction");
PORTABLE_REQUIRE(rho[m] > 0., "Non-positive density");
PORTABLE_REQUIRE(vfrac[m] > 0.,
"Non-positive volume fraction provided to PTE solver");
PORTABLE_REQUIRE(rho[m] > 0., "Non-positive density provided to PTE solver");
rhobar[m] = rho[m] * vfrac[m];
rho_total += rhobar[m];
}
Expand Down Expand Up @@ -353,12 +353,9 @@ class PTESolverBase {
Asum += vfrac[m] * robust::ratio(press[m], temp[m]);
rhoBsum += rho[m] * vfrac[m] * robust::ratio(sie[m], temp[m]);
}
PORTABLE_REQUIRE(Tnorm > 0., "Non-positive temperature guess");
Asum *= uscale / Tnorm;
rhoBsum /= Tnorm;
PORTABLE_REQUIRE(rhoBsum > 0., "Non-positive energy density");
Tideal = uscale / rhoBsum / Tnorm;
PORTABLE_REQUIRE(vfrac_total > 0., "Non-positive volume fraction sum");
Pideal = robust::ratio(Tnorm * Tideal * Asum, uscale * vfrac_total);
}

Expand Down Expand Up @@ -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];
Expand Down Expand Up @@ -671,6 +669,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
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) {
Expand All @@ -684,8 +683,15 @@ class PTESolverRhoT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
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
Expand All @@ -701,8 +707,10 @@ class PTESolverRhoT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
// 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];
Expand Down Expand Up @@ -879,6 +887,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
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) {
Expand All @@ -890,6 +899,13 @@ class PTESolverFixedT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
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]);
}
Expand All @@ -903,8 +919,10 @@ class PTESolverFixedT : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
// 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];
}
Expand Down Expand Up @@ -1093,6 +1111,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
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) {
Expand All @@ -1106,6 +1125,13 @@ class PTESolverFixedP : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
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]);
}
Expand All @@ -1123,8 +1149,10 @@ class PTESolverFixedP : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer>
// 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];
Expand Down Expand Up @@ -1332,6 +1360,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
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
Expand All @@ -1348,6 +1377,13 @@ class PTESolverRhoU : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
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]);
}
Expand All @@ -1361,8 +1397,10 @@ class PTESolverRhoU : public mix_impl::PTESolverBase<EOSIndexer, RealIndexer> {
// 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];
Expand Down Expand Up @@ -1417,25 +1455,28 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) {

// possibly scale the update to stay within reasonable bounds
Real scale = s.ScaleDx();
// const Real scale_save = scale;
PORTABLE_REQUIRE(scale <= 1.0, "PTE Solver is attempting to increase the step size");

// Line search
Real gradfdx = -2.0 * scale * err;
scale = 1.0;
scale = 1.0; // New scale for line search
Real err_old = err;
err = s.TestUpdate(scale);
// Test the update and reset the cache the current state
err = s.TestUpdate(scale, true /* cache_state */);
if (err > err_old + line_search_alpha * gradfdx) {
// backtrack
Real err_mid = s.TestUpdate(0.5);
// backtrack to middle of step
scale = 0.5;
Real err_mid = s.TestUpdate(scale);
if (err_mid < err && err_mid < err_old) {
// We know the half step is better than both the full step and the
// prior result, so try a pseudo parabolic fit to the error and find its
// minimum. The `scale` value is bound between 0.75 and 0.25.
scale = 0.75 + 0.5 * robust::ratio(err_mid - err, err - 2.0 * err_mid + err_old);
} else {
scale = line_search_fac;
}

for (int line_iter = 0; line_iter < line_search_max_iter; line_iter++) {
err = s.TestUpdate(scale);
if (err < err_old + line_search_alpha * scale * gradfdx) break;
// shrink the step if the error isn't reduced enough
scale *= line_search_fac;
}
}
Expand Down
19 changes: 13 additions & 6 deletions singularity-eos/eos/get_sg_eos_functors.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,32 +63,37 @@ struct init_functor {

PORTABLE_INLINE_FUNCTION
void operator()(const int i, const int tid, double &mass_sum, int &npte,
const Real t_mult, const Real s_mult, const Real p_mult) const {
double &vfrac_sum, const Real t_mult, const Real s_mult,
const Real p_mult) const {
/* normalize mass fractions */
/* first find the mass sum */
/* also set idxs as the decrement of the eos offsets */
/* to take into account 1 based indexing in fortran */
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) {
// participating materials are those with non-negligible mass fractions
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;
Expand All @@ -98,10 +103,12 @@ struct init_functor {
// NOTE: the volume fractions and densities need to be consistent with the
// total specific volume since they are used to calculate internal
// quantities for the PTE solver
vfrac_sum = 0.;
for (int mp = 0; mp < npte; ++mp) {
const int m = pte_mats(tid, mp);
// Need to guess volume fractions
vfrac_pte(tid, mp) = frac_mass_v(i, m);
vfrac_sum += vfrac_pte(tid, mp);
// Calculate densities to be consistent with these volume fractions
rho_pte(tid, mp) = frac_mass_v(i, m) / spvol_v(i) / vfrac_pte(tid, mp);
temp_pte(tid, mp) = temp_v(i) * ev2k * t_mult;
Expand Down
5 changes: 3 additions & 2 deletions singularity-eos/eos/get_sg_eos_rho_e.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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<singularity::EOSAccessor_, Real *, Real **> 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);
Expand Down
9 changes: 5 additions & 4 deletions singularity-eos/eos/get_sg_eos_rho_p.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand All @@ -54,9 +55,9 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v,
// eos accessor
singularity::EOSAccessor_ eos_inx(eos_v, &pte_idxs(tid, 0));
PTESolverFixedP<singularity::EOSAccessor_, Real *, Real *> method(
npte, eos_inx, 1.0, press_pte(tid, 0), &rho_pte(tid, 0), &vfrac_pte(tid, 0),
&sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache[0],
&solver_scratch(tid, 0));
npte, eos_inx, vfrac_sum, press_pte(tid, 0), &rho_pte(tid, 0),
&vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0),
cache[0], &solver_scratch(tid, 0));
pte_converged = PTESolver(method);
// calculate total sie
for (int mp = 0; mp < npte; ++mp) {
Expand Down
Loading

0 comments on commit f5e52ba

Please sign in to comment.