[fcdiamond][localresidual] Support time information in computeStorage

This allows for an ad-hoc fix of the current known bug in multidomain,
that does not correctly account for residual evaluations at different
time levels in a time integration scheme.

dumux/assembly/fcdiamondlocalresidual.hh

@@ -25,6 +25,7 @@
 #ifndef DUMUX_FC_DIAMOND_LOCAL_RESIDUAL_HH
 #define DUMUX_FC_DIAMOND_LOCAL_RESIDUAL_HH
 
+#include <dune/common/std/type_traits.hh>
 #include <dune/geometry/type.hh>
 #include <dune/istl/matrix.hh>
 
@@ -33,6 +34,22 @@
 #include <dumux/assembly/fvlocalresidual.hh>
 #include <dumux/discretization/extrusion.hh>
 
+namespace Dumux::Detail {
+
+template<class Imp, class P, class G, class S, class V>
+using TimeInfoInterfaceDiamondDetector = decltype(
+    std::declval<Imp>().computeStorage(
+        std::declval<P>(), std::declval<G>(), std::declval<S>(), std::declval<V>(), true
+    )
+);
+
+template<class Imp, class P, class G, class S, class V>
+constexpr inline bool hasTimeInfoInterfaceDiamond()
+{ return Dune::Std::is_detected<TimeInfoInterfaceDiamondDetector, Imp, P, G, S, V>::value; }
+
+} // end namespace Dumux::Detail
+
+
 namespace Dumux {
 
 /*!
@@ -45,6 +62,7 @@ template<class TypeTag>
 class FaceCenteredDiamondLocalResidual : public FVLocalResidual<TypeTag>
 {
     using ParentType = FVLocalResidual<TypeTag>;
+    using Implementation = GetPropType<TypeTag, Properties::LocalResidual>;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Problem = GetPropType<TypeTag, Properties::Problem>;
     using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
@@ -52,7 +70,9 @@ class FaceCenteredDiamondLocalResidual : public FVLocalResidual<TypeTag>
     using Element = typename GridView::template Codim<0>::Entity;
     using ElementBoundaryTypes = GetPropType<TypeTag, Properties::ElementBoundaryTypes>;
     using FVElementGeometry = typename GridGeometry::LocalView;
-    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
+    using GridVolumeVariables = GetPropType<TypeTag, Properties::GridVolumeVariables>;
+    using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
+    using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
     using SubControlVolume = typename FVElementGeometry::SubControlVolume;
     using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
     using ElementFluxVariablesCache = typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView;
@@ -129,6 +149,62 @@ public:
 
         return flux;
     }
+
+    using ParentType::evalStorage;
+    /*!
+     * \brief Compute the storage local residual, i.e. the deviation of the
+     *        storage term from zero for instationary problems.
+     *
+     * \param residual The residual vector to fill
+     * \param problem The problem to solve
+     * \param element The DUNE Codim<0> entity for which the residual
+     *                ought to be calculated
+     * \param fvGeometry The finite-volume geometry of the element
+     * \param prevElemVolVars The volume averaged variables for all
+     *                        sub-control volumes of the element at the previous time level
+     * \param curElemVolVars The volume averaged variables for all
+     *                       sub-control volumes of the element at the current  time level
+     * \param scv The sub control volume the storage term is integrated over
+     */
+    void evalStorage(ElementResidualVector& residual,
+                     const Problem& problem,
+                     const Element& element,
+                     const FVElementGeometry& fvGeometry,
+                     const ElementVolumeVariables& prevElemVolVars,
+                     const ElementVolumeVariables& curElemVolVars,
+                     const SubControlVolume& scv) const
+    {
+        const auto& curVolVars = curElemVolVars[scv];
+        const auto& prevVolVars = prevElemVolVars[scv];
+
+        // Compute storage with the model specific storage residual
+        // This addresses issues #792/#940 in ad-hoc way by additionally providing crude time level information (previous or current)
+        // to the low-level interfaces if this is supported by the LocalResidual implementation
+        NumEqVector prevStorage = computeStorageImpl_(problem, fvGeometry, scv, prevVolVars, /*previous time level?*/true);
+        NumEqVector storage = computeStorageImpl_(problem, fvGeometry, scv, curVolVars, /*previous time level?*/false);
+
+        prevStorage *= prevVolVars.extrusionFactor();
+        storage *= curVolVars.extrusionFactor();
+
+        storage -= prevStorage;
+        storage *= Extrusion::volume(scv);
+        storage /= this->timeLoop().timeStepSize();
+
+        residual[scv.localDofIndex()] += storage;
+    }
+
+private:
+    NumEqVector computeStorageImpl_(const Problem& problem,
+                                    const FVElementGeometry& fvGeometry,
+                                    const SubControlVolume& scv,
+                                    const VolumeVariables& volVars,
+                                    [[maybe_unused]] bool isPreviousTimeStep) const
+    {
+        if constexpr (Detail::hasTimeInfoInterfaceDiamond<Implementation, Problem, FVElementGeometry, SubControlVolume, VolumeVariables>())
+            return this->asImp().computeStorage(problem, fvGeometry, scv, volVars, isPreviousTimeStep);
+        else
+            return this->asImp().computeStorage(problem, scv, volVars);
+    }
 };
 
 } // end namespace Dumux