Commit a6e9a635 authored by Markus Wolff's avatar Markus Wolff
Browse files

Implementation of naming conventions



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@11554 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent b814e29a
......@@ -253,19 +253,19 @@ public:
}
case pn:
{
Scalar potNW = this->pressure()[globalIdx];
Scalar potNw = this->pressure()[globalIdx];
Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
Scalar potPc = cellData.capillaryPressure()
+ gravityDiff * (density_[nPhaseIdx] - density_[wPhaseIdx]);
cellData.setPotential(nPhaseIdx, potNW);
cellData.setPotential(wPhaseIdx, potNW - potPc);
cellData.setPotential(nPhaseIdx, potNw);
cellData.setPotential(wPhaseIdx, potNw - potPc);
Scalar pressNW = potNW - gravityDiff * density_[nPhaseIdx];
Scalar pressNw = potNw - gravityDiff * density_[nPhaseIdx];
cellData.setPressure(wPhaseIdx, pressNW - cellData.capillaryPressure());
cellData.setPressure(nPhaseIdx, pressNW);
cellData.setPressure(wPhaseIdx, pressNw - cellData.capillaryPressure());
cellData.setPressure(nPhaseIdx, pressNw);
break;
}
......@@ -1870,27 +1870,27 @@ void FvMpfaL2dPressure2p<TypeTag>::assemble()
//calculate potential gradients
Scalar potentialDiffW = 0;
Scalar potentialDiffNW = 0;
Scalar potentialDiffNw = 0;
switch (pressureType_)
{
case pw:
{
potentialBound += density_[wPhaseIdx]*gdeltaZ;
potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound) / dist;
potentialDiffNW = (cellData.potential(nPhaseIdx) - potentialBound - pcBound) / dist;
potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound - pcBound) / dist;
break;
}
case pn:
{
potentialBound += density_[nPhaseIdx]*gdeltaZ;
potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound + pcBound) / dist;
potentialDiffNW = (cellData.potential(nPhaseIdx) - potentialBound) / dist;
potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound) / dist;
break;
}
}
Scalar lambdaTotal = (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
lambdaTotal += (potentialDiffNW >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
lambdaTotal += (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
DimVector permTimesNormal(0);
permeability.mv(unitDistVec, permTimesNormal);
......@@ -2345,16 +2345,16 @@ void FvMpfaL2dPressure2p<TypeTag>::updateMaterialLaws()
// initialize mobilities
Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*eIt), satW)
/ viscosity_[wPhaseIdx];
Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), satW)
Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), satW)
/ viscosity_[nPhaseIdx];
// initialize mobilities
cellData.setMobility(wPhaseIdx, mobilityW);
cellData.setMobility(nPhaseIdx, mobilityNW);
cellData.setMobility(nPhaseIdx, mobilityNw);
//initialize fractional flow functions
cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNW));
cellData.setFracFlowFunc(nPhaseIdx, mobilityNW / (mobilityW + mobilityNW));
cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNw));
cellData.setFracFlowFunc(nPhaseIdx, mobilityNw / (mobilityW + mobilityNw));
}
return;
}
......
......@@ -264,19 +264,19 @@ public:
}
case pn:
{
Scalar potNW = this->pressure()[globalIdx];
Scalar potNw = this->pressure()[globalIdx];
Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
Scalar potPc = cellData.capillaryPressure()
+ gravityDiff * (density_[nPhaseIdx] - density_[wPhaseIdx]);
cellData.setPotential(nPhaseIdx, potNW);
cellData.setPotential(wPhaseIdx, potNW - potPc);
cellData.setPotential(nPhaseIdx, potNw);
cellData.setPotential(wPhaseIdx, potNw - potPc);
Scalar pressNW = potNW - gravityDiff * density_[nPhaseIdx];
Scalar pressNw = potNw - gravityDiff * density_[nPhaseIdx];
cellData.setPressure(wPhaseIdx, pressNW - cellData.capillaryPressure());
cellData.setPressure(nPhaseIdx, pressNW);
cellData.setPressure(wPhaseIdx, pressNw - cellData.capillaryPressure());
cellData.setPressure(nPhaseIdx, pressNw);
break;
}
......@@ -2686,27 +2686,27 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::assemble()
//calculate potential gradients
Scalar potentialDiffW = 0;
Scalar potentialDiffNW = 0;
Scalar potentialDiffNw = 0;
switch (pressureType_)
{
case pw:
{
potentialBound += density_[wPhaseIdx]*gdeltaZ;
potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound) / dist;
potentialDiffNW = (cellData.potential(nPhaseIdx) - potentialBound - pcBound) / dist;
potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound - pcBound) / dist;
break;
}
case pn:
{
potentialBound += density_[nPhaseIdx]*gdeltaZ;
potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBound + pcBound) / dist;
potentialDiffNW = (cellData.potential(nPhaseIdx) - potentialBound) / dist;
potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBound) / dist;
break;
}
}
Scalar lambdaTotal = (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
lambdaTotal += (potentialDiffNW >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
lambdaTotal += (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
DimVector permTimesNormal(0);
permeability.mv(unitDistVec, permTimesNormal);
......@@ -3392,16 +3392,16 @@ void FvMpfaL2dPressure2pAdaptive<TypeTag>::updateMaterialLaws()
// initialize mobilities
Scalar mobilityW = MaterialLaw::krw(problem_.spatialParams().materialLawParams(*eIt), satW)
/ viscosity_[wPhaseIdx];
Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), satW)
Scalar mobilityNw = MaterialLaw::krn(problem_.spatialParams().materialLawParams(*eIt), satW)
/ viscosity_[nPhaseIdx];
// initialize mobilities
cellData.setMobility(wPhaseIdx, mobilityW);
cellData.setMobility(nPhaseIdx, mobilityNW);
cellData.setMobility(nPhaseIdx, mobilityNw);
//initialize fractional flow functions
cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNW));
cellData.setFracFlowFunc(nPhaseIdx, mobilityNW / (mobilityW + mobilityNW));
cellData.setFracFlowFunc(wPhaseIdx, mobilityW / (mobilityW + mobilityNw));
cellData.setFracFlowFunc(nPhaseIdx, mobilityNw / (mobilityW + mobilityNw));
}
return;
}
......
......@@ -236,7 +236,7 @@ public:
CellData& cellData = problem_.variables().cellData(globalIdx);
Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
Dune::FieldVector < Scalar, 2 * dim > fluxNW(0);
Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
// run through all intersections with neighbors and boundary
IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
......@@ -245,7 +245,7 @@ public:
fluxW[isIndex] = isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(wPhaseIdx, isIndex));
fluxNW[isIndex] =
fluxNw[isIndex] =
isIt->geometry().volume()
* (isIt->centerUnitOuterNormal()
* cellData.fluxData().velocity(nPhaseIdx, isIndex));
......@@ -270,7 +270,7 @@ public:
refVelocity = 0;
for (int i = 0; i < dim; i++)
refVelocity[i] = 0.5 * (fluxNW[2*i + 1] - fluxNW[2*i]);
refVelocity[i] = 0.5 * (fluxNw[2*i + 1] - fluxNw[2*i]);
// calculate the element velocity by the Piola transformation
elementVelocity = 0;
......@@ -340,17 +340,17 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
// get pressure values
Dune::FieldVector < Scalar, 2 * dim > potW(0);
Dune::FieldVector < Scalar, 2 * dim > potNW(0);
Dune::FieldVector < Scalar, 2 * dim > potNw(0);
potW[0] = cellData1.potential(wPhaseIdx);
potW[1] = cellData2.potential(wPhaseIdx);
potW[2] = cellData3.potential(wPhaseIdx);
potW[3] = cellData4.potential(wPhaseIdx);
potNW[0] = cellData1.potential(nPhaseIdx);
potNW[1] = cellData2.potential(nPhaseIdx);
potNW[2] = cellData3.potential(nPhaseIdx);
potNW[3] = cellData4.potential(nPhaseIdx);
potNw[0] = cellData1.potential(nPhaseIdx);
potNw[1] = cellData2.potential(nPhaseIdx);
potNw[2] = cellData3.potential(nPhaseIdx);
potNw[3] = cellData4.potential(nPhaseIdx);
//get mobilities of the phases
Dune::FieldVector<Scalar, numPhases> lambda1(cellData1.mobility(wPhaseIdx));
......@@ -396,14 +396,14 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
Scalar potentialDiffW32 = 0;
Scalar potentialDiffW34 = 0;
Scalar potentialDiffNW12 = 0;
Scalar potentialDiffNW14 = 0;
Scalar potentialDiffNW32 = 0;
Scalar potentialDiffNW34 = 0;
Scalar potentialDiffNw12 = 0;
Scalar potentialDiffNw14 = 0;
Scalar potentialDiffNw32 = 0;
Scalar potentialDiffNw34 = 0;
//flux vector
Dune::FieldVector < Scalar, 2 * dim > fluxW(0);
Dune::FieldVector < Scalar, 2 * dim > fluxNW(0);
Dune::FieldVector < Scalar, 2 * dim > fluxNw(0);
Dune::FieldMatrix < Scalar, dim, 2 * dim - dim + 1 > T(0);
DimVector Tu(0);
......@@ -422,14 +422,14 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
fluxW[0] = Tu[1];
potentialDiffW12 = Tu[1];
u[0] = potNW[1];
u[1] = potNW[2];
u[2] = potNW[0];
u[0] = potNw[1];
u[1] = potNw[2];
u[2] = potNw[0];
T.mv(u, Tu);
fluxNW[0] = Tu[1];
potentialDiffNW12 = Tu[1];
fluxNw[0] = Tu[1];
potentialDiffNw12 = Tu[1];
}
else
{
......@@ -442,14 +442,14 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
fluxW[0] = Tu[1];
potentialDiffW12 = Tu[1];
u[0] = potNW[0];
u[1] = potNW[3];
u[2] = potNW[1];
u[0] = potNw[0];
u[1] = potNw[3];
u[2] = potNw[1];
T.mv(u, Tu);
fluxNW[0] = Tu[1];
potentialDiffNW12 = Tu[1];
fluxNw[0] = Tu[1];
potentialDiffNw12 = Tu[1];
}
rightTriangle = this->calculateTransmissibility(T, interactionVolume, lambda, 1, 2, 3, 0);
......@@ -465,14 +465,14 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
fluxW[1] = Tu[1];
potentialDiffW32 = -Tu[1];
u[0] = potNW[2];
u[1] = potNW[3];
u[2] = potNW[1];
u[0] = potNw[2];
u[1] = potNw[3];
u[2] = potNw[1];
T.mv(u, Tu);
fluxNW[1] = Tu[1];
potentialDiffNW32 = -Tu[1];
fluxNw[1] = Tu[1];
potentialDiffNw32 = -Tu[1];
}
else
{
......@@ -485,14 +485,14 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
fluxW[1] = Tu[1];
potentialDiffW32 = -Tu[1];
u[0] = potNW[1];
u[1] = potNW[0];
u[2] = potNW[2];
u[0] = potNw[1];
u[1] = potNw[0];
u[2] = potNw[2];
T.mv(u, Tu);
fluxNW[1] = Tu[1];
potentialDiffNW32 = -Tu[1];
fluxNw[1] = Tu[1];
potentialDiffNw32 = -Tu[1];
}
rightTriangle = this->calculateTransmissibility(T, interactionVolume, lambda, 2, 3, 0, 1);
......@@ -508,14 +508,14 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
fluxW[2] = Tu[1];
potentialDiffW34 = Tu[1];
u[0] = potNW[3];
u[1] = potNW[0];
u[2] = potNW[2];
u[0] = potNw[3];
u[1] = potNw[0];
u[2] = potNw[2];
T.mv(u, Tu);
fluxNW[2] = Tu[1];
potentialDiffNW34 = Tu[1];
fluxNw[2] = Tu[1];
potentialDiffNw34 = Tu[1];
}
else
{
......@@ -528,14 +528,14 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
fluxW[2] = Tu[1];
potentialDiffW34 = Tu[1];
u[0] = potNW[2];
u[1] = potNW[1];
u[2] = potNW[3];
u[0] = potNw[2];
u[1] = potNw[1];
u[2] = potNw[3];
T.mv(u, Tu);
fluxNW[2] = Tu[1];
potentialDiffNW34 = Tu[1];
fluxNw[2] = Tu[1];
potentialDiffNw34 = Tu[1];
}
rightTriangle = this->calculateTransmissibility(T, interactionVolume, lambda, 3, 0, 1, 2);
......@@ -551,14 +551,14 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
fluxW[3] = Tu[1];
potentialDiffW14 = -Tu[1];
u[0] = potNW[0];
u[1] = potNW[1];
u[2] = potNW[3];
u[0] = potNw[0];
u[1] = potNw[1];
u[2] = potNw[3];
T.mv(u, Tu);
fluxNW[3] = Tu[1];
potentialDiffNW14 = -Tu[1];
fluxNw[3] = Tu[1];
potentialDiffNw14 = -Tu[1];
}
else
{
......@@ -571,53 +571,53 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
fluxW[3] = Tu[1];
potentialDiffW14 = -Tu[1];
u[0] = potNW[3];
u[1] = potNW[2];
u[2] = potNW[0];
u[0] = potNw[3];
u[1] = potNw[2];
u[2] = potNw[0];
T.mv(u, Tu);
fluxNW[3] = Tu[1];
potentialDiffNW14 = -Tu[1];
fluxNw[3] = Tu[1];
potentialDiffNw14 = -Tu[1];
}
//store potentials for further calculations (saturation, ...)
cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffW12);
cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffNW12);
cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 0), potentialDiffNw12);
cellData1.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(0, 1), potentialDiffW14);
cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 1), potentialDiffNW14);
cellData1.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(0, 1), potentialDiffNw14);
cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 0), -potentialDiffW32);
cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 0), -potentialDiffNW32);
cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 0), -potentialDiffNw32);
cellData2.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffW12);
cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffNW12);
cellData2.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(1, 1), -potentialDiffNw12);
cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 0), potentialDiffW34);
cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 0), potentialDiffNW34);
cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 0), potentialDiffNw34);
cellData3.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(2, 1), potentialDiffW32);
cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 1), potentialDiffNW32);
cellData3.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(2, 1), potentialDiffNw32);
cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -potentialDiffW14);
cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -potentialDiffNW14);
cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 0), -potentialDiffNw14);
cellData4.fluxData().addUpwindPotential(wPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffW34);
cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffNW34);
cellData4.fluxData().addUpwindPotential(nPhaseIdx, interactionVolume.getIndexOnElement(3, 1), -potentialDiffNw34);
//compute mobilities of face 1
Dune::FieldVector<Scalar, numPhases> lambda12Upw(0.0);
lambda12Upw[wPhaseIdx] = (potentialDiffW12 >= 0) ? lambda1[wPhaseIdx] : lambda2[wPhaseIdx];
lambda12Upw[nPhaseIdx] = (potentialDiffNW12 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
lambda12Upw[nPhaseIdx] = (potentialDiffNw12 >= 0) ? lambda1[nPhaseIdx] : lambda2[nPhaseIdx];
//compute mobilities of face 4
Dune::FieldVector<Scalar, numPhases> lambda14Upw(0.0);
lambda14Upw[wPhaseIdx] = (potentialDiffW14 >= 0) ? lambda1[wPhaseIdx] : lambda4[wPhaseIdx];
lambda14Upw[nPhaseIdx] = (potentialDiffNW14 >= 0) ? lambda1[nPhaseIdx] : lambda4[nPhaseIdx];
lambda14Upw[nPhaseIdx] = (potentialDiffNw14 >= 0) ? lambda1[nPhaseIdx] : lambda4[nPhaseIdx];
//compute mobilities of face 2
Dune::FieldVector<Scalar, numPhases> lambda32Upw(0.0);
lambda32Upw[wPhaseIdx] = (potentialDiffW32 >= 0) ? lambda3[wPhaseIdx] : lambda2[wPhaseIdx];
lambda32Upw[nPhaseIdx] = (potentialDiffNW32 >= 0) ? lambda3[nPhaseIdx] : lambda2[nPhaseIdx];
lambda32Upw[nPhaseIdx] = (potentialDiffNw32 >= 0) ? lambda3[nPhaseIdx] : lambda2[nPhaseIdx];
//compute mobilities of face 3
Dune::FieldVector<Scalar, numPhases> lambda34Upw(0.0);
lambda34Upw[wPhaseIdx] = (potentialDiffW34 >= 0) ? lambda3[wPhaseIdx] : lambda4[wPhaseIdx];
lambda34Upw[nPhaseIdx] = (potentialDiffNW34 >= 0) ? lambda3[nPhaseIdx] : lambda4[nPhaseIdx];
lambda34Upw[nPhaseIdx] = (potentialDiffNw34 >= 0) ? lambda3[nPhaseIdx] : lambda4[nPhaseIdx];
for (int i = 0; i < numPhases; i++)
{
......@@ -641,7 +641,7 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
}
case nPhaseIdx:
{
flux = fluxNW;
flux = fluxNw;
break;
}
}
......@@ -827,14 +827,14 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
Scalar gdeltaZ = (problem_.bBoxMax()-globalPosFace) * gravity_;
Scalar potentialBoundW = interactionVolume.getDirichletValues(intVolFaceIdx)[pressureIdx] + density_[wPhaseIdx]*gdeltaZ;
Scalar potentialBoundNW = potentialBoundW;
Scalar potentialBoundNw = potentialBoundW;
//calculate potential gradients
switch (pressureType_)
{
case pw:
{
potentialBoundNW += pcBound;
potentialBoundNw += pcBound;
break;
}
case pn:
......@@ -846,15 +846,15 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
}
Scalar potentialDiffW = (cellData.potential(wPhaseIdx) - potentialBoundW) / dist;
Scalar potentialDiffNW = (cellData.potential(nPhaseIdx) - potentialBoundNW) / dist;
Scalar potentialDiffNw = (cellData.potential(nPhaseIdx) - potentialBoundNw) / dist;
//store potentials for further calculations (saturation, ...)
cellData.fluxData().addUpwindPotential(wPhaseIdx, boundaryFaceIdx, potentialDiffW);
cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, potentialDiffNW);
cellData.fluxData().addUpwindPotential(nPhaseIdx, boundaryFaceIdx, potentialDiffNw);
//calculated phase velocities from advective velocities -> capillary pressure velocity already added in pressure part!
DimVector velocityW(0);
DimVector velocityNW(0);
DimVector velocityNw(0);
// calculate capillary pressure gradient
DimVector pressGradient = unitDistVec;
......@@ -862,21 +862,21 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
permeability.mv(pressGradient, velocityW);
pressGradient = unitDistVec;
pressGradient *= (cellData.potential(nPhaseIdx) - potentialBoundNW) / dist;
permeability.mv(pressGradient, velocityNW);
pressGradient *= (cellData.potential(nPhaseIdx) - potentialBoundNw) / dist;
permeability.mv(pressGradient, velocityNw);
velocityW *= (potentialDiffW >= 0.) ? lambda[wPhaseIdx] : lambdaBound[wPhaseIdx];
velocityNW *= (potentialDiffNW >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
velocityNw *= (potentialDiffNw >= 0.) ? lambda[nPhaseIdx] : lambdaBound[nPhaseIdx];
//velocity is calculated from two vertices of one intersection!
velocityW *= 0.5;
velocityNW *= 0.5;
velocityNw *= 0.5;
//store velocities
velocityW += cellData.fluxData().velocity(wPhaseIdx, boundaryFaceIdx);
velocityNW += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
velocityNw += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
cellData.fluxData().setVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNW);
cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
}
......@@ -903,10 +903,10 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
boundValues[nPhaseIdx] /= density_[nPhaseIdx];
DimVector velocityW(unitDistVec);
DimVector velocityNW(unitDistVec);
DimVector velocityNw(unitDistVec);
velocityW *= boundValues[wPhaseIdx] / (2 * interactionVolume.getFaceArea(elemIdx, faceIdx));
velocityNW *= boundValues[nPhaseIdx]
velocityNw *= boundValues[nPhaseIdx]
/ (2 * interactionVolume.getFaceArea(elemIdx, faceIdx));
//store potentials for further calculations (saturation, ...)
......@@ -915,9 +915,9 @@ void FvMpfaL2dPressureVelocity2p<TypeTag>::calculateVelocity()
//store velocities
velocityW += cellData.fluxData().velocity(wPhaseIdx, boundaryFaceIdx);
velocityNW += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
velocityNw += cellData.fluxData().velocity(nPhaseIdx, boundaryFaceIdx);
cellData.fluxData().setVelocity(wPhaseIdx, boundaryFaceIdx, velocityW);
cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNW);
cellData.fluxData().setVelocity(nPhaseIdx, boundaryFaceIdx, velocityNw);
cellData.fluxData().setVelocityMarker(boundaryFaceIdx);
}
else
......
......@@ -284,19 +284,19 @@ public:
}
case pn:
{
Scalar potNW = this->pressure()[globalIdx];
Scalar potNw = this->pressure()[globalIdx];
Scalar gravityDiff = (problem_.bBoxMax() - element.geometry().center()) * gravity_;
Scalar potPc = cellData.capillaryPressure()
+ gravityDiff * (density_[nPhaseIdx] - density_[wPhaseIdx]);
cellData.setPotential(nPhaseIdx, potNW);
cellData.setPotential(wPhaseIdx, potNW - potPc);
cellData.setPotential(nPhaseIdx, potNw);
cellData.setPotential(wPhaseIdx, potNw - potPc);
Scalar pressNW = potNW - gravityDiff * density_[nPhaseIdx];
Scalar pressNw = potNw - gravityDiff * density_[nPhaseIdx];