[staggered] Fix caching of fluxVariablesCache

* empty cache filler used
* contains several todos
* fluxVarsCache not used anywehere, but might be required in future

dumux/discretization/staggered/elementfluxvariablescache.hh

@@ -53,6 +53,9 @@ public:
     //! export the type of the flux variables cache
     using FluxVariablesCache = typename GFVC::FluxVariablesCache;
 
+    //! export the type of the flux variables cache filler
+    using FluxVariablesCacheFiller = typename GFVC::FluxVariablesCacheFiller;
+
     StaggeredElementFluxVariablesCache(const GridFluxVariablesCache& global)
     : gridFluxVarsCachePtr_(&global) {}
 
@@ -75,6 +78,15 @@ public:
                   const ElementVolumeVariables& elemVolVars,
                   const typename FVElementGeometry::SubControlVolumeFace& scvf) {}
 
+    //! Specialization for the global caching being enabled - do nothing here
+    template<class FVElementGeometry, class ElementVolumeVariables>
+    void update(const typename FVElementGeometry::FVGridGeometry::GridView::template Codim<0>::Entity& element,
+                const FVElementGeometry& fvGeometry,
+                const ElementVolumeVariables& elemVolVars)
+    {
+        DUNE_THROW(Dune::InvalidStateException, "In case of enabled caching, the grid flux variables cache has to be updated");
+    }
+
     //! operators in the case of caching
     template<class SubControlVolumeFace>
     const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const
@@ -103,11 +115,17 @@ public:
     //! export the type of the flux variables cache
     using FluxVariablesCache = typename GFVC::FluxVariablesCache;
 
+    //! export the type of the flux variables cache filler
+    using FluxVariablesCacheFiller = typename GFVC::FluxVariablesCacheFiller;
+
     StaggeredElementFluxVariablesCache(const GridFluxVariablesCache& global)
     : gridFluxVarsCachePtr_(&global) {}
 
-    //! This function has to be called prior to flux calculations on the element.
-    //! Prepares the transmissibilities of the scv faces in an element. The FvGeometry is assumed to be bound.
+    /*!
+     * \brief Prepares the transmissibilities of the scv faces in an element
+     * \note the fvGeometry is assumed to be bound to the same element
+     * \note this function has to be called prior to flux calculations on the element.
+     */
     template<class FVElementGeometry, class ElementVolumeVariables>
     void bindElement(const typename FVElementGeometry::FVGridGeometry::GridView::template Codim<0>::Entity& element,
                      const FVElementGeometry& fvGeometry,
@@ -117,63 +135,55 @@ public:
         const auto numScvf = fvGeometry.numScvf();
         fluxVarsCache_.resize(numScvf);
         globalScvfIndices_.resize(numScvf);
-        std::size_t localScvfIdx = 0;
 
+        // instantiate helper class to fill the caches
+        // FluxVariablesCacheFiller filler(gridFluxVarsCache().problem()); TODO: use proper ctor
+        FluxVariablesCacheFiller filler;
+
+        std::size_t localScvfIdx = 0;
         // fill the containers
         for (auto&& scvf : scvfs(fvGeometry))
         {
-            fluxVarsCache_[localScvfIdx].update(gridFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf);
+            filler.fill(*this, fluxVarsCache_[localScvfIdx], element, fvGeometry, elemVolVars, scvf, true);
             globalScvfIndices_[localScvfIdx] = scvf.index();
             localScvfIdx++;
         }
     }
 
