Merge pull request #4048 from hakonhagland/krnum

Add support for directional relative permeabilities

ebos/eclfluxmodule.hh

@@ -36,9 +36,12 @@
 #include <opm/models/utils/signum.hh>
 
 #include <opm/material/common/Valgrind.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
 
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
+#include <fmt/format.h>
 
 namespace Opm {
 
@@ -227,7 +230,6 @@ public:
 
         unsigned I = stencil.globalSpaceIndex(interiorDofIdx);
         unsigned J = stencil.globalSpaceIndex(exteriorDofIdx);
-
         Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
         Scalar faceArea = scvf.area();
         Scalar thpres = problem.thresholdPressure(I, J);
@@ -265,7 +267,7 @@ public:
                                         intQuantsEx,
                                         phaseIdx,//input
                                         interiorDofIdx,//input
-                                        exteriorDofIdx,//intput
+                                        exteriorDofIdx,//input
                                         Vin,
                                         Vex,
                                         I,
@@ -283,12 +285,18 @@ public:
             // or averaged? all fluids should see the same compaction?!
             const Evaluation& transMult = up.rockCompTransMultiplier();
 
+            const auto& materialLawManager = problem.materialLawManager();
+            FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown;
+            if (materialLawManager->hasDirectionalRelperms()) {
+                facedir = scvf.faceDirFromDirId();  // direction (X, Y, or Z) of the face
+            }
             if (upwindIsInterior)
                 volumeFlux[phaseIdx] =
-                    pressureDifferences[phaseIdx]*up.mobility(phaseIdx)*transMult*(-trans/faceArea);
+                    pressureDifferences[phaseIdx]*up.mobility(phaseIdx, facedir)*transMult*(-trans/faceArea);
             else
                 volumeFlux[phaseIdx] =
-                    pressureDifferences[phaseIdx]*(Toolbox::value(up.mobility(phaseIdx))*Toolbox::value(transMult)*(-trans/faceArea));
+                    pressureDifferences[phaseIdx]*
+                        (Toolbox::value(up.mobility(phaseIdx, facedir))*Toolbox::value(transMult)*(-trans/faceArea));
         }
     }
 

ebos/eclgenericproblem.cc

@@ -322,6 +322,15 @@ updatePlmixnum_()
     updateNum("PLMIXNUM", plmixnum_);
 }
 
+template<class GridView, class FluidSystem, class Scalar>
+void EclGenericProblem<GridView,FluidSystem,Scalar>::
+updateKrnum_()
+{
+    updateNum("KRNUMX", krnumx_);
+    updateNum("KRNUMY", krnumy_);
+    updateNum("KRNUMZ", krnumz_);
+}
+
 template<class GridView, class FluidSystem, class Scalar>
 bool EclGenericProblem<GridView,FluidSystem,Scalar>::
 vapparsActive(int episodeIdx) const

ebos/eclgenericproblem.hh

@@ -306,6 +306,7 @@ protected:
     void updateSatnum_();
     void updateMiscnum_();
     void updatePlmixnum_();
+    void updateKrnum_();
 
     const EclipseState& eclState_;
     const Schedule& schedule_;
@@ -318,6 +319,9 @@ protected:
     std::vector<unsigned short> satnum_;
     std::vector<unsigned short> miscnum_;
     std::vector<unsigned short> plmixnum_;
+    std::vector<unsigned short> krnumx_;
+    std::vector<unsigned short> krnumy_;
+    std::vector<unsigned short> krnumz_;
 
     std::vector<RockParams> rockParams_;
     std::vector<unsigned short> rockTableIdx_;

ebos/eclproblem.hh

@@ -1505,6 +1505,12 @@ public:
     std::shared_ptr<const EclMaterialLawManager> materialLawManager() const
     { return materialLawManager_; }
 
+    // TODO: See discussion in multiphasebaseproblem.hh for the reason why we need this method
+    const EclMaterialLawManager* materialLawManagerPtr() const
+    {
+        return materialLawManager_.get();
+    }
+
     /*!
      * \copydoc materialLawManager()
      */
@@ -2223,6 +2229,8 @@ private:
         // the PLMIX region numbers (polymer model)
         this->updatePlmixnum_();
 
+        // directional relative permeabilities
+        this->updateKrnum_();
         ////////////////////////////////
         // porosity
         updateReferencePorosity_();
@@ -2849,6 +2857,8 @@ private:
                     case FaceDir::ZPlus:
                         data = &massratebcZ_;
                         break;
+                    case FaceDir::Unknown:
+                        throw std::runtime_error("Unexpected unknown face direction");
                     }
 
                     const Evaluation rate = bcface.rate;
@@ -2883,6 +2893,8 @@ private:
                     case FaceDir::ZPlus:
                         data = &freebcZ_;
                         break;
+                    case FaceDir::Unknown:
+                        throw std::runtime_error("Unexpected unknown face direction");
                     }
 
                     for (int i = bcface.i1; i <= bcface.i2; ++i) {

opm/simulators/utils/PartiallySupportedFlowKeywords.cpp

@@ -196,7 +196,6 @@ partiallySupported()
          {
             "SATOPTS",
             {
-               {1,{false, allow_values<std::string> {"HYSTER"}, "SATOPTS(DIRECT): directional relative permeability assignment option not supported - value ignored"}}, // options
                {2,{false, allow_values<std::string> {"HYSTER"}, "SATOPTS(IRREVERS): reversible directional relative permeability option not supported – value ignored"}}, // IRREVERS
                {3,{false, allow_values<std::string> {"HYSTER"}, "SATOPTS(HYSTER): hysteresis directional relative permeability option not supported - value ignored"}}, // HYSTER
                {4,{false, allow_values<std::string> {"HYSTER"}, "SATOPTS(SURFTENS): capillary pressure surface tension pressure dependency option not supported – value ignored"}}, // SURFTENS