Merge pull request #714 from hakonhagland/krnum

Add support for directional relative permeabilities

opm/models/blackoil/blackoilintensivequantities.hh

@@ -39,11 +39,19 @@
 #include "blackoilmicpmodules.hh"
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
 #include <opm/material/common/Valgrind.hpp>
+#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
+
+#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
+#include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/utility/CopyablePtr.hpp>
+
 #include <dune/common/fmatrix.hh>
 
 #include <cstring>
 #include <utility>
 
+#include <fmt/format.h>
+
 namespace Opm {
 /*!
  * \ingroup BlackOilModel
@@ -111,6 +119,32 @@ class BlackOilIntensiveQuantities
     using FluxIntensiveQuantities = typename FluxModule::FluxIntensiveQuantities;
     using DiffusionIntensiveQuantities = BlackOilDiffusionIntensiveQuantities<TypeTag, enableDiffusion>;
 
+    struct DirectionalMobility {
+        using array_type = std::array<Evaluation,numPhases>;
+        DirectionalMobility(const DirectionalMobility& other)
+            : mobilityX_{other.mobilityX_}, mobilityY_{other.mobilityY_}, mobilityZ_{other.mobilityZ_} {}
+        DirectionalMobility(const array_type& mX, const array_type& mY, const array_type& mZ)
+            : mobilityX_{mX}, mobilityY_{mY}, mobilityZ_{mZ} {}
+        DirectionalMobility() : mobilityX_{}, mobilityY_{}, mobilityZ_{} {}
+        array_type& getArray(int index) {
+            switch(index) {
+                case 0:
+                    return mobilityX_;
+                case 1:
+                    return mobilityY_;
+                case 2:
+                    return mobilityZ_;
+                default:
+                    throw std::runtime_error("Unexpected mobility array index");
+            }
+        }
+        array_type mobilityX_;
+        array_type mobilityY_;
+        array_type mobilityZ_;
+    };
+    using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility>;
+
+
 public:
     using FluidState = BlackOilFluidState<Evaluation,
                                           FluidSystem,
@@ -130,7 +164,6 @@ public:
             fluidState_.setRv(0.0);
         }
     }
-
     BlackOilIntensiveQuantities(const BlackOilIntensiveQuantities& other) = default;
 
     BlackOilIntensiveQuantities& operator=(const BlackOilIntensiveQuantities& other) = default;
@@ -144,7 +177,6 @@ public:
 
         const auto& problem = elemCtx.problem();
         const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
-
         const auto& linearizationType = problem.model().linearizer().getLinearizationType();
         unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
         Scalar RvMax = FluidSystem::enableVaporizedOil()
@@ -226,6 +258,7 @@ public:
         std::array<Evaluation, numPhases> pC;
         const auto& materialParams = problem.materialLawParams(globalSpaceIdx);
         MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
+        BlackOilIntensiveQuantities::updateRelperms(mobility_, dirMob_, fluidState_, problem, materialParams, globalSpaceIdx);
 
         // oil is the reference phase for pressure
         if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv || priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_pg_Rv) {
@@ -242,11 +275,6 @@ public:
                     fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
         }
 
-        // calculate relative permeabilities. note that we store the result into the
-        // mobility_ class attribute. the division by the phase viscosity happens later.
-        MaterialLaw::relativePermeabilities(mobility_, materialParams, fluidState_);
-        Valgrind::CheckDefined(mobility_);
-
         // update the Saturation functions for the blackoil solvent module.
         asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
 
@@ -395,20 +423,31 @@ public:
         paramCache.updateAll(fluidState_);
 
         // compute the phase densities and transform the phase permeabilities into mobilities
+        int nmobilities = 1;
+        std::vector<std::array<Evaluation,numPhases>*> mobilities = {&mobility_};
+        if (dirMob_) {
+            for (int i=0; i<3; i++) {
+                nmobilities += 1;
+                mobilities.push_back(&(dirMob_->getArray(i)));
+            }
+        }
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx))
                 continue;
-
             const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx);
             fluidState_.setInvB(phaseIdx, b);
-
             const auto& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
-            if (enableExtbo && phaseIdx == oilPhaseIdx)
+            for (int i = 0; i<nmobilities; i++) {
-              mobility_[phaseIdx] /= asImp_().oilViscosity();
+                if (enableExtbo && phaseIdx == oilPhaseIdx) {
-            else if (enableExtbo && phaseIdx == gasPhaseIdx)
+                    (*mobilities[i])[phaseIdx] /= asImp_().oilViscosity();
-              mobility_[phaseIdx] /= asImp_().gasViscosity();
+                }
-            else
+                else if (enableExtbo && phaseIdx == gasPhaseIdx) {
-              mobility_[phaseIdx] /= mu;
+                    (*mobilities[i])[phaseIdx] /= asImp_().gasViscosity();
+                }
+                else {
+                    (*mobilities[i])[phaseIdx] /= mu;
+                }
+            }
         }
         Valgrind::CheckDefined(mobility_);
 
@@ -532,6 +571,27 @@ public:
     const Evaluation& mobility(unsigned phaseIdx) const
     { return mobility_[phaseIdx]; }
 
+    const Evaluation& mobility(unsigned phaseIdx, FaceDir::DirEnum facedir) const
+    {
+        using Dir = FaceDir::DirEnum;
+        if (dirMob_) {
+            switch(facedir) {
+                case Dir::XPlus:
+                    return dirMob_->mobilityX_[phaseIdx];
+                case Dir::YPlus:
+                    return dirMob_->mobilityY_[phaseIdx];
+                case Dir::ZPlus:
+                    return dirMob_->mobilityZ_[phaseIdx];
+                default:
+                    throw std::runtime_error("Unexpected face direction");
+            }
+        }
+        else {
+            return mobility_[phaseIdx];
+        }
+
+    }
+
     /*!
      * \copydoc ImmiscibleIntensiveQuantities::porosity
      */
@@ -588,6 +648,46 @@ private:
     friend BlackOilBrineIntensiveQuantities<TypeTag>;
     friend BlackOilMICPIntensiveQuantities<TypeTag>;
 
+    template <class MaterialLawParams>
+    static void updateRelperms(
+        std::array<Evaluation,numPhases> &mobility,
+        DirectionalMobilityPtr &dirMob,
+        FluidState &fluidState,
+        const Problem& problem,
+        const MaterialLawParams& materialParams,
+        unsigned globalSpaceIdx)
+    {
+        // calculate relative permeabilities. note that we store the result into the
+        // mobility_ class attribute. the division by the phase viscosity happens later.
+        MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState);
+        Valgrind::CheckDefined(mobility);
+        const auto* materialLawManager = problem.materialLawManagerPtr();
+        if (materialLawManager && materialLawManager->hasDirectionalRelperms()) {
+            auto satnumIdx = materialLawManager->satnumRegionIdx(globalSpaceIdx);
+            using Dir = FaceDir::DirEnum;
+            constexpr int ndim = 3;
+            dirMob = std::make_unique<DirectionalMobility>();
+            Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus};
+            for (int i = 0; i<ndim; i++) {
+                auto krnumSatIdx = materialLawManager->getKrnumSatIdx(globalSpaceIdx, facedirs[i]);
+                auto& mob_array = dirMob->getArray(i);
+                if (krnumSatIdx != satnumIdx) {
+                    // This hack is also used by StandardWell_impl.hpp:getMobilityEval() to temporarily use a different
+                    // satnum index for a cell
+                    const auto& paramsCell = materialLawManager->connectionMaterialLawParams(krnumSatIdx, globalSpaceIdx);
+                    MaterialLaw::relativePermeabilities(mob_array, paramsCell, fluidState);
+                    // reset the cell's satnum index back to the original
+                    materialLawManager->connectionMaterialLawParams(satnumIdx, globalSpaceIdx);
+                }
+                else {
+                    // Copy the default (non-directional dependent) mobility
+                    for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+                        mob_array[phaseIdx] = mobility[phaseIdx];
+                    }
+                }
+            }
+        }
+    }
 
     Implementation& asImp_()
     { return *static_cast<Implementation*>(this); }