-    //! This function is called by the StaggeredLocalResidual before flux calculations during assembly.
-    //! Prepares the transmissibilities of the scv faces in the stencil. The FvGeometries are assumed to be bound.
+    /*!
+     * \brief Prepares the transmissibilities of the scv faces in the stencil of an element
+     * \note the fvGeometry is assumed to be bound to the same element
+     * \note this function has to be called prior to flux calculations on the element.
+     */
     template<class FVElementGeometry, class ElementVolumeVariables>
     void bind(const typename FVElementGeometry::FVGridGeometry::GridView::template Codim<0>::Entity& element,
               const FVElementGeometry& fvGeometry,
               const ElementVolumeVariables& elemVolVars)
     {
-        const auto globalI = gridFluxVarsCache().problem_().elementMapper().index(element);
-        const auto& neighborStencil = gridFluxVarsCache().problem_().model().stencils(element).neighborStencil();
-        const auto numNeighbors = neighborStencil.size();
+        // instantiate helper class to fill the caches
+        // FluxVariablesCacheFiller filler(problem); TODO: use proper ctor
+        FluxVariablesCacheFiller filler;
 
         // find the number of scv faces that need to be prepared
-        auto numScvf = fvGeometry.numScvf();
-        for (std::size_t localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ)
-        {
-            const auto& fluxVarIndicesJ = gridFluxVarsCache().problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ];
-            numScvf += fluxVarIndicesJ.size();
-        }
+        const auto numScvf = fvGeometry.numScvf();
 
         // fill the containers with the data on the scv faces inside the actual element
         fluxVarsCache_.resize(numScvf);
         globalScvfIndices_.resize(numScvf);
-        std::size_t localScvfIdx = 0;
+        unsigned int localScvfIdx = 0;
         for (auto&& scvf : scvfs(fvGeometry))
         {
-            fluxVarsCache_[localScvfIdx].update(gridFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf);
+            filler.fill(*this, fluxVarsCache_[localScvfIdx], element, fvGeometry, elemVolVars, scvf, true);
             globalScvfIndices_[localScvfIdx] = scvf.index();
             localScvfIdx++;
         }
-
-        // add required data on the scv faces in the neighboring elements
-        for (std::size_t localIdxJ = 0; localIdxJ < numNeighbors; ++localIdxJ)
-        {
-            const auto& fluxVarIndicesJ = gridFluxVarsCache().problem_().model().localJacobian().assemblyMap()[globalI][localIdxJ];
-            const auto elementJ = fvGeometry.fvGridGeometry().element(neighborStencil[localIdxJ]);
-
-            for (auto fluxVarIdx : fluxVarIndicesJ)
-            {
-                auto&& scvfJ = fvGeometry.scvf(fluxVarIdx);
-                fluxVarsCache_[localScvfIdx].update(gridFluxVarsCache().problem_(), elementJ, fvGeometry, elemVolVars, scvfJ);
-                globalScvfIndices_[localScvfIdx] = scvfJ.index();
-                localScvfIdx++;
-            }
-        }
     }
 
+    /*!
+     * \brief Prepares the transmissibilities of a single scv face
+     * \note the fvGeometry is assumed to be bound to the same element
+     * \note this function has to be called prior to flux calculations on the element.
+     */
     template<class FVElementGeometry, class ElementVolumeVariables>
     void bindScvf(const typename FVElementGeometry::FVGridGeometry::GridView::template Codim<0>::Entity& element,
                   const FVElementGeometry& fvGeometry,
@@ -183,10 +193,46 @@ public:
         fluxVarsCache_.resize(1);
         globalScvfIndices_.resize(1);
 
-        fluxVarsCache_[0].update(gridFluxVarsCache().problem_(), element, fvGeometry, elemVolVars, scvf);
+        // instantiate helper class to fill the caches
+        // FluxVariablesCacheFiller filler(gridFluxVarsCache().problem());
+        FluxVariablesCacheFiller filler; // TODO: use proper ctor
+
+        filler.fill(*this, fluxVarsCache_[0], element, fvGeometry, elemVolVars, scvf, true);
         globalScvfIndices_[0] = scvf.index();
     }
 
