Add more support for directional relperms

Adds support for directional relperms for the tpfa linearizer.

opm/models/blackoil/blackoillocalresidualtpfa.hh

@@ -38,6 +38,8 @@
 #include "blackoildiffusionmodule.hh"
 #include "blackoilmicpmodules.hh"
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
+
 
 namespace Opm {
 /*!
@@ -198,7 +200,8 @@ public:
                             const IntensiveQuantities& intQuantsIn,
                             const IntensiveQuantities& intQuantsEx,
                             const Scalar trans,
-                            const Scalar faceArea)
+                            const Scalar faceArea,
+                            const FaceDir::DirEnum facedir)
     {
         flux = 0.0;
         Scalar Vin = problem.model().dofTotalVolume(globalIndexIn);
@@ -233,7 +236,8 @@ public:
                          distZ * g,
                          thpres,
                          trans,
-                         faceArea);
+                         faceArea,
+                         facedir);
     }
 
     // This function demonstrates compatibility with the ElementContext-based interface.
@@ -264,6 +268,11 @@ public:
         const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx);
         Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
         Scalar faceArea = scvf.area();
+        const auto& materialLawManager = problem.materialLawManager();
+        FaceDir::DirEnum facedir = FaceDir::DirEnum::Unknown; // Use an arbitrary
+        if (materialLawManager->hasDirectionalRelperms()) {
+            facedir = scvf.dirId();
+        }
         Scalar thpres = problem.thresholdPressure(globalIndexIn, globalIndexEx);
 
         // estimate the gravity correction: for performance reasons we use a simplified
@@ -295,7 +304,8 @@ public:
                          distZ * g,
                          thpres,
                          trans,
-                         faceArea);
+                         faceArea,
+                         facedir);
     }
 
     static void calculateFluxes_(RateVector& flux,
@@ -308,7 +318,8 @@ public:
                                  const Scalar& distZg,
                                  const Scalar& thpres,
                                  const Scalar& trans,
-                                 const Scalar& faceArea)
+                                 const Scalar& faceArea,
+                                 const FaceDir::DirEnum facedir)
     {
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx))
@@ -346,9 +357,10 @@ public:
                 darcyFlux = 0.0; // NB maybe we could drop calculations
             } else {
                 if (globalUpIndex == globalIndexIn)
-                    darcyFlux = pressureDifference * up.mobility(phaseIdx) * transMult * (-trans / faceArea);
+                    darcyFlux = pressureDifference * up.mobility(phaseIdx, facedir) * transMult * (-trans / faceArea);
                 else
-                    darcyFlux = pressureDifference * (Toolbox::value(up.mobility(phaseIdx)) * Toolbox::value(transMult) * (-trans / faceArea));
+                    darcyFlux = pressureDifference *
+                       (Toolbox::value(up.mobility(phaseIdx, facedir)) * Toolbox::value(transMult) * (-trans / faceArea));
             }
 
             unsigned pvtRegionIdx = up.pvtRegionIndex();

opm/models/discretization/common/tpfalinearizer.hh

@@ -35,6 +35,7 @@
 
 #include <opm/material/common/Exceptions.hpp>
 #include <opm/grid/utility/SparseTable.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
 
 #include <dune/common/version.hh>
 #include <dune/common/fvector.hh>
@@ -314,7 +315,8 @@ private:
         unsigned numCells = model.numTotalDof();
         neighborInfo_.reserve(numCells, 6 * numCells);
         std::vector<NeighborInfo> loc_nbinfo;
-
+        const auto& materialLawManager = problem_().materialLawManager();
+        using FaceDirection = FaceDir::DirEnum;
         ElementIterator elemIt = gridView_().template begin<0>();
         const ElementIterator elemEndIt = gridView_().template end<0>();
         for (; elemIt != elemEndIt; ++elemIt) {
@@ -330,8 +332,14 @@ private:
                     sparsityPattern[myIdx].insert(neighborIdx);
                     if (dofIdx > 0) {
                         const double trans = problem_().transmissibility(myIdx, neighborIdx);
-                        const double area = stencil.interiorFace(dofIdx - 1).area();
+                        const auto scvfIdx = dofIdx - 1;
-                        loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, trans, area};
+                        const auto& scvf = stencil.interiorFace(scvfIdx);
+                        const double area = scvf.area();
+                        FaceDirection dirId = FaceDirection::Unknown;
+                        if (materialLawManager->hasDirectionalRelperms()) {
+                            dirId = scvf.dirId();
+                        }
+                        loc_nbinfo[dofIdx - 1] = NeighborInfo{neighborIdx, trans, area, dirId};
                     }
                 }
                 neighborInfo_.appendRow(loc_nbinfo.begin(), loc_nbinfo.end());
@@ -410,7 +418,8 @@ private:
                 }
                 const IntensiveQuantities& intQuantsEx = *intQuantsExP;
                 LocalResidual::computeFlux(
-                    adres, problem_(), globI, globJ, intQuantsIn, intQuantsEx, nbInfo.trans, nbInfo.faceArea);
+                       adres, problem_(), globI, globJ, intQuantsIn, intQuantsEx,
+                           nbInfo.trans, nbInfo.faceArea, nbInfo.faceDirection);
                 adres *= nbInfo.faceArea;
                 setResAndJacobi(res, bMat, adres);
                 residual_[globI] += res;
@@ -467,6 +476,7 @@ private:
         unsigned int neighbor;
         double trans;
         double faceArea;
+        FaceDir::DirEnum faceDirection;
     };
     SparseTable<NeighborInfo> neighborInfo_;
 };