+class GridVariablesFacade
+{
+ using P = typename ElemVolVars::GridVolumeVariables::Problem;
+public:
+ using GridVolumeVariables = typename ElemVolVars::GridVolumeVariables;
+ using GridFluxVariablesCache = typename ElemFluxVarsCache::GridFluxVariablesCache;
+ using GridGeometry = typename ProblemTraits::GridGeometry;
+};
+
+// implements necessary interfaces for compatibility
+template
+class ElemVariablesFacade
+{
+ const ElemVolVars* elemVolVars_;
+ const ElemFluxVarsCache* elemFluxVarsCache_;
+
+public:
+
+ using GridVariables = GridVariablesFacade;
+
+ ElemVariablesFacade(const ElemVolVars* elemVolVars,
+ const ElemFluxVarsCache* elemFluxVarsCache)
+ : elemVolVars_(elemVolVars)
+ , elemFluxVarsCache_(elemFluxVarsCache)
+ {}
+
+ const ElemVolVars& elemVolVars() const { return *elemVolVars_; }
+ const ElemFluxVarsCache& elemFluxVarsCache() const { return *elemFluxVarsCache_; }
+};
+
+// helper function to construct an element variables facade from given local views
+template
+ElemVariablesFacade
+makeElemVariablesFacade(const ElemVolVars& elemVolVars,
+ const ElemFluxVarsCache& elemFluxVarsCache)
+{ return {&elemVolVars, &elemFluxVarsCache}; }
+
+// get the Dirichlet values for a boundary entity from a user problem
+template
+decltype(auto) getDirichletValues(const Problem& problem,
+ const Element& element,
+ const BoundaryEntity& boundaryEntity,
+ const ElemSol& elemSol)
+{
+ // given element solution models the experimental elemsol concept
+ if constexpr (Dune::models())
+ {
+ if constexpr (ProblemTraits::hasTransientDirichletInterface)
+ return problem.dirichlet(element, boundaryEntity, elemSol.timeLevel());
+ else
+ return problem.dirichlet(element, boundaryEntity);
+ }
+
+ // current default style (elemSolState = element solution, not passed to user interface)
+ else
+ return problem.dirichlet(element, boundaryEntity);
+}
+
+} // end namespace CompatibilityHelpers
+} // end namespace Dumux::Experimental
+
+#endif
diff --git a/dumux/assembly/fv/CMakeLists.txt b/dumux/assembly/fv/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..890727a9263f546be0e2cc0938586ff7981885e1
--- /dev/null
+++ b/dumux/assembly/fv/CMakeLists.txt
@@ -0,0 +1,6 @@
+install(FILES
+boxlocalassembler.hh
+cclocalassembler.hh
+localoperator.hh
+operators.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/assembly/fv)
diff --git a/dumux/assembly/fv/boxlocalassembler.hh b/dumux/assembly/fv/boxlocalassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..633323c54688055373667d35405e2a3992600d3e
--- /dev/null
+++ b/dumux/assembly/fv/boxlocalassembler.hh
@@ -0,0 +1,500 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief An assembler for Jacobian and residual contribution
+ * per element for the box scheme.
+ */
+#ifndef DUMUX_BOX_LOCAL_ASSEMBLER_HH
+#define DUMUX_BOX_LOCAL_ASSEMBLER_HH
+
+#include
+#include
+
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+
+namespace Dumux::Experimental {
+
+/*!
+ * \ingroup Assembly
+ * \brief An assembler for Jacobian and residual contribution
+ * per element for the box scheme.
+ * \tparam LO The element-local operator
+ */
+template
+class BoxLocalAssembler
+{
+ using LocalContext = typename LO::LocalContext;
+ using ElementVariables = typename LocalContext::ElementVariables;
+ using ElementGridGeometry = typename LocalContext::ElementGridGeometry;
+
+ using GG = typename ElementGridGeometry::GridGeometry;
+ using GV = typename ElementVariables::GridVariables;
+
+ using FVElementGeometry = typename GG::LocalView;
+ using SubControlVolume = typename GG::SubControlVolume;
+ using Element = typename GG::GridView::template Codim<0>::Entity;
+
+ using PrimaryVariables = typename GV::PrimaryVariables;
+ using Scalar = typename GV::Scalar;
+
+ static constexpr int numEq = NumEqVectorTraits::numEq;
+
+public:
+
+ //! the parameters of a stage in time integration
+ using StageParams = MultiStageParams;
+
+ //! export the grid variables type on which to operate
+ using GridVariables = GV;
+
+ //! export the underlying local operator
+ using LocalOperator = LO;
+
+ //! export the vector storing the residuals of all dofs of the element
+ using ElementResidualVector = typename LO::ElementResidualVector;
+
+ /*!
+ * \brief Constructor for stationary problems.
+ */
+ explicit BoxLocalAssembler(const Element& element,
+ GridVariables& gridVariables,
+ DiffMethod dm = DiffMethod::numeric)
+ : diffMethod_(dm)
+ , gridVariables_(gridVariables)
+ , fvGeometry_(localView(gridVariables.gridGeometry()))
+ , elemVariables_(localView(gridVariables))
+ , prevElemVariables_()
+ , elementIsGhost_(element.partitionType() == Dune::GhostEntity)
+ , stageParams_(nullptr)
+ {
+ if (diffMethod_ != DiffMethod::numeric)
+ DUNE_THROW(Dune::NotImplemented, "Provided differentiation method");
+
+ fvGeometry_.bind(element);
+ elemVariables_.bind(element, fvGeometry_);
+ }
+
+ /*!
+ * \brief Constructor for instationary problems.
+ * \note Using this constructor, we assemble one stage within
+ * a time integration step using multi-stage methods.
+ */
+ explicit BoxLocalAssembler(const Element& element,
+ const std::vector& prevGridVars,
+ GridVariables& gridVariables,
+ std::shared_ptr stageParams,
+ DiffMethod dm = DiffMethod::numeric)
+ : diffMethod_(dm)
+ , gridVariables_(gridVariables)
+ , fvGeometry_(localView(gridVariables.gridGeometry()))
+ , elemVariables_(localView(gridVariables))
+ , prevElemVariables_()
+ , elementIsGhost_(element.partitionType() == Dune::GhostEntity)
+ , stageParams_(stageParams)
+ {
+ if (diffMethod_ != DiffMethod::numeric)
+ DUNE_THROW(Dune::NotImplemented, "Provided differentiation method");
+
+ if (prevGridVars.size() != stageParams->size() - 1)
+ DUNE_THROW(Dune::InvalidStateException, "Grid Variables for all stages needed");
+
+ fvGeometry_.bind(element);
+ for (const auto& gridVars : prevGridVars)
+ {
+ prevElemVariables_.emplace_back(localView(gridVars));
+ prevElemVariables_.back().bind(element, fvGeometry_);
+ }
+
+ elemVariables_.bind(element, 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,
+ const PartialReassembler* partialReassembler = nullptr)
+ {
+ const auto& element = fvGeometry_.element();
+ const auto eIdxGlobal = fvGeometry_.gridGeometry().elementMapper().index(element);
+
+ if (partialReassembler && partialReassembler->elementColor(eIdxGlobal) == EntityColor::green)
+ {
+ const auto residual = evalLocalResidual();
+ for (const auto& scv : scvs(fvGeometry_))
+ res[scv.dofIndex()] += residual[scv.localDofIndex()];
+ }
+ else if (!elementIsGhost_)
+ {
+ const auto residual = assembleJacobianAndResidual_(jac, partialReassembler);
+ for (const auto& scv : scvs(fvGeometry_))
+ res[scv.dofIndex()] += residual[scv.localDofIndex()];
+ }
+ else
+ {
+ static constexpr int dim = GG::GridView::dimension;
+
+ for (const auto& scv : scvs(fvGeometry_))
+ {
+ const auto vIdxLocal = scv.indexInElement();
+ const auto& v = fvGeometry_.element().template subEntity(vIdxLocal);
+
+ // do not change the non-ghost vertices
+ if (v.partitionType() == Dune::InteriorEntity ||
+ v.partitionType() == Dune::BorderEntity)
+ continue;
+
+ const auto dofIdx = scv.dofIndex();
+ res[dofIdx] = 0;
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+ jac[dofIdx][dofIdx][pvIdx][pvIdx] = 1.0;
+ }
+ }
+
+ auto applyDirichlet = [&] (const auto& scvI,
+ const auto& dirichletValues,
+ const auto eqIdx,
+ const auto pvIdx)
+ {
+ const auto& priVars = elemVariables_.elemVolVars()[scvI].priVars();
+ res[scvI.dofIndex()][eqIdx] = 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 (fvGeometry_.gridGeometry().dofOnPeriodicBoundary(scvI.dofIndex()))
+ {
+ const auto periodicDof = fvGeometry_.gridGeometry().periodicallyMappedDof(scvI.dofIndex());
+ res[periodicDof][eqIdx] = 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;
+ }
+ };
+
+ enforceDirichletConstraints(applyDirichlet);
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the dofs of the given
+ * element and adds them to the given jacobian matrix.
+ */
+ template
+ void assembleJacobian(JacobianMatrix& jac)
+ {
+ assembleJacobianAndResidual_(jac);
+
+ 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;
+ };
+
+ enforceDirichletConstraints(applyDirichlet);
+ }
+
+ //! Assemble the residual only
+ template
+ void assembleResidual(ResidualVector& res)
+ {
+ const auto residual = evalLocalResidual();
+ for (const auto& scv : scvs(fvGeometry_))
+ res[scv.dofIndex()] += residual[scv.localDofIndex()];
+
+ auto applyDirichlet = [&] (const auto& scvI,
+ const auto& dirichletValues,
+ const auto eqIdx,
+ const auto pvIdx)
+ {
+ const auto& priVars = elemVariables_.elemVolVars()[scvI].priVars();
+ res[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
+ };
+
+ enforceDirichletConstraints(applyDirichlet);
+ }
+
+ //! Enforce Dirichlet constraints
+ template
+ void enforceDirichletConstraints(const ApplyFunction& applyDirichlet)
+ {
+ // enforce Dirichlet boundary conditions
+ evalDirichletBoundaries(applyDirichlet);
+ // TODO: take care of internal Dirichlet constraints (if enabled)
+ // enforceInternalDirichletConstraints(applyDirichlet);
+ }
+
+ //! Evaluates Dirichlet boundaries
+ template< typename ApplyDirichletFunctionType >
+ void evalDirichletBoundaries(ApplyDirichletFunctionType applyDirichlet)
+ {
+ const auto& element = fvGeometry_.element();
+ const auto& problem = gridVariables_.gridVolVars().problem();
+
+ for (const auto& scvI : scvs(fvGeometry_))
+ {
+ const auto bcTypes = problem.boundaryTypes(element, scvI);
+ if (bcTypes.hasDirichlet())
+ {
+ const auto dirichletValues = getDirichletValues_(problem, element, scvI,
+ gridVariables_.timeLevel());
+
+ // 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);
+ }
+ }
+ }
+ }
+ }
+
+ //! Evaluate the complete local residual for the current element.
+ ElementResidualVector evalLocalResidual() const
+ {
+ if (isStationary())
+ {
+ const auto context = makeLocalContext(fvGeometry_, elemVariables_);
+ LocalOperator localOperator(context);
+ return elementIsGhost_ ? localOperator.getEmptyResidual()
+ : localOperator.evalFluxesAndSources();
+ }
+ else
+ {
+ ElementResidualVector residual(fvGeometry_.numScv());
+ residual = 0.0;
+
+ if (!elementIsGhost_)
+ {
+ // add the terms associated with previous stages
+ for (std::size_t k = 0; k < stageParams_->size()-1; ++k)
+ addStageTerms_(residual, k,
+ makeLocalContext(fvGeometry_, prevElemVariables_[k]));
+
+ // add the terms of the current stage
+ addStageTerms_(residual, stageParams_->size()-1,
+ makeLocalContext(fvGeometry_, elemVariables_));
+ }
+
+ return residual;
+ }
+ }
+
+ //! Return true if a stationary problem is assembled
+ bool isStationary() const
+ { return !stageParams_; }
+
+protected:
+
+ //! add the terms of a stage to the current element residual
+ void addStageTerms_(ElementResidualVector& r,
+ std::size_t stageIdx,
+ const LocalContext& context) const
+ {
+ LocalOperator localOperator(context);
+ if (!stageParams_->skipTemporal(stageIdx))
+ r.axpy(stageParams_->temporalWeight(stageIdx),
+ localOperator.evalStorage());
+ if (!stageParams_->skipSpatial(stageIdx))
+ r.axpy(stageParams_->spatialWeight(stageIdx),
+ localOperator.evalFluxesAndSources());
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the dofs of the given
+ * element and adds them to the global matrix.
+ * \return The element residual at the current solution.
+ */
+ template
+ ElementResidualVector assembleJacobianAndResidual_(JacobianMatrix& A,
+ const PartialReassembler* partialReassembler = nullptr)
+ {
+ // TODO: DO we need this to be constexpr?
+ if (diffMethod_ == DiffMethod::numeric)
+ return assembleJacobianAndResidualNumeric_(A, partialReassembler);
+ else
+ DUNE_THROW(Dune::NotImplemented, "Analytic assembler for box");
+ }
+
+ /*!
+ * \brief Computes the derivatives by means of numeric differentiation
+ * and adds them to the global matrix.
+ * \return The element residual at the current solution.
+ * \note This specialization is for the box scheme with numeric differentiation
+ */
+ template
+ ElementResidualVector assembleJacobianAndResidualNumeric_(JacobianMatrix& A,
+ const PartialReassembler* partialReassembler = nullptr)
+ {
+ // alias for the variables of the current stage
+ auto& curVariables = elemVariables_;
+ auto& curElemVolVars = curVariables.elemVolVars();
+ const auto& element = fvGeometry_.element();
+ const auto& problem = curElemVolVars.gridVolVars().problem();
+
+ const auto origResiduals = evalLocalResidual();
+ const Experimental::SolutionStateView solStateView(curVariables.gridVariables());
+ const auto origElemSol = elementSolution(element, solStateView, fvGeometry_.gridGeometry());
+ auto elemSol = origElemSol;
+
+ //////////////////////////////////////////////////////////////////////////////////////////////
+ // Calculate derivatives of the residual of all dofs in element with respect to themselves. //
+ //////////////////////////////////////////////////////////////////////////////////////////////
+
+ ElementResidualVector partialDerivs(fvGeometry_.numScv());
+ for (const auto& scvI : scvs(fvGeometry_))
+ {
+ // dof index and corresponding actual pri vars
+ const auto globalI = scvI.dofIndex();
+ const auto localI = scvI.localDofIndex();
+
+ const auto origCurVolVars = curElemVolVars[scvI];
+ auto& curVolVars = getVolVarAcess_(curElemVolVars, gridVariables_.gridVolVars(), scvI);
+
+ // 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[scvI.localDofIndex()][pvIdx] = priVar;
+ curVolVars.update(elemSol, problem, element, scvI);
+ return evalLocalResidual();
+ };
+
+ static const NumericEpsilon numEps{problem.paramGroup()};
+ static const int numDiffMethod = getParamFromGroup(problem.paramGroup(),
+ "Assembly.NumericDifferenceMethod");
+
+ // epsilon used for numeric differentiation
+ const auto eps = numEps(elemSol[localI][pvIdx], pvIdx);
+
+ // derive the residuals numerically
+ NumericDifferentiation::partialDerivative(evalResiduals, elemSol[localI][pvIdx], partialDerivs,
+ origResiduals, eps, numDiffMethod);
+
+ // TODO: Distinguish between implicit/explicit here. For explicit schemes,
+ // no entries between different scvs of an element are reserved. Thus,
+ // we currently get an error when using explicit schemes.
+ // TODO: Doesn't this mean we only have to evaluate the residual for a single
+ // scv instead of calling evalLocalResidual()? That computes the residuals
+ // and derivs for all other scvs of the element, too, which are never used.
+ // Note: this is the same in the current implementation of master.
+ // Should we try to optimize this for explicit schemes? Or adjust the Jacobian pattern?
+ // update the global stiffness matrix with the current partial derivatives
+ for (const auto& scvJ : scvs(fvGeometry_))
+ {
+ const auto globalJ = scvJ.dofIndex();
+ const auto localJ = scvJ.localDofIndex();
+
+ // don't add derivatives for green entities
+ if (!partialReassembler || partialReassembler->dofColor(globalJ) != EntityColor::green)
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ {
+ // A[i][col][eqIdx][pvIdx] is the rate of change of the
+ // the residual of equation 'eqIdx' at dof 'i'
+ // depending on the primary variable 'pvIdx' at dof 'col'
+ A[globalJ][globalI][eqIdx][pvIdx] += partialDerivs[localJ][eqIdx];
+ }
+ }
+ }
+
+ // restore the original element solution & volume variables
+ elemSol[localI][pvIdx] = origElemSol[localI][pvIdx];
+ curVolVars = origCurVolVars;
+
+ // TODO additional dof dependencies
+ }
+ }
+
+ return origResiduals;
+ }
+
+ //! Returns the volume variables of the grid vol vars (case of global caching)
+ template
+ auto& getVolVarAcess_(ElemVolVars& elemVolVars,
+ GridVolVars& gridVolVars,
+ const SubControlVolume& scv)
+ {
+ if constexpr (GridVolVars::cachingEnabled)
+ return gridVolVars.volVars(scv);
+ else
+ return elemVolVars[scv];
+ }
+
+ //! get the user-defined Dirichlet boundary conditions
+ template
+ auto getDirichletValues_(const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv,
+ const TimeLevel& timeLevel)
+ {
+ if constexpr(ProblemTraits::hasTransientDirichletInterface)
+ return problem.dirichlet(element, scv, timeLevel);
+ else
+ return problem.dirichlet(element, scv);
+ }
+
+ DiffMethod diffMethod_; //!< the type of differentiation method
+ GridVariables& gridVariables_; //!< reference to the grid variables
+ FVElementGeometry fvGeometry_; //!< element-local finite volume geometry
+ ElementVariables elemVariables_; //!< element variables of the current stage
+ std::vector prevElemVariables_; //!< element variables of prior stages
+
+ bool elementIsGhost_;
+ std::shared_ptr stageParams_;
+};
+
+} // end namespace Dumux::Experimental
+
+#endif
diff --git a/dumux/assembly/fv/cclocalassembler.hh b/dumux/assembly/fv/cclocalassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..396b732f7cf1a8ce6eaaf72e4f5c810c282477af
--- /dev/null
+++ b/dumux/assembly/fv/cclocalassembler.hh
@@ -0,0 +1,483 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief An assembler for Jacobian and residual contribution
+ * per element for cell-centered finite volume schemes.
+ */
+#ifndef DUMUX_CC_LOCAL_ASSEMBLER_HH
+#define DUMUX_CC_LOCAL_ASSEMBLER_HH
+
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+
+
+namespace Dumux::Experimental {
+
+/*!
+ * \ingroup Assembly
+ * \brief An assembler for Jacobian and residual contribution
+ * per element for cell-centered finite volume schemes.
+ * \tparam LO The element-local operator
+ */
+template
+class CCLocalAssembler
+{
+ using Operators = typename LO::Operators;
+ using LocalContext = typename LO::LocalContext;
+ using ElementVariables = typename LocalContext::ElementVariables;
+ using ElementGridGeometry = typename LocalContext::ElementGridGeometry;
+
+ using GG = typename ElementGridGeometry::GridGeometry;
+ using GV = typename ElementVariables::GridVariables;
+
+ using FVElementGeometry = typename GG::LocalView;
+ using SubControlVolume = typename GG::SubControlVolume;
+ using Element = typename GG::GridView::template Codim<0>::Entity;
+
+ using PrimaryVariables = typename GV::PrimaryVariables;
+ using NumEqVector = Dumux::NumEqVector;
+ using Scalar = typename GV::Scalar;
+
+ static constexpr int numEq = NumEqVectorTraits::numEq;
+ static constexpr int maxElementStencilSize = GG::maxElementStencilSize;
+ static constexpr bool enableGridFluxVarsCache = GV::GridFluxVariablesCache::cachingEnabled;
+
+public:
+
+ //! the parameters of a stage in time integration
+ using StageParams = MultiStageParams;
+
+ //! export the grid variables type on which to operate
+ using GridVariables = GV;
+
+ //! export the underlying local operator
+ using LocalOperator = LO;
+
+ //! export the vector storing the residuals of all dofs of the element
+ using ElementResidualVector = typename LO::ElementResidualVector;
+
+ /*!
+ * \brief Constructor for stationary problems.
+ */
+ explicit CCLocalAssembler(const Element& element,
+ GridVariables& gridVariables,
+ DiffMethod dm = DiffMethod::numeric)
+ : diffMethod_(dm)
+ , gridVariables_(gridVariables)
+ , fvGeometry_(localView(gridVariables.gridGeometry()))
+ , elemVariables_(localView(gridVariables))
+ , prevElemVariables_()
+ , elementIsGhost_(element.partitionType() == Dune::GhostEntity)
+ , stageParams_(nullptr)
+ {
+ if (diffMethod_ != DiffMethod::numeric)
+ DUNE_THROW(Dune::NotImplemented, "Provided differentiation method");
+
+ fvGeometry_.bind(element);
+ elemVariables_.bind(element, fvGeometry_);
+ }
+
+ /*!
+ * \brief Constructor for instationary problems.
+ * \note Using this constructor, we assemble one stage within
+ * a time integration step using multi-stage methods.
+ */
+ explicit CCLocalAssembler(const Element& element,
+ const std::vector& prevGridVars,
+ GridVariables& gridVariables,
+ std::shared_ptr stageParams,
+ DiffMethod dm = DiffMethod::numeric)
+ : diffMethod_(dm)
+ , gridVariables_(gridVariables)
+ , fvGeometry_(localView(gridVariables.gridGeometry()))
+ , elemVariables_(localView(gridVariables))
+ , prevElemVariables_()
+ , elementIsGhost_(element.partitionType() == Dune::GhostEntity)
+ , stageParams_(stageParams)
+ {
+ if (diffMethod_ != DiffMethod::numeric)
+ DUNE_THROW(Dune::NotImplemented, "Provided differentiation method");
+
+ if (prevGridVars.size() != stageParams->size() - 1)
+ DUNE_THROW(Dune::InvalidStateException, "Grid Variables for all stages needed");
+
+ fvGeometry_.bind(element);
+ for (const auto& gridVars : prevGridVars)
+ {
+ prevElemVariables_.emplace_back(localView(gridVars));
+ prevElemVariables_.back().bind(element, fvGeometry_);
+ }
+
+ elemVariables_.bind(element, 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,
+ const PartialReassembler* partialReassembler = nullptr)
+ {
+ const auto& element = fvGeometry_.element();
+ const auto globalI = fvGeometry_.gridGeometry().elementMapper().index(element);
+ if (partialReassembler
+ && partialReassembler->elementColor(globalI) == EntityColor::green)
+ {
+ res[globalI] = evalLocalResidual()[0]; // forward to the internal implementation
+ }
+ else
+ {
+ res[globalI] = assembleJacobianAndResidual_(jac); // forward to the internal implementation
+ }
+ }
+
+
+ /*!
+ * \brief Computes the derivatives with respect to the dofs of the given
+ * element and adds them to the given jacobian matrix.
+ */
+ template
+ void assembleJacobian(JacobianMatrix& jac)
+ {
+ assembleJacobianAndResidual_(jac);
+ }
+
+ //! Assemble the residual only
+ template
+ void assembleResidual(ResidualVector& res)
+ {
+ const auto residual = evalLocalResidual();
+ for (const auto& scv : scvs(fvGeometry_))
+ res[scv.dofIndex()] += residual[scv.localDofIndex()];
+ }
+
+ //! Evaluate the complete local residual for the current element.
+ ElementResidualVector evalLocalResidual() const
+ {
+ if (isStationary())
+ {
+ const auto context = makeLocalContext(fvGeometry_, elemVariables_);
+ LocalOperator localOperator(context);
+ return elementIsGhost_ ? localOperator.getEmptyResidual()
+ : localOperator.evalFluxesAndSources();
+ }
+ else
+ {
+ ElementResidualVector residual(fvGeometry_.numScv());
+ residual = 0.0;
+
+ if (!elementIsGhost_)
+ {
+ // add the terms associated with previous stages
+ for (std::size_t k = 0; k < stageParams_->size()-1; ++k)
+ addStageTerms_(residual, k,
+ makeLocalContext(fvGeometry_, prevElemVariables_[k]));
+
+ // add the terms of the current stage
+ addStageTerms_(residual, stageParams_->size()-1,
+ makeLocalContext(fvGeometry_, elemVariables_));
+ }
+
+ return residual;
+ }
+ }
+
+ //! Return true if a stationary problem is assembled
+ bool isStationary() const
+ { return !stageParams_; }
+
+protected:
+
+ //! add the terms of a stage to the current element residual
+ void addStageTerms_(ElementResidualVector& r,
+ std::size_t stageIdx,
+ const LocalContext& context) const
+ {
+ LocalOperator localOperator(context);
+ if (!stageParams_->skipTemporal(stageIdx))
+ r.axpy(stageParams_->temporalWeight(stageIdx),
+ localOperator.evalStorage());
+ if (!stageParams_->skipSpatial(stageIdx))
+ r.axpy(stageParams_->spatialWeight(stageIdx),
+ localOperator.evalFluxesAndSources());
+ }
+
+ /*!
+ * \brief Computes the derivatives with respect to the dofs of the given
+ * element and adds them to the global matrix.
+ * \return The element residual at the current solution.
+ */
+ template
+ NumEqVector assembleJacobianAndResidual_(JacobianMatrix& A,
+ const PartialReassembler* partialReassembler = nullptr)
+ {
+ // TODO: DO we need this to be constexpr?
+ if (diffMethod_ == DiffMethod::numeric)
+ return assembleJacobianAndResidualNumeric_(A, partialReassembler);
+ else
+ DUNE_THROW(Dune::NotImplemented, "Analytic assembler for box");
+ }
+
+ /*!
+ * \brief Computes the derivatives by means of numeric differentiation
+ * and adds them to the global matrix.
+ * \return The element residual at the current solution.
+ * \note This specialization is for the box scheme with numeric differentiation
+ */
+ template
+ NumEqVector assembleJacobianAndResidualNumeric_(JacobianMatrix& A,
+ const PartialReassembler* partialReassembler = nullptr)
+ {
+ // alias for the variables of the current stage
+ auto& curVariables = elemVariables_;
+ auto& curElemVolVars = curVariables.elemVolVars();
+ auto& curElemFluxVarsCache = curVariables.elemFluxVarsCache();
+ const auto& x = curVariables.gridVariables().dofs();
+
+ // get stencil informations
+ const auto& element = fvGeometry_.element();
+ const auto& gridGeometry = fvGeometry_.gridGeometry();
+ const auto& connectivityMap = gridGeometry.connectivityMap();
+ const auto globalI = gridGeometry.elementMapper().index(element);
+ const auto numNeighbors = connectivityMap[globalI].size();
+
+ // deduce the problem type
+ const auto& problem = curElemVolVars.gridVolVars().problem();
+ using Problem = std::decay_t;
+
+ // assemble the undeflected residual
+ using Residuals = ReservedBlockVector;
+ Residuals origResiduals(numNeighbors + 1); origResiduals = 0.0;
+ origResiduals[0] = evalLocalResidual()[0];
+
+ // lambda for convenient evaluation of the fluxes across scvfs in the neighbors
+ auto evalNeighborFlux = [&] (const auto& neighbor, const auto& scvf)
+ {
+ const auto context = makeLocalContext(fvGeometry_, elemVariables_);
+ return LocalOperator(context).computeFlux(neighbor, scvf);
+ };
+
+ // get the elements in which we need to evaluate the fluxes
+ // and calculate these in the undeflected state
+ Dune::ReservedVector neighborElements;
+ neighborElements.resize(numNeighbors);
+
+ unsigned int j = 1;
+ for (const auto& dataJ : connectivityMap[globalI])
+ {
+ neighborElements[j-1] = gridGeometry.element(dataJ.globalJ);
+
+ if (neighborElements[j-1].partitionType() != Dune::GhostEntity)
+ for (const auto scvfIdx : dataJ.scvfsJ)
+ origResiduals[j] += evalNeighborFlux(neighborElements[j-1],
+ fvGeometry_.scvf(scvfIdx));
+
+ ++j;
+ }
+
+ // reference to the element's scv (needed later) and corresponding vol vars
+ const auto& scv = fvGeometry_.scv(globalI);
+ auto& curVolVars = getVolVarAcess_(curElemVolVars, gridVariables_.gridVolVars(), scv);
+
+ // save a copy of the original privars and vol vars in order
+ // to restore the original solution after deflection
+ const auto origPriVars = x[globalI];
+ const auto origVolVars = curVolVars;
+
+ // element solution to be deflected
+ const Experimental::SolutionStateView solStateView(curVariables.gridVariables());
+ const auto origElemSol = elementSolution(element, solStateView, fvGeometry_.gridGeometry());
+ auto elemSol = origElemSol;
+
+ // derivatives in the neighbors with repect to the current elements
+ // in index 0 we save the derivative of the element residual with respect to it's own dofs
+ Residuals partialDerivs(numNeighbors + 1);
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+ {
+ partialDerivs = 0.0;
+
+ auto evalResiduals = [&](Scalar priVar)
+ {
+ Residuals partialDerivsTmp(numNeighbors + 1);
+ partialDerivsTmp = 0.0;
+
+ // update the volume variables and the flux var cache
+ elemSol[0][pvIdx] = priVar;
+ curVolVars.update(elemSol, problem, element, scv);
+ curElemFluxVarsCache.update(element, fvGeometry_, curElemVolVars);
+ if (enableGridFluxVarsCache)
+ gridVariables_.gridFluxVarsCache().updateElement(element, fvGeometry_, curElemVolVars);
+
+ // calculate the residual with the deflected primary variables
+ partialDerivsTmp[0] = evalLocalResidual()[0];
+
+ // calculate the fluxes in the neighbors with the deflected primary variables
+ for (std::size_t k = 0; k < numNeighbors; ++k)
+ if (neighborElements[k].partitionType() != Dune::GhostEntity)
+ for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
+ partialDerivsTmp[k+1] += evalNeighborFlux(neighborElements[k],
+ fvGeometry_.scvf(scvfIdx));
+
+ return partialDerivsTmp;
+ };
+
+ // derive the residuals numerically
+ static const NumericEpsilon numEps{problem.paramGroup()};
+ static const int numDiffMethod = getParamFromGroup(problem.paramGroup(), "Assembly.NumericDifferenceMethod");
+
+ const auto eps = numEps(elemSol[0][pvIdx], pvIdx);
+ NumericDifferentiation::partialDerivative(evalResiduals, elemSol[0][pvIdx],
+ partialDerivs, origResiduals,
+ eps, numDiffMethod);
+
+ // Correct derivative for ghost elements, i.e. set a 1 for the derivative w.r.t. the
+ // current primary variable and a 0 elsewhere. As we always solve for a delta of the
+ // solution with repect to the initial one, this results in a delta of 0 for ghosts.
+ if (elementIsGhost_)
+ {
+ partialDerivs[0] = 0.0;
+ partialDerivs[0][pvIdx] = 1.0;
+ }
+
+ // For instationary simulations, scale the coupling
+ // fluxes of the current stage correctly
+ if (stageParams_)
+ {
+ for (std::size_t k = 0; k < numNeighbors; ++k)
+ partialDerivs[k+1] *= stageParams_->spatialWeight(stageParams_->size()-1);
+ }
+
+ // add the current partial derivatives to the global jacobian matrix
+ // no special treatment is needed if globalJ is a ghost because then derivatives have been assembled to 0 above
+ if constexpr (Problem::enableInternalDirichletConstraints())
+ {
+ // check if own element has internal Dirichlet constraint
+ const auto internalDirichletConstraintsOwnElement = problem.hasInternalDirichletConstraint(element, scv);
+ const auto dirichletValues = problem.internalDirichlet(element, scv);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraintsOwnElement[eqIdx])
+ {
+ origResiduals[0][eqIdx] = origVolVars.priVars()[eqIdx] - dirichletValues[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] = (eqIdx == pvIdx) ? 1.0 : 0.0;
+ }
+ else
+ A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
+ }
+
+ // off-diagonal entries
+ j = 1;
+ for (const auto& dataJ : connectivityMap[globalI])
+ {
+ const auto& neighborElement = gridGeometry.element(dataJ.globalJ);
+ const auto& neighborScv = fvGeometry_.scv(dataJ.globalJ);
+ const auto internalDirichletConstraintsNeighbor = problem().hasInternalDirichletConstraint(neighborElement, neighborScv);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (internalDirichletConstraintsNeighbor[eqIdx])
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] = 0.0;
+ else
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j][eqIdx];
+ }
+
+ ++j;
+ }
+ }
+ else // no internal Dirichlet constraints specified
+ {
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ {
+ // the diagonal entries
+ A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
+
+ // off-diagonal entries
+ j = 1;
+ for (const auto& dataJ : connectivityMap[globalI])
+ A[dataJ.globalJ][globalI][eqIdx][pvIdx] += partialDerivs[j++][eqIdx];
+ }
+ }
+
+ // restore the original state of the scv's volume variables
+ curVolVars = origVolVars;
+
+ // restore the current element solution
+ elemSol[0][pvIdx] = origPriVars[pvIdx];
+ }
+
+ // restore original state of the flux vars cache in case of global caching.
+ // This has to be done in order to guarantee that everything is in an undeflected
+ // state before the assembly of another element is called. In the case of local caching
+ // this is obsolete because the elemFluxVarsCache used here goes out of scope after this.
+ // We only have to do this for the last primary variable, for all others the flux var cache
+ // is updated with the correct element volume variables before residual evaluations
+ curElemFluxVarsCache.update(element, fvGeometry_, curElemVolVars);
+ if (enableGridFluxVarsCache)
+ gridVariables_.gridFluxVarsCache().updateElement(element, fvGeometry_, curElemVolVars);
+
+ // return the original residual
+ return origResiduals[0];
+ }
+
+ //! Returns the volume variables of the local view (case of no caching)
+ template = 0>
+ auto& getVolVarAcess_(ElemVolVars& elemVolVars, GridVolVars& gridVolVars, const SubControlVolume& scv)
+ { return elemVolVars[scv]; }
+
+ //! Returns the volume variables of the grid vol vars (case of global caching)
+ template = 0>
+ auto& getVolVarAcess_(ElemVolVars& elemVolVars, GridVolVars& gridVolVars, const SubControlVolume& scv)
+ { return gridVolVars.volVars(scv); }
+
+ DiffMethod diffMethod_; //!< the type of differentiation method
+ GridVariables& gridVariables_; //!< reference to the grid variables
+ FVElementGeometry fvGeometry_; //!< element-local finite volume geometry
+ ElementVariables elemVariables_; //!< element variables of the current stage
+ std::vector prevElemVariables_; //!< element variables of prior stages
+
+ bool elementIsGhost_;
+ std::shared_ptr stageParams_;
+};
+
+} // end namespace Dumux::Experimental
+
+#endif
diff --git a/dumux/assembly/fv/localoperator.hh b/dumux/assembly/fv/localoperator.hh
new file mode 100644
index 0000000000000000000000000000000000000000..55cf9e9ad375339a384bb24e5e5e88f64a5c9e25
--- /dev/null
+++ b/dumux/assembly/fv/localoperator.hh
@@ -0,0 +1,271 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief An element-wise local operator for finite-volume schemes.
+ */
+#ifndef DUMUX_FV_LOCAL_OPERATOR_HH
+#define DUMUX_FV_LOCAL_OPERATOR_HH
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+
+namespace Dumux::Experimental {
+
+/*!
+ * \ingroup Assembly
+ * \brief The element-wise local operator for finite volume schemes.
+ * This allows for element-wise evaluation of individual terms
+ * of the equations to be solved.
+ * \tparam OP The model-specific operators
+ */
+template
+class FVLocalOperator
+{
+ using LC = typename OP::LocalContext;
+ using ElementVariables = typename LC::ElementVariables;
+ using GridVariables = typename ElementVariables::GridVariables;
+ using PrimaryVariables = typename GridVariables::PrimaryVariables;
+ using NumEqVector = Dumux::NumEqVector;
+
+ using FVElementGeometry = typename LC::ElementGridGeometry;
+ using GridGeometry = typename FVElementGeometry::GridGeometry;
+ using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using Extrusion = Extrusion_t;
+
+ static constexpr int numEq = NumEqVectorTraits::numEq();
+ static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box;
+
+public:
+ //! export the expected local context type
+ using LocalContext = LC;
+
+ //! export the underlying operators
+ using Operators = OP;
+
+ //! vector type storing the operator values on all dofs of an element
+ //! TODO: Use ReservedBlockVector
+ using ElementResidualVector = Dune::BlockVector;
+
+ //! Constructor from a local context
+ explicit FVLocalOperator(const LocalContext& context)
+ : context_(context)
+ {}
+
+ /*!
+ * \name Main interface
+ * \note Methods used by the assembler to compute derivatives and residual
+ */
+ // \{
+
+ /*!
+ * \brief Compute the terms of the local residual that do not appear in
+ * time derivatives. These are the sources and the fluxes.
+ */
+ ElementResidualVector evalFluxesAndSources() const
+ {
+ auto result = getEmptyResidual();
+ const auto& problem = context_.elementVariables().gridVariables().gridVolVars().problem();
+ const auto& fvGeometry = context_.elementGridGeometry();
+
+ // source term
+ for (const auto& scv : scvs(fvGeometry))
+ result[scv.localDofIndex()] -= Operators::source(problem, context_, scv);
+
+ // flux term
+ for (const auto& scvf : scvfs(fvGeometry))
+ addFlux_(result, scvf);
+
+ return result;
+ }
+
+ /*!
+ * \brief Compute the storage term, i.e. the term appearing in the time derivative.
+ */
+ ElementResidualVector evalStorage() const
+ {
+ const auto& problem = context_.elementVariables().gridVariables().gridVolVars().problem();
+ const auto& fvGeometry = context_.elementGridGeometry();
+
+ // TODO: Until now, FVLocalResidual defined storage as the entire
+ // time derivative. Now it means the term above the time derivative.
+ // We should think about the correct naming here...
+ // TODO: Should storage() NOT multiply with volume?? That was different until
+ // now but flux() also returns the flux multiplied with area so this should
+ // be more consistent
+ auto result = getEmptyResidual();
+ for (const auto& scv : scvs(fvGeometry))
+ result[scv.localDofIndex()] += Operators::storage(problem, context_, scv);
+
+ return result;
+ }
+
+ ElementResidualVector getEmptyResidual() const
+ {
+ ElementResidualVector res(context_.elementGridGeometry().numScv());
+ res = 0.0;
+ return res;
+ }
+
+ // \}
+
+ /*!
+ * \name Interfaces for analytic Jacobian computation
+ */
+ // \{
+
+ //! \todo TODO: Add interfaces. Or, should this be here at all!?
+
+ //\}
+
+ // \}
+
+ /*!
+ * \brief Compute the flux across a single face embedded in the given element.
+ * \note This function behaves different than the flux function in the
+ * operators, as it checks if the face is on the boundary. If so,
+ * it returns the flux resulting from the user-specified conditions.
+ * \note This assumes that the given element and face are contained within
+ * the context that was used to instantiate this class.
+ */
+ NumEqVector computeFlux(const Element& element,
+ const SubControlVolumeFace& scvf) const
+ {
+ const auto& evv = context_.elementVariables().elemVolVars();
+ const auto& problem = evv.gridVolVars().problem();
+
+ // interior faces
+ if (!scvf.boundary())
+ return Operators::flux(problem, context_, scvf);
+
+ const auto& insideScv = context_.elementGridGeometry().scv(scvf.insideScvIdx());
+
+ // for cell-centered schemes, evaluate fluxes also across Dirichlet boundaries
+ if constexpr (!isBox)
+ {
+ const auto& bcTypes = problem.boundaryTypes(element, scvf);
+ if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
+ return Operators::flux(problem, context_, scvf);
+ else if (bcTypes.hasNeumann() && bcTypes.hasDirichlet())
+ DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions for cell-centered schemes. " <<
+ "Use pure boundary conditions by converting Dirichlet BCs to Robin BCs.");
+ }
+ // for the box scheme, only pure Neumann boundaries are supported
+ else
+ {
+ if (!problem.boundaryTypes(element, insideScv).hasOnlyNeumann())
+ DUNE_THROW(Dune::NotImplemented, "For the box scheme only Neumann boundary fluxes can be computed.");
+ }
+
+ // compute and return the Neumann boundary fluxes
+ auto neumannFluxes = problem.neumann(context_, scvf);
+ neumannFluxes *= Extrusion::area(scvf)*evv[insideScv].extrusionFactor();
+ return neumannFluxes;
+ }
+
+protected:
+
+ //! compute and add the flux across the given face to the container (cc schemes)
+ template = 0>
+ void addFlux_(ElementResidualVector& r, const SubControlVolumeFace& scvf) const
+ {
+ const auto& fvGeometry = context_.elementGridGeometry();
+ const auto& evv = context_.elementVariables().elemVolVars();
+ const auto& problem = evv.gridVolVars().problem();
+
+ const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto localDofIdx = insideScv.localDofIndex();
+
+ if (!scvf.boundary())
+ r[localDofIdx] += Operators::flux(problem, context_, scvf);
+ else
+ {
+ const auto& bcTypes = problem.boundaryTypes(fvGeometry.element(), scvf);
+
+ // Dirichlet boundaries
+ if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
+ r[localDofIdx] += Operators::flux(problem, context_, scvf);
+
+ // Neumann and Robin ("solution dependent Neumann") boundary conditions
+ else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
+ {
+ auto neumannFluxes = problem.neumann(context_, scvf);
+ neumannFluxes *= Extrusion::area(scvf)*evv[insideScv].extrusionFactor();
+ r[localDofIdx] += neumannFluxes;
+ }
+
+ else
+ DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions for cell-centered schemes. " <<
+ "Use pure boundary conditions by converting Dirichlet BCs to Robin BCs");
+ }
+ }
+
+ //! compute and add the flux across the given face to the container (box scheme)
+ template = 0>
+ void addFlux_(ElementResidualVector& r, const SubControlVolumeFace& scvf) const
+ {
+ const auto& fvGeometry = context_.elementGridGeometry();
+ const auto& evv = context_.elementVariables().elemVolVars();
+ const auto& problem = evv.gridVolVars().problem();
+
+ // inner faces
+ if (!scvf.boundary())
+ {
+ const auto flux = Operators::flux(problem, context_, scvf);
+ r[fvGeometry.scv(scvf.insideScvIdx()).localDofIndex()] += flux;
+
+ if (scvf.numOutsideScvs() > 0)
+ r[fvGeometry.scv(scvf.outsideScvIdx()).localDofIndex()] -= flux;
+ }
+
+ // boundary faces
+ else
+ {
+ const auto& scv = fvGeometry.scv(scvf.insideScvIdx());
+ const auto& bcTypes = problem.boundaryTypes(fvGeometry.element(), scv);
+
+ // Treat Neumann and Robin ("solution dependent Neumann") boundary conditions.
+ if (bcTypes.hasNeumann())
+ {
+ const auto neumannFluxes = problem.neumann(context_, scvf);
+ const auto area = Extrusion::area(scvf)*evv[scv].extrusionFactor();
+
+ // only add fluxes to equations for which Neumann is set
+ for (int eqIdx = 0; eqIdx < NumEqVector::size(); ++eqIdx)
+ if (bcTypes.isNeumann(eqIdx))
+ r[scv.localDofIndex()][eqIdx] += neumannFluxes[eqIdx]*area;
+ }
+ }
+ }
+
+private:
+ const LocalContext& context_;
+};
+
+} // end namespace Dumux::Experimental
+
+#endif
diff --git a/dumux/assembly/fv/operators.hh b/dumux/assembly/fv/operators.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0fbc70812bed097609248d9878eeb510d2e2df8f
--- /dev/null
+++ b/dumux/assembly/fv/operators.hh
@@ -0,0 +1,155 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief The base class for the sub-control entity-local evaluation of
+ * the terms of equations in the context of finite-volume schemes
+ */
+#ifndef DUMUX_FV_OPERATORS_HH
+#define DUMUX_FV_OPERATORS_HH
+
+#include
+#include
+#include
+
+namespace Dumux::Experimental {
+
+/*!
+ * \ingroup Assembly
+ * \brief Convenience alias to define the context finite-volume operators work on.
+ * \tparam EV The element-(stencil)-local variables
+ */
+template
+using FVOperatorsContext = DefaultLocalContext;
+
+/*!
+ * \ingroup Assembly
+ * \brief The base class for the sub-control entity-local evaluation of
+ * the terms of equations in the context of finite-volume schemes
+ * \tparam EV The element-(stencil)-local variables
+ */
+template
+class FVOperators
+{
+ // context type on which to operate
+ using LC = FVOperatorsContext;
+
+ // The grid geometry on which the scheme operates
+ using FVElementGeometry = typename LC::ElementGridGeometry;
+ using GridGeometry = typename FVElementGeometry::GridGeometry;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+ using Extrusion = Extrusion_t;
+
+ // The variables (primary/secondary) on the grid
+ using GridVariables = typename LC::ElementVariables::GridVariables;
+ using PrimaryVariables = typename GridVariables::PrimaryVariables;
+
+public:
+ //! export the local context on which this operates
+ using LocalContext = LC;
+
+ //! export the type used to store scalar values for all equations
+ using NumEqVector = Dumux::NumEqVector;
+
+ // export the types of the flux/source/storage terms
+ using FluxTerm = NumEqVector;
+ using SourceTerm = NumEqVector;
+ using StorageTerm = NumEqVector;
+
+ /*!
+ * \name Model specific interfaces
+ * \note The following method are the model specific implementations of the
+ * operators appearing in the equation and should be overloaded by the
+ * implementations.
+ */
+ // \{
+
+ /*!
+ * \brief Compute the storage term of the equations for the given sub-control volume
+ * \param problem The problem to be solved (could store additionally required quantities)
+ * \param context The element-stencil-local required data
+ * \param scv The sub-control volume
+ * \note This must be overloaded by the implementation
+ * \note In cell-centered schemes, depending on the assembler or the policy used
+ * therein, the context may contain variables (e.g. from neighboring cells
+ * within the stencil) with respect to which the storage term is not
+ * differentiated. This means that if a Newton solver is used to solve
+ * the system of equations, one ends up with a quasi-Newton scheme. If
+ * this is undesireable, make sure to set up an assembler that takes into
+ * account the dependencies of the storage term w.r.t variables in the context.
+ */
+ template
+ static StorageTerm storage(const Problem& problem,
+ const LocalContext& context,
+ const SubControlVolume& scv)
+ { DUNE_THROW(Dune::NotImplemented, "Storage operator not implemented!"); }
+
+ /*!
+ * \brief Compute the flux term of the equations for the given sub-control volume face
+ * \param problem The problem to be solved (could store additionally required quantities)
+ * \param context The element-stencil-local required data
+ * \param scvf The sub-control volume face for which the flux term is to be computed
+ * \note This must be overloaded by the implementation
+ */
+ template
+ static FluxTerm flux(const Problem& problem,
+ const LocalContext& context,
+ const SubControlVolumeFace& scvf)
+ { DUNE_THROW(Dune::NotImplemented, "This model does not implement a flux method!"); }
+
+ /*!
+ * \brief Compute the source term of the equations for the given sub-control volume
+ * \param problem The problem to be solved (could store additionally required quantities)
+ * \param context The element-stencil-local required data
+ * \param scv The sub-control volume for which the source term is to be computed
+ * \note This is a default implementation forwarding to interfaces in the problem
+ * \note In cell-centered schemes, depending on the assembler or the policy used
+ * therein, the context may contain variables (e.g. from neighboring cells
+ * within the stencil) with respect to which the source term is not
+ * differentiated. This means that if a Newton solver is used to solve
+ * the system of equations, one ends up with a quasi-Newton scheme. If
+ * this is undesireable, make sure to set up an assembler that takes into
+ * account the dependencies of the storage term w.r.t variables in the context.
+ */
+ template
+ static SourceTerm source(const Problem& problem,
+ const LocalContext& context,
+ const SubControlVolume& scv)
+ {
+ SourceTerm source(0.0);
+
+ // add contributions from volume flux sources
+ source += problem.source(context, scv);
+
+ // add contribution from possible point sources
+ source += problem.scvPointSources(context, scv);
+
+ // multiply with scv volume
+ const auto& elemVolVars = context.elementVariables().elemVolVars();
+ source *= Extrusion::volume(scv)*elemVolVars[scv].extrusionFactor();
+
+ return source;
+ }
+};
+
+} // end namespace Dumux::Experimental
+
+#endif
diff --git a/dumux/assembly/localassembler.hh b/dumux/assembly/localassembler.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5e51dfd9f401d57920fe5c8d4679b0e34c1b6692
--- /dev/null
+++ b/dumux/assembly/localassembler.hh
@@ -0,0 +1,65 @@
+// -*- 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 . *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Assembly
+ * \brief Convenience alias to select a local assembler based on the discretization scheme.
+ */
+#ifndef DUMUX_LOCAL_ASSEMBLER_HH
+#define DUMUX_LOCAL_ASSEMBLER_HH
+
+#include
+#include "fv/boxlocalassembler.hh"
+#include "fv/cclocalassembler.hh"
+
+namespace Dumux::Experimental {
+namespace Detail {
+
+ template
+ struct LocalAssemblerChooser;
+
+ template
+ struct LocalAssemblerChooser
+ { using type = BoxLocalAssembler; };
+
+ template
+ struct LocalAssemblerChooser
+ { using type = CCLocalAssembler; };
+
+ template
+ struct LocalAssemblerChooser
+ { using type = CCLocalAssembler; };
+
+ template
+ using LocalAssemblerType
+ = typename LocalAssemblerChooser::type;
+
+} // end namespace Detail
+
+/*!
+ * \ingroup Assembly
+ * \brief Helper alias to select the local assembler type from a local operator.
+ */
+template
+using LocalAssembler = Detail::LocalAssemblerType;
+
+} // end namespace Dumux::Experimental
+
+#endif
diff --git a/dumux/common/fvproblem.hh b/dumux/common/fvproblem.hh
index f3063569475d9ea073ead553445ac21c90927da7..ed4d05e119f3ab4cece3b646103e5216abfa116a 100644
--- a/dumux/common/fvproblem.hh
+++ b/dumux/common/fvproblem.hh
@@ -601,6 +601,156 @@ private:
PointSourceMap pointSourceMap_;
};
+namespace Experimental {
+
+//! Experimental problem implementation compatible with new time
+//! integration schemes and corresponding assembly.
+template
+class FVProblem : public Dumux::FVProblem
+{
+ using ParentType = Dumux::FVProblem;
+
+ using GridGeometry = GetPropType;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace;
+
+ using GridView = typename GridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+ using PrimaryVariables = typename ParentType::Traits::PrimaryVariables;
+ using NumEqVector = typename ParentType::Traits::NumEqVector;
+ using Scalar = typename ParentType::Traits::Scalar;
+
+ static constexpr bool isBox = GridGeometry::discMethod == DiscretizationMethod::box;
+
+public:
+ // pull up constructor
+ using ParentType::ParentType;
+
+ /*!
+ * \brief Specify the kind of boundary conditions used on a discrete
+ * entity on the boundary.
+ *
+ * \param element The finite element
+ * \param boundaryEntity The boundary entity (scv/scvf)
+ * \note In cell-centered schemes, boundaryEntity is a sub-control
+ * volume face (scvf). In the box scheme, a sub-control volume (scv).
+ */
+ template
+ auto boundaryTypes(const Element& element,
+ const BoundaryEntity& boundaryEntity) const
+ {
+ if constexpr (isBox)
+ return this->asImp_().boundaryTypesAtPos(boundaryEntity.dofPosition());
+ else
+ return this->asImp_().boundaryTypesAtPos(boundaryEntity.ipGlobal());
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a discrete entity on the boundary.
+ *
+ * \param element The finite element
+ * \param boundaryEntity The boundary entity (scv/scvf)
+ * \param timeLevel The time level on which to evaluate the boundary conditions
+ * \note In cell-centered schemes, boundaryEntity is a sub-control
+ * volume face (scvf). In the box scheme, a sub-control volume (scv).
+ */
+ template
+ PrimaryVariables dirichlet(const Element& element,
+ const BoundaryEntity& boundaryEntity,
+ const TimeLevel& timeLevel) const
+ { return this->asImp_().dirichlet(element, boundaryEntity); }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a discrete entity on the boundary.
+ *
+ * \param element The finite element
+ * \param boundaryEntity The boundary entity (scv/scvf)
+ * \note In cell-centered schemes, boundaryEntity is a sub-control
+ * volume face (scvf). In the box scheme, a sub-control volume (scv).
+ */
+ template
+ PrimaryVariables dirichlet(const Element& element,
+ const BoundaryEntity& boundaryEntity) const
+ {
+ if constexpr (isBox)
+ return this->asImp_().dirichletAtPos(boundaryEntity.dofPosition());
+ else
+ return this->asImp_().dirichletAtPos(boundaryEntity.ipGlobal());
+ }
+
+ /*!
+ * \brief Evaluate the boundary conditions for a neumann
+ * boundary segment.
+ *
+ * This is the method for the case where the Neumann condition is
+ * potentially solution dependent
+ *
+ * \param context The element-local context
+ * \param scvf The sub control volume face
+ * \note The element-local context consists of an element and local
+ * geometric information, together with the primary/secondary
+ * variables in that local scope. Potentially, users can define
+ * the context to carry an additional object containing further
+ * locally required data.
+ */
+ template
+ NumEqVector neumann(const LocalContext& context,
+ const SubControlVolumeFace& scvf) const
+ { return this->asImp_().neumannAtPos(scvf.ipGlobal()); }
+
+ /*!
+ * \brief Evaluate the source term for all phases within a given
+ * sub-control-volume.
+ *
+ * This is the method for the case where the source term is
+ * potentially solution dependent and requires some quantities that
+ * are specific to the fully-implicit method.
+ *
+ * \param context The element-local context
+ * \param scv The sub control volume
+ * \note The element-local context consists of an element and local
+ * geometric information, together with the primary/secondary
+ * variables in that local scope. Potentially, users can define
+ * the context to carry an additional object containing further
+ * locally required data.
+ *
+ * For this method, the return parameter stores the conserved quantity rate
+ * generated or annihilate per volume unit. Positive values mean
+ * that the conserved quantity is created, negative ones mean that it vanishes.
+ * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
+ */
+ template
+ NumEqVector source(const LocalContext& context,
+ const SubControlVolume& scv) const
+ { return this->asImp_().sourceAtPos(scv.center()); }
+
+ /*!
+ * \brief Adds contribution of point sources for a specific sub control volume
+ * to the values.
+ *
+ * \param context The element-local context
+ * \param scv The sub control volume
+ * \note The element-local context consists of an element and local
+ * geometric information, together with the primary/secondary
+ * variables in that local scope. Potentially, users can define
+ * the context to carry an additional object containing further
+ * locally required data.
+ *
+ */
+ template
+ NumEqVector scvPointSources(const LocalContext& context,
+ const SubControlVolume& scv) const
+ {
+ const auto& elemVolVars = context.elementVariables().elemVolVars();
+ const auto& fvGeometry = context.elementGridGeometry();
+ const auto& element = fvGeometry.element();
+ return ParentType::scvPointSources(element, fvGeometry, elemVolVars, scv);
+ }
+};
+
+} // end namespace Experimental
} // end namespace Dumux
#endif
diff --git a/dumux/common/typetraits/problem.hh b/dumux/common/typetraits/problem.hh
index 0131254abd3831b190cb0ac61025ca26c5a5c4a3..4de22d482f792281589040b7d10bbefdd018bfb3 100644
--- a/dumux/common/typetraits/problem.hh
+++ b/dumux/common/typetraits/problem.hh
@@ -25,6 +25,9 @@
#define DUMUX_TYPETRAITS_PROBLEM_HH
#include
+
+#include
+#include
#include
namespace Dumux {
@@ -33,6 +36,29 @@ namespace Dumux {
namespace Detail {
template
struct ProblemTraits;
+
+template
+struct hasTransientDirichlet
+{
+private:
+ using TimeLevel = Dumux::Experimental::TimeLevel;
+ using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+ using BoundaryEntity = std::conditional_t;
+public:
+ template
+ auto operator()(const Problem& p)
+ -> decltype(p.dirichlet(std::declval(),
+ std::declval(),
+ std::declval()))
+ {}
+};
+
+template
+inline constexpr bool hasTransientDirichletInterface
+ = decltype(isValid(hasTransientDirichlet())(std::declval()))::value;
+
} // end namespace Detail
/*!
@@ -44,6 +70,9 @@ struct ProblemTraits
{
using GridGeometry = std::decay_t().gridGeometry())>;
using BoundaryTypes = typename Detail::template ProblemTraits::BoundaryTypes;
+
+ static constexpr bool hasTransientDirichletInterface
+ = Detail::hasTransientDirichletInterface;
};
} // end namespace Dumux
diff --git a/dumux/discretization/CMakeLists.txt b/dumux/discretization/CMakeLists.txt
index 3b145a92fb493e50c73f4afb2f7dfafc068b10b5..82f9d4dcf9ece340cf0b9bc966175ac67c7af98b 100644
--- a/dumux/discretization/CMakeLists.txt
+++ b/dumux/discretization/CMakeLists.txt
@@ -11,6 +11,7 @@ box.hh
ccmpfa.hh
cctpfa.hh
checkoverlapsize.hh
+concepts.hh
elementsolution.hh
evalgradients.hh
evalsolution.hh
@@ -20,6 +21,7 @@ functionspacebasis.hh
fvgridvariables.hh
fvproperties.hh
gridvariables.hh
+localcontext.hh
localview.hh
method.hh
rotationpolicy.hh
@@ -27,6 +29,7 @@ rotationsymmetricgridgeometrytraits.hh
rotationsymmetricscv.hh
rotationsymmetricscvf.hh
scvandscvfiterators.hh
+solutionstateview.hh
staggered.hh
subcontrolvolumebase.hh
subcontrolvolumefacebase.hh
diff --git a/dumux/discretization/box/elementvolumevariables.hh b/dumux/discretization/box/elementvolumevariables.hh
index 2ddb9591a67f91b90822bf71a67c3da68e8a3d4b..e072f7ff6629133da0a6f464e260f6a8b157bb7a 100644
--- a/dumux/discretization/box/elementvolumevariables.hh
+++ b/dumux/discretization/box/elementvolumevariables.hh
@@ -68,19 +68,19 @@ public:
// For compatibility reasons with the case of not storing the vol vars.
// function to be called before assembling an element, preparing the vol vars within the stencil
- template
+ template
void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
const FVElementGeometry& fvGeometry,
- const SolutionVector& sol)
+ const SolutionState& solState)
{
- bindElement(element, fvGeometry, sol);
+ bindElement(element, fvGeometry, solState);
}
// function to prepare the vol vars within the element
- template
+ template
void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
const FVElementGeometry& fvGeometry,
- const SolutionVector& sol)
+ const SolutionState& solState)
{
eIdx_ = fvGeometry.gridGeometry().elementMapper().index(element);
}
@@ -114,22 +114,21 @@ public:
: gridVolVarsPtr_(&gridVolVars) {}
// specialization for box models, simply forwards to the bindElement method
- template
+ template
void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
const FVElementGeometry& fvGeometry,
- const SolutionVector& sol)
+ const SolutionState& solState)
{
- bindElement(element, fvGeometry, sol);
+ bindElement(element, fvGeometry, solState);
}
// specialization for box models
- template
+ template
void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
const FVElementGeometry& fvGeometry,
- const SolutionVector& sol)
+ const SolutionState& solState)
{
- // get the solution at the dofs of the element
- auto elemSol = elementSolution(element, sol, fvGeometry.gridGeometry());
+ const auto elemSol = elementSolution(element, solState, fvGeometry.gridGeometry());
// resize volume variables to the required size
volumeVariables_.resize(fvGeometry.numScv());
diff --git a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
index 83b6b3a93acab074b3d263b98156c0c843a34969..1666053e045b026c6c89f9134b1ad0d47a439bcd 100644
--- a/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
+++ b/dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh
@@ -29,6 +29,7 @@
#include
#include
+#include
#include
namespace Dumux {
@@ -76,20 +77,23 @@ namespace CCMpfa {
* \param element The element to which the finite volume geometry is bound
* \param fvGeometry The element finite volume geometry
* \param nodalIndexSet The dual grid index set around a node
+ * \param solState The state of the discrete solution
*/
- template
+ template
void addBoundaryVolVarsAtNode(std::vector& volVars,
std::vector& volVarIndices,
const Problem& problem,
const typename FVElemGeom::GridGeometry::GridView::template Codim<0>::Entity& element,
const FVElemGeom& fvGeometry,
- const NodalIndexSet& nodalIndexSet)
+ const NodalIndexSet& nodalIndexSet,
+ const SolState& solState)
{
if (nodalIndexSet.numBoundaryScvfs() == 0)
return;
// index of the element the fvGeometry was bound to
const auto boundElemIdx = fvGeometry.gridGeometry().elementMapper().index(element);
+ auto elemSol = elementSolution(element, solState, fvGeometry.gridGeometry());
// check each scvf in the index set for boundary presence
for (auto scvfIdx : nodalIndexSet.gridScvfIndices())
@@ -108,11 +112,11 @@ namespace CCMpfa {
// boundaries the "outside" vol vars cannot be properly defined.
if (bcTypes.hasOnlyDirichlet())
{
+ using namespace Experimental::CompatibilityHelpers;
+ elemSol[0] = getDirichletValues(problem, insideElement, ivScvf, elemSol);
+
VolumeVariables dirichletVolVars;
- dirichletVolVars.update(elementSolution(problem.dirichlet(insideElement, ivScvf)),
- problem,
- insideElement,
- fvGeometry.scv(insideScvIdx));
+ dirichletVolVars.update(elemSol, problem, insideElement, fvGeometry.scv(insideScvIdx));
volVars.emplace_back(std::move(dirichletVolVars));
volVarIndices.push_back(ivScvf.outsideScvIdx());
@@ -130,15 +134,18 @@ namespace CCMpfa {
* \param problem The problem containing the Dirichlet boundary conditions
* \param element The element to which the finite volume geometry was bound
* \param fvGeometry The element finite volume geometry
+ * \param solState The state of the discrete solution
*/
- template
+ template
void addBoundaryVolVars(std::vector& volVars,
std::vector& volVarIndices,
const Problem& problem,
const typename FVElemGeom::GridGeometry::GridView::template Codim<0>::Entity& element,
- const FVElemGeom& fvGeometry)
+ const FVElemGeom& fvGeometry,
+ const SolState& solState)
{
const auto& gridGeometry = fvGeometry.gridGeometry();
+ auto elemSol = elementSolution(element, solState, fvGeometry.gridGeometry());
// treat the BCs inside the element
if (fvGeometry.hasBoundaryScvf())
@@ -155,11 +162,11 @@ namespace CCMpfa {
// boundaries the "outside" vol vars cannot be properly defined.
if (problem.boundaryTypes(element, scvf).hasOnlyDirichlet())
{
+ using namespace Experimental::CompatibilityHelpers;
+ elemSol[0] = getDirichletValues(problem, element, scvf, elemSol);
+
VolumeVariables dirichletVolVars;
- dirichletVolVars.update(elementSolution(problem.dirichlet(element, scvf)),
- problem,
- element,
- scvI);
+ dirichletVolVars.update(elemSol, problem, element, scvI);
volVars.emplace_back(std::move(dirichletVolVars));
volVarIndices.push_back(scvf.outsideScvIdx());
@@ -173,10 +180,10 @@ namespace CCMpfa {
{
if (!gridGeometry.vertexUsesSecondaryInteractionVolume(scvf.vertexIndex()))
addBoundaryVolVarsAtNode( volVars, volVarIndices, problem, element, fvGeometry,
- gridIvIndexSets.primaryIndexSet(scvf).nodalIndexSet() );
+ gridIvIndexSets.primaryIndexSet(scvf).nodalIndexSet(), solState );
else
addBoundaryVolVarsAtNode( volVars, volVarIndices, problem, element, fvGeometry,
- gridIvIndexSets.secondaryIndexSet(scvf).nodalIndexSet() );
+ gridIvIndexSets.secondaryIndexSet(scvf).nodalIndexSet(), solState );
}
}
} // end namespace CCMpfa
@@ -228,10 +235,10 @@ public:
}
//! precompute all volume variables in a stencil of an element - bind Dirichlet vol vars in the stencil
- template
+ template
void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
const FVElementGeometry& fvGeometry,
- const SolutionVector& sol)
+ const SolutionState& solState)
{
clear();
@@ -241,7 +248,8 @@ public:
{
boundaryVolVars_.reserve(maxNumBoundaryVolVars);
boundaryVolVarIndices_.reserve(maxNumBoundaryVolVars);
- CCMpfa::addBoundaryVolVars(boundaryVolVars_, boundaryVolVarIndices_, gridVolVars().problem(), element, fvGeometry);
+ CCMpfa::addBoundaryVolVars(boundaryVolVars_, boundaryVolVarIndices_,
+ gridVolVars().problem(), element, fvGeometry, solState);
}
}
@@ -299,10 +307,10 @@ public:
: gridVolVarsPtr_(&gridVolVars) {}
//! Prepares the volume variables within the element stencil
- template
+ template
void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
const FVElementGeometry& fvGeometry,
- const SolutionVector& sol)
+ const SolutionState& solState)
{
clear();
@@ -321,10 +329,8 @@ public:
VolumeVariables volVars;
const auto& scvI = fvGeometry.scv(globalI);
- volVars.update(elementSolution(element, sol, gridGeometry),
- problem,
- element,
- scvI);
+ const auto elemSol = elementSolution(element, solState, gridGeometry);
+ volVars.update(elemSol, problem, element, scvI);
volVarIndices_.push_back(scvI.dofIndex());
volumeVariables_.emplace_back(std::move(volVars));
@@ -333,12 +339,10 @@ public:
for (auto&& dataJ : assemblyMapI)
{
const auto& elementJ = gridGeometry.element(dataJ.globalJ);
+ const auto elemSolJ = elementSolution(elementJ, solState, gridGeometry);
const auto& scvJ = fvGeometry.scv(dataJ.globalJ);
VolumeVariables volVarsJ;
- volVarsJ.update(elementSolution(elementJ, sol, gridGeometry),
- problem,
- elementJ,
- scvJ);
+ volVarsJ.update(elemSolJ, problem, elementJ, scvJ);
volVarIndices_.push_back(scvJ.dofIndex());
volumeVariables_.emplace_back(std::move(volVarsJ));
@@ -346,7 +350,8 @@ public:
// maybe prepare boundary volume variables
if (maxNumBoundaryVolVars > 0)
- CCMpfa::addBoundaryVolVars(volumeVariables_, volVarIndices_, problem, element, fvGeometry);
+ CCMpfa::addBoundaryVolVars(volumeVariables_, volVarIndices_,
+ problem, element, fvGeometry, solState);
// //! TODO Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
// //! on additional DOFs not included in the discretization schemes' occupation pattern
@@ -373,10 +378,10 @@ public:
}
//! Prepares the volume variables of an element
- template
+ template
void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
const FVElementGeometry& fvGeometry,
- const SolutionVector& sol)
+ const SolutionState& solState)
{
clear();
@@ -387,7 +392,7 @@ public:
// update the volume variables of the element
const auto& scv = fvGeometry.scv(eIdx);
- volumeVariables_[0].update(elementSolution(element, sol, gridGeometry),
+ volumeVariables_[0].update(elementSolution(element, solState, gridGeometry),
gridVolVars().problem(),
element,
scv);
diff --git a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh
index e10d37ac5cbae1cc063cf1378144a11a855f0be6..19da0f37744dfd9b3ba07b1bb81dbe7ae7adc3cf 100644
--- a/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh
+++ b/dumux/discretization/cellcentered/tpfa/elementvolumevariables.hh
@@ -28,6 +28,7 @@
#include
#include
+#include
#include
namespace Dumux {
@@ -84,10 +85,10 @@ public:
}
//! precompute all boundary volume variables in a stencil of an element, the remaining ones are cached
- template
+ template
void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
const FVElementGeometry& fvGeometry,
- const SolutionVector& sol)
+ const SolutionState& solutionState)
{
if (!fvGeometry.hasBoundaryScvf())
return;
@@ -96,6 +97,7 @@ public:
boundaryVolVarIndices_.reserve(fvGeometry.numScvf());
boundaryVolumeVariables_.reserve(fvGeometry.numScvf());
+ auto elemSol = elementSolution(element, solutionState, fvGeometry.gridGeometry());
for (const auto& scvf : scvfs(fvGeometry))
{
if (!scvf.boundary())
@@ -106,14 +108,11 @@ public:
const auto bcTypes = problem.boundaryTypes(element, scvf);
if (bcTypes.hasOnlyDirichlet())
{
- const auto dirichletPriVars = elementSolution(problem.dirichlet(element, scvf));
- auto&& scvI = fvGeometry.scv(scvf.insideScvIdx());
-
+ using namespace Experimental::CompatibilityHelpers;
+ elemSol[0] = getDirichletValues(problem, element, scvf, elemSol);
+ const auto& scvI = fvGeometry.scv(scvf.insideScvIdx());
VolumeVariables volVars;
- volVars.update(dirichletPriVars,
- problem,
- element,
- scvI);
+ volVars.update(elemSol, problem, element, scvI);
boundaryVolumeVariables_.emplace_back(std::move(volVars));
boundaryVolVarIndices_.push_back(scvf.outsideScvIdx());
@@ -174,10 +173,10 @@ public:
: gridVolVarsPtr_(&gridVolVars) {}
//! Prepares the volume variables within the element stencil
- template
+ template
void bind(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
const FVElementGeometry& fvGeometry,
- const SolutionVector& sol)
+ const SolutionState& solState)
{
clear_();
@@ -193,8 +192,9 @@ public:
int localIdx = 0;
// update the volume variables of the element at hand
- auto&& scvI = fvGeometry.scv(globalI);
- volumeVariables_[localIdx].update(elementSolution(element, sol, gridGeometry),
+ auto elemSol = elementSolution(element, solState, gridGeometry);
+ const auto& scvI = fvGeometry.scv(globalI);
+ volumeVariables_[localIdx].update(elemSol,
problem,
element,
scvI);
@@ -206,7 +206,7 @@ public:
{
const auto& elementJ = gridGeometry.element(dataJ.globalJ);
auto&& scvJ = fvGeometry.scv(dataJ.globalJ);
- volumeVariables_[localIdx].update(elementSolution(elementJ, sol, gridGeometry),
+ volumeVariables_[localIdx].update(elementSolution(elementJ, solState, gridGeometry),
problem,
elementJ,
scvJ);
@@ -227,18 +227,16 @@ public:
const auto bcTypes = problem.boundaryTypes(element, scvf);
if (bcTypes.hasOnlyDirichlet())
{
- const auto dirichletPriVars = elementSolution(problem.dirichlet(element, scvf));
+ using namespace Experimental::CompatibilityHelpers;
+ elemSol[0] = getDirichletValues(problem, element, scvf, elemSol);
volumeVariables_.resize(localIdx+1);
volVarIndices_.resize(localIdx+1);
- volumeVariables_[localIdx].update(dirichletPriVars,
- problem,
- element,
- scvI);
- volVarIndices_[localIdx] = scvf.outsideScvIdx();
- ++localIdx;
- }
+ volumeVariables_[localIdx].update(elemSol, problem, element, scvI);
+ volVarIndices_[localIdx] = scvf.outsideScvIdx();
+ ++localIdx;
}
+ }
}
//! Check if user added additional DOF dependencies, i.e. the residual of DOF globalI depends
@@ -263,10 +261,10 @@ public:
// }
}
- template
+ template
void bindElement(const typename FVElementGeometry::GridGeometry::GridView::template Codim<0>::Entity& element,
const FVElementGeometry& fvGeometry,
- const SolutionVector& sol)
+ const SolutionState& solState)
{
clear_();
@@ -275,8 +273,9 @@ public:
volVarIndices_.resize(1);
// update the volume variables of the element
- auto&& scv = fvGeometry.scv(eIdx);
- volumeVariables_[0].update(elementSolution(element, sol, fvGeometry.gridGeometry()),
+ const auto elemSol = elementSolution(element, solState, fvGeometry.gridGeometry());
+ const auto& scv = fvGeometry.scv(eIdx);
+ volumeVariables_[0].update(elemSol,
gridVolVars().problem(),
element,
scv);
diff --git a/test/porousmediumflow/1p/compressible/instationary/assembler.hh b/dumux/discretization/concepts.hh
similarity index 52%
rename from test/porousmediumflow/1p/compressible/instationary/assembler.hh
rename to dumux/discretization/concepts.hh
index 115b25570c655bbae65f89bde25120ca7db1dc7e..65bdff436e218125ef874234c2bbc9da758ff396 100644
--- a/test/porousmediumflow/1p/compressible/instationary/assembler.hh
+++ b/dumux/discretization/concepts.hh
@@ -18,37 +18,45 @@
*****************************************************************************/
/*!
* \file
- * \ingroup OnePTests
- * \brief Dummy assembler that fulfills the new experimental assembly interface.
+ * \ingroup Discretization
+ * \brief Concepts related to discretization methods.
*/
-#ifndef DUMUX_COMPRESSIBLE_ONEP_TEST_ASSEMBLER_HH
-#define DUMUX_COMPRESSIBLE_ONEP_TEST_ASSEMBLER_HH
+#ifndef DUMUX_DISCRETIZATION_CONCEPT_HH
+#define DUMUX_DISCRETIZATION_CONCEPT_HH
-namespace Dumux::OnePCompressibleTest {
+#include
-// Custom assembler to test assembly with grid variables,
-// fulfilling the foreseen required interface
-template
-class GridVarsAssembler : public Assembler
-{
-public:
- using Assembler::Assembler;
- using typename Assembler::GridVariables;
- using typename Assembler::ResidualType;
-
- using Variables = GridVariables;
+namespace Dumux{
- void assembleJacobianAndResidual(const GridVariables& gridVars)
- { Assembler::assembleJacobianAndResidual(gridVars.dofs()); }
+// experimental concepts
+namespace Experimental::Concept {
- void assembleResidual(const GridVariables& gridVars)
- { Assembler::assembleResidual(gridVars.dofs()); }
+//! Concept for states of discrete numerical solutions
+struct SolutionState
+{
+ template
+ auto require(const T& t) -> decltype(
+ Dune::Concept::requireType(),
+ Dune::Concept::requireType(),
+ Dune::Concept::requireConvertible(t.dofs()),
+ Dune::Concept::requireConvertible(t.timeLevel())
+ );
+};
- //! Remove residualNorm once this is fixed in the solvers !2113
- auto residualNorm(const GridVariables& gridVars)
- { return Assembler::residualNorm(gridVars.dofs()); }
+//! Concept for element-local discrete numerical solutions
+struct ElementSolution
+{
+ template
+ auto require(const T& t) -> decltype(
+ Dune::Concept::requireType(),
+ Dune::Concept::requireType(),
+ Dune::Concept::requireConvertible(t.size()),
+ Dune::Concept::requireConvertible(t[0]),
+ Dune::Concept::requireConvertible(t.timeLevel())
+ );
};
-} // end namespace Dumux::OnePCompressibleTest
+} // end namespace Experimental::Concept
+} // end namespace Dumux
#endif
diff --git a/dumux/discretization/elementsolution.hh b/dumux/discretization/elementsolution.hh
index 3f9d62f11b956db440e2f603443177cf5c64e2cc..3e5e8bae55054b8059afc1711a22393639cd8d78 100644
--- a/dumux/discretization/elementsolution.hh
+++ b/dumux/discretization/elementsolution.hh
@@ -24,6 +24,11 @@
#ifndef DUMUX_DISCRETIZATION_ELEMENT_SOLUTION_HH
#define DUMUX_DISCRETIZATION_ELEMENT_SOLUTION_HH
+#include
+
+#include
+#include
+
#include
#include
#include
@@ -32,6 +37,64 @@ namespace Dumux {
struct EmptyElementSolution {};
+namespace Experimental {
+
+/*!
+ * \ingroup Discretization
+ * \brief State class to represent the element-local state of the primary
+ * variables of a numerical model, consisting of the spatial degrees
+ * of freedom and potentially a corresponding time level.
+ * \note Here, we extend the concept of elementSolution to additionally carry
+ * time information. Therefore, for now we inherit from the given element
+ * solution, extending its interface
+ */
+template
+class ElementSolution : public BaseElemSol
+{
+public:
+ using TimeLevel = TL;
+
+ //! Constructor
+ template, TL>, int> = 0>
+ ElementSolution(BaseElemSol&& base,
+ TheTimeLevel&& timeLevel)
+ : BaseElemSol(std::forward(base))
+ , timeLevel_(timeLevel)
+ {
+ static_assert(std::is_lvalue_reference_v,
+ "Time level must be an lvalue reference");
+ }
+
+ //! return the time level
+ const TimeLevel& timeLevel() const
+ { return timeLevel_; }
+
+private:
+ const TimeLevel& timeLevel_;
+};
+
+// hide deduction guides from doxygen
+#ifndef DOXYGEN
+template ElementSolution(const B& b, const TL& t) -> ElementSolution;
+template ElementSolution(B&& b, const TL& t) -> ElementSolution;
+#endif // DOXYGEN
+
+/*!
+ * \ingroup Discretization
+ * \brief Make an element solution from a global solution state.
+ */
+template