+    /*!
+     * \brief Update the transmissibilities if the volume variables have changed
+     * \note Results in undefined behaviour if called before bind() or with a different element
+     */
+    template<class FVElementGeometry, class ElementVolumeVariables>
+    void update(const typename FVElementGeometry::FVGridGeometry::GridView::template Codim<0>::Entity& element,
+                const FVElementGeometry& fvGeometry,
+                const ElementVolumeVariables& elemVolVars)
+    {
+        // if (FluxVariablesCacheFiller::isSolDependent) TODO
+        // {
+        //     const auto& problem = gridFluxVarsCache().problem();
+        //     const auto globalI = fvGeometry.fvGridGeometry().elementMapper().index(element);
+        //
+        //     // instantiate filler class
+        //     FluxVariablesCacheFiller filler(problem);
+        //
+        //     // let the filler class update the caches
+        //     for (unsigned int localScvfIdx = 0; localScvfIdx < fluxVarsCache_.size(); ++localScvfIdx)
+        //     {
+        //         const auto& scvf = fvGeometry.scvf(globalScvfIndices_[localScvfIdx]);
+        //
+        //         const auto scvfInsideScvIdx = scvf.insideScvIdx();
+        //         const auto& insideElement = scvfInsideScvIdx == globalI ?
+        //                                     element :
+        //                                     fvGeometry.fvGridGeometry().element(scvfInsideScvIdx);
+        //
+        //         filler.fill(*this, fluxVarsCache_[localScvfIdx], insideElement, fvGeometry, elemVolVars, scvf);
+        //     }
+        // }
+    }
+
     //! access operators in the case of no caching
     template<class SubControlVolumeFace>
     const FluxVariablesCache& operator [](const SubControlVolumeFace& scvf) const

dumux/discretization/staggered/gridfluxvariablescache.hh