@@ -597,6 +697,23 @@ private:
     Evaluation porosity_;
     Evaluation rockCompTransMultiplier_;
     std::array<Evaluation,numPhases> mobility_;
+
+    // Instead of writing a custom copy constructor and a custom assignment operator just to handle
+    // the dirMob_ unique ptr member variable when copying BlackOilIntensiveQuantites (see for example
+    // updateIntensitiveQuantities_() in fvbaseelementcontext.hh for a copy example) we write the below
+    // custom wrapper class CopyablePtr which wraps the unique ptr and makes it copyable.
+    //
+    // The advantage of this approach is that we avoid having to call all the base class copy constructors and
+    // assignment operators explicitly (which is needed when writing the custom copy constructor and assignment
+    // operators) which could become a maintenance burden. For example, when adding a new base class (if that should
+    // be needed sometime in the future) to BlackOilIntensiveQuantites we could forget to update the copy
+    // constructor and assignment operators.
+    //
+    // We want each copy of the BlackOilIntensiveQuantites to be unique, (TODO: why?) so we have to make a copy
+    // of the unique_ptr each time we copy construct or assign to it from another BlackOilIntensiveQuantites.
+    // (On the other hand, if a copy could share the ptr with the original, a shared_ptr could be used instead and the
+    // wrapper would not be needed)
+    DirectionalMobilityPtr dirMob_;
 };
 
 } // namespace Opm

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.faceDirFromDirId();
+        }
         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/common/multiphasebaseproblem.hh

