From 09d58b770a48d7ba9808702dcd6cca5749bc58c9 Mon Sep 17 00:00:00 2001
From: hnil <Halvor.M.Nilsen@sintef.no>
Date: Thu, 7 May 2020 20:44:56 +0200
Subject: [PATCH] Added option for projection saturation. It do not change the
 result for default setting. To be able to do this it is slightly unlogical
 for surfactant runs which do projection two times.

---
 opm/models/blackoil/blackoilnewtonmethod.hh   |  8 ++
 .../blackoil/blackoilprimaryvariables.hh      | 77 +++++++++++++++++++
 2 files changed, 85 insertions(+)

diff --git a/opm/models/blackoil/blackoilnewtonmethod.hh b/opm/models/blackoil/blackoilnewtonmethod.hh
index c9d583853..345a67d3b 100644
--- a/opm/models/blackoil/blackoilnewtonmethod.hh
+++ b/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
diff --git a/opm/models/blackoil/blackoilprimaryvariables.hh b/opm/models/blackoil/blackoilprimaryvariables.hh
index 044e17836..ca9e350e5 100644
--- a/opm/models/blackoil/blackoilprimaryvariables.hh
+++ b/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)
     {