diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
index d15a930b733733624d9406bb8c876ff6a6c0796c..dfb3ad9b6f8bb79fa3ce5b220c8af48df4d4716e 100644
--- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
+++ b/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);
diff --git a/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh b/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh
index c1cc25b6b088e56bd0b6191f481f54537418ace9..3244b8c74e3b2bf81d07001e2003fdaf3fe8b427 100644
--- a/dumux/discretization/cellcentered/mpfa/interiorboundarydata.hh
+++ b/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");
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index 7d3bb189bc6d6e1536f5dd0bced2d932b84c8f43..444fc6a823bf5866652f14127c7957da6f390c29 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/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);
diff --git a/dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh b/dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh
new file mode 100644
index 0000000000000000000000000000000000000000..548cd6598e072cd06ba515c80b4d60f6bff6d95a
--- /dev/null
+++ b/dumux/discretization/cellcentered/mpfa/tensorlambdafactory.hh
@@ -0,0 +1,149 @@
+// -*- 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