Ensure FixedT and FixedP PTE solvers use sensible uscale (#568)

* Ensure fixedT and fixedP uscale values make sense

* P -> Pequil

* Update changelog

* Add protections gainst divide by zero for uscale

Author: jhp-lanl

CHANGELOG.md

@@ -11,6 +11,7 @@
 - [[PR561]](https://github.com/lanl/singularity-eos/pull/561) Fix logic for kokkos-kernels in spackage so that it is only required for closure models on GPU
 - [[PR563]](https://github.com/lanl/singularity-eos/pull/563) Fixed DensityFromPressureTemperature for the Carnahan-Starling EOS.
 - [[PR564]](https://github.com/lanl/singularity-eos/pull/564) Fix logic for numerical vs type indices by adding safeGet(), safeSet(), safeMustGet(), and safeMustSet() helpers
+- [[PR568]](https://github.com/lanl/singularity-eos/pull/568) Fix bug in FixedT and FixedP PTE solvers when materials have negative energies at the initial guess
 
 ### Changed (changing behavior/API/variables/...)
 

singularity-eos/closure/mixed_cell_models.hpp

@@ -1428,10 +1428,12 @@ class PTESolverFixedT
       // this change in volume fraction
       rho[m] = robust::ratio(rhobar[m], vfrac[m]);
       sie[m] = eos[m].InternalEnergyFromDensityTemperature(rho[m], Tequil, lambda[m]);
-      uscale += sie[m] * rho[m];
+      // Guess uscale from initial density guess. Make sure it's positive
+      uscale += std::abs(sie[m]) * rho[m];
       // note the scaling of pressure
       press[m] = eos[m].PressureFromDensityTemperature(rho[m], Tequil, lambda[m]);
     }
+    uscale = std::max(uscale, 1.0e-14);
     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);
@@ -1641,13 +1643,12 @@ class PTESolverFixedP
     const Real Tguess = this->GetTguess();
     // set the temperature normalization
     Tnorm = Tguess;
-    // calculate u normalization as internal energy guess
-    uscale = 0.0;
+    // Assume pressure and energy scales are similar
+    uscale = std::max(std::abs(Pequil), 1.0e-14);
     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, lambda[m]);
-      uscale += sie[m] * rho[m];
     }
 
     // note the scaling of the material internal energy densities