diff --git a/dumux/mixeddimension/assembler.hh b/dumux/mixeddimension/assembler.hh
index 8c25ab1cbfc203d5a73ca6a5a004afb26e08cc04..b66efd7f4fb983ef36a1fff23473d4ea02aaaeff 100644
--- a/dumux/mixeddimension/assembler.hh
+++ b/dumux/mixeddimension/assembler.hh
@@ -182,7 +182,8 @@ protected:
residual_[lowDimIdx]);
}
- void buildBulkMatrixBlocksCC_(BulkMatrixBlock& m, BulkCouplingMatrixBlock& cm)
+ template <class T = BulkProblemTypeTag, typename std::enable_if<((!GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
+ void buildBulkMatrixBlocks_(BulkMatrixBlock& m, BulkCouplingMatrixBlock& cm)
{
// get occupation pattern of the matrix
Dune::MatrixIndexSet bulkPattern, bulkCouplingPattern;
@@ -213,7 +214,8 @@ protected:
bulkCouplingPattern.exportIdx(cm);
}
- void buildLowDimMatrixBlocksCC_(LowDimMatrixBlock& m, LowDimCouplingMatrixBlock& cm)
+ template <class T = LowDimProblemTypeTag, typename std::enable_if<((!GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
+ void buildLowDimMatrixBlocks_(LowDimMatrixBlock& m, LowDimCouplingMatrixBlock& cm)
{
// get occupation pattern of the matrix
Dune::MatrixIndexSet lowDimPattern, lowDimCouplingPattern;
@@ -244,63 +246,73 @@ protected:
lowDimCouplingPattern.exportIdx(cm);
}
- // void buildBulkMatrixBlocksBox_(BulkMatrixBlock& m, BulkCouplingMatrixBlock& cm)
- // {
- // // get occupation pattern of the matrix
- // Dune::MatrixIndexSet bulkPattern, bulkCouplingPattern;
- // bulkPattern.resize(m.N(), m.M());
- // bulkCouplingPattern.resize(cm.N(), cm.M());
-
- // for (const auto& element : elements(problem_().bulkGridView()))
- // {
- // const auto& stencil = problem_().couplingManager().stencil(element);
- // const auto& couplingStencil = problem_().couplingManager().couplingStencil(element);
-
- // for (unsigned int scvIdx = 0; scvIdx < element.subEntities(bulkDim); ++scvIdx)
- // {
- // auto globalI = problem_().model().bulkVertexMapper().subIndex(element, scvIdx, bulkDim);
-
- // for (auto&& globalJ : stencil)
- // bulkPattern.add(globalI, globalJ);
-
- // for (auto&& globalJ : couplingStencil)
- // bulkCouplingPattern.add(globalI, globalJ);
- // }
- // }
+ template <class T = BulkProblemTypeTag, typename std::enable_if<((GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
+ void buildBulkMatrixBlocks_(BulkMatrixBlock& m, BulkCouplingMatrixBlock& cm)
+ {
+ // get occupation pattern of the matrix
+ Dune::MatrixIndexSet bulkPattern, bulkCouplingPattern;
+ bulkPattern.resize(m.N(), m.M());
+ bulkCouplingPattern.resize(cm.N(), cm.M());
- // // export occupation patterns to matrices
- // bulkPattern.exportIdx(m);
- // bulkCouplingPattern.exportIdx(cm);
- // }
+ for (const auto& element : elements(problem_().bulkGridView()))
+ {
+ const auto& couplingStencil = problem_().couplingManager().couplingStencil(element);
- // void buildLowDimMatrixBlocksBox_(LowDimMatrixBlock& m, LowDimCouplingMatrixBlock& cm)
- // {
- // // get occupation pattern of the matrix
- // Dune::MatrixIndexSet lowDimPattern, lowDimCouplingPattern;
- // lowDimPattern.resize(m.N(), m.M());
- // lowDimCouplingPattern.resize(cm.N(), cm.M());
+ for (unsigned int vIdx = 0; vIdx < element.subEntities(bulkDim); ++vIdx)
+ {
+ const auto globalI = problem_().model().bulkVertexMapper().subIndex(element, vIdx, bulkDim);
+ for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(bulkDim); ++vIdx2)
+ {
+ const auto globalJ = problem_().model().bulkVertexMapper().subIndex(element, vIdx2, bulkDim);
+ bulkPattern.add(globalI, globalJ);
+ bulkPattern.add(globalJ, globalI);
+ }
+
+ for (auto&& globalJ : couplingStencil)
+ bulkCouplingPattern.add(globalI, globalJ);
+
+ //TODO: additionalDofDependencies
+ }
+ }
- // for (const auto& element : elements(problem_().lowDimGridView()))
- // {
- // const auto& stencil = problem_().couplingManager().stencil(element);
- // const auto& couplingStencil = problem_().couplingManager().couplingStencil(element);
+ // export occupation patterns to matrices
+ bulkPattern.exportIdx(m);
+ bulkCouplingPattern.exportIdx(cm);
+ }
- // for (unsigned int scvIdx = 0; scvIdx < element.subEntities(lowDimDim); ++scvIdx)
- // {
- // auto globalI = problem_().model().lowDimVertexMapper().subIndex(element, scvIdx, lowDimDim);
+ template <class T = LowDimProblemTypeTag, typename std::enable_if<((GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
+ void buildLowDimMatrixBlocks_(LowDimMatrixBlock& m, LowDimCouplingMatrixBlock& cm)
+ {
+ // get occupation pattern of the matrix
+ Dune::MatrixIndexSet lowDimPattern, lowDimCouplingPattern;
+ lowDimPattern.resize(m.N(), m.M());
+ lowDimCouplingPattern.resize(cm.N(), cm.M());
- // for (auto&& globalJ : stencil)
- // lowDimPattern.add(globalI, globalJ);
+ for (const auto& element : elements(problem_().lowDimGridView()))
+ {
+ const auto& couplingStencil = problem_().couplingManager().couplingStencil(element);
- // for (auto&& globalJ : couplingStencil)
- // lowDimCouplingPattern.add(globalI, globalJ);
- // }
- // }
+ for (unsigned int vIdx = 0; vIdx < element.subEntities(lowDimDim); ++vIdx)
+ {
+ const auto globalI = problem_().model().lowDimVertexMapper().subIndex(element, vIdx, lowDimDim);
+ for (unsigned int vIdx2 = vIdx; vIdx2 < element.subEntities(lowDimDim); ++vIdx2)
+ {
+ const auto globalJ = problem_().model().lowDimVertexMapper().subIndex(element, vIdx2, lowDimDim);
+ lowDimPattern.add(globalI, globalJ);
+ lowDimPattern.add(globalJ, globalI);
+ }
+
+ for (auto&& globalJ : couplingStencil)
+ lowDimCouplingPattern.add(globalI, globalJ);
+
+ //TODO: additionalDofDependencies
+ }
+ }
- // // export occupation patterns to matrices
- // lowDimPattern.exportIdx(m);
- // lowDimCouplingPattern.exportIdx(cm);
- // }
+ // export occupation patterns to matrices
+ lowDimPattern.exportIdx(m);
+ lowDimCouplingPattern.exportIdx(cm);
+ }
// Construct the multitype matrix for the global jacobian
void createMatrix_()
@@ -318,18 +330,8 @@ protected:
auto A22 = LowDimMatrixBlock(lowDimSize, lowDimSize, LowDimMatrixBlock::random);
auto A21 = LowDimCouplingMatrixBlock(lowDimSize, bulkSize, LowDimCouplingMatrixBlock::random);
- // cell-centered
- if (!bulkIsBox)
- buildBulkMatrixBlocksCC_(A11, A12);
- else{
- // buildBulkMatrixBlocksBox_(A11, A12);
- }
-
- if (!lowDimIsBox)
- buildLowDimMatrixBlocksCC_(A22, A21);
- else{
- // buildLowDimMatrixBlocksBox_(A22, A21);
- }
+ buildBulkMatrixBlocks_(A11, A12);
+ buildLowDimMatrixBlocks_(A22, A21);
(*matrix_)[bulkIdx][bulkIdx] = A11;
(*matrix_)[bulkIdx][lowDimIdx] = A12;
diff --git a/dumux/mixeddimension/subproblemlocaljacobian.hh b/dumux/mixeddimension/subproblemlocaljacobian.hh
index 605be0a6a62fcadc2969e7b0084ca537c25720b5..43ddbb99cdf8466a1248a791eaabdbcd28672e56 100644
--- a/dumux/mixeddimension/subproblemlocaljacobian.hh
+++ b/dumux/mixeddimension/subproblemlocaljacobian.hh
@@ -104,6 +104,8 @@ class SubProblemLocalJacobian : public GET_PROP_TYPE(SubProblemTypeTag, LocalJac
using Element = typename GridView::template Codim<0>::Entity;
using SolutionVector = typename GET_PROP_TYPE(SubProblemTypeTag, SolutionVector);
+ enum { isBox = GET_PROP_VALUE(SubProblemTypeTag, ImplicitIsBox) };
+
public:
// copying a local jacobian is not a good idea
@@ -138,8 +140,6 @@ public:
JacobianMatrixCoupling& couplingMatrix,
SolutionVector& residual)
{
- const bool isGhost = (element.partitionType() == Dune::GhostEntity);
-
// prepare the volvars/fvGeometries in case caching is disabled
auto fvGeometry = localView(this->model_().globalFvGeometry());
fvGeometry.bind(element);
@@ -153,13 +153,32 @@ public:
auto elemFluxVarsCache = localView(this->model_().globalFluxVarsCache());
elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
- // set the actual dof index
- this->globalI_ = this->problem_().elementMapper().index(element);
-
// check for boundaries on the element
ElementBoundaryTypes elemBcTypes;
elemBcTypes.update(this->problem_(), element, fvGeometry);
+ assemble_(element, fvGeometry, prevElemVolVars, curElemVolVars, elemFluxVarsCache, elemBcTypes, matrix, couplingMatrix, residual);
+ }
+
+protected:
+
+ // for cell-centered models
+ template<class JacobianMatrix, class JacobianMatrixCoupling, class T = SubProblemTypeTag, typename std::enable_if<((!GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
+ void assemble_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ ElementVolumeVariables& prevElemVolVars,
+ ElementVolumeVariables& curElemVolVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ JacobianMatrixCoupling& couplingMatrix,
+ SolutionVector& residual)
+ {
+ const bool isGhost = (element.partitionType() == Dune::GhostEntity);
+
+ // set the actual dof index
+ this->globalI_ = this->problem_().elementMapper().index(element);
+
// Evaluate the undeflected element local residual
this->localResidual().eval(element,
fvGeometry,
@@ -204,7 +223,73 @@ public:
residual);
}
-protected:
+ // for box models
+ template<class JacobianMatrix, class JacobianMatrixCoupling, class T = SubProblemTypeTag, typename std::enable_if<((GET_PROP_VALUE(T, ImplicitIsBox))), bool>::type = 0>
+ void assemble_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ ElementVolumeVariables& prevElemVolVars,
+ ElementVolumeVariables& curElemVolVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrix& matrix,
+ JacobianMatrixCoupling& couplingMatrix,
+ SolutionVector& residual)
+ {
+ constexpr auto numEq = GET_PROP_VALUE(T, NumEq);
+
+ // the element solution
+ auto curElemSol = this->problem_().model().elementSolution(element, this->problem_().model().curSol());
+
+ // calculate the actual element residual
+ this->localResidual().eval(element, fvGeometry, prevElemVolVars, curElemVolVars, elemBcTypes, elemFluxVarsCache);
+ this->residual_ = this->localResidual().residual();
+
+ // residual[this->globalI_] = this->localResidual().residual(0);
+
+ this->problem_().model().updatePVWeights(fvGeometry);
+
+ // calculation of the derivatives
+ for (auto&& scv : scvs(fvGeometry))
+ {
+ // dof index and corresponding actual pri vars
+ const auto dofIdx = scv.dofIndex();
+ auto& curVolVars = this->getCurVolVars(curElemVolVars, scv);
+ VolumeVariables origVolVars(curVolVars);
+
+ // add precalculated residual for this scv into the global container
+ residual[dofIdx] += this->residual_[scv.index()];
+
+ // calculate derivatives w.r.t to the privars at the dof at hand
+ for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ {
+ this->evalPartialDerivative_(matrix,
+ element,
+ fvGeometry,
+ prevElemVolVars,
+ curElemVolVars,
+ curElemSol,
+ scv,
+ elemBcTypes,
+ elemFluxVarsCache,
+ pvIdx);
+
+ // restore the original state of the scv's volume variables
+ curVolVars = origVolVars;
+
+ // restore the original element solution
+ curElemSol[scv.index()][pvIdx] = this->problem_().model().curSol()[scv.dofIndex()][pvIdx];
+ }
+ // TODO: what if we have an extended source stencil????
+ }
+ //TODO: is otherElement in this method required? is otherResidual correct?
+ evalPartialDerivativeCouplingBox_(element,
+ fvGeometry,
+ curElemVolVars,
+ elemFluxVarsCache,
+ elemBcTypes,
+ couplingMatrix,
+ residual);
+ }
/*!
* \brief Returns a reference to the problem.
@@ -234,6 +319,7 @@ protected:
decltype(auto) otherProblem_(typename std::enable_if<!std::is_same<typename GET_PROP_TYPE(T, LowDimProblemTypeTag), SubProblemTypeTag>::value, void>::type* x = nullptr)
{ return globalProblem_().lowDimProblem(); }
+ // cell-centered
template<class JacobianMatrixCoupling>
void evalPartialDerivativeCoupling_(const Element& element,
const FVElementGeometry& fvGeometry,
@@ -329,6 +415,102 @@ protected:
}
}
+ // box
+ template<class JacobianMatrixCoupling>
+ void evalPartialDerivativeCouplingBox_(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ ElementVolumeVariables& curElemVolVars,
+ ElementFluxVariablesCache& elemFluxVarsCache,
+ const ElementBoundaryTypes& elemBcTypes,
+ JacobianMatrixCoupling& couplingMatrix,
+ SolutionVector& residual)
+ {
+ const auto& couplingStencil = globalProblem_().couplingManager().couplingStencil(element);
+ constexpr auto numEq = GET_PROP_VALUE(SubProblemTypeTag, NumEq);
+
+ for (auto globalJ : couplingStencil)
+ {
+ const auto otherResidual = globalProblem_().couplingManager().evalCouplingResidual(element,
+ fvGeometry,
+ curElemVolVars,
+ elemBcTypes,
+ elemFluxVarsCache);
+
+ auto& otherPriVars = otherProblem_().model().curSol()[globalJ];
+ auto originalOtherPriVars = otherPriVars;
+
+ // derivatives in the neighbors with repect to the current elements
+ ElementSolutionVector partialDeriv(fvGeometry.numScv());
+ for (int pvIdx = 0; pvIdx < numEq; pvIdx++)
+ {
+ const Scalar eps = this->numericEpsilon(otherPriVars[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
+ otherPriVars[pvIdx] += eps;
+ delta += eps;
+
+ // calculate the residual with the deflected primary variables
+ partialDeriv = globalProblem_().couplingManager().evalCouplingResidual(element,
+ fvGeometry,
+ 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 = otherResidual;
+ }
+
+
+ if (numericDifferenceMethod_ <= 0)
+ {
+ // we are not using forward differences, i.e. we
+ // need to calculate f(x - \epsilon)
+
+ // deflect the primary variables
+ otherPriVars[pvIdx] -= 2*eps;
+ delta += eps;
+
+ // calculate the residual with the deflected primary variables
+ partialDeriv -= globalProblem_().couplingManager().evalCouplingResidual(element,
+ fvGeometry,
+ 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 -= otherResidual;
+ }
+
+ // divide difference in residuals by the magnitude of the
+ // deflections between the two function evaluation
+ partialDeriv /= delta;
+
+ // restore the original state of the element solution vector
+ otherPriVars = originalOtherPriVars;
+
+ // update the global jacobian matrix (coupling block)
+ // this->updateGlobalJacobian_(couplingMatrix, scv.dofIndex(), globalJ, pvIdx, partialDeriv[scv.dofIndex()]);
+
+ for (auto&& scv : scvs(fvGeometry))
+ this->updateGlobalJacobian_(couplingMatrix, scv.dofIndex(), globalJ, pvIdx, partialDeriv[scv.index()]);
+ }
+ }
+ }
+
// The problem we would like to solve
GlobalProblem *globalProblemPtr_;
// The type of the numeric difference method (forward, center, backward)