diff --git a/dumux/assembly/boxlocalassembler.hh b/dumux/assembly/boxlocalassembler.hh
index add56b133acbadefd1c7584b2ebe0870d8ccc091..9487428f95ef5381fd0ae9118cd6f782b80f256e 100644
--- a/dumux/assembly/boxlocalassembler.hh
+++ b/dumux/assembly/boxlocalassembler.hh
@@ -31,177 +31,134 @@
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/assembly/diffmethod.hh>
+#include <dumux/assembly/fvlocalassemblerbase.hh>
namespace Dumux {
/*!
* \ingroup Assembly
* \ingroup BoxDiscretization
- * \brief An assembler for Jacobian and residual contribution per element (box method)
- * \tparam TypeTag the TypeTag
- * \tparam DM the differentiation method to residual compute derivatives
- * \tparam implicit if to use an implicit or explicit time discretization
- */
-template<class TypeTag,
- DiffMethod DM = DiffMethod::numeric,
- bool implicit = true>
-class BoxLocalAssembler;
-
-/*!
- * \ingroup Assembly
- * \ingroup BoxDiscretization
- * \brief Box local assembler using numeric differentiation and implicit time discretization
+ * \brief A base class for all local box 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 BoxLocalAssembler<TypeTag,
- DiffMethod::numeric,
- /*implicit=*/true>
+template<class TypeTag, class Assembler, class Implementation, bool implicit>
+class BoxLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler, Implementation, implicit>
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementResidualVector = Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, NumEqVector)>;
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
+ using ParentType = FVLocalAssemblerBase<TypeTag, Assembler, Implementation, implicit>;
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
- static constexpr int dim = GridView::dimension;
-
public:
+ using ParentType::ParentType;
+
+ void bindLocalViews()
+ {
+ ParentType::bindLocalViews();
+ this->elemBcTypes().update(this->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 Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
+ void assembleJacobianAndResidual(JacobianMatrix& jac, SolutionVector& res, GridVariables& gridVariables)
{
- assemble_(assembler, jac, res, element, curSol);
+ this->asImp_().bindLocalViews();
+
+ const auto residual = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables);
+
+ for (const auto& scv : scvs(this->fvGeometry()))
+ res[scv.dofIndex()] += residual[scv.indexInElement()];
+
+ 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];
+ for (const auto& scvJ : scvs(this->fvGeometry()))
+ {
+ jac[scvI.dofIndex()][scvJ.dofIndex()][eqIdx] = 0.0;
+ if (scvI.indexInElement() == scvJ.indexInElement())
+ jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
+ }
+ };
+ this->asImp_().evalDirichletBoundaries(applyDirichlet);
}
/*!
* \brief Computes the derivatives with respect to the given element and adds them
* to the global matrix.
*/
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac,
- const Element& element, const SolutionVector& curSol)
+ void assembleJacobian(JacobianMatrix& jac, GridVariables& gridVariables)
{
- auto dummyresidual = curSol;
- assemble_(assembler, jac, dummyresidual, element, curSol);
- }
+ this->asImp_().bindLocalViews();
+ this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
- /*!
- * \brief Assemble the residual only
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, res, element, curSol);
- }
+ auto applyDirichlet = [&] (const auto& scvI,
+ const auto& dirichletValues,
+ const auto eqIdx,
+ const auto pvIdx)
+ {
+ for (const auto& scvJ : scvs(this->fvGeometry()))
+ {
+ jac[scvI.dofIndex()][scvJ.dofIndex()][eqIdx] = 0.0;
+ if (scvI.indexInElement() == scvJ.indexInElement())
+ jac[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
+ }
+ };
- /*!
- * \brief Computes the epsilon used for numeric differentiation
- * for a given value of a primary variable.
- *
- * \param priVar The value of the primary variable
- */
- static Scalar numericEpsilon(const Scalar priVar)
- {
- // define the base epsilon as the geometric mean of 1 and the
- // resolution of the scalar type. E.g. for standard 64 bit
- // floating point values, the resolution is about 10^-16 and
- // the base epsilon is thus approximately 10^-8.
- /*
- static const Scalar baseEps
- = Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
- */
- static const Scalar baseEps = 1e-10;
- assert(std::numeric_limits<Scalar>::epsilon()*1e4 < baseEps);
- // the epsilon value used for the numeric differentiation is
- // now scaled by the absolute value of the primary variable...
- return baseEps*(std::abs(priVar) + 1.0);
+ this->asImp_().evalDirichletBoundaries(applyDirichlet);
}
-private:
/*!
- * \brief Computes the residual
- *
- * \return The element residual at the current solution.
+ * \brief Assemble the residual only
*/
- template<class Assembler>
- static void assemble_(Assembler& assembler, SolutionVector& r,
- const Element& element, const SolutionVector& curSol)
+ void assembleResidual(SolutionVector& res)
{
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, curSol);
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+ this->asImp_().bindLocalViews();
+ const auto residual = this->evalLocalResidual();
- const bool isStationary = localResidual.isStationary();
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- if (!isStationary)
- prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
+ for (const auto& scv : scvs(this->fvGeometry()))
+ res[scv.dofIndex()] += residual[scv.indexInElement()];
- // the element boundary types
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // the actual element's current residual
- ElementResidualVector residual(0.0);
- if (isStationary)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
- }
- else
+ auto applyDirichlet = [&] (const auto& scvI,
+ const auto& dirichletValues,
+ const auto eqIdx,
+ const auto pvIdx)
{
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
- }
+ res[scvI.dofIndex()][eqIdx] = this->curElemVolVars()[scvI].priVars()[pvIdx] - dirichletValues[pvIdx];
+ };
- for (const auto& scv : scvs(fvGeometry))
- r[scv.dofIndex()] += residual[scv.indexInElement()];
+ this->asImp_().evalDirichletBoundaries(applyDirichlet);
+ }
- // enforce Dirichlet boundaries by setting the residual to (privar - dirichletvalue)
- if (elemBcTypes.hasDirichlet())
+ /*!
+ * \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(fvGeometry))
+ for (const auto& scvI : scvs(this->fvGeometry()))
{
- const auto bcTypes = elemBcTypes[scvI.indexInElement()];
+ const auto bcTypes = this->elemBcTypes()[scvI.indexInElement()];
if (bcTypes.hasDirichlet())
{
- const auto dirichletValues = problem.dirichlet(element, scvI);
+ const auto dirichletValues = this->problem().dirichlet(this->element(), scvI);
// set the dirichlet conditions in residual and jacobian
for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
@@ -210,9 +167,7 @@ private:
{
const auto pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
assert(0 <= pvIdx && pvIdx < numEq);
-
- const auto& priVars = curElemVolVars[scvI].priVars();
- r[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
+ applyDirichlet(scvI, dirichletValues, eqIdx, pvIdx);
}
}
}
@@ -220,65 +175,62 @@ private:
}
}
+};
+
+/*!
+ * \ingroup Assembly
+ * \ingroup BoxDiscretization
+ * \brief An assembler for Jacobian and residual contribution per element (box methods)
+ * \tparam TypeTag The TypeTag
+ * \tparam DM The differentiation method to residual compute derivatives
+ * \tparam implicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
+ */
+template<class TypeTag, class Assembler, DiffMethod DM = DiffMethod::numeric, bool implicit = true>
+class BoxLocalAssembler;
+
+/*!
+ * \ingroup Assembly
+ * \ingroup BoxDiscretization
+ * \brief Box local assembler using numeric differentiation and implicit time discretization
+ */
+template<class TypeTag, class Assembler>
+class BoxLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>
+: public BoxLocalAssemblerBase<TypeTag, Assembler,
+ BoxLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, true>, true>
+{
+ using ThisType = BoxLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, true>;
+ using ParentType = BoxLocalAssemblerBase<TypeTag, Assembler, ThisType, true>;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ enum { dim = GET_PROP_TYPE(TypeTag, GridView)::dimension };
+
+ static constexpr bool enableGridFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableGridFluxVariablesCache);
+
+public:
+
+ 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 Assembler>
- static void assemble_(Assembler& assembler, JacobianMatrix& A, SolutionVector& r,
- const Element& element, const SolutionVector& curSol)
+ auto assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
{
- // get some references for convenience
- const auto& problem = assembler.problem();
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // prepare the local views
- auto fvGeometry = localView(fvGridGeometry);
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, curSol);
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
-
- const bool isStationary = localResidual.isStationary();
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- if (!isStationary)
- prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
-
- // check for boundaries on the element
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // create the element solution
- ElementSolutionVector elemSol(element, curSol, fvGeometry);
+ // get some aliases for convenience
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ const auto& curSol = this->curSol();
+ auto&& curElemVolVars = this->curElemVolVars();
- // the actual element's current residual
- ElementResidualVector residual(0.0);
- if (isStationary)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
- }
- else
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
- }
+ // get the vecor of the acutal element residuals
+ const auto origResiduals = this->evalLocalResidual();
//////////////////////////////////////////////////////////////////////////////////////////////////
// //
@@ -288,99 +240,35 @@ private:
// //
//////////////////////////////////////////////////////////////////////////////////////////////////
- static const int numericDifferenceMethod = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.NumericDifferenceMethod");
+ // create the element solution
+ ElementSolutionVector elemSol(element, curSol, fvGeometry);
+
+ using ElementResidualVector = std::decay_t<decltype(origResiduals)>;
+ ElementResidualVector partialDerivs(element.subEntities(dim));
// calculation of the derivatives
for (auto&& scv : scvs(fvGeometry))
{
// dof index and corresponding actual pri vars
const auto dofIdx = scv.dofIndex();
- auto& volVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
- VolumeVariables origVolVars(volVars);
-
- // add precalculated residual for this scv into the global container
- r[dofIdx] += residual[scv.indexInElement()];
+ 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++)
{
- ElementResidualVector partialDeriv(element.subEntities(dim));
- Scalar eps = numericEpsilon(volVars.priVar(pvIdx));
- Scalar delta = 0;
+ partialDerivs = 0.0;
- // calculate the residual with the forward deflected primary variables
- if (numericDifferenceMethod >= 0)
+ auto evalResiduals = [&](Scalar priVar)
{
- // we are not using backward differences, i.e. we need to
- // calculate f(x + \epsilon)
-
- // deflect primary variables
- elemSol[scv.indexInElement()][pvIdx] += eps;
- delta += eps;
+ // update the volume variables and compute element residual
+ elemSol[scv.indexInElement()][pvIdx] = priVar;
+ curVolVars.update(elemSol, this->problem(), element, scv);
+ return this->evalLocalResidual();
+ };
- // update the volume variables connected to the dof
- volVars.update(elemSol, problem, element, scv);
-
- // store the deflected residual
- if (isStationary)
- {
- partialDeriv = localResidual.eval(problem, element, fvGeometry,
- curElemVolVars,
- elemBcTypes, elemFluxVarsCache);
- }
- else
- {
- partialDeriv = localResidual.eval(problem, element, fvGeometry,
- prevElemVolVars, curElemVolVars,
- elemBcTypes, elemFluxVarsCache);
- }
- }
- else
- {
- // we are using backward differences, i.e. we don't need
- // to calculate f(x + \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- partialDeriv = residual;
- }
-
- if (numericDifferenceMethod <= 0)
- {
- // we are not using forward differences, i.e. we
- // need to calculate f(x - \epsilon)
-
- // deflect the primary variables
- elemSol[scv.indexInElement()][pvIdx] -= delta + eps;
- delta += eps;
-
- // update the volume variables connected to the dof
- volVars.update(elemSol, problem, element, scv);
-
- // subtract the deflected residual from the derivative storage
- // store the deflected residual
- if (isStationary)
- {
- partialDeriv -= localResidual.eval(problem, element, fvGeometry,
- curElemVolVars,
- elemBcTypes, elemFluxVarsCache);
- }
- else
- {
- partialDeriv -= localResidual.eval(problem, element, fvGeometry,
- prevElemVolVars, curElemVolVars,
- elemBcTypes, elemFluxVarsCache);
- }
- }
- else
- {
- // we are using forward differences, i.e. we don't need to
- // calculate f(x - \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- partialDeriv -= residual;
- }
-
- // divide difference in residuals by the magnitude of the
- // deflections between the two function evaluation
- partialDeriv /= delta;
+ // derive the residuals numerically
+ NumericDifferentiation::partialDerivative(evalResiduals, elemSol[scv.indexInElement()][pvIdx], partialDerivs, origResiduals);
// update the global stiffness matrix with the current partial derivatives
for (auto&& scvJ : scvs(fvGeometry))
@@ -391,296 +279,64 @@ private:
// the residual of equation 'eqIdx' at dof 'i'
// depending on the primary variable 'pvIdx' at dof
// 'col'.
- A[scvJ.dofIndex()][dofIdx][eqIdx][pvIdx] += partialDeriv[scvJ.indexInElement()][eqIdx];
+ A[scvJ.dofIndex()][dofIdx][eqIdx][pvIdx] += partialDerivs[scvJ.indexInElement()][eqIdx];
}
}
// restore the original state of the scv's volume variables
- volVars = origVolVars;
+ curVolVars = origVolVars;
// restore the original element solution
elemSol[scv.indexInElement()][pvIdx] = curSol[scv.dofIndex()][pvIdx];
- }
-
- // TODO additional dof dependencies
- }
-
- // enforce Dirichlet boundaries by overwriting partial derivatives with 1 or 0
- // and set the residual to (privar - dirichletvalue)
- if (elemBcTypes.hasDirichlet())
- {
- for (const auto& scvI : scvs(fvGeometry))
- {
- const auto bcTypes = elemBcTypes[scvI.indexInElement()];
- if (bcTypes.hasDirichlet())
- {
- const auto dirichletValues = problem.dirichlet(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);
-
- const auto& priVars = curElemVolVars[scvI].priVars();
- r[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
- for (const auto& scvJ : scvs(fvGeometry))
- {
- A[scvI.dofIndex()][scvJ.dofIndex()][eqIdx] = 0.0;
- if (scvI.indexInElement() == scvJ.indexInElement())
- A[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
- }
- }
- }
- }
+ // TODO additional dof dependencies
}
}
+ return origResiduals;
}
-private:
- template<class T = TypeTag>
- static typename std::enable_if<!GET_PROP_VALUE(T, EnableGridVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return elemVolVars[scv]; }
-
- template<class T = TypeTag>
- static typename std::enable_if<GET_PROP_VALUE(T, EnableGridVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return gridVolVars.volVars(scv.elementIndex(), scv.indexInElement()); }
}; // implicit BoxAssembler with numeric Jacobian
-
/*!
* \ingroup Assembly
* \ingroup BoxDiscretization
* \brief Box local assembler using numeric differentiation and explicit time discretization
*/
-template<class TypeTag>
-class BoxLocalAssembler<TypeTag,
- DiffMethod::numeric,
- /*implicit=*/false>
+template<class TypeTag, class Assembler>
+class BoxLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false>
+: public BoxLocalAssemblerBase<TypeTag, Assembler,
+ BoxLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, false>, false>
{
+ using ThisType = BoxLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, false>;
+ using ParentType = BoxLocalAssemblerBase<TypeTag, Assembler, ThisType, false>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementResidualVector = Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, NumEqVector)>;
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
-
- static constexpr int dim = GridView::dimension;
+ enum { dim = GET_PROP_TYPE(TypeTag, GridView)::dimension };
public:
- /*!
- * \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 Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, jac, res, element, curSol);
-
- }
+ using ParentType::ParentType;
/*!
* \brief Computes the derivatives with respect to the given element and adds them
* to the global matrix.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac,
- const Element& element, const SolutionVector& curSol)
- {
- auto dummyresidual = curSol;
- assemble_(assembler, jac, dummyresidual, element, curSol);
- }
-
- /*!
- * \brief Assemble the residual only
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, res, element, curSol);
- }
-
- /*!
- * \brief Computes the epsilon used for numeric differentiation
- * for a given value of a primary variable.
- *
- * \param priVar The value of the primary variable
- */
- static Scalar numericEpsilon(const Scalar priVar)
- {
- // define the base epsilon as the geometric mean of 1 and the
- // resolution of the scalar type. E.g. for standard 64 bit
- // floating point values, the resolution is about 10^-16 and
- // the base epsilon is thus approximately 10^-8.
- /*
- static const Scalar baseEps
- = Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
- */
- static const Scalar baseEps = 1e-10;
- assert(std::numeric_limits<Scalar>::epsilon()*1e4 < baseEps);
- // the epsilon value used for the numeric differentiation is
- // now scaled by the absolute value of the primary variable...
- return baseEps*(std::abs(priVar) + 1.0);
- }
-
-private:
- /*!
- * \brief Computes the residual
*
* \return The element residual at the current solution.
*/
- template<class Assembler>
- static void assemble_(Assembler& assembler, SolutionVector& r,
- const Element& element, const SolutionVector& curSol)
+ auto assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
{
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
+ // get some aliases for convenience
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ const auto& curSol = this->curSol();
+ auto&& curElemVolVars = this->curElemVolVars();
- // the element boundary types
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // the actual element's current residual
- auto residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
-
- auto storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
- residual += storageResidual;
-
- for (const auto& scv : scvs(fvGeometry))
- r[scv.dofIndex()] += residual[scv.indexInElement()];
-
- // enforce Dirichlet boundaries by setting the residual to (privar - dirichletvalue)
- if (elemBcTypes.hasDirichlet())
- {
- for (const auto& scvI : scvs(fvGeometry))
- {
- const auto bcTypes = elemBcTypes[scvI.indexInElement()];
- if (bcTypes.hasDirichlet())
- {
- const auto dirichletValues = problem.dirichlet(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);
-
- const auto& priVars = curElemVolVars[scvI].priVars();
- r[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
- }
- }
- }
- }
- }
- }
-
- /*!
- * \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 Assembler>
- static void assemble_(Assembler& assembler, JacobianMatrix& A, SolutionVector& r,
- const Element& element, const SolutionVector& curSol)
- {
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
-
- // the element boundary types
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // get the element solution
- const auto numVert = element.subEntities(dim);
- ElementSolutionVector elemSol(numVert);
- for (const auto& scv : scvs(fvGeometry))
- elemSol[scv.indexInElement()] = localResidual.prevSol()[scv.dofIndex()];
-
- // the actual element's previous time step residual
- auto residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
-
- auto storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
-
- residual += storageResidual;
+ // get the vecor of the acutal element residuals
+ const auto origResiduals = this->evalLocalResidual();
//////////////////////////////////////////////////////////////////////////////////////////////////
// //
@@ -690,78 +346,35 @@ private:
// //
//////////////////////////////////////////////////////////////////////////////////////////////////
- static const int numericDifferenceMethod = getParamFromGroup<int>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.NumericDifferenceMethod");
+ // create the element solution
+ ElementSolutionVector elemSol(element, curSol, fvGeometry);
+
+ using ElementResidualVector = std::decay_t<decltype(origResiduals)>;
+ ElementResidualVector partialDerivs(element.subEntities(dim));
// calculation of the derivatives
for (auto&& scv : scvs(fvGeometry))
{
// dof index and corresponding actual pri vars
const auto dofIdx = scv.dofIndex();
- auto& volVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
- VolumeVariables origVolVars(volVars);
-
- // add precalculated residual for this scv into the global container
- r[dofIdx] += residual[scv.indexInElement()];
+ 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++)
{
- ElementSolutionVector partialDeriv(element.subEntities(dim));
- Scalar eps = numericEpsilon(volVars.priVar(pvIdx));
- Scalar delta = 0;
-
- // calculate the residual with the forward deflected primary variables
- if (numericDifferenceMethod >= 0)
- {
- // we are not using backward differences, i.e. we need to
- // calculate f(x + \epsilon)
-
- // deflect primary variables
- elemSol[scv.indexInElement()][pvIdx] += eps;
- delta += eps;
-
- // update the volume variables connected to the dof
- volVars.update(elemSol, problem, element, scv);
+ partialDerivs = 0.0;
- partialDeriv = localResidual.evalStorage(problem, element, fvGeometry,
- prevElemVolVars, curElemVolVars,
- elemBcTypes, elemFluxVarsCache);
- }
- else
- {
- // we are using backward differences, i.e. we don't need
- // to calculate f(x + \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- partialDeriv = residual;
- }
-
- if (numericDifferenceMethod <= 0)
+ auto evalStorage = [&](Scalar priVar)
{
- // we are not using forward differences, i.e. we
- // need to calculate f(x - \epsilon)
+ // auto partialDerivsTmp = partialDerivs;
+ elemSol[scv.indexInElement()][pvIdx] = priVar;
+ curVolVars.update(elemSol, this->problem(), element, scv);
+ return this->evalLocalStorageResidual();
+ };
- // deflect the primary variables
- elemSol[scv.indexInElement()][pvIdx] -= delta + eps;
- delta += eps;
-
- // update the volume variables connected to the dof
- volVars.update(elemSol, problem, element, scv);
-
- partialDeriv -= localResidual.evalStorage(problem, element, fvGeometry,
- prevElemVolVars, curElemVolVars,
- elemBcTypes, elemFluxVarsCache);
- }
- else
- {
- // we are using forward differences, i.e. we don't need to
- // calculate f(x - \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- partialDeriv -= residual;
- }
-
- // divide difference in residuals by the magnitude of the
- // deflections between the two function evaluation
- partialDeriv /= delta;
+ // derive the residuals numerically
+ NumericDifferentiation::partialDerivative(evalStorage, elemSol[scv.indexInElement()][pvIdx], partialDerivs, origResiduals);
// update the global stiffness matrix with the current partial derivatives
for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
@@ -770,210 +383,39 @@ private:
// the residual of equation 'eqIdx' at dof 'i'
// depending on the primary variable 'pvIdx' at dof
// 'col'.
- A[dofIdx][dofIdx][eqIdx][pvIdx] += partialDeriv[scv.indexInElement()][eqIdx];
+ A[dofIdx][dofIdx][eqIdx][pvIdx] += partialDerivs[scv.indexInElement()][eqIdx];
}
// restore the original state of the scv's volume variables
- volVars = origVolVars;
+ curVolVars = origVolVars;
// restore the original element solution
elemSol[scv.indexInElement()][pvIdx] = curSol[scv.dofIndex()][pvIdx];
- }
-
- // TODO additional dof dependencies
- }
-
- // enforce Dirichlet boundaries by overwriting partial derivatives with 1 or 0
- // and set the residual to (privar - dirichletvalue)
- if (elemBcTypes.hasDirichlet())
- {
- for (const auto& scvI : scvs(fvGeometry))
- {
- const auto bcTypes = elemBcTypes[scvI.indexInElement()];
- if (bcTypes.hasDirichlet())
- {
- const auto dirichletValues = problem.dirichlet(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);
-
- const auto& priVars = curElemVolVars[scvI].priVars();
- r[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
- A[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
- }
- }
- }
+ // TODO additional dof dependencies
}
}
+ return origResiduals;
}
-private:
- template<class T = TypeTag>
- static typename std::enable_if<!GET_PROP_VALUE(T, EnableGridVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return elemVolVars[scv]; }
-
- template<class T = TypeTag>
- static typename std::enable_if<GET_PROP_VALUE(T, EnableGridVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return gridVolVars.volVars(scv.elementIndex(), scv.indexInElement()); }
-
}; // explicit BoxAssembler with numeric Jacobian
-
/*!
* \ingroup Assembly
* \ingroup BoxDiscretization
- * \brief Box local assembler using analytic (hard coded) derivatives and implicit time discretization
+ * \brief Box local assembler using analytic differentiation and implicit time discretization
*/
-template<class TypeTag>
-class BoxLocalAssembler<TypeTag,
- DiffMethod::analytic,
- /*implicit=*/true>
+template<class TypeTag, class Assembler>
+class BoxLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/true>
+: public BoxLocalAssemblerBase<TypeTag, Assembler,
+ BoxLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, true>, true>
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementResidualVector = Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, NumEqVector)>;
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using ThisType = BoxLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, true>;
+ using ParentType = BoxLocalAssemblerBase<TypeTag, Assembler, ThisType, true>;
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
-
- static constexpr int dim = GridView::dimension;
-
public:
- /*!
- * \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 Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, jac, res, element, curSol);
-
- }
-
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac,
- const Element& element, const SolutionVector& curSol)
- {
- auto dummyresidual = curSol;
- assemble_(assembler, jac, dummyresidual, element, curSol);
- }
-
- /*!
- * \brief Assemble the residual only
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, res, element, curSol);
- }
-
-private:
- /*!
- * \brief Computes the residual
- *
- * \return The element residual at the current solution.
- */
- template<class Assembler>
- static void assemble_(Assembler& assembler, SolutionVector& r,
- const Element& element, const SolutionVector& curSol)
- {
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, curSol);
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
-
- const bool isStationary = localResidual.isStationary();
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- if (!isStationary)
- prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
-
- // the element boundary types
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // the actual element's current residual
- ElementResidualVector residual(0.0);
- if (isStationary)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
- }
- else
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
- }
-
- for (const auto& scv : scvs(fvGeometry))
- r[scv.dofIndex()] += residual[scv.indexInElement()];
-
- // enforce Dirichlet boundaries by setting the residual to (privar - dirichletvalue)
- if (elemBcTypes.hasDirichlet())
- {
- for (const auto& scvI : scvs(fvGeometry))
- {
- const auto bcTypes = elemBcTypes[scvI.indexInElement()];
- if (bcTypes.hasDirichlet())
- {
- const auto dirichletValues = problem.dirichlet(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);
-
- const auto& priVars = curElemVolVars[scvI].priVars();
- r[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
- }
- }
- }
- }
- }
- }
+ using ParentType::ParentType;
/*!
* \brief Computes the derivatives with respect to the given element and adds them
@@ -981,59 +423,17 @@ private:
*
* \return The element residual at the current solution.
*/
- template<class Assembler>
- static void assemble_(Assembler& assembler, JacobianMatrix& A, SolutionVector& r,
- const Element& element, const SolutionVector& curSol)
+ auto assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
{
- // get some references for convenience
- const auto& problem = assembler.problem();
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // prepare the local views
- auto fvGeometry = localView(fvGridGeometry);
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, curSol);
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
-
- const bool isStationary = localResidual.isStationary();
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- if (!isStationary)
- prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
-
- // check for boundaries on the element
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // the actual element's current residual
- ElementResidualVector residual(0.0);
- if (isStationary)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
- }
- else
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
- }
+ // get some aliases for convenience
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ const auto& problem = this->problem();
+ auto&& curElemVolVars = this->curElemVolVars();
+ auto&& elemFluxVarsCache = this->elemFluxVarsCache();
- for (auto&& scv : scvs(fvGeometry))
- r[scv.dofIndex()] += residual[scv.indexInElement()];
+ // get the vecor of the acutal element residuals
+ const auto origResiduals = this->evalLocalResidual();
//////////////////////////////////////////////////////////////////////////////////////////////////
// //
@@ -1052,22 +452,23 @@ private:
// 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
- if (!isStationary)
- localResidual.addStorageDerivatives(A[dofIdx][dofIdx],
- problem,
- element,
- fvGeometry,
- volVars,
- scv);
+ // 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
- localResidual.addSourceDerivatives(A[dofIdx][dofIdx],
- problem,
- element,
- fvGeometry,
- volVars,
- scv);
+ this->localResidual().addSourceDerivatives(A[dofIdx][dofIdx],
+ problem,
+ element,
+ fvGeometry,
+ volVars,
+ scv);
}
// localJacobian[scvIdx][otherScvIdx][eqIdx][priVarIdx] of the fluxes
@@ -1076,13 +477,13 @@ private:
if (!scvf.boundary())
{
// add flux term derivatives
- localResidual.addFluxDerivatives(A,
- problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- scvf);
+ this->localResidual().addFluxDerivatives(A,
+ problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ scvf);
}
// the boundary gets special treatment to simplify
@@ -1090,261 +491,60 @@ private:
else
{
const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
- if (elemBcTypes[insideScv.indexInElement()].hasNeumann())
+ if (this->elemBcTypes()[insideScv.indexInElement()].hasNeumann())
{
// add flux term derivatives
- localResidual.addRobinFluxDerivatives(A[insideScv.dofIndex()],
- problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- scvf);
+ this->localResidual().addRobinFluxDerivatives(A[insideScv.dofIndex()],
+ problem,
+ element,
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ scvf);
}
}
}
- // enforce Dirichlet boundaries by overwriting partial derivatives with 1 or 0
- // and set the residual to (privar - dirichletvalue)
- if (elemBcTypes.hasDirichlet())
- {
- for (const auto& scvI : scvs(fvGeometry))
- {
- const auto bcTypes = elemBcTypes[scvI.indexInElement()];
- if (bcTypes.hasDirichlet())
- {
- const auto dirichletValues = problem.dirichlet(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);
-
- const auto& priVars = curElemVolVars[scvI].priVars();
- r[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
- for (const auto& scvJ : scvs(fvGeometry))
- {
- A[scvI.dofIndex()][scvJ.dofIndex()][eqIdx] = 0.0;
- if (scvI.indexInElement() == scvJ.indexInElement())
- A[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
- }
- }
- }
- }
- }
- }
-
- // TODO additional dof dependencies
+ return origResiduals;
}
}; // implicit BoxAssembler with analytic Jacobian
-
/*!
* \ingroup Assembly
* \ingroup BoxDiscretization
- * \brief Box local assembler using analytic (hard coded) derivatives and explicit time discretization
+ * \brief Box local assembler using analytic differentiation and explicit time discretization
*/
-template<class TypeTag>
-class BoxLocalAssembler<TypeTag,
- DiffMethod::analytic,
- /*implicit=*/false>
+template<class TypeTag, class Assembler>
+class BoxLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/false>
+: public BoxLocalAssemblerBase<TypeTag, Assembler,
+ BoxLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, false>, false>
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementResidualVector = Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, NumEqVector)>;
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using ThisType = BoxLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, false>;
+ using ParentType = BoxLocalAssemblerBase<TypeTag, Assembler, ThisType, false>;
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
-
- static constexpr int dim = GridView::dimension;
-
public:
- /*!
- * \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 Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, jac, res, element, curSol);
-
- }
+ using ParentType::ParentType;
/*!
* \brief Computes the derivatives with respect to the given element and adds them
* to the global matrix.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac,
- const Element& element, const SolutionVector& curSol)
- {
- auto dummyresidual = curSol;
- assemble_(assembler, jac, dummyresidual, element, curSol);
- }
-
- /*!
- * \brief Assemble the residual only
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, res, element, curSol);
- }
-
-private:
- /*!
- * \brief Computes the residual
*
* \return The element residual at the current solution.
*/
- template<class Assembler>
- static void assemble_(Assembler& assembler, SolutionVector& r,
- const Element& element, const SolutionVector& curSol)
+ auto assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
{
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
+ // get some aliases for convenience
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ const auto& problem = this->problem();
+ auto&& curElemVolVars = this->curElemVolVars();
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
-
- // the element boundary types
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // the actual element's current residual
- auto residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
-
- auto storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
- residual += storageResidual;
-
- for (const auto& scv : scvs(fvGeometry))
- r[scv.dofIndex()] += residual[scv.indexInElement()];
-
- // enforce Dirichlet boundaries by setting the residual to (privar - dirichletvalue)
- if (elemBcTypes.hasDirichlet())
- {
- for (const auto& scvI : scvs(fvGeometry))
- {
- const auto bcTypes = elemBcTypes[scvI.indexInElement()];
- if (bcTypes.hasDirichlet())
- {
- const auto dirichletValues = problem.dirichlet(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);
-
- const auto& priVars = curElemVolVars[scvI].priVars();
- r[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
- }
- }
- }
- }
- }
- }
-
- /*!
- * \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 Assembler>
- static void assemble_(Assembler& assembler, JacobianMatrix& A, SolutionVector& r,
- const Element& element, const SolutionVector& curSol)
- {
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
-
- // the element boundary types
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // the actual element's current residual
- auto residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
-
- auto storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache);
- residual += storageResidual;
-
- for (const auto& scv : scvs(fvGeometry))
- r[scv.dofIndex()] += residual[scv.indexInElement()];
+ // get the vecor of the acutal element residuals
+ const auto origResiduals = this->evalLocalResidual();
//////////////////////////////////////////////////////////////////////////////////////////////////
// //
@@ -1363,43 +563,15 @@ private:
// 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
- localResidual.addStorageDerivatives(A[dofIdx][dofIdx],
- problem,
- element,
- fvGeometry,
- volVars,
- scv);
- }
-
- // enforce Dirichlet boundaries by overwriting partial derivatives with 1 or 0
- // and set the residual to (privar - dirichletvalue)
- if (elemBcTypes.hasDirichlet())
- {
- for (const auto& scvI : scvs(fvGeometry))
- {
- const auto bcTypes = elemBcTypes[scvI.indexInElement()];
- if (bcTypes.hasDirichlet())
- {
- const auto dirichletValues = problem.dirichlet(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);
-
- const auto& priVars = curElemVolVars[scvI].priVars();
- r[scvI.dofIndex()][eqIdx] = priVars[pvIdx] - dirichletValues[pvIdx];
- A[scvI.dofIndex()][scvI.dofIndex()][eqIdx][pvIdx] = 1.0;
- }
- }
- }
- }
+ this->localResidual().addStorageDerivatives(A[dofIdx][dofIdx],
+ problem,
+ element,
+ fvGeometry,
+ volVars,
+ scv);
}
- // TODO additional dof dependencies
+ return origResiduals;
}
}; // explicit BoxAssembler with analytic Jacobian
diff --git a/dumux/assembly/boxlocalresidual.hh b/dumux/assembly/boxlocalresidual.hh
index 7beb17dc02973953b3fa9d6b6b9f907e5911d56d..2858c752202145edb07110893489d12c8647f21f 100644
--- a/dumux/assembly/boxlocalresidual.hh
+++ b/dumux/assembly/boxlocalresidual.hh
@@ -53,9 +53,9 @@ class BoxLocalResidual : public FVLocalResidual<TypeTag>
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementResidualVector = Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, NumEqVector)>;
public:
+ using ElementResidualVector = typename ParentType::ElementResidualVector;
using ParentType::ParentType;
//! evaluate flux residuals for one sub control volume face and add to residual
diff --git a/dumux/assembly/cclocalassembler.hh b/dumux/assembly/cclocalassembler.hh
index 9d5be567a28ee279d7a5379c8ebc2116413defc3..00d9c51e2d82bf7b0d11df075241cce916f99710 100644
--- a/dumux/assembly/cclocalassembler.hh
+++ b/dumux/assembly/cclocalassembler.hh
@@ -21,180 +21,121 @@
* \ingroup Assembly
* \ingroup CCDiscretization
* \brief An assembler for Jacobian and residual contribution per element (cell-centered methods)
- * \tparam TypeTag the TypeTag
- * \tparam DM the differentiation method to residual compute derivatives
- * \tparam implicit if to use an implicit or explicit time discretization
*/
#ifndef DUMUX_CC_LOCAL_ASSEMBLER_HH
#define DUMUX_CC_LOCAL_ASSEMBLER_HH
+#include <dune/common/reservedvector.hh>
#include <dune/grid/common/gridenums.hh> // for GhostEntity
#include <dune/istl/matrixindexset.hh>
-#include <dune/istl/bvector.hh>
+#include <dumux/common/reservedblockvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/assembly/diffmethod.hh>
+#include <dumux/assembly/numericdifferentiation.hh>
+#include <dumux/assembly/fvlocalassemblerbase.hh>
+#include <dumux/discretization/fluxstencil.hh>
namespace Dumux {
/*!
* \ingroup Assembly
* \ingroup CCDiscretization
- * \brief An assembler for Jacobian and residual contribution per element (cell-centered methods)
- * \tparam TypeTag the TypeTag
- * \tparam DM the differentiation method to residual compute derivatives
- * \tparam implicit if to use an implicit or explicit time discretization
+ * \brief A base class for all local cell-centered assemblers
+ * \tparam TypeTag The TypeTag
+ * \tparam Assembler The assembler type
+ * \tparam Implementation The actual implementation
+ * \tparam implicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
*/
-template<class TypeTag,
- DiffMethod DM = DiffMethod::numeric,
- bool implicit = true>
-class CCLocalAssembler;
-
-/*!
- * \ingroup Assembly
- * \ingroup CCDiscretization
- * \brief Cell-centered scheme local assembler using numeric differentiation and implicit time discretization
- */
-template<class TypeTag>
-class CCLocalAssembler<TypeTag,
- DiffMethod::numeric,
- /*implicit=*/true>
+template<class TypeTag, class Assembler, class Implementation, bool implicit>
+class CCLocalAssemblerBase : public FVLocalAssemblerBase<TypeTag, Assembler, Implementation, implicit>
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
+ using ParentType = FVLocalAssemblerBase<TypeTag, Assembler, Implementation, implicit>;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
-
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
-
- static constexpr bool enableGridFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableGridFluxVariablesCache);
+ using LocalResidualValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
public:
+ using ParentType::ParentType;
+
/*!
* \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 Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
+ void assembleJacobianAndResidual(JacobianMatrix& jac, SolutionVector& res, GridVariables& gridVariables)
{
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, jac, element, curSol);
+ this->asImp_().bindLocalViews();
+ const auto globalI = this->assembler().fvGridGeometry().elementMapper().index(this->element());
+ res[globalI] = this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
}
/*!
* \brief Computes the derivatives with respect to the given element and adds them
* to the global matrix.
*/
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac,
- const Element& element, const SolutionVector& curSol)
+ void assembleJacobian(JacobianMatrix& jac, GridVariables& gridVariables)
{
- assemble_(assembler, jac, element, curSol);
+ this->asImp_().bindLocalViews();
+ this->asImp_().assembleJacobianAndResidualImpl(jac, gridVariables); // forward to the internal implementation
}
/*!
* \brief Assemble the residual only
*/
- template<class Assembler>
- static void assemble(Assembler& assembler, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, element, curSol);
- }
-
- /*!
- * \brief Computes the epsilon used for numeric differentiation
- * for a given value of a primary variable.
- *
- * \param priVar The value of the primary variable
- */
- static Scalar numericEpsilon(const Scalar priVar)
+ void assembleResidual(SolutionVector& res)
{
- // define the base epsilon as the geometric mean of 1 and the
- // resolution of the scalar type. E.g. for standard 64 bit
- // floating point values, the resolution is about 10^-16 and
- // the base epsilon is thus approximately 10^-8.
- /*
- static const Scalar baseEps
- = Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
- */
- static const Scalar baseEps = 1e-10;
- assert(std::numeric_limits<Scalar>::epsilon()*1e4 < baseEps);
- // the epsilon value used for the numeric differentiation is
- // now scaled by the absolute value of the primary variable...
- return baseEps*(std::abs(priVar) + 1.0);
+ this->asImp_().bindLocalViews();
+ const auto globalI = this->assembler().fvGridGeometry().elementMapper().index(this->element());
+ res[globalI] = this->asImp_().evalLocalResidual()[0]; // forward to the internal implementation
}
+};
-private:
- /*!
- * \brief Computes the residual
- *
- * \return The element residual at the current solution.
- */
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler,
- const Element& element, const SolutionVector& curSol)
- {
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- if (isGhost) return NumEqVector(0.0);
-
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, curSol);
+/*!
+ * \ingroup Assembly
+ * \ingroup CCDiscretization
+ * \brief An assembler for Jacobian and residual contribution per element (cell-centered methods)
+ * \tparam TypeTag The TypeTag
+ * \tparam DM The differentiation method to residual compute derivatives
+ * \tparam implicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
+ */
+template<class TypeTag, class Assembler, DiffMethod DM = DiffMethod::numeric, bool implicit = true>
+class CCLocalAssembler;
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+/*!
+ * \ingroup Assembly
+ * \ingroup CCDiscretization
+ * \brief Cell-centered scheme local assembler using numeric differentiation and implicit time discretization
+ */
+template<class TypeTag, class Assembler>
+class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/true>
+: public CCLocalAssemblerBase<TypeTag, Assembler,
+ CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, true>, true >
+{
+ using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, true>;
+ using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, true>;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using LocalResidualValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ enum { dim = GET_PROP_TYPE(TypeTag, GridView)::dimension };
- // for compatibility with box models
- ElementBoundaryTypes elemBcTypes;
+ using FluxStencil = Dumux::FluxStencil<TypeTag>;
+ static constexpr int maxNeighbors = FluxStencil::maxFluxStencilSize*FVElementGeometry::maxNumElementScvfs;
+ static constexpr bool enableGridFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableGridFluxVariablesCache);
- // the actual element's current residual
- NumEqVector residual(0.0);
- if (localResidual.isStationary())
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
+public:
- return residual;
- }
+ using ParentType::ParentType;
/*!
* \brief Computes the derivatives with respect to the given element and adds them
@@ -202,117 +143,64 @@ private:
*
* \return The element residual at the current solution.
*/
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
- const Element& element, const SolutionVector& curSol)
+ LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
{
- // get some references for convenience
- const auto& problem = assembler.problem();
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- const auto& connectivityMap = fvGridGeometry.connectivityMap();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, curSol);
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
-
- // the global dof of the actual element
- const auto globalI = fvGridGeometry.elementMapper().index(element);
-
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- // is the local residual stationary?
- const bool isStationary = localResidual.isStationary();
-
- // the actual element's current residual
- NumEqVector residual(0.0);
- if (!isGhost)
- {
- if (isStationary)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
-
-
- // TODO Do we really need this??????????
- // this->model_().updatePVWeights(fvGeometry);
-
//////////////////////////////////////////////////////////////////////////////////////////////////
- // //
// 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. //
- // //
//////////////////////////////////////////////////////////////////////////////////////////////////
- static const std::string group = GET_PROP_VALUE(TypeTag, ModelParameterGroup);
- static const int numericDifferenceMethod = getParamFromGroup<int>(group, "Implicit.NumericDifferenceMethod");
+ // get some aliases for convenience
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ const auto& fvGridGeometry = this->assembler().fvGridGeometry();
+ auto&& curElemVolVars = this->curElemVolVars();
+ auto&& elemFluxVarsCache = this->elemFluxVarsCache();
// get stencil informations
+ const auto globalI = fvGridGeometry.elementMapper().index(element);
+ const auto& connectivityMap = fvGridGeometry.connectivityMap();
const auto numNeighbors = connectivityMap[globalI].size();
// container to store the neighboring elements
- std::vector<Element> neighborElements;
- neighborElements.reserve(numNeighbors);
+ Dune::ReservedVector<Element, maxNeighbors+1> neighborElements;
+ neighborElements.resize(numNeighbors);
+
+ // assemble the undeflected residual
+ using Residuals = ReservedBlockVector<LocalResidualValues, maxNeighbors+1>;
+ Residuals origResiduals(numNeighbors + 1); origResiduals = 0.0;
+ origResiduals[0] = this->evalLocalResidual()[0];
+
+ // lambda for convenient evaluation of the fluxes across scvfs in the neighbors
+ auto evalNeighborFlux = [&] (const auto& neighbor, const auto& scvf)
+ {
+ return this->localResidual().evalFlux(this->problem(),
+ neighbor,
+ this->fvGeometry(),
+ this->curElemVolVars(),
+ this->elemFluxVarsCache(), scvf);
+ };
// get the elements in which we need to evaluate the fluxes
// and calculate these in the undeflected state
- Dune::BlockVector<NumEqVector> origFlux(numNeighbors);
- origFlux = 0.0;
- unsigned int j = 0;
+ unsigned int j = 1;
for (const auto& dataJ : connectivityMap[globalI])
{
- neighborElements.emplace_back(fvGridGeometry.element(dataJ.globalJ));
+ neighborElements[j-1] = fvGridGeometry.element(dataJ.globalJ);
for (const auto scvfIdx : dataJ.scvfsJ)
- {
- origFlux[j] += localResidual.evalFlux(problem,
- neighborElements.back(),
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- fvGeometry.scvf(scvfIdx));
- }
- // increment neighbor counter
+ 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 = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
+ auto& curVolVars = ParentType::getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
// save a copy of the original privars and vol vars in order
// to restore the original solution after deflection
+ const auto& curSol = this->curSol();
const auto origPriVars = curSol[globalI];
const auto origVolVars = curVolVars;
@@ -320,435 +208,102 @@ private:
ElementSolutionVector elemSol(origPriVars);
// derivatives in the neighbors with repect to the current elements
- Dune::BlockVector<NumEqVector> neighborDeriv(numNeighbors);
- for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
- {
- // reset derivatives of element dof with respect to itself
- // as well as neighbor derivatives
- NumEqVector partialDeriv(0.0);
- neighborDeriv = 0.0;
+ // in index 0 we save the derivative of the element residual with respect to it's own dofs
+ Residuals partialDerivs(numNeighbors + 1);
- if (isGhost)
- partialDeriv[pvIdx] = 1.0;
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
+ {
- Scalar eps = numericEpsilon(curVolVars.priVar(pvIdx));
- Scalar delta = 0;
+ // for ghost elements we assemble a 1.0 where the primary variable and zero everywhere else
+ // as we always solve for a delta of the solution with repect to the initial solution this
+ // results in a delta of zero for ghosts, we still need to do the neighbor derivatives though
+ // so we are not done yet here.
+ partialDerivs = 0.0;
+ if (this->elementIsGhost()) partialDerivs[0][pvIdx] = 1.0;
- if (numericDifferenceMethod >= 0)
+ auto evalResiduals = [&](Scalar priVar)
{
- // we are not using backward differences, i.e. we need to
- // calculate f(x + \epsilon)
-
- // deflect primary variables
- elemSol[0][pvIdx] += eps;
- delta += eps;
-
+ Residuals partialDerivsTmp(numNeighbors + 1);
+ partialDerivsTmp = 0.0;
// update the volume variables and the flux var cache
- curVolVars.update(elemSol, problem, element, scv);
+ elemSol[0][pvIdx] = priVar;
+ curVolVars.update(elemSol, this->problem(), element, scv);
if (enableGridFluxVarsCache)
gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
else
elemFluxVarsCache.update(element, fvGeometry, curElemVolVars);
// calculate the residual with the deflected primary variables
- if (!isGhost)
- {
- if (isStationary)
- {
- partialDeriv = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- partialDeriv = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
+ if (!this->elementIsGhost()) partialDerivsTmp[0] = this->evalLocalResidual()[0];
// calculate the fluxes in the neighbors with the deflected primary variables
for (std::size_t k = 0; k < numNeighbors; ++k)
for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
- {
- neighborDeriv[k] += localResidual.evalFlux(problem,
- neighborElements[k],
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- fvGeometry.scvf(scvfIdx));
- }
- }
- else
- {
- // we are using backward differences, i.e. we don't need
- // to calculate f(x + \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- if (!isGhost)
- partialDeriv = residual;
- neighborDeriv = origFlux;
- }
-
- if (numericDifferenceMethod <= 0)
- {
- // we are not using forward differences, i.e. we
- // need to calculate f(x - \epsilon)
-
- // deflect the primary variables
- elemSol[0][pvIdx] -= delta + eps;
- delta += eps;
-
- // update the volume variables and the flux var cache
- curVolVars.update(elemSol, problem, element, scv);
- if (enableGridFluxVarsCache)
- gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
- else
- elemFluxVarsCache.update(element, fvGeometry, curElemVolVars);
-
- // calculate the residual with the deflected primary variables and subtract it
- if (!isGhost)
- {
- if (isStationary)
- {
- partialDeriv -= localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- partialDeriv -= localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
-
- // calculate the fluxes into element with the deflected primary variables
- for (std::size_t k = 0; k < numNeighbors; ++k)
- for (auto scvfIdx : connectivityMap[globalI][k].scvfsJ)
- {
- neighborDeriv[k] -= localResidual.evalFlux(problem,
- neighborElements[k],
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- fvGeometry.scvf(scvfIdx));
- }
- }
- else
- {
- // we are using forward differences, i.e. we don't need to
- // calculate f(x - \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- if (!isGhost)
- partialDeriv -= residual;
- neighborDeriv -= origFlux;
- }
+ partialDerivsTmp[k+1] += evalNeighborFlux(neighborElements[k], fvGeometry.scvf(scvfIdx));
- // divide difference in residuals by the magnitude of the
- // deflections between the two function evaluation
- if (!isGhost)
- partialDeriv /= delta;
- neighborDeriv /= delta;
+ return partialDerivsTmp;
+ };
- // restore the original state of the scv's volume variables
- curVolVars = origVolVars;
-
- // restore the current element solution
- elemSol[0][pvIdx] = origPriVars[pvIdx];
+ // derive the residuals numerically
+ NumericDifferentiation::partialDerivative(evalResiduals, elemSol[0][pvIdx], partialDerivs, origResiduals);
// add the current partial derivatives to the global jacobian matrix
for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
{
// the diagonal entries
- A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
+ A[globalI][globalI][eqIdx][pvIdx] += partialDerivs[0][eqIdx];
// off-diagonal entries
- j = 0;
+ j = 1;
for (const auto& dataJ : connectivityMap[globalI])
- A[dataJ.globalJ][globalI][eqIdx][pvIdx] += neighborDeriv[j++][eqIdx];
+ 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.
+ // 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
if (enableGridFluxVarsCache)
gridVariables.gridFluxVarsCache().updateElement(element, fvGeometry, curElemVolVars);
- //////////////////////////////////////////////////////////////////////////////////////////////
- // //
- // Calculate derivatives of the dofs in the element with respect to user-defined additional //
- // dof dependencies. We do so by evaluating the change in the source term of the current //
- // element with respect to the primary variables at the given additional dofs. //
- // //
- //////////////////////////////////////////////////////////////////////////////////////////////
-
- // const auto& additionalDofDepedencies = problem.getAdditionalDofDependencies(globalI);
- // if (!additionalDofDepedencies.empty() && !isGhost)
- // {
- // // compute the source in the undeflected state
- // auto source = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
- // source *= -scv.volume()*curVolVarsI.extrusionFactor();
-
- // // deflect solution at given dofs and recalculate the source
- // for (auto globalJ : additionalDofDependencies)
- // {
- // const auto& scvJ = fvGeometry.scv(globalJ);
- // auto& curVolVarsJ = curElemVolVars[scv];
- // const auto& elementJ = fvGridGeometry.element(globalJ);
-
- // // save a copy of the original privars and volvars
- // // to restore original solution after deflection
- // const auto origPriVars = curSol[globalJ];
- // const auto origVolVarsJ = curVolVarsJ;
-
- // // derivatives with repect to the additional DOF we depend on
- // for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
- // {
- // // derivatives of element dof with respect to itself
- // NumEqVector partialDeriv(0.0);
- // const auto eps = numericEpsilon(curVolVarsJ.priVar(pvIdx));
- // Scalar delta = 0;
-
- // if (numericDifferenceMethod >= 0)
- // {
- // // we are not using backward differences, i.e. we need to
- // // calculate f(x + \epsilon)
-
- // // deflect primary variables
- // curSol[globalJ][pvIdx] += eps;
- // delta += eps;
-
- // // update the volume variables and the flux var cache
- // curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
-
- // // calculate the source with the deflected primary variables
- // auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
- // deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
- // partialDeriv = std::move(deflSource);
- // }
- // else
- // {
- // // we are using backward differences, i.e. we don't need
- // // to calculate f(x + \epsilon) and we can recycle the
- // // (already calculated) source f(x)
- // partialDeriv = source;
- // }
-
- // if (numericDifferenceMethod <= 0)
- // {
- // // we are not using forward differences, i.e. we
- // // need to calculate f(x - \epsilon)
-
- // // deflect the primary variables
- // curSol[globalJ][pvIdx] -= delta + eps;
- // delta += eps;
-
- // // update the volume variables and the flux var cache
- // curVolVarsJ.update(gridVariables.elementSolution(elementJ, curSol), problem, elementJ, scvJ);
-
- // // calculate the source with the deflected primary variables and subtract
- // auto deflSource = localResidual.computeSource(element, fvGeometry, curElemVolVars, scv);
- // deflSource *= -scv.volume()*curVolVarsI.extrusionFactor();
- // partialDeriv -= std::move(deflSource);
- // }
- // else
- // {
- // // we are using forward differences, i.e. we don't need to
- // // calculate f(x - \epsilon) and we can recycle the
- // // (already calculated) source f(x)
- // partialDeriv -= source;
- // }
-
- // // divide difference in residuals by the magnitude of the
- // // deflections between the two function evaluation
- // partialDeriv /= delta;
-
- // // restore the original state of the dofs privars and the volume variables
- // curSol[globalJ] = origPriVars;
- // curVolVarsJ = origVolVarsJ;
-
- // // add the current partial derivatives to the global jacobian matrix
- // for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- // A[globalI][globalJ][eqIdx][pvIdx] += partialDeriv[eqIdx];
- // }
- // }
- // }
-
// return the original residual
- return residual;
+ return origResiduals[0];
}
-private:
- template<class T = TypeTag>
- static typename std::enable_if<!GET_PROP_VALUE(T, EnableGridVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return elemVolVars[scv]; }
-
- template<class T = TypeTag>
- static typename std::enable_if<GET_PROP_VALUE(T, EnableGridVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return gridVolVars.volVars(scv); }
};
/*!
* \ingroup Assembly
* \ingroup CCDiscretization
- * \brief Cell-centered scheme local assembler using numeric differentiation and explicits time discretization
+ * \brief Cell-centered scheme local assembler using numeric differentiation and explicit time discretization
*/
-template<class TypeTag>
-class CCLocalAssembler<TypeTag,
- DiffMethod::numeric,
- /*implicit=*/false>
+template<class TypeTag, class Assembler>
+class CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, /*implicit=*/false>
+: public CCLocalAssemblerBase<TypeTag, Assembler,
+ CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, false>, false>
{
+ using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::numeric, false>;
+ using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, false>;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+ using LocalResidualValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
public:
-
- /*!
- * \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 Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, jac, element, curSol);
- }
-
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, jac, element, curSol);
- }
-
- /*!
- * \brief Assemble the residual only
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, element, curSol);
- }
-
- /*!
- * \brief Computes the epsilon used for numeric differentiation
- * for a given value of a primary variable.
- *
- * \param priVar The value of the primary variable
- */
- static Scalar numericEpsilon(const Scalar priVar)
- {
- // define the base epsilon as the geometric mean of 1 and the
- // resolution of the scalar type. E.g. for standard 64 bit
- // floating point values, the resolution is about 10^-16 and
- // the base epsilon is thus approximately 10^-8.
- /*
- static const Scalar baseEps
- = Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(), 1.0);
- */
- static const Scalar baseEps = 1e-10;
- assert(std::numeric_limits<Scalar>::epsilon()*1e4 < baseEps);
- // the epsilon value used for the numeric differentiation is
- // now scaled by the absolute value of the primary variable...
- return baseEps*(std::abs(priVar) + 1.0);
- }
-
-private:
-
- /*!
- * \brief Computes the residual
- *
- * \return The element residual at the current solution.
- */
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler,
- const Element& element, const SolutionVector& curSol)
- {
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- if (isGhost) return NumEqVector(0.0);
-
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
-
- // compatibility with box method
- ElementBoundaryTypes elemBcTypes;
-
- // the actual element's previous time step residual
- auto residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- auto storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- residual += storageResidual;
-
- return residual;
- }
+ using ParentType::ParentType;
/*!
* \brief Computes the derivatives with respect to the given element and adds them
@@ -756,65 +311,13 @@ private:
*
* \return The element residual at the current solution.
*/
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
- const Element& element, const SolutionVector& curSol)
+ LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrix& A, GridVariables& gridVariables)
{
- // get some references for convenience
- const auto& problem = assembler.problem();
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
+ if (this->assembler().isStationaryProblem())
DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
-
- // the global dof of the actual element
- const auto globalI = fvGridGeometry.elementMapper().index(element);
-
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
-
- // the actual element's previous time step residual
- NumEqVector residual(0.0), storageResidual(0.0);
- if (!isGhost)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- residual += storageResidual;
- }
+ // assemble the undeflected residual
+ const auto residual = this->evalLocalResidual()[0];
//////////////////////////////////////////////////////////////////////////////////////////////////
// Calculate derivatives of all dofs in stencil with respect to the dofs in the element. In the //
@@ -822,132 +325,66 @@ private:
// derivatives are non-zero. //
//////////////////////////////////////////////////////////////////////////////////////////////////
- static const std::string group = GET_PROP_VALUE(TypeTag, ModelParameterGroup);
- static const int numericDifferenceMethod = getParamFromGroup<int>(group, "Implicit.NumericDifferenceMethod");
+ // get some aliases for convenience
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ const auto& fvGridGeometry = this->assembler().fvGridGeometry();
+ auto&& curElemVolVars = this->curElemVolVars();
// reference to the element's scv (needed later) and corresponding vol vars
+ const auto globalI = fvGridGeometry.elementMapper().index(element);
const auto& scv = fvGeometry.scv(globalI);
- auto& curVolVars = getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
+ auto& curVolVars = ParentType::getVolVarAccess(gridVariables.curGridVolVars(), curElemVolVars, scv);
// save a copy of the original privars and vol vars in order
// to restore the original solution after deflection
+ const auto& curSol = this->curSol();
const auto origPriVars = curSol[globalI];
const auto origVolVars = curVolVars;
// element solution container to be deflected
- ElementSolutionVector elemSol({origPriVars});
+ ElementSolutionVector elemSol(origPriVars);
+ LocalResidualValues partialDeriv;
// derivatives in the neighbors with repect to the current elements
- for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ for (int pvIdx = 0; pvIdx < numEq; ++pvIdx)
{
// reset derivatives of element dof with respect to itself
- // as well as neighbor derivatives
- NumEqVector partialDeriv(0.0);
-
- if (isGhost)
- partialDeriv[pvIdx] = 1.0;
-
- Scalar eps = numericEpsilon(curVolVars.priVar(pvIdx));
- Scalar delta = 0;
-
- if (numericDifferenceMethod >= 0)
- {
- // we are not using backward differences, i.e. we need to
- // calculate f(x + \epsilon)
-
- // deflect primary variables
- elemSol[0][pvIdx] += eps;
- delta += eps;
-
- // update the volume variables and the flux var cache
- curVolVars.update(elemSol, problem, element, scv);
-
- // calculate the residual with the deflected primary variables
- if (!isGhost)
- {
- partialDeriv = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
- else
- {
- // we are using backward differences, i.e. we don't need
- // to calculate f(x + \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- if (!isGhost)
- partialDeriv = storageResidual;
- }
+ partialDeriv = 0.0;
- if (numericDifferenceMethod <= 0)
+ auto evalStorage = [&](Scalar priVar)
{
- // we are not using forward differences, i.e. we
- // need to calculate f(x - \epsilon)
+ // update the volume variables and calculate
+ // the residual with the deflected primary variables
+ elemSol[0][pvIdx] = priVar;
+ curVolVars.update(elemSol, this->problem(), element, scv);
+ return this->evalLocalStorageResidual()[0];
+ };
+
+ // for non-ghosts compute the derivative numerically
+ if (!this->elementIsGhost())
+ NumericDifferentiation::partialDerivative(evalStorage, elemSol[0][pvIdx], partialDeriv, residual);
+
+ // for ghost elements we assemble a 1.0 where the primary variable and zero everywhere else
+ // as we always solve for a delta of the solution with repect to the initial solution this
+ // results in a delta of zero for ghosts
+ else partialDeriv[pvIdx] = 1.0;
- // deflect the primary variables
- elemSol[0][pvIdx] -= delta + eps;
- delta += eps;
-
- // update the volume variables and the flux var cache
- curVolVars.update(elemSol, problem, element, scv);
-
- // calculate the residual with the deflected primary variables and subtract it
- if (!isGhost)
- {
- partialDeriv -= localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
- else
- {
- // we are using forward differences, i.e. we don't need to
- // calculate f(x - \epsilon) and we can recycle the
- // (already calculated) residual f(x)
- if (!isGhost)
- partialDeriv -= storageResidual;
- }
-
- // divide difference in residuals by the magnitude of the
- // deflections between the two function evaluation
- if (!isGhost)
- partialDeriv /= delta;
+ // add the current partial derivatives to the global jacobian matrix
+ // only diagonal entries for explicit jacobians
+ for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
+ A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
// restore the original state of the scv's volume variables
curVolVars = origVolVars;
// restore the current element solution
elemSol[0][pvIdx] = origPriVars[pvIdx];
-
- // add the current partial derivatives to the global jacobian matrix
- for (int eqIdx = 0; eqIdx < numEq; eqIdx++)
- {
- // the diagonal entries
- A[globalI][globalI][eqIdx][pvIdx] += partialDeriv[eqIdx];
- }
}
// return the original residual
return residual;
}
-private:
- template<class T = TypeTag>
- static typename std::enable_if<!GET_PROP_VALUE(T, EnableGridVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return elemVolVars[scv]; }
-
- template<class T = TypeTag>
- static typename std::enable_if<GET_PROP_VALUE(T, EnableGridVolumeVariablesCache), VolumeVariables&>::type
- getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
- { return gridVolVars.volVars(scv); }
};
/*!
@@ -955,123 +392,19 @@ private:
* \ingroup CCDiscretization
* \brief Cell-centered scheme local assembler using analytic (hand-coded) differentiation and implicit time discretization
*/
-template<class TypeTag>
-class CCLocalAssembler<TypeTag,
- DiffMethod::analytic,
- /*implicit=*/true>
+template<class TypeTag, class Assembler>
+class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/true>
+: public CCLocalAssemblerBase<TypeTag, Assembler,
+ CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, true>, true>
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
- using IndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, true>;
+ using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, true>;
+ using LocalResidualValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
- using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
-
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
public:
-
- /*!
- * \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 Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, jac, element, curSol);
- }
-
- /*!
- * \brief Computes the derivatives with respect to the given element and adds them
- * to the global matrix.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, jac, element, curSol);
- }
-
- /*!
- * \brief Assemble the residual only
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, element, curSol);
- }
-
-private:
-
- /*!
- * \brief Computes the residual
- *
- * \return The element residual at the current solution.
- */
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler,
- const Element& element, const SolutionVector& curSol)
- {
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- if (isGhost) return NumEqVector(0.0);
-
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, curSol);
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
-
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- NumEqVector residual(0.0);
- if (localResidual.isStationary())
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
-
- return residual;
- }
+ using ParentType::ParentType;
/*!
* \brief Computes the derivatives with respect to the given element and adds them
@@ -1079,137 +412,55 @@ private:
*
* \return The element residual at the current solution.
*/
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
- const Element& element, const SolutionVector& curSol)
+ LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
{
- // get some references for convenience
- const auto& problem = assembler.problem();
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bind(element, fvGeometry, curSol);
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
-
- // the global dof of the actual element
- const auto globalI = fvGridGeometry.elementMapper().index(element);
-
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- // is this a stationary simulation?
- const bool isStationary = localResidual.isStationary();
+ // assemble the undeflected residual
+ const auto residual = this->evalLocalResidual()[0];
- // the actual element's current residual (will be returned by this function)
- NumEqVector residual(0.0);
- if (!isGhost)
- {
- if (isStationary)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- else
- {
- prevElemVolVars.bindElement(element, fvGeometry, localResidual.prevSol());
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
- }
- }
+ // get some aliases for convenience
+ const auto& problem = this->problem();
+ const auto& element = this->element();
+ const auto& fvGeometry = this->fvGeometry();
+ auto&& curElemVolVars = this->curElemVolVars();
+ auto&& elemFluxVarsCache = this->elemFluxVarsCache();
// get reference to the element's current vol vars
+ const auto globalI = this->assembler().fvGridGeometry().elementMapper().index(element);
const auto& scv = fvGeometry.scv(globalI);
const auto& volVars = curElemVolVars[scv];
// if the problem is instationary, add derivative of storage term
- if (!isStationary)
- localResidual.addStorageDerivatives(A[globalI][globalI],
- problem,
- element,
- fvGeometry,
- volVars,
- scv);
+ if (!this->assembler().isStationaryProblem())
+ this->localResidual().addStorageDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv);
// add source term derivatives
- localResidual.addSourceDerivatives(A[globalI][globalI],
- problem,
- element,
- fvGeometry,
- volVars,
- scv);
+ this->localResidual().addSourceDerivatives(A[globalI][globalI], problem, element, fvGeometry, volVars, scv);
// add flux derivatives for each scvf
for (const auto& scvf : scvfs(fvGeometry))
{
+ // inner faces
if (!scvf.boundary())
- {
- localResidual.addFluxDerivatives(A[globalI],
- problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- scvf);
- }
+ this->localResidual().addFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
+
+ // boundary faces
else
{
const auto& bcTypes = problem.boundaryTypes(element, scvf);
// add Dirichlet boundary flux derivatives
if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
- {
- localResidual.addCCDirichletFluxDerivatives(A[globalI],
- problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- scvf);
- }
+ this->localResidual().addCCDirichletFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
+
// add Robin ("solution dependent Neumann") boundary flux derivatives
else if (bcTypes.hasNeumann() && !bcTypes.hasDirichlet())
- {
- localResidual.addRobinFluxDerivatives(A[globalI],
- problem,
- element,
- fvGeometry,
- curElemVolVars,
- elemFluxVarsCache,
- scvf);
- }
+ this->localResidual().addRobinFluxDerivatives(A[globalI], problem, element, fvGeometry, curElemVolVars, elemFluxVarsCache, scvf);
+
else
DUNE_THROW(Dune::NotImplemented, "Mixed boundary conditions. Use pure boundary conditions by converting Dirichlet BCs to Robin BCs");
}
}
- // TODO Do we really need this??????????
- // this->model_().updatePVWeights(fvGeometry);
-
- // TODO: Additional dof dependencies???
-
// return element residual
return residual;
}
@@ -1220,202 +471,38 @@ private:
* \ingroup CCDiscretization
* \brief Cell-centered scheme local assembler using analytic (hand-coded) differentiation and explicit time discretization
*/
-template<class TypeTag>
-class CCLocalAssembler<TypeTag,
- DiffMethod::analytic,
- /*implicit=*/false>
+template<class TypeTag, class Assembler>
+class CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, /*implicit=*/false>
+: public CCLocalAssemblerBase<TypeTag, Assembler,
+ CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, false>, false>
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using NumEqVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
- using IndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using ThisType = CCLocalAssembler<TypeTag, Assembler, DiffMethod::analytic, false>;
+ using ParentType = CCLocalAssemblerBase<TypeTag, Assembler, ThisType, false>;
+ using LocalResidualValues = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
- using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
-
- enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
public:
-
- /*!
- * \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 Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, jac, element, curSol);
- }
+ using ParentType::ParentType;
/*!
* \brief Computes the derivatives with respect to the given element and adds them
* to the global matrix.
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, JacobianMatrix& jac,
- const Element& element, const SolutionVector& curSol)
- {
- assemble_(assembler, jac, element, curSol);
- }
-
- /*!
- * \brief Assemble the residual only
- */
- template<class Assembler>
- static void assemble(Assembler& assembler, SolutionVector& res,
- const Element& element, const SolutionVector& curSol)
- {
- const auto globalI = assembler.fvGridGeometry().elementMapper().index(element);
- res[globalI] = assemble_(assembler, element, curSol);
- }
-
-private:
-
- /*!
- * \brief Computes the residual
*
* \return The element residual at the current solution.
*/
- template<class Assembler>
- static NumEqVector assemble_(Assembler& assembler,
- const Element& element, const SolutionVector& curSol)
+ LocalResidualValues assembleJacobianAndResidualImpl(JacobianMatrix& A, const GridVariables& gridVariables)
{
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
- if (isGhost) return NumEqVector(0.0);
-
- // get some references for convenience
- const auto& problem = assembler.problem();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
-
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // the actual element's previous time step residual
- auto residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- auto storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- residual += storageResidual;
-
- return residual;
- }
-
- /*!
- * \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 Assembler>
- static NumEqVector assemble_(Assembler& assembler, JacobianMatrix& A,
- const Element& element, const SolutionVector& curSol)
- {
- // get some references for convenience
- const auto& problem = assembler.problem();
- const auto& fvGridGeometry = assembler.fvGridGeometry();
- auto& localResidual = assembler.localResidual();
- auto& gridVariables = assembler.gridVariables();
-
- // using an explicit assembler doesn't make sense for stationary problems
- if (localResidual.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
-
- // prepare the local views
- auto fvGeometry = localView(assembler.fvGridGeometry());
- fvGeometry.bind(element);
-
- auto curElemVolVars = localView(gridVariables.curGridVolVars());
- curElemVolVars.bindElement(element, fvGeometry, curSol);
-
- auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
- prevElemVolVars.bind(element, fvGeometry, localResidual.prevSol());
-
- auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
- elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
-
- // the global dof of the actual element
- const auto globalI = fvGridGeometry.elementMapper().index(element);
-
- // check for boundaries on the element
- // TODO Do we need them for cell-centered models?
- ElementBoundaryTypes elemBcTypes;
- elemBcTypes.update(problem, element, fvGeometry);
-
- // is the actual element a ghost element?
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
-
- // the actual element's previous time step residual
- NumEqVector residual(0.0), storageResidual(0.0);
- if (!isGhost)
- {
- residual = localResidual.eval(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- storageResidual = localResidual.evalStorage(problem,
- element,
- fvGeometry,
- prevElemVolVars,
- curElemVolVars,
- elemBcTypes,
- elemFluxVarsCache)[0];
-
- residual += storageResidual;
- }
+ // assemble the undeflected residual
+ const auto residual = this->evalLocalResidual()[0];
// get reference to the element's current vol vars
- const auto& scv = fvGeometry.scv(globalI);
- const auto& volVars = curElemVolVars[scv];
+ const auto globalI = this->assembler().fvGridGeometry().elementMapper().index(this->element());
+ const auto& scv = this->fvGeometry().scv(globalI);
+ const auto& volVars = this->curElemVolVars()[scv];
- // add derivative of storage term
- localResidual.addStorageDerivatives(A[globalI][globalI],
- problem,
- element,
- fvGeometry,
- volVars,
- scv);
+ // add hand-code derivative of storage term
+ this->localResidual().addStorageDerivatives(A[globalI][globalI], this->problem(), this->element(), this->fvGeometry(), volVars, scv);
// return the original residual
return residual;
diff --git a/dumux/assembly/cclocalresidual.hh b/dumux/assembly/cclocalresidual.hh
index bc13dd2a8961fe643a804fab43c520c7b0f22aa7..c90d2d47765ea392936d19076fab772adf8488ad 100644
--- a/dumux/assembly/cclocalresidual.hh
+++ b/dumux/assembly/cclocalresidual.hh
@@ -25,8 +25,7 @@
#ifndef DUMUX_CC_LOCAL_RESIDUAL_HH
#define DUMUX_CC_LOCAL_RESIDUAL_HH
-#include <dune/istl/matrix.hh>
-
+#include <dumux/common/reservedblockvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/assembly/fvlocalresidual.hh>
@@ -43,7 +42,6 @@ class CCLocalResidual : public FVLocalResidual<TypeTag>
using ParentType = FVLocalResidual<TypeTag>;
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
- using ElementResidualVector = Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, NumEqVector)>;
using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
@@ -52,6 +50,7 @@ class CCLocalResidual : public FVLocalResidual<TypeTag>
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
public:
+ using ElementResidualVector = typename ParentType::ElementResidualVector;
using ParentType::ParentType;
//! evaluate the flux residual for a sub control volume face and add to residual
diff --git a/dumux/assembly/fvassembler.hh b/dumux/assembly/fvassembler.hh
index d8142264a3c4bbc22864c6633c7064cf14340a34..4310071b2a7d97005dc21e045ac3cd38a330c8a6 100644
--- a/dumux/assembly/fvassembler.hh
+++ b/dumux/assembly/fvassembler.hh
@@ -41,10 +41,10 @@ namespace Dumux {
/*!
* \ingroup Assembly
- * \brief A linear system assembler (residual and Jacobian) for finite volume schemes
- * \tparam TypeTag the TypeTag
- * \tparam diffMethod the differentiation method to residual compute derivatives
- * \tparam isImplicit if to use an implicit or explicit time discretization
+ * \brief A linear system assembler (residual and Jacobian) for finite volume schemes (box, tpfa, mpfa, ...)
+ * \tparam TypeTag The TypeTag
+ * \tparam diffMethod The differentiation method to residual compute derivatives
+ * \tparam isImplicit Specifies whether the time discretization is implicit or not not (i.e. explicit)
*/
template<class TypeTag, DiffMethod diffMethod, bool isImplicit = true>
class FVAssembler
@@ -52,44 +52,55 @@ class FVAssembler
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
using VertexMapper = typename GET_PROP_TYPE(TypeTag, VertexMapper);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+ using Element = typename GridView::template Codim<0>::Entity;
using TimeLoop = TimeLoopBase<Scalar>;
- static constexpr int dim = GridView::dimension;
static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
- using LocalAssembler = std::conditional_t<isBox, BoxLocalAssembler<TypeTag, diffMethod, isImplicit>,
- CCLocalAssembler<TypeTag, diffMethod, isImplicit>>;
+
+ using ThisType = FVAssembler<TypeTag, diffMethod, isImplicit>;
+ using LocalAssembler = std::conditional_t<isBox, BoxLocalAssembler<TypeTag, ThisType, diffMethod, isImplicit>,
+ CCLocalAssembler<TypeTag, ThisType, diffMethod, isImplicit>>;
public:
using ResidualType = SolutionVector;
- //! The constructor for stationary problems
+ /*!
+ * \brief The constructor for stationary problems
+ * \note the grid variables might be temporarily changed during assembly (if caching is enabled)
+ * it is however guaranteed that the state after assembly will be the same as before
+ */
FVAssembler(std::shared_ptr<const Problem> problem,
std::shared_ptr<const FVGridGeometry> fvGridGeometry,
std::shared_ptr<GridVariables> gridVariables)
: problem_(problem)
, fvGridGeometry_(fvGridGeometry)
, gridVariables_(gridVariables)
- , stationary_(true)
+ , timeLoop_()
+ , isStationaryProblem_(true)
{
static_assert(isImplicit, "Explicit assembler for stationary problem doesn't make sense!");
}
- //! The constructor for instationary problems
+ /*!
+ * \brief The constructor for instationary problems
+ * \note the grid variables might be temporarily changed during assembly (if caching is enabled)
+ * it is however guaranteed that the state after assembly will be the same as before
+ */
FVAssembler(std::shared_ptr<const Problem> problem,
std::shared_ptr<const FVGridGeometry> fvGridGeometry,
std::shared_ptr<GridVariables> gridVariables,
- std::shared_ptr<TimeLoop> timeLoop)
+ std::shared_ptr<const TimeLoop> timeLoop)
: problem_(problem)
, fvGridGeometry_(fvGridGeometry)
, gridVariables_(gridVariables)
- , localResidual_(timeLoop)
- , stationary_(false)
+ , timeLoop_(timeLoop)
+ , isStationaryProblem_(false)
{}
/*!
@@ -98,50 +109,15 @@ public:
*/
void assembleJacobianAndResidual(const SolutionVector& curSol)
{
- if (!stationary_ && localResidual_.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
-
- if(!jacobian_)
- {
- jacobian_ = std::make_shared<JacobianMatrix>();
- jacobian_->setBuildMode(JacobianMatrix::random);
- setJacobianPattern();
- }
-
- if(!residual_)
- {
- residual_ = std::make_shared<SolutionVector>();
- setResidualSize();
- }
-
+ checkAssemblerState_();
resetJacobian_();
resetResidual_();
- bool succeeded;
- // try assembling the global linear system
- try
+ assemble_([&](const Element& element)
{
- // let the local assembler add the element contributions
- for (const auto element : elements(gridView()))
- LocalAssembler::assemble(*this, *jacobian_, *residual_, element, curSol);
-
- // if we get here, everything worked well
- succeeded = true;
- if (gridView().comm().size() > 1)
- succeeded = gridView().comm().min(succeeded);
- }
- // throw exception if a problem ocurred
- catch (NumericalProblem &e)
- {
- std::cout << "rank " << gridView().comm().rank()
- << " caught an exception while assembling:" << e.what()
- << "\n";
- succeeded = false;
- if (gridView().comm().size() > 1)
- succeeded = gridView().comm().min(succeeded);
- }
- if (!succeeded)
- DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
+ LocalAssembler localAssembler(*this, element, curSol);
+ localAssembler.assembleJacobianAndResidual(*jacobian_, *residual_, *gridVariables_);
+ });
}
/*!
@@ -149,72 +125,39 @@ public:
*/
void assembleJacobian(const SolutionVector& curSol)
{
- if (!stationary_ && localResidual_.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
-
- if(!jacobian_)
- {
- jacobian_ = std::make_shared<JacobianMatrix>();
- jacobian_->setBuildMode(JacobianMatrix::random);
- setJacobianPattern();
- }
-
+ checkAssemblerState_();
resetJacobian_();
- bool succeeded;
- // try assembling the global linear system
- try
- {
- // let the local assembler add the element contributions
- for (const auto element : elements(gridView()))
- LocalAssembler::assemble(*this, *jacobian_, element, curSol);
-
- // if we get here, everything worked well
- succeeded = true;
- if (gridView().comm().size() > 1)
- succeeded = gridView().comm().min(succeeded);
- }
- // throw exception if a problem ocurred
- catch (NumericalProblem &e)
+ assemble_([&](const Element& element)
{
- std::cout << "rank " << gridView().comm().rank()
- << " caught an exception while assembling:" << e.what()
- << "\n";
- succeeded = false;
- if (gridView().comm().size() > 1)
- succeeded = gridView().comm().min(succeeded);
- }
- if (!succeeded)
- DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
+ LocalAssembler localAssembler(*this, element, curSol);
+ localAssembler.assembleJacobianAndResidual(*jacobian_, *gridVariables_);
+ });
}
- //! compute the residuals
+ //! compute the residuals using the internal residual
void assembleResidual(const SolutionVector& curSol)
{
- if(!residual_)
- {
- residual_ = std::make_shared<SolutionVector>();
- setResidualSize();
- }
-
+ resetResidual_();
assembleResidual(*residual_, curSol);
}
- //! compute the residuals
+ //! assemble a residual r
void assembleResidual(ResidualType& r, const SolutionVector& curSol) const
{
- if (!stationary_ && localResidual_.isStationary())
- DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
+ checkAssemblerState_();
// update the grid variables for the case of active caching
gridVariables_->update(curSol);
- // let the local assembler add the element contributions
- for (const auto element : elements(gridView()))
- LocalAssembler::assemble(*this, r, element, curSol);
+ assemble_([&](const Element& element)
+ {
+ LocalAssembler localAssembler(*this, element, curSol);
+ localAssembler.assembleResidual(r);
+ });
}
- //! computes the residual norm
+ //! compute the residual and return it's vector norm
Scalar residualNorm(const SolutionVector& curSol) const
{
ResidualType residual(numDofs());
@@ -262,170 +205,256 @@ public:
/*!
* \brief The version without arguments uses the default constructor to create
- * the jacobian and residual objects in this assembler.
+ * the jacobian and residual objects in this assembler if you don't need them outside this class
*/
void setLinearSystem()
{
jacobian_ = std::make_shared<JacobianMatrix>();
jacobian_->setBuildMode(JacobianMatrix::random);
-
residual_ = std::make_shared<SolutionVector>();
setJacobianPattern();
setResidualSize();
}
- /*!
- * \brief Sets the solution from which to start the time integration. Has to be
- * called prior to assembly for time-dependent problems.
- */
- void setPreviousSolution(const SolutionVector& u)
- { localResidual_.setPreviousSolution(u); }
-
- /*!
- * \brief Return the solution that has been set as the previous one.
- */
- const SolutionVector& prevSol() const
- { return localResidual_.prevSol(); }
-
/*!
* \brief Resizes the jacobian and sets the jacobian' sparsity pattern.
*/
-
void setJacobianPattern()
{
// resize the jacobian and the residual
const auto numDofs = this->numDofs();
jacobian_->setSize(numDofs, numDofs);
- // get occupation pattern of the jacobian
+ // create occupation pattern of the jacobian
Dune::MatrixIndexSet occupationPattern;
occupationPattern.resize(numDofs, numDofs);
- // matrix pattern for implicit jacobians
- if (isImplicit)
- setImplicitJacobianPattern_(occupationPattern, numDofs);
-
- // matrix pattern for explicit jacobians -> diagonal matrix
- else
- for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
- occupationPattern.add(globalI, globalI);
+ // set the jacobian pattern depending on space and time discretization
+ setJacobianPattern_(occupationPattern, numDofs);
// export pattern to jacobian
occupationPattern.exportIdx(*jacobian_);
}
- /*!
- * \brief Resizes the residual
- */
+ //! Resizes the residual
void setResidualSize()
{ residual_->resize(numDofs()); }
- //! cell-centered schemes have one dof per cell
+ //! Returns the number of degrees of freedom
std::size_t numDofs() const
{ return fvGridGeometry_->numDofs(); }
+ //! The problem
const Problem& problem() const
{ return *problem_; }
+ //! The global finite volume geometry
const FVGridGeometry& fvGridGeometry() const
{ return *fvGridGeometry_; }
+ //! The gridview
const GridView& gridView() const
{ return fvGridGeometry().gridView(); }
+ //! The global grid variables
GridVariables& gridVariables()
{ return *gridVariables_; }
+ //! The global grid variables
const GridVariables& gridVariables() const
{ return *gridVariables_; }
+ //! The jacobian matrix
JacobianMatrix& jacobian()
- {
- if (!residual_)
- DUNE_THROW(Dune::InvalidStateException, "No jacobian was set.");
- return *jacobian_;
- }
+ { return *jacobian_; }
+ //! The residual vector (rhs)
SolutionVector& residual()
- {
- if (!residual_)
- DUNE_THROW(Dune::InvalidStateException, "No residual was set.");
- return *residual_;
- }
+ { return *residual_; }
+
+ //! The solution of the previous time step
+ const SolutionVector& prevSol() const
+ { return *prevSol_; }
+
+ /*!
+ * \brief Set time loop for instationary problems
+ * \note calling this turns this into a stationary assembler
+ */
+ void setTimeManager(std::shared_ptr<const TimeLoop> timeLoop)
+ { timeLoop_ = timeLoop_; isStationaryProblem_ = true; }
- const LocalResidual& localResidual() const
- { return localResidual_; }
+ /*!
+ * \brief Sets the solution from which to start the time integration. Has to be
+ * called prior to assembly for time-dependent problems.
+ */
+ void setPreviousSolution(const SolutionVector& u)
+ { prevSol_ = &u; }
+
+ /*!
+ * \brief Whether we are assembling a stationary or instationary problem
+ */
+ bool isStationaryProblem() const
+ { return isStationaryProblem_; }
+
+ /*!
+ * \brief Create a local residual object (used by the local assembler)
+ */
+ LocalResidual localResidual() const
+ { return LocalResidual(problem_.get(), timeLoop_.get()); }
private:
- // reset the residual to 0.0
+ // reset the residual vector to 0.0
void resetResidual_()
{
+ if(!residual_)
+ {
+ residual_ = std::make_shared<SolutionVector>();
+ setResidualSize();
+ }
+
(*residual_) = 0.0;
}
- // reset the jacobian to 0.0
+ // reset the jacobian vector to 0.0
void resetJacobian_()
{
+ if(!jacobian_)
+ {
+ jacobian_ = std::make_shared<JacobianMatrix>();
+ jacobian_->setBuildMode(JacobianMatrix::random);
+ setJacobianPattern();
+ }
+
(*jacobian_) = 0.0;
}
+ // check if the assembler is in a correct state for assembly
+ void checkAssemblerState_() const
+ {
+ if (!isStationaryProblem_ && !prevSol_)
+ DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
+ }
+
//! Implicit jacobian pattern for cell-centered fv schemes
- template<typename T = TypeTag>
- std::enable_if_t<GET_PROP_VALUE(T, DiscretizationMethod) != DiscretizationMethods::Box, void>
- setImplicitJacobianPattern_(Dune::MatrixIndexSet& pattern, std::size_t numDofs)
+ template<class T = TypeTag, typename std::enable_if_t<GET_PROP_VALUE(T, DiscretizationMethod) != DiscretizationMethods::Box, int> = 0>
+ void setJacobianPattern_(Dune::MatrixIndexSet& pattern, std::size_t numDofs)
{
- for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
+ // matrix pattern for implicit Jacobians
+ if (isImplicit)
{
- pattern.add(globalI, globalI);
- for (const auto& dataJ : fvGridGeometry().connectivityMap()[globalI])
- pattern.add(dataJ.globalJ, globalI);
-
- // reserve index for additional user defined DOF dependencies
- // const auto& additionalDofDependencies = problem().getAdditionalDofDependencies(globalI);
- // for (auto globalJ : additionalDofDependencies)
- // pattern.add(globalI, globalJ);
+ for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
+ {
+ pattern.add(globalI, globalI);
+ for (const auto& dataJ : fvGridGeometry().connectivityMap()[globalI])
+ pattern.add(dataJ.globalJ, globalI);
+
+ // reserve index for additional user defined DOF dependencies
+ // const auto& additionalDofDependencies = problem().getAdditionalDofDependencies(globalI);
+ // for (auto globalJ : additionalDofDependencies)
+ // pattern.add(globalI, globalJ);
+ }
+ }
+
+ // matrix pattern for explicit Jacobians -> diagonal matrix
+ else
+ {
+ for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
+ pattern.add(globalI, globalI);
}
}
//! Implicit jacobian pattern for vertex-centered fv schemes
- template<typename T = TypeTag>
- std::enable_if_t<GET_PROP_VALUE(T, DiscretizationMethod) == DiscretizationMethods::Box, void>
- setImplicitJacobianPattern_(Dune::MatrixIndexSet& pattern, std::size_t numDofs)
+ template<class T = TypeTag, typename std::enable_if_t<GET_PROP_VALUE(T, DiscretizationMethod) == DiscretizationMethods::Box, int> = 0>
+ void setJacobianPattern_(Dune::MatrixIndexSet& pattern, std::size_t numDofs)
{
- for (const auto& element : elements(fvGridGeometry().gridView()))
+ // matrix pattern for implicit Jacobians
+ if (isImplicit)
{
- for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
+ static constexpr int dim = GridView::dimension;
+ for (const auto& element : elements(fvGridGeometry().gridView()))
{
- const auto globalI = fvGridGeometry().vertexMapper().subIndex(element, vIdx, dim);
- for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(dim); ++vIdx2)
+ for (unsigned int vIdx = 0; vIdx < element.subEntities(dim); ++vIdx)
{
- const auto globalJ = fvGridGeometry().vertexMapper().subIndex(element, vIdx2, dim);
- pattern.add(globalI, globalJ);
- pattern.add(globalJ, globalI);
+ const auto globalI = fvGridGeometry().vertexMapper().subIndex(element, vIdx, dim);
+ for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(dim); ++vIdx2)
+ {
+ const auto globalJ = fvGridGeometry().vertexMapper().subIndex(element, vIdx2, dim);
+ pattern.add(globalI, globalJ);
+ pattern.add(globalJ, globalI);
+ }
}
}
}
+
+ // matrix pattern for explicit Jacobians -> diagonal matrix
+ else
+ {
+ for (unsigned int globalI = 0; globalI < numDofs; ++globalI)
+ pattern.add(globalI, globalI);
+ }
+ }
+
+ /*!
+ * \brief A method assembling something per element
+ * \note Handles exceptions for parallel runs
+ * \throws NumericalProblem on all processes if something throwed during assembly
+ */
+ template<typename AssembleElementFunc>
+ void assemble_(AssembleElementFunc&& assembleElement) const
+ {
+ // a state that will be checked on all processes
+ bool succeeded = false;
+
+ // try assembling using the local assembly function
+ try
+ {
+ // let the local assembler add the element contributions
+ for (const auto& element : elements(gridView()))
+ assembleElement(element);
+
+ // if we get here, everything worked well on this process
+ succeeded = true;
+ }
+ // throw exception if a problem ocurred
+ catch (NumericalProblem &e)
+ {
+ std::cout << "rank " << gridView().comm().rank()
+ << " caught an exception while assembling:" << e.what()
+ << "\n";
+ succeeded = false;
+ }
+
+ // make sure everything worked well on all processes
+ if (gridView().comm().size() > 1)
+ succeeded = gridView().comm().min(succeeded);
+
+ // if not succeeded rethrow the error on all processes
+ if (!succeeded)
+ DUNE_THROW(NumericalProblem, "A process did not succeed in linearizing the system");
}
- // pointer to the problem to be solved
+ //! pointer to the problem to be solved
std::shared_ptr<const Problem> problem_;
- // the finite volume geometry of the grid
+ //! the finite volume geometry of the grid
std::shared_ptr<const FVGridGeometry> fvGridGeometry_;
- // the variables container for the grid
+ //! the variables container for the grid
std::shared_ptr<GridVariables> gridVariables_;
- // shared pointers to the jacobian matrix and residual
- std::shared_ptr<JacobianMatrix> jacobian_;
- std::shared_ptr<SolutionVector> residual_;
+ //! the time loop for instationary problem assembly
+ std::shared_ptr<const TimeLoop> timeLoop_;
+
+ //! an observing pointer to the previous solution for instationary problems
+ const SolutionVector* prevSol_ = nullptr;
- // class computing the residual of an element
- LocalResidual localResidual_;
+ //! if this assembler is assembling an instationary problem
+ bool isStationaryProblem_;
- // if this assembler is assembling a time dependent problem
- bool stationary_;
+ //! shared pointers to the jacobian matrix and residual
+ std::shared_ptr<JacobianMatrix> jacobian_;
+ std::shared_ptr<SolutionVector> residual_;
};
} // namespace Dumux
diff --git a/dumux/assembly/fvlocalassemblerbase.hh b/dumux/assembly/fvlocalassemblerbase.hh
new file mode 100644
index 0000000000000000000000000000000000000000..dbae670f0b5ca590c0bb6867d820e56afccde073
--- /dev/null
+++ b/dumux/assembly/fvlocalassemblerbase.hh
@@ -0,0 +1,290 @@
+// -*- 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 2 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
+ * \copydoc Dumux::FVLocalAssemblerBase
+ */
+#ifndef DUMUX_FV_LOCAL_ASSEMBLER_BASE_HH
+#define DUMUX_FV_LOCAL_ASSEMBLER_BASE_HH
+
+#include <dune/common/reservedvector.hh>
+#include <dune/grid/common/gridenums.hh> // for GhostEntity
+#include <dune/istl/matrixindexset.hh>
+
+#include <dumux/common/reservedblockvector.hh>
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/assembly/diffmethod.hh>
+#include <dumux/assembly/numericdifferentiation.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Assembly
+ * \brief A base class for all local assemblers
+ * \tparam TypeTag The TypeTag
+ * \tparam Assembler The assembler type
+ * \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 FVLocalAssemblerBase
+{
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using LocalResidual = typename GET_PROP_TYPE(TypeTag, LocalResidual);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using JacobianMatrix = typename GET_PROP_TYPE(TypeTag, JacobianMatrix);
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using GridVolumeVariables = typename GET_PROP_TYPE(TypeTag, GridVolumeVariables);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
+ using Element = typename GridView::template Codim<0>::Entity;
+
+public:
+
+ /*!
+ * \brief The constructor
+ */
+ explicit FVLocalAssemblerBase(const Assembler& assembler,
+ const Element& element,
+ const SolutionVector& curSol)
+ : assembler_(assembler)
+ , element_(element)
+ , curSol_(curSol)
+ , fvGeometry_(localView(assembler.fvGridGeometry()))
+ , curElemVolVars_(localView(assembler.gridVariables().curGridVolVars()))
+ , prevElemVolVars_(localView(assembler.gridVariables().prevGridVolVars()))
+ , elemFluxVarsCache_(localView(assembler.gridVariables().gridFluxVarsCache()))
+ , localResidual_(assembler.localResidual())
+ , elementIsGhost_((element.partitionType() == Dune::GhostEntity))
+ {}
+
+ /*!
+ * \brief Convenience function to evaluate the complete local residual for the current element. Automatically chooses the the appropriate
+ * element volume variables.
+ */
+ auto evalLocalResidual() const
+ {
+ if(this->assembler().isStationaryProblem() && !isImplicit())
+ DUNE_THROW(Dune::InvalidStateException, "Using explicit jacobian assembler with stationary local residual");
+
+ if(elementIsGhost())
+ {
+ using ResdiualType = decltype(evalLocalResidual(std::declval<ElementVolumeVariables>()));
+ return ResdiualType(0.0);
+ }
+
+ return isImplicit() ? evalLocalResidual(curElemVolVars())
+ : evalLocalResidual(prevElemVolVars());
+ }
+
+ /*!
+ * \brief Evaluates the complete local residual for the current element.
+ * \param elemVolVars The element volume variables
+ */
+ auto evalLocalResidual(const ElementVolumeVariables& elemVolVars) const
+ {
+ if (!assembler().isStationaryProblem())
+ {
+ auto residual = evalLocalFluxAndSourceResidual(elemVolVars);
+ residual += evalLocalStorageResidual();
+ return residual;
+ }
+ else
+ return evalLocalFluxAndSourceResidual(elemVolVars);
+ }
+
+ /*!
+ * \brief Convenience function to evaluate the flux and source terms (i.e, the terms without a time derivative)
+ * of the local residual for the current element. Automatically chooses the the appropriate
+ * element volume variables.
+ */
+ auto evalLocalFluxAndSourceResidual() const
+ {
+ return isImplicit() ? evalLocalFluxAndSourceResidual(curElemVolVars())
+ : evalLocalFluxAndSourceResidual(prevElemVolVars());
+ }
+
+ /*!
+ * \brief Evaluates the flux and source terms (i.e, the terms without a time derivative)
+ * of the local residual for the current element.
+ *
+ * \param elemVolVars The element volume variables
+ */
+ auto evalLocalFluxAndSourceResidual(const ElementVolumeVariables& elemVolVars) const
+ {
+ return localResidual_.evalFluxSource(element_, fvGeometry_, elemVolVars, elemFluxVarsCache_, elemBcTypes_);
+ }
+
+ /*!
+ * \brief Convenience function to evaluate storage term (i.e, the term with a time derivative)
+ * of the local residual for the current element. Automatically chooses the the appropriate
+ * element volume variables.
+ */
+ auto evalLocalStorageResidual() const
+ {
+ return localResidual_.evalStorage(element_, fvGeometry_, prevElemVolVars_, curElemVolVars_);
+ }
+
+ /*!
+ * \brief Convenience function bind and prepare all relevant variables required for the
+ * evaluation of the local residual.
+ */
+ void bindLocalViews()
+ {
+ // get some references for convenience
+ const auto& element = this->element();
+ const auto& curSol = this->curSol();
+ const auto& prevSol = this->assembler().prevSol();
+ auto&& fvGeometry = this->fvGeometry();
+ auto&& curElemVolVars = this->curElemVolVars();
+ auto&& prevElemVolVars = this->prevElemVolVars();
+ auto&& elemFluxVarsCache = this->elemFluxVarsCache();
+
+ // bind the caches
+ fvGeometry.bind(element);
+
+ if(isImplicit())
+ {
+ curElemVolVars.bind(element, fvGeometry, curSol);
+ elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
+ if (!this->assembler().isStationaryProblem())
+ prevElemVolVars.bindElement(element, fvGeometry, this->assembler().prevSol());
+ }
+ else
+ {
+ curElemVolVars.bindElement(element, fvGeometry, curSol);
+ prevElemVolVars.bind(element, fvGeometry, prevSol);
+ elemFluxVarsCache.bind(element, fvGeometry, prevElemVolVars);
+ }
+ }
+
+ //! The problem
+ const Problem& problem() const
+ { return assembler_.problem(); }
+
+ //! The assembler
+ const Assembler& assembler() const
+ { return assembler_; }
+
+ //! The current element
+ const Element& element() const
+ { return element_; }
+
+ //! Returns if element is a ghost entity
+ bool elementIsGhost() const
+ { return elementIsGhost_; }
+
+ //! The current solution
+ const SolutionVector& curSol() const
+ { return curSol_; }
+
+ //! The global finite volume geometry
+ FVElementGeometry& fvGeometry()
+ { return fvGeometry_; }
+
+ //! The current element volume variables
+ ElementVolumeVariables& curElemVolVars()
+ { return curElemVolVars_; }
+
+ //! The element volume variables of the provious time step
+ ElementVolumeVariables& prevElemVolVars()
+ { return prevElemVolVars_; }
+
+ //! The element flux variables cache
+ ElementFluxVariablesCache& elemFluxVarsCache()
+ { return elemFluxVarsCache_; }
+
+ //! The local residual for the current element
+ LocalResidual& localResidual()
+ { return localResidual_; }
+
+ //! The element's boundary types
+ ElementBoundaryTypes& elemBcTypes()
+ { return elemBcTypes_; }
+
+ //! The finite volume geometry
+ const FVElementGeometry& fvGeometry() const
+ { return fvGeometry_; }
+
+ //! The current element volume variables
+ const ElementVolumeVariables& curElemVolVars() const
+ { return curElemVolVars_; }
+
+ //! The element volume variables of the provious time step
+ const ElementVolumeVariables& prevElemVolVars() const
+ { return prevElemVolVars_; }
+
+ //! The element flux variables cache
+ const ElementFluxVariablesCache& elemFluxVarsCache() const
+ { return elemFluxVarsCache_; }
+
+ //! The element's boundary types
+ const ElementBoundaryTypes& elemBcTypes() const
+ { return elemBcTypes_; }
+
+ //! The local residual for the current element
+ const LocalResidual& localResidual() const
+ { return localResidual_; }
+
+ constexpr bool isImplicit() const
+ { return implicit; }
+
+protected:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+ template<class T = TypeTag, typename std::enable_if_t<!GET_PROP_VALUE(T, EnableGridVolumeVariablesCache), int> = 0>
+ VolumeVariables& getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return elemVolVars[scv]; }
+
+ template<class T = TypeTag, typename std::enable_if_t<GET_PROP_VALUE(T, EnableGridVolumeVariablesCache), int> = 0>
+ VolumeVariables& getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
+ { return gridVolVars.volVars(scv); }
+
+private:
+
+ const Assembler& assembler_; //!< access pointer to assembler instance
+ const Element& element_; //!< the element whose residual is assembled
+ const SolutionVector& curSol_; //!< the current solution
+
+ FVElementGeometry fvGeometry_;
+ ElementVolumeVariables curElemVolVars_;
+ ElementVolumeVariables prevElemVolVars_;
+ ElementFluxVariablesCache elemFluxVarsCache_;
+ ElementBoundaryTypes elemBcTypes_;
+
+ LocalResidual localResidual_; //!< the local residual evaluating the equations per element
+ bool elementIsGhost_; //!< whether the element's partitionType is ghost
+};
+
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/assembly/fvlocalresidual.hh b/dumux/assembly/fvlocalresidual.hh
index 67b4f2b5350d6eb7d4975d63bc733bd17e9b7421..3e06c94d388f9e506e94027724775678c5d5f0d1 100644
--- a/dumux/assembly/fvlocalresidual.hh
+++ b/dumux/assembly/fvlocalresidual.hh
@@ -24,11 +24,13 @@
#ifndef DUMUX_FV_LOCAL_RESIDUAL_HH
#define DUMUX_FV_LOCAL_RESIDUAL_HH
+#include <dune/common/deprecated.hh>
#include <dune/common/exceptions.hh>
#include <dune/istl/bvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/timeloop.hh>
+#include <dumux/common/reservedblockvector.hh>
#include <dumux/discretization/methods.hh>
namespace Dumux {
@@ -45,13 +47,13 @@ class FVLocalResidual
using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using ElementResidualVector = Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, NumEqVector)>;
using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
@@ -61,13 +63,14 @@ class FVLocalResidual
using TimeLoop = TimeLoopBase<Scalar>;
public:
- //! the constructor for stationary problems
- FVLocalResidual() : prevSol_(nullptr) {}
-
- //! the constructor for instationary problems
- FVLocalResidual(std::shared_ptr<TimeLoop> timeLoop)
- : timeLoop_(timeLoop)
- , prevSol_(nullptr)
+ //! the container storing all element residuals
+ using ElementResidualVector = ReservedBlockVector<ResidualVector, FVElementGeometry::maxNumElementScvs>;
+
+ //! the constructor
+ FVLocalResidual(const Problem* problem,
+ const TimeLoop* timeLoop = nullptr)
+ : problem_(problem)
+ , timeLoop_(timeLoop)
{}
/*!
@@ -141,6 +144,7 @@ public:
* \param bcTypes The types of the boundary conditions for all boundary entities of an element
* \param elemFluxVarsCache The flux variable caches for the element stencil
*/
+ DUNE_DEPRECATED_MSG("eval is deprecated because it doesn't allow to specify on which time level to evaluate. Use evalFluxSource, and evalStorage instead!")
ElementResidualVector eval(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
@@ -150,7 +154,6 @@ public:
const ElementFluxVariablesCache& elemFluxVarsCache) const
{
assert(timeLoop_ && "no time loop set for storage term evaluation");
- assert(prevSol_ && "no solution set for storage term evaluation");
// initialize the residual vector for all scvs in this element
ElementResidualVector residual(fvGeometry.numScv());
@@ -177,7 +180,6 @@ public:
* \brief Compute the storage local residual, i.e. the deviation of the
* storage term from zero for instationary problems.
*
- * \param problem The problem to solve
* \param element The DUNE Codim<0> entity for which the residual
* ought to be calculated
* \param fvGeometry The finite-volume geometry of the element
@@ -185,30 +187,21 @@ public:
* sub-control volumes of the element at the previous time level
* \param curElemVolVars The volume averaged variables for all
* sub-control volumes of the element at the current time level
- * \param bcTypes The types of the boundary conditions for all boundary entities of an element
- * \param elemFluxVarsCache The flux variable caches for the element stencil
*/
- ElementResidualVector evalStorage(const Problem& problem,
- const Element& element,
+ ElementResidualVector evalStorage(const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const ElementBoundaryTypes &bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache) const
+ const ElementVolumeVariables& curElemVolVars) const
{
assert(timeLoop_ && "no time loop set for storage term evaluation");
- assert(prevSol_ && "no solution set for storage term evaluation");
// initialize the residual vector for all scvs in this element
ElementResidualVector residual(fvGeometry.numScv());
- residual = 0.0;
// evaluate the volume terms (storage + source terms)
+ // forward to the local residual specialized for the discretization methods
for (auto&& scv : scvs(fvGeometry))
- {
- //! foward to the local residual specialized for the discretization methods
- asImp().evalStorage(residual, problem, element, fvGeometry, prevElemVolVars, curElemVolVars, scv);
- }
+ asImp().evalStorage(residual, this->problem(), element, fvGeometry, prevElemVolVars, curElemVolVars, scv);
return residual;
}
@@ -226,29 +219,34 @@ public:
* \param bcTypes The types of the boundary conditions for all boundary entities of an element
* \param elemFluxVarsCache The flux variable caches for the element stencil
*/
+ DUNE_DEPRECATED_MSG("Use evalFluxSource instead!")
ElementResidualVector eval(const Problem& problem,
const Element& element,
const FVElementGeometry& fvGeometry,
const ElementVolumeVariables& curElemVolVars,
const ElementBoundaryTypes &bcTypes,
const ElementFluxVariablesCache& elemFluxVarsCache) const
+ {
+ return evalFluxSource(element, fvGeometry, curElemVolVars, elemFluxVarsCache, bcTypes);
+ }
+
+ ElementResidualVector evalFluxSource(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes &bcTypes) const
{
// initialize the residual vector for all scvs in this element
ElementResidualVector residual(fvGeometry.numScv());
- residual = 0.0;
// evaluate the volume terms (storage + source terms)
+ // forward to the local residual specialized for the discretization methods
for (auto&& scv : scvs(fvGeometry))
- {
- //! foward to the local residual specialized for the discretization methods
- asImp().evalSource(residual, problem, element, fvGeometry, curElemVolVars, scv);
- }
+ asImp().evalSource(residual, this->problem(), element, fvGeometry, elemVolVars, scv);
+ // forward to the local residual specialized for the discretization methods
for (auto&& scvf : scvfs(fvGeometry))
- {
- //! foward to the local residual specialized for the discretization methods
- asImp().evalFlux(residual, problem, element, fvGeometry, curElemVolVars, bcTypes, elemFluxVarsCache, scvf);
- }
+ asImp().evalFlux(residual, this->problem(), element, fvGeometry, elemVolVars, bcTypes, elemFluxVarsCache, scvf);
return residual;
}
@@ -552,40 +550,22 @@ public:
//\}
/*!
- * \brief Sets the solution from which to start the time integration.
- * \note Has to be called prior to assembly for time-dependent problems.
- */
- void setPreviousSolution(const SolutionVector& u)
- { prevSol_ = &u; }
-
- /*!
- * \brief Return the solution of the previous time level that has been set
- */
- const SolutionVector& prevSol() const
- {
- assert(prevSol_ && "no solution set for storage term evaluation");
- return *prevSol_;
- }
-
- /*!
- * \brief if the residual is stationary (storage residual is strictly zero)
- * \note If no solution has been set, we treat the problem as stationary.
+ * \name Interfaces accessed by local residual implementations
*/
- bool isStationary() const
- { return !prevSol_; }
+ // \{
- // \}
-protected:
- //! the timeloop for instationary problems
- //! calling this for stationary leads to undefined behaviour
- TimeLoop& timeLoop()
- { return *timeLoop_; }
+ //! the problem
+ const Problem& problem() const
+ { return *problem_; }
//! the timeloop for instationary problems
//! calling this for stationary leads to undefined behaviour
const TimeLoop& timeLoop() const
{ return *timeLoop_; }
+ // \}
+protected:
+
Implementation &asImp()
{ return *static_cast<Implementation*>(this); }
@@ -593,8 +573,8 @@ protected:
{ return *static_cast<const Implementation*>(this); }
private:
- std::shared_ptr<TimeLoop> timeLoop_; //!< the timeloop for instationary problems
- const SolutionVector* prevSol_; //!< the solution of the previous time level for instationary problems
+ const Problem* problem_; //!< the problem we are assembling this residual for
+ const TimeLoop* timeLoop_; //!< the timeloop for instationary problems
};
} // end namespace Dumux
diff --git a/dumux/assembly/numericdifferentiation.hh b/dumux/assembly/numericdifferentiation.hh
new file mode 100644
index 0000000000000000000000000000000000000000..cdf0ee3f1bb810e5a3d1d4f0196e2baad6c87332
--- /dev/null
+++ b/dumux/assembly/numericdifferentiation.hh
@@ -0,0 +1,126 @@
+// -*- 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 2 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
+ * \brief A class for numeric differentiation
+ *
+ */
+#ifndef DUMUX_NUMERIC_DIFFERENTIATION_HH
+#define DUMUX_NUMERIC_DIFFERENTIATION_HH
+
+#include <cmath>
+#include <limits>
+#include <dumux/common/parameters.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Assembly
+ * \brief A class for numeric differentiation with respect to a scalar parameter
+ */
+class NumericDifferentiation
+{
+public:
+
+ /*!
+ * \brief Computes the epsilon used for numeric differentiation
+ * \param value The value of the variable with respect to which we are differentiating
+ */
+ template<class Scalar>
+ static Scalar epsilon(const Scalar value)
+ {
+ static constexpr Scalar baseEps = 1e-10;
+ assert(std::numeric_limits<Scalar>::epsilon()*1e4 < baseEps);
+ // the epsilon value used for the numeric differentiation is
+ // now scaled by the absolute value of the primary variable...
+ using std::abs;
+ return baseEps*(abs(value) + 1.0);
+ }
+
+ /*!
+ * \brief Computes the derivative of a function with repect to a function parameter
+ * \note Overload using default epsilon computation
+ */
+ template<class Function, class Scalar, class FunctionEvalType>
+ static void partialDerivative(const Function& function, Scalar x0,
+ FunctionEvalType& derivative,
+ const FunctionEvalType& fx0)
+ { partialDerivative(function, x0, derivative, fx0, epsilon(x0)); }
+
+ /*!
+ * \brief Computes the derivative of a function with repect to a function parameter
+ * \param function The function to derive
+ * \param x0 The parameter at which the derivative is ought to be evaluated
+ * \param derivative The partial derivative (output)
+ * \param fx0 The result of the function evaluated at x0
+ * \param eps The numeric epsilon used in the differentiation
+ */
+ template<class Function, class Scalar, class FunctionEvalType>
+ static void partialDerivative(const Function& function, Scalar x0,
+ FunctionEvalType& derivative,
+ const FunctionEvalType& fx0,
+ const Scalar eps)
+ {
+ static const int numericDifferenceMethod = numDiffMethod();
+
+ Scalar delta = 0.0;
+ // we are using forward or central differences, i.e. we need to calculate f(x + \epsilon)
+ if (numericDifferenceMethod >= 0)
+ {
+ delta += eps;
+ // calculate the function evaluated with the deflected variable
+ derivative = function(x0 + eps);
+ }
+
+ // we are using backward differences,
+ // i.e. we don't need to calculate f(x + \epsilon)
+ // we can recycle the (possibly cached) f(x)
+ else derivative = fx0;
+
+ // we are using backward or central differences,
+ // i.e. we need to calculate f(x - \epsilon)
+ if (numericDifferenceMethod <= 0)
+ {
+ delta += eps;
+ // subtract the function evaluated with the deflected variable
+ derivative -= function(x0 - eps);
+ }
+
+ // we are using forward differences,
+ // i.e. we don't need to calculate f(x - \epsilon)
+ // we can recycle the (possibly cached) f(x)
+ else derivative -= fx0;
+
+ // divide difference in residuals by the magnitude of the
+ // deflections between the two function evaluation
+ derivative /= delta;
+ }
+
+ static int numDiffMethod()
+ {
+ static const int numericDifferenceMethod
+ = getParam<int>("Assembly.NumericDifferenceMethod", 1);
+ return numericDifferenceMethod;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/assembly/staggeredfvassembler.hh b/dumux/assembly/staggeredfvassembler.hh
index 7042b115c2a7e310ff76590de67a02116d3cc2b3..4b144c880197922f9fdd7d743eae951f8da9f7a7 100644
--- a/dumux/assembly/staggeredfvassembler.hh
+++ b/dumux/assembly/staggeredfvassembler.hh
@@ -128,7 +128,7 @@ public:
try
{
// let the local assembler add the element contributions
- for (const auto element : elements(gridView()))
+ for (const auto& element : elements(gridView()))
LocalAssembler::assembleJacobianAndResidual(*this, *jacobian_, *residual_, element, curSol);
// if we get here, everything worked well
@@ -170,7 +170,7 @@ public:
DUNE_THROW(Dune::InvalidStateException, "Assembling instationary problem but previous solution was not set!");
// let the local assembler add the element contributions
- for (const auto element : elements(gridView()))
+ for (const auto& element : elements(gridView()))
LocalAssembler::assembleResidual(*this, r, element, curSol);
}
@@ -362,6 +362,9 @@ public:
const LocalResidual& localResidual() const
{ return localResidual_; }
+ bool isStationaryProblem() const
+ { return stationary_; }
+
private:
//! reset the residual to 0.0
void resetResidual_()
diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index 2cbc1b8ed7f371161403edffa0f0962cea284e1b..ca9b5c443d63cca5184e0e87b6fb27ce10ef43f3 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -98,6 +98,7 @@ NEW_PROP_TAG(ElementFluxVariablesCache); //!< A local vector of flux v
NEW_PROP_TAG(GridFluxVariablesCache); //!< The global vector of flux variable containers
NEW_PROP_TAG(EnableGridFluxVariablesCache); //!< specifies if data on flux vars should be saved (faster, but more memory consuming)
NEW_PROP_TAG(GridVariables); //!< The grid variables object managing variable data on the grid (volvars/fluxvars cache)
+NEW_PROP_TAG(MaxNumNeighborsPerScvf); //!< The maximum number of neighboring elements allowed per scvf (for static memory allocation)
/////////////////////////////////////////////////////////////////
// Additional properties used by the cell-centered mpfa schemes:
@@ -106,7 +107,7 @@ NEW_PROP_TAG(MpfaMethod); //!< Specifies the mpfa metho
NEW_PROP_TAG(MpfaHelper); //!< A Helper class depending on the mpfa method and grid dimension
NEW_PROP_TAG(PrimaryInteractionVolume); //!< The primary interaction volume type
NEW_PROP_TAG(SecondaryInteractionVolume); //!< The secondary interaction volume type used e.g. on the boundaries
-
+NEW_PROP_TAG(DualGridNodalIndexSet); //!< The type used for the nodal index sets of the dual grid
/////////////////////////////////////////////////////////////
// Properties used by models involving flow in porous media:
diff --git a/dumux/common/reservedblockvector.hh b/dumux/common/reservedblockvector.hh
new file mode 100644
index 0000000000000000000000000000000000000000..277cb10994b0a8972dfffde73fdbc3f58a5a071a
--- /dev/null
+++ b/dumux/common/reservedblockvector.hh
@@ -0,0 +1,96 @@
+// -*- 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 2 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 Common
+ * \brief A arithmetic block vector type based on DUNE's reserved vector
+ */
+#ifndef DUMUX_RESERVED_BLOCK_VECTOR_HH
+#define DUMUX_RESERVED_BLOCK_VECTOR_HH
+
+#include <algorithm>
+#include <dune/common/reservedvector.hh>
+
+namespace Dumux {
+
+template<class BlockType, int capacity>
+class ReservedBlockVector : public Dune::ReservedVector<BlockType, capacity>
+{
+ using Base = Dune::ReservedVector<BlockType, capacity>;
+public:
+
+ using size_type = typename Base::size_type;
+ using value_type = BlockType;
+
+ using Base::Base;
+
+ explicit ReservedBlockVector() : Base() {}
+ explicit ReservedBlockVector(size_type size) : Base() { this->resize(size); }
+
+ ReservedBlockVector(const ReservedBlockVector&) = default;
+ ReservedBlockVector(ReservedBlockVector&&) = default;
+
+ ReservedBlockVector& operator= (const ReservedBlockVector&) = default;
+ ReservedBlockVector& operator= (ReservedBlockVector&&) = default;
+
+ ~ReservedBlockVector() = default;
+
+ //! assigment from scalar
+ ReservedBlockVector& operator= (const typename BlockType::field_type& v)
+ {
+ std::fill(this->begin(), this->end(), v);
+ return *this;
+ }
+
+ //! vector space addition
+ ReservedBlockVector& operator+= (const ReservedBlockVector& other)
+ {
+ for (size_type i = 0; i < this->size(); ++i)
+ (*this)[i] += other[i];
+ return *this;
+ }
+
+ //! vector space subtraction
+ ReservedBlockVector& operator-= (const ReservedBlockVector& other)
+ {
+ for (size_type i = 0; i < this->size(); ++i)
+ (*this)[i] -= other[i];
+ return *this;
+ }
+
+ //! division by scalar
+ ReservedBlockVector& operator/= (const typename BlockType::field_type& v)
+ {
+ for (size_type i = 0; i < this->size(); ++i)
+ (*this)[i] /= v;
+ return *this;
+ }
+
+ //! multiplication by scalar
+ ReservedBlockVector& operator*= (const typename BlockType::field_type& v)
+ {
+ for (size_type i = 0; i < this->size(); ++i)
+ (*this)[i] *= v;
+ return *this;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/discretization/box/fvelementgeometry.hh b/dumux/discretization/box/fvelementgeometry.hh
index e0639746b5c43a45212ddb8ef6d6d29a2da5c1b2..393b8cd8b94d53f454a6289a4b182cef70993c72 100644
--- a/dumux/discretization/box/fvelementgeometry.hh
+++ b/dumux/discretization/box/fvelementgeometry.hh
@@ -54,28 +54,33 @@ class BoxFVElementGeometry
template<class TypeTag>
class BoxFVElementGeometry<TypeTag, true>
{
- using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using IndexType = typename GridView::IndexSet::IndexType;
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+public:
+ //! export type of subcontrol volume
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ //! export type of subcontrol volume face
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using Element = typename GridView::template Codim<0>::Entity;
+ //! export type of finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ //! the maximum number of scvs per element (2^dim for cubes)
+ static constexpr std::size_t maxNumElementScvs = (1<<dim);
+
+private:
+ using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using Element = typename GridView::template Codim<0>::Entity;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using CoordScalar = typename GridView::ctype;
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
-
using FeCache = Dune::PQkLocalFiniteElementCache<CoordScalar, Scalar, dim, 1>;
using FeLocalBasis = typename FeCache::FiniteElementType::Traits::LocalBasisType;
using ReferenceElements = typename Dune::ReferenceElements<CoordScalar, dim>;
- using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>;
- using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
-
public:
+
//! Constructor
BoxFVElementGeometry(const FVGridGeometry& fvGridGeometry)
: fvGridGeometryPtr_(&fvGridGeometry) {}
@@ -169,16 +174,22 @@ template<class TypeTag>
class BoxFVElementGeometry<TypeTag, false>
{
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using IndexType = typename GridView::IndexSet::IndexType;
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+public:
+ //! export type of subcontrol volume
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ //! export type of subcontrol volume face
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using LocalIndexType = typename SubControlVolumeFace::Traits::LocalIndexType;
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ //! export type of finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ //! the maximum number of scvs per element (2^dim for cubes)
+ static constexpr std::size_t maxNumElementScvs = (1<<dim);
- static const int dim = GridView::dimension;
- static const int dimWorld = GridView::dimensionworld;
-
+private:
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using LocalIndexType = typename SubControlVolumeFace::Traits::LocalIndexType;
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using Element = typename GridView::template Codim<0>::Entity;
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
diff --git a/dumux/discretization/box/gridvolumevariables.hh b/dumux/discretization/box/gridvolumevariables.hh
index 7875b88304c981892fe6a4fa40bbd779882c2770..b4aa07fd2645d209dfd3bbf60e01528926387d96 100644
--- a/dumux/discretization/box/gridvolumevariables.hh
+++ b/dumux/discretization/box/gridvolumevariables.hh
@@ -86,6 +86,12 @@ public:
friend inline ElementVolumeVariables localView(const BoxGridVolumeVariables& global)
{ return ElementVolumeVariables(global); }
+ const VolumeVariables& volVars(const SubControlVolume& scv) const
+ { return volumeVariables_[scv.elementIndex()][scv.indexInElement()]; }
+
+ VolumeVariables& volVars(const SubControlVolume& scv)
+ { return volumeVariables_[scv.elementIndex()][scv.indexInElement()]; }
+
const VolumeVariables& volVars(const IndexType eIdx, const IndexType scvIdx) const
{ return volumeVariables_[eIdx][scvIdx]; }
diff --git a/dumux/discretization/cellcentered/connectivitymap.hh b/dumux/discretization/cellcentered/connectivitymap.hh
index f675817c404424acc13afee6e9a6b060bfd20e1c..0f935b288ec93fa1040d5f4121963ba02a01cec2 100644
--- a/dumux/discretization/cellcentered/connectivitymap.hh
+++ b/dumux/discretization/cellcentered/connectivitymap.hh
@@ -30,6 +30,7 @@
#include <utility>
#include <algorithm>
+#include <dune/common/reservedvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/discretization/fluxstencil.hh>
@@ -49,6 +50,7 @@ template<class TypeTag>
class CCSimpleConnectivityMap
{
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
using FluxStencil = Dumux::FluxStencil<TypeTag>;
@@ -56,10 +58,10 @@ class CCSimpleConnectivityMap
struct DataJ
{
IndexType globalJ;
- std::vector<IndexType> scvfsJ;
+ Dune::ReservedVector<IndexType, FluxStencil::maxNumScvfJForI> scvfsJ;
// A list of additional scvfs is needed for compatibility
// reasons with more complex connectivity maps (see mpfa)
- std::vector<IndexType> additionalScvfs;
+ Dune::ReservedVector<IndexType, FluxStencil::maxNumScvfJForI> additionalScvfs;
};
using Map = std::vector<std::vector<DataJ>>;
@@ -75,16 +77,21 @@ public:
{
map_.clear();
map_.resize(fvGridGeometry.gridView().size(0));
+
+ // container to store for each element J the elements I that appear in J's flux stencils
+ static constexpr int maxNumJ = FluxStencil::maxFluxStencilSize*FVElementGeometry::maxNumElementScvfs;
+ Dune::ReservedVector<std::pair<IndexType, DataJ>, maxNumJ> dataJForI;
+
for (const auto& element : elements(fvGridGeometry.gridView()))
{
// We are looking for the elements I, for which this element J is in the flux stencil
- auto globalJ = fvGridGeometry.elementMapper().index(element);
+ const auto globalJ = fvGridGeometry.elementMapper().index(element);
auto fvGeometry = localView(fvGridGeometry);
fvGeometry.bindElement(element);
// obtain the data of J in elements I
- std::vector<std::pair<IndexType, DataJ>> dataJForI;
+ dataJForI.clear();
// loop over sub control faces
for (auto&& scvf : scvfs(fvGeometry))
@@ -97,16 +104,13 @@ public:
if (globalI == globalJ)
continue;
- auto it = std::find_if(dataJForI.begin(),
- dataJForI.end(),
+ auto it = std::find_if(dataJForI.begin(), dataJForI.end(),
[globalI](const auto& pair) { return pair.first == globalI; });
if (it != dataJForI.end())
it->second.scvfsJ.push_back(scvf.index());
else
- dataJForI.emplace_back(std::make_pair(globalI, DataJ({globalJ,
- std::vector<IndexType>({scvf.index()}),
- std::vector<IndexType>()})));
+ dataJForI.push_back(std::make_pair(globalI, DataJ({globalJ, {scvf.index()}, {}})));
}
}
diff --git a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
index 0d7f2a7be9c275cea7faf65ef00b5d8562e03ea2..9ff81414195549d3849e38f683cee513a8b94a82 100644
--- a/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/darcyslaw.hh
@@ -94,6 +94,9 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
static constexpr int dimWorld = GridView::dimensionworld;
static constexpr int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
+ using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
+ using Stencil = typename DualGridNodalIndexSet::GridStencilType;
+
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs();
@@ -109,7 +112,6 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
using PrimaryIvTij = typename PrimaryIvMatrix::row_type;
- using PrimaryIvStencil = typename PrimaryInteractionVolume::Traits::Stencil;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
@@ -117,7 +119,6 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
using SecondaryIvTij = typename SecondaryIvMatrix::row_type;
- using SecondaryIvStencil = typename SecondaryInteractionVolume::Traits::Stencil;
public:
// export the filler type
@@ -138,7 +139,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SubControlVolumeFace &scvf)
{
switchFluxSign_ = localFaceData.isOutside();
- primaryStencil_ = &iv.stencil();
+ stencil_ = &iv.stencil();
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
primaryPj_[pIdx] = &dataHandle.pressures(pIdx);
@@ -191,7 +192,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SubControlVolumeFace &scvf)
{
switchFluxSign_ = localFaceData.isOutside();
- secondaryStencil_ = &iv.stencil();
+ stencil_ = &iv.stencil();
for (unsigned int pIdx = 0; pIdx < numPhases; ++pIdx)
secondaryPj_[pIdx] = &dataHandle.pressures(pIdx);
@@ -229,10 +230,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
}
//! The stencil corresponding to the transmissibilities (primary type)
- const PrimaryIvStencil& advectionStencilPrimaryIv() const { return *primaryStencil_; }
-
- //! The stencil corresponding to the transmissibilities (secondary type)
- const SecondaryIvStencil& advectionStencilSecondaryIv() const { return *secondaryStencil_; }
+ const Stencil& advectionStencil() const { return *stencil_; }
//! Coefficients for the cell (& Dirichlet) unknowns in flux expressions (primary type)
const PrimaryIvTij& advectionTijPrimaryIv() const { return *primaryTij_; }
@@ -259,8 +257,7 @@ class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
bool switchFluxSign_;
//! The stencil, i.e. the grid indices j
- const PrimaryIvStencil* primaryStencil_;
- const SecondaryIvStencil* secondaryStencil_;
+ const Stencil* stencil_;
//! The transmissibilities such that f = Tij*pj
const PrimaryIvTij* primaryTij_;
diff --git a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
index 52817810a3787d72f3259f47afd5935c80cae7d2..a93c7ce63970a589ca92cdf883dd9547556e1ee2 100644
--- a/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
+++ b/dumux/discretization/cellcentered/mpfa/dualgridindexset.hh
@@ -28,6 +28,8 @@
#include <vector>
#include <algorithm>
+#include <dune/common/reservedvector.hh>
+
namespace Dumux
{
@@ -39,20 +41,40 @@ namespace Dumux
* \tparam GI The type used for indices on the grid
* \tparam LI The type used for indexing in interaction volumes
* \tparam dim The dimension of the grid
+ * \tparam maxE The maximum admissible number of elements around vertices.
+ * \tparam maxB The maximum admissible number of branches on intersections.
+ * This is only to be specified for network grids and defaults to 1
+ * for normal grids.
*/
-template< class GI, class LI, int dim >
-class DualGridNodalIndexSet
+template< class GI, class LI, int dim, int maxE, int maxB = 2 >
+class CCMpfaDualGridNodalIndexSet
{
+ using DimIndexVector = Dune::ReservedVector<GI, dim>;
+
public:
+ //! Export the used index types
using GridIndexType = GI;
using LocalIndexType = LI;
- // we use dynamic containers to store indices here
- using GridIndexContainer = std::vector<GridIndexType>;
- using LocalIndexContainer = std::vector<LocalIndexType>;
+ //! Export the specified maximum admissible sizes
+ static constexpr int dimension = dim;
+ static constexpr int maxBranches = maxB;
+ static constexpr int maxNumElementsAtNode = maxE*(maxBranches-1);
+ static constexpr int maxNumScvfsAtNode = maxNumElementsAtNode*dim;
+
+ //! Export the stencil types used
+ using GridStencilType = Dune::ReservedVector<GridIndexType, maxNumElementsAtNode>;
+ using LocalStencilType = Dune::ReservedVector<LocalIndexType, maxNumElementsAtNode>;
+
+ //! Export the type used for storing the global scvf indices at this node
+ using GridScvfStencilType = Dune::ReservedVector<GridIndexType, maxNumScvfsAtNode>;
+
+ //! Data structure to store the neighboring scv indices of an scvf (grid/local indices)
+ using ScvfNeighborIndexSet = Dune::ReservedVector<GridIndexType, maxBranches>;
+ using ScvfNeighborLocalIndexSet = Dune::ReservedVector<LocalIndexType, maxBranches>;
//! Constructor
- DualGridNodalIndexSet() : numBoundaryScvfs_(0) {}
+ CCMpfaDualGridNodalIndexSet() : numBoundaryScvfs_(0) {}
//! Inserts data for a given scvf
template<typename SubControlVolumeFace>
@@ -78,10 +100,10 @@ public:
const LocalIndexType curScvfLocalIdx = scvfIndices_.size();
// add grid scvf data
- GridIndexContainer scvIndices;
- scvIndices.reserve(outsideScvIndices.size()+1);
+ ScvfNeighborIndexSet scvIndices;
scvIndices.push_back(insideScvIdx);
- scvIndices.insert(scvIndices.end(), outsideScvIndices.begin(), outsideScvIndices.end());
+ for (const auto& outsideIdx : outsideScvIndices)
+ scvIndices.push_back(outsideIdx);
// if scvf is on boundary, increase counter
if (boundary)
@@ -90,7 +112,7 @@ public:
// insert data on the new scv
scvfIndices_.push_back(scvfIdx);
scvfIsOnBoundary_.push_back(boundary);
- scvfNeighborScvIndices_.emplace_back( std::move(scvIndices) );
+ scvfNeighborScvIndices_.push_back(scvIndices);
// if entry for the inside scv exists, append scvf local index, create otherwise
auto it = std::find( scvIndices_.begin(), scvIndices_.end(), insideScvIdx );
@@ -98,10 +120,8 @@ public:
localScvfIndicesInScv_[ std::distance(scvIndices_.begin(), it) ].push_back(curScvfLocalIdx);
else
{
- LocalIndexContainer localScvfIndices;
- localScvfIndices.reserve(dim);
- localScvfIndices.push_back(curScvfLocalIdx);
- localScvfIndicesInScv_.emplace_back( std::move(localScvfIndices) );
+ localScvfIndicesInScv_.push_back({});
+ localScvfIndicesInScv_.back().push_back(curScvfLocalIdx);
scvIndices_.push_back(insideScvIdx);
}
}
@@ -116,10 +136,10 @@ public:
std::size_t numBoundaryScvfs() const { return numBoundaryScvfs_; }
//! returns the global scv indices connected to this dual grid node
- const GridIndexContainer& globalScvIndices() const { return scvIndices_; }
+ const GridStencilType& globalScvIndices() const { return scvIndices_; }
//! returns the global scvf indices connected to this dual grid node
- const GridIndexContainer& globalScvfIndices() const { return scvfIndices_; }
+ const GridScvfStencilType& globalScvfIndices() const { return scvfIndices_; }
//! returns whether or not the i-th scvf is on a domain boundary
bool scvfIsOnBoundary(unsigned int i) const { return scvfIsOnBoundary_[i]; }
@@ -145,39 +165,40 @@ public:
}
//! returns the indices of the neighboring scvs of the i-th scvf
- const GridIndexContainer& neighboringScvIndices(unsigned int i) const
+ const ScvfNeighborIndexSet& neighboringScvIndices(unsigned int i) const
{ return scvfNeighborScvIndices_[i]; }
private:
- GridIndexContainer scvIndices_; //!< The indices of the scvs around a dual grid node
- std::vector<LocalIndexContainer> localScvfIndicesInScv_; //!< Maps to each scv a list of scvf indices embedded in it
+ GridStencilType scvIndices_; //!< The indices of the scvs around a dual grid node
+ Dune::ReservedVector<DimIndexVector, maxNumElementsAtNode> localScvfIndicesInScv_; //!< Maps to each scv a list of scvf indices embedded in it
+
+ GridScvfStencilType scvfIndices_; //!< the indices of the scvfs around a dual grid node
+ std::size_t numBoundaryScvfs_; //!< stores how many boundary scvfs are embedded in this dual grid node
+ Dune::ReservedVector<bool, maxNumScvfsAtNode> scvfIsOnBoundary_; //!< Maps to each scvf a boolean to indicate if it is on the boundary
- GridIndexContainer scvfIndices_; //!< the indices of the scvfs around a dual grid node
- std::vector<bool> scvfIsOnBoundary_; //!< Maps to each scvf a boolean to indicate if it is on the boundary
- std::vector<GridIndexContainer> scvfNeighborScvIndices_; //!< The neighboring scvs for the scvfs (order: 0 - inside, 1..n - outside)
- std::size_t numBoundaryScvfs_; //!< stores how many boundary scvfs are embedded in this dual grid node
+ //! The neighboring scvs for the scvfs (order: 0 - inside, 1..n - outside)
+ Dune::ReservedVector<ScvfNeighborIndexSet, maxNumScvfsAtNode> scvfNeighborScvIndices_;
};
/*!
* \ingroup CCMpfaDiscretization
* \brief Class for the index sets of the dual grid in mpfa schemes.
*
- * \tparam GridIndexType The type used for indices on the grid
- * \tparam LocalIndexType The type used for indexing in interaction volumes
- * \tparam dim The dimension of the grid
+ * \tparam NI The type used for the nodal index sets.
*/
-template< class GridIndexType, class LocalIndexType, int dim >
+template< class NI >
class CCMpfaDualGridIndexSet
{
public:
- using NodalIndexSet = DualGridNodalIndexSet< GridIndexType, LocalIndexType, dim >;
+ using NodalIndexSet = NI;
+ using GridIndexType = typename NodalIndexSet::GridIndexType;
//! Default constructor should not be used
CCMpfaDualGridIndexSet() = delete;
//! Constructor taking a grid view
template< class GridView >
- CCMpfaDualGridIndexSet(const GridView& gridView) : nodalIndexSets_(gridView.size(dim)) {}
+ CCMpfaDualGridIndexSet(const GridView& gridView) : nodalIndexSets_(gridView.size(NodalIndexSet::dimension)) {}
//! Access with an scvf
template< class SubControlVolumeFace >
diff --git a/dumux/discretization/cellcentered/mpfa/fickslaw.hh b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
index be10d54f81548e4e79bf883b6cffde8bf49ad9c2..1f3ea99206a341aff67ccd83bdc78c8ca4ba6d34 100644
--- a/dumux/discretization/cellcentered/mpfa/fickslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fickslaw.hh
@@ -90,6 +90,9 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
//! The cache used in conjunction with the mpfa Fick's Law
class MpfaFicksLawCache
{
+ using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
+ using Stencil = typename DualGridNodalIndexSet::GridStencilType;
+
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs();
@@ -105,7 +108,6 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
using PrimaryIvTij = typename PrimaryIvMatrix::row_type;
- using PrimaryIvStencil = typename PrimaryInteractionVolume::Traits::Stencil;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
@@ -113,7 +115,6 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
using SecondaryIvTij = typename SecondaryIvMatrix::row_type;
- using SecondaryIvStencil = typename SecondaryInteractionVolume::Traits::Stencil;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
@@ -138,7 +139,7 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SubControlVolumeFace &scvf,
unsigned int phaseIdx, unsigned int compIdx)
{
- primaryStencil_[phaseIdx][compIdx] = &iv.stencil();
+ stencil_[phaseIdx][compIdx] = &iv.stencil();
switchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutside();
// store pointer to the mole fraction vector of this iv
@@ -168,7 +169,7 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SubControlVolumeFace &scvf,
unsigned int phaseIdx, unsigned int compIdx)
{
- secondaryStencil_[phaseIdx][compIdx] = &iv.stencil();
+ stencil_[phaseIdx][compIdx] = &iv.stencil();
switchFluxSign_[phaseIdx][compIdx] = localFaceData.isOutside();
// store pointer to the mole fraction vector of this iv
@@ -205,20 +206,15 @@ class FicksLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SecondaryIvVector& moleFractionsSecondaryIv(unsigned int phaseIdx, unsigned int compIdx) const
{ return *secondaryXj_[phaseIdx][compIdx]; }
- //! The stencils corresponding to the transmissibilities (primary type)
- const PrimaryIvStencil& diffusionStencilPrimaryIv(unsigned int phaseIdx, unsigned int compIdx) const
- { return *primaryStencil_[phaseIdx][compIdx]; }
-
- //! The stencils corresponding to the transmissibilities (secondary type)
- const SecondaryIvStencil& diffusionStencilSecondaryIv(unsigned int phaseIdx, unsigned int compIdx) const
- { return *secondaryStencil_[phaseIdx][compIdx]; }
+ //! The stencils corresponding to the transmissibilities
+ const Stencil& diffusionStencil(unsigned int phaseIdx, unsigned int compIdx) const
+ { return *stencil_[phaseIdx][compIdx]; }
private:
std::array< std::array<bool, numComponents>, numPhases > switchFluxSign_;
//! The stencils, i.e. the grid indices j
- std::array< std::array<const PrimaryIvStencil*, numComponents>, numPhases > primaryStencil_;
- std::array< std::array<const SecondaryIvStencil*, numComponents>, numPhases > secondaryStencil_;
+ std::array< std::array<const Stencil*, numComponents>, numPhases > stencil_;
//! The transmissibilities such that f = Tij*xj
std::array< std::array<const PrimaryIvVector*, numComponents>, numPhases > primaryTij_;
diff --git a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
index 11e2a4accd03b1fa82106cc6efae2618ee8650fa..ebdb9d76aa4e3e286acd5d59142653bf0158e81f 100644
--- a/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
+++ b/dumux/discretization/cellcentered/mpfa/fourierslaw.hh
@@ -88,6 +88,9 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
//! The cache used in conjunction with the mpfa Fourier's Law
class MpfaFouriersLawCache
{
+ using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
+ using Stencil = typename DualGridNodalIndexSet::GridStencilType;
+
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
static constexpr bool considerSecondaryIVs = MpfaHelper::considerSecondaryIVs();
@@ -103,7 +106,6 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using PrimaryIvVector = typename PrimaryInteractionVolume::Traits::Vector;
using PrimaryIvMatrix = typename PrimaryInteractionVolume::Traits::Matrix;
using PrimaryIvTij = typename PrimaryIvMatrix::row_type;
- using PrimaryIvStencil = typename PrimaryInteractionVolume::Traits::Stencil;
using SecondaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, SecondaryInteractionVolume);
using SecondaryIvLocalFaceData = typename SecondaryInteractionVolume::Traits::LocalFaceData;
@@ -111,7 +113,6 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using SecondaryIvVector = typename SecondaryInteractionVolume::Traits::Vector;
using SecondaryIvMatrix = typename SecondaryInteractionVolume::Traits::Matrix;
using SecondaryIvTij = typename SecondaryIvMatrix::row_type;
- using SecondaryIvStencil = typename SecondaryInteractionVolume::Traits::Stencil;
static constexpr int dim = GridView::dimension;
static constexpr int dimWorld = GridView::dimensionworld;
@@ -134,7 +135,7 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const PrimaryIvDataHandle& dataHandle,
const SubControlVolumeFace &scvf)
{
- primaryStencil_ = &iv.stencil();
+ stencil_ = &iv.stencil();
switchFluxSign_ = localFaceData.isOutside();
// store pointer to the temperature vector of this iv
@@ -163,7 +164,7 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SecondaryIvDataHandle& dataHandle,
const SubControlVolumeFace &scvf)
{
- secondaryStencil_ = &iv.stencil();
+ stencil_ = &iv.stencil();
switchFluxSign_ = localFaceData.isOutside();
// store pointer to the temperature vector of this iv
@@ -184,10 +185,7 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
const SecondaryIvTij& heatConductionTijSecondaryIv() const { return *secondaryTij_; }
//! The stencil corresponding to the transmissibilities (primary type)
- const PrimaryIvStencil& heatConductionStencilPrimaryIv() const { return *primaryStencil_; }
-
- //! The stencil corresponding to the transmissibilities (secondary type)
- const SecondaryIvStencil& heatConductionStencilSecondaryIv() const { return *secondaryStencil_; }
+ const Stencil& heatConductionStencil() const { return *stencil_; }
//! The cell (& Dirichlet) temperatures within this interaction volume (primary type)
const PrimaryIvVector& temperaturesPrimaryIv() const { return *primaryTj_; }
@@ -205,8 +203,7 @@ class FouriersLawImplementation<TypeTag, DiscretizationMethods::CCMpfa>
bool switchFluxSign_;
//! The stencil, i.e. the grid indices j
- const PrimaryIvStencil* primaryStencil_;
- const SecondaryIvStencil* secondaryStencil_;
+ const Stencil* stencil_;
//! The transmissibilities such that f = Tij*Tj
const PrimaryIvTij* primaryTij_;
diff --git a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
index f3fbca56476ed92c90b1194dbb6e5c5f1e1e73ff..44b8f899fe6fdb14221087ea24a8bd92bf36fb6b 100644
--- a/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvelementgeometry.hh
@@ -59,16 +59,22 @@ class CCMpfaFVElementGeometry<TypeTag, true>
using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Element = typename GridView::template Codim<0>::Entity;
-
using GridIndexType = typename GridView::IndexSet::IndexType;
+
+ static constexpr int dim = GridView::dimension;
+
+public:
+ //! export type of subcontrol volume
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ //! export type of subcontrol volume face
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ //! export type of finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ //! the maximum number of scvs per element
+ static constexpr std::size_t maxNumElementScvs = 1;
+ //! the maximum number of scvfs per element (use cubes for maximum)
+ static constexpr std::size_t maxNumElementScvfs = dim == 3 ? 24 : 8;
- using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<GridIndexType>, ThisType>;
- using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType>;
-
-public:
//! Constructor
CCMpfaFVElementGeometry(const FVGridGeometry& fvGridGeometry)
: fvGridGeometryPtr_(&fvGridGeometry) {}
@@ -97,9 +103,10 @@ public:
//! This is a free function found by means of ADL
//! To iterate over all sub control volumes of this FVElementGeometry use
//! for (auto&& scv : scvs(fvGeometry))
- friend inline Dune::IteratorRange<ScvIterator>
+ friend inline Dune::IteratorRange< ScvIterator<SubControlVolume, std::array<GridIndexType, 1>, ThisType> >
scvs(const CCMpfaFVElementGeometry& fvGeometry)
{
+ using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::array<GridIndexType, 1>, ThisType>;
return Dune::IteratorRange<ScvIterator>(ScvIterator(fvGeometry.scvIndices_.begin(), fvGeometry),
ScvIterator(fvGeometry.scvIndices_.end(), fvGeometry));
}
@@ -109,11 +116,12 @@ public:
//! This is a free function found by means of ADL
//! To iterate over all sub control volume faces of this FVElementGeometry use
//! for (auto&& scvf : scvfs(fvGeometry))
- friend inline Dune::IteratorRange<ScvfIterator>
+ friend inline Dune::IteratorRange< ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType> >
scvfs(const CCMpfaFVElementGeometry& fvGeometry)
{
const auto& g = fvGeometry.fvGridGeometry();
const auto scvIdx = fvGeometry.scvIndices_[0];
+ using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType>;
return Dune::IteratorRange<ScvfIterator>(ScvfIterator(g.scvfIndicesOfScv(scvIdx).begin(), fvGeometry),
ScvfIterator(g.scvfIndicesOfScv(scvIdx).end(), fvGeometry));
}
@@ -139,7 +147,7 @@ public:
//! Bind only element-local
void bindElement(const Element& element)
{
- scvIndices_ = std::vector<GridIndexType>({fvGridGeometry().elementMapper().index(element)});
+ scvIndices_[0] = fvGridGeometry().elementMapper().index(element);
}
//! The global finite volume geometry we are a restriction of
@@ -148,7 +156,7 @@ public:
private:
- std::vector<GridIndexType> scvIndices_;
+ std::array<GridIndexType, 1> scvIndices_;
const FVGridGeometry* fvGridGeometryPtr_;
};
@@ -168,12 +176,6 @@ class CCMpfaFVElementGeometry<TypeTag, false>
using GridIndexType = typename GridView::IndexSet::IndexType;
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
-
- using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<GridIndexType>, ThisType>;
- using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<GridIndexType>, ThisType>;
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
@@ -182,6 +184,17 @@ class CCMpfaFVElementGeometry<TypeTag, false>
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
public:
+ //! export type of subcontrol volume
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ //! export type of subcontrol volume face
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ //! export type of finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ //! the maximum number of scvs per element
+ static constexpr std::size_t maxNumElementScvs = 1;
+ //! the maximum number of scvfs per element (use cubes for maximum)
+ static constexpr std::size_t maxNumElementScvfs = dim == 3 ? 24 : 8;
+
//! Constructor
CCMpfaFVElementGeometry(const FVGridGeometry& fvGridGeometry)
: fvGridGeometryPtr_(&fvGridGeometry) {}
@@ -234,11 +247,11 @@ public:
//! This is a free function found by means of ADL
//! To iterate over all sub control volumes of this FVElementGeometry use
//! for (auto&& scv : scvs(fvGeometry))
- friend inline Dune::IteratorRange<typename std::vector<SubControlVolume>::const_iterator>
+ friend inline Dune::IteratorRange<typename std::array<SubControlVolume, 1>::const_iterator>
scvs(const ThisType& g)
{
- using Iter = typename std::vector<SubControlVolume>::const_iterator;
- return Dune::IteratorRange<Iter>(g.scvs_.begin(), g.scvs_.end());
+ using IteratorType = typename std::array<SubControlVolume, 1>::const_iterator;
+ return Dune::IteratorRange<IteratorType>(g.scvs_.begin(), g.scvs_.end());
}
//! iterator range for sub control volumes faces. Iterates over
@@ -249,8 +262,8 @@ public:
friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator>
scvfs(const ThisType& g)
{
- using Iter = typename std::vector<SubControlVolumeFace>::const_iterator;
- return Dune::IteratorRange<Iter>(g.scvfs_.begin(), g.scvfs_.end());
+ using IteratorType = typename std::vector<SubControlVolumeFace>::const_iterator;
+ return Dune::IteratorRange<IteratorType>(g.scvfs_.begin(), g.scvfs_.end());
}
//! number of sub control volumes in this fv element geometry
@@ -336,8 +349,8 @@ private:
{
// make the scv
const auto eIdx = fvGridGeometry().elementMapper().index(element);
- scvs_.emplace_back(element.geometry(), eIdx);
- scvIndices_.emplace_back(eIdx);
+ scvs_[0] = SubControlVolume(element.geometry(), eIdx);
+ scvIndices_[0] = eIdx;
// get data on the scv faces
const auto& scvFaceIndices = fvGridGeometry().scvfIndicesOfScv(eIdx);
@@ -617,9 +630,7 @@ private:
//! clear all containers
void clear()
{
- scvIndices_.clear();
scvfIndices_.clear();
- scvs_.clear();
scvfs_.clear();
neighborScvIndices_.clear();
@@ -634,9 +645,9 @@ private:
const FVGridGeometry* fvGridGeometryPtr_;
// local storage after binding an element
- std::vector<GridIndexType> scvIndices_;
+ std::array<GridIndexType, 1> scvIndices_;
std::vector<GridIndexType> scvfIndices_;
- std::vector<SubControlVolume> scvs_;
+ std::array<SubControlVolume, 1> scvs_;
std::vector<SubControlVolumeFace> scvfs_;
std::vector<GridIndexType> neighborScvIndices_;
diff --git a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
index 0b1e09f5d9722e5a3babc41ab78ad4a2b3adb46d..dfba3592e209700e45d4642e397f517b5433f180 100644
--- a/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
+++ b/dumux/discretization/cellcentered/mpfa/fvgridgeometry.hh
@@ -64,6 +64,7 @@ class CCMpfaFVGridGeometry<TypeTag, true> : public BaseFVGridGeometry<TypeTag>
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
+ using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
@@ -164,7 +165,7 @@ public:
const auto isGhostVertex = MpfaHelper::findGhostVertices(this->gridView(), this->vertexMapper());
// instantiate the dual grid index set (to be used for construction of interaction volumes)
- CCMpfaDualGridIndexSet<GridIndexType, LocalIndexType, dim> dualIdSet(this->gridView());
+ CCMpfaDualGridIndexSet< DualGridNodalIndexSet > dualIdSet(this->gridView());
// Build the SCVs and SCV faces
GridIndexType scvfIdx = 0;
@@ -405,6 +406,7 @@ class CCMpfaFVGridGeometry<TypeTag, false> : public BaseFVGridGeometry<TypeTag>
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
using PrimaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
+ using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
@@ -505,7 +507,7 @@ public:
isGhostVertex_ = MpfaHelper::findGhostVertices(this->gridView(), this->vertexMapper());
// instantiate the dual grid index set (to be used for construction of interaction volumes)
- CCMpfaDualGridIndexSet<GridIndexType, LocalIndexType, dim> dualIdSet(this->gridView());
+ CCMpfaDualGridIndexSet< DualGridNodalIndexSet > dualIdSet(this->gridView());
// Build the SCVs and SCV faces
numScvf_ = 0;
diff --git a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
index 5a7a7882f0ddfcfe07e9b4ed0005244d9e0ef336..7ace724febac7f73d204653a76a66ea72fc4ac84 100644
--- a/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
+++ b/dumux/discretization/cellcentered/mpfa/gridinteractionvolumeindexsets.hh
@@ -43,6 +43,7 @@ class CCMpfaGridInteractionVolumeIndexSets
using GridIndexType = typename GridView::IndexSet::IndexType;
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using DualGridNodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
using PrimaryIV = typename GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume);
using PrimaryIVIndexSet = typename PrimaryIV::Traits::IndexSet;
@@ -51,7 +52,7 @@ class CCMpfaGridInteractionVolumeIndexSets
static constexpr int dim = GridView::dimension;
using LocalIndexType = typename PrimaryIV::Traits::LocalIndexType;
- using DualGridIndexSet = CCMpfaDualGridIndexSet< GridIndexType, LocalIndexType, dim>;
+ using DualGridIndexSet = CCMpfaDualGridIndexSet< DualGridNodalIndexSet >;
public:
/*!
diff --git a/dumux/discretization/cellcentered/mpfa/helper.hh b/dumux/discretization/cellcentered/mpfa/helper.hh
index f129138540b5be817f3789ac4880a4f75802f17c..5a32421a9ed42b371e90f08cf5fd5fae8155fb0c 100644
--- a/dumux/discretization/cellcentered/mpfa/helper.hh
+++ b/dumux/discretization/cellcentered/mpfa/helper.hh
@@ -654,8 +654,8 @@ public:
{ return !std::is_same<PrimaryInteractionVolume, SecondaryInteractionVolume>::value; }
//! returns whether or not a value exists in a vector
- template<typename V1, typename V2>
- static bool vectorContainsValue(const std::vector<V1>& vector, const V2 value)
+ template<typename V, typename T>
+ static bool vectorContainsValue(const V& vector, const T value)
{ return std::find(vector.begin(), vector.end(), value) != vector.end(); }
};
diff --git a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
index 51ac0be10eebe29b0c12ce5c8df5480725f1c059..d5e46c9b0ff4a71744be6ffe45a8dd158bcac447 100644
--- a/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
+++ b/dumux/discretization/cellcentered/mpfa/interactionvolumebase.hh
@@ -65,8 +65,6 @@ namespace Dumux
* using Matrix = ...;
* //! export the type used for iv-local vectors
* using Vector = ...;
- * //! export the type used for the iv-stencils
- * using Stencil = ...;
* //! export the data handle type for this iv
* using DataHandle = ...;
* \endcode
@@ -123,7 +121,7 @@ public:
{ DUNE_THROW(Dune::NotImplemented, "Interaction volume implementation does not provide a localFaceData() funtion"); }
//! returns the cell-stencil of this interaction volume
- const typename Traits::Stencil& stencil() const { return asImp().stencil(); }
+ const typename Traits::IndexSet::GridStencilType& stencil() const { return asImp().stencil(); }
//! returns the local scvf entity corresponding to a given iv-local scvf idx
const typename Traits::LocalScvfType& localScvf(typename Traits::LocalIndexType ivLocalScvfIdx) const
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
index 298ec506bd6461955c8b68fe6ee4c25d3417d8ab..06fefbd08a6d351b67beb5e264c78e17b3a76743 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh
@@ -51,7 +51,8 @@ template< class TypeTag >
struct CCMpfaODefaultInteractionVolumeTraits
{
private:
- using GridIndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
+ using NodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
+ using GridIndexType = typename NodalIndexSet::GridIndexType;
static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
@@ -69,9 +70,9 @@ public:
//! export the type of the mpfa helper class
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
//! export the type used for local indices
- using LocalIndexType = std::uint8_t;
+ using LocalIndexType = typename NodalIndexSet::LocalIndexType;
//! export the type for the interaction volume index set
- using IndexSet = CCMpfaOInteractionVolumeIndexSet< DualGridNodalIndexSet<GridIndexType, LocalIndexType, dim> >;
+ using IndexSet = CCMpfaOInteractionVolumeIndexSet< NodalIndexSet >;
//! export the type of interaction-volume local scvs
using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, ScalarType, dim, dimWorld >;
//! export the type of interaction-volume local scvfs
@@ -82,8 +83,6 @@ public:
using Matrix = Dune::DynamicMatrix< ScalarType >;
//! export the type used for iv-local vectors
using Vector = Dune::DynamicVector< ScalarType >;
- //! export the type used for the iv-stencils
- using Stencil = std::vector< GridIndexType >;
//! export the data handle type for this iv
using DataHandle = InteractionVolumeDataHandle< TypeTag, Matrix, Vector, LocalIndexType >;
};
@@ -119,15 +118,7 @@ class CCMpfaOInteractionVolume : public CCMpfaInteractionVolumeBase< CCMpfaOInte
using IndexSet = typename Traits::IndexSet;
using GridIndexType = typename GridView::IndexSet::IndexType;
using LocalIndexType = typename Traits::LocalIndexType;
- using Stencil = typename Traits::Stencil;
-
- //! For the o method, the interaction volume stencil can be taken directly
- //! from the nodal index set, which always uses dynamic types to be compatible
- //! on the boundaries and unstructured grids. Thus, we have to make sure that
- //! the type set for the stencils in the traits is castable.
- static_assert( std::is_convertible<Stencil*, typename IndexSet::GridIndexContainer*>::value,
- "The o-method uses the (dynamic) nodal index set's stencil as the interaction volume stencil. "
- "Using a different type is not permissive here." );
+ using Stencil = typename IndexSet::GridStencilType;
//! Data attached to scvf touching Dirichlet boundaries.
//! For the default o-scheme, we only store the corresponding vol vars index.
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
index 84a26f7a4e4a253f905b1ca6405f6b1e2a2bdd65..1d63d804f259675a3619111b2056aa24dcd9430f 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/interactionvolumeindexset.hh
@@ -24,6 +24,8 @@
#ifndef DUMUX_DISCRETIZATION_MPFA_O_INTERACTIONVOLUME_INDEXSET_HH
#define DUMUX_DISCRETIZATION_MPFA_O_INTERACTIONVOLUME_INDEXSET_HH
+#include <dune/common/reservedvector.hh>
+
#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh>
namespace Dumux
@@ -38,46 +40,27 @@ template< class DualGridNodalIndexSet >
class CCMpfaOInteractionVolumeIndexSet
{
public:
+ //! Export the type used for the nodal grid index sets
+ using NodalIndexSet = DualGridNodalIndexSet;
+
+ //! Export the types used for local/grid indices
using LocalIndexType = typename DualGridNodalIndexSet::LocalIndexType;
using GridIndexType = typename DualGridNodalIndexSet::GridIndexType;
- using LocalIndexContainer = typename DualGridNodalIndexSet::LocalIndexContainer;
- using GridIndexContainer = typename DualGridNodalIndexSet::GridIndexContainer;
-
- /*!
- * \brief The constructor
- * \note The actual type used for the nodal index sets might be different, as maybe
- * a different type for the local indexes is used. We therefore template this
- * constructor. However, a static assertion enforces you to use the same LocalIndexType
- * in the traits for both the secondary and the primary interaction volume traits.
- *
- * \tparam NodalIndexSet Possibly differing type for the DualGridNodalIndexSet
- */
- template< class NodalIndexSet >
- CCMpfaOInteractionVolumeIndexSet(const NodalIndexSet& nodalIndexSet)
- : nodalIndexSet_( static_cast<const DualGridNodalIndexSet&>(nodalIndexSet) )
- {
- // make sure the index types used are the same in order to avoid any losses due to type conversion
- static_assert(std::is_same<GridIndexType, typename NodalIndexSet::GridIndexType>::value,
- "Provided nodal index set does not use the same type for grid indices as the given template argument");
- static_assert(std::is_same<LocalIndexType, typename NodalIndexSet::LocalIndexType>::value,
- "Provided nodal index set does not use the same type for local indices as the given template argument");
+ // Export the types used for local/grid stencils
+ using LocalStencilType = typename DualGridNodalIndexSet::LocalStencilType;
+ using GridStencilType = typename DualGridNodalIndexSet::GridStencilType;
+ using GridScvfStencilType = typename DualGridNodalIndexSet::GridScvfStencilType;
+ //! Export the type used for the neighbor scv index sets of the scvfs
+ using ScvfNeighborLocalIndexSet = typename DualGridNodalIndexSet::ScvfNeighborLocalIndexSet;
+
+ //! The constructor
+ CCMpfaOInteractionVolumeIndexSet(const NodalIndexSet& nodalIndexSet) : nodalIndexSet_(nodalIndexSet)
+ {
// determine the number of iv-local faces for memory reservation
// note that this might be a vast overestimation on surface grids!
const auto numNodalScvfs = nodalIndexSet.numScvfs();
- const auto numBoundaryScvfs = nodalIndexSet.numBoundaryScvfs();
- const std::size_t numFaceEstimate = numBoundaryScvfs + (numNodalScvfs-numBoundaryScvfs)/2;
-
- // make sure we found a reasonable number of faces
- assert((numNodalScvfs-numBoundaryScvfs)%2 == 0);
-
- // index transformation from interaction-volume-local to node-local
- ivToNodeScvf_.reserve(numFaceEstimate);
- nodeToIvScvf_.resize(numNodalScvfs);
-
- // the local neighboring scv indices of the faces
- scvfNeighborScvLocalIndices_.reserve(numFaceEstimate);
// keeps track of which nodal scvfs have been handled already
std::vector<bool> isHandled(numNodalScvfs, false);
@@ -172,7 +155,7 @@ public:
{ return nodeToIvScvf_[ nodalIndexSet_.scvfIdxLocal(scvIdxLocal, i) ]; }
//! returns the local indices of the neighboring scvs of an scvf
- const LocalIndexContainer& neighboringLocalScvIndices(LocalIndexType ivLocalScvfIdx) const
+ const ScvfNeighborLocalIndexSet& neighboringLocalScvIndices(LocalIndexType ivLocalScvfIdx) const
{ return scvfNeighborScvLocalIndices_[ivLocalScvfIdx]; }
//! returns the number of faces in the interaction volume
@@ -182,10 +165,10 @@ public:
std::size_t numScvs() const { return nodalIndexSet_.numScvs(); }
//! returns the global scv indices connected to this dual grid node
- const GridIndexContainer& globalScvIndices() const { return nodalIndexSet_.globalScvIndices(); }
+ const GridStencilType& globalScvIndices() const { return nodalIndexSet_.globalScvIndices(); }
//! returns the global scvf indices connected to this dual grid node
- const GridIndexContainer& globalScvfIndices() const { return nodalIndexSet_.globalScvfIndices(); }
+ const GridScvfStencilType& globalScvfIndices() const { return nodalIndexSet_.globalScvfIndices(); }
private:
//! returns the local scv index to a given global scv index
@@ -199,11 +182,12 @@ private:
const DualGridNodalIndexSet& nodalIndexSet_;
std::size_t numFaces_;
- std::vector<LocalIndexType> ivToNodeScvf_;
- std::vector<LocalIndexType> nodeToIvScvf_;
+ Dune::ReservedVector< LocalIndexType, NodalIndexSet::maxNumScvfsAtNode > ivToNodeScvf_;
+ Dune::ReservedVector< LocalIndexType, NodalIndexSet::maxNumScvfsAtNode > nodeToIvScvf_;
+
// maps to each scvf a list of neighbouring scv indices
// ordering: 0 - inside scv idx; 1..n - outside scv indices
- std::vector< LocalIndexContainer > scvfNeighborScvLocalIndices_;
+ Dune::ReservedVector< ScvfNeighborLocalIndexSet, NodalIndexSet::maxNumScvfsAtNode > scvfNeighborScvLocalIndices_;
};
} // end namespace Dumux
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
index 4b5f11aa27cf41a64c3a2cb8b0124f6db21bfb4d..fd6700c7273d300f4433326e262320f653e77614 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localassembler.hh
@@ -627,7 +627,7 @@ private:
* \param getTensor Lambda to evaluate the scv-wise tensors
*/
template< class IV, class GetTensor >
- void assembleLocalMatrices_(Matrix& A, Matrix& B, Matrix& C, Matrix& D,
+ void assembleLocalMatrices_(Matrix& A, Matrix& B, Matrix& C, Matrix& D,
IV& iv, const GetTensor& getTensor)
{
// Matrix D is assumed to have the right size already
@@ -672,7 +672,7 @@ private:
// we require the matrices A,B,C to have the correct size already
assert(A.rows() == iv.numUnknowns() && A.cols() == iv.numUnknowns() && "Matrix A does not have the correct size");
assert(B.rows() == iv.numUnknowns() && B.cols() == iv.numKnowns() && "Matrix B does not have the correct size");
- assert(C.rows() == iv.numFaces() && C.cols() == iv.numKnowns() && "Matrix C does not have the correct size");
+ assert(C.rows() == iv.numFaces() && C.cols() == iv.numUnknowns() && "Matrix C does not have the correct size");
// reset matrices
A = 0.0;
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
index 7ddded18c3abba12cc578707e0358606da02597d..510ace43c52a0393a4acfa1003eb7d3482409b6a 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/localsubcontrolentities.hh
@@ -134,7 +134,7 @@ private:
template< class IvIndexSet >
struct CCMpfaOInteractionVolumeLocalScvf
{
- using LocalIndexContainer = typename IvIndexSet::LocalIndexContainer;
+ using ScvfNeighborLocalIndexSet = typename IvIndexSet::ScvfNeighborLocalIndexSet;
public:
// export index types
@@ -154,7 +154,7 @@ public:
*/
template< class SubControlVolumeFace >
CCMpfaOInteractionVolumeLocalScvf(const SubControlVolumeFace& scvf,
- const LocalIndexContainer& localScvIndices,
+ const ScvfNeighborLocalIndexSet& localScvIndices,
const LocalIndexType localDofIdx,
const bool isDirichlet)
: isDirichlet_(isDirichlet)
@@ -171,7 +171,7 @@ public:
GridIndexType globalScvfIndex() const { return scvfIdxGlobal_; }
//! Returns the local indices of the scvs neighboring this scvf
- const LocalIndexContainer& neighboringLocalScvIndices() const { return *neighborScvIndicesLocal_; }
+ const ScvfNeighborLocalIndexSet& neighboringLocalScvIndices() const { return *neighborScvIndicesLocal_; }
//! states if this is scvf is on a Dirichlet boundary
bool isDirichlet() const { return isDirichlet_; }
@@ -180,7 +180,7 @@ private:
bool isDirichlet_;
GridIndexType scvfIdxGlobal_;
LocalIndexType localDofIndex_;
- const LocalIndexContainer* neighborScvIndicesLocal_;
+ const ScvfNeighborLocalIndexSet* neighborScvIndicesLocal_;
};
} // end namespace Dumux
diff --git a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
index d4401e178678120e2d3c807aa683281c052f6a35..a800e18829a237495894504e040109f07ea70073 100644
--- a/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
+++ b/dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh
@@ -52,7 +52,8 @@ template< class TypeTag, int localSize >
struct CCMpfaODefaultStaticInteractionVolumeTraits
{
private:
- using GridIndexType = typename GET_PROP_TYPE(TypeTag, GridView)::IndexSet::IndexType;
+ using NodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
+ using GridIndexType = typename NodalIndexSet::GridIndexType;
static constexpr int dim = GET_PROP_TYPE(TypeTag, GridView)::dimension;
static constexpr int dimWorld = GET_PROP_TYPE(TypeTag, GridView)::dimensionworld;
@@ -70,9 +71,9 @@ public:
//! export the type of the mpfa helper class
using MpfaHelper = typename GET_PROP_TYPE(TypeTag, MpfaHelper);
//! export the type used for local indices
- using LocalIndexType = std::uint8_t;
+ using LocalIndexType = typename NodalIndexSet::LocalIndexType;
//! export the type for the interaction volume index set
- using IndexSet = CCMpfaOInteractionVolumeIndexSet< DualGridNodalIndexSet<GridIndexType, LocalIndexType, dim> >;
+ using IndexSet = CCMpfaOInteractionVolumeIndexSet< NodalIndexSet >;
//! export the type of interaction-volume local scvs
using LocalScvType = CCMpfaOInteractionVolumeLocalScv< IndexSet, ScalarType, dim, dimWorld >;
//! export the type of interaction-volume local scvfs
@@ -83,8 +84,6 @@ public:
using Matrix = Dune::FieldMatrix< ScalarType, localSize, localSize >;
//! export the type used for iv-local vectors
using Vector = Dune::FieldVector< ScalarType, localSize >;
- //! export the type used for the iv-stencils
- using Stencil = std::vector< GridIndexType >;
//! export the data handle type for this iv
using DataHandle = InteractionVolumeDataHandle< TypeTag, Matrix, Vector, LocalIndexType >;
};
@@ -126,15 +125,7 @@ class CCMpfaOStaticInteractionVolume : public CCMpfaInteractionVolumeBase< CCMpf
using IndexSet = typename Traits::IndexSet;
using GridIndexType = typename GridView::IndexSet::IndexType;
using LocalIndexType = typename Traits::LocalIndexType;
- using Stencil = typename Traits::Stencil;
-
- //! For the o method, the interaction volume stencil can be taken directly
- //! from the nodal index set, which always uses dynamic types to be compatible
- //! on the boundaries and unstructured grids. Thus, we have to make sure that
- //! the type set for the stencils in the traits is castable.
- static_assert( std::is_convertible<Stencil*, typename IndexSet::GridIndexContainer*>::value,
- "The o-method uses the (dynamic) nodal index set's stencil as the interaction volume stencil. "
- "Using a different type is not permissive here." );
+ using Stencil = typename IndexSet::GridStencilType;
//! Data attached to scvf touching Dirichlet boundaries.
//! For the default o-scheme, we only store the corresponding vol vars index.
diff --git a/dumux/discretization/cellcentered/mpfa/properties.hh b/dumux/discretization/cellcentered/mpfa/properties.hh
index 1a5e96d98a389f00c0976233008b749cc5f4e1d0..3a20ad88449e6cfb1d6809240a930c123f73a2a1 100644
--- a/dumux/discretization/cellcentered/mpfa/properties.hh
+++ b/dumux/discretization/cellcentered/mpfa/properties.hh
@@ -47,10 +47,10 @@
#include <dumux/discretization/cellcentered/mpfa/elementvolumevariables.hh>
#include <dumux/discretization/cellcentered/mpfa/elementfluxvariablescache.hh>
#include <dumux/discretization/cellcentered/mpfa/subcontrolvolumeface.hh>
+#include <dumux/discretization/cellcentered/mpfa/dualgridindexset.hh>
#include <dumux/discretization/cellcentered/mpfa/helper.hh>
#include <dumux/discretization/cellcentered/mpfa/omethod/interactionvolume.hh>
-#include <dumux/discretization/cellcentered/mpfa/omethod/staticinteractionvolume.hh>
namespace Dumux
{
@@ -71,6 +71,44 @@ SET_PROP(CCMpfaModel, MpfaMethod)
static const MpfaMethods value = GET_PROP_TYPE(TypeTag, PrimaryInteractionVolume)::MpfaMethod;
};
+//! Set the maximum admissible number of branches per scvf
+SET_PROP(CCMpfaModel, MaxNumNeighborsPerScvf)
+{
+private:
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+
+public:
+ // Per default, we allow for 8 neighbors on network/surface grids
+ static const std::size_t value = dim < dimWorld ? 9 : 2;
+};
+
+//! Set the index set type used on the dual grid nodes
+SET_PROP(CCMpfaModel, DualGridNodalIndexSet)
+{
+private:
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+
+ // Use the grid view's index type for grid indices
+ using GI = typename GridView::IndexSet::IndexType;
+
+ // per default, use uint8_t as iv-local index type
+ using LI = std::uint8_t;
+
+ // the specified maximum admissible number of branches per scvf
+ static constexpr int maxB = GET_PROP_VALUE(TypeTag, MaxNumNeighborsPerScvf);
+
+ // maximum admissible number of elements around a node
+ // if for a given grid this number is still not high enough,
+ // overwrite this property in your problem with a higher number
+ static constexpr int dim = GridView::dimension;
+ static constexpr int maxE = dim == 3 ? 45 : 15;
+
+public:
+ using type = CCMpfaDualGridNodalIndexSet<GI, LI, dim, maxE, maxB>;
+};
+
//! The mpfa helper class
SET_TYPE_PROP(CCMpfaModel, MpfaHelper, CCMpfaHelper<TypeTag>);
diff --git a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
index 12837c09b265a23fdba0819608df3a9daac31aeb..014c54320a1529052b78b67b2c8ef9654dd42195 100644
--- a/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
+++ b/dumux/discretization/cellcentered/tpfa/fvelementgeometry.hh
@@ -61,15 +61,20 @@ class CCTpfaFVElementGeometry<TypeTag, true>
using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+public:
+ //! export type of subcontrol volume
using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ //! export type of subcontrol volume face
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using Element = typename GridView::template Codim<0>::Entity;
+ //! export type of finite volume grid geometry
using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ //! the maximum number of scvs per element
+ static constexpr std::size_t maxNumElementScvs = 1;
+ //! the maximum number of scvfs per element (use cubes for maximum)
+ static constexpr std::size_t maxNumElementScvfs = 2*GridView::dimension;
- using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>;
- using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
-
-public:
//! Constructor
CCTpfaFVElementGeometry(const FVGridGeometry& fvGridGeometry)
: fvGridGeometryPtr_(&fvGridGeometry) {}
@@ -100,9 +105,10 @@ public:
//! This is a free function found by means of ADL
//! To iterate over all sub control volumes of this FVElementGeometry use
//! for (auto&& scv : scvs(fvGeometry))
- friend inline Dune::IteratorRange<ScvIterator>
+ friend inline Dune::IteratorRange< ScvIterator<SubControlVolume, std::array<IndexType, 1>, ThisType> >
scvs(const CCTpfaFVElementGeometry& fvGeometry)
{
+ using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::array<IndexType, 1>, ThisType>;
return Dune::IteratorRange<ScvIterator>(ScvIterator(fvGeometry.scvIndices_.begin(), fvGeometry),
ScvIterator(fvGeometry.scvIndices_.end(), fvGeometry));
}
@@ -112,11 +118,12 @@ public:
//! This is a free function found by means of ADL
//! To iterate over all sub control volume faces of this FVElementGeometry use
//! for (auto&& scvf : scvfs(fvGeometry))
- friend inline Dune::IteratorRange<ScvfIterator>
+ friend inline Dune::IteratorRange< ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType> >
scvfs(const CCTpfaFVElementGeometry& fvGeometry)
{
const auto& g = fvGeometry.fvGridGeometry();
const auto scvIdx = fvGeometry.scvIndices_[0];
+ using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
return Dune::IteratorRange<ScvfIterator>(ScvfIterator(g.scvfIndicesOfScv(scvIdx).begin(), fvGeometry),
ScvfIterator(g.scvfIndicesOfScv(scvIdx).end(), fvGeometry));
}
@@ -143,7 +150,7 @@ public:
void bindElement(const Element& element)
{
elementPtr_ = &element;
- scvIndices_ = std::vector<IndexType>({fvGridGeometry().elementMapper().index(*elementPtr_)});
+ scvIndices_[0] = fvGridGeometry().elementMapper().index(*elementPtr_);
}
//! The global finite volume geometry we are a restriction of
@@ -153,7 +160,7 @@ public:
private:
const Element* elementPtr_;
- std::vector<IndexType> scvIndices_;
+ std::array<IndexType, 1> scvIndices_;
const FVGridGeometry* fvGridGeometryPtr_;
};
@@ -168,18 +175,23 @@ class CCTpfaFVElementGeometry<TypeTag, false>
using ThisType = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using IndexType = typename GridView::IndexSet::IndexType;
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using Element = typename GridView::template Codim<0>::Entity;
- using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
static const int dim = GridView::dimension;
static const int dimWorld = GridView::dimensionworld;
- using ScvIterator = Dumux::ScvIterator<SubControlVolume, std::vector<IndexType>, ThisType>;
- using ScvfIterator = Dumux::ScvfIterator<SubControlVolumeFace, std::vector<IndexType>, ThisType>;
-
public:
+ //! export type of subcontrol volume
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ //! export type of subcontrol volume face
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ //! export type of finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ //! the maximum number of scvs per element
+ static constexpr std::size_t maxNumElementScvs = 1;
+ //! the maximum number of scvfs per element (use cubes for maximum)
+ static constexpr std::size_t maxNumElementScvfs = 2*dim;
+
//! Constructor
CCTpfaFVElementGeometry(const FVGridGeometry& fvGridGeometry)
: fvGridGeometryPtr_(&fvGridGeometry) {}
@@ -231,11 +243,11 @@ public:
//! This is a free function found by means of ADL
//! To iterate over all sub control volumes of this FVElementGeometry use
//! for (auto&& scv : scvs(fvGeometry))
- friend inline Dune::IteratorRange<typename std::vector<SubControlVolume>::const_iterator>
+ friend inline Dune::IteratorRange<typename std::array<SubControlVolume, 1>::const_iterator>
scvs(const ThisType& g)
{
- using Iter = typename std::vector<SubControlVolume>::const_iterator;
- return Dune::IteratorRange<Iter>(g.scvs_.begin(), g.scvs_.end());
+ using IteratorType = typename std::array<SubControlVolume, 1>::const_iterator;
+ return Dune::IteratorRange<IteratorType>(g.scvs_.begin(), g.scvs_.end());
}
//! iterator range for sub control volumes faces. Iterates over
@@ -246,8 +258,8 @@ public:
friend inline Dune::IteratorRange<typename std::vector<SubControlVolumeFace>::const_iterator>
scvfs(const ThisType& g)
{
- using Iter = typename std::vector<SubControlVolumeFace>::const_iterator;
- return Dune::IteratorRange<Iter>(g.scvfs_.begin(), g.scvfs_.end());
+ using IteratorType = typename std::vector<SubControlVolumeFace>::const_iterator;
+ return Dune::IteratorRange<IteratorType>(g.scvfs_.begin(), g.scvfs_.end());
}
//! number of sub control volumes in this fv element geometry
@@ -384,8 +396,8 @@ private:
void makeElementGeometries(const Element& element)
{
const auto eIdx = fvGridGeometry().elementMapper().index(element);
- scvs_.emplace_back(element.geometry(), eIdx);
- scvIndices_.emplace_back(eIdx);
+ scvs_[0] = SubControlVolume(element.geometry(), eIdx);
+ scvIndices_[0] = eIdx;
const auto& scvFaceIndices = fvGridGeometry().scvfIndicesOfScv(eIdx);
const auto& neighborVolVarIndices = fvGridGeometry().neighborVolVarIndices(eIdx);
@@ -523,9 +535,7 @@ private:
//! Clear all local data
void clear()
{
- scvIndices_.clear();
scvfIndices_.clear();
- scvs_.clear();
scvfs_.clear();
flippedScvfIndices_.clear();
@@ -536,21 +546,21 @@ private:
flippedNeighborScvfIndices_.clear();
}
- //! the bound element
- const Element* elementPtr_;
-
- const FVGridGeometry* fvGridGeometryPtr_;
+ const Element* elementPtr_; //!< the element to which this fvgeometry is bound
+ const FVGridGeometry* fvGridGeometryPtr_; //!< the grid fvgeometry
// local storage after binding an element
- std::vector<IndexType> scvIndices_;
+ std::array<IndexType, 1> scvIndices_;
+ std::array<SubControlVolume, 1> scvs_;
+
std::vector<IndexType> scvfIndices_;
- std::vector<SubControlVolume> scvs_;
std::vector<SubControlVolumeFace> scvfs_;
std::vector<std::vector<IndexType>> flippedScvfIndices_;
std::vector<IndexType> neighborScvIndices_;
- std::vector<IndexType> neighborScvfIndices_;
std::vector<SubControlVolume> neighborScvs_;
+
+ std::vector<IndexType> neighborScvfIndices_;
std::vector<SubControlVolumeFace> neighborScvfs_;
std::vector<std::vector<IndexType>> flippedNeighborScvfIndices_;
};
diff --git a/dumux/discretization/cellcentered/tpfa/properties.hh b/dumux/discretization/cellcentered/tpfa/properties.hh
index 8e2bca2a4857885649ade46a74e54528762004b1..3710cf0fa7f339df01da81fc79ae8b570dac3f81 100644
--- a/dumux/discretization/cellcentered/tpfa/properties.hh
+++ b/dumux/discretization/cellcentered/tpfa/properties.hh
@@ -82,6 +82,19 @@ SET_TYPE_PROP(CCTpfaModel, ElementFluxVariablesCache, CCTpfaElementFluxVariables
//! The global current volume variables vector class
SET_TYPE_PROP(CCTpfaModel, GridVolumeVariables, CCGridVolumeVariables<TypeTag, GET_PROP_VALUE(TypeTag, EnableGridVolumeVariablesCache)>);
+//! Set the maximum admissible number of branches per scvf
+SET_PROP(CCTpfaModel, MaxNumNeighborsPerScvf)
+{
+private:
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+
+public:
+ // Per default, we allow for 8 neighbors on network/surface grids
+ static const std::size_t value = dim < dimWorld ? 9 : 2;
+};
+
//! The sub control volume
SET_PROP(CCTpfaModel, SubControlVolume)
{
diff --git a/dumux/discretization/fluxstencil.hh b/dumux/discretization/fluxstencil.hh
index 0f24ee01654c42b9d7f7e074f9c0915a65b96ab1..a138efb0a86c1d47f8812807a75e088fcde77d3e 100644
--- a/dumux/discretization/fluxstencil.hh
+++ b/dumux/discretization/fluxstencil.hh
@@ -24,6 +24,7 @@
#ifndef DUMUX_DISCRETIZATION_FLUXSTENCIL_HH
#define DUMUX_DISCRETIZATION_FLUXSTENCIL_HH
+#include <dune/common/reservedvector.hh>
#include <dumux/common/properties.hh>
#include <dumux/discretization/methods.hh>
@@ -58,9 +59,17 @@ class FluxStencilImplementation<TypeTag, DiscretizationMethods::CCTpfa>
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
- using Stencil = std::vector<IndexType>;
public:
+ //! The maximum number of elements in a flux stencil
+ static constexpr int maxFluxStencilSize = GET_PROP_VALUE(TypeTag, MaxNumNeighborsPerScvf);
+
+ //! States how many scvfs of an element J might have an element I in the flux stencil
+ static constexpr int maxNumScvfJForI = 1;
+
+ //! The flux stencil type
+ using Stencil = Dune::ReservedVector<IndexType, maxFluxStencilSize>;
+
static Stencil stencil(const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
@@ -91,12 +100,36 @@ class FluxStencilImplementation<TypeTag, DiscretizationMethods::CCMpfa>
using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using Element = typename GridView::template Codim<0>::Entity;
using IndexType = typename GridView::IndexSet::IndexType;
- using Stencil = std::vector<IndexType>;
+ static constexpr int dim = GridView::dimension;
+
+ // Use the stencil type of the primary interaction volume
+ using NodalIndexSet = typename GET_PROP_TYPE(TypeTag, DualGridNodalIndexSet);
public:
- static Stencil stencil(const Element& element,
- const FVElementGeometry& fvGeometry,
- const SubControlVolumeFace& scvf)
+ //! The maximum number of elements in a flux stencil (equal to number of elements at node)
+ static constexpr int maxFluxStencilSize = NodalIndexSet::maxNumElementsAtNode;
+
+ //! States how many scvfs of an element J might have an element I in the flux stencil
+ //! We use cubes here for determining the maximum. Only basic geometries are not supported.
+ static constexpr int maxNumScvfJForI = dim == 3 ? 12 : 4;
+
+ //! The flux stencil type
+ using Stencil = typename NodalIndexSet::GridStencilType;
+
+ /*
+ * \brief Returns a set of grid element indices that participate in the
+ * flux calculations on a given scvf.
+ *
+ * \note The interaction volume index sets must use the same type for the
+ * stencils as the nodal index set. If not, the compiler will complain here.
+ *
+ * \param element The grid element
+ * \param fvGeometry The finite volume geometry of this element
+ * \param scvf The sub-control volume face embedded in this element
+ */
+ static const Stencil& stencil(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const SubControlVolumeFace& scvf)
{
const auto& fvGridGeometry = fvGeometry.fvGridGeometry();
diff --git a/dumux/discretization/upwindscheme.hh b/dumux/discretization/upwindscheme.hh
index 38eb68b2299b36883176f3e61143c08bb56ab213..e2056f47faf943f2889144b35f73fe47c2c2066a 100644
--- a/dumux/discretization/upwindscheme.hh
+++ b/dumux/discretization/upwindscheme.hh
@@ -56,14 +56,15 @@ public:
const UpwindTermFunction& upwindTerm,
Scalar flux, int phaseIdx)
{
- // retrieve the upwind weight for the mass conservation equations. Use the value
- // specified via the property system as default, and overwrite
- // it by the run-time parameter from the Dune::ParameterTree
- static const Scalar upwindWeight = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
-
- const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()];
- const auto& outsideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().outsideScvIdx()];
- if (std::signbit(flux))
+ static const Scalar upwindWeight = getParamFromGroup<Scalar>
+ (GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
+
+ const auto& scvf = fluxVars.scvFace();
+ const auto& elemVolVars = fluxVars.elemVolVars();
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+
+ if (std::signbit(flux)) // if sign of flux is negative
return flux*(upwindWeight*upwindTerm(outsideVolVars)
+ (1.0 - upwindWeight)*upwindTerm(insideVolVars));
else
@@ -91,17 +92,16 @@ public:
const UpwindTermFunction& upwindTerm,
Scalar flux, int phaseIdx)
{
- // retrieve the upwind weight for the mass conservation equations. Use the value
- // specified via the property system as default, and overwrite
- // it by the run-time parameter from the Dune::ParameterTree
- static const std::string modelParamGroup = GET_PROP_VALUE(TypeTag, ModelParameterGroup);
- static const Scalar upwindWeight = getParamFromGroup<Scalar>(modelParamGroup, "Implicit.UpwindWeight");
+ static const Scalar upwindWeight = getParamFromGroup<Scalar>
+ (GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
// the volume variables of the inside sub-control volume
- const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()];
+ const auto& scvf = fluxVars.scvFace();
+ const auto& elemVolVars = fluxVars.elemVolVars();
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
// check if this is a branching point
- if (fluxVars.scvFace().numOutsideScvs() > 1)
+ if (scvf.numOutsideScvs() > 1)
{
// more complicated upwind scheme
// we compute a flux-weighted average of all inflowing branches
@@ -114,23 +114,24 @@ public:
else
sumUpwindFluxes += flux;
- for (unsigned int i = 0; i < fluxVars.scvFace().numOutsideScvs(); ++i)
+ for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i)
{
- // compute the outside flux
- const auto outsideScvIdx = fluxVars.scvFace().outsideScvIdx(i);
- const auto outsideElement = fluxVars.fvGeometry().fvGridGeometry().element(outsideScvIdx);
- const auto& flippedScvf = fluxVars.fvGeometry().flipScvf(fluxVars.scvFace().index(), i);
+ // compute the outside flux
+ const auto& fvGeometry = fluxVars.fvGeometry();
+ const auto outsideScvIdx = scvf.outsideScvIdx(i);
+ const auto outsideElement = fvGeometry.fvGridGeometry().element(outsideScvIdx);
+ const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i);
const auto outsideFlux = AdvectionType::flux(fluxVars.problem(),
outsideElement,
- fluxVars.fvGeometry(),
- fluxVars.elemVolVars(),
+ fvGeometry,
+ elemVolVars,
flippedScvf,
phaseIdx,
fluxVars.elemFluxVarsCache());
if (!std::signbit(outsideFlux))
- branchingPointUpwindTerm += upwindTerm(fluxVars.elemVolVars()[outsideScvIdx])*outsideFlux;
+ branchingPointUpwindTerm += upwindTerm(elemVolVars[outsideScvIdx])*outsideFlux;
else
sumUpwindFluxes += outsideFlux;
}
@@ -153,7 +154,7 @@ public:
else
{
// upwind scheme
- const auto& outsideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().outsideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
if (std::signbit(flux))
return flux*(upwindWeight*upwindTerm(outsideVolVars)
+ (1.0 - upwindWeight)*upwindTerm(insideVolVars));
@@ -170,14 +171,15 @@ public:
const UpwindTermFunction& upwindTerm,
Scalar flux, int phaseIdx)
{
- // retrieve the upwind weight for the mass conservation equations. Use the value
- // specified via the property system as default, and overwrite
- // it by the run-time parameter from the Dune::ParameterTree
- static const Scalar upwindWeight = getParamFromGroup<Scalar>(GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
-
- const auto& insideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().insideScvIdx()];
- const auto& outsideVolVars = fluxVars.elemVolVars()[fluxVars.scvFace().outsideScvIdx()];
- if (std::signbit(flux))
+ static const Scalar upwindWeight = getParamFromGroup<Scalar>
+ (GET_PROP_VALUE(TypeTag, ModelParameterGroup), "Implicit.UpwindWeight");
+
+ const auto& scvf = fluxVars.scvFace();
+ const auto& elemVolVars = fluxVars.elemVolVars();
+ const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
+ const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
+
+ if (std::signbit(flux)) // if sign of flux is negative
return flux*(upwindWeight*upwindTerm(outsideVolVars)
+ (1.0 - upwindWeight)*upwindTerm(insideVolVars));
else
diff --git a/dumux/nonlinear/newtoncontroller.hh b/dumux/nonlinear/newtoncontroller.hh
index 39ff8a2ebb8a0a4349f02fad3ee217da4b192cd7..4a5513fb17152b067decd4ca522b97943fbe59d8 100644
--- a/dumux/nonlinear/newtoncontroller.hh
+++ b/dumux/nonlinear/newtoncontroller.hh
@@ -456,7 +456,7 @@ public:
template<class Assembler, class SolutionVector>
void newtonFail(Assembler& assembler, SolutionVector& u)
{
- if (!assembler.localResidual().isStationary())
+ if (!assembler.isStationaryProblem())
{
// set solution to previous solution
u = assembler.prevSol();
diff --git a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
index 9e0a78bb7a16aa579d48d212d11011b094491a2f..dfdc173f082d9f9aced2fd6b416775cca63f27c4 100644
--- a/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
+++ b/dumux/porousmediumflow/1p/incompressiblelocalresidual.hh
@@ -124,9 +124,9 @@ public:
/ curElemVolVars[scvf.insideScvIdx()].viscosity();
const auto& fluxVarsCache = elemFluxVarsCache[scvf];
+ const auto& stencil = fluxVarsCache.advectionStencil();
if (fluxVarsCache.usesSecondaryIv())
{
- const auto& stencil = fluxVarsCache.advectionStencilSecondaryIv();
const auto& tij = fluxVarsCache.advectionTijSecondaryIv();
assert(stencil.size() == tij.size());
@@ -141,7 +141,6 @@ public:
}
else
{
- const auto& stencil = fluxVarsCache.advectionStencilPrimaryIv();
const auto& tij = fluxVarsCache.advectionTijPrimaryIv();
assert(stencil.size() == tij.size());
diff --git a/dumux/porousmediumflow/mpnc/localresidual.hh b/dumux/porousmediumflow/mpnc/localresidual.hh
index 5bf07dc6ea32bc0eac0a81c242a6525ce41e5219..d982de82f85d55f6c8ae25fa820ff196f0967036 100644
--- a/dumux/porousmediumflow/mpnc/localresidual.hh
+++ b/dumux/porousmediumflow/mpnc/localresidual.hh
@@ -45,16 +45,10 @@ template<class TypeTag>
class MPNCLocalResidual : public CompositionalLocalResidual<TypeTag>
{
using ParentType = CompositionalLocalResidual<TypeTag>;
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
using Element = typename GET_PROP_TYPE(TypeTag, GridView)::template Codim<0>::Entity;
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
- using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
- using ElementResidualVector = Dune::BlockVector<typename GET_PROP_TYPE(TypeTag, NumEqVector)>;
using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
@@ -64,102 +58,22 @@ class MPNCLocalResidual : public CompositionalLocalResidual<TypeTag>
public:
using ParentType::ParentType;
+ using typename ParentType::ElementResidualVector;
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero for instationary problems.
- * \param problem The problem to solve
-
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param fvGeometry The finite-volume geometry of the element
- * \param prevVolVars The volume averaged variables for all
- * sub-control volumes of the element at the previous
- * time level
- * \param curVolVars The volume averaged variables for all
- * sub-control volumes of the element at the current
- * time level
- * \param bcTypes The types of the boundary conditions for all
- * vertices of the element
- */
- ElementResidualVector eval(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& prevElemVolVars,
- const ElementVolumeVariables& curElemVolVars,
- const ElementBoundaryTypes &bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache) const
- {
- // initialize the residual vector for all scvs in this element
- ElementResidualVector residual(fvGeometry.numScv());
- residual = 0.0;
-
- // evaluate the volume terms (storage + source terms)
- for (auto&& scv : scvs(fvGeometry))
- {
- //! foward to the local residual specialized for the discretization methods
- this->asImp().evalStorage(residual, problem, element, fvGeometry, prevElemVolVars, curElemVolVars, scv);
- this->asImp().evalSource(residual, problem, element, fvGeometry, curElemVolVars, scv);
- //here we need to set the constraints of the mpnc model into the residual
- const auto localScvIdx = scv.indexInElement();
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- residual[localScvIdx][phase0NcpIdx + phaseIdx] = curElemVolVars[scv].phaseNcp(phaseIdx);
- }
- }
-
- for (auto&& scvf : scvfs(fvGeometry))
- {
- //! foward to the local residual specialized for the discretization methods
- this->asImp().evalFlux(residual, problem, element, fvGeometry, curElemVolVars, bcTypes, elemFluxVarsCache, scvf);
- }
-
- return residual;
- }
-
- /*!
- * \brief Compute the local residual, i.e. the deviation of the
- * equations from zero for stationary problem.
- * \param problem The problem to solve
-
- * \param element The DUNE Codim<0> entity for which the residual
- * ought to be calculated
- * \param fvGeometry The finite-volume geometry of the element
- * \param curVolVars The volume averaged variables for all
- * sub-control volumes of the element at the current
- * time level
- * \param bcTypes The types of the boundary conditions for all
- * vertices of the element
- */
- ElementResidualVector eval(const Problem& problem,
- const Element& element,
- const FVElementGeometry& fvGeometry,
- const ElementVolumeVariables& curElemVolVars,
- const ElementBoundaryTypes &bcTypes,
- const ElementFluxVariablesCache& elemFluxVarsCache) const
+ ElementResidualVector evalFluxSource(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes &bcTypes) const
{
- // initialize the residual vector for all scvs in this element
- ElementResidualVector residual(fvGeometry.numScv());
- residual = 0.0;
+ ElementResidualVector residual = ParentType::evalFluxSource(element, fvGeometry, elemVolVars, elemFluxVarsCache, bcTypes);
- // evaluate the volume terms (storage + source terms)
for (auto&& scv : scvs(fvGeometry))
{
- //! foward to the local residual specialized for the discretization methods
- this->asImp().evalSource(residual, problem, element, fvGeometry, curElemVolVars, scv);
-
- //here we need to set the constraints of the mpnc model into the residual
- const auto localScvIdx = scv.indexInElement();
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- residual[localScvIdx][phase0NcpIdx + phaseIdx] = curElemVolVars[scv].phaseNcp(phaseIdx);
- }
- }
-
- for (auto&& scvf : scvfs(fvGeometry))
- {
- //! foward to the local residual specialized for the discretization methods
- this->asImp().evalFlux(residual, problem, element, fvGeometry, curElemVolVars, bcTypes, elemFluxVarsCache, scvf);
+ //here we need to set the constraints of the mpnc model into the residual
+ const auto localScvIdx = scv.indexInElement();
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ residual[localScvIdx][phase0NcpIdx + phaseIdx] = elemVolVars[scv].phaseNcp(phaseIdx);
}
return residual;