Commit dd683a8a authored by Thomas Fetzer's avatar Thomas Fetzer
Browse files

[multidomain]

- removed Beavers-Joseph dependency for the handling of the
couplingInflow and couplingOutflow boundaries
- removed todo and unused lines of code
- updated reference solutions for the multidomain problem (only small
deviations from reference)

reviewed by gruenich, bernd



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@14102 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent de6b47fc
......@@ -176,9 +176,7 @@ public:
cParams,
couplingRes1, couplingRes2);
// TODO: this should be done only once
const DimVector& globalPos2 = cParams.fvGeometry2.subContVol[vertInElem2].global;
// const DimVector& bfNormal2 = boundaryVars2.face().normal;
DimVector normalMassFlux2(0.);
// velocity*normal*area*rho
......
......@@ -196,7 +196,6 @@ class TwoCStokesTwoPTwoCLocalOperator :
const Dune::ReferenceElement<typename MDGrid::ctype,dim>& referenceElement2 =
Dune::ReferenceElements<typename MDGrid::ctype,dim>::general(mdElement2.type());
// TODO: assumes same number of vertices on a coupling face
for (int vertexInFace = 0; vertexInFace < numVerticesOfFace; ++vertexInFace)
{
const int vertInElem1 = referenceElement1.subEntity(faceIdx1, 1, vertexInFace, dim);
......@@ -343,17 +342,6 @@ class TwoCStokesTwoPTwoCLocalOperator :
}
if (cParams.boundaryTypes2.isCouplingOutflow(massBalanceIdx2))
{
//TODO what happens at the corners? Idea: set only pressure in corners, so that residual is not needed there
//pi-tau
// if (globalProblem_.sdProblem1().isCornerPoint(globalPos1))
// {
// couplingRes2[getIndex_<LFSU2,massBalanceIdx2> (lfsu_n,vertInElem2)] -= cParams.elemVolVarsCur1[vertInElem1].pressure;
// }
// else
// {
// // n.(pI-tau)n as dirichlet condition for Darcy p (set, if B&J defined as Dirichlet condition)
// set residualDarcy[massBalance] = p in 2p2clocalresidual.hh
couplingRes2.accumulate(lfsu_n.child(massBalanceIdx2), vertInElem2,
globalProblem_.localResidual1().residual(vertInElem1)[momentumYIdx1]
-cParams.elemVolVarsCur1[vertInElem1].pressure());
......@@ -513,13 +501,11 @@ class TwoCStokesTwoPTwoCLocalOperator :
couplingRes1.accumulate(lfsu_s.child(momentumYIdx1), vertInElem1,
globalProblem_.localResidual2().residual(vertInElem2)[massBalanceIdx2]
/ cParams.elemVolVarsCur1[vertInElem1].density());
// TODO: * bfNormal2.two_norm());
}
}
//coupling residual is added to "real" residual
//here each node is passed twice, hence only half of the dirichlet condition has to be set
//TODO what to do at boundary nodes which appear only once?
if (cParams.boundaryTypes1.isCouplingOutflow(transportEqIdx1))
{
// set residualStokes[transportEqIdx1] = x in stokes2clocalresidual.hh
......
......@@ -115,7 +115,7 @@ NEW_PROP_TAG(MultiDomainCoupling);
//! Property tag for the multidomain constraints transformation
NEW_PROP_TAG(MultiDomainConstraintsTrafo);
NEW_PROP_TAG(ConstraintsTrafo); // TODO: required?
NEW_PROP_TAG(ConstraintsTrafo);
//! Specifies the type of the jacobian matrix as used for the linear solver
NEW_PROP_TAG(JacobianMatrix);
......
......@@ -292,7 +292,7 @@ public:
if (this->bcTypes_(scvIdx).isCouplingOutflow(massBalanceIdx))
this->residual_[scvIdx][massBalanceIdx] = volVars.pressure(nPhaseIdx);
// set mass fraction; TODO: this is fixed to contiWEqIdx so far!
// set mass fraction;
if (this->bcTypes_(scvIdx).isCouplingOutflow(contiWEqIdx))
this->residual_[scvIdx][contiWEqIdx] = volVars.massFraction(nPhaseIdx, wCompIdx);
}
......
......@@ -300,7 +300,7 @@ public:
if (this->bcTypes_(scvIdx).isCouplingOutflow(massBalanceIdx))
this->residual_[scvIdx][massBalanceIdx] = volVars.pressure(nPhaseIdx);
// set mass fraction; TODO: this is fixed to contiWEqIdx so far!
// set mass fraction;
if (this->bcTypes_(scvIdx).isCouplingOutflow(contiWEqIdx))
this->residual_[scvIdx][contiWEqIdx] = volVars.massFraction(nPhaseIdx, wCompIdx);
......
......@@ -18,7 +18,8 @@
*****************************************************************************/
/*!
* \file
* \brief Calculates the residual of models based on the box scheme element-wise.
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the coupled box model.
*/
#ifndef DUMUX_BOX_COUPLING_LOCAL_RESIDUAL_HH
#define DUMUX_BOX_COUPLING_LOCAL_RESIDUAL_HH
......@@ -30,10 +31,8 @@ namespace Dumux
/*!
* \ingroup ImplicitLocalResidual
* \ingroup TwoPTwoCNIStokesTwoCNIModel
* \brief Element-wise calculation of the residual matrix for models
* based on the box scheme .
*
* \todo Please doc me more!
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the coupled box model.
*/
template<class TypeTag>
class BoxCouplingLocalResidual : public BoxLocalResidual<TypeTag>
......
......@@ -20,7 +20,7 @@
* \file
*
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the stokes box model.
* using the coupled compositional Stokes box model.
*/
#ifndef DUMUX_STOKESNC_COUPLING_LOCAL_RESIDUAL_HH
......@@ -34,13 +34,8 @@ namespace Dumux
* \ingroup ImplicitLocalResidual
* \ingroup TwoPTwoCStokesTwoCModel
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the Stokes box model.
*
* This class is also used for the stokes transport
* model, which means that it uses static polymorphism.
*
* \todo Please doc me more!
* This file should contain a more detailed description of the coupling conditions.
* using the coupled compositional Stokes box model.
* It is derived from the compositional Stokes box model.
*/
template<class TypeTag>
class StokesncCouplingLocalResidual : public StokesncLocalResidual<TypeTag>
......@@ -272,17 +267,8 @@ protected:
const VolumeVariables &volVars = this->curVolVars_()[scvIdx];
const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
// TODO: beaversJosephCoeff is used to determine, if we are on a coupling face
// this is important at the corners. However, a better way should be found.
// check if BJ-coeff is not zero to decide, if coupling condition
// for the momentum balance (Dirichlet vor v.n) has to be set;
// may be important at corners
const GlobalPosition &globalPos = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeffAtPos(globalPos);
// set velocity normal to the interface
if (bcTypes.isCouplingInflow(momentumYIdx) && beaversJosephCoeff)
if (bcTypes.isCouplingInflow(momentumYIdx))
this->residual_[scvIdx][momentumYIdx] =
volVars.velocity() *
boundaryVars.face().normal /
......@@ -291,10 +277,10 @@ protected:
// add pressure correction - required for pressure coupling,
// if p.n comes from the pm
if ((bcTypes.isCouplingOutflow(momentumYIdx) && beaversJosephCoeff) || bcTypes.isMortarCoupling(momentumYIdx))
if (bcTypes.isCouplingOutflow(momentumYIdx) || bcTypes.isMortarCoupling(momentumYIdx))
{
DimVector pressureCorrection(isIt->centerUnitOuterNormal());
pressureCorrection *= volVars.pressure(); // TODO: 3D
pressureCorrection *= volVars.pressure();
this->residual_[scvIdx][momentumYIdx] += pressureCorrection[momentumYIdx]*
boundaryVars.face().area;
}
......@@ -319,69 +305,65 @@ protected:
void evalBeaversJoseph_(const IntersectionIterator &isIt,
const int scvIdx,
const int boundaryFaceIdx,
const FluxVariables &boundaryVars) //TODO: required
const FluxVariables &boundaryVars)
{
const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
const GlobalPosition &globalPos = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeffAtPos(globalPos);
// only enter here, if we are on a coupling boundary (see top)
// and the BJ coefficient is not zero
if (beaversJosephCoeff)
const Scalar Kxx = this->problem_().permeability(this->element_(),
this->fvGeometry_(),
scvIdx);
beaversJosephCoeff /= std::sqrt(Kxx);
const DimVector& elementUnitNormal = isIt->centerUnitOuterNormal();
// implementation as NEUMANN condition /////////////////////////////////////////////
// (v.n)n
if (bcTypes.isCouplingOutflow(momentumXIdx))
{
const Scalar Kxx = this->problem_().permeability(this->element_(),
this->fvGeometry_(),
scvIdx);
beaversJosephCoeff /= sqrt(Kxx);
const DimVector& elementUnitNormal = isIt->centerUnitOuterNormal();
// implementation as NEUMANN condition /////////////////////////////////////////////
// (v.n)n
if (bcTypes.isCouplingOutflow(momentumXIdx))
{
const Scalar normalComp = boundaryVars.velocity()*elementUnitNormal;
DimVector normalV = elementUnitNormal;
normalV *= normalComp; // v*n*n
// v_tau = v - (v.n)n
const DimVector tangentialV = boundaryVars.velocity() - normalV;
const Scalar boundaryFaceArea = boundaryVars.face().area;
for (int dimIdx=0; dimIdx < dim; ++dimIdx)
this->residual_[scvIdx][dimIdx] += beaversJosephCoeff *
boundaryFaceArea *
tangentialV[dimIdx] *
(boundaryVars.dynamicViscosity()
+ boundaryVars.dynamicEddyViscosity());
////////////////////////////////////////////////////////////////////////////////////
//normal component has only to be set if no coupling conditions are defined
//set NEUMANN flux (set equal to pressure in problem)
// PrimaryVariables priVars(0.0);
// this->problem_().neumann(priVars, this->element_(), this->fvGeometry_(),
// *isIt, scvIdx, boundaryFaceIdx);
// for (int dimIdx=0; dimIdx < dim; ++dimIdx)
// this->residual_[scvIdx][dimIdx] += priVars[dimIdx]*
// boundaryFaceArea;
}
if (bcTypes.isCouplingInflow(momentumXIdx)) //TODO: 3D
{
///////////////////////////////////////////////////////////////////////////////////////////
//IMPLEMENTATION AS DIRICHLET CONDITION
//tangential componment: vx = sqrt K /alpha * (grad v n(unity))t
DimVector tangentialVelGrad(0);
boundaryVars.velocityGrad().umv(elementUnitNormal, tangentialVelGrad);
tangentialVelGrad /= -beaversJosephCoeff; // was - before
this->residual_[scvIdx][momentumXIdx] =
tangentialVelGrad[momentumXIdx] - this->curPriVars_(scvIdx)[momentumXIdx];
////////////////////////////////////////////////////////////////////////////////////
//for testing Beavers and Joseph without adjacent porous medium set vy to 0
const Scalar normalComp = boundaryVars.velocity()*elementUnitNormal;
DimVector normalV = elementUnitNormal;
normalV *= normalComp; // v*n*n
// v_tau = v - (v.n)n
const DimVector tangentialV = boundaryVars.velocity() - normalV;
const Scalar boundaryFaceArea = boundaryVars.face().area;
for (int dimIdx=0; dimIdx < dim; ++dimIdx)
this->residual_[scvIdx][dimIdx] += beaversJosephCoeff *
boundaryFaceArea *
tangentialV[dimIdx] *
(boundaryVars.dynamicViscosity()
+ boundaryVars.dynamicEddyViscosity());
////////////////////////////////////////////////////////////////////////////////////
//normal component has only to be set if no coupling conditions are defined
//set NEUMANN flux (set equal to pressure in problem)
// PrimaryVariables priVars(0.0);
// this->problem_().neumann(priVars, this->element_(), this->fvGeometry_(),
// *isIt, scvIdx, boundaryFaceIdx);
// for (int dimIdx=0; dimIdx < dim; ++dimIdx)
// this->residual_[scvIdx][dimIdx] += priVars[dimIdx]*
// boundaryFaceArea;
}
if (bcTypes.isCouplingInflow(momentumXIdx))
{
assert(beaversJosephCoeff > 0);
///////////////////////////////////////////////////////////////////////////////////////////
//IMPLEMENTATION AS DIRICHLET CONDITION
//tangential componment: vx = sqrt K /alpha * (grad v n(unity))t
DimVector tangentialVelGrad(0);
boundaryVars.velocityGrad().umv(elementUnitNormal, tangentialVelGrad);
tangentialVelGrad /= -beaversJosephCoeff; // was - before
this->residual_[scvIdx][momentumXIdx] =
tangentialVelGrad[momentumXIdx] - this->curPriVars_(scvIdx)[momentumXIdx];
////////////////////////////////////////////////////////////////////////////////////
//for testing Beavers and Joseph without adjacent porous medium set vy to 0
// Scalar normalVel(0);
// this->residual_[scvIdx][momentumYIdx] = this->curPrimaryVars_(scvIdx)[momentumYIdx] - normalVel;
////////////////////////////////////////////////////////////////////////////////////
}
////////////////////////////////////////////////////////////////////////////////////
}
}
......
......@@ -20,7 +20,7 @@
* \file
*
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the stokes box model.
* using the coupled compositional non-isothermal Stokes box model.
*/
#ifndef DUMUX_STOKESNCNI_COUPLING_LOCAL_RESIDUAL_HH
......@@ -34,14 +34,9 @@ namespace Dumux
/*!
* \ingroup ImplicitLocalResidual
* \ingroup TwoPTwoCNIStokesTwoCNIModel
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the Stokes box model.
*
* This class is also used for the stokes transport
* model, which means that it uses static polymorphism.
*
* \todo Please doc me more!
* This file should contain a more detailed description of the coupling conditions.
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the coupled compositional non-isothermal Stokes box model.
* It is derived from the compositional non-isothermal Stokes box model.
*/
template<class TypeTag>
class StokesncniCouplingLocalResidual : public StokesncniLocalResidual<TypeTag>
......@@ -277,18 +272,9 @@ namespace Dumux
{
const VolumeVariables &volVars = this->curVolVars_()[scvIdx];
const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
// TODO: beaversJosephCoeff is used to determine, if we are on a coupling face
// this is important at the corners. However, a better way should be found.
// check if BJ-coeff is not zero to decide, if coupling condition
// for the momentum balance (Dirichlet vor v.n) has to be set;
// may be important at corners
const GlobalPosition &globalPos = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeffAtPos(globalPos);
// set velocity normal to the interface
if (bcTypes.isCouplingInflow(momentumYIdx) && beaversJosephCoeff)
if (bcTypes.isCouplingInflow(momentumYIdx))
this->residual_[scvIdx][momentumYIdx] =
volVars.velocity() *
boundaryVars.face().normal /
......@@ -297,10 +283,10 @@ namespace Dumux
// add pressure correction - required for pressure coupling,
// if p.n comes from the pm
if ((bcTypes.isCouplingOutflow(momentumYIdx) && beaversJosephCoeff) || bcTypes.isMortarCoupling(momentumYIdx))
if (bcTypes.isCouplingOutflow(momentumYIdx) || bcTypes.isMortarCoupling(momentumYIdx))
{
DimVector pressureCorrection(isIt->centerUnitOuterNormal());
pressureCorrection *= volVars.pressure(); // TODO: 3D
pressureCorrection *= volVars.pressure();
this->residual_[scvIdx][momentumYIdx] += pressureCorrection[momentumYIdx]*
boundaryVars.face().area;
}
......@@ -331,69 +317,66 @@ namespace Dumux
void evalBeaversJoseph_(const IntersectionIterator &isIt,
const int scvIdx,
const int boundaryFaceIdx,
const FluxVariables &boundaryVars) //TODO: required
const FluxVariables &boundaryVars)
{
const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
const GlobalPosition &globalPos = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeffAtPos(globalPos);
const GlobalPosition &globalPos = this->fvGeometry_().boundaryFace[boundaryFaceIdx].ipGlobal;
Scalar beaversJosephCoeff = this->problem_().beaversJosephCoeffAtPos(globalPos);
// only enter here, if we are on a coupling boundary (see top)
// and the BJ coefficient is not zero
if (beaversJosephCoeff)
{
const Scalar Kxx = this->problem_().permeability(this->element_(),
this->fvGeometry_(),
scvIdx);
beaversJosephCoeff /= sqrt(Kxx);
const DimVector& elementUnitNormal = isIt->centerUnitOuterNormal();
// implementation as NEUMANN condition /////////////////////////////////////////////
// (v.n)n
if (bcTypes.isCouplingOutflow(momentumXIdx))
{
const Scalar normalComp = boundaryVars.velocity()*elementUnitNormal;
DimVector normalV = elementUnitNormal;
normalV *= normalComp; // v*n*n
// v_tau = v - (v.n)n
const DimVector tangentialV = boundaryVars.velocity() - normalV;
const Scalar boundaryFaceArea = boundaryVars.face().area;
for (int dimIdx=0; dimIdx < dim; ++dimIdx)
this->residual_[scvIdx][dimIdx] += beaversJosephCoeff *
boundaryFaceArea *
tangentialV[dimIdx] *
(boundaryVars.dynamicViscosity() + boundaryVars.dynamicEddyViscosity());
////////////////////////////////////////////////////////////////////////////////////
//normal component has only to be set if no coupling conditions are defined
//set NEUMANN flux (set equal to pressure in problem)
// PrimaryVariables priVars(0.0);
// this->problem_().neumann(priVars, this->element_(), this->fvGeometry_(),
// *isIt, scvIdx, boundaryFaceIdx);
// for (int dimIdx=0; dimIdx < dim; ++dimIdx)
// this->residual_[scvIdx][dimIdx] += priVars[dimIdx]*
// boundaryFaceArea;
}
if (bcTypes.isCouplingInflow(momentumXIdx)) //TODO: 3D
{
///////////////////////////////////////////////////////////////////////////////////////////
//IMPLEMENTATION AS DIRICHLET CONDITION
//tangential componment: vx = sqrt K /alpha * (grad v n(unity))t
DimVector tangentialVelGrad(0);
boundaryVars.velocityGrad().umv(elementUnitNormal, tangentialVelGrad);
tangentialVelGrad /= -beaversJosephCoeff; // was - before
this->residual_[scvIdx][momentumXIdx] =
tangentialVelGrad[momentumXIdx] - this->curPriVars_(scvIdx)[momentumXIdx];
////////////////////////////////////////////////////////////////////////////////////
//for testing Beavers and Joseph without adjacent porous medium set vy to 0
// Scalar normalVel(0);
// this->residual_[scvIdx][momentumYIdx] = this->curPrimaryVars_(scvIdx)[momentumYIdx] - normalVel;
////////////////////////////////////////////////////////////////////////////////////
}
}
const Scalar Kxx = this->problem_().permeability(this->element_(),
this->fvGeometry_(),
scvIdx);
beaversJosephCoeff /= sqrt(Kxx);
assert(beaversJosephCoeff > 0.0);
const DimVector& elementUnitNormal = isIt->centerUnitOuterNormal();
// implementation as NEUMANN condition /////////////////////////////////////////////
// (v.n)n
if (bcTypes.isCouplingOutflow(momentumXIdx))
{
const Scalar normalComp = boundaryVars.velocity()*elementUnitNormal;
DimVector normalV = elementUnitNormal;
normalV *= normalComp; // v*n*n
// v_tau = v - (v.n)n
const DimVector tangentialV = boundaryVars.velocity() - normalV;
const Scalar boundaryFaceArea = boundaryVars.face().area;
for (int dimIdx=0; dimIdx < dim; ++dimIdx)
this->residual_[scvIdx][dimIdx] += beaversJosephCoeff *
boundaryFaceArea *
tangentialV[dimIdx] *
(boundaryVars.dynamicViscosity() + boundaryVars.dynamicEddyViscosity());
////////////////////////////////////////////////////////////////////////////////////
//normal component has only to be set if no coupling conditions are defined
//set NEUMANN flux (set equal to pressure in problem)
// PrimaryVariables priVars(0.0);
// this->problem_().neumann(priVars, this->element_(), this->fvGeometry_(),
// *isIt, scvIdx, boundaryFaceIdx);
// for (int dimIdx=0; dimIdx < dim; ++dimIdx)
// this->residual_[scvIdx][dimIdx] += priVars[dimIdx]*
// boundaryFaceArea;
}
if (bcTypes.isCouplingInflow(momentumXIdx))
{
assert(beaversJosephCoeff > 0);
///////////////////////////////////////////////////////////////////////////////////////////
//IMPLEMENTATION AS DIRICHLET CONDITION
//tangential componment: vx = sqrt K /alpha * (grad v n(unity))t
DimVector tangentialVelGrad(0);
boundaryVars.velocityGrad().umv(elementUnitNormal, tangentialVelGrad);
tangentialVelGrad /= -beaversJosephCoeff; // was - before
this->residual_[scvIdx][momentumXIdx] =
tangentialVelGrad[momentumXIdx] - this->curPriVars_(scvIdx)[momentumXIdx];
////////////////////////////////////////////////////////////////////////////////////
//for testing Beavers and Joseph without adjacent porous medium set vy to 0
// Scalar normalVel(0);
// this->residual_[scvIdx][momentumYIdx] = this->curPrimaryVars_(scvIdx)[momentumYIdx] - normalVel;
////////////////////////////////////////////////////////////////////////////////////
}
}
// return true, if at least one equation on the boundary has a coupling condition
......
......@@ -357,10 +357,7 @@ public:
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
{
if (onLowerBoundary_(globalPos))
return alphaBJ_;
else
return 0.0;
return alphaBJ_;
}
/*!
......
......@@ -325,10 +325,7 @@ public:
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
{
if (onLowerBoundary_(globalPos))
return alphaBJ_;
else
return 0.0;
return alphaBJ_;
}
/*!
......
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