diff --git a/dumux/assembly/cvfelocalassembler.hh b/dumux/assembly/cvfelocalassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..44ef08fcc45bd3795222f5b316e066cbcefac422
--- /dev/null
+++ b/dumux/assembly/cvfelocalassembler.hh
@@ -0,0 +1,778 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Assembly
+ * \ingroup CVFEDiscretization
+ * \brief An assembler for Jacobian and residual contribution per element (CVFE methods)
+ */
+#ifndef DUMUX_CVFE_LOCAL_ASSEMBLER_HH
+#define DUMUX_CVFE_LOCAL_ASSEMBLER_HH
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/hybridutilities.hh>
+#include <dune/common/reservedvector.hh>
+#include <dune/grid/common/gridenums.hh>
+#include <dune/istl/matrixindexset.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/numericdifferentiation.hh>
+
+#include <dumux/assembly/numericepsilon.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/assembly/fvlocalassemblerbase.hh>
+#include <dumux/assembly/partialreassembler.hh>
+#include <dumux/assembly/entitycolor.hh>
+
+#include <dumux/discretization/cvfe/elementsolution.hh>
+
+#include "volvardeflectionhelper_.hh"
+
+namespace Dumux {
+
+#ifndef DOXYGEN
+namespace Detail::CVFE {
+
+struct NoOperator
+{
+ template<class... Args>
+ constexpr void operator()(Args&&...) const {}
+};
+
+template<class X, class Y>
+using Impl = std::conditional_t<!std::is_same_v<X, void>, X, Y>;
+
+} // end namespace Detail
+#endif // DOXYGEN
+
+/*!
+ * \ingroup Assembly
+ * \ingroup CVFEDiscretization
+ * \brief A base class for all local CVFE 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 not (i.e. explicit)
+ */
+template<class TypeTag, class Assembler, class Implementation, bool implicit>
+class CVFELocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler, Implementation, implicit>
+{
+ using ParentType = FVLocalAssemblerBase<TypeTag, Assembler, Implementation, implicit>;
+ using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using ElementVolumeVariables = typename GridVariables::GridVolumeVariables::LocalView;
+ using SolutionVector = typename Assembler::SolutionVector;
+
+ static constexpr int numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq();
+ static constexpr int dim = GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension;
+
+public:
+
+ using ParentType::ParentType;
+
+ void bindLocalViews()
+ {
+ std::cout << "hi" << std::endl;
+ 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 <class ResidualVector, class PartialReassembler = DefaultPartialReassembler, class CouplingFunction = Detail::CVFE::NoOperator>
+ void assembleJacobianAndResidual(JacobianMatrix& jac, ResidualVector& res, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler = nullptr,
+ const CouplingFunction& maybeAssembleCouplingBlocks = {})
+ {
+ this->asImp_().bindLocalViews();
+ const auto eIdxGlobal = this->asImp_().problem().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
+ 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
+ {
+ // Treatment of ghost elements
+ assert(this->elementIsGhost());
+
+ // handle dofs per codimension
+ const auto& gridGeometry = this->asImp_().problem().gridGeometry();
+ Dune::Hybrid::forEach(std::make_integer_sequence<int, dim>{}, [&](auto d)
+ {
+ constexpr int codim = dim - d;
+ const auto& localCoeffs = gridGeometry.feCache().get(this->element().type()).localCoefficients();
+ for (int idx = 0; idx < localCoeffs.size(); ++idx)
+ {
+ const auto& localKey = localCoeffs.localKey(idx);
+
+ // skip if we are not handling this codim right now
+ if (localKey.codim() != codim)
+ continue;
+
+ // do not change the non-ghost entities
+ auto entity = this->element().template subEntity<codim>(localKey.subEntity());
+ if (entity.partitionType() == Dune::InteriorEntity || entity.partitionType() == Dune::BorderEntity)
+ continue;
+
+ // WARNING: this only works if the mapping from codim+subEntity to
+ // global dofIndex is unique (on dof per entity of this codim).
+ // For more general mappings, we should use a proper local-global mapping here.
+ // For example through dune-functions.
+ const auto dofIndex = gridGeometry.dofMapper().index(entity);
+
+ // this might be a vector-valued dof
+ using BlockType = typename JacobianMatrix::block_type;
+ BlockType &J = jac[dofIndex][dofIndex];
+ for (int j = 0; j < BlockType::rows; ++j)
+ J[j][j] = 1.0;
+
+ // set residual for the ghost dof
+ res[dofIndex] = 0;
+ }
+ });
+ }
+
+ 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->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
+ const auto end = jac[periodicDof].end();
+ for (auto it = jac[periodicDof].begin(); it != end; ++it)
+ (*it) = periodicDof != it.index() ? 0.0 : 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 <class ResidualVector>
+ 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);
+ }
+
+ //! Enforce Dirichlet constraints
+ template<typename ApplyFunction>
+ 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& scvI : scvs(this->fvGeometry()))
+ {
+ const auto bcTypes = this->elemBcTypes().get(this->fvGeometry(), scvI);
+ if (bcTypes.hasDirichlet())
+ {
+ const auto dirichletValues = this->asImp_().problem().dirichlet(this->element(), scvI);
+
+ // set the Dirichlet conditions in residual and jacobian
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (bcTypes.isDirichlet(eqIdx))
+ {
+ const auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
+ assert(0 <= pvIdx && pvIdx < numEq);
+ applyDirichlet(scvI, dirichletValues, eqIdx, pvIdx);
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /*!
+ * \brief Update the coupling context for coupled models.
+ * \note This does nothing per default (not a coupled model).
+ */
+ template<class... Args>
+ void maybeUpdateCouplingContext(Args&&...) {}
+
+ /*!
+ * \brief Update the additional domain derivatives for coupled models.
+ * \note This does nothing per default (not a coupled model).
+ */
+ template<class... Args>
+ void maybeEvalAdditionalDomainDerivatives(Args&&...) {}
+};
+
+/*!
+ * \ingroup Assembly
+ * \ingroup CVFEDiscretization
+ * \brief An assembler for Jacobian and residual contribution per element (CVFE 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 not (i.e. explicit)
+ * \tparam Implementation via CRTP
+ */
+template<class TypeTag, class Assembler, DiffMethod diffMethod = DiffMethod::numeric, bool implicit = true, class Implementation = void>
+class CVFELocalAssembler;
+
+/*!
+ * \ingroup Assembly
+ * \ingroup CVFEDiscretization
+ * \brief Control volume finite element local assembler using numeric differentiation and implicit time discretization
+ */
+template<class TypeTag, class Assembler, class Implementation>
+class CVFELocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true, Implementation>
+: public CVFELocalAssemblerBase<TypeTag, Assembler,
+ Detail::CVFE::Impl<Implementation, CVFELocalAssembler<TypeTag, Assembler, DiffMethod::numeric, true, Implementation>>,
+ true>
+{
+ using ThisType = CVFELocalAssembler<TypeTag, Assembler, DiffMethod::numeric, true, Implementation>;
+ using ParentType = CVFELocalAssemblerBase<TypeTag, Assembler, Detail::CVFE::Impl<Implementation, ThisType>, true>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
+
+ static constexpr int numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq();
+ static constexpr int dim = GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension;
+
+ static constexpr bool enableGridFluxVarsCache
+ = GridVariables::GridFluxVariablesCache::cachingEnabled;
+ static constexpr bool solutionDependentFluxVarsCache
+ = GridVariables::GridFluxVariablesCache::FluxVariablesCache::isSolDependent;
+
+public:
+
+ using LocalResidual = GetPropType<TypeTag, Properties::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 <class PartialReassembler = DefaultPartialReassembler>
+ ElementResidualVector assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables,
+ const PartialReassembler* partialReassembler = nullptr)
+ {
+ // get some aliases for convenience
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ const auto& curSol = this->asImp_().curSol();
+
+ auto&& curElemVolVars = this->curElemVolVars();
+ auto&& elemFluxVarsCache = this->elemFluxVarsCache();
+
+ // 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. //
+ // //
+ //////////////////////////////////////////////////////////////////////////////////////////////////
+
+ // if all volvars in the stencil have to be updated or if it's enough to only update the
+ // volVars for the scv whose associated dof has been deflected
+ static const bool updateAllVolVars = getParamFromGroup<bool>(
+ this->asImp_().problem().paramGroup(), "Assembly.BoxVolVarsDependOnAllElementDofs", false
+ );
+
+ // create the element solution
+ auto elemSol = elementSolution(element, curSol, fvGeometry.gridGeometry());
+
+ // create the vector storing the partial derivatives
+ ElementResidualVector partialDerivs(fvGeometry.numScv());
+
+ Detail::VolVarsDeflectionHelper deflectionHelper(
+ [&] (const auto& scv) -> VolumeVariables& {
+ return this->getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
+ },
+ fvGeometry,
+ updateAllVolVars
+ );
+
+ // calculation of the derivatives
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ // dof index and corresponding actual pri vars
+ const auto dofIdx = scv.dofIndex();
+ deflectionHelper.setCurrent(scv);
+
+ // calculate derivatives w.r.t to the privars at the dof at hand
+ for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ {
+ partialDerivs = 0.0;
+
+ auto evalResiduals = [&](Scalar priVar)
+ {
+ // update the volume variables and compute element residual
+ elemSol[scv.localDofIndex()][pvIdx] = priVar;
+ deflectionHelper.deflect(elemSol, scv, this->asImp_().problem());
+ if constexpr (solutionDependentFluxVarsCache)
+ {
+ elemFluxVarsCache.update(element, fvGeometry, curElemVolVars);
+ if constexpr (enableGridFluxVarsCache)
+ gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
+ }
+ this->asImp_().maybeUpdateCouplingContext(scv, elemSol, pvIdx);
+ return this->evalLocalResidual();
+ };
+
+ // derive the residuals numerically
+ static const NumericEpsilon<Scalar, numEq> eps_{this->asImp_().problem().paramGroup()};
+ static const int numDiffMethod = getParamFromGroup<int>(this->asImp_().problem().paramGroup(), "Assembly.NumericDifferenceMethod");
+ NumericDifferentiation::partialDerivative(evalResiduals, elemSol[scv.localDofIndex()][pvIdx], partialDerivs, origResiduals,
+ eps_(elemSol[scv.localDofIndex()][pvIdx], pvIdx), numDiffMethod);
+
+ // update the global stiffness matrix with the current partial derivatives
+ for (const auto& scvJ : scvs(fvGeometry))
+ {
+ // don't add derivatives for green dofs
+ if (!partialReassembler
+ || partialReassembler->dofColor(scvJ.dofIndex()) != EntityColor::green)
+ {
+ 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[scvJ.dofIndex()][dofIdx][eqIdx][pvIdx] += partialDerivs[scvJ.localDofIndex()][eqIdx];
+ }
+ }
+ }
+
+ // restore the original state of the scv's volume variables
+ deflectionHelper.restore(scv);
+
+ // restore the original element solution
+ elemSol[scv.localDofIndex()][pvIdx] = curSol[scv.dofIndex()][pvIdx];
+ this->asImp_().maybeUpdateCouplingContext(scv, elemSol, pvIdx);
+ }
+ }
+
+ // restore original state of the flux vars cache in case of global caching.
+ // In the case of local caching this is obsolete because the elemFluxVarsCache used here goes out of scope after this.
+ if constexpr (enableGridFluxVarsCache)
+ gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
+
+ // evaluate additional derivatives that might arise from the coupling (no-op if not coupled)
+ this->asImp_().maybeEvalAdditionalDomainDerivatives(origResiduals, A, gridVariables);
+
+ return origResiduals;
+ }
+
+}; // implicit CVFEAssembler with numeric Jacobian
+
+/*!
+ * \ingroup Assembly
+ * \ingroup CVFEDiscretization
+ * \brief Control volume finite element local assembler using numeric differentiation and explicit time discretization
+ */
+template<class TypeTag, class Assembler, class Implementation>
+class CVFELocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false, Implementation>
+: public CVFELocalAssemblerBase<TypeTag, Assembler,
+ Detail::CVFE::Impl<Implementation, CVFELocalAssembler<TypeTag, Assembler, DiffMethod::numeric, false, Implementation>>,
+ false>
+{
+ using ThisType = CVFELocalAssembler<TypeTag, Assembler, DiffMethod::numeric, false, Implementation>;
+ using ParentType = CVFELocalAssemblerBase<TypeTag, Assembler, Detail::CVFE::Impl<Implementation, ThisType>, false>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using VolumeVariables = GetPropType<TypeTag, Properties::VolumeVariables>;
+ using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
+
+ static constexpr int numEq = GetPropType<TypeTag, Properties::ModelTraits>::numEq();
+ static constexpr int dim = GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension;
+
+public:
+
+ using LocalResidual = GetPropType<TypeTag, Properties::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 <class PartialReassembler = DefaultPartialReassembler>
+ 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& curSol = this->asImp_().curSol();
+ auto&& curElemVolVars = this->curElemVolVars();
+
+ // get the vecor 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, this->asImp_().problem(), element, scv);
+ return this->evalLocalStorageResidual();
+ };
+
+ // derive the residuals numerically
+ static const NumericEpsilon<Scalar, numEq> eps_{this->asImp_().problem().paramGroup()};
+ static const int numDiffMethod = getParamFromGroup<int>(this->asImp_().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];
+ // TODO additional dof dependencies
+ }
+ }
+
+ return origResiduals;
+ }
+}; // explicit CVFEAssembler with numeric Jacobian
+
+/*!
+ * \ingroup Assembly
+ * \ingroup CVFEDiscretization
+ * \brief Control volume finite element local assembler using analytic differentiation and implicit time discretization
+ */
+template<class TypeTag, class Assembler, class Implementation>
+class CVFELocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/true, Implementation>
+: public CVFELocalAssemblerBase<TypeTag, Assembler,
+ Detail::CVFE::Impl<Implementation, CVFELocalAssembler<TypeTag, Assembler, DiffMethod::analytic, true, Implementation>>,
+ true>
+{
+ using ThisType = CVFELocalAssembler<TypeTag, Assembler, DiffMethod::analytic, true, Implementation>;
+ using ParentType = CVFELocalAssemblerBase<TypeTag, Assembler, Detail::CVFE::Impl<Implementation, ThisType>, true>;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
+
+public:
+
+ using LocalResidual = GetPropType<TypeTag, Properties::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 <class PartialReassembler = DefaultPartialReassembler>
+ 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().get(fvGeometry, insideScv).hasNeumann())
+ {
+ // add flux term derivatives
+ this->localResidual().addRobinFluxDerivatives(A[insideScv.dofIndex()],
+ problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ scvf);
+ }
+ }
+ }
+
+ return origResiduals;
+ }
+
+}; // implicit CVFEAssembler with analytic Jacobian
+
+/*!
+ * \ingroup Assembly
+ * \ingroup CVFEDiscretization
+ * \brief Control volume finite element local assembler using analytic differentiation and explicit time discretization
+ */
+template<class TypeTag, class Assembler, class Implementation>
+class CVFELocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/false, Implementation>
+: public CVFELocalAssemblerBase<TypeTag, Assembler,
+ Detail::CVFE::Impl<Implementation, CVFELocalAssembler<TypeTag, Assembler, DiffMethod::analytic, false, Implementation>>,
+ false>
+{
+ using ThisType = CVFELocalAssembler<TypeTag, Assembler, DiffMethod::analytic, false, Implementation>;
+ using ParentType = CVFELocalAssemblerBase<TypeTag, Assembler, Detail::CVFE::Impl<Implementation, ThisType>, false>;
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ using JacobianMatrix = GetPropType<TypeTag, Properties::JacobianMatrix>;
+
+public:
+
+ using LocalResidual = GetPropType<TypeTag, Properties::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 <class PartialReassembler = DefaultPartialReassembler>
+ 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 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
+ this->localResidual().addStorageDerivatives(A[dofIdx][dofIdx],
+ problem,
+ element,
+ fvGeometry,
+ volVars,
+ scv);
+ }
+
+ return origResiduals;
+ }
+
+}; // explicit CVFEAssembler with analytic Jacobian
+
+} // end namespace Dumux
+
+#endif