Commit a78975bb authored by Bernd Flemisch's avatar Bernd Flemisch
Browse files

[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
parent d27d40a5
......@@ -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);
}
}
......@@ -186,9 +208,10 @@ 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,
void evalDirichletSegmentPure_(const IntersectionIterator &isIt,
const BoundaryTypes &bcTypes)
{
// temporary vector to store the Dirichlet boundary fluxes
......@@ -197,34 +220,42 @@ protected:
unsigned bfIdx = isIt->indexInInside();
// check for mixed Dirichlet/Neumann conditions
if (bcTypes.hasNeumann())
{
this->problem_().dirichlet(values, *isIt);
Valgrind::CheckDefined(values);
this->asImp_().computeFlux(values, bfIdx, true);
values *= this->curVolVars_(0).extrusionFactor();
// set Dirichlet conditions in a strong sense
// add fluxes to the residual
Valgrind::CheckDefined(values);
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];
}
this->residual_[0][eqIdx] += values[eqIdx];
}
}
else // pure Dirichlet conditions
/*!
* \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)
{
this->asImp_().computeFlux(values, bfIdx, true);
values *= this->curVolVars_(0).extrusionFactor();
// temporary vector to store the Dirichlet boundary fluxes
PrimaryVariables values;
Valgrind::SetUndefined(values);
// add fluxes to the residual
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))
this->residual_[0][eqIdx] += values[eqIdx];
{
int pvIdx = bcTypes.eqToDirichletIndex(eqIdx);
this->residual_[0][eqIdx]
= this->curPriVar_(0, pvIdx) - values[pvIdx];
}
}
}
......
......@@ -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
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment