Move updateRelperms() to the problem class

Moves the updateRelperms() method in BlackOilInstensiveQuantities to the
Problem class. This is a more natural place for this method and it avoids
including EclMaterialManager into the  BlackOilInstensiveQuantities. The
DirectionalMobility struct is moved to a separate file such that it can be
include from both the Problem files and the BlackOilInstensiveQuanMove
updateRelperms() to the problem class

opm/models/blackoil/blackoilintensivequantities.hh

@@ -39,7 +39,6 @@
 #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>
@@ -119,30 +118,7 @@ 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>;
+    using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
 
 
 public:
@@ -258,7 +234,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);
+        problem.template updateRelperms<FluidState>(mobility_, dirMob_, fluidState_, globalSpaceIdx);
 
         // oil is the reference phase for pressure
         if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv || priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_pg_Rv) {
@@ -648,47 +624,6 @@ 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); }
 

opm/models/common/directionalmobility.hh

@@ -0,0 +1,67 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  Copyright 2022 Equinor ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+/*!
+ * \file
+ *
+ * \brief This file contains definitions related to directional mobilities
+ */
+#ifndef OPM_MODELS_DIRECTIONAL_MOBILITY_HH
+#define OPM_MODELS_DIRECTIONAL_MOBILITY_HH
+
+#include <opm/models/utils/propertysystem.hh>
+#include <opm/material/densead/Evaluation.hpp>
+
+#include <array>
+#include <stdexcept>
+
+namespace Opm {
+template <class TypeTag, class Evaluation>
+struct DirectionalMobility {
+    enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
+    // TODO: This (line below) did not work. I get error: ‘Evaluation’ is not a member of ‘Opm::Properties’
+    //  when compiling the tracer model (eclgenerictracermodel.cc). 
+    //  QuickFix: I am adding Evaluation as a class template parameter..
+    //using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+
+    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_;
+};
+} // namespace Opm
+#endif

opm/models/common/multiphasebaseproblem.hh

@@ -30,23 +30,20 @@
 
 #include "multiphasebaseproperties.hh"
 
+#include <opm/models/common/directionalmobility.hh>
 #include <opm/models/discretization/common/fvbaseproblem.hh>
 #include <opm/models/discretization/common/fvbaseproperties.hh>
 
 #include <opm/material/fluidmatrixinteractions/NullMaterial.hpp>
 #include <opm/material/common/Means.hpp>
+#include <opm/material/densead/Evaluation.hpp>
+
+#include <opm/utility/CopyablePtr.hpp>
 
 #include <dune/common/fvector.hh>
 #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
  *
@@ -70,6 +67,7 @@ class MultiPhaseBaseProblem
     using SolidEnergyLawParams = GetPropType<TypeTag, Properties::SolidEnergyLawParams>;
     using ThermalConductionLawParams = GetPropType<TypeTag, Properties::ThermalConductionLawParams>;
     using MaterialLawParams = typename GetPropType<TypeTag, Properties::MaterialLaw>::Params;
+    using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
 
     enum { dimWorld = GridView::dimensionworld };
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
@@ -251,12 +249,14 @@ public:
         static MaterialLawParams dummy;
         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
+    template <class FluidState>
+    void updateRelperms(
+        std::array<Evaluation,numPhases> &mobility,
+        DirectionalMobilityPtr &dirMob,
+        FluidState &fluidState,
+        unsigned globalSpaceIdx) const
     {
-        return nullptr;
+        return;
     }
 
     /*!