diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
index 9d1a9a9937c2690f932d6add1b8a964724a92dda..09bee3c3f5bb6255dd48276a95a59f01a702ff3f 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
@@ -657,29 +657,26 @@ void FVPressure2P<TypeTag>::getFlux(EntryType& entry, const Intersection& inters
density_[nPhaseIdx] = (potentialNW == 0) ? rhoMeanNW : density_[nPhaseIdx];
}
+ Scalar scalarPerm = permeability.two_norm();
//calculate current matrix entry
- entry[matrix] = (lambdaW + lambdaNW) * ((permeability * unitOuterNormal) / dist) * faceArea;
+ entry[matrix] = (lambdaW + lambdaNW) * scalarPerm / dist * faceArea;
//calculate right hand side
- entry[rhs] = (lambdaW * density_[wPhaseIdx] + lambdaNW * density_[nPhaseIdx]) * (permeability * gravity_) * faceArea;
+ //calculate unit distVec
+ distVec /= dist;
+ Scalar areaScaling = (unitOuterNormal * distVec);
+ //this treatment of g allows to account for gravity flux through faces where the face normal has no z component (e.g. parallelepiped grids)
+ entry[rhs] = (lambdaW * density_[wPhaseIdx] + lambdaNW * density_[nPhaseIdx]) * scalarPerm * (gravity_ * distVec) * faceArea * areaScaling;
if (pressureType_ == pw)
{
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
- pCGradient *= (pcI - pcJ) / dist;
-
//add capillary pressure term to right hand side
- entry[rhs] += 0.5 * (lambdaNWI + lambdaNWJ) * (permeability * pCGradient) * faceArea;
+ entry[rhs] += 0.5 * (lambdaNWI + lambdaNWJ) * scalarPerm * (pcI - pcJ) / dist * faceArea;
}
else if (pressureType_ == pn)
{
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
- pCGradient *= (pcI - pcJ) / dist;
-
//add capillary pressure term to right hand side
- entry[rhs] -= 0.5 * (lambdaWI + lambdaWJ) * (permeability * pCGradient) * faceArea;
+ entry[rhs] -= 0.5 * (lambdaWI + lambdaWJ) * scalarPerm * (pcI - pcJ) / dist * faceArea;
}
return;
@@ -888,31 +885,28 @@ const Intersection& intersection, const CellData& cellData, const bool first)
density_[nPhaseIdx] = (potentialNW == 0) ? rhoMeanNW : density_[nPhaseIdx];
}
+ Scalar scalarPerm = permeability.two_norm();
//calculate current matrix entry
- entry[matrix] = (lambdaW + lambdaNW) * ((permeability * unitOuterNormal) / dist) * faceArea;
+ entry[matrix] = (lambdaW + lambdaNW) * scalarPerm / dist * faceArea;
entry[rhs] = entry[matrix] * pressBound;
//calculate right hand side
- entry[rhs] -= (lambdaW * density_[wPhaseIdx] + lambdaNW * density_[nPhaseIdx]) * (permeability * gravity_)
- * faceArea;
+ //calculate unit distVec
+ distVec /= dist;
+ Scalar areaScaling = (unitOuterNormal * distVec);
+ //this treatment of g allows to account for gravity flux through faces where the face normal has no z component (e.g. parallelepiped grids)
+ entry[rhs] -= (lambdaW * density_[wPhaseIdx] + lambdaNW * density_[nPhaseIdx]) * scalarPerm * (gravity_ * distVec)
+ * faceArea * areaScaling;
if (pressureType_ == pw)
{
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
- pCGradient *= (pcI - pcBound) / dist;
-
//add capillary pressure term to right hand side
- entry[rhs] -= 0.5 * (lambdaNWI + lambdaNWBound) * (permeability * pCGradient) * faceArea;
+ entry[rhs] -= 0.5 * (lambdaNWI + lambdaNWBound) * scalarPerm * (pcI - pcBound) / dist * faceArea;
}
else if (pressureType_ == pn)
{
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
- pCGradient *= (pcI - pcBound) / dist;
-
//add capillary pressure term to right hand side
- entry[rhs] += 0.5 * (lambdaWI + lambdaWBound) * (permeability * pCGradient) * faceArea;
+ entry[rhs] += 0.5 * (lambdaWI + lambdaWBound) * scalarPerm * (pcI - pcBound) / dist * faceArea;
}
}
//set neumann boundary condition
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
index c593168b945e1da3b432b1806dd47167da972459..5a1cd7804d6b152ca7e3f6ff8bdc4633d03ff7f2 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
@@ -390,41 +390,40 @@ void FVVelocity2P<TypeTag>::calculateVelocity(const Intersection& intersection,
density_[nPhaseIdx];
}
+ Scalar scalarPerm = permeability.two_norm();
+
//calculate the gravity term
- Dune::FieldVector < Scalar, dimWorld > velocityW(permeability);
- Dune::FieldVector < Scalar, dimWorld > velocityNW(permeability);
- Dune::FieldVector < Scalar, dimWorld > gravityTermW(unitOuterNormal);
- Dune::FieldVector < Scalar, dimWorld > gravityTermNW(unitOuterNormal);
+ Dune::FieldVector < Scalar, dimWorld > velocityW(unitOuterNormal);
+ Dune::FieldVector < Scalar, dimWorld > velocityNW(unitOuterNormal);
- gravityTermW *= (gravity_ * permeability) * (lambdaW * density_[wPhaseIdx]);
- gravityTermNW *= (gravity_ * permeability) * (lambdaNW * density_[nPhaseIdx]);
+ //calculate unit distVec
+ distVec /= dist;
+ Scalar areaScaling = (unitOuterNormal * distVec);
+ //this treatment of g allows to account for gravity flux through faces where the face normal has no z component (e.g. parallelepiped grids)
+ Scalar gravityTermW = (gravity_* distVec) * density_[wPhaseIdx] * areaScaling;
+ Scalar gravityTermNW = (gravity_ * distVec) * density_[nPhaseIdx] * areaScaling;
//calculate velocity depending on the pressure used -> use pc = pn - pw
switch (pressureType_)
{
case pw:
{
- velocityW *= lambdaW * (cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist;
- velocityNW *= lambdaNW * (cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist
- + 0.5 * (lambdaNWI + lambdaNWJ) * (pcI - pcJ) / dist;
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
+ velocityW *= lambdaW * scalarPerm * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermW);
+ velocityNW *= lambdaNW * scalarPerm * ((cellData.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx)) / dist + gravityTermNW)
+ + 0.5 * (lambdaNWI + lambdaNWJ) * scalarPerm * (pcI - pcJ) / dist;
break;
}
case pn:
{
- velocityW *= lambdaW * (cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist
- - 0.5 * (lambdaWI + lambdaWJ) * (pcI - pcJ) / dist;
- velocityNW *= lambdaNW * (cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist;
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
+ velocityW *= lambdaW * scalarPerm * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermW)
+ - 0.5 * (lambdaWI + lambdaWJ) * scalarPerm * (pcI - pcJ) / dist;
+ velocityNW *= lambdaNW * scalarPerm * ((cellData.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx)) / dist + gravityTermNW);
break;
}
case pglobal:
{
- velocityW *= (lambdaW + lambdaNW) * (cellData.globalPressure() - cellDataJ.globalPressure()) / dist;
- velocityW += gravityTermW;
- velocityW += gravityTermNW;
+ velocityW *= (lambdaW + lambdaNW) * scalarPerm * (cellData.globalPressure() - cellDataJ.globalPressure()) / dist +
+ scalarPerm * (lambdaW * gravityTermW + lambdaNW * gravityTermNW);
velocityNW = 0;
break;
}
@@ -619,41 +618,40 @@ void FVVelocity2P<TypeTag>::calculateVelocityOnBoundary(const Intersection& inte
density_[nPhaseIdx] = (potentialNW == 0) ? 0.5 * (cellData.density(nPhaseIdx) + densityNWBound) : density_[nPhaseIdx];
}
+ Scalar scalarPerm = permeability.two_norm();
+
//calculate the gravity term
- Dune::FieldVector < Scalar, dimWorld > velocityW(permeability);
- Dune::FieldVector < Scalar, dimWorld > velocityNW(permeability);
- Dune::FieldVector < Scalar, dimWorld > gravityTermW(unitOuterNormal);
- Dune::FieldVector < Scalar, dimWorld > gravityTermNW(unitOuterNormal);
+ Dune::FieldVector < Scalar, dimWorld > velocityW(unitOuterNormal);
+ Dune::FieldVector < Scalar, dimWorld > velocityNW(unitOuterNormal);
- gravityTermW *= (gravity_ * permeability) * (lambdaW * density_[wPhaseIdx]);
- gravityTermNW *= (gravity_ * permeability) * (lambdaNW * density_[nPhaseIdx]);
+ //calculate unit distVec
+ distVec /= dist;
+ Scalar areaScaling = (unitOuterNormal * distVec);
+ //this treatment of g allows to account for gravity flux through faces where the face normal has no z component (e.g. parallelepiped grids)
+ Scalar gravityTermW = (gravity_ * distVec) * density_[wPhaseIdx] * areaScaling;
+ Scalar gravityTermNW = (gravity_ * distVec) * density_[nPhaseIdx] * areaScaling;
//calculate velocity depending on the pressure used -> use pc = pn - pw
switch (pressureType_)
{
case pw:
{
- velocityW *= lambdaW * (cellData.pressure(wPhaseIdx) - pressBound) / dist;
- velocityNW *= lambdaNW * (cellData.pressure(wPhaseIdx) - pressBound) / dist
- + 0.5 * (lambdaNWI + lambdaNWBound) * (pcI - pcBound) / dist;
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
+ velocityW *= lambdaW * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermW);
+ velocityNW *= lambdaNW * scalarPerm * ((cellData.pressure(wPhaseIdx) - pressBound) / dist + gravityTermNW)
+ + 0.5 * (lambdaNWI + lambdaNWBound) * scalarPerm * (pcI - pcBound) / dist;
break;
}
case pn:
{
- velocityW *= lambdaW * (cellData.pressure(nPhaseIdx) - pressBound) / dist
- - 0.5 * (lambdaWI + lambdaWBound) * (pcI - pcBound) / dist;
- velocityNW *= lambdaNW * (cellData.pressure(nPhaseIdx) - pressBound) / dist;
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
+ velocityW *= lambdaW * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermW)
+ - 0.5 * (lambdaWI + lambdaWBound) * scalarPerm * (pcI - pcBound) / dist;
+ velocityNW *= lambdaNW * scalarPerm * ((cellData.pressure(nPhaseIdx) - pressBound) / dist + gravityTermNW);
break;
}
case pglobal:
{
- velocityW *= (lambdaW + lambdaNW) * (cellData.globalPressure() - pressBound) / dist;
- velocityW += gravityTermW;
- velocityW += gravityTermNW;
+ velocityW *= (lambdaW + lambdaNW) * scalarPerm * (cellData.globalPressure() - pressBound) / dist +
+ scalarPerm * (lambdaW * gravityTermW + lambdaNW * gravityTermNW);
velocityNW = 0;
break;
}