diff --git a/dumux/assembly/jacobianpattern.hh b/dumux/assembly/jacobianpattern.hh
index 977b6d01374b78e91d4326819729f20286bb2f59..3829a66b73a91210da7c05062dc495ec91b0633b 100644
--- a/dumux/assembly/jacobianpattern.hh
+++ b/dumux/assembly/jacobianpattern.hh
@@ -125,7 +125,6 @@ auto getJacobianPattern(const GridGeometry& gridGeometry)
const auto numDofs = gridGeometry.numDofs();
Dune::MatrixIndexSet pattern(numDofs, numDofs);
-
const auto& connectivityMap = gridGeometry.connectivityMap();
// evaluate the acutal pattern
@@ -136,7 +135,7 @@ auto getJacobianPattern(const GridGeometry& gridGeometry)
// the global index of the element at hand
static constexpr auto cellCenterIdx = GridGeometry::cellCenterIdx();
const auto ccGlobalI = gridGeometry.elementMapper().index(element);
-
+ pattern.add(ccGlobalI, ccGlobalI);
for (auto&& ccGlobalJ : connectivityMap(cellCenterIdx, cellCenterIdx, ccGlobalI))
pattern.add(ccGlobalI, ccGlobalJ);
}
@@ -150,6 +149,7 @@ auto getJacobianPattern(const GridGeometry& gridGeometry)
for (auto&& scvf : scvfs(fvGeometry))
{
const auto faceGlobalI = scvf.dofIndex();
+ pattern.add(faceGlobalI, faceGlobalI);
for (auto&& faceGlobalJ : connectivityMap(faceIdx, faceIdx, scvf.index()))
pattern.add(faceGlobalI, faceGlobalJ);
}
diff --git a/dumux/discretization/staggered/facesolution.hh b/dumux/discretization/staggered/facesolution.hh
index 71a486ec462afc104d742ca6631f47eb4254f26a..29d3bad86ca5629f76b8d810f9378597d0447fa3 100644
--- a/dumux/discretization/staggered/facesolution.hh
+++ b/dumux/discretization/staggered/facesolution.hh
@@ -50,9 +50,11 @@ public:
const auto& connectivityMap = fvGridGeometry.connectivityMap();
const auto& stencil = connectivityMap(FVGridGeometry::faceIdx(), FVGridGeometry::faceIdx(), scvf.index());
- facePriVars_.reserve(stencil.size());
- map_.reserve(stencil.size());
+ facePriVars_.reserve(stencil.size()+1);
+ map_.reserve(stencil.size()+1);
+ map_.push_back(scvf.dofIndex());
+ facePriVars_.push_back(sol[scvf.dofIndex()]);
for(const auto dofJ : stencil)
{
map_.push_back(dofJ);
diff --git a/dumux/discretization/staggered/freeflow/connectivitymap.hh b/dumux/discretization/staggered/freeflow/connectivitymap.hh
index 1aa82d24fe1e852e2225648e240e1eb4ac79be2b..573af548418aa8435c554ce165e0e9b321e51ad2 100644
--- a/dumux/discretization/staggered/freeflow/connectivitymap.hh
+++ b/dumux/discretization/staggered/freeflow/connectivitymap.hh
@@ -132,9 +132,6 @@ private:
const FVElementGeometry& fvGeometry,
const SubControlVolumeFace& scvf)
{
- // the first entry is always the cc dofIdx itself
- if(stencil.empty())
- stencil.push_back(scvf.insideScvIdx());
if(!scvf.boundary())
stencil.push_back(scvf.outsideScvIdx());
}
@@ -184,7 +181,6 @@ private:
{
if(stencil.empty())
{
- stencil.push_back(scvf.axisData().selfDof);
stencil.push_back(scvf.axisData().oppositeDof);
addHigherOrderInAxisDofs_(scvf, stencil, std::integral_constant<bool, useHigherOrder>{});
}
diff --git a/dumux/discretization/staggered/freeflow/elementvolumevariables.hh b/dumux/discretization/staggered/freeflow/elementvolumevariables.hh
index ccb55ea3f9a193ef6a8f2b5868c6dc76904822da..4da5e6a857379d1de801e0f45ef506f5c9f5505b 100644
--- a/dumux/discretization/staggered/freeflow/elementvolumevariables.hh
+++ b/dumux/discretization/staggered/freeflow/elementvolumevariables.hh
@@ -219,16 +219,15 @@ public:
auto&& scvI = fvGeometry.scv(globalI);
// resize local containers to the required size (for internal elements)
- volumeVariables_.resize(numDofs);
- volVarIndices_.resize(numDofs);
+ volumeVariables_.resize(numDofs+1);
+ volVarIndices_.resize(numDofs+1);
int localIdx = 0;
- // Update the volume variables of the element at hand and the neighboring elements
- for (auto globalJ : connectivityMapI)
+ // Lambda to update the volume variables of the given index
+ auto doVolVarUpdate = [&](int globalJ)
{
const auto& elementJ = fvGridGeometry.element(globalJ);
auto&& scvJ = fvGeometry.scv(globalJ);
-
const auto elemSol = makeElementSolutionFromCellCenterPrivars<PrimaryVariables>(sol[globalJ]);
volumeVariables_[localIdx].update(elemSol,
problem,
@@ -236,7 +235,14 @@ public:
scvJ);
volVarIndices_[localIdx] = scvJ.dofIndex();
++localIdx;
- }
+ };
+
+ // Update the volume variables of the element at hand
+ doVolVarUpdate(globalI);
+
+ // Update the volume variables of the neighboring elements
+ for (const auto& globalJ : connectivityMapI)
+ doVolVarUpdate(globalJ);
if (fvGeometry.hasBoundaryScvf())
{
diff --git a/dumux/multidomain/subdomainstaggeredlocalassembler.hh b/dumux/multidomain/subdomainstaggeredlocalassembler.hh
index c8d7db523a40d0b693294a182363b85b0194cb77..3b1e741274d42310573826fed18f0a4dabc3c528 100644
--- a/dumux/multidomain/subdomainstaggeredlocalassembler.hh
+++ b/dumux/multidomain/subdomainstaggeredlocalassembler.hh
@@ -538,11 +538,6 @@ public:
CellCenterResidualValue assembleCellCenterJacobianAndResidualImpl(JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
{
assert(domainI == cellCenterId);
- //////////////////////////////////////////////////////////////////////////////////////////////////
- // 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. //
- //////////////////////////////////////////////////////////////////////////////////////////////////
// get some aliases for convenience
const auto& element = this->element();
@@ -554,18 +549,12 @@ public:
const auto cellCenterGlobalI = fvGridGeometry.elementMapper().index(element);
const auto origResidual = this->evalLocalResidualForCellCenter();
- //////////////////////////////////////////////////////////////////////////////////////////////////
- // //
- // 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. //
- // //
- //////////////////////////////////////////////////////////////////////////////////////////////////
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Calculate derivatives of all cell center residuals in the element w.r.t. to other cell center dofs. //
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////
- // build derivatives with for cell center dofs w.r.t. cell center dofs
- const auto& connectivityMap = fvGridGeometry.connectivityMap();
-
- for (const auto& globalJ : connectivityMap(cellCenterId, cellCenterId, cellCenterGlobalI))
+ // lambda to evaluate the derivatives for cell center dofs with respect to neighbor cells
+ auto evaluateCellCenterDerivatives = [&](const std::size_t globalJ)
{
// get the volVars of the element with respect to which we are going to build the derivative
auto&& scvJ = fvGeometry.scv(globalJ);
@@ -615,7 +604,17 @@ public:
// restore the undeflected state of the coupling context
this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalJ, curSol[globalJ], pvIdx);
}
- }
+ };
+
+ // get the list of cell center dofs that have an influence on the cell center resdiual of the current element
+ const auto& connectivityMap = fvGridGeometry.connectivityMap();
+
+ // evaluate derivatives w.r.t. own dof
+ evaluateCellCenterDerivatives(cellCenterGlobalI);
+
+ // evaluate derivatives w.r.t. all other related cell center dofs
+ for (const auto& globalJ : connectivityMap(cellCenterId, cellCenterId, cellCenterGlobalI))
+ evaluateCellCenterDerivatives(globalJ);
return origResidual;
}
@@ -630,11 +629,6 @@ public:
auto assembleFaceJacobianAndResidualImpl(JacobianMatrixDiagBlock& A, GridVariables& gridVariables)
{
assert(domainI == faceId);
- //////////////////////////////////////////////////////////////////////////////////////////////////
- // 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. //
- //////////////////////////////////////////////////////////////////////////////////////////////////
// get some aliases for convenience
const auto& problem = this->problem();
@@ -652,16 +646,12 @@ public:
for (auto&& scvf : scvfs(fvGeometry))
origResiduals[scvf.localFaceIdx()] = this->evalLocalResidualForFace(scvf);
- //////////////////////////////////////////////////////////////////////////////////////////////////
- // //
- // 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. //
- // //
- //////////////////////////////////////////////////////////////////////////////////////////////////
-
- // build derivatives with for cell center dofs w.r.t. cell center dofs
- const auto& connectivityMap = fvGridGeometry.connectivityMap();
+ ///////////////////////////////////////////////////////////////////////////////////////////////////
+ // Calculate derivatives of all face residuals in the element w.r.t. to other face dofs. //
+ // Note that we do an element-wise assembly, therefore this is only the contribution of the //
+ // current element while the contribution of the element on the opposite side of the scvf will //
+ // be added separately. //
+ ///////////////////////////////////////////////////////////////////////////////////////////////////
for (auto&& scvf : scvfs(fvGeometry))
{
@@ -671,8 +661,8 @@ public:
using FaceSolution = GetPropType<TypeTag, Properties::StaggeredFaceSolution>;
const auto origFaceSolution = FaceSolution(scvf, curSol, fvGridGeometry);
- // build derivatives with for face dofs w.r.t. cell center dofs
- for (const auto& globalJ : connectivityMap(faceId, faceId, scvf.index()))
+ // Lambda to evaluate the derivatives for faces
+ auto evaluateFaceDerivatives = [&](const std::size_t globalJ)
{
// get the faceVars of the face with respect to which we are going to build the derivative
auto& faceVars = getFaceVarAccess(gridVariables.curGridFaceVars(), this->curElemFaceVars(), scvf);
@@ -711,9 +701,19 @@ public:
// restore the undeflected state of the coupling context
this->couplingManager().updateCouplingContext(domainI, *this, domainI, globalJ, origFaceSolution[globalJ], pvIdx);
- }
- }
- }
+ }
+ };
+
+ // evaluate derivatives w.r.t. own dof
+ evaluateFaceDerivatives(scvf.dofIndex());
+
+ // get the list of face dofs that have an influence on the resdiual of the current face
+ const auto& connectivityMap = fvGridGeometry.connectivityMap();
+
+ // evaluate derivatives w.r.t. all other related face dofs
+ for (const auto& globalJ : connectivityMap(faceId, faceId, scvf.index()))
+ evaluateFaceDerivatives(globalJ);
+ }
return origResiduals;
}
@@ -728,9 +728,9 @@ public:
void assembleJacobianCellCenterCoupling(Dune::index_constant<1> domainJ, JacobianBlock& A,
const CellCenterResidualValue& origResidual, GridVariables& gridVariables)
{
- // ////////////////////////////////////////////////////////////////////////////////////////////////////////
- // // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. //
- // ////////////////////////////////////////////////////////////////////////////////////////////////////////
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Calculate derivatives of all cell center residuals in the element w.r.t. to all coupled faces dofs //
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////
// get some aliases for convenience
const auto& element = this->element();
@@ -791,9 +791,9 @@ public:
void assembleJacobianCellCenterCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A,
const CellCenterResidualValue& res, GridVariables& gridVariables)
{
- ////////////////////////////////////////////////////////////////////////////////////////////////////////
- // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. //
- ////////////////////////////////////////////////////////////////////////////////////////////////////////
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Calculate derivatives of all cell center residuals in the element w.r.t. all dofs in the coupling stencil. //
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// get some aliases for convenience
const auto& element = this->element();
@@ -850,9 +850,9 @@ public:
void assembleJacobianFaceCoupling(Dune::index_constant<0> domainJ, JacobianBlock& A,
const ElementResidualVector& origResiduals, GridVariables& gridVariables)
{
- ////////////////////////////////////////////////////////////////////////////////////////////////////////
- // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. //
- ////////////////////////////////////////////////////////////////////////////////////////////////////////
+ /////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Calculate derivatives of all face residuals in the element w.r.t. all coupled cell center dofs //
+ /////////////////////////////////////////////////////////////////////////////////////////////////////
// get some aliases for convenience
const auto& problem = this->problem();
@@ -925,9 +925,9 @@ public:
void assembleJacobianFaceCoupling(Dune::index_constant<otherId> domainJ, JacobianBlock& A,
const ElementResidualVector& res, GridVariables& gridVariables)
{
- ////////////////////////////////////////////////////////////////////////////////////////////////////////
- // Calculate derivatives of all dofs in the element with respect to all dofs in the coupling stencil. //
- ////////////////////////////////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////
+ // Calculate derivatives of all face residuals in the element w.r.t. all dofs in the coupling stencil. //
+ //////////////////////////////////////////////////////////////////////////////////////////////////////////
// get some aliases for convenience
const auto& fvGeometry = this->fvGeometry();