From 232c202e87d4ebba68ac9e10558707a6b4504560 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 26 Jul 2022 13:35:25 +0200
Subject: [PATCH] [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 | 78 +++++++++++++++++++++++-
1 file changed, 77 insertions(+), 1 deletion(-)
diff --git a/dumux/assembly/fcdiamondlocalresidual.hh b/dumux/assembly/fcdiamondlocalresidual.hh
index c122e09e63..3e6e1da4b6 100644
--- a/dumux/assembly/fcdiamondlocalresidual.hh
+++ b/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
--
GitLab