().computeStorage(
+ std::declval(), std::declval(), std::declval(), std::declval(), true
+ )
+);
+
+template
+constexpr inline bool hasTimeInfoInterfaceCVFE()
+{ return Dune::Std::is_detected::value; }
+
+template
+using SCVFIsOverlappingDetector = decltype(
+ std::declval().isOverlapping()
+);
+
+template
+constexpr inline bool hasScvfIsOverlapping()
+{ return Dune::Std::is_detected::value; }
+
+} // end namespace Dumux::Detail
+
+
+namespace Dumux {
+
+/*!
+ * \ingroup Assembly
+ * \ingroup CVFEDiscretization
+ * \brief The element-wise residual for control-volume finite element schemes
+ * \tparam TypeTag the TypeTag
+ */
+template
+class CVFELocalResidual : public FVLocalResidual
+{
+ using ParentType = FVLocalResidual;
+ using Implementation = GetPropType;
+ using Scalar = GetPropType;
+ using Problem = GetPropType;
+ using GridGeometry = GetPropType;
+ using GridView = typename GridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementBoundaryTypes = GetPropType;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using GridVolumeVariables = GetPropType;
+ using ElementVolumeVariables = typename GridVolumeVariables::LocalView;
+ using VolumeVariables = typename GridVolumeVariables::VolumeVariables;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using ElementFluxVariablesCache = typename GetPropType::LocalView;
+ using PrimaryVariables = GetPropType;
+ using NumEqVector = Dumux::NumEqVector;
+ using Extrusion = Extrusion_t;
+
+public:
+ using ElementResidualVector = typename ParentType::ElementResidualVector;
+ using ParentType::ParentType;
+
+ //! evaluate flux residuals for one sub control volume face and add to residual
+ void evalFlux(ElementResidualVector& residual,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementBoundaryTypes& elemBcTypes,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ const auto flux = evalFlux(problem, element, fvGeometry, elemVolVars, elemBcTypes, elemFluxVarsCache, scvf);
+ if (!scvf.boundary())
+ {
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+ residual[insideScv.localDofIndex()] += flux;
+
+ // for control-volume finite element schemes with overlapping control volumes
+ if constexpr (Detail::hasScvfIsOverlapping())
+ {
+ if (!scvf.isOverlapping())
+ residual[outsideScv.localDofIndex()] -= flux;
+ }
+ else
+ residual[outsideScv.localDofIndex()] -= flux;
+ }
+ else
+ {
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ residual[insideScv.localDofIndex()] += flux;
+ }
+ }
+
+ //! evaluate flux residuals for one sub control volume face
+ NumEqVector evalFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementBoundaryTypes& elemBcTypes,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const SubControlVolumeFace& scvf) const
+ {
+ NumEqVector flux(0.0);
+
+ // inner faces
+ if (!scvf.boundary())
+ flux += this->asImp().computeFlux(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+
+ // boundary faces
+ else
+ {
+ const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& bcTypes = elemBcTypes.get(fvGeometry, scv);
+
+ // Treat Neumann and Robin ("solution dependent Neumann") boundary conditions.
+ // For Dirichlet there is no addition to the residual here but they
+ // are enforced strongly by replacing the residual entry afterwards.
+ if (bcTypes.hasNeumann())
+ {
+ auto neumannFluxes = problem.neumann(element, fvGeometry, elemVolVars, elemFluxVarsCache, scvf);
+
+ // multiply neumann fluxes with the area and the extrusion factor
+ neumannFluxes *= Extrusion::area(fvGeometry, scvf)*elemVolVars[scv].extrusionFactor();
+
+ // only add fluxes to equations for which Neumann is set
+ for (int eqIdx = 0; eqIdx < NumEqVector::dimension; ++eqIdx)
+ if (bcTypes.isNeumann(eqIdx))
+ flux[eqIdx] += neumannFluxes[eqIdx];
+ }
+ }
+
+ 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(fvGeometry, 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::hasTimeInfoInterfaceCVFE())
+ return this->asImp().computeStorage(problem, fvGeometry, scv, volVars, isPreviousTimeStep);
+ else
+ return this->asImp().computeStorage(problem, scv, volVars);
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/diffmethod.hh b/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/diffmethod.hh
new file mode 100644
index 0000000..86295c6
--- /dev/null
+++ b/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/diffmethod.hh
@@ -0,0 +1,31 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief An enum class to define various differentiation methods available in order to compute
+ the derivatives of the residual i.e. the entries in the jacobian matrix.
+ */
+#ifndef DUMUX_JACOBIAN_DIFFERENTIATION_METHODS_HH
+#define DUMUX_JACOBIAN_DIFFERENTIATION_METHODS_HH
+
+namespace Dumux {
+
+/*!
+ * \ingroup Assembly
+ * \brief Differentiation methods in order to compute the derivatives
+ * of the residual i.e. the entries in the jacobian matrix.
+ * \todo automatic differentation is not yet implemented
+ */
+enum class DiffMethod
+{
+ numeric, analytic, automatic
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/entitycolor.hh b/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/entitycolor.hh
new file mode 100644
index 0000000..33964c7
--- /dev/null
+++ b/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/entitycolor.hh
@@ -0,0 +1,45 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief An enum class to define the colors of elements and vertices required
+ * for partial Jacobian reassembly.
+ */
+#ifndef DUMUX_ENTITY_COLOR_HH
+#define DUMUX_ENTITY_COLOR_HH
+
+namespace Dumux {
+
+/*!
+ * \ingroup Assembly
+ * \brief The colors of elements and vertices required for partial
+ * Jacobian reassembly.
+ */
+enum class EntityColor {
+ //! distance from last linearization is above the tolerance
+ red,
+
+ //! neighboring entity is red
+ yellow,
+
+ /*!
+ * A yellow entity that has only non-green neighbor elements.
+ *
+ * This means that its relative error is below the tolerance,
+ * but its defect can be linearized without any additional
+ * cost.
+ */
+ orange,
+
+ //! does not need to be reassembled
+ green
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/fclocalassembler.hh b/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/fclocalassembler.hh
new file mode 100644
index 0000000..148dd47
--- /dev/null
+++ b/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/fclocalassembler.hh
@@ -0,0 +1,781 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup Assembly
+ * \ingroup FaceCenteredStaggeredDiscretization
+ * \brief An assembler for Jacobian and residual contribution per element (face-centered staggered methods)
+ */
+#ifndef DUMUX_FC_LOCAL_ASSEMBLER_HH
+#define DUMUX_FC_LOCAL_ASSEMBLER_HH
+
+#include
+
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+#include
+
+namespace Dumux {
+
+namespace Detail {
+
+struct NoOpFunctor
+{
+ template
+ constexpr void operator()(Args&&...) const {}
+};
+
+template
+using NonVoidOrDefault_t = std::conditional_t, T, Default>;
+
+} // end namespace Detail
+
+/*!
+ * \ingroup Assembly
+ * \ingroup FaceCenteredStaggeredDiscretization
+ * \brief A base class for all local cell-centered assemblers
+ * \tparam TypeTag The TypeTag
+ * \tparam Assembler The assembler type
+ * \tparam Implementation The actual implementation
+ * \tparam implicit Specifies whether the time discretization is implicit or not (i.e. explicit)
+ */
+template
+class FaceCenteredLocalAssemblerBase : public FVLocalAssemblerBase
+{
+ using ParentType = FVLocalAssemblerBase;
+ using GridView = typename GetPropType::GridView;
+ using JacobianMatrix = GetPropType;
+ using GridVariables = GetPropType;
+ using PrimaryVariables = GetPropType;
+ using Scalar = GetPropType;
+
+ static constexpr auto numEq = GetPropType::numEq();
+
+public:
+
+ using ParentType::ParentType;
+
+ void bindLocalViews()
+ {
+ ParentType::bindLocalViews();
+ this->elemBcTypes().update(this->asImp_().problem(), this->element(), this->fvGeometry());
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix. The element residual is written into the right hand side.
+ */
+ template
+ void assembleJacobianAndResidual(JacobianMatrix& jac, ResidualVector& res, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler,
+ const CouplingFunction& maybeAssembleCouplingBlocks = CouplingFunction{})
+ {
+ static_assert(!std::decay_tasImp_().problem())>::enableInternalDirichletConstraints(),
+ "Internal Dirichlet constraints are currently not implemented for face-centered staggered models!");
+
+ this->asImp_().bindLocalViews();
+ const auto& gridGeometry = this->asImp_().problem().gridGeometry();
+ const auto eIdxGlobal = gridGeometry.elementMapper().index(this->element());
+ if (partialReassembler
+ && partialReassembler->elementColor(eIdxGlobal) == EntityColor::green)
+ {
+ const auto residual = this->asImp_().evalLocalResidual(); // forward to the internal implementation
+ for (const auto& scv : scvs(this->fvGeometry()))
+ res[scv.dofIndex()] += residual[scv.localDofIndex()];
+
+ // assemble the coupling blocks for coupled models (does nothing if not coupled)
+ maybeAssembleCouplingBlocks(residual);
+ }
+ else if (!this->elementIsGhost())
+ {
+ const auto residual = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables, partialReassembler); // forward to the internal implementation
+
+ if (this->element().partitionType() == Dune::InteriorEntity)
+ {
+ for (const auto& scv : scvs(this->fvGeometry()))
+ res[scv.dofIndex()] += residual[scv.localDofIndex()];
+ }
+ else
+ {
+ // handle residual and matrix entries for parallel runs
+ for (const auto& scv : scvs(this->fvGeometry()))
+ {
+ const auto& facet = this->element().template subEntity <1> (scv.indexInElement());
+ // make sure that the residual at border entities is consistent by adding the
+ // the contribution from the neighboring overlap element's scv
+ if (facet.partitionType() == Dune::BorderEntity)
+ res[scv.dofIndex()] += residual[scv.localDofIndex()];
+
+ // set the matrix entries of all DOFs within the overlap region (except the border DOF)
+ // to 1.0 and the residual entries to 0.0
+ else
+ {
+ const auto idx = scv.dofIndex();
+ jac[idx][idx] = 0.0;
+ for (int i = 0; i < jac[idx][idx].size(); ++i)
+ jac[idx][idx][i][i] = 1.0;
+ res[idx] = 0;
+ }
+ }
+ }
+
+ // assemble the coupling blocks for coupled models (does nothing if not coupled)
+ maybeAssembleCouplingBlocks(residual);
+ }
+ else
+ DUNE_THROW(Dune::NotImplemented, "Ghost elements not supported");
+
+
+ auto applyDirichlet = [&] (const auto& scvI,
+ const auto& dirichletValues,
+ const auto eqIdx,
+ const auto pvIdx)
+ {
+ res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
+
+ auto& row = jac[scvI.dofIndex()];
+ for (auto col = row.begin(); col != row.end(); ++col)
+ row[col.index()][eqIdx] = 0.0;
+
+ jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
+
+ // if a periodic dof has Dirichlet values also apply the same Dirichlet values to the other dof
+ if (this->asImp_().problem().gridGeometry().dofOnPeriodicBoundary(scvI.dofIndex()))
+ {
+ const auto periodicDof = this->asImp_().problem().gridGeometry().periodicallyMappedDof(scvI.dofIndex());
+ res[periodicDof][eqIdx] = this->asImp_().curSol()[periodicDof][pvIdx] - dirichletValues[pvIdx];
+
+ auto& rowP = jac[periodicDof];
+ for (auto col = rowP.begin(); col != rowP.end(); ++col)
+ row[col.index()][eqIdx] = 0.0;
+
+ rowP[periodicDof][eqIdx][pvIdx] = 1.0;
+ }
+ };
+
+ this->asImp_().enforceDirichletConstraints(applyDirichlet);
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ */
+ void assembleJacobian(JacobianMatrix& jac, GridVariables& gridVariables)
+ {
+ this->asImp_().bindLocalViews();
+ this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
+
+ auto applyDirichlet = [&] (const auto& scvI,
+ const auto& dirichletValues,
+ const auto eqIdx,
+ const auto pvIdx)
+ {
+ auto& row = jac[scvI.dofIndex()];
+ for (auto col = row.begin(); col != row.end(); ++col)
+ row[col.index()][eqIdx] = 0.0;
+
+ jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
+ };
+
+ this->asImp_().enforceDirichletConstraints(applyDirichlet);
+ }
+
+ /*!
+ * \brief Assemble the residual only
+ */
+ template
+ void assembleResidual(ResidualVector& res)
+ {
+ this->asImp_().bindLocalViews();
+ const auto residual = this->evalLocalResidual();
+
+ for (const auto& scv : scvs(this->fvGeometry()))
+ res[scv.dofIndex()] += residual[scv.localDofIndex()];
+
+ auto applyDirichlet = [&] (const auto& scvI,
+ const auto& dirichletValues,
+ const auto eqIdx,
+ const auto pvIdx)
+ {
+ res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
+ };
+
+ this->asImp_().enforceDirichletConstraints(applyDirichlet);
+ }
+
+ /*!
+ * \brief Enforce Dirichlet constraints
+ */
+ template
+ void enforceDirichletConstraints(const ApplyFunction& applyDirichlet)
+ {
+ // enforce Dirichlet boundary conditions
+ this->asImp_().evalDirichletBoundaries(applyDirichlet);
+ // take care of internal Dirichlet constraints (if enabled)
+ this->asImp_().enforceInternalDirichletConstraints(applyDirichlet);
+ }
+
+ /*!
+ * \brief Evaluates Dirichlet boundaries
+ */
+ template< typename ApplyDirichletFunctionType >
+ void evalDirichletBoundaries(ApplyDirichletFunctionType applyDirichlet)
+ {
+ // enforce Dirichlet boundaries by overwriting partial derivatives with 1 or 0
+ // and set the residual to (privar - dirichletvalue)
+ if (this->elemBcTypes().hasDirichlet())
+ {
+ for (const auto& scvf : scvfs(this->fvGeometry()))
+ {
+ if (scvf.isFrontal() && scvf.boundary())
+ {
+ const auto bcTypes = this->elemBcTypes()[scvf.localIndex()];
+ if (bcTypes.hasDirichlet())
+ {
+ const auto& scv = this->fvGeometry().scv(scvf.insideScvIdx());
+ const auto dirichletValues = this->asImp_().problem().dirichlet(this->element(), scvf);
+
+ // set the Dirichlet conditions in residual and jacobian
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ static_assert(numEq == 1, "Not yet implemented for more than one vector-valued primary variable");
+ const int pvIdx = eqIdx;
+ const int componentIdx = scv.dofAxis();
+ if (bcTypes.isDirichlet(componentIdx))
+ applyDirichlet(scv, std::array{{dirichletValues[componentIdx]}}, eqIdx, pvIdx);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /*!
+ * \brief Update the coupling context for coupled models.
+ * \note This does nothing per default (not a coupled model).
+ */
+ template
+ void maybeUpdateCouplingContext(Args&&...) {}
+
+ /*!
+ * \brief Update the additional domain derivatives for coupled models.
+ * \note This does nothing per default (not a coupled model).
+ */
+ template
+ void maybeEvalAdditionalDomainDerivatives(Args&&...) {}
+};
+
+/*!
+ * \ingroup Assembly
+ * \ingroup FaceCenteredStaggeredDiscretization
+ * \brief An assembler for Jacobian and residual contribution per element (Face-centered methods)
+ * \tparam TypeTag The TypeTag
+ * \tparam diffMethod The differentiation method to residual compute derivatives
+ * \tparam implicit Specifies whether the time discretization is implicit or not (i.e. explicit)
+ * \tparam Implementation The actual implementation, if void this class is the actual implementation
+ */
+template
+class FaceCenteredLocalAssembler;
+
+/*!
+ * \ingroup Assembly
+ * \ingroup FaceCenteredStaggeredDiscretization
+ * \brief Face-centered scheme local assembler using numeric differentiation and implicit time discretization
+ */
+template
+class FaceCenteredLocalAssembler
+: public FaceCenteredLocalAssemblerBase<
+ TypeTag, Assembler,
+ Detail::NonVoidOrDefault_t>,
+ /*implicit=*/true
+>
+{
+ using ThisType = FaceCenteredLocalAssembler;
+ using ParentType = FaceCenteredLocalAssemblerBase, true>;
+ using Scalar = GetPropType;
+ using Element = typename GetPropType::GridView::template Codim<0>::Entity;
+ using GridGeometry = GetPropType;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using GridVariables = GetPropType;
+ using JacobianMatrix = GetPropType;
+ using VolumeVariables = GetPropType;
+
+ static constexpr auto numEq = GetPropType::numEq();
+ static constexpr bool enableGridFluxVarsCache = GetPropType::GridFluxVariablesCache::cachingEnabled;
+
+public:
+
+ using LocalResidual = typename ParentType::LocalResidual;
+ using ElementResidualVector = typename LocalResidual::ElementResidualVector;
+ using ParentType::ParentType;
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ *
+ * \return The element residual at the current solution.
+ */
+ template
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler = nullptr)
+ {
+ // get some aliases for convenience
+ const auto& problem = this->asImp_().problem();
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ const auto& curSol = this->asImp_().curSol();
+ auto&& curElemVolVars = this->curElemVolVars();
+
+ // get the vector of the actual element residuals
+ const auto origResiduals = this->evalLocalResidual();
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+ // //
+ // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
+ // neighboring elements we do so by computing the derivatives of the fluxes which depend on the //
+ // actual element. In the actual element we evaluate the derivative of the entire residual. //
+ // //
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+
+ // one residual per element facet
+ const auto numElementResiduals = fvGeometry.numScv();
+
+ // create the vector storing the partial derivatives
+ ElementResidualVector partialDerivs(numElementResiduals);
+
+ const auto evalSource = [&](ElementResidualVector& residual, const SubControlVolume& scv)
+ {
+ this->localResidual().evalSource(residual, problem, element, fvGeometry, curElemVolVars, scv);
+ };
+
+ const auto evalStorage = [&](ElementResidualVector& residual, const SubControlVolume& scv)
+ {
+ this->localResidual().evalStorage(residual, problem, element, fvGeometry, this->prevElemVolVars(), curElemVolVars, scv);
+ };
+
+ const auto evalFlux = [&](ElementResidualVector& residual, const SubControlVolumeFace& scvf)
+ {
+ if (!scvf.processorBoundary())
+ this->localResidual().evalFlux(residual, problem, element, fvGeometry, curElemVolVars, this->elemBcTypes(), this->elemFluxVarsCache(), scvf);
+ };
+
+ const auto evalDerivative = [&] (const auto& scvI, const auto& scvJ)
+ {
+ // derivative w.r.t. own DOF
+ for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ {
+ partialDerivs = 0.0;
+ const auto& otherElement = fvGeometry.gridGeometry().element(scvJ.elementIndex());
+ auto otherElemSol = elementSolution(otherElement, curSol, fvGeometry.gridGeometry()); // TODO allow selective creation of elemsol (for one scv)
+ auto& curOtherVolVars = this->getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scvJ);
+ const VolumeVariables origOtherVolVars(curOtherVolVars);
+
+ auto evalResiduals = [&](Scalar priVar)
+ {
+ // update the volume variables and compute element residual
+ otherElemSol[scvJ.localDofIndex()][pvIdx] = priVar;
+ curOtherVolVars.update(otherElemSol, problem, otherElement, scvJ);
+ this->asImp_().maybeUpdateCouplingContext(scvJ, otherElemSol, pvIdx);
+
+ ElementResidualVector residual(numElementResiduals);
+ residual = 0;
+
+ evalSource(residual, scvI);
+
+ if (!this->assembler().isStationaryProblem())
+ evalStorage(residual, scvI);
+
+ for (const auto& scvf : scvfs(fvGeometry, scvI))
+ evalFlux(residual, scvf);
+
+ return residual;
+ };
+
+ // derive the residuals numerically
+ static const NumericEpsilon eps_{this->asImp_().problem().paramGroup()};
+ static const int numDiffMethod = getParamFromGroup(this->asImp_().problem().paramGroup(), "Assembly.NumericDifferenceMethod");
+ NumericDifferentiation::partialDerivative(evalResiduals, otherElemSol[scvJ.localDofIndex()][pvIdx], partialDerivs, origResiduals,
+ eps_(otherElemSol[scvJ.localDofIndex()][pvIdx], pvIdx), numDiffMethod);
+
+ const auto updateJacobian = [&]()
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ {
+ // A[i][col][eqIdx][pvIdx] is the rate of change of
+ // the residual of equation 'eqIdx' at dof 'i'
+ // depending on the primary variable 'pvIdx' at dof
+ // 'col'.
+ A[scvI.dofIndex()][scvJ.dofIndex()][eqIdx][pvIdx] += partialDerivs[scvI.localDofIndex()][eqIdx];
+ }
+ };
+
+ if (element.partitionType() == Dune::InteriorEntity)
+ updateJacobian();
+ else
+ {
+ const auto localIdxI = scvI.indexInElement();
+ const auto& facetI = element.template subEntity <1> (localIdxI);
+
+ if (facetI.partitionType() == Dune::BorderEntity)
+ updateJacobian();
+ }
+
+ // restore the original state of the scv's volume variables
+ curOtherVolVars = origOtherVolVars;
+
+ // restore the original element solution
+ otherElemSol[scvJ.localDofIndex()][pvIdx] = curSol[scvJ.dofIndex()][pvIdx];
+ this->asImp_().maybeUpdateCouplingContext(scvJ, otherElemSol, pvIdx);
+ // TODO additional dof dependencies
+ }
+ };
+
+ // calculation of the derivatives
+ for (const auto& scvI : scvs(fvGeometry))
+ {
+ // derivative w.r.t. own DOFs
+ evalDerivative(scvI, scvI);
+
+ // derivative w.r.t. other DOFs
+ const auto& otherScvIndices = fvGeometry.gridGeometry().connectivityMap()[scvI.index()];
+ for (const auto globalJ : otherScvIndices)
+ evalDerivative(scvI, fvGeometry.scv(globalJ));
+ }
+
+ // evaluate additional derivatives that might arise from the coupling
+ this->asImp_().maybeEvalAdditionalDomainDerivatives(origResiduals, A, gridVariables);
+
+ return origResiduals;
+ }
+};
+
+
+/*!
+ * \ingroup Assembly
+ * \ingroup FaceCenteredStaggeredDiscretization
+ * \brief TODO docme
+ */
+template
+class FaceCenteredLocalAssembler
+: public FaceCenteredLocalAssemblerBase<
+ TypeTag, Assembler,
+ Detail::NonVoidOrDefault_t>,
+ /*implicit=*/false
+>
+{
+ using ThisType = FaceCenteredLocalAssembler;
+ using ParentType = FaceCenteredLocalAssemblerBase, false>;
+ using Scalar = GetPropType;
+ using Element = typename GetPropType::GridView::template Codim<0>::Entity;
+ using GridVariables = GetPropType;
+ using JacobianMatrix = GetPropType;
+ using VolumeVariables = GetPropType;
+
+ enum { numEq = GetPropType::numEq() };
+
+public:
+ using LocalResidual = typename ParentType::LocalResidual;
+ using ElementResidualVector = typename LocalResidual::ElementResidualVector;
+ using ParentType::ParentType;
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ *
+ * \return The element residual at the current solution.
+ */
+ template
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler = nullptr)
+ {
+ if (partialReassembler)
+ DUNE_THROW(Dune::NotImplemented, "partial reassembly for explicit time discretization");
+
+ // get some aliases for convenience
+ const auto& problem = this->asImp_().problem();
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ const auto& curSol = this->asImp_().curSol();
+ auto&& curElemVolVars = this->curElemVolVars();
+
+ // get the vector of the actual element residuals
+ const auto origResiduals = this->evalLocalResidual();
+ const auto origStorageResiduals = this->evalLocalStorageResidual();
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+ // //
+ // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
+ // neighboring elements we do so by computing the derivatives of the fluxes which depend on the //
+ // actual element. In the actual element we evaluate the derivative of the entire residual. //
+ // //
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+
+ // create the element solution
+ auto elemSol = elementSolution(element, curSol, fvGeometry.gridGeometry());
+
+ // create the vector storing the partial derivatives
+ ElementResidualVector partialDerivs(fvGeometry.numScv());
+
+ // calculation of the derivatives
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ // dof index and corresponding actual pri vars
+ const auto dofIdx = scv.dofIndex();
+ auto& curVolVars = this->getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
+ const VolumeVariables origVolVars(curVolVars);
+
+ // calculate derivatives w.r.t to the privars at the dof at hand
+ for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ {
+ partialDerivs = 0.0;
+
+ auto evalStorage = [&](Scalar priVar)
+ {
+ // auto partialDerivsTmp = partialDerivs;
+ elemSol[scv.localDofIndex()][pvIdx] = priVar;
+ curVolVars.update(elemSol, problem, element, scv);
+ return this->evalLocalStorageResidual();
+ };
+
+ // derive the residuals numerically
+ static const NumericEpsilon eps_{problem.paramGroup()};
+ static const int numDiffMethod = getParamFromGroup(problem.paramGroup(), "Assembly.NumericDifferenceMethod");
+ NumericDifferentiation::partialDerivative(evalStorage, elemSol[scv.localDofIndex()][pvIdx], partialDerivs, origStorageResiduals,
+ eps_(elemSol[scv.localDofIndex()][pvIdx], pvIdx), numDiffMethod);
+
+ // update the global stiffness matrix with the current partial derivatives
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ {
+ // A[i][col][eqIdx][pvIdx] is the rate of change of
+ // the residual of equation 'eqIdx' at dof 'i'
+ // depending on the primary variable 'pvIdx' at dof
+ // 'col'.
+ A[dofIdx][dofIdx][eqIdx][pvIdx] += partialDerivs[scv.localDofIndex()][eqIdx];
+ }
+
+ // restore the original state of the scv's volume variables
+ curVolVars = origVolVars;
+
+ // restore the original element solution
+ elemSol[scv.localDofIndex()][pvIdx] = curSol[scv.dofIndex()][pvIdx];
+ }
+ }
+ return origResiduals;
+ }
+};
+
+/*!
+ * \ingroup Assembly
+ * \ingroup FaceCenteredStaggeredDiscretization
+ * \brief TODO docme
+ */
+template
+class FaceCenteredLocalAssembler
+: public FaceCenteredLocalAssemblerBase<
+ TypeTag, Assembler,
+ FaceCenteredLocalAssembler,
+ /*implicit=*/true
+>
+{
+ using ThisType = FaceCenteredLocalAssembler;
+ using ParentType = FaceCenteredLocalAssemblerBase, true>;
+ using JacobianMatrix = GetPropType;
+ using GridVariables = GetPropType;
+ using VolumeVariables = GetPropType;
+
+ enum { numEq = GetPropType::numEq() };
+
+public:
+ using LocalResidual = typename ParentType::LocalResidual;
+ using ElementResidualVector = typename LocalResidual::ElementResidualVector;
+ using ParentType::ParentType;
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ *
+ * \return The element residual at the current solution.
+ */
+ template
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler = nullptr)
+ {
+ if (partialReassembler)
+ DUNE_THROW(Dune::NotImplemented, "partial reassembly for analytic differentiation");
+
+ // get some aliases for convenience
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ const auto& problem = this->asImp_().problem();
+ const auto& curElemVolVars = this->curElemVolVars();
+ const auto& elemFluxVarsCache = this->elemFluxVarsCache();
+
+ // get the vecor of the actual element residuals
+ const auto origResiduals = this->evalLocalResidual();
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+ // //
+ // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
+ // neighboring elements we do so by computing the derivatives of the fluxes which depend on the //
+ // actual element. In the actual element we evaluate the derivative of the entire residual. //
+ // //
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+
+ // calculation of the source and storage derivatives
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ // dof index and corresponding actual pri vars
+ const auto dofIdx = scv.dofIndex();
+ const auto& volVars = curElemVolVars[scv];
+
+ // derivative of this scv residual w.r.t the d.o.f. of the same scv (because of mass lumping)
+ // only if the problem is instationary we add derivative of storage term
+ // TODO if e.g. porosity depends on all dofs in the element, we would have off-diagonal matrix entries!?
+ if (!this->assembler().isStationaryProblem())
+ this->localResidual().addStorageDerivatives(A[dofIdx][dofIdx],
+ problem,
+ element,
+ fvGeometry,
+ volVars,
+ scv);
+
+ // derivative of this scv residual w.r.t the d.o.f. of the same scv (because of mass lumping)
+ // add source term derivatives
+ this->localResidual().addSourceDerivatives(A[dofIdx][dofIdx],
+ problem,
+ element,
+ fvGeometry,
+ volVars,
+ scv);
+ }
+
+ // localJacobian[scvIdx][otherScvIdx][eqIdx][priVarIdx] of the fluxes
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ if (!scvf.boundary())
+ {
+ // add flux term derivatives
+ this->localResidual().addFluxDerivatives(A,
+ problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ scvf);
+ }
+
+ // the boundary gets special treatment to simplify
+ // for the user
+ else
+ {
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ if (this->elemBcTypes()[insideScv.localDofIndex()].hasNeumann())
+ {
+ // add flux term derivatives
+ this->localResidual().addRobinFluxDerivatives(A[insideScv.dofIndex()],
+ problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ scvf);
+ }
+ }
+ }
+
+ return origResiduals;
+ }
+};
+
+/*!
+ * \ingroup Assembly
+ * \ingroup FaceCenteredStaggeredDiscretization
+ * \brief TODO docme
+ */
+template
+class FaceCenteredLocalAssembler
+: public FaceCenteredLocalAssemblerBase<
+ TypeTag, Assembler,
+ FaceCenteredLocalAssembler,
+ /*implicit=*/false
+>
+{
+ using ThisType = FaceCenteredLocalAssembler;
+ using ParentType = FaceCenteredLocalAssemblerBase, false>;
+ using JacobianMatrix = GetPropType;
+ using GridVariables = GetPropType;
+ using VolumeVariables = GetPropType;
+
+ enum { numEq = GetPropType::numEq() };
+
+public:
+ using LocalResidual = typename ParentType::LocalResidual;
+ using ElementResidualVector = typename LocalResidual::ElementResidualVector;
+ using ParentType::ParentType;
+
+ /*!
+ * \brief Computes the derivatives with respect to the given element and adds them
+ * to the global matrix.
+ *
+ * \return The element residual at the current solution.
+ */
+ template
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler = nullptr)
+ {
+ if (partialReassembler)
+ DUNE_THROW(Dune::NotImplemented, "partial reassembly for explicit time discretization");
+
+ // get some aliases for convenience
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ const auto& problem = this->asImp_().problem();
+ const auto& curElemVolVars = this->curElemVolVars();
+
+ // get the vector of the actual element residuals
+ const auto origResiduals = this->evalLocalResidual();
+
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+ // //
+ // Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
+ // neighboring elements we do so by computing the derivatives of the fluxes which depend on the //
+ // actual element. In the actual element we evaluate the derivative of the entire residual. //
+ // //
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+
+ // calculation of the source and storage derivatives
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ // dof index and corresponding actual pri vars
+ const auto dofIdx = scv.dofIndex();
+ const auto& volVars = curElemVolVars[scv];
+
+ // derivative of this scv residual w.r.t the d.o.f. of the same scv (because of mass lumping)
+ // only if the problem is instationary we add derivative of storage term
+ this->localResidual().addStorageDerivatives(A[dofIdx][dofIdx],
+ problem,
+ element,
+ fvGeometry,
+ volVars,
+ scv);
+ }
+
+ return origResiduals;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/fvlocalassemblerbase.hh b/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/fvlocalassemblerbase.hh
new file mode 100644
index 0000000..0307d07
--- /dev/null
+++ b/benchmarks/rotating-cylinders/dumux/dumuxModule/dumux/assembly/fvlocalassemblerbase.hh
@@ -0,0 +1,330 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+//
+// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
+// SPDX-License-Identifier: GPL-3.0-or-later
+//
+/*!
+ * \file
+ * \ingroup Assembly
+ * \copydoc Dumux::FVLocalAssemblerBase
+ */
+#ifndef DUMUX_FV_LOCAL_ASSEMBLER_BASE_HH
+#define DUMUX_FV_LOCAL_ASSEMBLER_BASE_HH
+
+#include
+#include // for GhostEntity
+#include
+
+#include
+#include
+#include
+#include
+
+namespace Dumux {
+
+/*!
+ * \ingroup Assembly
+ * \brief A base class for all local assemblers
+ * \tparam TypeTag The TypeTag
+ * \tparam Assembler The assembler type
+ * \tparam Implementation The assembler implementation
+ * \tparam useImplicitAssembly Specifies whether the time discretization is implicit or not not (i.e. explicit)
+ */
+template
+class FVLocalAssemblerBase
+{
+ using Problem = GetPropType;
+ using GridView = typename GetPropType