From 2686c4fec1df2e276c0189e7167deb23e24b422c Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 15 Oct 2024 11:19:49 -0600 Subject: [PATCH 1/7] convice compiler there's no local memory addresses being returned --- singularity-eos/eos/eos_base.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 8b305dba54..daa55bbdb3 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -845,7 +845,9 @@ class EosBase { inline constexpr decltype(auto) GetUnmodifiedObject() { if constexpr (CRTP::IsModified()) { - return ((static_cast(this))->UnmodifyOnce()).GetUnmodifiedObject(); + auto unmodified = + ((static_cast(this))->UnmodifyOnce()).GetUnmodifiedObject(); + return unmodified; } else { return *static_cast(this); } From 8dceffcf2d6ac1c506f23be0046963780b2a6989 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 15 Oct 2024 11:20:37 -0600 Subject: [PATCH 2/7] update PTE solver to report instrumentation and expose parameters --- singularity-eos/closure/mixed_cell_models.hpp | 441 +++++++++--------- test/CMakeLists.txt | 1 + test/test_closure_pte.cpp | 3 +- test/test_eos_modifiers.cpp | 113 ----- test/test_eos_modifiers_minimal.cpp | 161 +++++++ test/test_pte.cpp | 4 +- test/test_pte_2phase.cpp | 4 +- test/test_pte_3phase.cpp | 4 +- 8 files changed, 403 insertions(+), 328 deletions(-) create mode 100644 test/test_eos_modifiers_minimal.cpp diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index ace00b8acb..69cdf80807 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -37,27 +37,32 @@ namespace singularity { // Implementation details below // ====================================================================== -// TODO(JCD): some of these should be exposed to consumers to allow changes to defaults -namespace mix_params { -constexpr Real derivative_eps = 3.0e-6; -constexpr Real pte_rel_tolerance_p = 1.e-6; -constexpr Real pte_rel_tolerance_e = 1.e-6; -constexpr Real pte_rel_tolerance_t = 1.e-4; -constexpr Real pte_abs_tolerance_p = 0.0; -constexpr Real pte_abs_tolerance_e = 1.e-4; -constexpr Real pte_abs_tolerance_t = 0.0; -constexpr Real pte_residual_tolerance = 1.e-8; -constexpr int pte_max_iter_per_mat = 128; -constexpr Real line_search_alpha = 1.e-2; -constexpr int line_search_max_iter = 6; -constexpr Real line_search_fac = 0.5; -constexpr Real vfrac_safety_fac = 0.95; -constexpr Real minimum_temperature = 1.e-9; -constexpr Real maximum_temperature = 1.e9; -constexpr Real temperature_limit = 1.0e15; -constexpr Real default_tguess = 300.; -constexpr Real min_dtde = 1.0e-16; -} // namespace mix_params +struct MixParams { + bool verbose = false; + Real derivative_eps = 3.0e-6; + Real pte_rel_tolerance_p = 1.e-6; + Real pte_rel_tolerance_e = 1.e-6; + Real pte_rel_tolerance_t = 1.e-4; + Real pte_abs_tolerance_p = 0.0; + Real pte_abs_tolerance_e = 1.e-4; + Real pte_abs_tolerance_t = 0.0; + Real pte_residual_tolerance = 1.e-8; + std::size_t pte_max_iter_per_mat = 128; + Real line_search_alpha = 1.e-2; + std::size_t line_search_max_iter = 6; + Real line_search_fac = 0.5; + Real vfrac_safety_fac = 0.95; + Real temperature_limit = 1.0e15; + Real default_tguess = 300.; + Real min_dtde = 1.0e-16; +}; + +struct SolverStatus { + bool converged = false; + std::size_t max_niter = 0; + std::size_t max_line_niter = 0; + Real residual; +}; namespace mix_impl { template class PTESolverBase { public: PTESolverBase() = delete; - PORTABLE_INLINE_FUNCTION int Nmat() const { return nmat; } - PORTABLE_INLINE_FUNCTION int &Niter() { return niter; } + PORTABLE_INLINE_FUNCTION std::size_t Nmat() const { return nmat; } + PORTABLE_INLINE_FUNCTION std::size_t &Niter() { return niter; } // Fixup is meant to be a hook for derived classes to provide arbitrary manipulations // after each iteration of the Newton solver. This version just renormalizes the // volume fractions, which is useful to deal with roundoff error. PORTABLE_INLINE_FUNCTION virtual void Fixup() const { Real vsum = 0; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vsum += vfrac[m]; } - for (int m = 0; m < nmat; ++m) + for (std::size_t m = 0; m < nmat; ++m) vfrac[m] *= robust::ratio(vfrac_total, vsum); } // Finalize restores the temperatures, energies, and pressures to unscaled values from // the internally scaled quantities used by the solvers PORTABLE_INLINE_FUNCTION void Finalize() { - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { temp[m] *= Tnorm; u[m] *= uscale; press[m] *= uscale; @@ -189,20 +194,24 @@ class PTESolverBase { // Solve the linear system for the update dx PORTABLE_INLINE_FUNCTION bool Solve() const { - for (int m = 0; m < neq; m++) + for (std::size_t m = 0; m < neq; m++) dx[m] = residual[m]; return solve_Ax_b_wscr(neq, jacobian, dx, sol_scratch); } + // Parameters for the solver + PORTABLE_INLINE_FUNCTION + const MixParams &GetParams() const { return params_; } protected: PORTABLE_INLINE_FUNCTION - PTESolverBase(int nmats, int neqs, const EOSIndexer &eos_, const Real vfrac_tot, - const Real sie_tot, const RealIndexer &rho_, const RealIndexer &vfrac_, - const RealIndexer &sie_, const RealIndexer &temp_, - const RealIndexer &press_, Real *&scratch, Real Tguess) + PTESolverBase(std::size_t nmats, std::size_t neqs, const EOSIndexer &eos_, + const Real vfrac_tot, const Real sie_tot, const RealIndexer &rho_, + const RealIndexer &vfrac_, const RealIndexer &sie_, + const RealIndexer &temp_, const RealIndexer &press_, Real *&scratch, + Real Tnorm, const MixParams ¶ms = MixParams()) : nmat(nmats), neq(neqs), niter(0), eos(eos_), vfrac_total(vfrac_tot), sie_total(sie_tot), rho(rho_), vfrac(vfrac_), sie(sie_), temp(temp_), - press(press_), Tnorm(Tguess) { + press(press_), Tnorm(Tnorm), params_(params) { jacobian = AssignIncrement(scratch, neq * neq); dx = AssignIncrement(scratch, neq); sol_scratch = AssignIncrement(scratch, 2 * neq); @@ -217,7 +226,7 @@ class PTESolverBase { // rhobar is a fixed quantity: the average density of // material m averaged over the full PTE volume rho_total = 0.0; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { 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"); @@ -230,14 +239,14 @@ class PTESolverBase { void SetVfracFromT(const Real T) { Real vsum = 0.0; // set volume fractions - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { const Real rho_min = eos[m].RhoPmin(T); const Real vmax = std::min(0.9 * robust::ratio(rhobar[m], rho_min), 1.0); vfrac[m] = (vfrac[m] > 0.0 ? std::min(vmax, vfrac[m]) : vmax); vsum += vfrac[m]; } // Normalize vfrac - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vfrac[m] *= robust::ratio(vfrac_total, vsum); } } @@ -250,22 +259,22 @@ class PTESolverBase { if (Tnorm > 0.0) { Tguess = Tnorm; } else { - Tguess = mix_params::default_tguess; - for (int m = 0; m < nmat; ++m) + Tguess = params_.default_tguess; + for (std::size_t m = 0; m < nmat; ++m) Tguess = std::max(Tguess, temp[m]); } PORTABLE_REQUIRE(Tguess > 0., "Non-positive temperature guess for PTE"); // check for sanity. basically checks that the input temperatures weren't garbage - PORTABLE_REQUIRE(Tguess < mix_params::temperature_limit, + PORTABLE_REQUIRE(Tguess < params_.temperature_limit, "Very large input temperature or temperature guess"); // iteratively increase temperature guess until all rho's are above rho_at_pmin const Real Tfactor = 10.0; bool rho_fail; - for (int i = 0; i < 3; i++) { + for (std::size_t i = 0; i < 3; i++) { SetVfracFromT(Tguess); // check to make sure the normalization didn't put us below rho_at_pmin rho_fail = false; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { const Real rho_min = eos[m].RhoPmin(Tguess); rho[m] = robust::ratio(rhobar[m], vfrac[m]); if (rho[m] < rho_min) { @@ -277,8 +286,9 @@ class PTESolverBase { if (!rho_fail) break; } - if (rho_fail) { - printf("rho < rho_min in PTE initialization! Solver may not converge."); + if (rho_fail && params_.verbose) { + PORTABLE_ALWAYS_WARN( + "rho < rho_min in PTE initialization! Solver may not converge.\n"); } return Tguess; } @@ -320,7 +330,7 @@ class PTESolverBase { // set the temperature normalization Tnorm = Tguess; - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { // scaled initial guess for temperature is just 1 temp[m] = 1.0; sie[m] = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tguess, Cache[m]); @@ -331,20 +341,22 @@ class PTESolverBase { } // note the scaling of the material internal energy densities - for (int m = 0; m < nmat; ++m) + for (std::size_t m = 0; m < nmat; ++m) u[m] = sie[m] * robust::ratio(rhobar[m], uscale); } PORTABLE_INLINE_FUNCTION Real ResidualNorm() const { Real norm = 0.0; - for (int m = 0; m < neq; m++) + for (std::size_t m = 0; m < neq; m++) norm += residual[m] * residual[m]; return 0.5 * norm; } PORTABLE_FORCEINLINE_FUNCTION - int MatIndex(const int &i, const int &j) const { return i * neq + j; } + std::size_t MatIndex(const std::size_t &i, const std::size_t &j) const { + return i * neq + j; + } // Compute the equilibrium pressure and temperature assuming an ideal EOS for each // material. Set Pideal and Tideal to the ***scaled*** solution. @@ -352,7 +364,7 @@ class PTESolverBase { void GetIdealPTE(Real &Pideal, Real &Tideal) const { Real rhoBsum = 0.0; Real Asum = 0.0; - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { Asum += vfrac[m] * robust::ratio(press[m], temp[m]); rhoBsum += rho[m] * vfrac[m] * robust::ratio(sie[m], temp[m]); } @@ -376,31 +388,31 @@ class PTESolverBase { Real *rtemp = jacobian + 3 * nmat; Real *res = jacobian + 4 * nmat; // copy out the initial guess - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { etemp[m] = u[m]; ptemp[m] = press[m]; vtemp[m] = vfrac[m]; rtemp[m] = rho[m]; } Real res_norm_old = 0.0; - for (int m = 0; m < neq; ++m) { + for (std::size_t m = 0; m < neq; ++m) { res[m] = residual[m]; res_norm_old += res[m] * res[m]; } // check if the volume fractions are reasonable const Real alpha = robust::ratio(Pideal, Tideal); - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vfrac[m] *= robust::ratio(press[m], (temp[m] * alpha)); 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--) + for (std::size_t n = m; n >= 0; n--) vfrac[n] = vtemp[n]; return; } } // fill in the rest of the state - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { rho[m] = robust::ratio(rhobar[m], vfrac[m]); const Real sie_m = @@ -414,23 +426,23 @@ class PTESolverBase { // fill in the residual solver->Residual(); Real res_norm_new = 0.0; - for (int m = 0; m < neq; m++) + for (std::size_t m = 0; m < neq; m++) res_norm_new += residual[m] * residual[m]; if (res_norm_new > res_norm_old) { // this didn't work out so reset everything - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m]; rho[m] = rtemp[m]; u[m] = etemp[m]; press[m] = ptemp[m]; } - for (int m = 0; m < neq; ++m) { + for (std::size_t m = 0; m < neq; ++m) { residual[m] = res[m]; } } else { // did work, fill in temp and energy density - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { temp[m] = Tideal; sie[m] = uscale * robust::ratio(u[m], rhobar[m]); } @@ -438,14 +450,15 @@ class PTESolverBase { } PORTABLE_INLINE_FUNCTION - Real *AssignIncrement(Real *&scratch, const int size) const { + Real *AssignIncrement(Real *&scratch, const std::size_t size) const { Real *p = scratch; scratch += size; return p; } - const int nmat, neq; - int niter; + const MixParams params_; + const std::size_t nmat, neq; + std::size_t niter; const Real vfrac_total, sie_total; const EOSIndexer &eos; const RealIndexer ρ @@ -462,28 +475,31 @@ class PTESolverBase { template PORTABLE_INLINE_FUNCTION Real ApproxTemperatureFromRhoMatU( - const int nmat, EOSIndexer &&eos, const Real u_tot, RealIndexer &&rho, + const std::size_t nmat, EOSIndexer &&eos, const Real u_tot, RealIndexer &&rho, RealIndexer &&vfrac, const Real Tguess = 0.0) { + // should these be passed in? + constexpr Real minimum_temperature = 1.e-9; + constexpr Real maximum_temperature = 1.e9; + // given material microphysical densities, volume fractions, and a total internal energy // density (rho e => erg/cm^3), solve for the temperature that gives the right sum // of material energies. this should only be used for a rough guess since it has a // hard coded and fairly loose tolerance auto ufunc = [&](const Real T) { Real usum = 0.0; - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { usum += rho[m] * vfrac[m] * eos[m].InternalEnergyFromDensityTemperature(rho[m], T); } return usum; }; - Real ulo = ufunc(mix_params::minimum_temperature); - if (u_tot < ulo) return mix_params::minimum_temperature; - Real uhi = ufunc(mix_params::maximum_temperature); - if (u_tot > uhi) return mix_params::maximum_temperature; - Real lTlo = FastMath::lg(mix_params::minimum_temperature); - Real lThi = FastMath::lg(mix_params::maximum_temperature); - if (Tguess > mix_params::minimum_temperature && - Tguess < mix_params::maximum_temperature) { + Real ulo = ufunc(minimum_temperature); + if (u_tot < ulo) return minimum_temperature; + Real uhi = ufunc(maximum_temperature); + if (u_tot > uhi) return maximum_temperature; + Real lTlo = FastMath::lg(minimum_temperature); + Real lThi = FastMath::lg(maximum_temperature); + if (Tguess > minimum_temperature && Tguess < maximum_temperature) { const Real ug = ufunc(Tguess); if (ug < u_tot) { lTlo = FastMath::lg(Tguess); @@ -493,8 +509,8 @@ PORTABLE_INLINE_FUNCTION Real ApproxTemperatureFromRhoMatU( uhi = ug; } } - int iter = 0; - constexpr int max_iter = 10; + std::size_t iter = 0; + constexpr std::size_t max_iter = 10; while (lThi - lTlo > 0.01 && iter < max_iter) { // apply bisection which is much better behaved // for materials that have a flat sie at low temperatures @@ -514,14 +530,14 @@ PORTABLE_INLINE_FUNCTION Real ApproxTemperatureFromRhoMatU( return FastMath::pow2((1.0 - alpha) * lTlo + alpha * lThi); } -inline int PTESolverRhoTRequiredScratch(const int nmat) { - int neq = nmat + 1; +inline int PTESolverRhoTRequiredScratch(const std::size_t nmat) { + std::size_t neq = nmat + 1; return neq * neq // jacobian + 4 * neq // dx, residual, and sol_scratch + 6 * nmat // all the nmat sized arrays + MAX_NUM_LAMBDAS * nmat; // the cache } -inline size_t PTESolverRhoTRequiredScratchInBytes(const int nmat) { +inline size_t PTESolverRhoTRequiredScratchInBytes(const std::size_t nmat) { return PTESolverRhoTRequiredScratch(nmat) * sizeof(Real); } @@ -552,17 +568,19 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { using mix_impl::PTESolverBase::rhobar; using mix_impl::PTESolverBase::Cache; using mix_impl::PTESolverBase::Tnorm; + using mix_impl::PTESolverBase::params_; public: // template the ctor to get type deduction/universal references prior to c++17 template PORTABLE_INLINE_FUNCTION - PTESolverRhoT(const int nmat, EOS_t &&eos, const Real vfrac_tot, const Real sie_tot, - Real_t &&rho, Real_t &&vfrac, Real_t &&sie, Real_t &&temp, Real_t &&press, - Lambda_t &&lambda, Real *scratch, const Real Tguess = 0.0) + PTESolverRhoT(const std::size_t nmat, EOS_t &&eos, const Real vfrac_tot, + const Real sie_tot, Real_t &&rho, Real_t &&vfrac, Real_t &&sie, + Real_t &&temp, Real_t &&press, Lambda_t &&lambda, Real *scratch, + const Real Tnorm = 0.0, const MixParams ¶ms = MixParams()) : mix_impl::PTESolverBase(nmat, nmat + 1, eos, vfrac_tot, sie_tot, rho, vfrac, sie, temp, - press, scratch, Tguess) { + press, scratch, Tnorm, params) { dpdv = AssignIncrement(scratch, nmat); dedv = AssignIncrement(scratch, nmat); dpdT = AssignIncrement(scratch, nmat); @@ -589,44 +607,43 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { void Residual() const { Real vsum = 0.0; Real esum = 0.0; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vsum += vfrac[m]; esum += u[m]; } residual[0] = vfrac_total - vsum; residual[1] = utotal_scale - esum; - for (int m = 0; m < nmat - 1; ++m) { + for (std::size_t m = 0; m < nmat - 1; ++m) { residual[2 + m] = press[m + 1] - press[m]; } } PORTABLE_INLINE_FUNCTION bool CheckPTE() const { - using namespace mix_params; Real mean_p = vfrac[0] * press[0]; Real error_p = 0.0; - for (int m = 1; m < nmat; ++m) { + for (std::size_t m = 1; m < nmat; ++m) { mean_p += vfrac[m] * press[m]; error_p += residual[m + 1] * residual[m + 1]; } error_p = std::sqrt(error_p); Real error_u = std::abs(residual[1]); // Check for convergence - bool converged_p = (error_p < pte_rel_tolerance_p * std::abs(mean_p) || - error_p < pte_abs_tolerance_p); - bool converged_u = (error_u < pte_rel_tolerance_e || error_u < pte_abs_tolerance_e); + bool converged_p = (error_p < params_.pte_rel_tolerance_p * std::abs(mean_p) || + error_p < params_.pte_abs_tolerance_p); + bool converged_u = + (error_u < params_.pte_rel_tolerance_e || error_u < params_.pte_abs_tolerance_e); return converged_p && converged_u; } PORTABLE_INLINE_FUNCTION void Jacobian() const { - using namespace mix_params; Real dedT_sum = 0.0; - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { ////////////////////////////// // perturb volume fractions ////////////////////////////// - Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; + Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * params_.derivative_eps; const Real vf_pert = vfrac[m] + dv; const Real rho_pert = robust::ratio(rhobar[m], vf_pert); @@ -642,7 +659,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { ////////////////////////////// // perturb temperature ////////////////////////////// - Real dT = Tequil * derivative_eps; + Real dT = Tequil * params_.derivative_eps; e_pert = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tnorm * (Tequil + dT), Cache[m]); p_pert = robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], @@ -654,15 +671,15 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { } // Fill in the Jacobian - for (int i = 0; i < neq * neq; ++i) + for (std::size_t i = 0; i < neq * neq; ++i) jacobian[i] = 0.0; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { jacobian[m] = 1.0; jacobian[neq + m] = dedv[m]; } jacobian[neq + nmat] = dedT_sum; - for (int m = 0; m < nmat - 1; m++) { - const int ind = MatIndex(2 + m, m); + for (std::size_t m = 0; m < nmat - 1; m++) { + const std::size_t ind = MatIndex(2 + m, m); jacobian[ind] = dpdv[m]; jacobian[ind + 1] = -dpdv[m + 1]; jacobian[MatIndex(2 + m, nmat)] = dpdT[m] - dpdT[m + 1]; @@ -671,18 +688,17 @@ 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) { - if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); + for (std::size_t m = 0; m < nmat; ++m) { + if (scale * dx[m] < -params_.vfrac_safety_fac * vfrac[m]) { + scale = -params_.vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); } } const Real Tnew = Tequil + scale * dx[nmat]; // control how big of a step toward rho = rho(Pmin) is allowed - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { 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); @@ -702,7 +718,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { scale = robust::ratio(-0.95 * Tequil, dx[nmat]); } // Now apply the overall scaling - for (int i = 0; i < neq; ++i) + for (std::size_t i = 0; i < neq; ++i) dx[i] *= scale; return scale; } @@ -715,11 +731,11 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { // Store the current state in temp variables for first iteration of line // search Ttemp = Tequil; - for (int m = 0; m < nmat; ++m) + for (std::size_t m = 0; m < nmat; ++m) vtemp[m] = vfrac[m]; } Tequil = Ttemp + scale * dx[nmat]; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = rhobar[m] * eos[m].InternalEnergyFromDensityTemperature( @@ -742,15 +758,15 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { }; // fixed temperature solver -inline int PTESolverFixedTRequiredScratch(const int nmat) { - int neq = nmat; +inline std::size_t PTESolverFixedTRequiredScratch(const std::size_t nmat) { + std::size_t neq = nmat; return neq * neq // jacobian + 4 * neq // dx, residual, and sol_scratch + 2 * nmat // rhobar and u in base + 2 * nmat // nmat sized arrays in fixed T solver + MAX_NUM_LAMBDAS * nmat; // the cache } -inline size_t PTESolverFixedTRequiredScratchInBytes(const int nmat) { +inline size_t PTESolverFixedTRequiredScratchInBytes(const std::size_t nmat) { return PTESolverFixedTRequiredScratch(nmat) * sizeof(Real); } @@ -777,18 +793,20 @@ class PTESolverFixedT : public mix_impl::PTESolverBase using mix_impl::PTESolverBase::rhobar; using mix_impl::PTESolverBase::Cache; using mix_impl::PTESolverBase::Tnorm; + using mix_impl::PTESolverBase::params_; public: // template the ctor to get type deduction/universal references prior to c++17 // allow the type of the temperature array to be different, potentially a const Real* template PORTABLE_INLINE_FUNCTION - PTESolverFixedT(const int nmat, EOS_t &&eos, const Real vfrac_tot, const Real T_true, - Real_t &&rho, Real_t &&vfrac, Real_t &&sie, CReal_t &&temp, - Real_t &&press, Lambda_t &&lambda, Real *scratch) + PTESolverFixedT(const std::size_t nmat, EOS_t &&eos, const Real vfrac_tot, + const Real T_true, Real_t &&rho, Real_t &&vfrac, Real_t &&sie, + CReal_t &&temp, Real_t &&press, Lambda_t &&lambda, Real *scratch, + const MixParams ¶ms = MixParams()) : mix_impl::PTESolverBase(nmat, nmat, eos, vfrac_tot, 1.0, rho, vfrac, sie, temp, press, - scratch, T_true) { + scratch, T_true, params) { dpdv = AssignIncrement(scratch, nmat); vtemp = AssignIncrement(scratch, nmat); Tequil = T_true; @@ -810,7 +828,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase this->InitRhoBarandRho(); this->SetVfracFromT(Tequil); uscale = 0.0; - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { // volume fractions have been potentially reset to ensure densitites are // larger than rho(Pmin(Tequil)); set the physical density to reflect // this change in volume fraction @@ -820,7 +838,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase // note the scaling of pressure press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tequil, Cache[m]); } - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { press[m] = robust::ratio(press[m], uscale); u[m] = sie[m] * robust::ratio(rhobar[m], uscale); } @@ -831,41 +849,40 @@ class PTESolverFixedT : public mix_impl::PTESolverBase PORTABLE_INLINE_FUNCTION void Residual() const { Real vsum = 0.0; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vsum += vfrac[m]; } residual[0] = vfrac_total - vsum; - for (int m = 0; m < nmat - 1; ++m) { + for (std::size_t m = 0; m < nmat - 1; ++m) { residual[1 + m] = press[m] - press[m + 1]; } } PORTABLE_INLINE_FUNCTION bool CheckPTE() const { - using namespace mix_params; Real mean_p = vfrac[0] * press[0]; Real error_p = 0; - for (int m = 1; m < nmat; ++m) { + for (std::size_t m = 1; m < nmat; ++m) { mean_p += vfrac[m] * press[m]; error_p += residual[m + 1] * residual[m + 1]; } error_p = std::sqrt(error_p); Real error_v = std::abs(residual[0]); // Check for convergence - bool converged_p = (error_p < pte_rel_tolerance_p * std::abs(mean_p) || - error_p < pte_abs_tolerance_p); - bool converged_v = (error_v < pte_rel_tolerance_e || error_v < pte_abs_tolerance_e); + bool converged_p = (error_p < params_.pte_rel_tolerance_p * std::abs(mean_p) || + error_p < params_.pte_abs_tolerance_p); + bool converged_v = + (error_v < params_.pte_rel_tolerance_e || error_v < params_.pte_abs_tolerance_e); return converged_p && converged_v; } PORTABLE_INLINE_FUNCTION void Jacobian() const { - using namespace mix_params; - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { ////////////////////////////// // perturb volume fractions ////////////////////////////// - Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; + Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * params_.derivative_eps; const Real vf_pert = vfrac[m] + dv; const Real rho_pert = robust::ratio(rhobar[m], vf_pert); @@ -875,12 +892,12 @@ class PTESolverFixedT : public mix_impl::PTESolverBase } // Fill in the Jacobian - for (int i = 0; i < neq * neq; ++i) + for (std::size_t i = 0; i < neq * neq; ++i) jacobian[i] = 0.0; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { jacobian[m] = 1.0; } - for (int m = 0; m < nmat - 1; m++) { + for (std::size_t m = 0; m < nmat - 1; m++) { jacobian[MatIndex(m + 1, m)] = -dpdv[m]; jacobian[MatIndex(m + 1, m + 1)] = dpdv[m + 1]; } @@ -888,17 +905,16 @@ 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) { - if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = robust::ratio(-vfrac_safety_fac * vfrac[m], dx[m]); + for (std::size_t m = 0; m < nmat; ++m) { + if (scale * dx[m] < -params_.vfrac_safety_fac * vfrac[m]) { + scale = robust::ratio(-params_.vfrac_safety_fac * vfrac[m], dx[m]); } } // control how big of a step toward rho = rho(Pmin) is allowed - for (int m = 0; m < nmat; m++) { + for (std::size_t 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]) { @@ -913,7 +929,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase } } // Now apply the overall scaling - for (int i = 0; i < neq; ++i) + for (std::size_t i = 0; i < neq; ++i) dx[i] *= scale; return scale; } @@ -925,10 +941,10 @@ class PTESolverFixedT : public mix_impl::PTESolverBase if (cache_state) { // Store the current state in temp variables for first iteration of line // search - for (int m = 0; m < nmat; ++m) + for (std::size_t m = 0; m < nmat; ++m) vtemp[m] = vfrac[m]; } - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = rhobar[m] * @@ -948,15 +964,15 @@ class PTESolverFixedT : public mix_impl::PTESolverBase }; // fixed P solver -inline int PTESolverFixedPRequiredScratch(const int nmat) { - int neq = nmat + 1; +inline std::size_t PTESolverFixedPRequiredScratch(const std::size_t nmat) { + std::size_t neq = nmat + 1; return neq * neq // jacobian + 4 * neq // dx, residual, and sol_scratch + 2 * nmat // all the nmat sized arrays in base + 3 * nmat // all the nmat sized arrays in fixedP + MAX_NUM_LAMBDAS * nmat; // the cache } -inline size_t PTESolverFixedPRequiredScratchInBytes(const int nmat) { +inline size_t PTESolverFixedPRequiredScratchInBytes(const std::size_t nmat) { return PTESolverFixedPRequiredScratch(nmat) * sizeof(Real); } @@ -986,17 +1002,19 @@ class PTESolverFixedP : public mix_impl::PTESolverBase using mix_impl::PTESolverBase::rhobar; using mix_impl::PTESolverBase::Cache; using mix_impl::PTESolverBase::Tnorm; + using mix_impl::PTESolverBase::params_; public: // template the ctor to get type deduction/universal references prior to c++17 template PORTABLE_INLINE_FUNCTION - PTESolverFixedP(const int nmat, EOS_t &&eos, const Real vfrac_tot, const Real P, + PTESolverFixedP(const std::size_t nmat, EOS_t &&eos, const Real vfrac_tot, const Real P, Real_t &&rho, Real_t &&vfrac, Real_t &&sie, Real_t &&temp, - CReal_t &&press, Lambda_t &&lambda, Real *scratch) + CReal_t &&press, Lambda_t &&lambda, Real *scratch, + const MixParams ¶ms = MixParams()) : mix_impl::PTESolverBase(nmat, nmat + 1, eos, vfrac_tot, 1.0, rho, vfrac, sie, temp, - press, scratch, 0.0) { + press, scratch, 0.0, params) { dpdv = AssignIncrement(scratch, nmat); dpdT = AssignIncrement(scratch, nmat); vtemp = AssignIncrement(scratch, nmat); @@ -1021,7 +1039,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase Tnorm = Tguess; // calculate u normalization as internal energy guess uscale = 0.0; - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { // scaled initial guess for temperature is just 1 temp[m] = 1.0; sie[m] = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tguess, Cache[m]); @@ -1029,7 +1047,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase } // note the scaling of the material internal energy densities - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { u[m] = robust::ratio(sie[m] * rhobar[m], uscale); press[m] = robust::ratio( eos[m].PressureFromDensityTemperature(rho[m], Tguess, Cache[m]), uscale); @@ -1045,7 +1063,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase PORTABLE_INLINE_FUNCTION void Residual() const { Real vsum = 0.0; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vsum += vfrac[m]; residual[m] = robust::ratio(Pequil, uscale) - press[m]; } @@ -1054,28 +1072,28 @@ class PTESolverFixedP : public mix_impl::PTESolverBase PORTABLE_INLINE_FUNCTION bool CheckPTE() const { - using namespace mix_params; Real error_p = 0; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { error_p += residual[m] * residual[m]; } error_p = std::sqrt(error_p); error_p *= uscale; Real error_v = std::abs(residual[neq - 1]); // Check for convergence - bool converged_p = (error_p < pte_rel_tolerance_p || error_p < pte_abs_tolerance_p); - bool converged_v = (error_v < pte_rel_tolerance_e || error_v < pte_abs_tolerance_e); + bool converged_p = + (error_p < params_.pte_rel_tolerance_p || error_p < params_.pte_abs_tolerance_p); + bool converged_v = + (error_v < params_.pte_rel_tolerance_e || error_v < params_.pte_abs_tolerance_e); return converged_p && converged_v; } PORTABLE_INLINE_FUNCTION void Jacobian() const { - using namespace mix_params; - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { ////////////////////////////// // perturb volume fractions ////////////////////////////// - Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; + Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * params_.derivative_eps; const Real vf_pert = vfrac[m] + dv; const Real rho_pert = robust::ratio(rhobar[m], vf_pert); @@ -1089,7 +1107,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase ////////////////////////////// // perturb temperature ////////////////////////////// - Real dT = Tequil * derivative_eps; + Real dT = Tequil * params_.derivative_eps; p_pert = robust::ratio(this->GetPressureFromPreferred(eos[m], rho[m], Tnorm * (Tequil + dT), e_pert, @@ -1099,31 +1117,30 @@ class PTESolverFixedP : public mix_impl::PTESolverBase } // Fill in the Jacobian - for (int i = 0; i < neq * neq; ++i) + for (std::size_t i = 0; i < neq * neq; ++i) jacobian[i] = 0.0; - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { jacobian[MatIndex(m, m)] = dpdv[m]; jacobian[MatIndex(m, nmat)] = dpdT[m]; } - for (int m = 0; m < neq - 1; ++m) { + for (std::size_t m = 0; m < neq - 1; ++m) { jacobian[MatIndex(neq - 1, m)] = 1.0; } } 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) { - if (scale * dx[m] < -vfrac_safety_fac * vfrac[m]) { - scale = -vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); + for (std::size_t m = 0; m < nmat; ++m) { + if (scale * dx[m] < -params_.vfrac_safety_fac * vfrac[m]) { + scale = -params_.vfrac_safety_fac * robust::ratio(vfrac[m], dx[m]); } } const Real Tnew = Tequil + scale * dx[nmat]; // control how big of a step toward rho = rho(Pmin) is allowed - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { 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); @@ -1143,7 +1160,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase scale = -0.95 * robust::ratio(Tequil, dx[nmat]); } // Now apply the overall scaling - for (int i = 0; i < neq; ++i) + for (std::size_t i = 0; i < neq; ++i) dx[i] *= scale; return scale; } @@ -1156,11 +1173,11 @@ class PTESolverFixedP : public mix_impl::PTESolverBase // Store the current state in temp variables for first iteration of line // search Ttemp = Tequil; - for (int m = 0; m < nmat; ++m) + for (std::size_t m = 0; m < nmat; ++m) vtemp[m] = vfrac[m]; } Tequil = Ttemp + scale * dx[nmat]; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = rhobar[m] * eos[m].InternalEnergyFromDensityTemperature( @@ -1181,14 +1198,14 @@ class PTESolverFixedP : public mix_impl::PTESolverBase Real Tequil, Ttemp, Pequil; }; -inline int PTESolverRhoURequiredScratch(const int nmat) { - int neq = 2 * nmat; +inline std::size_t PTESolverRhoURequiredScratch(const std::size_t nmat) { + std::size_t neq = 2 * nmat; return neq * neq // jacobian + 4 * neq // dx, residual, and sol_scratch + 8 * nmat // all the nmat sized arrays + MAX_NUM_LAMBDAS * nmat; // the cache } -inline size_t PTESolverRhoURequiredScratchInBytes(const int nmat) { +inline size_t PTESolverRhoURequiredScratchInBytes(const std::size_t nmat) { return PTESolverRhoURequiredScratch(nmat) * sizeof(Real); } @@ -1219,18 +1236,19 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { using mix_impl::PTESolverBase::rhobar; using mix_impl::PTESolverBase::Cache; using mix_impl::PTESolverBase::Tnorm; + using mix_impl::PTESolverBase::params_; public: // template the ctor to get type deduction/universal references prior to c++17 template - PORTABLE_INLINE_FUNCTION PTESolverRhoU(const int nmat, const EOS_t &&eos, - const Real vfrac_tot, const Real sie_tot, - Real_t &&rho, Real_t &&vfrac, Real_t &&sie, - Real_t &&temp, Real_t &&press, Lambda_t &&lambda, - Real *scratch, const Real Tguess = 0.0) + PORTABLE_INLINE_FUNCTION + PTESolverRhoU(const std::size_t nmat, const EOS_t &&eos, const Real vfrac_tot, + const Real sie_tot, Real_t &&rho, Real_t &&vfrac, Real_t &&sie, + Real_t &&temp, Real_t &&press, Lambda_t &&lambda, Real *scratch, + const Real Tnorm = 0.0, const MixParams ¶ms = MixParams()) : mix_impl::PTESolverBase(nmat, 2 * nmat, eos, vfrac_tot, sie_tot, rho, vfrac, sie, temp, - press, scratch, Tguess) { + press, scratch, Tnorm, params) { dpdv = AssignIncrement(scratch, nmat); dtdv = AssignIncrement(scratch, nmat); dpde = AssignIncrement(scratch, nmat); @@ -1255,28 +1273,27 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { void Residual() const { Real vsum = 0.0; Real esum = 0.0; - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vsum += vfrac[m]; esum += u[m]; } residual[0] = vfrac_total - vsum; residual[1] = utotal_scale - esum; - for (int m = 0; m < nmat - 1; ++m) { + for (std::size_t m = 0; m < nmat - 1; ++m) { residual[2 + m] = press[m + 1] - press[m]; } - for (int m = nmat + 1; m < neq; m++) { + for (std::size_t m = nmat + 1; m < neq; m++) { residual[m] = temp[m - nmat] - temp[m - nmat - 1]; } } PORTABLE_INLINE_FUNCTION bool CheckPTE() const { - using namespace mix_params; Real mean_p = vfrac[0] * press[0]; Real mean_t = rhobar[0] * temp[0]; Real error_p = 0.0; Real error_t = 0.0; - for (int m = 1; m < nmat; ++m) { + for (std::size_t m = 1; m < nmat; ++m) { mean_p += vfrac[m] * press[m]; mean_t += rhobar[m] * temp[m]; error_p += residual[m + 1] * residual[m + 1]; @@ -1286,21 +1303,20 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { error_p = std::sqrt(error_p); error_t = std::sqrt(error_t); // Check for convergence - bool converged_p = (error_p < pte_rel_tolerance_p * std::abs(mean_p) || - error_p < pte_abs_tolerance_p); - bool converged_t = - (error_t < pte_rel_tolerance_t * mean_t || error_t < pte_abs_tolerance_t); + bool converged_p = (error_p < params_.pte_rel_tolerance_p * std::abs(mean_p) || + error_p < params_.pte_abs_tolerance_p); + bool converged_t = (error_t < params_.pte_rel_tolerance_t * mean_t || + error_t < params_.pte_abs_tolerance_t); return (converged_p && converged_t); } PORTABLE_INLINE_FUNCTION void Jacobian() const { - using namespace mix_params; - for (int m = 0; m < nmat; m++) { + for (std::size_t m = 0; m < nmat; m++) { ////////////////////////////// // perturb volume fractions ////////////////////////////// - Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * derivative_eps; + Real dv = (vfrac[m] < 0.5 ? 1.0 : -1.0) * vfrac[m] * params_.derivative_eps; const Real vf_pert = vfrac[m] + dv; const Real rho_pert = robust::ratio(rhobar[m], vf_pert); @@ -1316,7 +1332,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { ////////////////////////////// // perturb energies ////////////////////////////// - const Real de = std::abs(u[m]) * derivative_eps; + const Real de = std::abs(u[m]) * params_.derivative_eps; Real e_pert = robust::ratio(u[m] + de, rhobar[m]); t_pert = robust::ratio( @@ -1328,14 +1344,14 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { uscale); dpde[m] = robust::ratio(p_pert - press[m], de); dtde[m] = robust::ratio(t_pert - temp[m], de); - if (std::abs(dtde[m]) < mix_params::min_dtde) { // must be on the cold curve - dtde[m] = derivative_eps; + if (std::abs(dtde[m]) < params_.min_dtde) { // must be on the cold curve + dtde[m] = params_.derivative_eps; } } - for (int i = 0; i < neq * neq; ++i) + for (std::size_t i = 0; i < neq * neq; ++i) jacobian[i] = 0.0; // TODO(JCD): clean all this up with MatIndex - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { jacobian[m] = 1.0; jacobian[2 * nmat + nmat + m] = 1.0; } @@ -1347,7 +1363,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { jacobian[nmat * 2 * nmat + 2 * nmat - 1] = -dpde[nmat - 1]; jacobian[(nmat + 1) * 2 * nmat + nmat] = dtde[0]; jacobian[(2 * nmat - 1) * 2 * nmat + 2 * nmat - 1] = -dtde[nmat - 1]; - for (int m = 1; m < nmat - 1; ++m) { + for (std::size_t m = 1; m < nmat - 1; ++m) { jacobian[(1 + m) * 2 * nmat + m] = -dpdv[m]; jacobian[(2 + m) * 2 * nmat + m] = dpdv[m]; jacobian[(nmat + m) * 2 * nmat + m] = -dtdv[m]; @@ -1361,10 +1377,9 @@ 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) { + for (std::size_t 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]); @@ -1391,7 +1406,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { } } // Now apply the overall scaling - for (int i = 0; i < neq; ++i) + for (std::size_t i = 0; i < neq; ++i) dx[i] *= scale; return scale; } @@ -1403,12 +1418,12 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { if (cache_state) { // Store the current state in temp variables for first iteration of line // search - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vtemp[m] = vfrac[m]; utemp[m] = u[m]; } } - for (int m = 0; m < nmat; ++m) { + for (std::size_t m = 0; m < nmat; ++m) { vfrac[m] = vtemp[m] + scale * dx[m]; rho[m] = robust::ratio(rhobar[m], vfrac[m]); u[m] = utemp[m] + scale * dx[nmat + m]; @@ -1429,16 +1444,24 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { }; template -PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { - using namespace mix_params; +PORTABLE_INLINE_FUNCTION SolverStatus PTESolver(System &s) { + SolverStatus status; + Real &err = status.residual; + bool &converged = status.converged; + // initialize the system, fill in residual, and get its norm - Real err = s.Init(); + err = s.Init(); - bool converged = false; - const int pte_max_iter = s.Nmat() * pte_max_iter_per_mat; - const Real residual_tol = s.Nmat() * pte_residual_tolerance; + // Pull out params + const MixParams ¶ms = s.GetParams(); + + converged = false; + const std::size_t pte_max_iter = s.Nmat() * params.pte_max_iter_per_mat; + const Real residual_tol = s.Nmat() * params.pte_residual_tolerance; auto &niter = s.Niter(); for (niter = 0; niter < pte_max_iter; ++niter) { + status.max_niter = std::max(status.max_niter, niter); + // Check for convergence converged = s.CheckPTE(); if (converged) break; @@ -1465,7 +1488,7 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { Real err_old = err; // 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) { + if (err > err_old + params.line_search_alpha * gradfdx) { // backtrack to middle of step scale = 0.5; Real err_mid = s.TestUpdate(scale); @@ -1475,11 +1498,13 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { // 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++) { + for (std::size_t line_iter = 0; line_iter < params.line_search_max_iter; + line_iter++) { err = s.TestUpdate(scale); - if (err < err_old + line_search_alpha * scale * gradfdx) break; + if (err < err_old + params.line_search_alpha * scale * gradfdx) break; // shrink the step if the error isn't reduced enough - scale *= line_search_fac; + scale *= params.line_search_fac; + status.max_line_niter = std::max(line_iter, status.max_line_niter); } } @@ -1498,7 +1523,7 @@ PORTABLE_INLINE_FUNCTION bool PTESolver(System &s) { if (!converged && err < residual_tol) converged = true; // undo any scaling that was applied internally for the solver s.Finalize(); - return converged; + return status; } } // namespace singularity diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index da04000dfc..747cbf5c57 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -31,6 +31,7 @@ add_executable( catch2_define.cpp eos_unit_test_helpers.hpp test_eos_modifiers.cpp + test_eos_modifiers_minimal.cpp test_eos_vector.cpp test_math_utils.cpp test_variadic_utils.cpp diff --git a/test/test_closure_pte.cpp b/test/test_closure_pte.cpp index 3db2d322ed..a9c734a3e4 100644 --- a/test/test_closure_pte.cpp +++ b/test/test_closure_pte.cpp @@ -125,7 +125,8 @@ bool run_PTE_from_state(const int num_pte, EOS *v_EOS, const Real spvol_bulk, PTESolverRhoT method( num_pte, v_EOS, vfrac_sum, sie_bulk, v_densities, v_vol_frac, v_sies, v_temperatures, v_pressures, lambdas, scratch); - pte_converged_d[0] = PTESolver(method); + auto status = PTESolver(method); + pte_converged_d[0] = status.converged; pte_u_out_d[0] = v_densities[0] * v_vol_frac[0] * v_sies[0]; for (int m = 1; m < num_pte; ++m) { pte_u_out_d[0] += v_densities[m] * v_vol_frac[m] * v_sies[m]; diff --git a/test/test_eos_modifiers.cpp b/test/test_eos_modifiers.cpp index e9ac142bd8..b75bfca16e 100644 --- a/test/test_eos_modifiers.cpp +++ b/test/test_eos_modifiers.cpp @@ -242,116 +242,3 @@ SCENARIO("EOS Builder and Modifiers", "[EOSBuilder][Modifiers][IdealGas]") { #endif // SINGULARITY_BUILD_CLOSURE } } - -SCENARIO("Relativistic EOS", "[EOSBuilder][RelativisticEOS][IdealGas]") { - GIVEN("Parameters for an ideal gas") { - constexpr Real Cv = 2.0; - constexpr Real gm1 = 0.5; - WHEN("We construct a relativistic IdealGas with EOSBuilder") { - constexpr Real cl = 1; - EOS eos = IdealGas(gm1, Cv); - eos = Modify(eos, cl); - THEN("The EOS has finite sound speeds") { - constexpr Real rho = 1e3; - constexpr Real sie = 1e3; - Real bmod = eos.BulkModulusFromDensityInternalEnergy(rho, sie); - Real cs2 = bmod / rho; - REQUIRE(cs2 < 1); - } - } - } -} - -SCENARIO("EOS Unit System", "[EOSBuilder][UnitSystem][IdealGas]") { - GIVEN("Parameters for an ideal gas") { - constexpr Real Cv = 2.0; - constexpr Real gm1 = 0.5; - GIVEN("Units with a thermal unit system") { - constexpr Real rho_unit = 1e1; - constexpr Real sie_unit = 1e-1; - constexpr Real temp_unit = 123; - WHEN("We construct an IdealGas with EOSBuilder") { - EOS eos = IdealGas(gm1, Cv); - eos = Modify(eos, rho_unit, sie_unit, temp_unit); - THEN("Units cancel out for an ideal gas") { - Real rho = 1e3; - Real sie = 1e3; - Real P = eos.PressureFromDensityInternalEnergy(rho, sie); - Real Ptrue = gm1 * rho * sie; - REQUIRE(std::abs(P - Ptrue) / Ptrue < 1e-3); - } - } - } - GIVEN("Units with length and time units") { - constexpr Real time_unit = 456; - constexpr Real length_unit = 1e2; - constexpr Real mass_unit = 1e6; - constexpr Real temp_unit = 789; - WHEN("We construct an IdealGas with EOSBuilder") { - EOS eos = IdealGas(gm1, Cv); - eos = Modify(eos, eos_units_init::length_time_units_init_tag, - time_unit, mass_unit, length_unit, temp_unit); - THEN("Units cancel out for an ideal gas") { - constexpr Real rho = 1e3; - constexpr Real sie = 1e3; - constexpr Real Ptrue = gm1 * rho * sie; - Real P = eos.PressureFromDensityInternalEnergy(rho, sie); - REQUIRE(std::abs(P - Ptrue) / Ptrue < 1e-3); - } - } - } - } -} - -SCENARIO("Serialization of modified EOSs preserves their properties", - "[ScaledEOS][IdealGas][Serialization]") { - GIVEN("A scaled ideal gas object") { - constexpr Real Cv = 2.0; - constexpr Real gm1 = 0.5; - constexpr Real scale = 2.0; - - constexpr Real rho_test = 1.0; - constexpr Real sie_test = 1.0; - constexpr Real temp_trivial = sie_test / (Cv); // = 1./2. - constexpr Real temp_test = sie_test / (Cv * scale); // = 1./4. - constexpr Real EPS = 10 * std::numeric_limits::epsilon(); - - ScaledEOS eos(IdealGas(gm1, Cv), scale); - REQUIRE(isClose(eos.TemperatureFromDensityInternalEnergy(rho_test, sie_test), - temp_test, EPS)); - EOS eos_scaled = eos; - - EOS eos_trivial = ScaledEOS(IdealGas(gm1, Cv), 1.0); - REQUIRE(isClose(eos_trivial.TemperatureFromDensityInternalEnergy(rho_test, sie_test), - temp_trivial, EPS)); - - THEN("The size of the object is larger than just the ideal gas by itself") { - REQUIRE(eos.SerializedSizeInBytes() > sizeof(IdealGas)); - } - - WHEN("We serialize the EOS") { - singularity::VectorSerializer serializer({eos_scaled, eos_trivial}); - auto [size, data] = serializer.Serialize(); - REQUIRE(size == serializer.SerializedSizeInBytes()); - REQUIRE(size > 0); - REQUIRE(data != nullptr); - - THEN("We can de-serialize the EOS") { - singularity::VectorSerializer deserializer; - deserializer.DeSerialize(data); - REQUIRE(deserializer.Size() == serializer.Size()); - - AND_THEN("The de-serialized EOS still evaluates properly") { - auto eos_new = deserializer.eos_objects[0]; - REQUIRE( - isClose(eos_new.TemperatureFromDensityInternalEnergy(rho_test, sie_test), - temp_test, EPS)); - } - } - - free(data); - } - - eos.Finalize(); - } -} diff --git a/test/test_eos_modifiers_minimal.cpp b/test/test_eos_modifiers_minimal.cpp new file mode 100644 index 0000000000..7cb3aefd41 --- /dev/null +++ b/test/test_eos_modifiers_minimal.cpp @@ -0,0 +1,161 @@ +//------------------------------------------------------------------------------ +// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// program was produced under U.S. Government contract 89233218CNA000001 +// for Los Alamos National Laboratory (LANL), which is operated by Triad +// National Security, LLC for the U.S. Department of Energy/National +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +#include + +namespace eos_units_init = singularity::eos_units_init; + +namespace EOSBuilder = singularity::EOSBuilder; +namespace thermalqs = singularity::thermalqs; +namespace variadic_utils = singularity::variadic_utils; + +using EOSBuilder::Modify; +using singularity::BilinearRampEOS; +using singularity::IdealGas; +using singularity::RelativisticEOS; +using singularity::ScaledEOS; +using singularity::ShiftedEOS; +using singularity::UnitSystem; +using singularity::Variant; + +SCENARIO("Relativistic EOS", "[EOSBuilder][RelativisticEOS][IdealGas]") { + using EOS = Variant>; + GIVEN("Parameters for an ideal gas") { + constexpr Real Cv = 2.0; + constexpr Real gm1 = 0.5; + WHEN("We construct a relativistic IdealGas with EOSBuilder") { + constexpr Real cl = 1; + EOS eos = IdealGas(gm1, Cv); + eos = Modify(eos, cl); + THEN("The EOS has finite sound speeds") { + constexpr Real rho = 1e3; + constexpr Real sie = 1e3; + Real bmod = eos.BulkModulusFromDensityInternalEnergy(rho, sie); + Real cs2 = bmod / rho; + REQUIRE(cs2 < 1); + } + } + } +} + +SCENARIO("EOS Unit System", "[EOSBuilder][UnitSystem][IdealGas]") { + using EOS = Variant>; + GIVEN("Parameters for an ideal gas") { + constexpr Real Cv = 2.0; + constexpr Real gm1 = 0.5; + GIVEN("Units with a thermal unit system") { + constexpr Real rho_unit = 1e1; + constexpr Real sie_unit = 1e-1; + constexpr Real temp_unit = 123; + WHEN("We construct an IdealGas with EOSBuilder") { + EOS eos = IdealGas(gm1, Cv); + eos = Modify(eos, rho_unit, sie_unit, temp_unit); + THEN("Units cancel out for an ideal gas") { + Real rho = 1e3; + Real sie = 1e3; + Real P = eos.PressureFromDensityInternalEnergy(rho, sie); + Real Ptrue = gm1 * rho * sie; + REQUIRE(std::abs(P - Ptrue) / Ptrue < 1e-3); + } + } + } + GIVEN("Units with length and time units") { + constexpr Real time_unit = 456; + constexpr Real length_unit = 1e2; + constexpr Real mass_unit = 1e6; + constexpr Real temp_unit = 789; + WHEN("We construct an IdealGas with EOSBuilder") { + EOS eos = IdealGas(gm1, Cv); + eos = Modify(eos, eos_units_init::length_time_units_init_tag, + time_unit, mass_unit, length_unit, temp_unit); + THEN("Units cancel out for an ideal gas") { + constexpr Real rho = 1e3; + constexpr Real sie = 1e3; + constexpr Real Ptrue = gm1 * rho * sie; + Real P = eos.PressureFromDensityInternalEnergy(rho, sie); + REQUIRE(std::abs(P - Ptrue) / Ptrue < 1e-3); + } + } + } + } +} + +SCENARIO("Serialization of modified EOSs preserves their properties", + "[ScaledEOS][IdealGas][Serialization]") { + using EOS = singularity::Variant, IdealGas>; + GIVEN("A scaled ideal gas object") { + constexpr Real Cv = 2.0; + constexpr Real gm1 = 0.5; + constexpr Real scale = 2.0; + + constexpr Real rho_test = 1.0; + constexpr Real sie_test = 1.0; + constexpr Real temp_trivial = sie_test / (Cv); // = 1./2. + constexpr Real temp_test = sie_test / (Cv * scale); // = 1./4. + constexpr Real EPS = 10 * std::numeric_limits::epsilon(); + + ScaledEOS eos(IdealGas(gm1, Cv), scale); + REQUIRE(isClose(eos.TemperatureFromDensityInternalEnergy(rho_test, sie_test), + temp_test, EPS)); + EOS eos_scaled = eos; + + EOS eos_trivial = ScaledEOS(IdealGas(gm1, Cv), 1.0); + REQUIRE(isClose(eos_trivial.TemperatureFromDensityInternalEnergy(rho_test, sie_test), + temp_trivial, EPS)); + + THEN("The size of the object is larger than just the ideal gas by itself") { + REQUIRE(eos.SerializedSizeInBytes() > sizeof(IdealGas)); + } + + WHEN("We serialize the EOS") { + singularity::VectorSerializer serializer({eos_scaled, eos_trivial}); + auto [size, data] = serializer.Serialize(); + REQUIRE(size == serializer.SerializedSizeInBytes()); + REQUIRE(size > 0); + REQUIRE(data != nullptr); + + THEN("We can de-serialize the EOS") { + singularity::VectorSerializer deserializer; + deserializer.DeSerialize(data); + REQUIRE(deserializer.Size() == serializer.Size()); + + AND_THEN("The de-serialized EOS still evaluates properly") { + auto eos_new = deserializer.eos_objects[0]; + REQUIRE( + isClose(eos_new.TemperatureFromDensityInternalEnergy(rho_test, sie_test), + temp_test, EPS)); + } + } + + free(data); + } + + eos.Finalize(); + } +} diff --git a/test/test_pte.cpp b/test/test_pte.cpp index 87f9e5bac1..52a8c63dd9 100644 --- a/test/test_pte.cpp +++ b/test/test_pte.cpp @@ -172,8 +172,8 @@ int main(int argc, char *argv[]) { PTESolverRhoT, decltype(lambda)>( NMAT, eos, 1.0, sie_tot, rho, vfrac, sie, temp, press, lambda, &scratch_d(t * nscratch_vars), Tguess); - bool success = PTESolver(method); - if (success) { + auto status = PTESolver(method); + if (status.converged) { nsuccess_d() += 1; } hist_d[method.Niter()] += 1; diff --git a/test/test_pte_2phase.cpp b/test/test_pte_2phase.cpp index ddb601e4d2..a50ec5532b 100644 --- a/test/test_pte_2phase.cpp +++ b/test/test_pte_2phase.cpp @@ -192,8 +192,8 @@ int main(int argc, char *argv[]) { PTESolverRhoT, decltype(lambda)>( NMAT, eos, 1.0, sie_tot, rho, vfrac, sie, temp, press, lambda, &scratch_d(t * nscratch_vars), Tguess); - bool success = PTESolver(method); - if (success) { + auto status = PTESolver(method); + if (status.converged) { nsuccess_d() += 1; } hist_d[method.Niter()] += 1; diff --git a/test/test_pte_3phase.cpp b/test/test_pte_3phase.cpp index 07addef331..83c3b2acb7 100644 --- a/test/test_pte_3phase.cpp +++ b/test/test_pte_3phase.cpp @@ -192,8 +192,8 @@ int main(int argc, char *argv[]) { PTESolverRhoT, decltype(lambda)>( NMAT, eos, 1.0, sie_tot, rho, vfrac, sie, temp, press, lambda, &scratch_d(t * nscratch_vars), Tguess); - bool success = PTESolver(method); - if (success) { + auto status = PTESolver(method); + if (status.converged) { nsuccess_d() += 1; } hist_d[method.Niter()] += 1; From 19af488264ba4425d653ce0666e6f0625aadeb49 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 15 Oct 2024 11:20:51 -0600 Subject: [PATCH 3/7] documentation --- doc/sphinx/src/using-closures.rst | 114 ++++++++++++++++++++++++------ 1 file changed, 94 insertions(+), 20 deletions(-) diff --git a/doc/sphinx/src/using-closures.rst b/doc/sphinx/src/using-closures.rst index 913926b7f8..d38c4ec72d 100644 --- a/doc/sphinx/src/using-closures.rst +++ b/doc/sphinx/src/using-closures.rst @@ -450,9 +450,11 @@ The constructor for the ``PTESolverRhoT`` is of the form .. code-block:: cpp template - PTESolverRhoT(const int nmat, EOS_t &&eos, const Real vfrac_tot, const Real sie_tot, + PORTABLE_INLINE_FUNCTION + PTESolverRhoT(const std::size_t nmat, EOS_t &&eos, const Real vfrac_tot, const Real sie_tot, Real_t &&rho, Real_t &&vfrac, Real_t &&sie, Real_t &&temp, Real_t &&press, - Lambda_t &&lambda, Real *scratch, const Real Tguess = 0); + Lambda_t &&lambda, Real *scratch, const Real Tnorm = 0.0, + const MixParams ¶ms = MixParams()) where ``nmat`` is the number of materials, ``eos`` is an indexer over equation of state objects, one per material, and ``vfrac_tot`` is a @@ -466,9 +468,64 @@ one per material. ``press`` is an indexer over pressures, one per material. ``lambda`` is an indexer over lambda arrays, one per material. ``scratch`` is a pointer to pre-allocated scratch memory, as described above. It is assumed enough scratch has been allocated. -Finally, the optional argument ``Tguess`` allows for host codes to -pass in an initial temperature guess for the solver. For more -information on initial guesses, see the section below. +The optional argument ``Tnorm`` allows for host codes to pass in a +normalization for the temperature scale. Initial guesses for density +and temperature may be passed in through the ``rho`` and ``sie`` input +parameters. + +The optional ``MixParams`` input contains a struct of runtime +parameters that may be used by the various PTE solvers. This struct +contains the following member fields, with default values: + +.. code-block:: cpp + + struct MixParams { + bool verbose = false; // verbosity + Real derivative_eps = 3.0e-6; + Real pte_rel_tolerance_p = 1.e-6; + Real pte_rel_tolerance_e = 1.e-6; + Real pte_rel_tolerance_t = 1.e-4; + Real pte_abs_tolerance_p = 0.0; + Real pte_abs_tolerance_e = 1.e-4; + Real pte_abs_tolerance_t = 0.0; + Real pte_residual_tolerance = 1.e-8; + std::size_t pte_max_iter_per_mat = 128; + Real line_search_alpha = 1.e-2; + std::size_t line_search_max_iter = 6; + Real line_search_fac = 0.5; + Real vfrac_safety_fac = 0.95; + Real temperature_limit = 1.0e15; + Real default_tguess = 300.; + Real min_dtde = 1.0e-16; + }; + +where here ``verbose`` specifies how erbose PTE output is, +``derivative_eps`` is the spacing used for finite differences +evaluations of equations of state when building a jacobian. The +``pte_rel_tolerance_p``, ``pte_rel_tolerance_e``, and +``pte_rel_tolerance_t`` variables are relative tolerances for the +error in the pressure, energy, temperature respectively. The +``pte_abs_tolerance_*`` variables are the same but are absolute, +rather than relative tolerances. ``pte_residual_tolerance`` is the +absolute tolerance for the residual. + +The maximum number of iterations the solver is allowed to take before +giving up is ``pte_max_iter_per_mat`` multiplied by the number of +materials used. ``line_search_alpha`` is used as a safety factor in +the line search. ``line_search_max_iter`` is the maximum number of +iterations the solver is allowed to take in the line +search. ``line_search_fac`` is the step size in the line +search. ``vfrac_safety_fac`` limites the relative amount the volume +fraction can take in a given iteration. ``temperature_limit`` is the +maximum temperature allowed by the solver. ``default_tguess`` is used +as an initial guess for temperature if a better guess is not passed in +or cannot be inferred. ``min_dtde`` is the minmum that temperature is +allowed to change with respect to energy when computing Jacobians. + +.. note:: + + If ``MixParams`` are not provided, the default values are used. Not + all ``MixParams`` are used by every solver. The constructor for the ``PTESolverRhoU`` has the same structure: @@ -478,7 +535,7 @@ The constructor for the ``PTESolverRhoU`` has the same structure: PTESolverRhoU(const int nmat, const EOS_t &&eos, const Real vfrac_tot, const Real sie_tot, Real_t &&rho, Real_t &&vfrac, Real_t &&sie, Real_t &&temp, Real_t &&press, Lambda_t &&lambda, Real *scratch, - const Real Tguess = 0); + const Real Tnorm = 0, const MixParams ¶ms = MixParams()); Both constructors are callable on host or device. In gerneral, densities and internal energies are the required inputs. However, all @@ -489,16 +546,31 @@ thermodynamically consistent with the equilibrium solution. Once a PTE solver has been constructed, one performs the solve with the ``PTESolver`` function, which takes a ``PTESolver`` object as -input and returns a boolean status of either success or failure. For -example: +input and returns a ``SolverStatus`` struct: .. code-block:: cpp auto method = PTESolverRhoT(NMAT, eos, 1.0, sie_tot, rho, vfrac, sie, temp, press, lambda, scratch); - bool success = PTESolver(method); + auto status = PTESolver(method); + +The status struct is of the form: + +.. code-block:: cpp -For an example of the PTE solver machinery in use, see the -``test_pte.cpp`` file in the tests directory. + struct SolverStatus { + bool converged; + std::size_t max_niter; + std::size_t max_line_niter; + Real residual; + }; + +where ``converged`` will report whether or not the solver successfully +converged, ``residual`` will report the final value of the residual, +``max_niter`` will report the total number of iterations that the +solver performed and ``max_line_niter`` will report the maximum number +of iterations within a line search that the solver performed. For an +example of the PTE solver machinery in use, see the ``test_pte.cpp`` +file in the tests directory. Initial Guesses for PTE Solvers ''''''''''''''''''''''''''''''' @@ -531,12 +603,14 @@ used to provide an initial guess. This function takes the form where ``nmat`` is the number of materials, ``eos`` is an indexer over equation of state objects, ``u_tot`` is the total material internal energy density (energy per unit volume), ``rho`` is an indexer over -material density, ``vfrac`` is an indexer over material volume fractions, -and the optional argument ``Tguess`` allows for callers to pass in a guess -that could accelerate finding a solution. This function does a 1-D root find -to find the temperature at which the material internal energies sum to the -total. The root find does not have a tight tolerance -- instead the -hard-coded tolerance was selected to balance performance with the accuracy -desired for an initial guess in a PTE solve. If a previous temperature value -is unavailable or some other process may have significantly modified the -temperature since it was last updated, this function can be quite effective. +material density, ``vfrac`` is an indexer over material volume +fractions, and the optional argument ``Tguess`` allows for callers to +pass in an initial guess that could accelerate finding a solution. +This function does a 1-D root find to find the temperature at which +the material internal energies sum to the total. The root find does +not have a tight tolerance -- instead the hard-coded tolerance was +selected to balance performance with the accuracy desired for an +initial guess in a PTE solve. If a previous temperature value is +unavailable or some other process may have significantly modified the +temperature since it was last updated, this function can be quite +effective. From 56ea7b024f633eaf2a4c8aae312a5cfa45bd1818 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 15 Oct 2024 11:52:08 -0600 Subject: [PATCH 4/7] fix calls to PTESolver in get_sg_eos --- singularity-eos/eos/get_sg_eos_rho_e.cpp | 3 ++- singularity-eos/eos/get_sg_eos_rho_p.cpp | 3 ++- singularity-eos/eos/get_sg_eos_rho_t.cpp | 3 ++- 3 files changed, 6 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/get_sg_eos_rho_e.cpp b/singularity-eos/eos/get_sg_eos_rho_e.cpp index 9518b3176b..adfa23ec81 100644 --- a/singularity-eos/eos/get_sg_eos_rho_e.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_e.cpp @@ -58,7 +58,8 @@ void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v, 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); + auto status = PTESolver(method); + pte_converged = status.converged; } else { // pure cell (nmat = 1) temp_pte(tid, 0) = eos_v(pte_idxs(tid, 0)) diff --git a/singularity-eos/eos/get_sg_eos_rho_p.cpp b/singularity-eos/eos/get_sg_eos_rho_p.cpp index f735ba83ca..32aef2772c 100644 --- a/singularity-eos/eos/get_sg_eos_rho_p.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_p.cpp @@ -58,7 +58,8 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v, 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); + auto status = PTESolver(method); + pte_converged = pte_status.converged; // calculate total sie for (int mp = 0; mp < npte; ++mp) { const int m = pte_mats(tid, mp); diff --git a/singularity-eos/eos/get_sg_eos_rho_t.cpp b/singularity-eos/eos/get_sg_eos_rho_t.cpp index cd292c811a..adedbea921 100644 --- a/singularity-eos/eos/get_sg_eos_rho_t.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_t.cpp @@ -59,7 +59,8 @@ void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v, 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); + auto status = PTESolver(method); + pte_converged = status.converged; // calculate total internal energy for (int mp = 0; mp < npte; ++mp) { const int m = pte_mats(tid, mp); From 304e78cd87d1fb3e78ffee6eb5c3aaf2df7483bb Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 15 Oct 2024 15:18:54 -0600 Subject: [PATCH 5/7] status --- singularity-eos/eos/get_sg_eos_rho_p.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/get_sg_eos_rho_p.cpp b/singularity-eos/eos/get_sg_eos_rho_p.cpp index 32aef2772c..24ab5119eb 100644 --- a/singularity-eos/eos/get_sg_eos_rho_p.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_p.cpp @@ -59,7 +59,7 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v, &vfrac_pte(tid, 0), &sie_pte(tid, 0), &temp_pte(tid, 0), &press_pte(tid, 0), cache[0], &solver_scratch(tid, 0)); auto status = PTESolver(method); - pte_converged = pte_status.converged; + pte_converged = status.converged; // calculate total sie for (int mp = 0; mp < npte; ++mp) { const int m = pte_mats(tid, mp); From fa8638b4ace55b2791b0ffd62d66399fc14308ea Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 15 Oct 2024 15:19:40 -0600 Subject: [PATCH 6/7] Update doc/sphinx/src/using-closures.rst Co-authored-by: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com> --- doc/sphinx/src/using-closures.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/sphinx/src/using-closures.rst b/doc/sphinx/src/using-closures.rst index d38c4ec72d..6900b448f5 100644 --- a/doc/sphinx/src/using-closures.rst +++ b/doc/sphinx/src/using-closures.rst @@ -470,7 +470,7 @@ material. ``scratch`` is a pointer to pre-allocated scratch memory, as described above. It is assumed enough scratch has been allocated. The optional argument ``Tnorm`` allows for host codes to pass in a normalization for the temperature scale. Initial guesses for density -and temperature may be passed in through the ``rho`` and ``sie`` input +and temperature may be passed in through the ``rho`` and ``temp`` input parameters. The optional ``MixParams`` input contains a struct of runtime From 15ae0161802e99352567797553572f3571c02af1 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 15 Oct 2024 15:20:02 -0600 Subject: [PATCH 7/7] Update doc/sphinx/src/using-closures.rst Co-authored-by: Jeff Peterson <83598606+jhp-lanl@users.noreply.github.com> --- doc/sphinx/src/using-closures.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/sphinx/src/using-closures.rst b/doc/sphinx/src/using-closures.rst index 6900b448f5..a3289d26b8 100644 --- a/doc/sphinx/src/using-closures.rst +++ b/doc/sphinx/src/using-closures.rst @@ -499,7 +499,7 @@ contains the following member fields, with default values: Real min_dtde = 1.0e-16; }; -where here ``verbose`` specifies how erbose PTE output is, +where here ``verbose`` enables verbose output in the PTE solve is, ``derivative_eps`` is the spacing used for finite differences evaluations of equations of state when building a jacobian. The ``pte_rel_tolerance_p``, ``pte_rel_tolerance_e``, and