BlackOilPrimaryVariables: use if constexpr

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -84,10 +84,10 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
     enum { saltConcentrationIdx  = Indices::saltConcentrationIdx };
 
-    static const bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
+    static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
-    static const bool waterEnabled = Indices::waterEnabled;
+    static constexpr bool waterEnabled = Indices::waterEnabled;
-    static const bool gasEnabled = Indices::gasEnabled;
+    static constexpr bool gasEnabled = Indices::gasEnabled;
-    static const bool oilEnabled = Indices::oilEnabled;
+    static constexpr bool oilEnabled = Indices::oilEnabled;
 
     // phase indices from the fluid system
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
@@ -293,7 +293,7 @@ public:
         bool gasPresent = FluidSystem::phaseIsActive(gasPhaseIdx)?(fluidState.saturation(gasPhaseIdx) > 0.0):false;
         bool oilPresent = FluidSystem::phaseIsActive(oilPhaseIdx)?(fluidState.saturation(oilPhaseIdx) > 0.0):false;
         bool waterPresent = FluidSystem::phaseIsActive(waterPhaseIdx)?(fluidState.saturation(waterPhaseIdx) > 0.0):false;
-        static const Scalar thresholdWaterFilledCell = 1.0 - 1e-6;
+        constexpr const Scalar thresholdWaterFilledCell = 1.0 - 1e-6;
         bool onlyWater = FluidSystem::phaseIsActive(waterPhaseIdx)?(fluidState.saturation(waterPhaseIdx) > thresholdWaterFilledCell):false;
         const auto& saltSaturation = BlackOil::getSaltSaturation_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
         bool precipitatedSaltPresent = enableSaltPrecipitation?(saltSaturation > 0.0):false;
@@ -333,7 +333,7 @@ public:
                 primaryVarsMeaning_ = Sw_po_Sg;
         }
 
