some tidying up in compostional

opm/input/eclipse/EclipseState/Compositional/CompositionalConfig.cpp

@@ -234,14 +234,19 @@ namespace {
     omegaDefaultsForEosType(Opm::CompositionalConfig::EOSType eos)
     {
         using EOSType = Opm::CompositionalConfig::EOSType;
+        // Omega_a and Omega_b of the two-parameter cubic equations of state.
+        constexpr double pengRobinsonOmegaA = 0.457235529;
+        constexpr double pengRobinsonOmegaB = 0.077796074;
+        constexpr double redlichKwongOmegaA = 0.4274802;
+        constexpr double redlichKwongOmegaB = 0.08664035;
         switch (eos) {
             case EOSType::PR:
             case EOSType::PRCORR:
-                return {0.457235529, 0.077796074};
+                return {pengRobinsonOmegaA, pengRobinsonOmegaB};
             case EOSType::RK:
             case EOSType::SRK:
             case EOSType::ZJ:
-                return {0.4274802, 0.08664035};
+                return {redlichKwongOmegaA, redlichKwongOmegaB};
         }
         throw std::invalid_argument("Unknown EOSType for OMEGAA/OMEGAB defaults");
     }

opm/material/constraintsolvers/PTFlash.hpp

@@ -77,6 +77,18 @@ class PTFlash
 
     using EOSType = CompositionalConfig::EOSType;
 
+    // Tangent-plane distance below which a trial phase does not split off.
+    static constexpr Scalar stabilityTolerance = 1e-5;
+    // Successive substitution has converged, or has collapsed onto the
+    // trivial solution with all equilibrium ratios at one.
+    static constexpr Scalar substitutionTolerance = 1e-10;
+    static constexpr Scalar trivialSolutionTolerance = 1e-5;
+    // Rachford-Rice bisection stops on the residual or the interval width.
+    static constexpr Scalar bisectionResidualTolerance = 1e-16;
+    static constexpr Scalar bisectionWidthTolerance = 1e-10;
+    // Slope of the Wilson correlation for the initial equilibrium ratios.
+    static constexpr Scalar wilsonSlope = 5.3727;
+
 public:
     /*!
      * \brief Calculates the fluid state from the global mole fractions of the components and the phase pressures
@@ -291,7 +303,7 @@ class PTFlash
         constexpr int max_it = 10000;
 
         auto closeLmaxLmin = [](double max_v, double min_v) {
-            return Opm::abs(max_v - min_v) / 2. < 1e-10;
+            return Opm::abs(max_v - min_v) / 2. < bisectionWidthTolerance;
             // what if max_v < min_v?
         };
 
@@ -308,7 +320,7 @@ class PTFlash
             }
 
             // Check if midpoint fulfills g=0 or L - Lmin is sufficiently small
-            if (Opm::abs(gMid) < 1e-16 || closeLmaxLmin(Lmax, Lmin)){
+            if (Opm::abs(gMid) < bisectionResidualTolerance || closeLmaxLmin(Lmax, Lmin)){
                 return L;
             }
             // Else we repeat with midpoint being either Lmin og Lmax (depending on the signs).
@@ -394,14 +406,14 @@ class PTFlash
             OpmLog::debug("Stability test for vapor phase:");
         }
         checkStability_(fluid_state, isTrivialV, K0, y, S_v, z, /*isGas=*/true, eos_type, verbosity);
-        bool V_unstable = (S_v < (1.0 + 1e-5)) || isTrivialV;
+        bool V_unstable = (S_v < (1.0 + stabilityTolerance)) || isTrivialV;
 
         // Check for liquids stable phase
         if (verbosity == 3 || verbosity == 4) {
             OpmLog::debug("Stability test for liquid phase:");
         }
         checkStability_(fluid_state, isTrivialL, K1, x, S_l, z, /*isGas=*/false, eos_type, verbosity);
-        bool L_stable = (S_l < (1.0 + 1e-5)) || isTrivialL;
+        bool L_stable = (S_l < (1.0 + stabilityTolerance)) || isTrivialL;
 
         // L-stable means success in making liquid, V-unstable means no success in making vapour
         isStable = L_stable && V_unstable;
@@ -432,7 +444,7 @@ class PTFlash
         const auto& p_crit = FluidSystem::criticalPressure(compIdx);
         const auto& p = fluid_state.pressure(0); //for now assume no capillary pressure
 
-        const auto& tmp = Opm::exp(5.3727 * (1+acf) * (1-T_crit/T)) * (p_crit/p);
+        const auto& tmp = Opm::exp(wilsonSlope * (1+acf) * (1-T_crit/T)) * (p_crit/p);
         return tmp;
     }
 
@@ -554,8 +566,8 @@ class PTFlash
             }
 
             // Check convergence
-            isTrivial = (K_norm < 1e-5);
-            if (isTrivial || R_norm < 1e-10)
+            isTrivial = (K_norm < trivialSolutionTolerance);
+            if (isTrivial || R_norm < substitutionTolerance)
                 return;
             //todo: make sure that no mole fraction is smaller than 1e-8 ?
             //todo: take care of water!

opm/material/eos/CubicEOS.hpp

@@ -36,6 +36,13 @@ class CubicEOS
 
     static constexpr Scalar R = Constants<Scalar>::R;
 
+    // The fugacity coefficient is clamped so a vanishing mole fraction still
+    // yields a usable fugacity; the molar volume is floored against a root
+    // that has collapsed.
+    static constexpr Scalar maxFugacityCoefficient = 1e10;
+    static constexpr Scalar minFugacityCoefficient = 1e-10;
+    static constexpr Scalar minMolarVolume = 1e-7;
+
 public:
     template <class FluidState, class Params, class LhsEval = typename FluidState::ValueType>
     static LhsEval computeFugacityCoefficient(const FluidState& fs,
@@ -90,12 +97,12 @@ class CubicEOS
         // on one side, we want the mole fraction to be at
         // least 10^-3 if the fugacity is at the current pressure
         //
-        fugCoeff = min(1e10, fugCoeff);
+        fugCoeff = min(maxFugacityCoefficient, fugCoeff);
         //
         // on the other hand, if the mole fraction of the component is 100%, we want the
         // fugacity to be at least 10^-3 Pa
         //
-        fugCoeff = max(1e-10, fugCoeff);
+        fugCoeff = max(minFugacityCoefficient, fugCoeff);
         ///////////
 
         return fugCoeff;
@@ -143,15 +150,15 @@ class CubicEOS
             // i.e. the molar volume of gas is the largest one and the
             // molar volume of liquid is the smallest one
             if (isGasPhase)
-                Vm = max(1e-7, Z[2] * RT_p);
+                Vm = max(minMolarVolume, Z[2] * RT_p);
             else
-                Vm = max(1e-7, Z[0] * RT_p);
+                Vm = max(minMolarVolume, Z[0] * RT_p);
         }
         else if (numSol == 1) {
             // the EOS only has one intersection with the pressure,
             // for the other phase, we take the extremum of the EOS
             // with the largest distance from the intersection.
-            Vm = max(1e-7, Z[0] * RT_p);
+            Vm = max(minMolarVolume, Z[0] * RT_p);
         }
 
         Valgrind::CheckDefined(Vm);