From ec4382eb9fcb2a7ef2b28bbcbdae46aefba02cb8 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Tue, 26 Jul 2022 15:14:11 +0200
Subject: [PATCH] [fcdiamond][assembler] Handle ghosts in local assembler
---
dumux/assembly/fcdiamondlocalassembler.hh | 40 +++++++++++++++++------
1 file changed, 30 insertions(+), 10 deletions(-)
diff --git a/dumux/assembly/fcdiamondlocalassembler.hh b/dumux/assembly/fcdiamondlocalassembler.hh
index 77ec484651..e8fd8879fb 100644
--- a/dumux/assembly/fcdiamondlocalassembler.hh
+++ b/dumux/assembly/fcdiamondlocalassembler.hh
@@ -144,8 +144,27 @@ public:
maybeAssembleCouplingBlocks(residual);
}
else
- DUNE_THROW(Dune::NotImplemented, "Ghost elements not supported");
+ {
+ const auto numLocalFaces = this->element().subEntities(1);
+ for (int i = 0; i < numLocalFaces; ++i)
+ {
+ // do not change the non-ghost vertices
+ const auto face = this->element().template subEntity<1>(i);
+ if (face.partitionType() == Dune::InteriorEntity || face.partitionType() == Dune::BorderEntity)
+ continue;
+
+ // set main diagonal entries for the vertex
+ const auto dofIndex = gridGeometry.dofMapper().index(face);
+ typedef typename JacobianMatrix::block_type BlockType;
+ BlockType &J = jac[dofIndex][dofIndex];
+ for (int j = 0; j < BlockType::rows; ++j)
+ J[j][j] = 1.0;
+
+ // set residual for the vertex
+ res[dofIndex] = 0;
+ }
+ }
auto applyDirichlet = [&] (const auto& scvI,
const auto& dirichletValues,
@@ -390,14 +409,14 @@ public:
{
// update the volume variables and compute element residual
elemSol[scv.localDofIndex()][pvIdx] = priVar;
- deflectionHelper.deflect(elemSol, scv, this->problem());
+ deflectionHelper.deflect(elemSol, scv, problem);
this->asImp_().maybeUpdateCouplingContext(scv, elemSol, pvIdx);
return this->evalLocalResidual();
};
// derive the residuals numerically
- static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()};
- static const int numDiffMethod = getParamFromGroup<int>(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod");
+ static const NumericEpsilon<Scalar, numEq> eps_{problem.paramGroup()};
+ static const int numDiffMethod = getParamFromGroup<int>(problem.paramGroup(), "Assembly.NumericDifferenceMethod");
NumericDifferentiation::partialDerivative(evalResiduals, elemSol[scv.localDofIndex()][pvIdx], partialDerivs, origResiduals,
eps_(elemSol[scv.localDofIndex()][pvIdx], pvIdx), numDiffMethod);
@@ -468,9 +487,10 @@ public:
DUNE_THROW(Dune::NotImplemented, "partial reassembly for explicit time discretization");
// get some aliases for convenience
+ const auto& problem = this->asImp_().problem();
const auto& element = this->element();
const auto& fvGeometry = this->fvGeometry();
- const auto& curSol = this->curSol();
+ const auto& curSol = this->asImp_().curSol();
auto&& curElemVolVars = this->curElemVolVars();
// get the vecor of the acutal element residuals
@@ -508,13 +528,13 @@ public:
{
// auto partialDerivsTmp = partialDerivs;
elemSol[scv.localDofIndex()][pvIdx] = priVar;
- curVolVars.update(elemSol, this->problem(), element, scv);
+ curVolVars.update(elemSol, problem, element, scv);
return this->evalLocalStorageResidual();
};
// derive the residuals numerically
- static const NumericEpsilon<Scalar, numEq> eps_{this->problem().paramGroup()};
- static const int numDiffMethod = getParamFromGroup<int>(this->problem().paramGroup(), "Assembly.NumericDifferenceMethod");
+ static const NumericEpsilon<Scalar, numEq> eps_{problem.paramGroup()};
+ static const int numDiffMethod = getParamFromGroup<int>(problem.paramGroup(), "Assembly.NumericDifferenceMethod");
NumericDifferentiation::partialDerivative(evalStorage, elemSol[scv.localDofIndex()][pvIdx], partialDerivs, origStorageResiduals,
eps_(elemSol[scv.localDofIndex()][pvIdx], pvIdx), numDiffMethod);
@@ -579,7 +599,7 @@ public:
// get some aliases for convenience
const auto& element = this->element();
const auto& fvGeometry = this->fvGeometry();
- const auto& problem = this->problem();
+ const auto& problem = this->asImp_().problem();
const auto& curElemVolVars = this->curElemVolVars();
const auto& elemFluxVarsCache = this->elemFluxVarsCache();
@@ -699,7 +719,7 @@ public:
// get some aliases for convenience
const auto& element = this->element();
const auto& fvGeometry = this->fvGeometry();
- const auto& problem = this->problem();
+ const auto& problem = this->asImp_().problem();
const auto& curElemVolVars = this->curElemVolVars();
// get the vector of the actual element residuals
--
GitLab