diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
index 12265729e9f1d96e2cd2b7d643076475c39869cc..ea0d8d1242eb95229a7927aa1660c172c09cd4ce 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
@@ -758,7 +758,7 @@ void FVPressure2P<TypeTag>::assemble(bool first)
else
{
- std::vector<Scalar> J(problem_.neumannPress(globalPosFace, *isIt));
+ std::vector<Scalar> J(problem_.neumann(globalPosFace, *isIt));
if (!compressibility)
{
J[wPhaseIdx] /= densityWI;
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
index b37922e7336d0b8ee4177b33c64a21fd077303f4..ac7b4498721b461ba43d75dfe6431e90ab9b9889 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh
@@ -656,7 +656,7 @@ void FVVelocity2P<TypeTag>::calculateVelocity()
else
{
- std::vector<Scalar> J = this->problem().neumannPress(globalPosFace, *isIt);
+ std::vector<Scalar> J = this->problem().neumann(globalPosFace, *isIt);
Dune::FieldVector<Scalar,dimWorld> velocityW(unitOuterNormal);
Dune::FieldVector<Scalar,dimWorld> velocityNW(unitOuterNormal);
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh
index 1ce4f8feb96240345f0b1e43913da1a96553afb4..18bb03bbf928407f64d72ea647596e7eb1c392be 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaopressure2p.hh
@@ -951,7 +951,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble()
if (nextisItbctype == BoundaryConditions::neumann)
{
// get Neumann boundary value of 'nextisIt'
- std::vector<Scalar> J(problem_.neumannPress(globalPosFace13, *nextisIt));
+ std::vector<Scalar> J(problem_.neumann(globalPosFace13, *nextisIt));
double J3 = (J[wPhaseIdx]/densityW+J[nPhaseIdx]/densityNW);
// get boundary condition for boundary face (isIt24) center
@@ -962,7 +962,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble()
if (isIt24bctype == BoundaryConditions::neumann)
{
// get neumann boundary value of 'isIt24'
- std::vector<Scalar> J(problem_.neumannPress(globalPosFace24, *isIt24));
+ std::vector<Scalar> J(problem_.neumann(globalPosFace24, *isIt24));
double J4 = (J[wPhaseIdx]/densityW+J[nPhaseIdx]/densityNW);
// compute normal vectors nu11,nu21; nu12, nu22;
@@ -1224,7 +1224,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble()
if (isIt24bctype == BoundaryConditions::neumann)
{
// get Neumann boundary value of 'isIt24'
- std::vector<Scalar> J(problem_.neumannPress(globalPosFace24, *isIt24));
+ std::vector<Scalar> J(problem_.neumann(globalPosFace24, *isIt24));
double J4 = (J[wPhaseIdx]/densityW+J[nPhaseIdx]/densityNW);
// compute normal vectors nu11,nu21; nu12, nu22;
@@ -1423,7 +1423,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble()
if (isItbctype == BoundaryConditions::neumann)
{
// get Neumann boundary value
- std::vector<Scalar> J(problem_.neumannPress(globalPosFace12, *isIt));
+ std::vector<Scalar> J(problem_.neumann(globalPosFace12, *isIt));
double J1 = (J[wPhaseIdx]/densityW+J[nPhaseIdx]/densityNW);
// evaluate right hand side
@@ -1583,7 +1583,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble()
if (isIt34bctype == BoundaryConditions::neumann)
{
// get Neumann boundary value
- std::vector<Scalar> J(problem_.neumannPress(globalPosFace34, *isIt34));
+ std::vector<Scalar> J(problem_.neumann(globalPosFace34, *isIt34));
double J2 = (J[wPhaseIdx]/densityW+J[nPhaseIdx]/densityNW);
// compute normal vectors nu11,nu21; nu13, nu23;
@@ -1963,7 +1963,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble()
else
{
// get Neumann boundary value of 'nextisIt'
- std::vector<Scalar> J(problem_.neumannPress(globalPosFace13, *nextisIt));
+ std::vector<Scalar> J(problem_.neumann(globalPosFace13, *nextisIt));
double J3 = (J[wPhaseIdx]/densityW+J[nPhaseIdx]/densityNW);
// compute normal vectors nu11,nu21;
@@ -2177,7 +2177,7 @@ void FVMPFAOPressure2P<TypeTag>::assemble()
else
{
// get Neumann boundary value of 'isIt34'
- std::vector<Scalar> J(problem_.neumannPress(globalPosFace34, *isIt34));
+ std::vector<Scalar> J(problem_.neumann(globalPosFace34, *isIt34));
double J2 = (J[wPhaseIdx]/densityW+J[nPhaseIdx]/densityNW);
// compute normal vectors nu11,nu21; nu13, nu23;
diff --git a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh
index 824ed06e8cdee4d830d0e327c50743966cdaead9..5ad2b8b703ebbdbf30748e684726d80ff19b39dc 100644
--- a/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh
+++ b/dumux/decoupled/2p/diffusion/fvmpfa/fvmpfaovelocity2p.hh
@@ -638,7 +638,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity()
if (nextisItbctype == BoundaryConditions::neumann)
{
// get Neumann boundary value of 'nextisIt'
- std::vector<Scalar> J(this->problem().neumannPress(globalPosFace13, *nextisIt));
+ std::vector<Scalar> J(this->problem().neumann(globalPosFace13, *nextisIt));
double J3 = (J[wPhaseIdx]/densityW + J[nPhaseIdx]/densityNW);
// get boundary condition for boundary face (isIt24) center
@@ -648,7 +648,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity()
if (isIt24bctype == BoundaryConditions::neumann)
{
// get neumann boundary value of 'isIt24'
- std::vector<Scalar> J(this->problem().neumannPress(globalPosFace24, *isIt24));
+ std::vector<Scalar> J(this->problem().neumann(globalPosFace24, *isIt24));
double J4 = (J[wPhaseIdx]/densityW + J[nPhaseIdx]/densityNW);
// compute normal vectors nu11,nu21; nu12, nu22;
@@ -914,7 +914,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity()
if (isIt24bctype == BoundaryConditions::neumann)
{
// get Neumann boundary value of 'isIt24'
- std::vector<Scalar> J(this->problem().neumannPress(globalPosFace24, *isIt24));
+ std::vector<Scalar> J(this->problem().neumann(globalPosFace24, *isIt24));
double J4 = (J[wPhaseIdx]/densityW + J[nPhaseIdx]/densityNW);
// compute normal vectors nu11,nu21; nu12, nu22;
@@ -1132,7 +1132,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity()
if (isItbctype == BoundaryConditions::neumann)
{
// get Neumann boundary value
- std::vector<Scalar> J(this->problem().neumannPress(globalPosFace12, *isIt));
+ std::vector<Scalar> J(this->problem().neumann(globalPosFace12, *isIt));
double J1 = (J[wPhaseIdx]/densityW + J[nPhaseIdx]/densityNW);
// evaluate velocity of facet 'isIt'
@@ -1298,7 +1298,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity()
if (isIt34bctype == BoundaryConditions::neumann)
{
// get Neumann boundary value
- std::vector<Scalar> J(this->problem().neumannPress(globalPosFace34, *isIt34));
+ std::vector<Scalar> J(this->problem().neumann(globalPosFace34, *isIt34));
double J2 = (J[wPhaseIdx]/densityW + J[nPhaseIdx]/densityNW);
// compute normal vectors nu11,nu21; nu13, nu23;
@@ -1692,7 +1692,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity()
else
{
// get Neumann boundary value of 'nextisIt'
- std::vector<Scalar> J(this->problem().neumannPress(globalPosFace13, *nextisIt));
+ std::vector<Scalar> J(this->problem().neumann(globalPosFace13, *nextisIt));
double J3 = (J[wPhaseIdx]/densityW + J[nPhaseIdx]/densityNW);
// compute normal vectors nu11,nu21;
@@ -1919,7 +1919,7 @@ void FVMPFAOVelocity2P<TypeTag>::calculateVelocity()
else
{
// get Neumann boundary value of 'isIt34'
- std::vector<Scalar> J(this->problem().neumannPress(globalPosFace34, *isIt34));
+ std::vector<Scalar> J(this->problem().neumann(globalPosFace34, *isIt34));
double J2 = (J[wPhaseIdx]/densityW + J[nPhaseIdx]/densityNW);
// compute normal vectors nu11,nu21; nu13, nu23;
diff --git a/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh b/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh
index 83b901d524e65df65e9b059fc8f47efb08844306..f54c93ab9722e5506f12ff11da267b70138ad475 100644
--- a/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh
+++ b/dumux/decoupled/2p/diffusion/mimetic/mimeticgroundwater.hh
@@ -382,7 +382,7 @@ private:
if (bctypeface == BoundaryConditions::neumann)
{
- std::vector<Scalar> neumannFlux(problem_.neumannPress(faceGlobal, *it));
+ std::vector<Scalar> neumannFlux(problem_.neumann(faceGlobal, *it));
Scalar J = (neumannFlux[wetting]+neumannFlux[nonwetting]);
this->b[faceIndex] -= J*it->geometry().volume();
}
diff --git a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
index 24f0c1e0d5e2fd148e866a02d527e5fbd89bbdc4..888fbcae44ff9d2bb70562077d84dadd461c3750 100644
--- a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
+++ b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
@@ -61,14 +61,19 @@ class FVSaturation2P
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Variables)) Variables;
- typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElementContainer;
- typedef Dune::GenericReferenceElements<Scalar, dim - 1> ReferenceElementFaceContainer;
+ typedef Dune::GenericReferenceElements<Scalar, dim>
+ ReferenceElementContainer;
+ typedef Dune::GenericReferenceElements<Scalar, dim - 1>
+ ReferenceElementFaceContainer;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(DiffusivePart)) DiffusivePart;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConvectivePart)) ConvectivePart;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(EvalCflFluxFunction)) EvalCflFluxFunction;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConvectivePart))
+ ConvectivePart;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(EvalCflFluxFunction))
+ EvalCflFluxFunction;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters))
+ SpatialParameters;
typedef typename SpatialParameters::MaterialLaw MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
@@ -150,7 +155,8 @@ public:
* Additionally to the \a update vector, the recommended time step size \a dt is calculated
* employing a CFL condition.
*/
- virtual int update(const Scalar t, Scalar& dt, RepresentationType& updateVec, bool impes);
+ virtual int update(const Scalar t, Scalar& dt,
+ RepresentationType& updateVec, bool impes);
//! Sets the initial solution \f$S_0\f$.
void initialize();
@@ -167,8 +173,9 @@ public:
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
- typename Variables::ScalarSolutionType *saturation = writer.template createField<Scalar, 1> (
- problem_.gridView().size(0));
+ typename Variables::ScalarSolutionType *saturation =
+ writer.template createField<Scalar, 1> (
+ problem_.gridView().size(0));
*saturation = problem_.variables().saturation();
@@ -207,14 +214,16 @@ public:
{
if (compressibility_ && velocityType_ == vt)
{
- DUNE_THROW(Dune::NotImplemented,
+ DUNE_THROW(
+ Dune::NotImplemented,
"Total velocity - global pressure - model cannot be used with compressible fluids!");
}
if (saturationType_ != Sw && saturationType_ != Sn)
{
DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
}
- if (pressureType_ != pw && pressureType_ != pn && pressureType_ != pglobal)
+ if (pressureType_ != pw && pressureType_ != pn && pressureType_
+ != pglobal)
{
DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
}
@@ -241,16 +250,21 @@ private:
ConvectivePart* convectivePart_;
EvalCflFluxFunction* evalCflFluxFunction_;
- static const bool compressibility_ = GET_PROP_VALUE(TypeTag, PTAG(EnableCompressibility));
+ static const bool compressibility_ =
+ GET_PROP_VALUE(TypeTag, PTAG(EnableCompressibility));
;
- static const int saturationType_ = GET_PROP_VALUE(TypeTag, PTAG(SaturationFormulation));
- static const int velocityType_ = GET_PROP_VALUE(TypeTag, PTAG(VelocityFormulation));
- static const int pressureType_ = GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation));
+ static const int saturationType_ =
+ GET_PROP_VALUE(TypeTag, PTAG(SaturationFormulation));
+ static const int velocityType_ =
+ GET_PROP_VALUE(TypeTag, PTAG(VelocityFormulation));
+ static const int pressureType_ =
+ GET_PROP_VALUE(TypeTag, PTAG(PressureFormulation));
bool switchNormals_;
};
template<class TypeTag>
-int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationType& updateVec, bool impes = false)
+int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt,
+ RepresentationType& updateVec, bool impes = false)
{
if (!impes)
{
@@ -268,7 +282,8 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
// compute update vector
ElementIterator eItEnd = problem_.gridView().template end<0> ();
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt
+ != eItEnd; ++eIt)
{
//
GlobalPosition globalPos = eIt->geometry().center();
@@ -283,10 +298,12 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
// problem_.variables().storeSrn(residualSatNW, globalIdxI);
//for benchmark only!
- Scalar porosity = problem_.spatialParameters().porosity(globalPos, *eIt);
+ Scalar porosity =
+ problem_.spatialParameters().porosity(globalPos, *eIt);
Scalar viscosityWI = problem_.variables().viscosityWetting(globalIdxI);
- Scalar viscosityNWI = problem_.variables().viscosityNonwetting(globalIdxI);
+ Scalar viscosityNWI = problem_.variables().viscosityNonwetting(
+ globalIdxI);
Scalar densityWI = problem_.variables().densityWetting(globalIdxI);
Scalar densityNWI = problem_.variables().densityNonwetting(globalIdxI);
@@ -296,12 +313,14 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
// run through all intersections with neighbors and boundary
IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
- for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
+ for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt
+ != isItEnd; ++isIt)
{
// local number of facet
int isIndex = isIt->indexInInside();
- Dune::FieldVector<Scalar, dimWorld> unitOuterNormal = isIt->centerUnitOuterNormal();
+ Dune::FieldVector<Scalar, dimWorld> unitOuterNormal =
+ isIt->centerUnitOuterNormal();
if (switchNormals_)
unitOuterNormal *= -1.0;
@@ -318,24 +337,32 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
int globalIdxJ = problem_.variables().index(*neighborPointer);
// neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
+ const GlobalPosition& globalPosNeighbor =
+ neighborPointer->geometry().center();
// distance vector between barycenters
- Dune::FieldVector<Scalar, dimWorld> distVec = globalPosNeighbor - globalPos;
+ Dune::FieldVector<Scalar, dimWorld> distVec = globalPosNeighbor
+ - globalPos;
// compute distance between cell centers
Scalar dist = distVec.two_norm();
//get phase potentials
- Scalar potentialW = problem_.variables().potentialWetting(globalIdxI, isIndex);
- Scalar potentialNW = problem_.variables().potentialNonwetting(globalIdxI, isIndex);
+ Scalar potentialW = problem_.variables().potentialWetting(
+ globalIdxI, isIndex);
+ Scalar potentialNW = problem_.variables().potentialNonwetting(
+ globalIdxI, isIndex);
- Scalar densityWJ = problem_.variables().densityWetting(globalIdxJ);
- Scalar densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
+ Scalar densityWJ = problem_.variables().densityWetting(
+ globalIdxJ);
+ Scalar densityNWJ = problem_.variables().densityNonwetting(
+ globalIdxJ);
//get velocity*normalvector*facearea/(volume*porosity)
- factor = (problem_.variables().velocity()[globalIdxI][isIndex] * unitOuterNormal) * (faceArea);
- factorSecondPhase = (problem_.variables().velocitySecondPhase()[globalIdxI][isIndex] * unitOuterNormal)
- * (faceArea);
+ factor = (problem_.variables().velocity()[globalIdxI][isIndex]
+ * unitOuterNormal) * (faceArea);
+ factorSecondPhase
+ = (problem_.variables().velocitySecondPhase()[globalIdxI][isIndex]
+ * unitOuterNormal) * (faceArea);
Scalar lambdaW = 0;
Scalar lambdaNW = 0;
@@ -359,7 +386,8 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
{
lambdaW /= densityWJ;
}//divide by density because lambda is saved as lambda*density
- viscosityW = problem_.variables().viscosityWetting(globalIdxJ);
+ viscosityW = problem_.variables().viscosityWetting(
+ globalIdxJ);
}
if (potentialNW >= 0.)
@@ -373,12 +401,14 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
}
else
{
- lambdaNW = problem_.variables().mobilityNonwetting(globalIdxJ);
+ lambdaNW = problem_.variables().mobilityNonwetting(
+ globalIdxJ);
if (compressibility_)
{
lambdaNW /= densityNWJ;
}//divide by density because lambda is saved as lambda*density
- viscosityNW = problem_.variables().viscosityNonwetting(globalIdxJ);
+ viscosityNW = problem_.variables().viscosityNonwetting(
+ globalIdxJ);
}
switch (velocityType_)
@@ -386,7 +416,8 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
case vt:
{
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factor, *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt);
//determine phase saturations from primary saturation variable
Scalar satWI = 0;
@@ -401,30 +432,38 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
}
case Sn:
{
- satWI = 1 - problem_.variables().saturation()[globalIdxI];
- satWJ = 1 - problem_.variables().saturation()[globalIdxJ];
+ satWI = 1
+ - problem_.variables().saturation()[globalIdxI];
+ satWJ = 1
+ - problem_.variables().saturation()[globalIdxJ];
break;
}
}
- Scalar pcI = problem_.variables().capillaryPressure(globalIdxI);
- Scalar pcJ = problem_.variables().capillaryPressure(globalIdxJ);
+ Scalar pcI = problem_.variables().capillaryPressure(
+ globalIdxI);
+ Scalar pcJ = problem_.variables().capillaryPressure(
+ globalIdxJ);
// calculate the saturation gradient
- Dune::FieldVector<Scalar, dimWorld> pcGradient = unitOuterNormal;
+ Dune::FieldVector<Scalar, dimWorld> pcGradient =
+ unitOuterNormal;
pcGradient *= (pcI - pcJ) / dist;
// get the diffusive part -> give 1-sat because sat = S_n and lambda = lambda(S_w) and pc = pc(S_w)
- Scalar diffPart = diffusivePart()(*eIt, isIndex, satWI, satWJ, pcGradient) * unitOuterNormal
- * faceArea;
+ Scalar diffPart = diffusivePart()(*eIt, isIndex, satWI,
+ satWJ, pcGradient) * unitOuterNormal * faceArea;
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, 10 * diffPart, *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, 10 * diffPart, *isIt);
- Scalar convPart = convectivePart()(*eIt, isIndex, satWI, satWJ) * unitOuterNormal * faceArea;
+ Scalar convPart = convectivePart()(*eIt, isIndex, satWI,
+ satWJ) * unitOuterNormal * faceArea;
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, 10 * convPart, *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, 10 * convPart, *isIt);
switch (saturationType_)
{
@@ -454,9 +493,10 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
}
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factor, *isIt,
- wPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factorSecondPhase,
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt, wPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factorSecondPhase,
*isIt, nPhaseIdx);
break;
@@ -470,9 +510,10 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
}
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factor, *isIt,
- nPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factorSecondPhase,
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt, nPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factorSecondPhase,
*isIt, wPhaseIdx);
break;
@@ -487,24 +528,32 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
GlobalPosition globalPosFace = isIt->geometry().center();
// distance vector between barycenters
- Dune::FieldVector<Scalar, dimWorld> distVec = globalPosFace - globalPos;
+ Dune::FieldVector<Scalar, dimWorld> distVec = globalPosFace
+ - globalPos;
// compute distance between cell centers
Scalar dist = distVec.two_norm();
//get boundary type
- BoundaryConditions::Flags bcTypeSat = problem_.bctypeSat(globalPosFace, *isIt);
+ BoundaryConditions::Flags bcTypeSat = problem_.bctypeSat(
+ globalPosFace, *isIt);
if (bcTypeSat == BoundaryConditions::dirichlet)
{
- Scalar satBound = problem_.dirichletSat(globalPosFace, *isIt);
+ Scalar satBound = problem_.dirichletSat(globalPosFace,
+ *isIt);
//get velocity*normalvector*facearea/(volume*porosity)
- factor = (problem_.variables().velocity()[globalIdxI][isIndex] * unitOuterNormal) * (faceArea);
- factorSecondPhase = (problem_.variables().velocitySecondPhase()[globalIdxI][isIndex]
- * unitOuterNormal) * (faceArea);
-
- Scalar pressBound = problem_.variables().pressure()[globalIdxI];
- Scalar temperature = problem_.temperature(globalPosFace, *eIt);
+ factor
+ = (problem_.variables().velocity()[globalIdxI][isIndex]
+ * unitOuterNormal) * (faceArea);
+ factorSecondPhase
+ = (problem_.variables().velocitySecondPhase()[globalIdxI][isIndex]
+ * unitOuterNormal) * (faceArea);
+
+ Scalar pressBound =
+ problem_.variables().pressure()[globalIdxI];
+ Scalar temperature = problem_.temperature(globalPosFace,
+ *eIt);
//determine phase saturations from primary saturation variable
Scalar satWI = 0;
@@ -519,14 +568,16 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
}
case Sn:
{
- satWI = 1 - problem_.variables().saturation()[globalIdxI];
+ satWI = 1
+ - problem_.variables().saturation()[globalIdxI];
satWBound = 1 - satBound;
break;
}
}
- Scalar pcBound = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(globalPos, *eIt),
- satBound);
+ Scalar pcBound = MaterialLaw::pC(
+ problem_.spatialParameters().materialLawParams(
+ globalPos, *eIt), satBound);
//determine phase pressures from primary pressure variable
Scalar pressW = 0;
@@ -551,8 +602,11 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
fluidState.update(satBound, pressW, pressNW, temperature);
//get phase potentials
- Scalar potentialW = problem_.variables().potentialWetting(globalIdxI, isIndex);
- Scalar potentialNW = problem_.variables().potentialNonwetting(globalIdxI, isIndex);
+ Scalar potentialW = problem_.variables().potentialWetting(
+ globalIdxI, isIndex);
+ Scalar potentialNW =
+ problem_.variables().potentialNonwetting(
+ globalIdxI, isIndex);
Scalar lambdaW = 0;
Scalar lambdaNW = 0;
@@ -570,14 +624,21 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
{
if (compressibility_)
{
- lambdaW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(globalPos, *eIt),
- satWBound)
- / FluidSystem::phaseViscosity(wPhaseIdx, temperature, pressW, fluidState);
+ lambdaW
+ = MaterialLaw::krw(
+ problem_.spatialParameters().materialLawParams(
+ globalPos, *eIt), satWBound)
+ / FluidSystem::phaseViscosity(
+ wPhaseIdx, temperature,
+ pressW, fluidState);
}
else
{
- lambdaW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(globalPos, *eIt),
- satWBound) / viscosityWI;
+ lambdaW
+ = MaterialLaw::krw(
+ problem_.spatialParameters().materialLawParams(
+ globalPos, *eIt), satWBound)
+ / viscosityWI;
}
}
@@ -593,15 +654,21 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
{
if (compressibility_)
{
- lambdaNW = MaterialLaw::krn(
- problem_.spatialParameters().materialLawParams(globalPos, *eIt), satWBound)
- / FluidSystem::phaseViscosity(nPhaseIdx, temperature, pressNW, fluidState);
+ lambdaNW
+ = MaterialLaw::krn(
+ problem_.spatialParameters().materialLawParams(
+ globalPos, *eIt), satWBound)
+ / FluidSystem::phaseViscosity(
+ nPhaseIdx, temperature,
+ pressNW, fluidState);
}
else
{
- lambdaNW = MaterialLaw::krn(
- problem_.spatialParameters().materialLawParams(globalPos, *eIt), satWBound)
- / viscosityNWI;
+ lambdaNW
+ = MaterialLaw::krn(
+ problem_.spatialParameters().materialLawParams(
+ globalPos, *eIt), satWBound)
+ / viscosityNWI;
}
}
// std::cout<<lambdaW<<" "<<lambdaNW<<std::endl;
@@ -610,28 +677,32 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
{
case vt:
{
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factor, *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt);
- Scalar pcI = problem_.variables().capillaryPressure(globalIdxI);
+ Scalar pcI = problem_.variables().capillaryPressure(
+ globalIdxI);
// calculate the saturation gradient
- Dune::FieldVector<Scalar, dimWorld> pcGradient = unitOuterNormal;
+ Dune::FieldVector<Scalar, dimWorld> pcGradient =
+ unitOuterNormal;
pcGradient *= (pcI - pcBound) / dist;
// get the diffusive part -> give 1-sat because sat = S_n and lambda = lambda(S_w) and pc = pc(S_w)
- Scalar diffPart = diffusivePart()(*eIt, isIndex, satWI, satWBound, pcGradient)
- * unitOuterNormal * faceArea;
+ Scalar diffPart = diffusivePart()(*eIt, isIndex, satWI,
+ satWBound, pcGradient) * unitOuterNormal
+ * faceArea;
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, 10 * diffPart,
- *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, 10 * diffPart, *isIt);
- Scalar convPart = convectivePart()(*eIt, isIndex, satWI, satWBound) * unitOuterNormal
- * faceArea;
+ Scalar convPart = convectivePart()(*eIt, isIndex,
+ satWI, satWBound) * unitOuterNormal * faceArea;
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, 10 * convPart,
- *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, 10 * convPart, *isIt);
switch (saturationType_)
{
@@ -663,9 +734,11 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
}
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factor, *isIt,
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt,
wPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factorSecondPhase,
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factorSecondPhase,
*isIt, nPhaseIdx);
break;
@@ -681,9 +754,11 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
}
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factor, *isIt,
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt,
nPhaseIdx);
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factorSecondPhase,
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factorSecondPhase,
*isIt, wPhaseIdx);
break;
@@ -708,72 +783,117 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
lambdaNW /= densityNWI;
}
- //get velocity*normalvector*facearea/(volume*porosity)
- factor = (problem_.variables().velocity()[globalIdxI][isIndex] * unitOuterNormal) * faceArea;
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ factor = problem_.neumann(globalPosFace, *isIt)[wPhaseIdx];
+ factor /= densityWI * faceArea;
+ break;
+ }
+ case Sn:
+ {
+ factor = problem_.neumann(globalPosFace, *isIt)[nPhaseIdx];
+ factor /= densityNWI * faceArea;
+ break;
+ }
+ }
+
switch (velocityType_)
{
case vt:
{
- switch (saturationType_)
- {
- case Sw:
- {
- //vt*fw
- factor *= lambdaW / (lambdaW + lambdaNW);
- break;
- }
- case Sn:
- {
- //vt*fn
- factor *= lambdaNW / (lambdaW + lambdaNW);
- break;
- }
- }
+ //add cflFlux for time-stepping
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt);
+ }
+ case vw:
+ {
+ //add cflFlux for time-stepping
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt,
+ wPhaseIdx);
break;
}
+ case vn:
+ {
+ //add cflFlux for time-stepping
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt,
+ nPhaseIdx);
+ break;
+ }
+ }
+
+ }//end neumann boundary
+ if (bcTypeSat == BoundaryConditions::outflow)
+ {
+ //get mobilities
+ Scalar lambdaW, lambdaNW;
+
+ lambdaW = lambdaWI;
+ if (compressibility_)
+ {
+ lambdaW /= densityWI;
+ }
+
+ lambdaNW = lambdaNWI;
+ if (compressibility_)
+ {
+ lambdaNW /= densityNWI;
}
- Scalar boundaryFactor = problem_.neumannSat(globalPosFace, *isIt, factor);
- if (factor != boundaryFactor)
+ //get velocity*normalvector*facearea/(volume*porosity)
+ factor
+ = (problem_.variables().velocity()[globalIdxI][isIndex]
+ * unitOuterNormal) * faceArea;
+
+ if (velocityType_ == vt)
{
switch (saturationType_)
{
case Sw:
{
- factor = boundaryFactor / densityWI * faceArea;
+ //vt*fw
+ factor *= lambdaW / (lambdaW + lambdaNW);
break;
}
case Sn:
{
- factor = boundaryFactor / densityNWI * faceArea;
+ //vt*fn
+ factor *= lambdaNW / (lambdaW + lambdaNW);
break;
}
}
}
+
switch (velocityType_)
{
case vt:
{
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factor, *isIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt);
}
case vw:
{
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factor, *isIt,
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt,
wPhaseIdx);
break;
}
case vn:
{
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, factor, *isIt,
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW,
+ viscosityWI, viscosityNWI, factor, *isIt,
nPhaseIdx);
break;
}
}
- }//end neumann boundary
+ }//end outflow boundary
}//end boundary
// add to update vector
updateVec[globalIdxI] -= factor / (volume * porosity);//-:v>0, if flow leaves the cell
@@ -834,30 +954,33 @@ int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationTy
case vw:
{
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, source * volume, *eIt,
- wPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI,
+ viscosityNWI, source * volume, *eIt, wPhaseIdx);
break;
}
case vn:
{
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, source * volume, *eIt,
- nPhaseIdx);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI,
+ viscosityNWI, source * volume, *eIt, nPhaseIdx);
break;
}
case vt:
{
//add cflFlux for time-stepping
- evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI, viscosityNWI, source * volume, *eIt);
+ evalCflFluxFunction().addFlux(lambdaW, lambdaNW, viscosityWI,
+ viscosityNWI, source * volume, *eIt);
break;
}
}
}
//calculate time step
- dt = std::min(dt, evalCflFluxFunction().getCflFluxFunction(globalPos, *eIt) * (porosity * volume));
+ dt = std::min(dt, evalCflFluxFunction().getCflFluxFunction(globalPos,
+ *eIt) * (porosity * volume));
- problem_.variables().volumecorrection(globalIdxI) = updateVec[globalIdxI];
+ problem_.variables().volumecorrection(globalIdxI)
+ = updateVec[globalIdxI];
} // end grid traversal
return 0;
@@ -868,26 +991,30 @@ void FVSaturation2P<TypeTag>::initialize()
{
// iterate through leaf grid an evaluate c0 at cell center
ElementIterator eItEnd = problem_.gridView().template end<0> ();
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt
+ != eItEnd; ++eIt)
{
// get global coordinate of cell center
GlobalPosition globalPos = eIt->geometry().center();
// initialize cell concentration
- problem_.variables().saturation()[problem_.variables().index(*eIt)] = problem_.initSat(globalPos, *eIt);
+ problem_.variables().saturation()[problem_.variables().index(*eIt)]
+ = problem_.initSat(globalPos, *eIt);
}
return;
}
template<class TypeTag>
-void FVSaturation2P<TypeTag>::updateMaterialLaws(RepresentationType& saturation = *(new RepresentationType(0)),
+void FVSaturation2P<TypeTag>::updateMaterialLaws(
+ RepresentationType& saturation = *(new RepresentationType(0)),
bool iterate = false)
{
FluidState fluidState;
// iterate through leaf grid an evaluate c0 at cell center
ElementIterator eItEnd = problem_.gridView().template end<0> ();
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt
+ != eItEnd; ++eIt)
{
// get global coordinate of cell center
GlobalPosition globalPos = eIt->geometry().center();
@@ -918,33 +1045,48 @@ void FVSaturation2P<TypeTag>::updateMaterialLaws(RepresentationType& saturation
Scalar temperature = problem_.temperature(globalPos, *eIt);
Scalar referencePressure = problem_.referencePressure(globalPos, *eIt);
- fluidState.update(satW, referencePressure, referencePressure, temperature);
-
- problem_.variables().densityWetting(globalIdx) = FluidSystem::phaseDensity(wPhaseIdx, temperature,
- referencePressure, fluidState);
- problem_.variables().densityNonwetting(globalIdx) = FluidSystem::phaseDensity(nPhaseIdx, temperature,
- referencePressure, fluidState);
- problem_.variables().viscosityWetting(globalIdx) = FluidSystem::phaseViscosity(wPhaseIdx, temperature,
- referencePressure, fluidState);
- problem_.variables().viscosityNonwetting(globalIdx) = FluidSystem::phaseViscosity(nPhaseIdx, temperature,
- referencePressure, fluidState);
+ fluidState.update(satW, referencePressure, referencePressure,
+ temperature);
+
+ problem_.variables().densityWetting(globalIdx)
+ = FluidSystem::phaseDensity(wPhaseIdx, temperature,
+ referencePressure, fluidState);
+ problem_.variables().densityNonwetting(globalIdx)
+ = FluidSystem::phaseDensity(nPhaseIdx, temperature,
+ referencePressure, fluidState);
+ problem_.variables().viscosityWetting(globalIdx)
+ = FluidSystem::phaseViscosity(wPhaseIdx, temperature,
+ referencePressure, fluidState);
+ problem_.variables().viscosityNonwetting(globalIdx)
+ = FluidSystem::phaseViscosity(nPhaseIdx, temperature,
+ referencePressure, fluidState);
// initialize mobilities
- problem_.variables().mobilityWetting(globalIdx) = MaterialLaw::krw(
- problem_.spatialParameters().materialLawParams(globalPos, *eIt), satW)
- / problem_.variables().viscosityWetting(globalIdx);
- problem_.variables().mobilityNonwetting(globalIdx) = MaterialLaw::krn(
- problem_.spatialParameters().materialLawParams(globalPos, *eIt), satW)
- / problem_.variables().viscosityNonwetting(globalIdx);
- problem_.variables().capillaryPressure(globalIdx) = MaterialLaw::pC(
- problem_.spatialParameters().materialLawParams(globalPos, *eIt), satW);
+ problem_.variables().mobilityWetting(globalIdx)
+ = MaterialLaw::krw(
+ problem_.spatialParameters().materialLawParams(
+ globalPos, *eIt), satW)
+ / problem_.variables().viscosityWetting(globalIdx);
+ problem_.variables().mobilityNonwetting(globalIdx)
+ = MaterialLaw::krn(
+ problem_.spatialParameters().materialLawParams(
+ globalPos, *eIt), satW)
+ / problem_.variables().viscosityNonwetting(globalIdx);
+ problem_.variables().capillaryPressure(globalIdx)
+ = MaterialLaw::pC(
+ problem_.spatialParameters().materialLawParams(
+ globalPos, *eIt), satW);
problem_.variables().fracFlowFuncWetting(globalIdx)
- = problem_.variables().mobilityWetting(globalIdx) / (problem_.variables().mobilityWetting(globalIdx)
- + problem_.variables().mobilityNonwetting(globalIdx));
+ = problem_.variables().mobilityWetting(globalIdx)
+ / (problem_.variables().mobilityWetting(globalIdx)
+ + problem_.variables().mobilityNonwetting(
+ globalIdx));
problem_.variables().fracFlowFuncNonwetting(globalIdx)
- = problem_.variables().mobilityNonwetting(globalIdx) / (problem_.variables().mobilityWetting(globalIdx)
- + problem_.variables().mobilityNonwetting(globalIdx));
+ = problem_.variables().mobilityNonwetting(globalIdx)
+ / (problem_.variables().mobilityWetting(globalIdx)
+ + problem_.variables().mobilityNonwetting(
+ globalIdx));
problem_.spatialParameters().update(satW, *eIt);
}
diff --git a/test/decoupled/1p/test_diffusion.cc b/test/decoupled/1p/test_diffusion.cc
index d67efc30f3b4dd21610cce3408bbea3374384c57..d557532db2a07e19e279175e616a1f9cc5a96172 100644
--- a/test/decoupled/1p/test_diffusion.cc
+++ b/test/decoupled/1p/test_diffusion.cc
@@ -87,6 +87,7 @@ int main(int argc, char** argv)
fvProblem.init();
fvProblem.model().calculateVelocity();
double fvTime = timer.elapsed();
+ fvProblem.writeOutput();
Dumux::ResultEvaluation fvResult;
fvResult.evaluate(grid.leafView(), fvProblem, fvProblem.variables().pressure(), fvProblem.variables().velocity(), consecutiveNumbering);
@@ -96,6 +97,7 @@ int main(int argc, char** argv)
mpfaProblem.init();
mpfaProblem.model().calculateVelocity();
double mpfaTime = timer.elapsed();
+ mpfaProblem.writeOutput();
Dumux::ResultEvaluation mpfaResult;
mpfaResult.evaluate(grid.leafView(), mpfaProblem, mpfaProblem.variables().pressure(), mpfaProblem.variables().velocity(), consecutiveNumbering);
@@ -105,6 +107,7 @@ int main(int argc, char** argv)
mimeticProblem.init();
mimeticProblem.model().calculateVelocity();
double mimeticTime = timer.elapsed();
+ mimeticProblem.writeOutput();
Dumux::ResultEvaluation mimeticResult;
mimeticResult.evaluate(grid.leafView(), mimeticProblem, mimeticProblem.variables().pressure(), mimeticProblem.variables().velocity(), consecutiveNumbering);
diff --git a/test/decoupled/1p/test_diffusion_problem.hh b/test/decoupled/1p/test_diffusion_problem.hh
index ba20b597a4e55ac8ef5328c38f7a3d72467630f7..04a0a4835e6546f9c4c500893a3ff914a9f5d971 100644
--- a/test/decoupled/1p/test_diffusion_problem.hh
+++ b/test/decoupled/1p/test_diffusion_problem.hh
@@ -269,7 +269,7 @@ public:
return 1.0;
}
- std::vector<Scalar> neumannPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+ std::vector<Scalar> neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
{
std::vector<Scalar> neumannFlux(2, 0.0);
diff --git a/test/decoupled/2p/test_impes_problem.hh b/test/decoupled/2p/test_impes_problem.hh
index 9ee86f535dabe826267e9486d1984ef942ad8f9b..002d78d55dedbd4b140d57127c4e2f5d9db36793 100644
--- a/test/decoupled/2p/test_impes_problem.hh
+++ b/test/decoupled/2p/test_impes_problem.hh
@@ -222,9 +222,10 @@ typename BoundaryConditions::Flags bctypePress(const GlobalPosition& globalPos,
BoundaryConditions::Flags bctypeSat(const GlobalPosition& globalPos, const Intersection& intersection) const
{
- if (globalPos[0] > (upperRight_[0] - eps_) || globalPos[0] < eps_)
-// if (globalPos[0] < eps_)
+ if (globalPos[0] < eps_)
return Dumux::BoundaryConditions::dirichlet;
+ else if (globalPos[0] > upperRight_[0] - eps_)
+ return Dumux::BoundaryConditions::outflow;
else
return Dumux::BoundaryConditions::neumann;
}
@@ -260,7 +261,7 @@ Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& interse
return 0.2;
}
-std::vector<Scalar> neumannPress(const GlobalPosition& globalPos, const Intersection& intersection) const
+std::vector<Scalar> neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
{
std::vector<Scalar> neumannFlux(2, 0.0);
if (globalPos[0] > upperRight_[0] - eps_)
@@ -270,13 +271,6 @@ std::vector<Scalar> neumannPress(const GlobalPosition& globalPos, const Intersec
return neumannFlux;
}
-Scalar neumannSat(const GlobalPosition& globalPos, const Intersection& intersection, Scalar factor) const
-{
- if (globalPos[0] > upperRight_[0] - eps_)
- return factor;
- return 0;
-}
-
Scalar initSat(const GlobalPosition& globalPos, const Element& element) const
{
return 0.2;
diff --git a/test/decoupled/2p/test_transport_problem.hh b/test/decoupled/2p/test_transport_problem.hh
index b153bdb6a984bce1cbd652ec7f0ca36aae4e56ca..40532d451e428211e9c1334a7102097183eca45a 100644
--- a/test/decoupled/2p/test_transport_problem.hh
+++ b/test/decoupled/2p/test_transport_problem.hh
@@ -189,17 +189,18 @@ public:
}
std::vector<Scalar> source(const GlobalPosition& globalPos, const Element& element)
- {
- return std::vector<Scalar>(2, 0.0);
- }
+ {
+ return std::vector<Scalar>(2, 0.0);
+ }
BoundaryConditions::Flags bctypeSat(const GlobalPosition& globalPos, const Intersection& intersection) const
{
- if (globalPos[0] > (upperRight_[0] - eps_) || globalPos[0] < eps_)
- // if (globalPos[0] < eps_)
- return Dumux::BoundaryConditions::dirichlet;
- else
- return Dumux::BoundaryConditions::neumann;
+ if (globalPos[0] < eps_)
+ return Dumux::BoundaryConditions::dirichlet;
+ else if (globalPos[0] > upperRight_[0] - eps_)
+ return Dumux::BoundaryConditions::outflow;
+ else
+ return Dumux::BoundaryConditions::neumann;
}
Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& intersection) const
@@ -210,11 +211,9 @@ public:
return 0.0;
}
- Scalar neumannSat(const GlobalPosition& globalPos, const Intersection& intersection, Scalar factor) const
+ std::vector<Scalar> neumann(const GlobalPosition& globalPos, const Intersection& intersection) const
{
- if (globalPos[0] > upperRight_[0] - eps_)
- return factor;
- return 0;
+ return std::vector<Scalar>(2, 0.0);
}
Scalar initSat(const GlobalPosition& globalPos, const Element& element) const