From 177c9fce2f023c08548280236f15aac0f5a77dc1 Mon Sep 17 00:00:00 2001
From: Bernd Flemisch <bernd@iws.uni-stuttgart.de>
Date: Thu, 12 Aug 2010 14:27:08 +0000
Subject: [PATCH] finished transport test
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4062 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/common/timemanager.hh | 3 +-
.../2p/transport/fv/fvsaturation2p.hh | 2085 +++++++++--------
.../2p/transport/transportproperties.hh | 35 +-
dumux/decoupled/2p/variableclass2p.hh | 2 +-
dumux/decoupled/common/onemodelproblem.hh | 20 +-
dumux/decoupled/common/variableclass.hh | 57 +-
test/decoupled/2p/grids/test_transport.dgf | 11 +
test/decoupled/2p/test_transport.cc | 23 +-
test/decoupled/2p/test_transport_problem.hh | 6 +-
9 files changed, 1168 insertions(+), 1074 deletions(-)
create mode 100644 test/decoupled/2p/grids/test_transport.dgf
diff --git a/dumux/common/timemanager.hh b/dumux/common/timemanager.hh
index 16d31364dc..ac118299e1 100644
--- a/dumux/common/timemanager.hh
+++ b/dumux/common/timemanager.hh
@@ -322,8 +322,7 @@ public:
if (problem_->shouldWriteOutput())
problem_->writeOutput();
- // advance the simulated time by the current time step
- // size
+ // advance the simulated time by the current time step size
Scalar dt = timeStepSize();
time_ += timeStepSize_;
++timeStepIdx_;
diff --git a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
index fca60f982d..1fd9020beb 100644
--- a/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
+++ b/dumux/decoupled/2p/transport/fv/fvsaturation2p.hh
@@ -57,1073 +57,1086 @@ namespace Dumux
template<class TypeTag>
class FVSaturation2P
{
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
- typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) ReferenceElements;
- typedef typename ReferenceElements::Container ReferenceElementContainer;
- typedef typename ReferenceElements::ContainerFaces 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(SpatialParameters)) SpatialParameters;
- typedef typename SpatialParameters::MaterialLaw MaterialLaw;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
-
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
-
- enum
- {
- dim = GridView::dimension, dimWorld = GridView::dimensionworld
- };
- enum
- {
- pw = Indices::pressureW,
- pn = Indices::pressureNW,
- pglobal = Indices::pressureGlobal,
- vw = Indices::velocityW,
- vn = Indices::velocityNW,
- vt = Indices::velocityTotal,
- Sw = Indices::saturationW,
- Sn = Indices::saturationNW,
- };
- enum
- {
- wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
- };
-
- typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
- typedef typename SolutionTypes::ScalarSolution RepresentationType;
- typedef typename GridView::Traits::template Codim<0>::Entity Element;
- typedef typename GridView::Grid Grid;
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
- typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
- typedef typename GridView::IntersectionIterator IntersectionIterator;
-
- typedef Dune::FieldVector<Scalar, dim> LocalPosition;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
- DiffusivePart& diffusivePart()
- {
- return *diffusivePart_;
- }
-
- const DiffusivePart& diffusivePart() const
- {
- return *diffusivePart_;
- }
-
- ConvectivePart& convectivPart()
- {
- return *convectivePart_;
- }
-
- const ConvectivePart& convectivePart() const
- {
- return *convectivePart_;
- }
-
- //function to calculate the time step if a non-wetting phase velocity is used
- Scalar evaluateTimeStepPhaseFlux(Scalar timestepFactorIn, Scalar timestepFactorOutNW, Scalar& residualSatW,
- Scalar& residualSatNW, int globalIdxI);
-
- //function to calculate the time step if a total velocity is used
- Scalar evaluateTimeStepTotalFlux(Scalar timestepFactorIn, Scalar timestepFactorOut, Scalar diffFactorIn,
- Scalar diffFactorOut, Scalar& residualSatW, Scalar& residualSatNW);
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) ReferenceElements;
+ typedef typename ReferenceElements::Container ReferenceElementContainer;
+ typedef typename ReferenceElements::ContainerFaces 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(SpatialParameters)) SpatialParameters;
+ typedef typename SpatialParameters::MaterialLaw MaterialLaw;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+ enum
+ {
+ pw = Indices::pressureW,
+ pn = Indices::pressureNW,
+ pglobal = Indices::pressureGlobal,
+ vw = Indices::velocityW,
+ vn = Indices::velocityNW,
+ vt = Indices::velocityTotal,
+ Sw = Indices::saturationW,
+ Sn = Indices::saturationNW,
+ };
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+ };
+
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::ScalarSolution RepresentationType;
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+ typedef typename GridView::Grid Grid;
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+ DiffusivePart& diffusivePart()
+ {
+ return *diffusivePart_;
+ }
+
+ const DiffusivePart& diffusivePart() const
+ {
+ return *diffusivePart_;
+ }
+
+ ConvectivePart& convectivePart()
+ {
+ return *convectivePart_;
+ }
+
+ const ConvectivePart& convectivePart() const
+ {
+ return *convectivePart_;
+ }
+
+ //function to calculate the time step if a non-wetting phase velocity is used
+ Scalar evaluateTimeStepPhaseFlux(Scalar timestepFactorIn, Scalar timestepFactorOutNW, Scalar& residualSatW,
+ Scalar& residualSatNW, int globalIdxI);
+
+ //function to calculate the time step if a total velocity is used
+ Scalar evaluateTimeStepTotalFlux(Scalar timestepFactorIn, Scalar timestepFactorOut, Scalar diffFactorIn,
+ Scalar diffFactorOut, Scalar& residualSatW, Scalar& residualSatNW);
public:
- //! Calculate the update vector.
- /*!
- * \param[in] t time
- * \param[in] dt time step size
- * \param[in] updateVec vector for the update values
- * \param[in] CLFFac security factor for the time step criterion (0 < CLFFac <= 1)
- * \param[in] impes variable is true if an impes algorithm is used and false if the transport part is solved independently
- *
- * This method calculates the update vector \f$ u \f$ of the discretized equation
- * \f[
- * S_{n_{new}} = S_{n_{old}} - u,
- * \f]
- * where \f$ u = \sum_{element faces} \boldsymbol{v}_n * \boldsymbol{n} * A_{element face}\f$, \f$\boldsymbol{n}\f$ is the face normal and \f$A_{element face}\f$ is the face area.
- *
- * 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);
-
- //! Sets the initial solution \f$S_0\f$.
- void initialize();
-
- //! Update the values of the material laws and constitutive relations.
- /*!
- * Constitutive relations like capillary pressure-saturation relationships, mobility-saturation relationships... are updated and stored in the variable class
- * of type Dumux::VariableClass2P. The update has to be done when new saturation are available.
- */
- void updateMaterialLaws(RepresentationType& saturation, bool iterate);
-
- //! Constructs a FVSaturationNonWetting2P object
- /**
- * \param gridView gridView object of type GridView
- * \param problem a problem class object
- * \param velocityType a string giving the type of velocity used (could be: vn, vt)
- * \param diffPart a object of class Dune::DiffusivePart or derived from Dune::DiffusivePart (only used with vt)
- */
-
- FVSaturation2P(Problem& problem) :
- problem_(problem), switchNormals_(false)
- {
- if (compressibility_ && velocityType_ == vt)
- {
- 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)
- {
- DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
- }
- if (velocityType_ != vw && velocityType_ != vn && velocityType_ != vt)
- {
- DUNE_THROW(Dune::NotImplemented, "Velocity type not supported!");
- }
-
-
- diffusivePart_ = new DiffusivePart(problem);
- convectivePart_ = new ConvectivePart(problem);
- }
-
- ~FVSaturation2P()
- {
- delete diffusivePart_;
- delete convectivePart_;
- }
+ //! Calculate the update vector.
+ /*!
+ * \param[in] t time
+ * \param[in] dt time step size
+ * \param[in] updateVec vector for the update values
+ * \param[in] CLFFac security factor for the time step criterion (0 < CLFFac <= 1)
+ * \param[in] impes variable is true if an impes algorithm is used and false if the transport part is solved independently
+ *
+ * This method calculates the update vector \f$ u \f$ of the discretized equation
+ * \f[
+ * S_{n_{new}} = S_{n_{old}} - u,
+ * \f]
+ * where \f$ u = \sum_{element faces} \boldsymbol{v}_n * \boldsymbol{n} * A_{element face}\f$, \f$\boldsymbol{n}\f$ is the face normal and \f$A_{element face}\f$ is the face area.
+ *
+ * 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);
+
+ //! Sets the initial solution \f$S_0\f$.
+ void initialize();
+
+ //! Update the values of the material laws and constitutive relations.
+ /*!
+ * Constitutive relations like capillary pressure-saturation relationships, mobility-saturation relationships... are updated and stored in the variable class
+ * of type Dumux::VariableClass2P. The update has to be done when new saturation are available.
+ */
+ void updateMaterialLaws(RepresentationType& saturation, bool iterate);
+
+ template<class MultiWriter>
+ void addOutputVtkFields(MultiWriter &writer)
+ {
+ problem_.variables().addOutputVtkFields(writer);
+ return;
+ }
+
+ // serialization methods
+ template<class Restarter>
+ void serialize(Restarter &res)
+ {
+ problem_.variables().serialize<Restarter> (res);
+ }
+ template<class Restarter>
+ void deserialize(Restarter &res)
+ {
+ problem_.variables().deserialize<Restarter> (res);
+ }
+
+ //! Constructs a FVSaturationNonWetting2P object
+ /**
+ * \param gridView gridView object of type GridView
+ * \param problem a problem class object
+ * \param velocityType a string giving the type of velocity used (could be: vn, vt)
+ * \param diffPart a object of class Dune::DiffusivePart or derived from Dune::DiffusivePart (only used with vt)
+ */
+
+ FVSaturation2P(Problem& problem) :
+ problem_(problem), switchNormals_(false)
+ {
+ if (compressibility_ && velocityType_ == vt)
+ {
+ 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)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
+ }
+ if (velocityType_ != vw && velocityType_ != vn && velocityType_ != vt)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Velocity type not supported!");
+ }
+
+
+ diffusivePart_ = new DiffusivePart(problem);
+ convectivePart_ = new ConvectivePart(problem);
+ }
+
+ ~FVSaturation2P()
+ {
+ delete diffusivePart_;
+ delete convectivePart_;
+ }
private:
- Problem& problem_;
- DiffusivePart* diffusivePart_;
- ConvectivePart* convectivePart_;
-
- 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));
- bool switchNormals_;
+ Problem& problem_;
+ DiffusivePart* diffusivePart_;
+ ConvectivePart* convectivePart_;
+
+ 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));
+ bool switchNormals_;
};
template<class TypeTag>
int FVSaturation2P<TypeTag>::update(const Scalar t, Scalar& dt, RepresentationType& updateVec, bool impes = false)
{
- if (!impes)
- {
- updateMaterialLaws();
- }
+ if (!impes)
+ {
+ updateMaterialLaws();
+ }
- // initialize dt very large
- dt = 1E100;
+ // initialize dt very large
+ dt = 1E100;
- // set update vector to zero
- updateVec = 0;
+ // set update vector to zero
+ updateVec = 0;
- // some phase properties
- Dune::FieldVector<Scalar, dimWorld> gravity = problem_.gravity();
+ // some phase properties
+ Dune::FieldVector<Scalar, dimWorld> gravity = problem_.gravity();
- // compute update vector
- ElementIterator eItEnd = problem_.gridView().template end<0> ();
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
- {
- // cell geometry type
- Dune::GeometryType gt = eIt->geometry().type();
+ // compute update vector
+ ElementIterator eItEnd = problem_.gridView().template end<0> ();
+ for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ {
+ // cell geometry type
+ Dune::GeometryType gt = eIt->geometry().type();
- // cell center in reference element
- const LocalPosition &localPos = ReferenceElementContainer::general(gt).position(0, 0);
+ // cell center in reference element
+ const LocalPosition &localPos = ReferenceElementContainer::general(gt).position(0, 0);
- //
- GlobalPosition globalPos = eIt->geometry().global(localPos);
+ //
+ GlobalPosition globalPos = eIt->geometry().global(localPos);
- // cell volume, assume linear map here
- Scalar volume = eIt->geometry().volume();
+ // cell volume, assume linear map here
+ Scalar volume = eIt->geometry().volume();
- // cell index
- int globalIdxI = problem_.variables().index(*eIt);
+ // cell index
+ int globalIdxI = problem_.variables().index(*eIt);
- Scalar residualSatW = problem_.spatialParameters().materialLawParams(globalPos, *eIt).Swr();
- Scalar residualSatNW = problem_.spatialParameters().materialLawParams(globalPos, *eIt).Snr();
-
- //for benchmark only!
- // problem_.variables().storeSrn(residualSatNW, globalIdxI);
- //for benchmark only!
-
- Scalar porosity = problem_.spatialParameters().porosity(globalPos, *eIt);
-
- Scalar viscosityWI = problem_.variables().viscosityWetting(globalIdxI);
- Scalar viscosityNWI = problem_.variables().viscosityNonwetting(globalIdxI);
- Scalar viscosityRatio = 1 - fabs(0.5 - viscosityNWI / (viscosityWI + viscosityNWI));//1 - fabs(viscosityWI-viscosityNWI)/(viscosityWI+viscosityNWI);
-
- Scalar densityWI = problem_.variables().densityWetting(globalIdxI);
- Scalar densityNWI = problem_.variables().densityNonwetting(globalIdxI);
-
- Scalar lambdaWI = problem_.variables().mobilityWetting(globalIdxI);
- Scalar lambdaNWI = problem_.variables().mobilityNonwetting(globalIdxI);
-
- Scalar timestepFactorIn = 0;
- Scalar timestepFactorOut = 0;
- Scalar diffFactorIn = 0;
- Scalar diffFactorOut = 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)
- {
- // local number of facet
- int isIndex = isIt->indexInInside();
-
- // get geometry type of face
- Dune::GeometryType faceGT = isIt->geometryInInside().type();
-
- // center in face's reference element
- const Dune::FieldVector<Scalar, dim - 1>& faceLocal =
- ReferenceElementFaceContainer::general(faceGT).position(0, 0);
-
- // center of face inside volume reference element
- const LocalPosition localPosFace(0);
-
- Dune::FieldVector<Scalar, dimWorld> unitOuterNormal = isIt->unitOuterNormal(faceLocal);
- if (switchNormals_)
- unitOuterNormal *= -1.0;
-
- Scalar faceArea = isIt->geometry().volume();
-
- Scalar factor = 0;
- Scalar factorSecondPhase = 0;
-
- // handle interior face
- if (isIt->neighbor())
- {
- // access neighbor
- ElementPointer neighborPointer = isIt->outside();
- int globalIdxJ = problem_.variables().index(*neighborPointer);
-
- // compute factor in neighbor
- Dune::GeometryType neighborGT = neighborPointer->geometry().type();
- const LocalPosition& localPosNeighbor = ReferenceElementContainer::general(neighborGT).position(0, 0);
-
- // cell center in global coordinates
- const GlobalPosition& globalPos = eIt->geometry().global(localPos);
-
- // neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().global(localPosNeighbor);
-
- // distance vector between barycenters
- Dune::FieldVector<Scalar, dimWorld> distVec = globalPosNeighbor - globalPos;
- // compute distance between cell centers
- Scalar dist = distVec.two_norm();
-
- Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
- unitDistVec /= dist;
-
- //get phase potentials
- 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);
-
- //get velocity*normalvector*facearea/(volume*porosity)
- factor = (problem_.variables().velocity()[globalIdxI][isIndex] * unitDistVec) * (faceArea
- * (unitOuterNormal * unitDistVec)) / (volume * porosity);
- factorSecondPhase = (problem_.variables().velocitySecondPhase()[globalIdxI][isIndex] * unitDistVec)
- * (faceArea * (unitOuterNormal * unitDistVec)) / (volume * porosity);
-
- Scalar lambdaW = 0;
- Scalar lambdaNW = 0;
- Scalar viscosityW = 0;
- Scalar viscosityNW = 0;
-
- //upwinding of lambda dependend on the phase potential gradients
- if (potentialW >= 0.)
- {
- lambdaW = lambdaWI;
- if (compressibility_)
- {
- lambdaW /= densityWI;
- }//divide by density because lambda is saved as lambda*density
- viscosityW = viscosityWI;
- }
- else
- {
- lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
- if (compressibility_)
- {
- lambdaW /= densityWJ;
- }//divide by density because lambda is saved as lambda*density
- viscosityW = problem_.variables().viscosityWetting(globalIdxJ);
- }
-
- if (potentialNW >= 0.)
- {
- lambdaNW = lambdaNWI;
- if (compressibility_)
- {
- lambdaNW /= densityNWI;
- }//divide by density because lambda is saved as lambda*density
- viscosityNW = viscosityNWI;
- }
- else
- {
- lambdaNW = problem_.variables().mobilityNonwetting(globalIdxJ);
- if (compressibility_)
- {
- lambdaNW /= densityNWJ;
- }//divide by density because lambda is saved as lambda*density
- viscosityNW = problem_.variables().viscosityNonwetting(globalIdxJ);
- }
- Scalar krSum = lambdaW * viscosityW + lambdaNW * viscosityNW;
-
- switch (velocityType_)
- {
- case vt:
- {
- //for time step criterion
-
- if (factor >= 0)
- {
- timestepFactorOut += factor / (krSum * viscosityRatio);
- }
- if (factor < 0)
- {
- timestepFactorIn -= factor / (krSum * viscosityRatio);
- }
-
- //determine phase saturations from primary saturation variable
- Scalar satWI = 0;
- Scalar satWJ = 0;
- switch (saturationType_)
- {
- case Sw:
- {
- satWI = problem_.variables().saturation()[globalIdxI];
- satWJ = problem_.variables().saturation()[globalIdxJ];
- break;
- }
- case Sn:
- {
- 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);
-
- // calculate the saturation gradient
- Dune::FieldVector<Scalar, dimWorld> pcGradient = unitDistVec;
- 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) * unitDistVec * faceArea
- / (volume * porosity) * (unitOuterNormal * unitDistVec);
- Scalar convPart = convectivePart() (*eIt, isIndex, satWI, satWJ) * unitDistVec * faceArea / (volume * porosity) * (unitOuterNormal * unitDistVec);
-
- //for time step criterion
- if (diffPart >= 0)
- {
- diffFactorIn += diffPart / (krSum * viscosityRatio);
- }
- if (diffPart < 0)
- {
- diffFactorOut -= diffPart / (krSum * viscosityRatio);
- }
- if (convPart >= 0)
- {
- timestepFactorOut += convPart / (krSum * viscosityRatio);
- }
- if (convPart < 0)
- {
- timestepFactorIn -= convPart / (krSum * viscosityRatio);
- }
-
- switch (saturationType_)
- {
- case Sw:
- {
- //vt*fw
- factor *= lambdaW / (lambdaW + lambdaNW);
- break;
- }
- case Sn:
- {
- //vt*fn
- factor *= lambdaNW / (lambdaW + lambdaNW);
- break;
- }
- }
- factor -= diffPart;
- factor += convPart;
- break;
- }
- case vw:
- {
- if (compressibility_)
- {
- factor /= densityWI;
- factorSecondPhase /= densityNWI;
- }
-
- //for time step criterion
- if (factor >= 0)
- {
- timestepFactorOut += factor / (krSum * viscosityRatio);
- }
- if (factor < 0)
- {
- timestepFactorIn -= factor / (krSum * viscosityRatio);
- }
- if (factorSecondPhase < 0)
- {
- timestepFactorIn -= factorSecondPhase / (krSum * viscosityRatio);
- }
-
- if (std::isnan(timestepFactorIn) || std::isinf(timestepFactorIn))
- {
- timestepFactorIn = 1e-100;
- }
- break;
- }
-
- //for time step criterion if the non-wetting phase velocity is used
- case vn:
- {
- if (compressibility_)
- {
- factor /= densityNWI;
- factorSecondPhase /= densityWI;
- }
-
- //for time step criterion
- if (factor >= 0)
- {
- timestepFactorOut += factor / (krSum * viscosityRatio);
- }
- if (factor < 0)
- {
- timestepFactorIn -= factor / (krSum * viscosityRatio);
- }
- if (factorSecondPhase < 0)
- {
- timestepFactorIn -= factorSecondPhase / (krSum * viscosityRatio);
- }
-
- if (std::isnan(timestepFactorIn) || std::isinf(timestepFactorIn))
- {
- timestepFactorIn = 1e-100;
- }
- break;
- }
- }
- }//end intersection with neighbor element
-
- // handle boundary face
- if (isIt->boundary())
- {
- // center of face in global coordinates
- GlobalPosition globalPosFace = isIt->geometry().global(faceLocal);
-
- // cell center in global coordinates
- GlobalPosition globalPos = eIt->geometry().global(localPos);
-
- // distance vector between barycenters
- Dune::FieldVector<Scalar, dimWorld> distVec = globalPosFace - globalPos;
- // compute distance between cell centers
- Scalar dist = distVec.two_norm();
-
- Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
- unitDistVec /= dist;
-
- //get boundary type
- BoundaryConditions::Flags bcTypeSat = problem_.bctypeSat(globalPosFace, *isIt);
-
- if (bcTypeSat == BoundaryConditions::dirichlet)
- {
- Scalar satBound = problem_.dirichletSat(globalPosFace, *isIt);
-
- //get velocity*normalvector*facearea/(volume*porosity)
- factor = (problem_.variables().velocity()[globalIdxI][isIndex] * unitDistVec) * (faceArea
- * (unitOuterNormal * unitDistVec)) / (volume * porosity);
- factorSecondPhase = (problem_.variables().velocitySecondPhase()[globalIdxI][isIndex] * unitDistVec)
- * (faceArea * (unitOuterNormal * unitDistVec)) / (volume * porosity);
-
- Scalar pressBound = problem_.variables().pressure()[globalIdxI];
- Scalar temperature = problem_.temperature(globalPosFace, *eIt);
-
- //determine phase saturations from primary saturation variable
- Scalar satWI = 0;
- Scalar satWBound = 0;
- switch (saturationType_)
- {
- case Sw:
- {
- satWI = problem_.variables().saturation()[globalIdxI];
- satWBound = satBound;
- break;
- }
- case Sn:
- {
- satWI = 1 - problem_.variables().saturation()[globalIdxI];
- satWBound = 1 - satBound;
- break;
- }
- }
-
- Scalar pcBound = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(globalPos, *eIt),
- satBound);
-
- //determine phase pressures from primary pressure variable
- Scalar pressW = 0;
- Scalar pressNW = 0;
- switch (pressureType_)
- {
- case pw:
- {
- pressW = pressBound;
- pressNW = pressBound + pcBound;
- break;
- }
- case pn:
- {
- pressW = pressBound - pcBound;
- pressNW = pressBound;
- break;
- }
- }
-
- FluidState fluidState;
- fluidState.update(satBound, pressW, pressNW, temperature);
-
- //get phase potentials
- Scalar potentialW = problem_.variables().potentialWetting(globalIdxI, isIndex);
- Scalar potentialNW = problem_.variables().potentialNonwetting(globalIdxI, isIndex);
-
- Scalar lambdaW = 0;
- Scalar lambdaNW = 0;
-
- //upwinding of lambda dependend on the phase potential gradients
- if (potentialW > 0.)
- {
- lambdaW = lambdaWI;
- if (compressibility_)
- {
- lambdaW /= densityWI;
- }//divide by density because lambda is saved as lambda*density
- }
- else
- {
- if (compressibility_)
- {
- 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;
- }
- }
-
- if (potentialNW >= 0.)
- {
- lambdaNW = lambdaNWI;
- if (compressibility_)
- {
- lambdaNW /= densityNWI;
- }//divide by density because lambda is saved as lambda*density
- }
- else
- {
- if (compressibility_)
- {
- 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;
- }
- }
- // std::cout<<lambdaW<<" "<<lambdaNW<<std::endl;
-
- Scalar krSum = lambdaW * viscosityWI + lambdaNW * viscosityNWI;
-
- switch (velocityType_)
- {
- case vt:
- {
- //for time step criterion
-
- if (factor >= 0)
- {
- timestepFactorOut += factor / (krSum * viscosityRatio);
- }
- if (factor < 0)
- {
- timestepFactorIn -= factor / (krSum * viscosityRatio);
- }
-
- Scalar pcI = problem_.variables().capillaryPressure(globalIdxI);
-
- // calculate the saturation gradient
- Dune::FieldVector<Scalar, dimWorld> pcGradient = unitDistVec;
- 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) * unitDistVec
- * faceArea / (volume * porosity) * (unitOuterNormal * unitDistVec);
- Scalar convPart = convectivePart()(*eIt, isIndex, satWI, satWBound) * unitDistVec * faceArea / (volume * porosity) * (unitOuterNormal * unitDistVec);
-
- //for time step criterion
- if (diffPart >= 0)
- {
- diffFactorIn += diffPart / (krSum * viscosityRatio);
- }
- if (diffPart < 0)
- {
- diffFactorOut -= diffPart / (krSum * viscosityRatio);
- }
- if (convPart >= 0)
- {
- timestepFactorOut += convPart / (krSum * viscosityRatio);
- }
- if (convPart < 0)
- {
- timestepFactorIn -= convPart / (krSum * viscosityRatio);
- }
-
- switch (saturationType_)
- {
- case Sw:
- {
- //vt*fw
- factor *= lambdaW / (lambdaW + lambdaNW);
- break;
- }
- case Sn:
- {
- //vt*fn
- factor *= lambdaNW / (lambdaW + lambdaNW);
- break;
- }
- }
- //vt*fw
- factor -= diffPart;
- factor += convPart;
- break;
- }
-
- case vw:
- {
- if (compressibility_)
- {
- factor /= densityWI;
- factorSecondPhase /= densityNWI;
- }
-
- //for time step criterion
- if (factor >= 0)
- {
- timestepFactorOut += factor / (krSum * viscosityRatio);
- }
- if (factor < 0)
- {
- timestepFactorIn -= factor / (krSum * viscosityRatio);
- }
- if (factorSecondPhase < 0)
- {
- timestepFactorIn -= factorSecondPhase / (krSum * viscosityRatio);
- }
-
- if (std::isnan(timestepFactorIn) || std::isinf(timestepFactorIn))
- {
- timestepFactorIn = 1e-100;
- }
- break;
- }
-
- //for time step criterion if the non-wetting phase velocity is used
- case vn:
- {
- if (compressibility_)
- {
- factor /= densityNWI;
- factorSecondPhase /= densityWI;
- }
-
- //for time step criterion
- if (factor >= 0)
- {
- timestepFactorOut += factor / (krSum * viscosityRatio);
- }
- if (factor < 0)
- {
- timestepFactorIn -= factor / (krSum * viscosityRatio);
- }
- if (factorSecondPhase < 0)
- {
- timestepFactorIn -= factorSecondPhase / (krSum * viscosityRatio);
- }
-
- if (std::isnan(timestepFactorIn) || std::isinf(timestepFactorIn))
- {
- timestepFactorIn = 1e-100;
- }
- break;
- }
- }
- }//end dirichlet boundary
-
- if (bcTypeSat == BoundaryConditions::neumann)
- {
- //get mobilities
- Scalar lambdaW, lambdaNW;
-
- lambdaW = lambdaWI;
- if (compressibility_)
- {
- lambdaW /= densityWI;
- }
-
- lambdaNW = lambdaNWI;
- if (compressibility_)
- {
- lambdaNW /= densityNWI;
- }
-
- Scalar krSum = lambdaW * viscosityWI + lambdaNW * viscosityNWI;
-
- //get velocity*normalvector*facearea/(volume*porosity)
- factor = (problem_.variables().velocity()[globalIdxI][isIndex] * unitOuterNormal) * faceArea
- / (volume * porosity);
- switch (velocityType_)
- {
- case vt:
- {
- switch (saturationType_)
- {
- case Sw:
- {
- //vt*fw
- factor *= lambdaW / (lambdaW + lambdaNW);
- break;
- }
- case Sn:
- {
- //vt*fn
- factor *= lambdaNW / (lambdaW + lambdaNW);
- break;
- }
- }
- break;
- }
- }
- Scalar boundaryFactor = problem_.neumannSat(globalPosFace, *isIt, factor);
-
- if (factor != boundaryFactor)
- {
- switch (saturationType_)
- {
- case Sw:
- {
- factor = boundaryFactor / densityWI * faceArea / (volume * porosity);
- break;
- }
- case Sn:
- {
- factor = boundaryFactor / densityNWI * faceArea / (volume * porosity);
- break;
- }
- }
- }
-
- //for time step criterion
-
- if (factor >= 0)
- {
- timestepFactorOut += factor / (krSum * viscosityRatio);
- }
- if (factor < 0)
- {
- timestepFactorIn -= factor / (krSum * viscosityRatio);
- }
- }//end neumann boundary
- }//end boundary
- // add to update vector
- updateVec[globalIdxI] -= factor;//-:v>0, if flow leaves the cell
- }// end all intersections
- Scalar source = 0;
- switch (velocityType_)
- {
- case vw:
- {
- source = problem_.source(globalPos, *eIt)[wPhaseIdx] / densityWI;
- break;
- }
- case vn:
- {
- source = problem_.source(globalPos, *eIt)[nPhaseIdx] / densityNWI;
- break;
- }
- case vt:
- {
- source = problem_.source(globalPos, *eIt)[wPhaseIdx] / densityWI + problem_.source(globalPos,
- *eIt)[nPhaseIdx] / densityNWI;
- break;
- }
- }
- if (source)
- {
- //get mobilities
- Scalar lambdaW = 0;
- Scalar lambdaNW = 0;
-
- lambdaW = lambdaWI;
- if (compressibility_)
- {
- lambdaW /= densityWI;
- }
- lambdaNW = lambdaNWI;
- if (compressibility_)
- {
- lambdaNW /= densityNWI;
- }
-
- Scalar krSum = lambdaW * viscosityWI + lambdaNW * viscosityNWI;
- switch (saturationType_)
- {
- case Sw:
- {
- updateVec[globalIdxI] += source / porosity;
- break;
- }
- case Sn:
- {
- updateVec[globalIdxI] += source / porosity;
- break;
- }
- }
- if (source >= 0)
- {
- timestepFactorIn += source / (porosity * viscosityRatio * krSum);
- }
- else
- {
- timestepFactorOut -= source / (porosity * viscosityRatio * krSum);
- }
- }
-
- //calculate time step
- if (velocityType_ == vw || velocityType_ == vn)
- {
- dt = std::min(dt, evaluateTimeStepPhaseFlux(timestepFactorIn, timestepFactorOut, residualSatW,
- residualSatNW, globalIdxI));
- }
- if (velocityType_ == vt)
- {
- dt = std::min(dt, evaluateTimeStepTotalFlux(timestepFactorIn, timestepFactorOut, diffFactorIn,
- diffFactorOut, residualSatW, residualSatNW));
- }
-
- problem_.variables().volumecorrection(globalIdxI) = updateVec[globalIdxI];
- } // end grid traversal
-
- return 0;
+ Scalar residualSatW = problem_.spatialParameters().materialLawParams(globalPos, *eIt).Swr();
+ Scalar residualSatNW = problem_.spatialParameters().materialLawParams(globalPos, *eIt).Snr();
+
+ //for benchmark only!
+ // problem_.variables().storeSrn(residualSatNW, globalIdxI);
+ //for benchmark only!
+
+ Scalar porosity = problem_.spatialParameters().porosity(globalPos, *eIt);
+
+ Scalar viscosityWI = problem_.variables().viscosityWetting(globalIdxI);
+ Scalar viscosityNWI = problem_.variables().viscosityNonwetting(globalIdxI);
+ Scalar viscosityRatio = 1 - fabs(0.5 - viscosityNWI / (viscosityWI + viscosityNWI));//1 - fabs(viscosityWI-viscosityNWI)/(viscosityWI+viscosityNWI);
+
+ Scalar densityWI = problem_.variables().densityWetting(globalIdxI);
+ Scalar densityNWI = problem_.variables().densityNonwetting(globalIdxI);
+
+ Scalar lambdaWI = problem_.variables().mobilityWetting(globalIdxI);
+ Scalar lambdaNWI = problem_.variables().mobilityNonwetting(globalIdxI);
+
+ Scalar timestepFactorIn = 0;
+ Scalar timestepFactorOut = 0;
+ Scalar diffFactorIn = 0;
+ Scalar diffFactorOut = 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)
+ {
+ // local number of facet
+ int isIndex = isIt->indexInInside();
+
+ // get geometry type of face
+ Dune::GeometryType faceGT = isIt->geometryInInside().type();
+
+ // center in face's reference element
+ const Dune::FieldVector<Scalar, dim - 1>& faceLocal =
+ ReferenceElementFaceContainer::general(faceGT).position(0, 0);
+
+ // center of face inside volume reference element
+ const LocalPosition localPosFace(0);
+
+ Dune::FieldVector<Scalar, dimWorld> unitOuterNormal = isIt->unitOuterNormal(faceLocal);
+ if (switchNormals_)
+ unitOuterNormal *= -1.0;
+
+ Scalar faceArea = isIt->geometry().volume();
+
+ Scalar factor = 0;
+ Scalar factorSecondPhase = 0;
+
+ // handle interior face
+ if (isIt->neighbor())
+ {
+ // access neighbor
+ ElementPointer neighborPointer = isIt->outside();
+ int globalIdxJ = problem_.variables().index(*neighborPointer);
+
+ // compute factor in neighbor
+ Dune::GeometryType neighborGT = neighborPointer->geometry().type();
+ const LocalPosition& localPosNeighbor = ReferenceElementContainer::general(neighborGT).position(0, 0);
+
+ // cell center in global coordinates
+ const GlobalPosition& globalPos = eIt->geometry().global(localPos);
+
+ // neighbor cell center in global coordinates
+ const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().global(localPosNeighbor);
+
+ // distance vector between barycenters
+ Dune::FieldVector<Scalar, dimWorld> distVec = globalPosNeighbor - globalPos;
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
+
+ Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ //get phase potentials
+ 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);
+
+ //get velocity*normalvector*facearea/(volume*porosity)
+ factor = (problem_.variables().velocity()[globalIdxI][isIndex] * unitOuterNormal)
+ * (faceArea) / (volume * porosity);
+ factorSecondPhase = (problem_.variables().velocitySecondPhase()[globalIdxI][isIndex] * unitOuterNormal)
+ * (faceArea) / (volume * porosity);
+
+ Scalar lambdaW = 0;
+ Scalar lambdaNW = 0;
+ Scalar viscosityW = 0;
+ Scalar viscosityNW = 0;
+
+ //upwinding of lambda dependend on the phase potential gradients
+ if (potentialW >= 0.)
+ {
+ lambdaW = lambdaWI;
+ if (compressibility_)
+ {
+ lambdaW /= densityWI;
+ }//divide by density because lambda is saved as lambda*density
+ viscosityW = viscosityWI;
+ }
+ else
+ {
+ lambdaW = problem_.variables().mobilityWetting(globalIdxJ);
+ if (compressibility_)
+ {
+ lambdaW /= densityWJ;
+ }//divide by density because lambda is saved as lambda*density
+ viscosityW = problem_.variables().viscosityWetting(globalIdxJ);
+ }
+
+ if (potentialNW >= 0.)
+ {
+ lambdaNW = lambdaNWI;
+ if (compressibility_)
+ {
+ lambdaNW /= densityNWI;
+ }//divide by density because lambda is saved as lambda*density
+ viscosityNW = viscosityNWI;
+ }
+ else
+ {
+ lambdaNW = problem_.variables().mobilityNonwetting(globalIdxJ);
+ if (compressibility_)
+ {
+ lambdaNW /= densityNWJ;
+ }//divide by density because lambda is saved as lambda*density
+ viscosityNW = problem_.variables().viscosityNonwetting(globalIdxJ);
+ }
+ Scalar krSum = lambdaW * viscosityW + lambdaNW * viscosityNW;
+
+ switch (velocityType_)
+ {
+ case vt:
+ {
+ //for time step criterion
+
+ if (factor >= 0)
+ {
+ timestepFactorOut += factor / (krSum * viscosityRatio);
+ }
+ if (factor < 0)
+ {
+ timestepFactorIn -= factor / (krSum * viscosityRatio);
+ }
+
+ //determine phase saturations from primary saturation variable
+ Scalar satWI = 0;
+ Scalar satWJ = 0;
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ satWI = problem_.variables().saturation()[globalIdxI];
+ satWJ = problem_.variables().saturation()[globalIdxJ];
+ break;
+ }
+ case Sn:
+ {
+ 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);
+
+ // calculate the saturation gradient
+ Dune::FieldVector<Scalar, dimWorld> pcGradient = unitDistVec;
+ 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) * unitDistVec * faceArea
+ / (volume * porosity) * (unitOuterNormal * unitDistVec);
+ Scalar convPart = convectivePart() (*eIt, isIndex, satWI, satWJ) * unitDistVec * faceArea / (volume * porosity) * (unitOuterNormal * unitDistVec);
+
+ //for time step criterion
+ if (diffPart >= 0)
+ {
+ diffFactorIn += diffPart / (krSum * viscosityRatio);
+ }
+ if (diffPart < 0)
+ {
+ diffFactorOut -= diffPart / (krSum * viscosityRatio);
+ }
+ if (convPart >= 0)
+ {
+ timestepFactorOut += convPart / (krSum * viscosityRatio);
+ }
+ if (convPart < 0)
+ {
+ timestepFactorIn -= convPart / (krSum * viscosityRatio);
+ }
+
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ //vt*fw
+ factor *= lambdaW / (lambdaW + lambdaNW);
+ break;
+ }
+ case Sn:
+ {
+ //vt*fn
+ factor *= lambdaNW / (lambdaW + lambdaNW);
+ break;
+ }
+ }
+ factor -= diffPart;
+ factor += convPart;
+ break;
+ }
+ case vw:
+ {
+ if (compressibility_)
+ {
+ factor /= densityWI;
+ factorSecondPhase /= densityNWI;
+ }
+
+ //for time step criterion
+ if (factor >= 0)
+ {
+ timestepFactorOut += factor / (krSum * viscosityRatio);
+ }
+ if (factor < 0)
+ {
+ timestepFactorIn -= factor / (krSum * viscosityRatio);
+ }
+ if (factorSecondPhase < 0)
+ {
+ timestepFactorIn -= factorSecondPhase / (krSum * viscosityRatio);
+ }
+
+ if (std::isnan(timestepFactorIn) || std::isinf(timestepFactorIn))
+ {
+ timestepFactorIn = 1e-100;
+ }
+ break;
+ }
+
+ //for time step criterion if the non-wetting phase velocity is used
+ case vn:
+ {
+ if (compressibility_)
+ {
+ factor /= densityNWI;
+ factorSecondPhase /= densityWI;
+ }
+
+ //for time step criterion
+ if (factor >= 0)
+ {
+ timestepFactorOut += factor / (krSum * viscosityRatio);
+ }
+ if (factor < 0)
+ {
+ timestepFactorIn -= factor / (krSum * viscosityRatio);
+ }
+ if (factorSecondPhase < 0)
+ {
+ timestepFactorIn -= factorSecondPhase / (krSum * viscosityRatio);
+ }
+
+ if (std::isnan(timestepFactorIn) || std::isinf(timestepFactorIn))
+ {
+ timestepFactorIn = 1e-100;
+ }
+ break;
+ }
+ }
+ }//end intersection with neighbor element
+
+ // handle boundary face
+ if (isIt->boundary())
+ {
+ // center of face in global coordinates
+ GlobalPosition globalPosFace = isIt->geometry().global(faceLocal);
+
+ // cell center in global coordinates
+ GlobalPosition globalPos = eIt->geometry().global(localPos);
+
+ // distance vector between barycenters
+ Dune::FieldVector<Scalar, dimWorld> distVec = globalPosFace - globalPos;
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
+
+ Dune::FieldVector<Scalar, dimWorld> unitDistVec(distVec);
+ unitDistVec /= dist;
+
+ //get boundary type
+ BoundaryConditions::Flags bcTypeSat = problem_.bctypeSat(globalPosFace, *isIt);
+
+ if (bcTypeSat == BoundaryConditions::dirichlet)
+ {
+ Scalar satBound = problem_.dirichletSat(globalPosFace, *isIt);
+
+ //get velocity*normalvector*facearea/(volume*porosity)
+ factor = (problem_.variables().velocity()[globalIdxI][isIndex] * unitOuterNormal)
+ * (faceArea) / (volume * porosity);
+ factorSecondPhase = (problem_.variables().velocitySecondPhase()[globalIdxI][isIndex] * unitOuterNormal)
+ * (faceArea) / (volume * porosity);
+
+ Scalar pressBound = problem_.variables().pressure()[globalIdxI];
+ Scalar temperature = problem_.temperature(globalPosFace, *eIt);
+
+ //determine phase saturations from primary saturation variable
+ Scalar satWI = 0;
+ Scalar satWBound = 0;
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ satWI = problem_.variables().saturation()[globalIdxI];
+ satWBound = satBound;
+ break;
+ }
+ case Sn:
+ {
+ satWI = 1 - problem_.variables().saturation()[globalIdxI];
+ satWBound = 1 - satBound;
+ break;
+ }
+ }
+
+ Scalar pcBound = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(globalPos, *eIt),
+ satBound);
+
+ //determine phase pressures from primary pressure variable
+ Scalar pressW = 0;
+ Scalar pressNW = 0;
+ switch (pressureType_)
+ {
+ case pw:
+ {
+ pressW = pressBound;
+ pressNW = pressBound + pcBound;
+ break;
+ }
+ case pn:
+ {
+ pressW = pressBound - pcBound;
+ pressNW = pressBound;
+ break;
+ }
+ }
+
+ FluidState fluidState;
+ fluidState.update(satBound, pressW, pressNW, temperature);
+
+ //get phase potentials
+ Scalar potentialW = problem_.variables().potentialWetting(globalIdxI, isIndex);
+ Scalar potentialNW = problem_.variables().potentialNonwetting(globalIdxI, isIndex);
+
+ Scalar lambdaW = 0;
+ Scalar lambdaNW = 0;
+
+ //upwinding of lambda dependend on the phase potential gradients
+ if (potentialW > 0.)
+ {
+ lambdaW = lambdaWI;
+ if (compressibility_)
+ {
+ lambdaW /= densityWI;
+ }//divide by density because lambda is saved as lambda*density
+ }
+ else
+ {
+ if (compressibility_)
+ {
+ 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;
+ }
+ }
+
+ if (potentialNW >= 0.)
+ {
+ lambdaNW = lambdaNWI;
+ if (compressibility_)
+ {
+ lambdaNW /= densityNWI;
+ }//divide by density because lambda is saved as lambda*density
+ }
+ else
+ {
+ if (compressibility_)
+ {
+ 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;
+ }
+ }
+ // std::cout<<lambdaW<<" "<<lambdaNW<<std::endl;
+
+ Scalar krSum = lambdaW * viscosityWI + lambdaNW * viscosityNWI;
+
+ switch (velocityType_)
+ {
+ case vt:
+ {
+ //for time step criterion
+
+ if (factor >= 0)
+ {
+ timestepFactorOut += factor / (krSum * viscosityRatio);
+ }
+ if (factor < 0)
+ {
+ timestepFactorIn -= factor / (krSum * viscosityRatio);
+ }
+
+ Scalar pcI = problem_.variables().capillaryPressure(globalIdxI);
+
+ // calculate the saturation gradient
+ Dune::FieldVector<Scalar, dimWorld> pcGradient = unitDistVec;
+ 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) * unitDistVec
+ * faceArea / (volume * porosity) * (unitOuterNormal * unitDistVec);
+ Scalar convPart = convectivePart()(*eIt, isIndex, satWI, satWBound) * unitDistVec * faceArea / (volume * porosity) * (unitOuterNormal * unitDistVec);
+
+ //for time step criterion
+ if (diffPart >= 0)
+ {
+ diffFactorIn += diffPart / (krSum * viscosityRatio);
+ }
+ if (diffPart < 0)
+ {
+ diffFactorOut -= diffPart / (krSum * viscosityRatio);
+ }
+ if (convPart >= 0)
+ {
+ timestepFactorOut += convPart / (krSum * viscosityRatio);
+ }
+ if (convPart < 0)
+ {
+ timestepFactorIn -= convPart / (krSum * viscosityRatio);
+ }
+
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ //vt*fw
+ factor *= lambdaW / (lambdaW + lambdaNW);
+ break;
+ }
+ case Sn:
+ {
+ //vt*fn
+ factor *= lambdaNW / (lambdaW + lambdaNW);
+ break;
+ }
+ }
+ //vt*fw
+ factor -= diffPart;
+ factor += convPart;
+ break;
+ }
+
+ case vw:
+ {
+ if (compressibility_)
+ {
+ factor /= densityWI;
+ factorSecondPhase /= densityNWI;
+ }
+
+ //for time step criterion
+ if (factor >= 0)
+ {
+ timestepFactorOut += factor / (krSum * viscosityRatio);
+ }
+ if (factor < 0)
+ {
+ timestepFactorIn -= factor / (krSum * viscosityRatio);
+ }
+ if (factorSecondPhase < 0)
+ {
+ timestepFactorIn -= factorSecondPhase / (krSum * viscosityRatio);
+ }
+
+ if (std::isnan(timestepFactorIn) || std::isinf(timestepFactorIn))
+ {
+ timestepFactorIn = 1e-100;
+ }
+ break;
+ }
+
+ //for time step criterion if the non-wetting phase velocity is used
+ case vn:
+ {
+ if (compressibility_)
+ {
+ factor /= densityNWI;
+ factorSecondPhase /= densityWI;
+ }
+
+ //for time step criterion
+ if (factor >= 0)
+ {
+ timestepFactorOut += factor / (krSum * viscosityRatio);
+ }
+ if (factor < 0)
+ {
+ timestepFactorIn -= factor / (krSum * viscosityRatio);
+ }
+ if (factorSecondPhase < 0)
+ {
+ timestepFactorIn -= factorSecondPhase / (krSum * viscosityRatio);
+ }
+
+ if (std::isnan(timestepFactorIn) || std::isinf(timestepFactorIn))
+ {
+ timestepFactorIn = 1e-100;
+ }
+ break;
+ }
+ }
+ }//end dirichlet boundary
+
+ if (bcTypeSat == BoundaryConditions::neumann)
+ {
+ //get mobilities
+ Scalar lambdaW, lambdaNW;
+
+ lambdaW = lambdaWI;
+ if (compressibility_)
+ {
+ lambdaW /= densityWI;
+ }
+
+ lambdaNW = lambdaNWI;
+ if (compressibility_)
+ {
+ lambdaNW /= densityNWI;
+ }
+
+ Scalar krSum = lambdaW * viscosityWI + lambdaNW * viscosityNWI;
+
+ //get velocity*normalvector*facearea/(volume*porosity)
+ factor = (problem_.variables().velocity()[globalIdxI][isIndex] * unitOuterNormal) * faceArea
+ / (volume * porosity);
+ switch (velocityType_)
+ {
+ case vt:
+ {
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ //vt*fw
+ factor *= lambdaW / (lambdaW + lambdaNW);
+ break;
+ }
+ case Sn:
+ {
+ //vt*fn
+ factor *= lambdaNW / (lambdaW + lambdaNW);
+ break;
+ }
+ }
+ break;
+ }
+ }
+ Scalar boundaryFactor = problem_.neumannSat(globalPosFace, *isIt, factor);
+
+ if (factor != boundaryFactor)
+ {
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ factor = boundaryFactor / densityWI * faceArea / (volume * porosity);
+ break;
+ }
+ case Sn:
+ {
+ factor = boundaryFactor / densityNWI * faceArea / (volume * porosity);
+ break;
+ }
+ }
+ }
+
+ //for time step criterion
+
+ if (factor >= 0)
+ {
+ timestepFactorOut += factor / (krSum * viscosityRatio);
+ }
+ if (factor < 0)
+ {
+ timestepFactorIn -= factor / (krSum * viscosityRatio);
+ }
+ }//end neumann boundary
+ }//end boundary
+ // add to update vector
+ updateVec[globalIdxI] -= factor;//-:v>0, if flow leaves the cell
+ }// end all intersections
+ Scalar source = 0;
+ switch (velocityType_)
+ {
+ case vw:
+ {
+ source = problem_.source(globalPos, *eIt)[wPhaseIdx] / densityWI;
+ break;
+ }
+ case vn:
+ {
+ source = problem_.source(globalPos, *eIt)[nPhaseIdx] / densityNWI;
+ break;
+ }
+ case vt:
+ {
+ source = problem_.source(globalPos, *eIt)[wPhaseIdx] / densityWI + problem_.source(globalPos,
+ *eIt)[nPhaseIdx] / densityNWI;
+ break;
+ }
+ }
+ if (source)
+ {
+ //get mobilities
+ Scalar lambdaW = 0;
+ Scalar lambdaNW = 0;
+
+ lambdaW = lambdaWI;
+ if (compressibility_)
+ {
+ lambdaW /= densityWI;
+ }
+ lambdaNW = lambdaNWI;
+ if (compressibility_)
+ {
+ lambdaNW /= densityNWI;
+ }
+
+ Scalar krSum = lambdaW * viscosityWI + lambdaNW * viscosityNWI;
+ switch (saturationType_)
+ {
+ case Sw:
+ {
+ updateVec[globalIdxI] += source / porosity;
+ break;
+ }
+ case Sn:
+ {
+ updateVec[globalIdxI] += source / porosity;
+ break;
+ }
+ }
+ if (source >= 0)
+ {
+ timestepFactorIn += source / (porosity * viscosityRatio * krSum);
+ }
+ else
+ {
+ timestepFactorOut -= source / (porosity * viscosityRatio * krSum);
+ }
+ }
+
+ //calculate time step
+ if (velocityType_ == vw || velocityType_ == vn)
+ {
+ dt = std::min(dt, evaluateTimeStepPhaseFlux(timestepFactorIn, timestepFactorOut, residualSatW,
+ residualSatNW, globalIdxI));
+ }
+ if (velocityType_ == vt)
+ {
+ dt = std::min(dt, evaluateTimeStepTotalFlux(timestepFactorIn, timestepFactorOut, diffFactorIn,
+ diffFactorOut, residualSatW, residualSatNW));
+ }
+
+ problem_.variables().volumecorrection(globalIdxI) = updateVec[globalIdxI];
+ } // end grid traversal
+
+ return 0;
}
template<class TypeTag>
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)
- {
- // get geometry type
- Dune::GeometryType gt = eIt->geometry().type();
-
- // get cell center in reference element
- const LocalPosition &localPos = ReferenceElementContainer::general(gt).position(0, 0);
-
- // get global coordinate of cell center
- GlobalPosition globalPos = eIt->geometry().global(localPos);
-
- // initialize cell concentration
- problem_.variables().saturation()[problem_.variables().index(*eIt)] = problem_.initSat(globalPos, *eIt);
- }
- return;
+ // 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)
+ {
+ // get geometry type
+ Dune::GeometryType gt = eIt->geometry().type();
+
+ // get cell center in reference element
+ const LocalPosition &localPos = ReferenceElementContainer::general(gt).position(0, 0);
+
+ // get global coordinate of cell center
+ GlobalPosition globalPos = eIt->geometry().global(localPos);
+
+ // initialize cell concentration
+ problem_.variables().saturation()[problem_.variables().index(*eIt)] = problem_.initSat(globalPos, *eIt);
+ }
+ return;
}
template<class TypeTag>
typename FVSaturation2P<TypeTag>::Scalar FVSaturation2P<TypeTag>::evaluateTimeStepTotalFlux(Scalar timestepFactorIn,
- Scalar timestepFactorOut, Scalar diffFactorIn, Scalar diffFactorOut, Scalar& residualSatW,
- Scalar& residualSatNW)
-{
- // compute volume correction
- Scalar volumeCorrectionFactor = (1 - residualSatW - residualSatNW);
-
- //make sure correction is in the right range. If not: force dt to be not min-dt!
- if (timestepFactorIn <= 0)
- {
- timestepFactorIn = 1e-100;
- }
- if (timestepFactorOut <= 0)
- {
- timestepFactorOut = 1e-100;
- }
-
- Scalar sumFactor = std::min(volumeCorrectionFactor / timestepFactorIn, volumeCorrectionFactor / timestepFactorOut);
-
- //make sure that diffFactor > 0
- if (diffFactorIn <= 0)
- {
- diffFactorIn = 1e-100;
- }
- if (diffFactorOut <= 0)
- {
- diffFactorOut = 1e-100;
- }
-
- Scalar minDiff = std::min(volumeCorrectionFactor / diffFactorIn, volumeCorrectionFactor / diffFactorOut);
-
- //determine time step
- sumFactor = std::min(sumFactor, 0.1 * minDiff);
-
- return sumFactor;
-}
+ Scalar timestepFactorOut, Scalar diffFactorIn, Scalar diffFactorOut, Scalar& residualSatW,
+ Scalar& residualSatNW)
+ {
+ // compute volume correction
+ Scalar volumeCorrectionFactor = (1 - residualSatW - residualSatNW);
+
+ //make sure correction is in the right range. If not: force dt to be not min-dt!
+ if (timestepFactorIn <= 0)
+ {
+ timestepFactorIn = 1e-100;
+ }
+ if (timestepFactorOut <= 0)
+ {
+ timestepFactorOut = 1e-100;
+ }
+
+ Scalar sumFactor = std::min(volumeCorrectionFactor / timestepFactorIn, volumeCorrectionFactor / timestepFactorOut);
+
+ //make sure that diffFactor > 0
+ if (diffFactorIn <= 0)
+ {
+ diffFactorIn = 1e-100;
+ }
+ if (diffFactorOut <= 0)
+ {
+ diffFactorOut = 1e-100;
+ }
+
+ Scalar minDiff = std::min(volumeCorrectionFactor / diffFactorIn, volumeCorrectionFactor / diffFactorOut);
+
+ //determine time step
+ sumFactor = std::min(sumFactor, 0.1 * minDiff);
+
+ return sumFactor;
+ }
template<class TypeTag>
typename FVSaturation2P<TypeTag>::Scalar FVSaturation2P<TypeTag>::evaluateTimeStepPhaseFlux(
- Scalar timestepFactorIn, Scalar timestepFactorOut, Scalar& residualSatW, Scalar& residualSatNW, int globalIdxI)
-{
- // compute dt restriction
- Scalar volumeCorrectionFactorIn = 1 - residualSatW - residualSatNW;
- Scalar volumeCorrectionFactorOut = 0;
- if (saturationType_ == Sw)
- {
- Scalar satI = problem_.variables().saturation()[globalIdxI];
- volumeCorrectionFactorOut = std::max((satI - residualSatW), 1e-2);
- }
- if (saturationType_ == Sn)
- {
- Scalar satI = problem_.variables().saturation()[globalIdxI];
- volumeCorrectionFactorOut = std::max((satI - residualSatNW), 1e-2);
- }
-
- //make sure correction is in the right range. If not: force dt to be not min-dt!
- if (volumeCorrectionFactorOut <= 0)
- {
- volumeCorrectionFactorOut = 1e100;
- }
-
- //make sure correction is in the right range. If not: force dt to be not min-dt!
- if (timestepFactorIn <= 0)
- {
- timestepFactorIn = 1e-100;
- }
- if (timestepFactorOut <= 0)
- {
- timestepFactorOut = 1e-100;
- }
-
- //correct volume
- timestepFactorIn = volumeCorrectionFactorIn / timestepFactorIn;
- timestepFactorOut = volumeCorrectionFactorOut / timestepFactorOut;
-
- //determine timestep
- Scalar timestepFactor = std::min(timestepFactorIn, timestepFactorOut);
-
- return timestepFactor;
-}
+ Scalar timestepFactorIn, Scalar timestepFactorOut, Scalar& residualSatW, Scalar& residualSatNW, int globalIdxI)
+ {
+ // compute dt restriction
+ Scalar volumeCorrectionFactorIn = 1 - residualSatW - residualSatNW;
+ Scalar volumeCorrectionFactorOut = 0;
+ if (saturationType_ == Sw)
+ {
+ Scalar satI = problem_.variables().saturation()[globalIdxI];
+ volumeCorrectionFactorOut = std::max((satI - residualSatW), 1e-2);
+ }
+ if (saturationType_ == Sn)
+ {
+ Scalar satI = problem_.variables().saturation()[globalIdxI];
+ volumeCorrectionFactorOut = std::max((satI - residualSatNW), 1e-2);
+ }
+
+ //make sure correction is in the right range. If not: force dt to be not min-dt!
+ if (volumeCorrectionFactorOut <= 0)
+ {
+ volumeCorrectionFactorOut = 1e100;
+ }
+
+ //make sure correction is in the right range. If not: force dt to be not min-dt!
+ if (timestepFactorIn <= 0)
+ {
+ timestepFactorIn = 1e-100;
+ }
+ if (timestepFactorOut <= 0)
+ {
+ timestepFactorOut = 1e-100;
+ }
+
+ //correct volume
+ timestepFactorIn = volumeCorrectionFactorIn / timestepFactorIn;
+ timestepFactorOut = volumeCorrectionFactorOut / timestepFactorOut;
+
+ //determine timestep
+ Scalar timestepFactor = std::min(timestepFactorIn, timestepFactorOut);
+
+ return timestepFactor;
+ }
template<class TypeTag>
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)
- {
- // get geometry type
- Dune::GeometryType gt = eIt->geometry().type();
-
- // get cell center in reference element
- const LocalPosition &localPos = ReferenceElementContainer::general(gt).position(0, 0);
-
- // get global coordinate of cell center
- GlobalPosition globalPos = eIt->geometry().global(localPos);
-
- int globalIdx = problem_.variables().index(*eIt);
-
- Scalar sat = 0;
- if (!iterate)
- {
- sat = problem_.variables().saturation()[globalIdx];
- }
- else
- {
- sat = saturation[globalIdx];
- }
-
- //determine phase saturations from primary saturation variable
- Scalar satW = 0;
- if (saturationType_ == Sw)
- {
- satW = sat;
- }
- if (saturationType_ == Sn)
- {
- satW = 1 - sat;
- }
-
- 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);
-
- // 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().fracFlowFuncWetting(globalIdx) = problem_.variables().densityWetting(globalIdx)
- * problem_.variables().mobilityWetting(globalIdx) / (problem_.variables().densityWetting(globalIdx)
- * problem_.variables().mobilityWetting(globalIdx) + problem_.variables().densityNonwetting(globalIdx)
- * problem_.variables().mobilityNonwetting(globalIdx));
- problem_.variables().fracFlowFuncNonwetting(globalIdx) = problem_.variables().densityNonwetting(globalIdx)
- * problem_.variables().mobilityNonwetting(globalIdx) / (problem_.variables().densityWetting(globalIdx)
- * problem_.variables().mobilityWetting(globalIdx) + problem_.variables().densityNonwetting(globalIdx)
- * problem_.variables().mobilityNonwetting(globalIdx));
-
- problem_.spatialParameters().update(satW, *eIt);
- }
- return;
-}
+ 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)
+ {
+ // get geometry type
+ Dune::GeometryType gt = eIt->geometry().type();
+
+ // get cell center in reference element
+ const LocalPosition &localPos = ReferenceElementContainer::general(gt).position(0, 0);
+
+ // get global coordinate of cell center
+ GlobalPosition globalPos = eIt->geometry().global(localPos);
+
+ int globalIdx = problem_.variables().index(*eIt);
+
+ Scalar sat = 0;
+ if (!iterate)
+ {
+ sat = problem_.variables().saturation()[globalIdx];
+ }
+ else
+ {
+ sat = saturation[globalIdx];
+ }
+
+ //determine phase saturations from primary saturation variable
+ Scalar satW = 0;
+ if (saturationType_ == Sw)
+ {
+ satW = sat;
+ }
+ if (saturationType_ == Sn)
+ {
+ satW = 1 - sat;
+ }
+
+ 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);
+
+ // 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().fracFlowFuncWetting(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_.spatialParameters().update(satW, *eIt);
+ }
+ return;
+ }
}
#endif
diff --git a/dumux/decoupled/2p/transport/transportproperties.hh b/dumux/decoupled/2p/transport/transportproperties.hh
index 0d7c4c1722..b5d2aaff09 100644
--- a/dumux/decoupled/2p/transport/transportproperties.hh
+++ b/dumux/decoupled/2p/transport/transportproperties.hh
@@ -31,6 +31,18 @@ class DiffusivePart;
template<class TypeTag>
class ConvectivePart;
+template<class TypeTag>
+class TwoPCommonIndices;
+
+template<class TypeTag>
+class FluidSystem2P;
+
+template<class TypeTag>
+class TwoPFluidState;
+
+template<class TypeTag>
+class VariableClass2P;
+
namespace Properties
{
/*!
@@ -51,10 +63,31 @@ NEW_TYPE_TAG(Transport);
NEW_PROP_TAG( DiffusivePart ); //!< The type of the diffusive part in a transport equation
NEW_PROP_TAG( ConvectivePart ); //!< The type of a convective part in a transport equation
+NEW_PROP_TAG( Variables );
+NEW_PROP_TAG( NumPhases );
+NEW_PROP_TAG( NumComponents );
+NEW_PROP_TAG( TwoPIndices );
+NEW_PROP_TAG( FluidSystem );
+NEW_PROP_TAG( FluidState );
+NEW_PROP_TAG( EnableCompressibility );
+NEW_PROP_TAG( PressureFormulation );
+NEW_PROP_TAG( SaturationFormulation );
+NEW_PROP_TAG( VelocityFormulation );
+NEW_PROP_TAG( CFLFactor );
SET_TYPE_PROP(Transport, DiffusivePart, DiffusivePart<TypeTag>);
SET_TYPE_PROP(Transport, ConvectivePart, ConvectivePart<TypeTag>);
-
+SET_TYPE_PROP(Transport, Variables, VariableClass2P<TypeTag>);
+SET_INT_PROP(Transport, NumPhases, 2);
+SET_INT_PROP(Transport, NumComponents, 1);
+SET_TYPE_PROP(Transport, TwoPIndices, TwoPCommonIndices<TypeTag>);
+SET_TYPE_PROP(Transport, FluidSystem, FluidSystem2P<TypeTag>);
+SET_TYPE_PROP(Transport, FluidState, TwoPFluidState<TypeTag>);
+SET_BOOL_PROP(Transport, EnableCompressibility, false);
+SET_INT_PROP(Transport, PressureFormulation, TwoPCommonIndices<TypeTag>::pressureW);
+SET_INT_PROP(Transport, SaturationFormulation, TwoPCommonIndices<TypeTag>::saturationW);
+SET_INT_PROP(Transport, VelocityFormulation, TwoPCommonIndices<TypeTag>::velocityTotal);
+SET_SCALAR_PROP(Transport, CFLFactor, 1.0);
}
}
diff --git a/dumux/decoupled/2p/variableclass2p.hh b/dumux/decoupled/2p/variableclass2p.hh
index 0ab82fe317..63588d9641 100644
--- a/dumux/decoupled/2p/variableclass2p.hh
+++ b/dumux/decoupled/2p/variableclass2p.hh
@@ -82,7 +82,7 @@ public:
typedef typename SolutionTypes::PhasePropertyElemFace PhasePropertyElemFaceType;//!<type for vector of vectors (of size 2 x dimension) of scalars
typedef typename SolutionTypes::DimVecElemFace DimVecElemFaceType;//!<type for vector of vectors (of size 2 x dimension) of vector (of size dimension) of scalars
-private:
+public:
const int codim_;
ScalarSolutionType saturation_;
diff --git a/dumux/decoupled/common/onemodelproblem.hh b/dumux/decoupled/common/onemodelproblem.hh
index cf59b35378..3265a85711 100644
--- a/dumux/decoupled/common/onemodelproblem.hh
+++ b/dumux/decoupled/common/onemodelproblem.hh
@@ -33,6 +33,11 @@
namespace Dumux
{
+namespace Properties
+{
+NEW_PROP_TAG(CFLFactor);
+}
+
/*! \ingroup diffusion
* @brief base class that defines the parameters of loosely coupled diffusion and transport equations
*
@@ -67,6 +72,9 @@ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP(TypeTag, PTAG(SolutionTypes)) SolutionTypes;
+ typedef typename SolutionTypes::ScalarSolution Solution;
+
enum
{
dim = GridView::dimension,
@@ -108,6 +116,7 @@ public:
// bboxMax_[i] = std::max(bboxMax_[i], vIt->geometry().corner(0)[i]);
// }
// }
+ cFLFactor_ = GET_PROP_VALUE(TypeTag, PTAG(CFLFactor));
model_ = new Model(asImp_()) ;
}
@@ -126,7 +135,7 @@ public:
void init()
{
// set the initial condition of the model
- model().initial();
+ model().initialize();
}
/*!
@@ -150,7 +159,6 @@ public:
void timeIntegration()
{
// allocate temporary vectors for the updates
- typedef typename Model::SolutionType Solution;
Solution k1 = asImp_().variables().saturation();
dt_ = 1e100;
@@ -160,7 +168,7 @@ public:
model().update(t, dt_, k1);
//make sure t_old + dt is not larger than tend
- dt_ = std::min(dt_, timeManager().episodeMaxTimeStepSize());
+ dt_ = std::min(dt_*cFLFactor_, timeManager().episodeMaxTimeStepSize());
timeManager().setTimeStepSize(dt_);
// explicit Euler: Sat <- Sat + dt*N(Sat)
@@ -226,7 +234,7 @@ public:
if (gridView().comm().rank() == 0)
std::cout << "Writing result file for current time step\n";
- resultWriter_.beginTimestep(timeManager_.time(), gridView());
+ resultWriter_.beginTimestep(timeManager_.time() + timeManager_.timeStepSize(), gridView());
model().addOutputVtkFields(resultWriter_);
resultWriter_.endTimestep();
}
@@ -338,7 +346,7 @@ public:
typedef Dumux::Restart Restarter;
Restarter res;
- res.serializeBegin(*asImp_());
+ res.serializeBegin(asImp_());
std::cerr << "Serialize to file " << res.fileName() << "\n";
timeManager_.serialize(res);
@@ -405,6 +413,8 @@ private:
Scalar dt_;
+ Scalar cFLFactor_;
+
Model* model_;
VtkMultiWriter resultWriter_;
diff --git a/dumux/decoupled/common/variableclass.hh b/dumux/decoupled/common/variableclass.hh
index 7dc0e3adaf..3b4bf93d3e 100644
--- a/dumux/decoupled/common/variableclass.hh
+++ b/dumux/decoupled/common/variableclass.hh
@@ -141,19 +141,6 @@ public:
instream >> pressure_[globalIdx];
}
-private:
- void initializeGlobalVariables(Dune::FieldVector<Scalar, dim>& initialVel)
- {
- //resize to grid size
- pressure_.resize(gridSize_);
- velocity_.resize(gridSize_);//depends on pressure
- potential_.resize(gridSize_);//depends on pressure
-
- //initialise variables
- pressure_ = 0;
- velocity_ = initialVel;
- }
-
void initializePotentials(Dune::FieldVector<Scalar, dim>& initialVel)
{
if (initialVel.two_norm())
@@ -185,6 +172,50 @@ private:
return;
}
+private:
+ void initializeGlobalVariables(Dune::FieldVector<Scalar, dim>& initialVel)
+ {
+ //resize to grid size
+ pressure_.resize(gridSize_);
+ velocity_.resize(gridSize_);//depends on pressure
+ potential_.resize(gridSize_);//depends on pressure
+
+ //initialise variables
+ pressure_ = 0;
+ velocity_ = initialVel;
+ }
+
+// void initializePotentials(Dune::FieldVector<Scalar, dim>& initialVel)
+// {
+// if (initialVel.two_norm())
+// {
+// // compute update vector
+// ElementIterator eItEnd = gridView_.template end<0> ();
+// for (ElementIterator eIt = gridView_.template begin<0> (); eIt != eItEnd; ++eIt)
+// {
+// // cell index
+// int globalIdxI = elementMapper_.map(*eIt);
+//
+// // run through all intersections with neighbors and boundary
+// IntersectionIterator isItEnd = gridView_.iend(*eIt);
+// for (IntersectionIterator isIt = gridView_.ibegin(*eIt); isIt != isItEnd; ++isIt)
+// {
+// // local number of facet
+// int indexInInside = isIt->indexInInside();
+//
+// Dune::FieldVector<Scalar, dimWorld> unitOuterNormal = isIt->centerUnitOuterNormal();
+//
+// for (int i = 0; i < numPhase; i++) {potential_[globalIdxI][indexInInside][i] = initialVel * unitOuterNormal;}
+// }
+// }
+// }
+// else
+// {
+// potential_ = Dune::FieldVector<Scalar, numPhase> (0);
+// }
+// return;
+// }
+
//Write saturation and pressure into file
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
diff --git a/test/decoupled/2p/grids/test_transport.dgf b/test/decoupled/2p/grids/test_transport.dgf
new file mode 100644
index 0000000000..0321c37549
--- /dev/null
+++ b/test/decoupled/2p/grids/test_transport.dgf
@@ -0,0 +1,11 @@
+DGF
+Interval
+0 0 % first corner
+1 1 % second corner
+10 10 % cells in x and y direction
+#
+
+BOUNDARYDOMAIN
+default 1 % all boundaries have id 1
+#BOUNDARYDOMAIN
+# unitcube.dgf
diff --git a/test/decoupled/2p/test_transport.cc b/test/decoupled/2p/test_transport.cc
index 71482ac270..f51a3df403 100644
--- a/test/decoupled/2p/test_transport.cc
+++ b/test/decoupled/2p/test_transport.cc
@@ -33,7 +33,7 @@
////////////////////////
void usage(const char *progname)
{
- std::cout << boost::format("usage: %s [--restart restartTime] tEnd\n")%progname;
+ std::cout << boost::format("usage: %s [--restart restartTime] gridFile.dgf tEnd\n")%progname;
exit(1);
}
@@ -54,7 +54,7 @@ int main(int argc, char** argv)
////////////////////////////////////////////////////////////
// parse the command line arguments
////////////////////////////////////////////////////////////
- if (argc < 2)
+ if (argc < 3)
usage(argv[0]);
// deal with the restart stuff
@@ -68,28 +68,27 @@ int main(int argc, char** argv)
std::istringstream(argv[argPos++]) >> restartTime;
}
- if (argc - argPos != 1) {
+ if (argc - argPos != 2) {
usage(argv[0]);
}
+ ////////////////////////////////////////////////////////////
+ // create the grid
+ ////////////////////////////////////////////////////////////
+ const char *dgfFileName = argv[argPos++];
+ Dune::GridPtr<Grid> gridPtr(dgfFileName);
+
// read the initial time step and the end time
double tEnd, dt;
std::istringstream(argv[argPos++]) >> tEnd;
dt = tEnd;
////////////////////////////////////////////////////////////
- // create the grid
+ // instantiate and run the concrete problem
////////////////////////////////////////////////////////////
- Dune::FieldVector<int,dim> N(4); N[0] = 4;
Dune::FieldVector<double ,dim> L(0);
Dune::FieldVector<double,dim> H(1);
- Grid grid(N,L,H);
-
- ////////////////////////////////////////////////////////////
- // instantiate and run the concrete problem
- ////////////////////////////////////////////////////////////
-
- Problem problem(grid.leafView(), L, H);
+ Problem problem(gridPtr->leafView(), L, H);
// load restart file if necessarry
if (restart)
diff --git a/test/decoupled/2p/test_transport_problem.hh b/test/decoupled/2p/test_transport_problem.hh
index 6a3fe9cd04..e055f93fe1 100644
--- a/test/decoupled/2p/test_transport_problem.hh
+++ b/test/decoupled/2p/test_transport_problem.hh
@@ -22,8 +22,7 @@
#include <dune/grid/uggrid.hh>
#endif
-#include <dune/grid/yaspgrid.hh>
-#include <dune/grid/sgrid.hh>
+#include <dune/grid/io/file/dgfparser/dgfug.hh>
#include <dumux/material/fluidsystems/liquidphase.hh>
#include <dumux/material/components/unit.hh>
@@ -51,8 +50,7 @@ NEW_TYPE_TAG(TransportTestProblem, INHERITS_FROM(DecoupledModel, Transport));
// Set the grid type
SET_PROP(TransportTestProblem, Grid)
{
- // typedef Dune::YaspGrid<2> type;
- typedef Dune::SGrid<2, 2> type;
+ typedef Dune::UGGrid<2> type;
};
// Set the problem property
--
GitLab