[mpfa] introduce TensorLambdaFactory and unique interface for lambdas

the lambdas are used to obtain the tensors related to a law, i.e. an
expression for a flux (advective, diffusive, conductive). We now introduce
a lambda factory that passes lambdas for the different physical processes.

dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh

@@ -25,6 +25,7 @@
 
 #include <dumux/implicit/properties.hh>
 #include <dumux/discretization/methods.hh>
+#include <dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh>
 
 namespace Dumux
 {
@@ -253,15 +254,12 @@ private:
         using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
         using AdvectionFiller = typename AdvectionType::CacheFiller;
 
-        static constexpr bool usesMpfa = AdvectionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
+        static constexpr auto AdvectionMethod = AdvectionType::myDiscretizationMethod;
+        using LambdaFactory = TensorLambdaFactory<TypeTag, AdvectionMethod>;
 
-        // maybe solve the local system subject to K
-        if (usesMpfa)
-        {
-            auto K = [] (const Element& element, const auto& volVars, const SubControlVolume& scv)
-                        { return volVars.permeability(); };
-            iv.solveLocalSystem(K);
-        }
+        // maybe solve the local system subject to K (if AdvectionType uses mpfa)
+        if (AdvectionMethod == DiscretizationMethods::CCMpfa)
+            iv.solveLocalSystem(LambdaFactory::getAdvectionLambda());
 
         // fill the caches of all scvfs within this interaction volume
         AdvectionFiller::fill(scvfFluxVarsCache, problem(), element(), fvGeometry(), elemVolVars(), scvFace(), *this);
@@ -308,9 +306,11 @@ private:
         using DiffusionType = typename GET_PROP_TYPE(TypeTag, MolecularDiffusionType);
         using DiffusionFiller = typename DiffusionType::CacheFiller;
 
+        static constexpr auto DiffusionMethod = DiffusionType::myDiscretizationMethod;
+        using LambdaFactory = TensorLambdaFactory<TypeTag, DiffusionMethod>;
+
         static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
         static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
-        static constexpr bool usesMpfa = DiffusionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
 
         for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
         {
@@ -320,12 +320,8 @@ private:
                     continue;
 
                 // solve the local system subject to the diffusion tensor (if uses mpfa)
-                if (usesMpfa)
-                {
-                    auto D = [phaseIdx, compIdx](const Element& element, const auto& volVars, const SubControlVolume& scv)
-                                                { return volVars.diffusionCoefficient(phaseIdx, compIdx); };
-                    iv.solveLocalSystem(D);
-                }
+                if (DiffusionMethod == DiscretizationMethods::CCMpfa)
+                    iv.solveLocalSystem(LambdaFactory::getDiffusionLambda(phaseIdx, compIdx));
 
                 // fill the caches of all scvfs within this interaction volume
                 DiffusionFiller::fill(scvfFluxVarsCache, phaseIdx, compIdx, problem(), element(), fvGeometry(), elemVolVars(), scvFace(), *this);
@@ -375,23 +371,13 @@ private:
     {
         using HeatConductionType = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
         using HeatConductionFiller = typename HeatConductionType::CacheFiller;
-        using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
 
-        static constexpr bool usesMpfa = HeatConductionType::myDiscretizationMethod == DiscretizationMethods::CCMpfa;
+        static constexpr auto HeatConductionMethod = HeatConductionType::myDiscretizationMethod;
+        using LambdaFactory = TensorLambdaFactory<TypeTag, HeatConductionMethod>;
 
-        // maybe solve the local system subject to K
-        if (usesMpfa)
-        {
-            const auto& prob = problem();
-            const auto& fvGeom = fvGeometry();
-            auto lambda = [&prob, &fvGeom](const Element& element, const auto& volVars, const SubControlVolume& scv)
-            { return ThermalConductivityModel::effectiveThermalConductivity(volVars,
-                                                                            prob.spatialParams(),
-                                                                            element,
-                                                                            fvGeom,
-                                                                            scv); };
-            iv.solveLocalSystem(lambda);
-        }
+        // maybe solve the local system subject to fourier coefficient
+        if (HeatConductionMethod == DiscretizationMethods::CCMpfa)
+            iv.solveLocalSystem(LambdaFactory::getHeatConductionLambda());
 
         // fill the caches of all scvfs within this interaction volume
         HeatConductionFiller::fill(scvfFluxVarsCache, problem(), element(), fvGeometry(), elemVolVars(), scvFace(), *this);

dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh

@@ -55,6 +55,11 @@ class InteriorBoundaryData
     //! Note that the return types are also "wrong" (Here just to satisfy the compiler)
     struct CompleteCoupledFacetData
     {
+        const Problem& problem() const
+        {
+            DUNE_THROW(Dune::NotImplemented, "No implementation of the CompleteCoupledFacetData class provided");
+        }
+
         const SpatialParams& spatialParams() const
         {
             DUNE_THROW(Dune::NotImplemented, "No implementation of the CompleteCoupledFacetData class provided");

dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh

@@ -485,7 +485,7 @@ private:
             const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
             const auto& posVolVars = elemVolVars_()[posGlobalScv];
             const auto& element = localElement_(posLocalScvIdx);
-            const auto tensor = getTensor(element, posVolVars, posGlobalScv);
+            const auto tensor = getTensor(problem_(), element, posVolVars, fvGeometry_(), posGlobalScv);
 
             // the omega factors of the "positive" sub volume
             auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor);
@@ -523,23 +523,28 @@ private:
                                 // get interior boundary data
                                 const auto& data = interiorBoundaryData_[this->findIndexInVector(interiorBoundaryScvfIndexSet_(), curLocalScvfIdx)];
 
-                                // get the volvars on the actual interior boundary face
-                                const auto facetVolVars = data.facetVolVars(fvGeometry_());
+                                // obtain the complete data on the facet element
+                                const auto completeFacetData = data.completeCoupledFacetData();
 
                                 // calculate "leakage factor"
                                 const auto n = curLocalScvf.unitOuterNormal();
                                 const auto v = [&] ()
                                         {
                                             auto res = n;
-                                            res *= -0.5*facetVolVars.extrusionFactor();
+                                            res *= -0.5*completeFacetData.volVars().extrusionFactor();
                                             res -= curLocalScvf.ip();
                                             res += curLocalScvf.globalScvf().facetCorner();
                                             res /= res.two_norm2();
                                             return res;
                                         } ();
 
-                                // substract from matrix entry
-                                A[idxInFluxFaces][curIdxInFluxFaces] -= MpfaHelper::nT_M_v(n, posVolVars.permeability(), v);
+                                // substract n*T*v from diagonal matrix entry
+                                const auto facetTensor = getTensor(completeFacetData.problem(),
+                                                                   completeFacetData.element(),
+                                                                   completeFacetData.volVars(),
+                                                                   completeFacetData.fvGeometry(),
+                                                                   completeFacetData.scv());
+                                A[idxInFluxFaces][curIdxInFluxFaces] -= MpfaHelper::nT_M_v(n, facetTensor, v);
                             }
                         }
                         // this means we are on an interior face
@@ -593,7 +598,7 @@ private:
                     const auto& negGlobalScv = fvGeometry_().scv(negLocalScv.globalIndex());
                     const auto& negVolVars = elemVolVars_()[negGlobalScv];
                     const auto& negElement = localElement_(negLocalScvIdx);
-                    const auto negTensor = getTensor(negElement, negVolVars, negGlobalScv);
+                    const auto negTensor = getTensor(problem_(), negElement, negVolVars, fvGeometry_(), negGlobalScv);
 
                     // the omega factors of the "negative" sub volume
                     DimVector negWijk;
@@ -681,7 +686,7 @@ private:
             const auto& posGlobalScv = fvGeometry_().scv(posLocalScv.globalIndex());
             const auto& posVolVars = elemVolVars_()[posGlobalScv];
             const auto element = localElement_(posLocalScvIdx);
-            const auto tensor = getTensor(element, posVolVars, posGlobalScv);
+            const auto tensor = getTensor(problem_(), element, posVolVars, fvGeometry_(), posGlobalScv);
 
             // the omega factors of the "positive" sub volume
             auto posWijk = calculateOmegas_(posLocalScv, localScvf.unitOuterNormal(), localScvf.area(), tensor);

dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program 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 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief Helper class to be used to obtain lambda functions for the tensors
+ *        involved in the laws that describe the different kind of fluxes that
+ *        occur in DuMuX models (i.e. advective, diffusive and heat conduction fluxes).
+ *        The local systems appearing in Mpfa methods have to be solved subject to
+ *        the different tensors. This class returns lambda expressions to be used in the
+ *        local systems. The specialization for other discretization methods allows
+ *        compatibility with the TPFA scheme, which could be used for one or more of the tensors.
+ */
+#ifndef DUMUX_DISCRETIZATION_MPFA_TENSOR_LAMBDA_FACTORY_HH
+#define DUMUX_DISCRETIZATION_MPFA_TENSOR_LAMBDA_FACTORY_HH
+
+#include <dumux/implicit/properties.hh>
+#include <dumux/discretization/methods.hh>
+#include <dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh>
+
+namespace Dumux
+{
+
+//! forward declaration of properties
+namespace Properties
+{
+NEW_PROP_TAG(ThermalConductivityModel);
+};
+
+/*!
+ * \ingroup MpfaModel
+ * \brief Helper class to be used to obtain lambda functions for the tensors
+ *        involved in the laws that describe the different kind of fluxes that
+ *        occur in DuMuX models (i.e. advective, diffusive and heat conduction fluxes).
+ *        The local systems appearing in Mpfa methods have to be solved subject to
+ *        the different tensors. This class returns lambda expressions to be used in the
+ *        local systems. The specialization for discretization methods other than mpfa allows
+ *        compatibility with the TPFA scheme, which could be used for one or more of the tensors.
+ *        The interfaces of the lambdas are chosen such that all involved tensors can be extracted
+ *        with the given arguments.
+ */
+template<class TypeTag, DiscretizationMethods Method>
+class TensorLambdaFactory
+{
+public:
+
+    //! We return zero scalars here in the functions below.
+    //! We have to return something as the local systems expect a type
+    //! to perform actions on. We return 0.0 as this call should never happen
+    //! for a tensor which is not treated by an mpfa method anyway.
+
+    //! lambda for the law describing the advective term
+    static auto getAdvectionLambda()
+    {
+        return [] (const auto& problem,
+                   const auto& element,
+                   const auto& volVars,
+                   const auto& fvGeometry,
+                   const auto& scv)
+               { return 0.0; };
+    }
+
+    //! lambda for the diffusion coefficient/tensor
+    static auto getDiffusionLambda(unsigned int phaseIdx, unsigned int compIdx)
+    {
+        return [] (const auto& problem,
+                   const auto& element,
+                   const auto& volVars,
+                   const auto& fvGeometry,
+                   const auto& scv)
+               { return 0.0; };
+    }
+
+    //! lambda for the fourier coefficient
+    static auto getHeatConductionLambda()
+    {
+        return [] (const auto& problem,
+                   const auto& element,
+                   const auto& volVars,
+                   const auto& fvGeometry,
+                   const auto& scv)
+               { return 0.0; };
+    }
+};
+
+//! Specialization for mpfa schemes
+template<class TypeTag>
+class TensorLambdaFactory<TypeTag, DiscretizationMethods::CCMpfa>
+{
+public:
+
+    //! return void lambda here
+    static auto getAdvectionLambda()
+    {
+        return [] (const auto& problem,
+                   const auto& element,
+                   const auto& volVars,
+                   const auto& fvGeometry,
+                   const auto& scv)
+               { return volVars.permeability(); };
+    }
+
+    //! return void lambda here
+    static auto getDiffusionLambda(unsigned int phaseIdx, unsigned int compIdx)
+    {
+        return [phaseIdx, compIdx] (const auto& problem,
+                                    const auto& element,
+                                    const auto& volVars,
+                                    const auto& fvGeometry,
+                                    const auto& scv)
+               { return volVars.diffusionCoefficient(phaseIdx, compIdx); };
+    }
+
+    //! return void lambda here
+    static auto getHeatConductionLambda()
+    {
+        using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
+
+        return [] (const auto& problem,
+                   const auto& element,
+                   const auto& volVars,
+                   const auto& fvGeometry,
+                   const auto& scv)
+               { return ThermalConductivityModel::effectiveThermalConductivity(volVars,
+                                                                               problem.spatialParams(),
+                                                                               element,
+                                                                               fvGeometry,
+                                                                               scv); };
+    }
+};
+
+} // end namespace
+
+#endif