From a78975bb738947b41cdf1f035cee61582330b425 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Fri, 28 Nov 2014 10:06:10 +0000
Subject: [PATCH] [cell-centered implicit] fix a bug in the handling of mixed
boundary conditions
Mixed boundary conditions here means that a part of the
equations/primary variables gets Dirichlet conditions, while the rest
gets Neumann conditions.
While in the cell-centered method pure Dirichlet conditions for all
equations are handled by calculating the resulting fluxes and adding
them to the cell residual, such a flux cannot be (easily?) calculated
for an equation that gets a Dirichlet condition as part of mixed
conditions. Therefore, such a Dirichlet condition is implemented in a
strong way by _replacing_ the corresponding cell residual.
It is important that this replacement is done at the very end of the
residual calculation. However, for corner cells this has not been
guaranteed so far. Therefore, fluxed resulting from the other boundary
parts of a corner cell could have been added to the replaced residual,
obviously leading to a wrong boundary condition treatment.
This patch resolves the issue by guaranteeing that the residual
replacement is done at the end of the residual calculation.
Brought to attention and reviewed by Thomas. Thanks.
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@13823 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
.../implicit/cellcentered/cclocalresidual.hh | 101 ++++++++++++------
.../implicit/common/implicitlocalresidual.hh | 4 +
2 files changed, 70 insertions(+), 35 deletions(-)
diff --git a/dumux/implicit/cellcentered/cclocalresidual.hh b/dumux/implicit/cellcentered/cclocalresidual.hh
index 242a7a7493..aa5d3048dd 100644
--- a/dumux/implicit/cellcentered/cclocalresidual.hh
+++ b/dumux/implicit/cellcentered/cclocalresidual.hh
@@ -68,8 +68,8 @@ public:
protected:
/*!
- * \brief Add all Neumann and outflow boundary conditions to the local
- * residual.
+ * \brief Add all fluxes resulting from Neumann, outflow and pure Dirichlet
+ * boundary conditions to the local residual.
*/
void evalBoundaryFluxes_()
{
@@ -92,9 +92,31 @@ protected:
if (bcTypes.hasOutflow())
this->asImp_().evalOutflowSegment_(isIt, bcTypes);
- // evaluate the Dirichlet conditions at the boundary face
- if (bcTypes.hasDirichlet())
- this->asImp_().evalDirichletSegment_(isIt, bcTypes);
+ // evaluate the pure Dirichlet conditions at the boundary face
+ if (bcTypes.hasDirichlet() && !bcTypes.hasNeumann())
+ this->asImp_().evalDirichletSegmentPure_(isIt, bcTypes);
+ }
+ }
+
+ /*!
+ * \brief Evaluate Dirichlet conditions on faces that have mixed
+ * Dirichlet/Neumann boundary conditions.
+ */
+ void evalDirichlet_()
+ {
+ IntersectionIterator isIt = this->gridView_().ibegin(this->element_());
+ const IntersectionIterator &isEndIt = this->gridView_().iend(this->element_());
+ for (; isIt != isEndIt; ++isIt)
+ {
+ // handle only faces on the boundary
+ if (!isIt->boundary())
+ continue;
+
+ BoundaryTypes bcTypes;
+ this->problem_().boundaryTypes(bcTypes, *isIt);
+
+ if (bcTypes.hasDirichlet() && bcTypes.hasNeumann())
+ this->asImp_().evalDirichletSegmentMixed_(isIt, bcTypes);
}
}
@@ -151,7 +173,7 @@ protected:
// manipulate the corresponding subcontrolvolume face
SCVFace& boundaryFace = fvGeometry.boundaryFace[bfIdx];
-
+
// set the second flux approximation index for the boundary face
for (int nIdx = 0; nIdx < fvGeometry.numNeighbors-1; nIdx++)
{
@@ -186,45 +208,54 @@ protected:
}
/*!
- * \brief Add Dirichlet boundary conditions for a single intersection
+ * \brief Treat Dirichlet boundary conditions in a weak sense for a single
+ * intersection that only has Dirichlet boundary conditions
*/
- void evalDirichletSegment_(const IntersectionIterator &isIt,
- const BoundaryTypes &bcTypes)
+ void evalDirichletSegmentPure_(const IntersectionIterator &isIt,
+ const BoundaryTypes &bcTypes)
{
// temporary vector to store the Dirichlet boundary fluxes
PrimaryVariables values;
Valgrind::SetUndefined(values);
unsigned bfIdx = isIt->indexInInside();
-
- // check for mixed Dirichlet/Neumann conditions
- if (bcTypes.hasNeumann())
+
+ this->asImp_().computeFlux(values, bfIdx, true);
+ values *= this->curVolVars_(0).extrusionFactor();
+
+ // add fluxes to the residual
+ Valgrind::CheckDefined(values);
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
{
- this->problem_().dirichlet(values, *isIt);
- Valgrind::CheckDefined(values);
-
- // set Dirichlet conditions in a strong sense
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
- {
- if (bcTypes.isDirichlet(eqIdx))
- {
- int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
- this->residual_[0][eqIdx]
- = this->curPriVar_(0, pvIdx) - values[pvIdx];
- }
- }
+ if (bcTypes.isDirichlet(eqIdx))
+ this->residual_[0][eqIdx] += values[eqIdx];
}
- else // pure Dirichlet conditions
- {
- this->asImp_().computeFlux(values, bfIdx, true);
- values *= this->curVolVars_(0).extrusionFactor();
+ }
- // add fluxes to the residual
- Valgrind::CheckDefined(values);
- for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ /*!
+ * \brief Treat Dirichlet boundary conditions in a strong sense for a
+ * single intersection that has mixed D/N boundary conditions
+ */
+ void evalDirichletSegmentMixed_(const IntersectionIterator &isIt,
+ const BoundaryTypes &bcTypes)
+ {
+ // temporary vector to store the Dirichlet boundary fluxes
+ PrimaryVariables values;
+ Valgrind::SetUndefined(values);
+
+ unsigned bfIdx = isIt->indexInInside();
+
+ this->problem_().dirichlet(values, *isIt);
+ Valgrind::CheckDefined(values);
+
+ // set Dirichlet conditions in a strong sense
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (bcTypes.isDirichlet(eqIdx))
{
- if (bcTypes.isDirichlet(eqIdx))
- this->residual_[0][eqIdx] += values[eqIdx];
+ int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
+ this->residual_[0][eqIdx]
+ = this->curPriVar_(0, pvIdx) - values[pvIdx];
}
}
}
@@ -244,7 +275,7 @@ protected:
continue;
fIdx++;
- PrimaryVariables flux;
+ PrimaryVariables flux;
Valgrind::SetUndefined(flux);
this->asImp_().computeFlux(flux, fIdx);
diff --git a/dumux/implicit/common/implicitlocalresidual.hh b/dumux/implicit/common/implicitlocalresidual.hh
index 5f50c79825..cb77d3eda9 100644
--- a/dumux/implicit/common/implicitlocalresidual.hh
+++ b/dumux/implicit/common/implicitlocalresidual.hh
@@ -311,6 +311,10 @@ protected:
else
{
asImp_().evalBoundaryFluxes_();
+
+ // additionally treat mixed D/N conditions in a strong sense
+ if (bcTypes_().hasDirichlet())
+ asImp_().evalDirichlet_();
}
#if !defined NDEBUG && HAVE_VALGRIND
--
GitLab