Catch singular explosions in saturation_PHSU_pure (#2250)

* Catch singular explosions in saturation_PHSU_pure

Sometimes, the Akasaka solver has issues. This is usually caught and
then it is retried with a new omega. However, sometimes it goes bad
because the J matrix is singular, and this is not caught because the
error is not recalculated.

This commit recalculates the error term to prevent a bad result.

* Add comment to explain error check in saturation_PHSU_pure

* Add extra update check after saturation_PHSU_pure

During saturated PHSU flash calculations, SatL and SatV states have an
imposed phase. This is good for stability, but there is a small chance
that they can both up up with a matching third variable (e.g.,
pressure) that is not actually at the saturation point. This commit
forces a final DT update without this requirement. If an actual solution
has been found, the the error term will still be small. If not, then we
throw an exception and try again.

This continues work on #2245.

* Ensure that saturated phase is specified

In saturation_PHSU_pure, we unspecify the phase of SatL and SatV to
perform a final check. However, for any future updates, these states
*must* be set with specified phase. This commit ensures that no matter
what happens (exception, etc.) the phase is always specified again.

* Revert "Ensure that saturated phase is specified"

This reverts commit c6b650b.

The commit caused potential recursive lookups and did not solve the
issue at hand.

* Update the ammonio saturated rhoV ancillary

The rhoV ancillary gave somewhat wrong results. This commit provides a
closer fit that prevents errors downstream.

dev/fluids/Ammonia.json

@@ -99,23 +99,31 @@
       "Tmax": 405.56,
       "Tmin": 195.495,
       "description": "rho'' = rhoc*exp(Tc/T*sum(n_i*theta^t_i))",
-      "max_abserror_percentage": 0.802513430755436,
+      "max_abserror_percentage": 0.312705,
       "n": [
-        -0.053296,
-        -3.4589,
-        -6.7572,
-        -17.26,
-        -43.12,
-        -115.18
+        0.0,
+        -3.6322306612603197,
+        43.53295636848281,
+        -396.7025942753399,
+        1946.6788208368914,
+        -5739.851852849688,
+        10341.987533476391,
+        -11168.484405912883,
+        6642.700525911098,
+        -1675.050456540271
       ],
       "reducing_value": 13696,
       "t": [
-        0.14,
-        0.44,
-        1.314,
-        3.225,
-        6.4,
-        14.0
+        0.0,
+        0.3333333333333333,
+        0.5,
+        0.6666666666666666,
+        0.8333333333333334,
+        1.0,
+        1.1666666666666667,
+        1.3333333333333333,
+        1.5,
+        1.6666666666666667
       ],
       "type": "rhoV",
       "using_tau_r": true

src/Backends/Helmholtz/VLERoutines.cpp

@@ -554,6 +554,52 @@ void SaturationSolvers::saturation_PHSU_pure(HelmholtzEOSMixtureBackend& HEOS, C
                                        info.c_str(), specified_value, error));
         }
     } while (error > 1e-9);
+    // Recalculate error
+    // The result has changed since the last error calculation.
+    // In rare scenarios, the final step can become unstable due to solving a singular
+    // J matrix. This final error check verifies that the solution is still good.
+    // Furthermore, the forced phase of SatL and SatV may have caused errors. We will recalculate them without this assumption.
+    SatL->unspecify_phase();
+    SatV->unspecify_phase();
+    SatL->update(DmolarT_INPUTS, rhoL, T);
+    SatV->update(DmolarT_INPUTS, rhoV, T);
+    negativer[0] = -(deltaV * (1 + deltaV * SatV->dalphar_dDelta()) - deltaL * (1 + deltaL * SatL->dalphar_dDelta()));
+    negativer[1] = -(deltaV * SatV->dalphar_dDelta() + SatV->alphar() + log(deltaV) - deltaL * SatL->dalphar_dDelta() - SatL->alphar() - log(deltaL));
+    switch (options.specified_variable) {
+        case saturation_PHSU_pure_options::IMPOSED_PL:
+            // -r_3 (equate calculated pressure and specified liquid pressure)
+            negativer[2] = -(SatL->p() / specified_value - 1);
+            break;
+        case saturation_PHSU_pure_options::IMPOSED_PV:
+            // -r_3 (equate calculated pressure and specified vapor pressure)
+            negativer[2] = -(SatV->p() / specified_value - 1);
+            break;
+        case saturation_PHSU_pure_options::IMPOSED_HL:
+            // -r_3 (equate calculated liquid enthalpy and specified liquid enthalpy)
+            negativer[2] = -(SatL->hmolar() - specified_value);
+            break;
+        case saturation_PHSU_pure_options::IMPOSED_HV:
+            // -r_3 (equate calculated vapor enthalpy and specified vapor enthalpy)
+            negativer[2] = -(SatV->hmolar() - specified_value);
+            break;
+        case saturation_PHSU_pure_options::IMPOSED_SL:
+            // -r_3 (equate calculated liquid entropy and specified liquid entropy)
+            negativer[2] = -(SatL->smolar() - specified_value);
+            break;
+        case saturation_PHSU_pure_options::IMPOSED_SV:
+            // -r_3 (equate calculated vapor entropy and specified vapor entropy)
+            negativer[2] = -(SatV->smolar() - specified_value);
+            break;
+        default:
+            throw ValueError(format("options.specified_variable to saturation_PHSU_pure [%d] is invalid", options.specified_variable));
+    }
+    error = sqrt(pow(negativer[0], 2) + pow(negativer[1], 2) + pow(negativer[2], 2));
+    // reset the phase for the next update.
+    SatL->specify_phase(iphase_liquid);
+    SatV->specify_phase(iphase_gas);
+    if (error > 1e-8){
+        throw SolutionError(format("saturation_PHSU_pure solver was good, but went bad. Current error is %Lg", error));
+    }
 }
 void SaturationSolvers::saturation_D_pure(HelmholtzEOSMixtureBackend& HEOS, CoolPropDbl rhomolar, saturation_D_pure_options& options)
 {