@@ -40,7 +40,13 @@
 #include <dune/common/fmatrix.hh>
 
 namespace Opm {
-
+// TODO: This hack is used to be able compile blackoilintensitivequantities.hh (see the function updateRelperms()) when
+//   the problem is not an EclProblem. For example if the problem is a ReservoirBlackOilVcfvProblem, the problem will not
+//   have a materialLawManager pointer (as the EclProblem has). Since this class MuitPhaseBaseProblem (see below) is a parent
+//   class for both those problem types, we can solve this problem by forward declaring EclMaterialLawManager<Traits> here
+//   and defining a method materialLawManagerPtr() here that returns a nullptr, but is overridden in EclProblem to
+//   return the real EclMaterialManager pointer.
+template <class TraitsT> class EclMaterialLawManager;
 /*!
  * \ingroup Discretization
  *
@@ -246,6 +252,13 @@ public:
         return dummy;
     }
 
+    // TODO: See the comment at the top of this file for the reason why we need this method
+    template <class TraitsT>
+    const ::Opm::EclMaterialLawManager<TraitsT>* materialLawManagerPtr() const
+    {
+        return nullptr;
+    }
+
     /*!
      * \brief Returns the temperature \f$\mathrm{[K]}\f$ within a control volume.
      *

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.faceDirFromDirId();
+                        }
+                        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_;
 };

opm/models/discretization/ecfv/ecfvstencil.hh

@@ -37,6 +37,8 @@
 #include <dune/geometry/type.hh>
 #include <dune/common/fvector.hh>
 #include <dune/common/version.hh>
+#include <opm/common/ErrorMacros.hpp>
+#include <opm/input/eclipse/EclipseState/Grid/FaceDir.hpp>
 
 #include <vector>
 
@@ -157,11 +159,15 @@ public:
 
             if (needNormal)
                 (*normal_) = intersection.centerUnitOuterNormal();
-
             const auto& geometry = intersection.geometry();
             if (needIntegrationPos)
                 (*integrationPos_) = geometry.center();
             area_ = geometry.volume();
+            // TODO: facedir_ = intersection.boundaryId();  // This did not compile for some reason
+            // using indexInInside() as in ecltransmissibility.cc instead
+            // Note that indexInInside() returns -1 for NNC faces, that's the reason we
+            // cannot convert directly to a FaceDir::DirEnum (see FaceDir.hpp)
+            dirId_  = intersection.indexInInside(); // Legal values: -1, 0, 1, 2, 3, 4, 5
         }
 
         /*!
@@ -205,10 +211,43 @@ public:
         Scalar area() const
         { return area_; }
 
+        /*!
+         * \brief Returns the direction of the face
+         */
+        int dirId() const
+        {
+            return dirId_;
+        }
+
+        /*!
+         * \brief Returns the direction of the face
+         */
+        FaceDir::DirEnum faceDirFromDirId() const
+        {
+            using Dir = FaceDir::DirEnum;
+            if (dirId_ == -1) {
+                OPM_THROW(std::runtime_error, "NNC faces does not have a face id");
+            }
+            switch(dirId_) {
+                case 0:
+                case 1:
+                    return Dir::XPlus;
+                case 2:
+                case 3:
+                    return Dir::YPlus;
+                case 4:
+                case 5:
+                    return Dir::ZPlus;
+                default:
+                    OPM_THROW(std::runtime_error, "Unexpected face id" << dirId_);
+            }
+        }
+
     private:
         ConditionalStorage<needIntegrationPos, GlobalPosition> integrationPos_;
         ConditionalStorage<needNormal, WorldVector> normal_;
         Scalar area_;
+        int dirId_;
 
         unsigned short exteriorIdx_;
     };