@@ -36,11 +36,12 @@ namespace Dumux {
  * \tparam P The problem type
  * \tparam FVC The flux variables cache type
  */
-template<class P, class FVC>
+template<class P, class FVC, class FVCF>
 struct StaggeredDefaultGridFluxVariablesCacheTraits
 {
     using Problem = P;
     using FluxVariablesCache = FVC;
+    using FluxVariablesCacheFiller = FVCF;
 
     template<class GridFluxVariablesCache, bool cachingEnabled>
     using LocalView = StaggeredElementFluxVariablesCache<GridFluxVariablesCache, cachingEnabled>;
@@ -52,8 +53,9 @@ struct StaggeredDefaultGridFluxVariablesCacheTraits
  */
 template<class Problem,
          class FluxVariablesCache,
-         bool cachingEnabled = false,
-         class Traits = StaggeredDefaultGridFluxVariablesCacheTraits<Problem, FluxVariablesCache> >
+         class FluxVariablesCacheFiller,
+         bool EnableGridFluxVariablesCache = false,
+         class Traits = StaggeredDefaultGridFluxVariablesCacheTraits<Problem, FluxVariablesCache, FluxVariablesCacheFiller> >
 class StaggeredGridFluxVariablesCache;
 
 /*!
@@ -61,16 +63,19 @@ class StaggeredGridFluxVariablesCache;
  * \brief Flux variables cache class for staggered models.
           Specialization in case of storing the flux cache.
  */
-template<class P, class FVC, class Traits>
-class StaggeredGridFluxVariablesCache<P, FVC, true, Traits>
+template<class P, class FVC, class FVCF, class Traits>
+class StaggeredGridFluxVariablesCache<P, FVC, FVCF, true, Traits>
 {
-    using ThisType = StaggeredGridFluxVariablesCache<P, FVC, true, Traits>;
     using Problem = typename Traits::Problem;
+    using ThisType = StaggeredGridFluxVariablesCache<P, FVC, FVCF, true, Traits>;
 
 public:
     //! export the flux variable cache type
     using FluxVariablesCache = typename Traits::FluxVariablesCache;
 
+    //! export the flux variable cache filler type
+    using FluxVariablesCacheFiller = typename Traits::FluxVariablesCacheFiller;
+
     //! make it possible to query if caching is enabled
     static constexpr bool cachingEnabled = true;
 
@@ -86,21 +91,31 @@ public:
                 const SolutionVector& sol,
                 bool forceUpdate = false)
     {
-        // fluxVarsCache_.resize(fvGridGeometry.numScvf());
-        // for (const auto& element : elements(fvGridGeometry.gridView()))
-        // {
-        //     // Prepare the geometries within the elements of the stencil
-        //     auto fvGeometry = localView(fvGridGeometry);
-        //     fvGeometry.bind(element);
-        //
-        //     auto elemVolVars = localView(gridVolVars);
-        //     elemVolVars.bind(element, fvGeometry, sol);
-        //
-        //     for (auto&& scvf : scvfs(fvGeometry))
-        //     {
-        //         fluxVarsCache_[scvf.index()].update(problem, element, fvGeometry, elemVolVars, scvf);
-        //     }
-        // }
+        // only do the update if fluxes are solution dependent or if update is forced
+        // TODO: so far, the staggered models do not use any fluxVar caches, therefore an empty cache filler
+        // is used which does not implement isSolDependent
+        if (/*FluxVariablesCacheFiller::isSolDependent ||*/ forceUpdate)
+        {
+            // instantiate helper class to fill the caches
+            // FluxVariablesCacheFiller filler(problem()); TODO: use proper ctor
+            FluxVariablesCacheFiller filler;
+
+            fluxVarsCache_.resize(fvGridGeometry.numScvf());
+            for (const auto& element : elements(fvGridGeometry.gridView()))
+            {
+                // Prepare the geometries within the elements of the stencil
+                auto fvGeometry = localView(fvGridGeometry);
+                fvGeometry.bind(element);
+
+                auto elemVolVars = localView(gridVolVars);
+                elemVolVars.bind(element, fvGeometry, sol);
+
+                for (auto&& scvf : scvfs(fvGeometry))
+                {
+                    filler.fill(*this, fluxVarsCache_[scvf.index()], element, fvGeometry, elemVolVars, scvf, forceUpdate);
+                }
+            }
+        }
     }
 
     const Problem& problem() const
@@ -125,31 +140,40 @@ private:
  * \brief Flux variables cache class for staggered models.
           Specialization in case of not storing the flux cache.
  */
-template<class P, class FVC, class Traits>
-class StaggeredGridFluxVariablesCache<P, FVC, false, Traits>
+template<class P, class FVC, class FVCF, class Traits>
+class StaggeredGridFluxVariablesCache<P, FVC, FVCF, false, Traits>
 {
-    using ThisType = StaggeredGridFluxVariablesCache<P, FVC, false, Traits>;
     using Problem = typename Traits::Problem;
+    using ThisType = StaggeredGridFluxVariablesCache<P, FVC, FVCF, false, Traits>;
 
 public:
     //! export the flux variable cache type
     using FluxVariablesCache = typename Traits::FluxVariablesCache;
 
+    //! export the flux variable cache filler type
+    using FluxVariablesCacheFiller = typename Traits::FluxVariablesCacheFiller;
+
     //! make it possible to query if caching is enabled
-    static constexpr bool cachingEnabled = true;
+    static constexpr bool cachingEnabled = false;
 
     //! export the type of the local view
     using LocalView = typename Traits::template LocalView<ThisType, cachingEnabled>;
 
-    // When global flux variables caching is disabled, we don't need to update the cache
-    void update(Problem& problem)
-    { problemPtr_ = &problem; }
+    StaggeredGridFluxVariablesCache(const Problem& problem) : problemPtr_(&problem) {}
 
-private:
+    // When global caching is enabled, precompute transmissibilities and stencils for all the scv faces
+    template<class FVGridGeometry, class GridVolumeVariables, class SolutionVector>
+    void update(const FVGridGeometry& fvGridGeometry,
+                const GridVolumeVariables& gridVolVars,
+                const SolutionVector& sol,
+                bool forceUpdate = false) {}
 
-    const Problem& problem_() const
+    const Problem& problem() const
     { return *problemPtr_; }
 
+private:
+
+
     const Problem* problemPtr_;
 };
 

dumux/discretization/staggered/properties.hh

@@ -32,6 +32,7 @@
 
 #include <dumux/discretization/methods.hh>
 #include <dumux/discretization/fvproperties.hh>
+#include <dumux/discretization/fluxvariablescaching.hh>
 
 #include <dumux/discretization/cellcentered/elementboundarytypes.hh>
 #include <dumux/assembly/staggeredlocalresidual.hh>
@@ -78,8 +79,9 @@ private:
     static constexpr auto enableCache = GET_PROP_VALUE(TypeTag, EnableGridFluxVariablesCache);
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
+    using FluxVariablesCacheFiller = FluxVariablesCaching::EmptyCacheFiller;
 public:
-    using type = StaggeredGridFluxVariablesCache<Problem, FluxVariablesCache, enableCache>;
+    using type = StaggeredGridFluxVariablesCache<Problem, FluxVariablesCache, FluxVariablesCacheFiller, enableCache>;
 };
 
 //! Set the face solution type