From 536ec739503fbb7a764d85c7df5c4a90b24ee452 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= <hakon.hagland@gmail.com>
Date: Sun, 8 Jan 2023 12:33:16 +0100
Subject: [PATCH] Add support for IMBNUMX, IMBNUMY, IMBNUMZ

---
 ebos/eclgenericproblem.cc |  9 +++++++++
 ebos/eclgenericproblem.hh |  4 ++++
 ebos/eclproblem.hh        | 39 +++++++++++++++++----------------------
 3 files changed, 30 insertions(+), 22 deletions(-)

diff --git a/ebos/eclgenericproblem.cc b/ebos/eclgenericproblem.cc
index a6eabb5f0..8e9c69531 100644
--- a/ebos/eclgenericproblem.cc
+++ b/ebos/eclgenericproblem.cc
@@ -367,6 +367,15 @@ updateKrnum_()
     updateNum("KRNUMZ", krnumz_);
 }
 
+template<class GridView, class FluidSystem, class Scalar>
+void EclGenericProblem<GridView,FluidSystem,Scalar>::
+updateImbnum_()
+{
+    updateNum("IMBNUMX", imbnumx_);
+    updateNum("IMBNUMY", imbnumy_);
+    updateNum("IMBNUMZ", imbnumz_);
+}
+
 template<class GridView, class FluidSystem, class Scalar>
 bool EclGenericProblem<GridView,FluidSystem,Scalar>::
 vapparsActive(int episodeIdx) const
diff --git a/ebos/eclgenericproblem.hh b/ebos/eclgenericproblem.hh
index 15e6e870e..f9037df4f 100644
--- a/ebos/eclgenericproblem.hh
+++ b/ebos/eclgenericproblem.hh
@@ -305,6 +305,7 @@ protected:
     void updateMiscnum_();
     void updatePlmixnum_();
     void updateKrnum_();
+    void updateImbnum_();
 
     const EclipseState& eclState_;
     const Schedule& schedule_;
@@ -320,6 +321,9 @@ protected:
     std::vector<unsigned short> krnumx_;
     std::vector<unsigned short> krnumy_;
     std::vector<unsigned short> krnumz_;
+    std::vector<unsigned short> imbnumx_;
+    std::vector<unsigned short> imbnumy_;
+    std::vector<unsigned short> imbnumz_;
 
     std::vector<RockParams> rockParams_;
     std::vector<unsigned short> rockTableIdx_;
diff --git a/ebos/eclproblem.hh b/ebos/eclproblem.hh
index d4d99bd06..90fcc4568 100644
--- a/ebos/eclproblem.hh
+++ b/ebos/eclproblem.hh
@@ -1415,6 +1415,11 @@ public:
         return materialLawManager_->materialLawParams(globalDofIdx);
     }
 
+    const MaterialLawParams& materialLawParams(unsigned globalDofIdx, FaceDir::DirEnum facedir) const
+    {
+        return materialLawManager_->materialLawParams(globalDofIdx, facedir);
+    }
+
     /*!
      * \brief Return the parameters for the energy storage law of the rock
      */
@@ -1455,34 +1460,24 @@ public:
         FluidState &fluidState,
         unsigned globalSpaceIdx) const
     {
-        // calculate relative permeabilities. note that we store the result into the
-        // mobility_ class attribute. the division by the phase viscosity happens later.
-        const auto& materialParams = materialLawParams(globalSpaceIdx);
-        MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState);
-        Valgrind::CheckDefined(mobility);
-        if (materialLawManager_->hasDirectionalRelperms()) {
-            auto satnumIdx = materialLawManager_->satnumRegionIdx(globalSpaceIdx);
+        {
+            // calculate relative permeabilities. note that we store the result into the
+            // mobility_ class attribute. the division by the phase viscosity happens later.
+            const auto& materialParams = materialLawParams(globalSpaceIdx);
+            MaterialLaw::relativePermeabilities(mobility, materialParams, fluidState);
+            Valgrind::CheckDefined(mobility);
+        }
+        if (materialLawManager_->hasDirectionalRelperms()
+               || materialLawManager_->hasDirectionalImbnum())
+        {
             using Dir = FaceDir::DirEnum;
             constexpr int ndim = 3;
             dirMob = std::make_unique<DirectionalMobility<TypeTag, Evaluation>>();
             Dir facedirs[ndim] = {Dir::XPlus, Dir::YPlus, Dir::ZPlus};
             for (int i = 0; i<ndim; i++) {
-                auto krnumSatIdx = materialLawManager_->getKrnumSatIdx(globalSpaceIdx, facedirs[i]);
+                const auto& materialParams = materialLawParams(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];
-                    }
-                }
+                MaterialLaw::relativePermeabilities(mob_array, materialParams, fluidState);
             }
         }
     }