diff --git a/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh b/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
index 7675b76053ebacb9afdf9757add290bce7964cee..1f402e26b1a4142163c16b0102bcbbfabe665bce 100644
--- a/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
+++ b/dumux/multidomain/2cnistokes2p2cni/2cnistokes2p2cnilocaloperator.hh
@@ -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
diff --git a/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh b/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
index 5c750244cad3ce5e9b8085c3ff62af098978a776..2f1304ef46296cc3d4352b8e8320ad8095532e5c 100644
--- a/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
+++ b/dumux/multidomain/2cstokes2p2c/2cstokes2p2clocaloperator.hh
@@ -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
diff --git a/dumux/multidomain/common/multidomainproperties.hh b/dumux/multidomain/common/multidomainproperties.hh
index 650aa46c3ef93ea1b34c7f55e13246d59cb0af5e..373062fc7e514aef02728bc2c05685cc4890b493 100644
--- a/dumux/multidomain/common/multidomainproperties.hh
+++ b/dumux/multidomain/common/multidomainproperties.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);
diff --git a/dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh
index 5fb0885a98a9b03e2341903643be9aa275cdf099..c2b0d168e3a8c25507a0ddd2cbe65fe1a30621b3 100644
--- a/dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/2p2ccouplinglocalresidual.hh
@@ -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);
}
diff --git a/dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh
index 036ddcd9a4f5c678d13f0c91b38905499f982d11..3b182f912fe9fdc2860389b072d67e4cffa6d5dd 100644
--- a/dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/2p2cnicouplinglocalresidual.hh
@@ -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);
diff --git a/dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh
index c1083413753a5da41aec9a00fb3fd9f25d6e40f2..68b59f9e4abf0ab3d8fc5c9d35a950cf2b197aef 100644
--- a/dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/boxcouplinglocalresidual.hh
@@ -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>
diff --git a/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
index d64e47a13529f17c88dd40ad161d61d3b88aee97..b4efd07c6677b26ae90946198e5353a6b7e6077c 100644
--- a/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/stokesnccouplinglocalresidual.hh
@@ -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;
- ////////////////////////////////////////////////////////////////////////////////////
- }
+ ////////////////////////////////////////////////////////////////////////////////////
}
}
diff --git a/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh b/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh
index 373546f19b10db8628ce38f8ba562902a051a7c1..dd9c3add7b3005122d098eb40f1f9ac181f8045a 100644
--- a/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh
+++ b/dumux/multidomain/couplinglocalresiduals/stokesncnicouplinglocalresidual.hh
@@ -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
diff --git a/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh b/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
index 1b698e0c8569373fc085266fa04c3028401aecd1..319a1b894283bd56054355604ef7ec7ea66d2812 100644
--- a/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
+++ b/test/multidomain/2cnistokes2p2cni/stokes2cnisubproblem.hh
@@ -357,10 +357,7 @@ public:
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
{
- if (onLowerBoundary_(globalPos))
- return alphaBJ_;
- else
- return 0.0;
+ return alphaBJ_;
}
/*!
diff --git a/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh b/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
index 2414382d0a09f3c17050cfa089f11343c57001b7..41232124c9fb68accefcab61dbd95fe82e036312 100644
--- a/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
+++ b/test/multidomain/2cstokes2p2c/stokes2csubproblem.hh
@@ -325,10 +325,7 @@ public:
*/
Scalar beaversJosephCoeffAtPos(const GlobalPosition &globalPos) const
{
- if (onLowerBoundary_(globalPos))
- return alphaBJ_;
- else
- return 0.0;
+ return alphaBJ_;
}
/*!