[mpfa] make fick's/fourier's law work & cleanup

dumux/discretization/cellcentered/mpfa/darcyslaw.hh

@@ -18,44 +18,35 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief This file contains the data which is required to calculate
+ * \brief This file contains the class which is required to calculate
  *        volume and mass fluxes of fluid phases over a face of a finite volume by means
- *        of the Darcy approximation. Specializations are provided for the different discretization methods.
+ *        of the Darcy approximation. This specializations is for cell-centered schemes
+ *        using multi-point flux approximation.
  */
 #ifndef DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
 #define DUMUX_DISCRETIZATION_CC_MPFA_DARCYS_LAW_HH
 
-namespace Dumux
-{
+#include <dumux/common/properties.hh>
 
-namespace Properties
+namespace Dumux
 {
-// forward declaration of properties
-NEW_PROP_TAG(MpfaHelper);
-}
 
 /*!
- * \ingroup DarcysLaw
+ * \ingroup Mpfa
  * \brief Specialization of Darcy's Law for the CCMpfa method.
  */
 template <class TypeTag>
 class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
 {
-    using Implementation = typename GET_PROP_TYPE(TypeTag, AdvectionType);
-
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using Element = typename GridView::template Codim<0>::Entity;
-    using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
-    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
-    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
     using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
-    using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
 
     // Always use the dynamic type for vectors (compatibility with the boundary)
     using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
@@ -125,23 +116,19 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
         }
 
         //! Returns the stencil for advective scvf flux computation
-        const Stencil& advectionVolVarsStencil() const
-        { return *advectionVolVarsStencil_; }
+        const Stencil& advectionVolVarsStencil() const { return *advectionVolVarsStencil_; }
 
         //! Returns the transmissibilities associated with the volume variables
         //! All phases flow through the same rock, thus, tij are equal for all phases
-        const CoefficientVector& advectionTij() const
-        { return *advectionTij_; }
+        const CoefficientVector& advectionTij() const { return *advectionTij_; }
 
         //! On faces that are "outside" w.r.t. a face in the interaction volume,
         //! we have to take the negative value of the fluxes, i.e. multiply by -1.0
-        bool advectionSwitchFluxSign() const
-        { return advectionSwitchFluxSign_; }
+        bool advectionSwitchFluxSign() const { return advectionSwitchFluxSign_; }
 
         //! Returns the data on dirichlet boundary conditions affecting
         //! the flux computation on this face
-        const DirichletDataContainer& advectionDirichletData() const
-        { return *advectionDirichletData_; }
+        const DirichletDataContainer& advectionDirichletData() const { return *advectionDirichletData_; }
 
     private:
         bool advectionSwitchFluxSign_;
@@ -167,17 +154,14 @@ public:
     {
         static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
 
-        const auto& fluxVarsCache = elemFluxVarsCache[scvf];
-        const auto& tij = fluxVarsCache.advectionTij();
-
         // Calculate the interface density for gravity evaluation
-        const auto rho = interpolateDensity(fvGeometry, elemVolVars, scvf, fluxVarsCache, phaseIdx);
-
-        // Variable for the flux to be computed
-        Scalar scvfFlux(0.0);
+        const auto rho = interpolateDensity(elemVolVars, scvf, phaseIdx);
 
-        // index counter to get corresponding transmissibilities
+        // prepare computations
         unsigned int i = 0;
+        Scalar scvfFlux(0.0);
+        const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+        const auto& tij = fluxVarsCache.advectionTij();
 
         // add contributions from cell-centered unknowns
         for (const auto volVarIdx : fluxVarsCache.advectionVolVarsStencil())
@@ -219,10 +203,8 @@ public:
     }
 
 private:
-    static Scalar interpolateDensity(const FVElementGeometry& fvGeometry,
-                                     const ElementVolumeVariables& elemVolVars,
+    static Scalar interpolateDensity(const ElementVolumeVariables& elemVolVars,
                                      const SubControlVolumeFace& scvf,
-                                     const FluxVariablesCache& fluxVarsCache,
                                      const unsigned int phaseIdx)
     {
         static const bool gravity = getParamFromGroup<bool>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Problem.EnableGravity");
@@ -235,7 +217,7 @@ private:
             if (!scvf.boundary())
             {
                 Scalar rho = elemVolVars[scvf.insideScvIdx()].density(phaseIdx);
-                for (auto outsideIdx : scvf.outsideScvIndices())
+                for (const auto outsideIdx : scvf.outsideScvIndices())
                     rho += elemVolVars[outsideIdx].density(phaseIdx);
                 return rho/(scvf.outsideScvIndices().size()+1);
             }

dumux/discretization/cellcentered/mpfa/dualgridindexset.hh

@@ -65,7 +65,7 @@ public:
                 const IndexType insideScvIdx,
                 const std::vector<IndexType>& outsideScvIndices)
     {
-        // check if it this really hasn't been called yet
+        // this should always be called only once
         assert( std::find(scvfIndices_.begin(), scvfIndices_.end(), scvfIdx ) == scvfIndices_.end() && "scvf has already been inserted!");
 
         // the local index of the scvf data about to be inserted
@@ -78,7 +78,8 @@ public:
         scvIndices.insert(scvIndices.end(), outsideScvIndices.begin(), outsideScvIndices.end());
 
         // if scvf is on boundary, increase counter
-        if (boundary) numBoundaryScvfs_++;
+        if (boundary)
+            numBoundaryScvfs_++;
 
         // insert data on the new scv
         scvfIndices_.push_back(scvfIdx);
@@ -88,10 +89,7 @@ public:
         // if entry for the inside scv exists append scvf local index, create otherwise
         auto it = std::find( scvIndices_.begin(), scvIndices_.end(), insideScvIdx );
         if (it != scvIndices_.end())
-        {
-            const auto localScvIdx = std::distance(scvIndices_.begin(), it);
-            localScvfIndicesInScv_[localScvIdx].push_back(curScvfLocalIdx);
-        }
+            localScvfIndicesInScv_[ std::distance(scvIndices_.begin(), it) ].push_back(curScvfLocalIdx);
         else
         {
             LocalIndexContainer localScvfIndices;
@@ -103,28 +101,25 @@ public:
     }
 
     //! returns the number of scvs
-    std::size_t numScvs() const
-    { return scvIndices_.size(); }
+    std::size_t numScvs() const { return scvIndices_.size(); }
 
     //! returns the number of scvfs
-    std::size_t numScvfs() const
-    { return scvfIndices_.size(); }
+    std::size_t numScvfs() const { return scvfIndices_.size(); }
 
     //! returns the number of boundary scvfs
-    std::size_t numBoundaryScvfs() const
-    { return numBoundaryScvfs_; }
+    std::size_t numBoundaryScvfs() const { return numBoundaryScvfs_; }
 
     //! returns the global scv indices connected to this dual grid node
-    const GlobalIndexContainer& globalScvIndices() const
-    { return scvIndices_; }
+    const GlobalIndexContainer& globalScvIndices() const { return scvIndices_; }
 
     //! returns the global scvf indices connected to this dual grid node
-    const GlobalIndexContainer& globalScvfIndices() const
-    { return scvfIndices_; }
+    const GlobalIndexContainer& globalScvfIndices() const { return scvfIndices_; }
 
     //! returns the global scv idx of the i-th scv
-    IndexType scvIdxGlobal(unsigned int i) const
-    { return scvIndices_[i]; }
+    IndexType scvIdxGlobal(unsigned int i) const { return scvIndices_[i]; }
+
+    //! returns the index of the i-th scvf
+    IndexType scvfIdxGlobal(unsigned int i) const { return scvfIndices_[i]; }
 
     //! returns the global index of the j-th scvf embedded in the i-th scv
     IndexType scvfIdxGlobal(unsigned int i, unsigned int j) const
@@ -140,13 +135,8 @@ public:
         return localScvfIndicesInScv_[i][j];
     }
 
-    //! returns the index of the i-th scvf
-    IndexType scvfIdxGlobal(unsigned int i) const
-    { return scvfIndices_[i]; }
-
     //! returns whether or not the i-th scvf touches the boundary
-    bool scvfIsOnBoundary(unsigned int i) const
-    { return scvfIsOnBoundary_[i]; }
+    bool scvfIsOnBoundary(unsigned int i) const { return scvfIsOnBoundary_[i]; }
 
     //! returns the indices of the neighboring scvs of the i-th scvf
     const GlobalIndexContainer& neighboringScvIndices(unsigned int i) const
@@ -174,7 +164,7 @@ private:
  * \brief Class for the index sets of the dual grid in mpfa schemes.
  */
 template<class TypeTag>
-class CCMpfaDualGridIndexSet : public std::vector<DualGridNodalIndexSet<TypeTag>>
+class CCMpfaDualGridIndexSet
 {
     using ParentType = std::vector<DualGridNodalIndexSet<TypeTag>>;
 
@@ -187,19 +177,26 @@ class CCMpfaDualGridIndexSet : public std::vector<DualGridNodalIndexSet<TypeTag>
 public:
     using NodalIndexSet = DualGridNodalIndexSet<TypeTag>;
 
-    CCMpfaDualGridIndexSet(const GridView& gridView)
-    {
-        (*this).resize(gridView.size(dim));
-    }
+    //! default constructor
+    CCMpfaDualGridIndexSet() = delete;
+
+    //! constructor
+    explicit CCMpfaDualGridIndexSet(const GridView& gridView) : nodalIndexSets_(gridView.size(dim)) {}
 
     //! access with an scvf
     const NodalIndexSet& operator[] (const SubControlVolumeFace& scvf) const
-    { return (*this)[scvf.vertexIndex()]; }
+    { return nodalIndexSets_[scvf.vertexIndex()]; }
     NodalIndexSet& operator[] (const SubControlVolumeFace& scvf)
-    { return (*this)[scvf.vertexIndex()]; }
+    { return nodalIndexSets_[scvf.vertexIndex()]; }
 
     //! access with an index
-    using ParentType::operator[];
+    const NodalIndexSet& operator[] (IndexType i) const
+    { return nodalIndexSets_[i]; }
+    NodalIndexSet& operator[] (IndexType i)
+    { return nodalIndexSets_[i]; }
+
+private:
+    std::vector<NodalIndexSet> nodalIndexSets_;
 };
 } // end namespace
 

dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh

@@ -29,13 +29,6 @@
 namespace Dumux
 {
 
-//! forward declaration of properties
-namespace Properties
-{
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(NumComponents);
-};
-
 /*!
  * \ingroup ImplicitModel
  * \brief Helper class to fill the flux var caches

dumux/discretization/cellcentered/mpfa/fourierslaw.hh

@@ -18,34 +18,28 @@
  *****************************************************************************/
 /*!
  * \file
- * \brief This file contains the data which is required to calculate
+ * \brief This file contains the class which is required to calculate
  *        heat conduction fluxes with Fourier's law for cell-centered MPFA models.
  */
 #ifndef DUMUX_DISCRETIZATION_CC_MPFA_FOURIERS_LAW_HH
 #define DUMUX_DISCRETIZATION_CC_MPFA_FOURIERS_LAW_HH
 
-#include <dune/common/float_cmp.hh>
-
-#include <dumux/common/math.hh>
-#include <dumux/common/parameters.hh>
-#include <dumux/implicit/properties.hh>
-#include <dumux/discretization/methods.hh>
+#include <dumux/common/properties.hh>
 
 namespace Dumux
 {
 
 /*!
- * \ingroup FouriersLaw
+ * \ingroup Mpfa
  * \brief Specialization of Fourier's Law for the CCMpfa method.
  */
 template <class TypeTag>
 class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
 {
-    using Implementation = typename GET_PROP_TYPE(TypeTag, HeatConductionType);
-
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
@@ -54,16 +48,12 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
     using ThermalConductivityModel = typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel);
 
     // Always use the dynamic type for vectors (compatibility with the boundary)
-    using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
-    using CoefficientVector = typename BoundaryInteractionVolume::Vector;
-
-    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-    using Element = typename GridView::template Codim<0>::Entity;
-    using IndexType = typename GridView::IndexSet::IndexType;
-
-    static constexpr bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
-    static constexpr bool enableInteriorBoundaries = GET_PROP_VALUE(TypeTag, EnableInteriorBoundaries);
+    using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
+    using CoefficientVector = typename PrimaryInteractionVolume::Traits::DynamicVector;
+    using DataHandle = typename PrimaryInteractionVolume::Traits::DataHandle;
 
+    static constexpr int dim = GridView::dimension;
+    static constexpr int dimWorld = GridView::dimensionworld;
     static constexpr int energyEqIdx = GET_PROP_TYPE(TypeTag, Indices)::energyEqIdx;
 
     //! Class that fills the cache corresponding to mpfa Darcy's Law
@@ -81,53 +71,67 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
                          const SubControlVolumeFace& scvf,
                          const FluxVariablesCacheFiller& fluxVarsCacheFiller)
         {
-            // get interaction volume from the flux vars cache filler & upate the cache
-            if (problem.model().fvGridGeometry().isInBoundaryInteractionVolume(scvf))
-                scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.boundaryInteractionVolume(), scvf);
-            else
-                scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.interactionVolume(), scvf);
+          // get interaction volume from the flux vars cache filler & upate the cache
+          if (fvGeometry.fvGridGeometry().vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
+              scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.secondaryInteractionVolume(),
+                                                     fluxVarsCacheFiller.dataHandle(),
+                                                     scvf);
+          else
+              scvfFluxVarsCache.updateHeatConduction(fluxVarsCacheFiller.primaryInteractionVolume(),
+                                                     fluxVarsCacheFiller.dataHandle(),
+                                                     scvf);
         }
     };
 
     //! The cache used in conjunction with the mpfa Fourier's Law
     class MpfaFouriersLawCache
     {
-        using Stencil = typename BoundaryInteractionVolume::GlobalIndexSet;
+        // We always use the dynamic types here to be compatible on the boundary
+        using Stencil = typename PrimaryInteractionVolume::Traits::DynamicGlobalIndexContainer;
+        using DirichletDataContainer = typename PrimaryInteractionVolume::DirichletDataContainer;
     public:
         // export filler type
         using Filler = MpfaFouriersLawCacheFiller;
 
         // update cached objects for heat conduction
-        template<typename InteractionVolume>
-        void updateHeatConduction(const InteractionVolume& iv, const SubControlVolumeFace &scvf)
+        template<class InteractionVolume>
+        void updateHeatConduction(const InteractionVolume& iv, const DataHandle& dataHandle, const SubControlVolumeFace &scvf)
         {
             const auto& localFaceData = iv.getLocalFaceData(scvf);
-            heatConductionVolVarsStencil_ = iv.volVarsStencil();
-            heatConductionTij_ = iv.getTransmissibilities(localFaceData);
-            heatNeumannFlux_ = iv.getNeumannFlux(localFaceData, energyEqIdx);
+
+            // update the quantities that are equal for all phases
+            heatConductionSwitchFluxSign_ = localFaceData.isOutside();
+            heatConductionVolVarsStencil_ = &dataHandle.volVarsStencil();
+            heatConductionDirichletData_ = &dataHandle.dirichletData();
+
+            // the transmissibilities on surface grids have to be obtained from the outside
+            if (dim == dimWorld)
+                heatConductionTij_ = &dataHandle.T()[localFaceData.ivLocalScvfIndex()];
+            else
+                heatConductionTij_ = localFaceData.isOutside() ?
+                                     &dataHandle.outsideTij()[localFaceData.ivLocalOutsideScvfIndex()] :
+                                     &dataHandle.T()[localFaceData.ivLocalScvfIndex()];
         }
 
-        //! Returns the volume variables indices necessary for heat conduction flux
-        //! computation. This includes all participating boundary volume variables
-        //! and it can be different for the phases & components.
-        const Stencil& heatConductionVolVarsStencil() const
-        { return heatConductionVolVarsStencil_; }
+        //! Returns the stencil for heat conduction flux computation on an scvf
+        const Stencil& heatConductionVolVarsStencil() const { return *heatConductionVolVarsStencil_; }
 
         //! Returns the transmissibilities associated with the volume variables
-        //! This can be different for the phases & components.
-        const CoefficientVector& heatConductionTij() const
-        { return heatConductionTij_; }
+        const CoefficientVector& heatConductionTij() const { return *heatConductionTij_; }
 
-        //! If the useTpfaBoundary property is set to false, the boundary conditions
-        //! are put into the local systems leading to possible contributions on all faces
-        Scalar heatNeumannFlux() const
-        { return heatNeumannFlux_; }
+        //! On faces that are "outside" w.r.t. a face in the interaction volume,
+        //! we have to take the negative value of the fluxes, i.e. multiply by -1.0
+        bool heatConductionSwitchFluxSign() const { return heatConductionSwitchFluxSign_; }
+
+        //! Returns the data on dirichlet boundary conditions affecting
+        //! the flux computation on this face
+        const DirichletDataContainer& heatConductionDirichletData() const { return *heatConductionDirichletData_; }
 
     private:
-        // Quantities associated with heat conduction
-        Stencil heatConductionVolVarsStencil_;
-        CoefficientVector heatConductionTij_;
-        Scalar heatNeumannFlux_;
+        bool heatConductionSwitchFluxSign_;
+        const Stencil* heatConductionVolVarsStencil_;
+        const CoefficientVector* heatConductionTij_;
+        const DirichletDataContainer* heatConductionDirichletData_;
     };
 
 public:
@@ -144,54 +148,22 @@ public:
                        const SubControlVolumeFace& scvf,
                        const ElementFluxVarsCache& elemFluxVarsCache)
     {
+        // prepare computations
+        Scalar flux(0.0);
+        unsigned int i = 0;
         const auto& fluxVarsCache = elemFluxVarsCache[scvf];
-        const auto& volVarsStencil = fluxVarsCache.heatConductionVolVarsStencil();
         const auto& tij = fluxVarsCache.heatConductionTij();
 
-        const bool isInteriorBoundary = enableInteriorBoundaries && fluxVarsCache.isInteriorBoundary();
-        // For interior Neumann boundaries when using Tpfa on boundaries, return the user-specified flux
-        if (isInteriorBoundary
-            && useTpfaBoundary
-            && fluxVarsCache.interiorBoundaryDataSelf().faceType() == MpfaFaceTypes::interiorNeumann)
-            return scvf.area()*
-                   elemVolVars[scvf.insideScvIdx()].extrusionFactor()*
-                   problem.neumann(element, fvGeometry, elemVolVars, scvf)[energyEqIdx];
-
         // calculate Tij*tj