-        if (enableSaltPrecipitation){
+        if constexpr (enableSaltPrecipitation){
             if (precipitatedSaltPresent)
                 primaryVarsMeaningBrine_ = Sp;
             else
@@ -342,76 +342,76 @@ public:
 
         // assign the actual primary variables
         if (primaryVarsMeaning() == OnePhase_p) {
-            if (waterEnabled) {
+            if constexpr (waterEnabled) {
                 (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(waterPhaseIdx));
             } else {
                 throw std::logic_error("For single-phase runs, only pure water is presently allowed.");
             }
         }
         else if (primaryVarsMeaning() == Sw_po_Sg) {
-            if (waterEnabled)
+            if constexpr (waterEnabled)
                 (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
-            if (gasEnabled && waterEnabled && !oilEnabled) {
+            if constexpr (gasEnabled && waterEnabled && !oilEnabled) {
                 //-> water-gas system
                 (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(gasPhaseIdx));
             }
-            else if (oilEnabled) {
+            else if constexpr (oilEnabled) {
                 (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(oilPhaseIdx));
             }
-            if( compositionSwitchEnabled )
+            if constexpr (compositionSwitchEnabled)
                 (*this)[compositionSwitchIdx] = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
         }
         else if (primaryVarsMeaning() == Sw_po_Rs) {
             const auto& Rs = BlackOil::getRs_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
 
-            if (waterEnabled)
+            if constexpr (waterEnabled)
                 (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
             (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(oilPhaseIdx));
-            if( compositionSwitchEnabled )
+            if constexpr (compositionSwitchEnabled)
                 (*this)[compositionSwitchIdx] = Rs;
         }
         else if (primaryVarsMeaning() == Rvw_po_Sg && FluidSystem::enableVaporizedWater()) {
             const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
-            if (waterEnabled)
+            if constexpr (waterEnabled)
                 (*this)[waterSaturationIdx] = Rvw; //waterSaturationIdx becomes a switching idx
-            if (gasEnabled && waterEnabled && !oilEnabled) {
+            if constexpr (gasEnabled && waterEnabled && !oilEnabled) {
                 //-> water-gas system
                 (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(gasPhaseIdx));
             }
-            else if (oilEnabled) {
+            else if constexpr (oilEnabled) {
                 (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(oilPhaseIdx));
             }
-            if( compositionSwitchEnabled )
+            if constexpr (compositionSwitchEnabled)
                 (*this)[compositionSwitchIdx] = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
         }
          else if (primaryVarsMeaning() == Rvw_pg_Rv && FluidSystem::enableVaporizedWater()) {
             const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
             const auto& Rv = BlackOil::getRv_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
-            if (waterEnabled)
+            if constexpr (waterEnabled)
                 (*this)[waterSaturationIdx] = Rvw; //waterSaturationIdx becomes a switching idx
-            if (gasEnabled && waterEnabled && !oilEnabled) {
+            if constexpr (gasEnabled && waterEnabled && !oilEnabled) {
                 //-> water-gas system
                 (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(gasPhaseIdx));
             }
-            else if (oilEnabled) {
+            else if constexpr (oilEnabled) {
                 (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(oilPhaseIdx));
             }
-            if( compositionSwitchEnabled )
+            if constexpr (compositionSwitchEnabled)
                 (*this)[compositionSwitchIdx] = Rv;
         }
         else {
             assert(primaryVarsMeaning() == Sw_pg_Rv);
 
             const auto& Rv = BlackOil::getRv_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
-            if (waterEnabled)
+            if constexpr (waterEnabled)
                 (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
 
             (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(gasPhaseIdx));
-            if( compositionSwitchEnabled )
+            if constexpr (compositionSwitchEnabled)
                 (*this)[compositionSwitchIdx] = Rv;
         }
 
-        if (enableSaltPrecipitation) {
+        if constexpr (enableSaltPrecipitation) {
             if (primaryVarsMeaningBrine() == Sp) {
                 (*this)[saltConcentrationIdx] = FsToolbox::value(saltSaturation);
             }
@@ -453,7 +453,7 @@ public:
         if (waterEnabled && primaryVarsMeaning() != Rvw_po_Sg && primaryVarsMeaning() != Rvw_pg_Rv )
             Sw = (*this)[Indices::waterSaturationIdx];
 
-        if (enableSaltPrecipitation) {
+        if constexpr (enableSaltPrecipitation) {
             Scalar saltSolubility = BrineModule::saltSol(pvtRegionIndex());
             if (primaryVarsMeaningBrine() == Sp) {
                 Scalar saltSat = (*this)[saltConcentrationIdx];
@@ -482,10 +482,10 @@ public:
             if (Sw >= thresholdWaterFilledCell) {
 
                 // make sure water saturations does not exceed 1.0
-                if (waterEnabled)
+                if constexpr (waterEnabled)
                     (*this)[Indices::waterSaturationIdx] = 1.0;
                 // the hydrocarbon gas saturation is set to 0.0
-                if (compositionSwitchEnabled)
+                if constexpr (compositionSwitchEnabled)
                     (*this)[Indices::compositionSwitchIdx] = 0.0;
 
                 return false;
@@ -493,7 +493,7 @@ public:
 
             // phase equilibrium, i.e., both hydrocarbon phases are present.
             Scalar Sg = 0.0;
-            if (compositionSwitchEnabled)
+            if constexpr (compositionSwitchEnabled)
                 Sg = (*this)[Indices::compositionSwitchIdx];
 
             Scalar So = 1.0 - Sw - Sg - solventSaturation_();
@@ -526,6 +526,10 @@ public:
                                                                                       T,
                                                                                       pg);
 
+                setPrimaryVarsMeaning(Rvw_pg_Rv);
+                (*this)[Indices::pressureSwitchIdx] = pg;
+                (*this)[Indices::waterSaturationIdx] = RvwSat; //primary variable becomes Rvw
+                if constexpr (compositionSwitchEnabled) {
                     Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
                     Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
                     Scalar RvSat = enableExtbo ? ExtboModule::rv(pvtRegionIndex(),
@@ -536,11 +540,8 @@ public:
                                                                                         pg,
                                                                                         Scalar(0),
                                                                                         SoMax);
-                setPrimaryVarsMeaning(Rvw_pg_Rv);
-                (*this)[Indices::pressureSwitchIdx] = pg;
-                (*this)[Indices::waterSaturationIdx] = RvwSat; //primary variable becomes Rvw 
-                if (compositionSwitchEnabled)
                     (*this)[Indices::compositionSwitchIdx] = std::min(RvMax, RvSat);
+                }
 
                 return true;
             }
@@ -554,6 +555,8 @@ public:
                 // represents the oil phase pressure, so we do not need to change
                 // this. For the gas dissolution factor, we use the low-level blackoil
                 // PVT objects to calculate the mole fraction of gas saturated oil.
+                setPrimaryVarsMeaning(Sw_po_Rs);
+                if constexpr (compositionSwitchEnabled) {
                     Scalar po = (*this)[Indices::pressureSwitchIdx];
                     Scalar T = asImp_().temperature_();
                     Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
@@ -566,11 +569,9 @@ public:
                                                                                        po,
                                                                                        So2,
                                                                                        SoMax);
-
-                setPrimaryVarsMeaning(Sw_po_Rs);
-                if (compositionSwitchEnabled)
                     (*this)[Indices::compositionSwitchIdx] =
                         std::min(RsMax, RsSat);
+                }
 
                 return true;
             }
@@ -590,6 +591,9 @@ public:
 
                 // we start at the Rv value that corresponds to that of oil-saturated
                 // hydrocarbon gas
+                setPrimaryVarsMeaning(Sw_pg_Rv);
+                (*this)[Indices::pressureSwitchIdx] = pg;
+                if constexpr (compositionSwitchEnabled) {
                     Scalar T = asImp_().temperature_();
                     Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
                     Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
@@ -601,10 +605,8 @@ public:
                                                                                         pg,
                                                                                         Scalar(0),
                                                                                         SoMax);
-                setPrimaryVarsMeaning(Sw_pg_Rv);
-                (*this)[Indices::pressureSwitchIdx] = pg;
-                if (compositionSwitchEnabled)
                     (*this)[Indices::compositionSwitchIdx] = std::min(RvMax, RvSat);
+                }
 
                 return true;
             }
@@ -619,7 +621,7 @@ public:
                 // switch back to phase equilibrium mode if the oil phase vanishes (i.e.,
                 // the water-only case)
                 setPrimaryVarsMeaning(Sw_po_Sg);
-                if (waterEnabled)
+                if constexpr (waterEnabled)
                     (*this)[Indices::waterSaturationIdx] = 1.0; // water saturation
 
                 (*this)[Indices::compositionSwitchIdx] = 0.0; // hydrocarbon gas saturation
@@ -657,11 +659,10 @@ public:
             return false;
         }
         else if (primaryVarsMeaning() == Rvw_po_Sg) {
-            
             Scalar po = (*this)[Indices::pressureSwitchIdx];
             Scalar T = asImp_().temperature_();
             Scalar Sg = 0.0;
-            if (compositionSwitchEnabled)
+            if constexpr (compositionSwitchEnabled)
                 Sg = (*this)[Indices::compositionSwitchIdx];
 
             Scalar So = 1.0 - Sg - solventSaturation_();
@@ -684,6 +685,9 @@ public:
                 //oil phase dissappears
                 computeCapillaryPressures_(pC, /*So=*/ 0.0, Sg + solventSaturation_(), /*Sw=*/ 0.0, matParams);
                 pg = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
+                setPrimaryVarsMeaning(Rvw_pg_Rv);
+                (*this)[Indices::pressureSwitchIdx] = pg;
+                if constexpr (compositionSwitchEnabled) {
                     Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
                     Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
                     Scalar RvSat = enableExtbo ? ExtboModule::rv(pvtRegionIndex(),
@@ -694,10 +698,8 @@ public:
                                                                                         pg,
                                                                                         Scalar(0),
                                                                                         SoMax);
-                setPrimaryVarsMeaning(Rvw_pg_Rv);
-                (*this)[Indices::pressureSwitchIdx] = pg;
-                if (compositionSwitchEnabled)
                     (*this)[Indices::compositionSwitchIdx] = std::min(RvMax, RvSat);
+                }
 
                 return true;
             }
@@ -775,7 +777,7 @@ public:
                 Scalar po = pg + (pC[oilPhaseIdx] - pC[gasPhaseIdx]);
 
                 setPrimaryVarsMeaning(Sw_po_Sg);
-                if (waterEnabled)
+                if constexpr (waterEnabled)
                     (*this)[Indices::waterSaturationIdx] = 1.0;
 
                 (*this)[Indices::pressureSwitchIdx] = po;
@@ -848,10 +850,10 @@ public:
 
         if (primaryVarsMeaning() == Sw_po_Sg) {
             Scalar Sg = 0.0;
-            if (compositionSwitchEnabled)
+            if constexpr (compositionSwitchEnabled)
                 Sg = (*this)[Indices::compositionSwitchIdx];
             Scalar Ssol = 0.0;
-            if (enableSolvent)
+            if constexpr (enableSolvent)
                 Ssol =(*this) [Indices::solventSaturationIdx];
 
             Scalar So = 1.0 - Sw - Sg - Ssol;
@@ -865,17 +867,17 @@ public:
             Sg = Sg/St;
             Ssol = Ssol/St;
             assert(St>0.5);
-            if (waterEnabled)
+            if constexpr (waterEnabled)
                 (*this)[Indices::waterSaturationIdx]= Sw;
-            if (compositionSwitchEnabled)
+            if constexpr (compositionSwitchEnabled)
                 (*this)[Indices::compositionSwitchIdx] = Sg;
-            if (enableSolvent)
+            if constexpr (enableSolvent)
                 (*this)[Indices::solventSaturationIdx] = Ssol;
             return not(St==1);
 
         }else if (primaryVarsMeaning() == Sw_po_Rs) {
             Scalar Ssol = 0.0;
-            if (enableSolvent)
+            if constexpr (enableSolvent)
                 Ssol = (*this)[Indices::solventSaturationIdx];
             Scalar So = 1.0 - Sw  - Ssol;
             Sw = std::min(std::max(Sw,0.0),1.0);
@@ -886,16 +888,16 @@ public:
             assert(St>0.5);
             Sw=Sw/St;
             Ssol=Ssol/St;
-            if (waterEnabled)
+            if constexpr (waterEnabled)
                 (*this)[Indices::waterSaturationIdx]= Sw;
-            if (enableSolvent)
+            if constexpr (enableSolvent)
                 (*this)[Indices::solventSaturationIdx] = Ssol;
             return not(St==1);
         }else{
             assert(primaryVarsMeaning() == Sw_pg_Rv);
             assert(compositionSwitchEnabled);
             Scalar Ssol=0.0;
-            if (enableSolvent)
+            if constexpr (enableSolvent)
                 Ssol = (*this)[Indices::solventSaturationIdx];
             Scalar Sg = 1.0 - Sw - Ssol;
             //So = 0.0;
@@ -906,9 +908,9 @@ public:
             assert(St>0.5);
             Sw=Sw/St;
             Ssol=Ssol/St;
-            if (waterEnabled)
+            if constexpr (waterEnabled)
                 (*this)[Indices::waterSaturationIdx]= Sw;
-            if (enableSolvent)
+            if constexpr (enableSolvent)
                 (*this)[Indices::solventSaturationIdx] = Ssol;
             return not(St==1);
         }
@@ -953,90 +955,90 @@ private:
 
     Scalar solventSaturation_() const
     {
-        if (!enableSolvent)
+        if constexpr (enableSolvent)
-            return 0.0;
-
             return (*this)[Indices::solventSaturationIdx];
+        else
+            return 0.0;
     }
 
     Scalar zFraction_() const
     {
-        if (!enableExtbo)
+        if constexpr (enableExtbo)
-            return 0.0;
-
             return (*this)[Indices::zFractionIdx];
+        else
+            return 0.0;
     }
 
     Scalar polymerConcentration_() const
     {
-        if (!enablePolymer)
+        if constexpr (enablePolymer)
-            return 0.0;
-
             return (*this)[Indices::polymerConcentrationIdx];
+        else
+            return 0.0;
     }
 
     Scalar foamConcentration_() const
     {
-        if (!enableFoam)
+        if constexpr (enableFoam)
-            return 0.0;
-
             return (*this)[Indices::foamConcentrationIdx];
+        else
+            return 0.0;
     }
 
     Scalar saltConcentration_() const
     {
-        if (!enableBrine)
+        if constexpr (enableBrine)
-            return 0.0;
-
             return (*this)[Indices::saltConcentrationIdx];
+        else
+            return 0.0;
     }
 
     Scalar temperature_() const
     {
-        if (!enableEnergy)
+        if constexpr (enableEnergy)
-            return FluidSystem::reservoirTemperature();
-
             return (*this)[Indices::temperatureIdx];
+        else
+            return FluidSystem::reservoirTemperature();
     }
 
     Scalar microbialConcentration_() const
     {
-        if (!enableMICP)
+        if constexpr (enableMICP)
-            return 0.0;
-
             return (*this)[Indices::microbialConcentrationIdx];
+        else
+            return 0.0;
     }
 
     Scalar oxygenConcentration_() const
     {
-        if (!enableMICP)
+        if constexpr (enableMICP)
-            return 0.0;
-
             return (*this)[Indices::oxygenConcentrationIdx];
+        else
+            return 0.0;
     }
 
     Scalar ureaConcentration_() const
     {
-        if (!enableMICP)
+        if constexpr (enableMICP)
-            return 0.0;
-
             return (*this)[Indices::ureaConcentrationIdx];
+        else
+            return 0.0;
     }
 
     Scalar biofilmConcentration_() const
     {
-        if (!enableMICP)
+        if constexpr (enableMICP)
-            return 0.0;
-
             return (*this)[Indices::biofilmConcentrationIdx];
+        else
+            return 0.0;
     }
 
     Scalar calciteConcentration_() const
     {
-        if (!enableMICP)
+        if constexpr (enableMICP)
-            return 0.0;
-
             return (*this)[Indices::calciteConcentrationIdx];
+        else
+            return 0.0;
     }
 
     template <class Container>