From 9c35a24d05810d5374a586da6b595eae807f5bf5 Mon Sep 17 00:00:00 2001 From: Patrick Mullen Date: Tue, 6 Aug 2024 11:30:15 -0600 Subject: [PATCH] Update changelog for internal energy scaling fix --- CHANGELOG.md | 3 ++- singularity-eos/closure/mixed_cell_models.hpp | 22 +++++++++++++------ 2 files changed, 17 insertions(+), 8 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 9340c409c0..6c4142d09c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,7 @@ ### Added (new features/APIs/variables/...) ### Fixed (Repair bugs, etc) +- [[PR401]](https://github.com/lanl/singularity-eos/pull/401) Fix for internal energy scaling in PTE closure ### Changed (changing behavior/API/variables/...) @@ -44,7 +45,7 @@ Date: 7/29/2024 ### Changed (changing behavior/API/variables/...) - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls - [[PR372]](https://github.com/lanl/singularity-eos/pull/372) Removed E0 from Davis Products EOS in favor of using the shifted EOS modifier. CHANGES API! -- [[PR#382]](https://github.com/lanl/singularity-eos/pull/382) Changed `get_sg_eos()` API to allow optionally specifying the mass fraction cutoff for materials to participate in the PTE solver +- [[PR#382]](https://github.com/lanl/singularity-eos/pull/382) Changed `get_sg_eos()` API to allow optionally specifying the mass fraction cutoff for materials to participate in the PTE solver ### Infrastructure (changes irrelevant to downstream codes) - [[PR329]](https://github.com/lanl/singularity-eos/pull/329) Move vinet tests into analytic test suite diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index f8c0335e38..ffbc105c25 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -328,8 +328,11 @@ class PTESolverBase { } // note the scaling of the material internal energy densities - for (int m = 0; m < nmat; ++m) + uscale_total = 0.0; + for (int m = 0; m < nmat; ++m) { u[m] = sie[m] * robust::ratio(rhobar[m], uscale); + uscale_total += u[m]; + } } PORTABLE_INLINE_FUNCTION @@ -452,7 +455,7 @@ class PTESolverBase { const RealIndexer &press; Real *jacobian, *dx, *sol_scratch, *residual, *u, *rhobar; CacheAccessor Cache; - Real rho_total, uscale, Tnorm; + Real rho_total, uscale, uscale_total, Tnorm; }; } // namespace mix_impl @@ -539,6 +542,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { using mix_impl::PTESolverBase::press; using mix_impl::PTESolverBase::rho_total; using mix_impl::PTESolverBase::uscale; + using mix_impl::PTESolverBase::uscale_total; using mix_impl::PTESolverBase::MatIndex; using mix_impl::PTESolverBase::TryIdealPTE; using mix_impl::PTESolverBase::jacobian; @@ -590,8 +594,7 @@ class PTESolverRhoT : public mix_impl::PTESolverBase { esum += u[m]; } residual[0] = vfrac_total - vsum; - // the 1 here is the scaled total internal energy density - residual[1] = 1.0 - esum; + residual[1] = uscale_total - esum; for (int m = 0; m < nmat - 1; ++m) { residual[2 + m] = press[m + 1] - press[m]; } @@ -766,6 +769,7 @@ class PTESolverFixedT : public mix_impl::PTESolverBase using mix_impl::PTESolverBase::press; using mix_impl::PTESolverBase::rho_total; using mix_impl::PTESolverBase::uscale; + using mix_impl::PTESolverBase::uscale_total; using mix_impl::PTESolverBase::MatIndex; using mix_impl::PTESolverBase::jacobian; using mix_impl::PTESolverBase::residual; @@ -817,9 +821,11 @@ class PTESolverFixedT : public mix_impl::PTESolverBase // note the scaling of pressure press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tequil, Cache[m]); } + uscale_total = 0.0; for (int m = 0; m < nmat; ++m) { press[m] = robust::ratio(press[m], uscale); u[m] = sie[m] * robust::ratio(rhobar[m], uscale); + uscale_total += u[m]; } Residual(); return ResidualNorm(); @@ -831,7 +837,6 @@ class PTESolverFixedT : public mix_impl::PTESolverBase for (int m = 0; m < nmat; ++m) { vsum += vfrac[m]; } - // the 1 here is the scaled volume fraction residual[0] = vfrac_total - vsum; for (int m = 0; m < nmat - 1; ++m) { residual[1 + m] = press[m] - press[m + 1]; @@ -975,6 +980,7 @@ class PTESolverFixedP : public mix_impl::PTESolverBase using mix_impl::PTESolverBase::press; using mix_impl::PTESolverBase::rho_total; using mix_impl::PTESolverBase::uscale; + using mix_impl::PTESolverBase::uscale_total; using mix_impl::PTESolverBase::MatIndex; using mix_impl::PTESolverBase::TryIdealPTE; using mix_impl::PTESolverBase::jacobian; @@ -1027,10 +1033,12 @@ class PTESolverFixedP : public mix_impl::PTESolverBase } // note the scaling of the material internal energy densities + uscale_total = 0.0; for (int 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); + uscale_total += u[m]; } Residual(); // Set the current guess for the equilibrium temperature. Note that this is already @@ -1207,6 +1215,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { using mix_impl::PTESolverBase::press; using mix_impl::PTESolverBase::rho_total; using mix_impl::PTESolverBase::uscale; + using mix_impl::PTESolverBase::uscale_total; using mix_impl::PTESolverBase::TryIdealPTE; using mix_impl::PTESolverBase::MatIndex; using mix_impl::PTESolverBase::jacobian; @@ -1257,8 +1266,7 @@ class PTESolverRhoU : public mix_impl::PTESolverBase { esum += u[m]; } residual[0] = vfrac_total - vsum; - // the 1 here is the scaled total internal energy density - residual[1] = 1.0 - esum; + residual[1] = uscale_total - esum; for (int m = 0; m < nmat - 1; ++m) { residual[2 + m] = press[m + 1] - press[m]; }