Merge pull request #610 from hnil/pure_saturation_chopping

Pure saturation chopping

opm/models/blackoil/blackoilnewtonmethod.hh

@@ -39,10 +39,12 @@ BEGIN_PROPERTIES
 NEW_PROP_TAG(DpMaxRel);
 NEW_PROP_TAG(DsMax);
 NEW_PROP_TAG(PriVarOscilationThreshold);
+NEW_PROP_TAG(ProjectSaturations);
 
 SET_SCALAR_PROP(NewtonMethod, DpMaxRel, 0.3);
 SET_SCALAR_PROP(NewtonMethod, DsMax, 0.2);
 SET_SCALAR_PROP(NewtonMethod, PriVarOscilationThreshold, 1e-5);
+SET_BOOL_PROP(NewtonMethod, ProjectSaturations, false);
 
 END_PROPERTIES
 
@@ -74,6 +76,7 @@ public:
         priVarOscilationThreshold_ = EWOMS_GET_PARAM(TypeTag, Scalar, PriVarOscilationThreshold);
         dpMaxRel_ = EWOMS_GET_PARAM(TypeTag, Scalar, DpMaxRel);
         dsMax_ = EWOMS_GET_PARAM(TypeTag, Scalar, DsMax);
+        projectSaturations_ = EWOMS_GET_PARAM(TypeTag, bool, ProjectSaturations);
     }
 
     /*!
@@ -98,6 +101,7 @@ public:
         EWOMS_REGISTER_PARAM(TypeTag, Scalar, DsMax, "Maximum absolute change of any saturation in a single iteration");
         EWOMS_REGISTER_PARAM(TypeTag, Scalar, PriVarOscilationThreshold,
                              "The threshold value for the primary variable switching conditions after its meaning has switched to hinder oscilations");
+        EWOMS_REGISTER_PARAM(TypeTag,bool, ProjectSaturations, "Option for doing saturation projection");
     }
 
     /*!
@@ -309,6 +313,9 @@ protected:
 
         if (wasSwitched_[globalDofIdx])
             ++ numPriVarsSwitched_;
+        if(projectSaturations_){
+            nextValue.chopAndNormalizeSaturations();
+        }
 
         nextValue.checkDefined();
     }
@@ -319,6 +326,7 @@ private:
     Scalar priVarOscilationThreshold_;
     Scalar dpMaxRel_;
     Scalar dsMax_;
+    bool projectSaturations_;
 
     // keep track of cells where the primary variable meaning has changed
     // to detect and hinder oscillations

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -565,6 +565,83 @@ public:
         return false;
     }
 
+    bool chopAndNormalizeSaturations(){
+        if (primaryVarsMeaning() == OnePhase_p){
+            return false;
+        }
+        Scalar Sw = 0.0;
+        if (waterEnabled)
+            Sw = (*this)[Indices::waterSaturationIdx];
+
+        if (primaryVarsMeaning() == Sw_po_Sg) {
+            Scalar Sg = 0.0;
+            if (compositionSwitchEnabled)
+                Sg = (*this)[Indices::compositionSwitchIdx];
+            Scalar Ssol = 0.0;
+            if (enableSolvent)
+                Ssol =(*this) [Indices::solventSaturationIdx];
+            
+            Scalar So = 1.0 - Sw - Sg - Ssol;
+            Sw = std::min(std::max(Sw,0.0),1.0);
+            So = std::min(std::max(So,0.0),1.0);
+            Sg = std::min(std::max(Sg,0.0),1.0);
+            Ssol = std::min(std::max(Ssol,0.0),1.0);
+            
+            Scalar St = Sw + So + Sg + Ssol;
+            Sw = Sw/St;
+            Sg = Sg/St;
+            Ssol = Ssol/St;
+            assert(St>0.5);
+            if (waterEnabled)
+                (*this)[Indices::waterSaturationIdx]= Sw;
+            if (compositionSwitchEnabled)
+                (*this)[Indices::compositionSwitchIdx] = Sg;
+            if (enableSolvent)
+                (*this)[Indices::solventSaturationIdx] = Ssol;
+            return not(St==1);
+            
+        }else if (primaryVarsMeaning() == Sw_po_Rs) { 
+            Scalar Ssol = 0.0;
+            if (enableSolvent)
+                Ssol = (*this)[Indices::solventSaturationIdx];
+            Scalar So = 1.0 - Sw  - Ssol;
+            Sw = std::min(std::max(Sw,0.0),1.0);
+            So = std::min(std::max(So,0.0),1.0);
+            Ssol = std::min(std::max(Ssol,0.0),1.0);
+            //Sg = 0.0;
+            Scalar St = Sw + So + Ssol;
+            assert(St>0.5);
+            Sw=Sw/St;
+            Ssol=Ssol/St;
+            if (waterEnabled)
+                (*this)[Indices::waterSaturationIdx]= Sw;
+            if (enableSolvent)
+                (*this)[Indices::solventSaturationIdx] = Ssol;            
+            return not(St==1);
+        }else{
+            assert(primaryVarsMeaning() == Sw_pg_Rv);
+            assert(compositionSwitchEnabled);
+            Scalar Ssol=0.0;
+            if (enableSolvent)
+                Ssol = (*this)[Indices::solventSaturationIdx];
+            Scalar Sg = 1.0 - Sw - Ssol;
+            //So = 0.0;
+            Sw = std::min(std::max(Sw,0.0),1.0);
+            Sg = std::min(std::max(Sg,0.0),1.0);
+            Ssol = std::min(std::max(Ssol,0.0),1.0);
+            Scalar St = Sw + Sg + Ssol;
+            assert(St>0.5);
+            Sw=Sw/St;
+            Ssol=Ssol;
+            if (waterEnabled)
+                (*this)[Indices::waterSaturationIdx]= Sw;
+            if (enableSolvent)
+                (*this)[Indices::solventSaturationIdx] = Ssol;
+            return not(St==1);            
+        }
+        assert(false); 
+    }
+    
     BlackOilPrimaryVariables& operator=(const BlackOilPrimaryVariables& other) = default;
     BlackOilPrimaryVariables& operator=(Scalar value)
     {