diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
index c5d04d53eeaaf1c3d5035995dcc0636194572903..77c046ede075d4aa627c92f0008c84b19138be63 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh
@@ -99,11 +99,6 @@ template<class TypeTag> class FVPressure2P: public FVPressure<TypeTag>
wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
};
- enum
- {
- rhs = ParentType::rhs, matrix = ParentType::matrix,
- };
-
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
@@ -113,6 +108,12 @@ template<class TypeTag> class FVPressure2P: public FVPressure<TypeTag>
typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix;
public:
+
+ enum
+ {
+ rhs = ParentType::rhs, matrix = ParentType::matrix,
+ };
+
void getSource(Dune::FieldVector<Scalar, 2>&, const Element&, const CellData&, const bool);
void getStorage(Dune::FieldVector<Scalar, 2>&, const Element&, const CellData&, const bool);
@@ -201,8 +202,6 @@ public:
}
}
- updateMaterialLaws();
-
ParentType::update();
storePressureSolution();
@@ -217,6 +216,7 @@ public:
{
CellData& cellData = problem_.variables().cellData(i);
storePressureSolution(i, cellData);
+ cellData.fluxData().resetVelocity();
}
}
diff --git a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
index 177cd765206c6ce6ea6a0bf24cd83202d7168eeb..36921788bac1075dfb5b435eb691e9b608128cd6 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh
@@ -26,7 +26,7 @@
// dumux environment
#include <dumux/decoupled/2p/2pproperties.hh>
-
+#include <dumux/decoupled/2p/diffusion/fv/fvpressure2p.hh>
/**
* @file
@@ -55,12 +55,13 @@ namespace Dumux
*
* \tparam TypeTag The Type Tag
*/
-template<class TypeTag> class FVPressure2Padaptive
+template<class TypeTag> class FVPressure2Padaptive: public FVPressure2P<TypeTag>
{
+ typedef FVPressure2P<TypeTag> ParentType;
+
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Variables) Variables;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters;
typedef typename SpatialParameters::MaterialLaw MaterialLaw;
@@ -73,6 +74,7 @@ template<class TypeTag> class FVPressure2Padaptive
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
enum
{
@@ -92,164 +94,37 @@ template<class TypeTag> class FVPressure2Padaptive
};
enum
{
- wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
+ };
+
+ enum
+ {
+ rhs = ParentType::rhs, matrix = ParentType::matrix,
};
typedef typename GridView::Traits::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<0>::Iterator ElementIterator;
- typedef typename GridView::Grid Grid;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::Intersection Intersection;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix;
- typedef typename GET_PROP_TYPE(TypeTag, PressureCoefficientMatrix) Matrix;
- typedef typename GET_PROP_TYPE(TypeTag, PressureRHSVector) Vector;
-
- //initializes the matrix to store the system of equations
- void initializeMatrix();
-
- //function which assembles the system of equations to be solved
- void assemble(bool first);
-
- //solves the system of equations to get the spatial distribution of the pressure
- void solve();
-
-protected:
- Problem& problem()
- {
- return problem_;
- }
-
- const Problem& problem() const
- {
- return problem_;
- }
-
public:
- //! updates and stores constitutive relations
- void updateMaterialLaws();
- //! \copydoc Dumux::FVPressure1P::initialize()
- void initialize(bool solveTwice = true)
- {
- updateMaterialLaws();
+ void getFlux(Dune::FieldVector<Scalar, 2>&, const Intersection&, const CellData&, const bool);
- assemble(true);
- solve();
- if (solveTwice)
- {
- Dune::BlockVector<Dune::FieldVector<Scalar, 1> > pressureOld(this->problem().variables().pressure());
-
- assemble(false);
- solve();
-
- Dune::BlockVector<Dune::FieldVector<Scalar, 1> > pressureDiff(pressureOld);
- pressureDiff -= this->problem().variables().pressure();
- pressureOld = this->problem().variables().pressure();
- Scalar pressureNorm = pressureDiff.infinity_norm();
- pressureNorm /= pressureOld.infinity_norm();
- int numIter = 0;
- while (pressureNorm > 1e-5 && numIter < 10)
- {
- updateMaterialLaws();
- assemble(false);
- solve();
-
- pressureDiff = pressureOld;
- pressureDiff -= this->problem().variables().pressure();
- pressureNorm = pressureDiff.infinity_norm();
- pressureOld = this->problem().variables().pressure();
- pressureNorm /= pressureOld.infinity_norm();
-
- numIter++;
- }
- // std::cout<<"Pressure defect = "<<pressureNorm<<"; "<<numIter<<" Iterations needed for initial pressure field"<<std::endl;
- }
- return;
- }
- //! \copydoc Dumux::FVPressure1P::pressure()
- void pressure(bool solveTwice = true)
- {
- assemble(false);
- solve();
-
- return;
- }
- //! updates the pressure field (analog to update function in Dumux::IMPET)
+ //pressure solution routine: update estimate for secants, assemble, solve.
void update()
{
- updateMaterialLaws();
-
- pressure(false);
-
- return;
- }
-
- // serialization methods
- //! \copydoc Dumux::FVPressure1P::serialize(Restarter &res)
- template<class Restarter>
- void serialize(Restarter &res)
- {
- return;
- }
-
- //! \copydoc Dumux::FVPressure1P::deserialize(Restarter &res)
- template<class Restarter>
- void deserialize(Restarter &res)
- {
- return;
- }
-
- //! \copydoc Dumux::FVPressure1P::addOutputVtkFields(MultiWriter &writer)
- template<class MultiWriter>
- void addOutputVtkFields(MultiWriter &writer)
- {
- typename Variables::ScalarSolutionType *pressure = writer.allocateManagedBuffer (problem_.gridView().size(0));
-
- *pressure = problem_.variables().pressure();
-
- if (pressureType == pw)
- {
- writer.attachCellData(*pressure, "wetting pressure");
- }
-
- if (pressureType == pn)
- {
- writer.attachCellData(*pressure, "nonwetting pressure");
- }
-
- if (pressureType == pglobal)
- {
- writer.attachCellData(*pressure, "global pressure");
- }
+ int gridSize = problem_.gridView().size(0);
+ // update RHS vector, matrix
+ this->A_.setSize(gridSize, gridSize); //
+ this->f_.resize(gridSize);
+ this->pressure().resize(gridSize);
- // output phase-dependent stuff
- typename Variables::ScalarSolutionType *pC = writer.allocateManagedBuffer (problem_.gridView().size(0));
- *pC = problem_.variables().capillaryPressure();
- writer.attachCellData(*pC, "capillary pressure");
+ ParentType::initializeMatrix();
- typename Variables::ScalarSolutionType *densityWetting = writer.allocateManagedBuffer (problem_.gridView().size(0));
- *densityWetting = problem_.variables().densityWetting();
- writer.attachCellData(*densityWetting, "wetting density");
-
- typename Variables::ScalarSolutionType *densityNonwetting = writer.allocateManagedBuffer (problem_.gridView().size(0));
- *densityNonwetting = problem_.variables().densityNonwetting();
- writer.attachCellData(*densityNonwetting, "nonwetting density");
-
- typename Variables::ScalarSolutionType *viscosityWetting = writer.allocateManagedBuffer (problem_.gridView().size(0));
- *viscosityWetting = problem_.variables().viscosityWetting();
- writer.attachCellData(*viscosityWetting, "wetting viscosity");
-
- typename Variables::ScalarSolutionType *viscosityNonwetting = writer.allocateManagedBuffer (problem_.gridView().size(0));
- *viscosityNonwetting = problem_.variables().viscosityNonwetting();
- writer.attachCellData(*viscosityNonwetting, "nonwetting viscosity");
-
-// typename Variables::ScalarSolutionType *saturation = writer.allocateManagedBuffer (problem_.gridView().size(0));
-//
-// *saturation = problem_.variables().saturation();
-//
-// writer.attachCellData(*saturation, "wetting saturation");
+ ParentType::update();
return;
}
@@ -259,966 +134,318 @@ public:
* \param problem a problem class object
*/
FVPressure2Padaptive(Problem& problem) :
- problem_(problem), A_(problem.variables().gridSize(), problem.variables().gridSize(), (4 * dim + 1)
- * problem.variables().gridSize(), Matrix::random), f_(problem.variables().gridSize()),
- gravity(problem.gravity())
-// problem_(problem), A_(problem.variables().gridSize(), problem.variables().gridSize(), (4 * dim + 1)
-// * problem.variables().gridSize(), Matrix::random), f_(problem.variables().gridSize()),
-// gravity(problem.gravity())
- {// Comment: the term (4 * dim + 1) only works in the 2D case (or less). For 3D applications, it has to be (8*dim+1)
- if (pressureType != pw && pressureType != pn && pressureType != pglobal)
+ ParentType(problem), problem_(problem), gravity_(problem.gravity())
+ {
+ if (pressureType_ != pw && pressureType_ != pn && pressureType_ != pglobal)
{
DUNE_THROW(Dune::NotImplemented, "Pressure type not supported!");
}
- if (saturationType != Sw && saturationType != Sn)
+ if (pressureType_ == pglobal && compressibility_)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Compressibility not supported for global pressure!");
+ }
+ if (saturationType_ != Sw && saturationType_ != Sn)
{
DUNE_THROW(Dune::NotImplemented, "Saturation type not supported!");
}
- initializeMatrix();
+ if (!compressibility_)
+ {
+ const Element& element = *(problem_.gridView().template begin<0>());
+ FluidState fluidState;
+ fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
+ fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
+ fluidState.setTemperature(problem_.temperature(element));
+ fluidState.setSaturation(wPhaseIdx, 1.);
+ fluidState.setSaturation(nPhaseIdx, 0.);
+ density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
+ density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
+ viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
+ }
}
private:
Problem& problem_;
- Matrix A_;
- Dune::BlockVector<Dune::FieldVector<Scalar, 1> > f_;
-protected:
- const GlobalPosition& gravity; //!< vector including the gravity constant
- static const bool compressibility = GET_PROP_VALUE(TypeTag, EnableCompressibility);
- static const int pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$p_w\f$, \f$p_n\f$, \f$p_{global}\f$)
- static const int saturationType = GET_PROP_VALUE(TypeTag, SaturationFormulation); //!< gives kind of saturation used (\f$S_w\f$, \f$S_n\f$)
+ const GlobalPosition& gravity_; //!< vector including the gravity constant
+
+ Scalar density_[numPhases];
+ Scalar viscosity_[numPhases];
+
+ static const bool compressibility_ = GET_PROP_VALUE(TypeTag, EnableCompressibility);
+ static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$p_w\f$, \f$p_n\f$, \f$p_{global}\f$)
+ static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation); //!< gives kind of saturation used (\f$S_w\f$, \f$S_n\f$)
};
-//initializes the matrix to store the system of equations
+//!function which assembles the system of equations to be solved. Only works for 2D.
template<class TypeTag>
-void FVPressure2Padaptive<TypeTag>::initializeMatrix()
+void FVPressure2Padaptive<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries, const Intersection& intersection
+ , const CellData& cellDataI, const bool first)
{
- // determine matrix row sizes
- ElementIterator eItEnd = problem_.gridView().template end<0> ();
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ ElementPointer elementI = intersection.inside();
+ ElementPointer elementJ = intersection.outside();
+
+ if (elementI->level() == elementJ->level())
{
- // cell index
- int globalIdxI = problem_.variables().index(*eIt);
-
- // initialize row size
- int rowSize = 1;
-
- // run through all intersections with neighbors
- IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
- for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
- if (isIt->neighbor())
- rowSize++;
- A_.setrowsize(globalIdxI, rowSize);
+ ParentType::getFlux(entries, intersection, cellDataI, first);
}
- A_.endrowsizes();
-
- // determine position of matrix entries
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
+ // hanging node situation: neighbor has higher level
+ else if (elementJ->level() == elementI->level() + 1)
{
- // cell index
- int globalIdxI = problem_.variables().index(*eIt);
+ const CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(*elementJ));
+
+ // get global coordinates of cell centers
+ const GlobalPosition& globalPosI = elementI->geometry().center();
+ const GlobalPosition& globalPosJ = elementJ->geometry().center();
- // add diagonal index
- A_.addindex(globalIdxI, globalIdxI);
+ int globalIdxI = problem_.variables().index(*elementI);
+ int globalIdxJ = problem_.variables().index(*elementJ);
- // run through all intersections with neighbors
- IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
- for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
- if (isIt->neighbor())
+ // get mobilities and fractional flow factors
+ Scalar lambdaWI = cellDataI.mobility(wPhaseIdx);
+ Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx);
+ Scalar fractionalWI = cellDataI.fracFlowFunc(wPhaseIdx);
+ Scalar fractionalNWI = cellDataI.fracFlowFunc(nPhaseIdx);
+ Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx);
+ Scalar lambdaNWJ = cellDataJ.mobility(nPhaseIdx);
+ Scalar fractionalWJ = cellDataJ.fracFlowFunc(wPhaseIdx);
+ Scalar fractionalNWJ = cellDataJ.fracFlowFunc(nPhaseIdx);
+
+ // get capillary pressure
+ Scalar pcI = cellDataI.capillaryPressure();
+ Scalar pcJ = cellDataJ.capillaryPressure();
+
+ //get face index
+ int isIndexI = intersection.indexInInside();
+
+ //get face normal
+ const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
+
+ // get face area
+ Scalar faceArea = intersection.geometry().volume();
+
+ // distance vector between barycenters
+ GlobalPosition distVecIJ = globalPosJ - globalPosI;
+
+ // compute distance between cell centers
+ Scalar distIJ = distVecIJ.two_norm();
+
+ // Count number of hanging nodes
+ // not really necessary
+ int isIndexJ = intersection.indexInOutside();
+ int globalIdxK = 0;
+ ElementPointer elementK = intersection.outside();
+ // We are looking for two things:
+ // IsIndexJ, the index of the interface from the neighbor-cell point of view
+ // GlobalIdxK, the index of the third cell
+ // Intersectioniterator around cell I
+ IntersectionIterator isItEndI = problem_.gridView().iend(*elementI);
+ for (IntersectionIterator isItI = problem_.gridView().ibegin(*elementI); isItI != isItEndI; ++isItI)
+ {
+ if (isItI->neighbor())
{
- // access neighbor
- ElementPointer outside = isIt->outside();
- int globalIdxJ = problem_.variables().index(*outside);
+ ElementPointer neighborPointer2 = isItI->outside();
- // add off diagonal index
- A_.addindex(globalIdxI, globalIdxJ);
+ // make sure we do not chose elemntI as third element
+ // -> faces with hanging node have more than one intersection but only one face index!
+ if (neighborPointer2 != elementJ && isItI->indexInInside() == isIndexI)
+ {
+ globalIdxK = problem_.variables().index(*neighborPointer2);
+ elementK = neighborPointer2;
+ break;
+ }
}
- }
- A_.endindices();
+ }
- return;
-}
+ const CellData& cellDataK = problem_.variables().cellData(globalIdxK);
-//!function which assembles the system of equations to be solved. Only works for 2D.
-template<class TypeTag>
-void FVPressure2Padaptive<TypeTag>::assemble(bool first)
-{
- updateMaterialLaws();
- // Adapt size of matrix A_ and Vector f_
- {int size = problem_.variables().gridSize();
- A_.setSize(size,size);
- f_.resize(size);}
- initializeMatrix(); // define non-zero entries in matrix
+ // neighbor cell center in global coordinates
+ const GlobalPosition& globalPosInterface = intersection.geometry().center();
- // initialization: set matrix A_ to zero
- A_ = 0;
- f_ = 0;
+ Dune::FieldVector < Scalar, dimWorld > distVec = globalPosInterface - globalPosI;
+ Scalar lI = distVec * unitOuterNormal;
+ distVec = globalPosJ - globalPosInterface;
+ Scalar lJ = distVec * unitOuterNormal;
+ Scalar l = lI + lJ;
- int nodecount = 0;
+ FieldMatrix permeabilityI(0);
+ FieldMatrix permeabilityJ(0);
+ FieldMatrix permeabilityK(0);
- BoundaryTypes bcType;
+ problem_.spatialParameters().meanK(permeabilityI,
+ problem_.spatialParameters().intrinsicPermeability(*elementI));
+ problem_.spatialParameters().meanK(permeabilityJ,
+ problem_.spatialParameters().intrinsicPermeability(*elementJ));
+ problem_.spatialParameters().meanK(permeabilityK,
+ problem_.spatialParameters().intrinsicPermeability(*elementK));
- ElementIterator eItEnd = problem_.gridView().template end<0> ();
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
- {
- // cell index
- int globalIdxI = problem_.variables().index(*eIt);
+ // Calculate permeablity component normal to interface
+ Scalar kI, kJ, kK, kMean, ng;
+ Dune::FieldVector < Scalar, dim > permI(0);
+ Dune::FieldVector < Scalar, dim > permJ(0);
+ Dune::FieldVector < Scalar, dim > permK(0);
- // get global coordinate of cell center
- const GlobalPosition& globalPos = eIt->geometry().center();
+ permeabilityI.mv(unitOuterNormal, permI);
+ permeabilityJ.mv(unitOuterNormal, permJ);
+ permeabilityK.mv(unitOuterNormal, permK);
+ // kI,kJ,kK=(n^T)Kn
+ kI = unitOuterNormal * permI;
+ kJ = unitOuterNormal * permJ;
+ kK = unitOuterNormal * permK;
- // cell volume, assume linear map here
- Scalar volume = eIt->geometry().volume();
+ // See Diplomarbeit Michael Sinsbeck
+ kMean = kI * kJ * kK * l / (kJ * kK * lI + kI * (kJ + kK) / 2 * lJ);
- Scalar densityWI = problem_.variables().densityWetting(globalIdxI);
- Scalar densityNWI = problem_.variables().densityNonwetting(globalIdxI);
+ ng = gravity_ * unitOuterNormal;
- // set sources
- PrimaryVariables source(0.0);
- problem().source(source, *eIt);
- if (!compressibility)
- {
- source[wPhaseIdx] /= densityWI;
- source[nPhaseIdx] /= densityNWI;
- }
- f_[globalIdxI] += volume * (source[wPhaseIdx] + source[nPhaseIdx]);
+ Scalar lambdaWK = cellDataK.mobility(wPhaseIdx);
+ Scalar lambdaNWK = cellDataK.mobility(nPhaseIdx);
+ Scalar fractionalWK = cellDataK.fracFlowFunc(wPhaseIdx);
+ Scalar fractionalNWK = cellDataK.fracFlowFunc(nPhaseIdx);
- Scalar porosity = problem_.spatialParameters().porosity(*eIt);
+ Scalar pcK = cellDataK.capillaryPressure();
+ Scalar pcJK = (pcJ + pcK) / 2;
- // get mobilities and fractional flow factors
- Scalar lambdaWI = problem_.variables().mobilityWetting(globalIdxI);
- Scalar lambdaNWI = problem_.variables().mobilityNonwetting(globalIdxI);
- Scalar fractionalWI = problem_.variables().fracFlowFuncWetting(globalIdxI);
- Scalar fractionalNWI = problem_.variables().fracFlowFuncNonwetting(globalIdxI);
- Scalar pcI = problem_.variables().capillaryPressure(globalIdxI);
-
- IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
- for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
- {
- int isIndex = isIt->indexInInside();
+ // Potentials from previous time step are not available.
+ // Instead, calculate mean density, then find potentials,
+ // then upwind density.
+ // pressure from previous time step might also be incorrect.
- // get normal vector
- const GlobalPosition& unitOuterNormal = isIt->centerUnitOuterNormal();
+ Scalar rhoMeanWIJ = density_[wPhaseIdx];
+ Scalar rhoMeanNWIJ = density_[nPhaseIdx];
+ Scalar rhoMeanWIK = density_[wPhaseIdx];
+ Scalar rhoMeanNWIK = density_[nPhaseIdx];
- // get face volume
- Scalar faceArea = isIt->geometry().volume();
+ if (compressibility_)
+ {
+ rhoMeanWIJ = (lJ * cellDataI.density(wPhaseIdx) + lI * cellDataJ.density(wPhaseIdx)) / l;
+ rhoMeanNWIJ = (lJ * cellDataI.density(nPhaseIdx) + lI * cellDataJ.density(nPhaseIdx)) / l;
+ rhoMeanWIK = (lJ * cellDataI.density(wPhaseIdx) + lI * cellDataK.density(wPhaseIdx)) / l;
+ rhoMeanNWIK = (lJ * cellDataI.density(nPhaseIdx) + lI * cellDataK.density(nPhaseIdx)) / l;
+ }
- // handle interior face
- if (isIt->neighbor())
- {
- // access neighbor
- ElementPointer neighborPointer = isIt->outside();
- int globalIdxJ = problem_.variables().index(*neighborPointer);
+ Scalar densityWIJ = density_[wPhaseIdx];
+ Scalar densityNWIJ = density_[nPhaseIdx];
+ Scalar densityWIK = density_[wPhaseIdx];
+ Scalar densityNWIK = density_[nPhaseIdx];
- // "normal" case: neighbor has same level
- if (neighborPointer->level()==eIt.level())
- {
- // neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
-
- // distance vector between barycenters
- GlobalPosition distVec = globalPosNeighbor - globalPos;
-
- // compute distance between cell centers
-// Scalar dist = distVec.two_norm();
- Scalar dist = distVec*unitOuterNormal;
-
- // compute vectorized permeabilities
- FieldMatrix meanPermeability(0);
-
- problem_.spatialParameters().meanK(meanPermeability,
- problem_.spatialParameters().intrinsicPermeability(*eIt),
- problem_.spatialParameters().intrinsicPermeability(*neighborPointer));
-
- Dune::FieldVector<Scalar, dim> permeability(0);
- meanPermeability.mv(unitOuterNormal, permeability);
-
- // get mobilities and fractional flow factors
- Scalar lambdaWJ = problem_.variables().mobilityWetting(globalIdxJ);
- Scalar lambdaNWJ = problem_.variables().mobilityNonwetting(globalIdxJ);
- Scalar fractionalWJ = problem_.variables().fracFlowFuncWetting(globalIdxJ);
- Scalar fractionalNWJ = problem_.variables().fracFlowFuncNonwetting(globalIdxJ);
- Scalar densityWJ = problem_.variables().densityWetting(globalIdxJ);
- Scalar densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
-
- Scalar pcJ = problem_.variables().capillaryPressure(globalIdxJ);
-
- Scalar rhoMeanW = 0.5 * (densityWI + densityWJ);
- Scalar rhoMeanNW = 0.5 * (densityNWI + densityNWJ);
- Scalar fMeanW = 0.5 * (fractionalWI + fractionalWJ);
- Scalar fMeanNW = 0.5 * (fractionalNWI + fractionalNWJ);
-
- // update diagonal entry
- Scalar entry;
-
- //calculate potential gradients
- Scalar potentialW = 0;
- Scalar potentialNW = 0;
-
- Scalar densityW = 0;
- Scalar densityNW = 0;
-
- //if we are at the very first iteration we can't calculate phase potentials
- if (!first)
- {
- // After grid adaption, the potentials from the previous time step
- // can not be used for upwinding. The mean values is assumed instead.
- densityW = rhoMeanW;
- densityNW = rhoMeanNW;
-
- switch (pressureType)
- {
- case pw:
- {
- potentialW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ]);
- potentialNW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] + pcI
- - pcJ);
- break;
- }
- case pn:
- {
- potentialW
- = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] - pcI + pcJ);
- potentialNW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ]);
- break;
- }
- case pglobal:
- {
- potentialW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] - fMeanNW
- * (pcI - pcJ));
- potentialNW = (problem_.variables().pressure()[globalIdxI] - problem_.variables().pressure()[globalIdxJ] + fMeanW
- * (pcI - pcJ));
- break;
- }
- }
-
- potentialW += densityW * (distVec * gravity);
- potentialNW += densityNW * (distVec * gravity);
-
- //store potentials for further calculations (velocity, saturation, ...)
- problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
- problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
- }
-
- //do the upwinding of the mobility depending on the phase potentials
- Scalar lambdaW = (potentialW > 0.) ? lambdaWI : lambdaWJ;
- lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
- Scalar lambdaNW = (potentialNW > 0) ? lambdaNWI : lambdaNWJ;
- lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNW;
-
- densityW = (potentialW > 0.) ? densityWI : densityWJ;
- densityNW = (potentialNW > 0.) ? densityNWI : densityNWJ;
-
- densityW = (potentialW == 0) ? rhoMeanW : densityW;
- densityNW = (potentialNW == 0) ? rhoMeanNW : densityNW;
-
- //calculate current matrix entry
- entry = (lambdaW + lambdaNW) * ((permeability * unitOuterNormal) / dist) * faceArea;
-
- //calculate right hand side
- Scalar rightEntry = (lambdaW * densityW + lambdaNW * densityNW) * (permeability * gravity) * faceArea;
-
- switch (pressureType)
- {
- case pw:
- {
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
- pCGradient *= (pcI - pcJ) / dist;
-
- //add capillary pressure term to right hand side
- rightEntry += 0.5 * (lambdaNWI + lambdaNWJ) * (permeability * pCGradient) * faceArea;
- break;
- }
- case pn:
- {
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
- pCGradient *= (pcI - pcJ) / dist;
-
- //add capillary pressure term to right hand side
- rightEntry -= 0.5 * (lambdaWI + lambdaWJ) * (permeability * pCGradient) * faceArea;
- break;
- }
- }
-
- //set right hand side
- f_[globalIdxI] -= rightEntry;
-
- // set diagonal entry
- A_[globalIdxI][globalIdxI] += entry;
-
- // set off-diagonal entry
- A_[globalIdxI][globalIdxJ] -= entry;
-
- } // end of case "neighbor has same level"
-
- // hanging node situation: neighbor has higher level
- if (neighborPointer->level()==eIt.level()+1)
- {
- // Count number of hanging nodes
- // nodecount counts each hanging node twice, we have to divide by 2 in the end
- // not really necessary
- nodecount++;
-
- int globalIdxK = 0;
- ElementPointer thirdCellPointer = isIt->outside();
- bool foundK=false;
- bool foundIJ=false;
- // We are looking for two things:
- // IsIndexJ, the index of the interface from the neighbor-cell point of view
- // GlobalIdxK, the index of the third cell
- // for efficienty this is done in one IntersectionIterator-Loop
-
- // Intersectioniterator around cell J
- IntersectionIterator isItEndJ = this->problem().gridView().iend(*neighborPointer);
-
- for (IntersectionIterator isItJ = this->problem().gridView().ibegin(*neighborPointer); isItJ != isItEndJ; ++isItJ)
- {
- if (isItJ->neighbor())
- {
- ElementPointer neighborPointer2 = isItJ->outside();
-
- // Neighbor of neighbor is Cell I?
- if (this->problem().variables().index(*neighborPointer2)==globalIdxI)
- {
- foundIJ=true;
- }
- else
- {
- if (neighborPointer2->level()==eIt.level()+1)
- {
- // To verify, we found the correct Cell K, we check
- // - is level(K)=level(J)?
- // - is distance(IJ)=distance(IK)?
- // - K is neighbor of J.
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
- const GlobalPosition& globalPosThirdCell = neighborPointer2->geometry().center();
- Dune::FieldVector<Scalar, dimWorld> distVecIJ = globalPosNeighbor - globalPos;
- Dune::FieldVector<Scalar, dimWorld> distVecIK = globalPosThirdCell - globalPos;
- if (((distVecIK.two_norm()-distVecIJ.two_norm()))/distVecIJ.two_norm() < 0.001)
- {
- globalIdxK= this->problem().variables().index(*neighborPointer2);
- thirdCellPointer = neighborPointer2;
- foundK=true;
- }
-
- }
- }
- }
- if (foundIJ && foundK) break;
- }
-
- int isIndexJ = isIt->indexInOutside();
-
- // neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
-// const GlobalPosition& globalPosThirdCell = thirdCellPointer->geometry().center();
- const GlobalPosition& globalPosInterface = isIt->geometry().center();
-
- Dune::FieldVector<Scalar,dimWorld> distVec = globalPosInterface - globalPos;
- Scalar lI= distVec*unitOuterNormal;
- distVec = globalPosNeighbor - globalPosInterface;
- Scalar lJ= distVec*unitOuterNormal;
- Scalar l=lI+lJ;
-
- FieldMatrix permeabilityI(0);
- FieldMatrix permeabilityJ(0);
- FieldMatrix permeabilityK(0);
-
- problem_.spatialParameters().meanK(permeabilityI,
- problem_.spatialParameters().intrinsicPermeability(*eIt));
- problem_.spatialParameters().meanK(permeabilityJ,
- problem_.spatialParameters().intrinsicPermeability(*neighborPointer));
- problem_.spatialParameters().meanK(permeabilityK,
- problem_.spatialParameters().intrinsicPermeability(*thirdCellPointer));
-
- // Calculate permeablity component normal to interface
- Scalar kI, kJ, kK, kMean, ng;
- Dune::FieldVector<Scalar, dim> permI(0);
- Dune::FieldVector<Scalar, dim> permJ(0);
- Dune::FieldVector<Scalar, dim> permK(0);
-
- permeabilityI.mv(unitOuterNormal, permI);
- permeabilityJ.mv(unitOuterNormal, permJ);
- permeabilityK.mv(unitOuterNormal, permK);
-
- // kI,kJ,kK=(n^T)Kn
- kI=unitOuterNormal*permI;
- kJ=unitOuterNormal*permJ;
- kK=unitOuterNormal*permK;
-
- // See Diplomarbeit Michael Sinsbeck
- kMean=kI*kJ*kK*l/(kJ*kK*lI+kI*(kJ+kK)/2*lJ);
-
- ng=this->gravity*unitOuterNormal;
-
- // get mobilities
- // fractional flow function is not evaluated, since we do not use p_global
- Scalar lambdaWJ = problem_.variables().mobilityWetting(globalIdxJ);
- Scalar lambdaNWJ = problem_.variables().mobilityNonwetting(globalIdxJ);
- Scalar densityWJ = problem_.variables().densityWetting(globalIdxJ);
- Scalar densityNWJ = problem_.variables().densityNonwetting(globalIdxJ);
- Scalar fractionalWJ = problem_.variables().fracFlowFuncWetting(globalIdxJ);
- Scalar fractionalNWJ = problem_.variables().fracFlowFuncNonwetting(globalIdxJ);
-
-// Scalar lambdaWK = problem_.variables().mobilityWetting(globalIdxK);
-// Scalar lambdaNWK = problem_.variables().mobilityNonwetting(globalIdxK);
- Scalar densityWK = problem_.variables().densityWetting(globalIdxK);
- Scalar densityNWK = problem_.variables().densityNonwetting(globalIdxK);
- Scalar fractionalWK = problem_.variables().fracFlowFuncWetting(globalIdxK);
- Scalar fractionalNWK = problem_.variables().fracFlowFuncNonwetting(globalIdxK);
-
- Scalar pcJ = problem_.variables().capillaryPressure(globalIdxJ);
- Scalar pcK = problem_.variables().capillaryPressure(globalIdxK);
- Scalar pcJK=(pcJ+pcK)/2;
-
- // Potentials from previous time step are not available.
- // Instead, calculate mean density, then find potentials,
- // then upwind density.
- // pressure from previous time step might also be incorrect.
-
- Scalar rhoMeanWIJ = (lJ*densityWI + lI*densityWJ)/l;
- Scalar rhoMeanNWIJ = (lJ*densityNWI + lI*densityNWJ)/l;
- Scalar rhoMeanWIK = (lJ*densityWI + lI*densityWK)/l;
- Scalar rhoMeanNWIK = (lJ*densityNWI + lI*densityNWK)/l;
-
- Scalar densityWIJ = 0;
- Scalar densityNWIJ = 0;
- Scalar densityWIK = 0;
- Scalar densityNWIK = 0;
-
- Scalar fMeanWIJ = (lJ*fractionalWI+lI*fractionalWJ)/l;
- Scalar fMeanNWIJ = (lJ*fractionalNWI+lI*fractionalNWJ)/l;
- Scalar fMeanWIK = (lJ*fractionalWI+lI*fractionalWK)/l;
- Scalar fMeanNWIK = (lJ*fractionalNWI+lI*fractionalNWK)/l;
-
- Scalar potentialWIJ = 0;
- Scalar potentialNWIJ = 0;
- Scalar potentialWIK = 0;
- Scalar potentialNWIK = 0;
-
- //if we are at the very first iteration we can't calculate phase potentials
- if (!first)
- {
- // potentials from previous time step no available.
- densityWIJ = rhoMeanWIJ;
- densityNWIJ = rhoMeanNWIJ;
- densityWIK = rhoMeanWIK;
- densityNWIK = rhoMeanNWIK;
-
- // Old pressure field
- Scalar pressI=problem_.variables().pressure()[globalIdxI];
- Scalar pressJ=problem_.variables().pressure()[globalIdxJ];
- Scalar pressK=problem_.variables().pressure()[globalIdxK];
- Scalar pressJK=(pressJ+pressK)/2;
- Scalar pcJK=(pcJ+pcK)/2;
-
- switch (pressureType)
- {
- case pw:
- {
- potentialWIJ = (pressI-pressJK)/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI+pcI-(pressJK+pcJK))/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-pressJK)/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI+pcI-(pressJK+pcJK))/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- case pn:
- {
- potentialWIJ = (pressI-pcI-(pressJK-pcJK))/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI-pressJK)/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-pcI-(pressJK-pcJK))/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI-pressJK)/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- case pglobal:
- {
- potentialWIJ = (pressI-fMeanNWIJ*pcI-(pressJK-fMeanNWIJ*pcJK))/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI+fMeanWIJ*pcI-(pressJK+fMeanWIJ*pcJK))/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-fMeanNWIK*pcI-(pressJK-fMeanNWIK*pcJK))/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI+fMeanWIK*pcI-(pressJK+fMeanWIK*pcJK))/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- }
-
- //store potentials for further calculations (velocity, saturation, ...)
-
- problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialWIJ;
- problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNWIJ;
- problem_.variables().potentialWetting(globalIdxJ, isIndexJ) = -potentialWIJ;
- problem_.variables().potentialNonwetting(globalIdxJ, isIndexJ) = -potentialNWIJ;
- }
-
- //do the upwinding of the mobility depending on the phase potentials
- Scalar lambdaWIJ = (potentialWIJ > 0.) ? lambdaWI : lambdaWJ;
- lambdaWIJ = (potentialWIJ == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
- Scalar lambdaNWIJ = (potentialNWIJ > 0) ? lambdaNWI : lambdaNWJ;
- lambdaNWIJ = (potentialNWIJ == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNWIJ;
- densityWIJ = (potentialWIJ > 0.) ? densityWI : densityWJ;
- densityNWIJ = (potentialNWIJ > 0.) ? densityNWI : densityNWJ;
- densityWIJ = (potentialWIJ == 0) ? rhoMeanWIJ : densityWIJ;
- densityNWIJ = (potentialNWIJ == 0) ? rhoMeanNWIJ : densityNWIJ;
-
- densityWIK = (potentialWIK > 0.) ? densityWI : densityWK;
- densityNWIK = (potentialNWIK > 0.) ? densityNWI : densityNWK;
- densityWIK = (potentialWIK == 0) ? rhoMeanWIK : densityWIK;
- densityNWIK = (potentialNWIK == 0) ? rhoMeanNWIK : densityNWIK;
-
-
- // update diagonal entry and right hand side
-
- Scalar entry;
- Scalar rightEntry;
-
- entry = (lambdaWIJ + lambdaNWIJ)*kMean/l*faceArea;
- rightEntry = faceArea*lambdaWIJ*kMean*ng*((densityWIJ)-(lJ/l)*(kI+kK)/kI*(densityWIK-densityWIJ)/2);
- rightEntry += faceArea*lambdaNWIJ*kMean*ng*((densityNWIJ)-(lJ/l)*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2);
-
- switch (pressureType)
- {
- case pw:
- {
- rightEntry += faceArea*lambdaNWIJ*kMean*(pcJK-pcI)/l;
- break;
- }
- case pn:
- {
- rightEntry -= faceArea*lambdaWIJ*kMean*(pcJK-pcI)/l;
- break;
- }
- }
-
- f_[globalIdxI] -= rightEntry;
- f_[globalIdxJ] += rightEntry;
-
- // set diagonal entry
- A_[globalIdxI][globalIdxI] += entry;
- A_[globalIdxI][globalIdxJ] -= 0.5*entry;
- A_[globalIdxI][globalIdxK] -= 0.5*entry;
-
- // set off-diagonal entry
- A_[globalIdxJ][globalIdxI] -= entry;
- A_[globalIdxJ][globalIdxJ] += 0.5*entry;
- A_[globalIdxJ][globalIdxK] += 0.5*entry;
+ Scalar fMeanWIJ = (lJ * fractionalWI + lI * fractionalWJ) / l;
+ Scalar fMeanNWIJ = (lJ * fractionalNWI + lI * fractionalNWJ) / l;
+ Scalar fMeanWIK = (lJ * fractionalWI + lI * fractionalWK) / l;
+ Scalar fMeanNWIK = (lJ * fractionalNWI + lI * fractionalNWK) / l;
- }
- }
+ Scalar potentialWIJ = 0;
+ Scalar potentialNWIJ = 0;
+ Scalar potentialWIK = 0;
+ Scalar potentialNWIK = 0;
- // boundary face
-
- else
+ //if we are at the very first iteration we can't calculate phase potentials
+ if (!first)
+ {
+ // potentials from previous time step no available.
+ if (compressibility_)
{
- // center of face in global coordinates
- const GlobalPosition& globalPosFace = isIt->geometry().center();
-
- problem().boundaryTypes(bcType, *isIt);
- PrimaryVariables boundValues(0.0);
-
- GlobalPosition distVec(globalPosFace - globalPos);
-// Scalar dist = distVec.two_norm();
- Scalar dist = distVec * unitOuterNormal;
-
- if (bcType.isDirichlet(eqIdxPress))
- {
- problem().dirichlet(boundValues, *isIt);
-
- //permeability vector at boundary
- // compute vectorized permeabilities
- FieldMatrix meanPermeability(0);
-
- problem_.spatialParameters().meanK(meanPermeability,
- problem_.spatialParameters().intrinsicPermeability(*eIt));
-
- Dune::FieldVector<Scalar, dim> permeability(0);
- meanPermeability.mv(unitOuterNormal, permeability);
-
- //determine saturation at the boundary -> if no saturation is known directly at the boundary use the cell saturation
- Scalar satBound;
- if (bcType.isDirichlet(eqIdxSat))
- {
- satBound = boundValues[saturationIdx];
- }
- else
- {
- satBound = problem_.variables().saturation()[globalIdxI];
- }
- Scalar temperature = problem_.temperature(*eIt);
-
- //get dirichlet pressure boundary condition
- Scalar pressBound = boundValues[pressureIdx];
-
- //calculate consitutive relations depending on the kind of saturation used
- //determine phase saturations from primary saturation variable
- Scalar satW = 0;
- //Scalar satNW = 0;
- switch (saturationType)
- {
- case Sw:
- {
- satW = satBound;
- //satNW = 1 - satBound;
- break;
- }
- case Sn:
- {
- satW = 1 - satBound;
- //satNW = satBound;
- break;
- }
- }
-
- Scalar pcI = problem_.variables().capillaryPressure(globalIdxI);
- Scalar pcBound = MaterialLaw::pC(problem_.spatialParameters().materialLawParams(*eIt), satW);
-
- //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;
- }
- }
-
- Scalar densityWBound = 0;
- Scalar densityNWBound = 0;
- Scalar lambdaWBound = 0;
- Scalar lambdaNWBound = 0;
-
- if (compressibility)
- {
- FluidState fluidState;
- fluidState.update(satW, pressW, pressNW, temperature);
- densityWBound = FluidSystem::phaseDensity(wPhaseIdx, temperature, pressW, fluidState);
- densityNWBound = FluidSystem::phaseDensity(nPhaseIdx, temperature, pressNW, fluidState);
- Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx, temperature, pressW, fluidState);
- Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx, temperature, pressNW, fluidState);
- lambdaWBound = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), satW)
- / viscosityWBound * densityWBound;
- lambdaNWBound = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), satW)
- / viscosityNWBound * densityNWBound;
- }
- else
- {
- Scalar referencePressure = problem_.referencePressure(*eIt);
- FluidState fluidState;
- fluidState.update(satW, referencePressure, referencePressure, temperature);
-
- densityWBound = FluidSystem::phaseDensity(wPhaseIdx, temperature, referencePressure, fluidState);
- densityNWBound = FluidSystem::phaseDensity(nPhaseIdx, temperature, referencePressure, fluidState);
- Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx, temperature, referencePressure, fluidState);
- Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx, temperature, referencePressure, fluidState);
- lambdaWBound = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), satW)
- / viscosityWBound;
- lambdaNWBound = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), satW)
- / viscosityNWBound;
- }
- Scalar fractionalWBound = lambdaWBound / (lambdaWBound + lambdaNWBound);
- Scalar fractionalNWBound = lambdaNWBound / (lambdaWBound + lambdaNWBound);
-
- Scalar rhoMeanW = 0.5 * (densityWI + densityWBound);
- Scalar rhoMeanNW = 0.5 * (densityNWI + densityNWBound);
- Scalar fMeanW = 0.5 * (fractionalWI + fractionalWBound);
- Scalar fMeanNW = 0.5 * (fractionalNWI + fractionalNWBound);
-
- Scalar potentialW = 0;
- Scalar potentialNW = 0;
-
- Scalar densityW = 0;
- Scalar densityNW = 0;
-
- if (!first)
- {
- // Comment: In the non-adaptive case, one can upwind the density using the
- // old potential. Here, we use the mean density instead.
- // For the incompressible case, it does not make a difference, anyway.
- densityW = rhoMeanW;
- densityNW = rhoMeanNW;
-
- //calculate potential gradient
- switch (pressureType)
- {
- case pw:
- {
- potentialW = (problem_.variables().pressure()[globalIdxI] - pressBound);
- potentialNW = (problem_.variables().pressure()[globalIdxI] + pcI - pressBound - pcBound);
- break;
- }
- case pn:
- {
- potentialW = (problem_.variables().pressure()[globalIdxI] - pcI - pressBound + pcBound);
- potentialNW = (problem_.variables().pressure()[globalIdxI] - pressBound);
- break;
- }
- case pglobal:
- {
- potentialW = (problem_.variables().pressure()[globalIdxI] - pressBound - fMeanNW * (pcI - pcBound));
- potentialNW = (problem_.variables().pressure()[globalIdxI] - pressBound + fMeanW * (pcI - pcBound));
- break;
- }
- }
-
- potentialW += densityW * (distVec * gravity);
- potentialNW += densityNW * (distVec * gravity);
-
- //store potential gradients for further calculations
- problem_.variables().potentialWetting(globalIdxI, isIndex) = potentialW;
- problem_.variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
- }
-
- //do the upwinding of the mobility depending on the phase potentials
- Scalar lambdaW = (potentialW > 0.) ? lambdaWI : lambdaWBound;
- lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
- Scalar lambdaNW = (potentialNW > 0.) ? lambdaNWI : lambdaNWBound;
- lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWBound) : lambdaNW;
-
- densityW = (potentialW > 0.) ? densityWI : densityWBound;
- densityW = (potentialW == 0) ? rhoMeanW : densityW;
- densityNW = (potentialNW > 0.) ? densityNWI : densityNWBound;
- densityNW = (potentialNW == 0) ? rhoMeanNW : densityNW;
-
- //calculate current matrix entry
- Scalar entry = (lambdaW + lambdaNW) * ((permeability * unitOuterNormal) / dist) * faceArea;
-
- //calculate right hand side
- Scalar rightEntry = (lambdaW * densityW + lambdaNW * densityNW) * (permeability * gravity) * faceArea;
-
- switch (pressureType)
- {
- case pw:
- {
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
- pCGradient *= (pcI - pcBound) / dist;
-
- //add capillary pressure term to right hand side
- rightEntry += 0.5 * (lambdaNWI + lambdaNWBound) * (permeability * pCGradient) * faceArea;
- break;
- }
- case pn:
- {
- // calculate capillary pressure gradient
- Dune::FieldVector<Scalar, dim> pCGradient = unitOuterNormal;
- pCGradient *= (pcI - pcBound) / dist;
-
- //add capillary pressure term to right hand side
- rightEntry -= 0.5 * (lambdaWI + lambdaWBound) * (permeability * pCGradient) * faceArea;
- break;
- }
- }
-
- // set diagonal entry and right hand side entry
- A_[globalIdxI][globalIdxI] += entry;
- f_[globalIdxI] += entry * pressBound;
- f_[globalIdxI] -= rightEntry;
- }
- //set neumann boundary condition
-
- else if (bcType.isNeumann(eqIdxPress))
- {
- problem().neumann(boundValues, *isIt);
-
- if (!compressibility)
- {
- boundValues[wPhaseIdx] /= densityWI;
- boundValues[nPhaseIdx] /= densityNWI;
- }
- f_[globalIdxI] -= (boundValues[wPhaseIdx] + boundValues[nPhaseIdx]) * faceArea;
-
- //Assumes that the phases flow in the same direction at the neumann boundary, which is the direction of the total flux!!!
- //needed to determine the upwind direction in the saturation equation
- problem_.variables().potentialWetting(globalIdxI, isIndex) = boundValues[wPhaseIdx];
- problem_.variables().potentialNonwetting(globalIdxI, isIndex) = boundValues[nPhaseIdx];
- }
- else
- {
- DUNE_THROW(Dune::NotImplemented, "No valid boundary condition type defined for pressure equation!");
- }
+ densityWIJ = rhoMeanWIJ;
+ densityNWIJ = rhoMeanNWIJ;
+ densityWIK = rhoMeanWIK;
+ densityNWIK = rhoMeanNWIK;
}
- } // end all intersections
- //volume correction due to density differences
- if (compressibility)
- {
- switch (saturationType)
+ switch (pressureType_)
{
- case Sw:
+ case pglobal:
{
- f_[globalIdxI] -= problem_.variables().volumecorrection(globalIdxI) * porosity * volume * (densityWI - densityNWI);
+ Scalar pressJK = (cellDataJ.globalPressure() + cellDataK.globalPressure()) / 2;
+
+ potentialWIJ = (cellDataI.globalPressure() - fMeanNWIJ * pcI - (pressJK - fMeanNWIJ * pcJK)) / l
+ + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng;
+ potentialNWIJ = (cellDataI.globalPressure() + fMeanWIJ * pcI - (pressJK + fMeanWIJ * pcJK)) / l
+ + (densityNWIJ - lJ / l * (kI + kK) / kI * (densityNWIK - densityNWIJ) / 2) * ng;
+ potentialWIK = (cellDataI.globalPressure() - fMeanNWIK * pcI - (pressJK - fMeanNWIK * pcJK)) / l
+ + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng;
+ potentialNWIK = (cellDataI.globalPressure() + fMeanWIK * pcI - (pressJK + fMeanWIK * pcJK)) / l
+ + (densityNWIK - lJ / l * (kI + kK) / kI * (densityNWIJ - densityNWIK) / 2) * ng;
break;
}
- case Sn:
+ default:
{
- f_[globalIdxI] -= problem_.variables().volumecorrection(globalIdxI) * porosity * volume * (densityNWI - densityWI);
+ potentialWIJ = (cellDataI.pressure(wPhaseIdx)
+ - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l
+ + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng;
+ potentialNWIJ = (cellDataI.pressure(nPhaseIdx)
+ - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l
+ + (densityNWIJ - lJ / l * (kI + kK) / kI * (densityNWIK - densityNWIJ) / 2) * ng;
+ potentialWIK = (cellDataI.pressure(wPhaseIdx)
+ - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l
+ + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng;
+ potentialNWIK = (cellDataI.pressure(nPhaseIdx)
+ - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l
+ + (densityNWIK - lJ / l * (kI + kK) / kI * (densityNWIJ - densityNWIK) / 2) * ng;
break;
}
}
}
-// Dune::printmatrix(std::cout,A_,"title","-");
- } // end grid traversal
- return;
-}
-
-//!solves the system of equations to get the spatial distribution of the pressure
-template<class TypeTag>
-void FVPressure2Padaptive<TypeTag>::solve()
-{
- typedef typename GET_PROP_TYPE(TypeTag, LinearSolver) Solver;
-
- int verboseLevelSolver = GET_PARAM(TypeTag, int, LinearSolver, Verbosity);
-
- if (verboseLevelSolver)
- std::cout << "FVPressure2Padaptive: solve for pressure" << std::endl;
-
-// printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
-// printvector(std::cout, f_, "right hand side", "row", 200, 1, 3);
-
- Solver solver(problem_);
- solver.solve(A_, problem_.variables().pressure(), f_);
-
-// printvector(std::cout, (problem_.variables().pressure()), "pressure", "row", 200, 1, 3);
+ //do the upwinding of the mobility depending on the phase potentials
+ Scalar lambdaWIJ = (potentialWIJ > 0.) ? lambdaWI : lambdaWJ;
+ lambdaWIJ = (potentialWIJ == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
+ Scalar lambdaNWIJ = (potentialNWIJ > 0) ? lambdaNWI : lambdaNWJ;
+ lambdaNWIJ = (potentialNWIJ == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNWIJ;
- //for parallel use
-// problem_.variables().communicatePressure();
-
- return;
-}
-//!constitutive functions are updated once if new saturations are calculated and stored in the variables object
-template<class TypeTag>
-void FVPressure2Padaptive<TypeTag>::updateMaterialLaws()
-{
-// //for parallel use
-//// printvector(std::cout, (problem_.variables().pressure()), "pressure", "row", 200, 1, 3);
-//// problem_.variables().communicatePressure();
-//// printvector(std::cout, (problem_.variables().pressure()), "pressureComm", "row", 200, 1, 3);
-//
-//// printvector(std::cout, (problem_.variables().saturation()), "sat", "row", 200, 1, 3);
-// problem_.variables().communicateTransportedQuantity();
-//// printvector(std::cout, (problem_.variables().saturation()), "satComm", "row", 200, 1, 3);
-
- 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)
- {
- int globalIdx = problem_.variables().index(*eIt);
-
- Scalar temperature = problem_.temperature(*eIt);
-
- //determine phase saturations from primary saturation variable
- Scalar satW = 0;
- //Scalar satNW = 0;
- switch (saturationType)
- {
- case Sw:
+ if (compressibility_)
{
- satW = problem_.variables().saturation()[globalIdx];
- //satNW = 1 - problem_.variables().saturation()[globalIdx];
- break;
- }
- case Sn:
- {
- satW = 1 - problem_.variables().saturation()[globalIdx];
- //satNW = problem_.variables().saturation()[globalIdx];
- break;
- }
+ densityWIJ = (potentialWIJ > 0.) ? cellDataI.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
+ densityNWIJ = (potentialNWIJ > 0.) ? cellDataI.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
+ densityWIJ = (potentialWIJ == 0) ? rhoMeanWIJ : densityWIJ;
+ densityNWIJ = (potentialNWIJ == 0) ? rhoMeanNWIJ : densityNWIJ;
+ densityWIK = (potentialWIK > 0.) ? cellDataI.density(wPhaseIdx) : cellDataK.density(wPhaseIdx);
+ densityNWIK = (potentialNWIK > 0.) ? cellDataI.density(nPhaseIdx) : densityNWIK;
+ densityWIK = (potentialWIK == 0) ? rhoMeanWIK : densityWIK;
+ densityNWIK = (potentialNWIK == 0) ? rhoMeanNWIK : densityNWIK;
}
- problem_.variables().capillaryPressure(globalIdx) = MaterialLaw::pC(
- problem_.spatialParameters().materialLawParams(*eIt), satW);
+ // update diagonal entry and right hand side
+ entries[matrix] = (lambdaWIJ + lambdaNWIJ) * kMean / l * faceArea;
+ entries[rhs] = faceArea * lambdaWIJ * kMean * ng
+ * ((densityWIJ) - (lJ / l) * (kI + kK) / kI * (densityWIK - densityWIJ) / 2);
+ entries[rhs] += faceArea * lambdaNWIJ * kMean * ng
+ * ((densityNWIJ) - (lJ / l) * (kI + kK) / kI * (densityNWIK - densityNWIJ) / 2);
- //determine phase pressures from primary pressure variable
- Scalar pressW = 0;
- Scalar pressNW = 0;
- switch (pressureType)
+ switch (pressureType_)
{
case pw:
{
- pressW = problem_.variables().pressure()[globalIdx];
- pressNW = problem_.variables().pressure()[globalIdx] + problem_.variables().capillaryPressure(globalIdx);
+ entries[rhs] += faceArea * lambdaNWIJ * kMean * (pcJK - pcI) / l;
break;
}
case pn:
{
- pressW = problem_.variables().pressure()[globalIdx] - problem_.variables().capillaryPressure(globalIdx);
- pressNW = problem_.variables().pressure()[globalIdx];
+ entries[rhs] -= faceArea * lambdaWIJ * kMean * (pcJK - pcI) / l;
break;
}
}
- Scalar densityW = 0;
- Scalar densityNW = 0;
- Scalar viscosityW = 0;
- Scalar viscosityNW = 0;
+ //write hanging-node-specific stuff directly into matrix and rhs!
+ this->f_[globalIdxI] -= entries[rhs];
+ this->f_[globalIdxJ] += entries[rhs];
- if (compressibility)
- {
- fluidState.update(satW, pressW, pressNW, temperature);
- }
- else
- {
- pressW = problem_.referencePressure(*eIt);
- pressNW = problem_.referencePressure(*eIt);
- fluidState.update(satW, pressW, pressNW, temperature);
- }
+ // set diagonal entry
+ this->A_[globalIdxI][globalIdxI] += entries[matrix];
+ this->A_[globalIdxI][globalIdxJ] -= 0.5 * entries[matrix];
+ this->A_[globalIdxI][globalIdxK] -= 0.5 * entries[matrix];
- densityW = FluidSystem::phaseDensity(wPhaseIdx, temperature, pressW, fluidState);
- densityNW = FluidSystem::phaseDensity(nPhaseIdx, temperature, pressNW, fluidState);
+ // set off-diagonal entry
+ this->A_[globalIdxJ][globalIdxI] -= entries[matrix];
+ this->A_[globalIdxJ][globalIdxJ] += 0.5 * entries[matrix];
+ this->A_[globalIdxJ][globalIdxK] += 0.5 * entries[matrix];
- viscosityW = FluidSystem::phaseViscosity(wPhaseIdx, temperature, pressW, fluidState);
- viscosityNW = FluidSystem::phaseViscosity(nPhaseIdx, temperature, pressNW, fluidState);
+ //set entries to zero -> matrix already written!
+ entries = 0.;
- // initialize mobilities
- Scalar mobilityW = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosityW;
- Scalar mobilityNW = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*eIt), satW) / viscosityNW;
-
- if (compressibility)
- {
- mobilityW *= densityW;
- mobilityNW *= densityNW;
- }
-
- // initialize mobilities
- problem_.variables().mobilityWetting(globalIdx) = mobilityW;
- problem_.variables().mobilityNonwetting(globalIdx) = mobilityNW;
-
- // initialize densities
- problem_.variables().densityWetting(globalIdx) = densityW;
- problem_.variables().densityNonwetting(globalIdx) = densityNW;
-
- // initialize viscosities
- problem_.variables().viscosityWetting(globalIdx) = viscosityW;
- problem_.variables().viscosityNonwetting(globalIdx) = viscosityNW;
-
- //initialize fractional flow functions
- problem_.variables().fracFlowFuncWetting(globalIdx) = mobilityW / (mobilityW + mobilityNW);
- problem_.variables().fracFlowFuncNonwetting(globalIdx) = mobilityNW / (mobilityW + mobilityNW);
+// std::cout<<"finished hanging node!\n";
}
+
return;
}
diff --git a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
index 0ebc8dfaabd4c29169853515422e931ca6910099..9063b2871ba0f1c110f959aa8f4156a938f74742 100644
--- a/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
+++ b/dumux/decoupled/2p/diffusion/fv/fvvelocity2padaptive.hh
@@ -28,8 +28,7 @@
* @author Markus Wolff
*/
-#include <dumux/decoupled/2p/diffusion/fv/fvpressure2padaptive.hh>
-#include <dune/grid/common/gridenums.hh>
+#include <dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh>
namespace Dumux
{
@@ -51,14 +50,14 @@ namespace Dumux
*/
template<class TypeTag>
-class FVVelocity2Padaptive: public FVPressure2Padaptive<TypeTag>
+class FVVelocity2Padaptive: public FVVelocity2P<TypeTag>
{
- typedef FVVelocity2Padaptive<TypeTag> ThisType;
- typedef FVPressure2Padaptive<TypeTag> ParentType;
+ typedef FVVelocity2P<TypeTag> ParentType;
+
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, Variables) Variables;
+
typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters;
typedef typename SpatialParameters::MaterialLaw MaterialLaw;
@@ -71,13 +70,13 @@ class FVVelocity2Padaptive: public FVPressure2Padaptive<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
- typedef typename GridView::Grid Grid;
- typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
- typedef typename Grid::template Codim<0>::EntityPointer ElementPointer;
+ typedef typename GridView::Intersection Intersection;
+ typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
enum
{
@@ -85,8 +84,6 @@ typedef typename GridView::Traits::template Codim<0>::Entity Element;
};
typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElementContainer;
- typedef Dune::GenericReferenceElements<Scalar, dim-1> ReferenceElementFaceContainer;
- typedef Dune::GenericReferenceElement<Scalar, dim> ReferenceElement;
enum
{
@@ -105,7 +102,7 @@ typedef typename GridView::Traits::template Codim<0>::Entity Element;
};
enum
{
- wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx, numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
};
typedef Dune::FieldVector<Scalar,dimWorld> GlobalPosition;
@@ -117,7 +114,7 @@ public:
* \param problem a problem class object
*/
FVVelocity2Padaptive(Problem& problem)
- : FVPressure2Padaptive<TypeTag>(problem)
+ : ParentType(problem), problem_(problem), gravity_(problem.gravity())
{
// todo: kompatibilität prüfen
if (GET_PROP_VALUE(TypeTag, EnableCompressibility) && velocityType_ == vt)
@@ -128,6 +125,21 @@ public:
{
DUNE_THROW(Dune::NotImplemented, "Velocity type not supported!");
}
+
+ if (!compressibility_)
+ {
+ const Element& element = *(problem_.gridView().template begin<0> ());
+ FluidState fluidState;
+ fluidState.setPressure(wPhaseIdx, problem_.referencePressure(element));
+ fluidState.setPressure(nPhaseIdx, problem_.referencePressure(element));
+ fluidState.setTemperature(problem_.temperature(element));
+ fluidState.setSaturation(wPhaseIdx, 1.);
+ fluidState.setSaturation(nPhaseIdx, 0.);
+ density_[wPhaseIdx] = FluidSystem::density(fluidState, wPhaseIdx);
+ density_[nPhaseIdx] = FluidSystem::density(fluidState, nPhaseIdx);
+ viscosity_[wPhaseIdx] = FluidSystem::viscosity(fluidState, wPhaseIdx);
+ viscosity_[nPhaseIdx] = FluidSystem::viscosity(fluidState, nPhaseIdx);
+ }
}
//! Calculate the velocity.
@@ -138,954 +150,324 @@ public:
* The method is needed in the IMPES (Implicit Pressure Explicit Saturation) algorithm which is used for a fractional flow formulation
* to provide the velocity field required for the solution of the saturation equation.
*/
- void calculateVelocity();
+ void calculateVelocity(const Intersection& intersection, CellData& cellDataI);
- void update()
- {
- ParentType::update();
+private:
+ Problem& problem_;
+ const GlobalPosition& gravity_; //!< vector including the gravity constant
+ Scalar density_[numPhases];
+ Scalar viscosity_[numPhases];
- calculateVelocity();
+ static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation); //!< gives kind of velocity used (\f$ 0 = v_w\f$, \f$ 1 = v_n\f$, \f$ 2 = v_t\f$)
+ static const bool compressibility_ = GET_PROP_VALUE(TypeTag, EnableCompressibility);
+ static const int pressureType_ = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$p_w\f$, \f$p_n\f$, \f$p_{global}\f$)
+ static const int saturationType_ = GET_PROP_VALUE(TypeTag, SaturationFormulation); //!< gives kind of saturation used (\f$S_w\f$, \f$S_n\f$)
+};
- return;
- }
+template<class TypeTag>
+void FVVelocity2Padaptive<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellDataI)
+{
+ ElementPointer elementI = intersection.inside();
+ ElementPointer elementJ = intersection.outside();
- //! \brief Write data files
- /* \param name file name */
- template<class MultiWriter>
- void addOutputVtkFields(MultiWriter &writer)
+ if (elementI->level() == elementJ->level())
{
- ParentType::addOutputVtkFields(writer);
-
- Dune::BlockVector<Dune::FieldVector<Scalar, dim> > &velocity = *(writer.template allocateManagedBuffer<Scalar, dim> (
- this->problem().gridView().size(0)));
- Dune::BlockVector<Dune::FieldVector<Scalar, dim> > &velocitySecondPhase = *(writer.template allocateManagedBuffer<Scalar, dim> (
- this->problem().gridView().size(0)));
-
- // compute update vector
- ElementIterator eItEnd = this->problem().gridView().template end<0>();
- for (ElementIterator eIt = this->problem().gridView().template begin<0>(); eIt != eItEnd; ++eIt)
+ ParentType::calculateVelocity(intersection, cellDataI);
+ }
+ else if (elementJ->level() == elementI->level() + 1)
+ {
+ CellData& cellDataJ = problem_.variables().cellData(problem_.variables().index(*elementJ));
+
+ // get global coordinates of cell centers
+ const GlobalPosition& globalPosI = elementI->geometry().center();
+ const GlobalPosition& globalPosJ = elementJ->geometry().center();
+
+ // get mobilities and fractional flow factors
+ Scalar lambdaWI = cellDataI.mobility(wPhaseIdx);
+ Scalar lambdaNWI = cellDataI.mobility(nPhaseIdx);
+ Scalar fractionalWI = cellDataI.fracFlowFunc(wPhaseIdx);
+ Scalar fractionalNWI = cellDataI.fracFlowFunc(nPhaseIdx);
+ Scalar lambdaWJ = cellDataJ.mobility(wPhaseIdx);
+ Scalar lambdaNWJ = cellDataJ.mobility(nPhaseIdx);
+ Scalar fractionalWJ = cellDataJ.fracFlowFunc(wPhaseIdx);
+ Scalar fractionalNWJ = cellDataJ.fracFlowFunc(nPhaseIdx);
+
+ // get capillary pressure
+ Scalar pcI = cellDataI.capillaryPressure();
+ Scalar pcJ = cellDataJ.capillaryPressure();
+
+ //get face index
+ int isIndexI = intersection.indexInInside();
+
+ //get face normal
+ const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
+
+ // get face area
+ Scalar faceArea = intersection.geometry().volume();
+
+ // distance vector between barycenters
+ GlobalPosition distVecIJ = globalPosJ - globalPosI;
+
+ // compute distance between cell centers
+ Scalar distIJ = distVecIJ.two_norm();
+
+ // Count number of hanging nodes
+ // not really necessary
+ int isIndexJ = intersection.indexInOutside();
+
+ int globalIdxK = 0;
+ ElementPointer elementK = intersection.outside();
+ // We are looking for two things:
+ // IsIndexJ, the index of the interface from the neighbor-cell point of view
+ // GlobalIdxK, the index of the third cell
+ // for efficienty this is done in one IntersectionIterator-Loop
+
+ // Intersectioniterator around cell I
+ IntersectionIterator isItEndI = problem_.gridView().iend(*elementI);
+ for (IntersectionIterator isItI = problem_.gridView().ibegin(*elementI); isItI != isItEndI; ++isItI)
{
- // cell index
- int globalIdx = this->problem().variables().index(*eIt);
-
- Dune::FieldVector<Scalar, 2*dim> fluxW(0);
- Dune::FieldVector<Scalar, 2*dim> fluxNW(0);
- // run through all intersections with neighbors and boundary
- IntersectionIterator
- isItEnd = this->problem().gridView().iend(*eIt);
- for (IntersectionIterator
- isIt = this->problem().gridView().ibegin(*eIt); isIt
- !=isItEnd; ++isIt)
+ if (isItI->neighbor())
{
- int isIndex = isIt->indexInInside();
+ ElementPointer neighborPointer2 = isItI->outside();
- fluxW[isIndex] = isIt->geometry().volume() * (isIt->centerUnitOuterNormal() * this->problem().variables().velocityElementFace(*eIt, isIndex));
- fluxNW[isIndex] = isIt->geometry().volume() * (isIt->centerUnitOuterNormal() * this->problem().variables().velocitySecondPhase()[this->problem().variables().index(*eIt)][isIndex]);
+ // make sure we do not chose elemntI as third element
+ // -> faces with hanging node have more than one intersection but only one face index!
+ if (neighborPointer2 != elementJ && isItI->indexInInside() == isIndexI)
+ {
+ globalIdxK = problem_.variables().index(*neighborPointer2);
+ elementK = neighborPointer2;
+ break;
+ }
}
+ }
- Dune::FieldVector<Scalar, dim> refVelocity(0);
- refVelocity[0] = 0.5 * (fluxW[1] - fluxW[0]);
- refVelocity[1] = 0.5 * (fluxW[3] - fluxW[2]);
-
- const Dune::FieldVector<Scalar, dim>& localPos = ReferenceElementContainer::general(eIt->geometry().type()).position(0,
- 0);
-
- // get the transposed Jacobian of the element mapping
- const FieldMatrix& jacobianInv = eIt->geometry().jacobianInverseTransposed(localPos);
- FieldMatrix jacobianT(jacobianInv);
- jacobianT.invert();
+ CellData& cellDataK = problem_.variables().cellData(globalIdxK);
+
+ // face global coordinates
+ const GlobalPosition& globalPosInterface = intersection.geometry().center();
+
+ Dune::FieldVector < Scalar, dimWorld > distVec = globalPosInterface - globalPosI;
+ Scalar lI = distVec * unitOuterNormal;
+ distVec = globalPosJ - globalPosInterface;
+ Scalar lJ = distVec * unitOuterNormal;
+ Scalar l = lI + lJ;
+
+ FieldMatrix permeabilityI(0);
+ FieldMatrix permeabilityJ(0);
+ FieldMatrix permeabilityK(0);
+
+ problem_.spatialParameters().meanK(permeabilityI,
+ problem_.spatialParameters().intrinsicPermeability(*elementI));
+ problem_.spatialParameters().meanK(permeabilityJ,
+ problem_.spatialParameters().intrinsicPermeability(*elementJ));
+ problem_.spatialParameters().meanK(permeabilityK,
+ problem_.spatialParameters().intrinsicPermeability(*elementK));
+
+ // Calculate permeablity component normal to interface
+ Scalar kI, kJ, kK, ng, kMean; //, gI, gJ, gK;
+ Dune::FieldVector < Scalar, dim > permI(0);
+ Dune::FieldVector < Scalar, dim > permJ(0);
+ Dune::FieldVector < Scalar, dim > permK(0);
+
+ permeabilityI.mv(unitOuterNormal, permI);
+ permeabilityJ.mv(unitOuterNormal, permJ);
+ permeabilityK.mv(unitOuterNormal, permK);
+
+ // kI,kJ,kK = (n^T)Kn
+ kI = unitOuterNormal * permI;
+ kJ = unitOuterNormal * permJ;
+ kK = unitOuterNormal * permK;
+ kMean = kI * kJ * kK * l / (kJ * kK * lI + kI * (kJ + kK) / 2 * lJ);
+
+ ng = gravity_ * unitOuterNormal;
+
+ Scalar lambdaWK = cellDataK.mobility(wPhaseIdx);
+ Scalar lambdaNWK = cellDataK.mobility(nPhaseIdx);
+ Scalar fractionalWK = cellDataK.fracFlowFunc(wPhaseIdx);
+ Scalar fractionalNWK = cellDataK.fracFlowFunc(nPhaseIdx);
+
+ Scalar pcK = cellDataK.capillaryPressure();
+ Scalar pcJK = (pcJ + pcK) / 2;
+
+ // calculate potential gradients
+ // reuse potentials from fvpressure2padaptive
+ Scalar potentialWIJ = cellDataI.fluxData().potential(wPhaseIdx, isIndexI);
+ Scalar potentialNWIJ = cellDataI.fluxData().potential(nPhaseIdx, isIndexI);
+ Scalar potentialWIK = potentialWIJ;
+ Scalar potentialNWIK = potentialNWIJ;
+ // preliminary potential. The "real" ones are found below
+
+ // Comment: reversed weighting is plausible, too (swap lJ and lI)
+ Scalar rhoMeanWIJ = density_[wPhaseIdx];
+ Scalar rhoMeanNWIJ = density_[nPhaseIdx];
+ Scalar rhoMeanWIK = density_[wPhaseIdx];
+ Scalar rhoMeanNWIK = density_[nPhaseIdx];
+
+ if (compressibility_)
+ {
+ rhoMeanWIJ = (lJ * cellDataI.density(wPhaseIdx) + lI * cellDataJ.density(wPhaseIdx)) / l;
+ rhoMeanNWIJ = (lJ * cellDataI.density(nPhaseIdx) + lI * cellDataJ.density(nPhaseIdx)) / l;
+ rhoMeanWIK = (lJ * cellDataI.density(wPhaseIdx) + lI * cellDataK.density(wPhaseIdx)) / l;
+ rhoMeanNWIK = (lJ * cellDataI.density(nPhaseIdx) + lI * cellDataK.density(nPhaseIdx)) / l;
+ }
- // calculate the element velocity by the Piola transformation
- Dune::FieldVector<Scalar, dim> elementVelocity(0);
- jacobianT.umtv(refVelocity, elementVelocity);
- elementVelocity /= eIt->geometry().integrationElement(localPos);
+ Scalar fMeanWIJ = (lJ * fractionalWI + lI * fractionalWJ) / l;
+ Scalar fMeanNWIJ = (lJ * fractionalNWI + lI * fractionalNWJ) / l;
+ Scalar fMeanWIK = (lJ * fractionalWI + lI * fractionalWK) / l;
+ Scalar fMeanNWIK = (lJ * fractionalNWI + lI * fractionalNWK) / l;
- velocity[globalIdx] = elementVelocity;
+ Scalar densityWIJ = density_[wPhaseIdx];
+ Scalar densityNWIJ = density_[nPhaseIdx];
+ Scalar densityWIK = density_[wPhaseIdx];
+ Scalar densityNWIK = density_[nPhaseIdx];
- refVelocity = 0;
- refVelocity[0] = 0.5 * (fluxNW[1] - fluxNW[0]);
- refVelocity[1] = 0.5 * (fluxNW[3] - fluxNW[2]);
+ if (compressibility_)
+ {
+ // Upwinding for finding the upwinding direction
+ densityWIJ = (potentialWIJ > 0.) ? cellDataI.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
+ densityNWIJ = (potentialNWIJ > 0.) ? cellDataI.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
+ densityWIK = (potentialWIK > 0.) ? cellDataI.density(wPhaseIdx) : cellDataK.density(wPhaseIdx);
+ densityNWIK = (potentialNWIK > 0.) ? cellDataI.density(nPhaseIdx) : cellDataK.density(nPhaseIdx);
+
+ densityWIJ = (potentialWIJ == 0.) ? rhoMeanWIJ : densityWIJ;
+ densityNWIJ = (potentialNWIJ == 0.) ? rhoMeanNWIJ : densityNWIJ;
+ densityWIK = (potentialWIK == 0.) ? rhoMeanWIK : densityWIK;
+ densityNWIK = (potentialNWIK == 0.) ? rhoMeanNWIK : densityNWIK;
+ }
- // calculate the element velocity by the Piola transformation
- elementVelocity = 0;
- jacobianT.umtv(refVelocity, elementVelocity);
- elementVelocity /= eIt->geometry().integrationElement(localPos);
+ Scalar fractionalWIJ = (potentialWIJ > 0.) ? fractionalWI : fractionalWJ;
+ Scalar fractionalNWIJ = (potentialNWIJ > 0.) ? fractionalNWI : fractionalNWJ;
+ Scalar fractionalWIK = (potentialWIK > 0.) ? fractionalWI : fractionalWK;
+ Scalar fractionalNWIK = (potentialNWIK > 0.) ? fractionalNWI : fractionalNWK;
- velocitySecondPhase[globalIdx] = elementVelocity;
- }
+ fractionalWIJ = (potentialWIJ == 0.) ? fMeanWIJ : fractionalWIJ;
+ fractionalNWIJ = (potentialNWIJ == 0.) ? fMeanNWIJ : fractionalNWIJ;
+ fractionalWIK = (potentialWIK == 0.) ? fMeanWIK : fractionalWIK;
+ fractionalNWIK = (potentialNWIK == 0.) ? fMeanNWIK : fractionalNWIK;
- //switch velocities
- switch (velocityType_)
+ switch (pressureType_)
{
- case vw:
+ case pglobal:
{
- writer.attachCellData(velocity, "wetting-velocity", dim);
- writer.attachCellData(velocitySecondPhase, "non-wetting-velocity", dim);
+ Scalar pressJK = (cellDataJ.globalPressure() + cellDataK.globalPressure()) / 2;
+
+ potentialWIJ = (cellDataI.globalPressure() - fMeanNWIJ * pcI - (pressJK - fMeanNWIJ * pcJK)) / l
+ + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng;
+ potentialNWIJ = (cellDataI.globalPressure() + fMeanWIJ * pcI - (pressJK + fMeanWIJ * pcJK)) / l
+ + (densityNWIJ - lJ / l * (kI + kK) / kI * (densityNWIK - densityNWIJ) / 2) * ng;
+ potentialWIK = (cellDataI.globalPressure() - fMeanNWIK * pcI - (pressJK - fMeanNWIK * pcJK)) / l
+ + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng;
+ potentialNWIK = (cellDataI.globalPressure() + fMeanWIK * pcI - (pressJK + fMeanWIK * pcJK)) / l
+ + (densityNWIK - lJ / l * (kI + kK) / kI * (densityNWIJ - densityNWIK) / 2) * ng;
break;
}
- case vn:
+ default:
{
- writer.attachCellData(velocity, "non-wetting-velocity", dim);
- writer.attachCellData(velocitySecondPhase, "wetting-velocity", dim);
+ potentialWIJ = (cellDataI.pressure(wPhaseIdx)
+ - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l
+ + (densityWIJ - lJ / l * (kI + kK) / kI * (densityWIK - densityWIJ) / 2) * ng;
+ potentialNWIJ = (cellDataI.pressure(nPhaseIdx)
+ - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l
+ + (densityNWIJ - lJ / l * (kI + kK) / kI * (densityNWIK - densityNWIJ) / 2) * ng;
+ potentialWIK = (cellDataI.pressure(wPhaseIdx)
+ - 0.5 * (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx))) / l
+ + (densityWIK - lJ / l * (kI + kK) / kI * (densityWIJ - densityWIK) / 2) * ng;
+ potentialNWIK = (cellDataI.pressure(nPhaseIdx)
+ - (0.5 * (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)))) / l
+ + (densityNWIK - lJ / l * (kI + kK) / kI * (densityNWIJ - densityNWIK) / 2) * ng;
break;
}
- case vt:
- {
- writer.attachCellData(velocity, "total velocity", dim);
- break;
}
+
+ //store potentials for further calculations (velocity, saturation, ...)
+ // these quantities only have correct sign (needed for upwinding)
+ // potentials are defined slightly different for adaptive scheme
+ cellDataI.fluxData().setPotential(wPhaseIdx, isIndexI, potentialWIJ);
+ cellDataI.fluxData().setPotential(nPhaseIdx, isIndexI, potentialNWIJ);
+ cellDataJ.fluxData().setPotential(wPhaseIdx, isIndexJ, -potentialWIJ);
+ cellDataJ.fluxData().setPotential(nPhaseIdx, isIndexJ, -potentialNWIJ);
+
+ //do the upwinding of the mobility depending on the phase potentials
+ Scalar lambdaWIJ = (potentialWIJ > 0.) ? lambdaWI : lambdaWJ;
+ lambdaWIJ = (potentialWIJ == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
+ Scalar lambdaNWIJ = (potentialNWIJ > 0.) ? lambdaNWI : lambdaNWJ;
+ lambdaNWIJ = (potentialNWIJ == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNWIJ;
+
+ if (compressibility_)
+ {
+ densityWIJ = (potentialWIJ > 0.) ? cellDataI.density(wPhaseIdx) : cellDataJ.density(wPhaseIdx);
+ densityNWIJ = (potentialNWIJ > 0.) ? cellDataI.density(nPhaseIdx) : cellDataJ.density(nPhaseIdx);
+ densityWIJ = (potentialWIJ == 0) ? rhoMeanWIJ : densityWIJ;
+ densityNWIJ = (potentialNWIJ == 0) ? rhoMeanNWIJ : densityNWIJ;
+ densityWIK = (potentialWIK > 0.) ? cellDataI.density(wPhaseIdx) : cellDataK.density(wPhaseIdx);
+ densityNWIK = (potentialNWIK > 0.) ? cellDataI.density(nPhaseIdx) : cellDataK.density(nPhaseIdx);
+ densityWIK = (potentialWIK == 0) ? rhoMeanWIK : densityWIK;
+ densityNWIK = (potentialNWIK == 0) ? rhoMeanNWIK : densityNWIK;
}
- return;
- }
+ fractionalWIJ = (potentialWIJ > 0.) ? fractionalWI : fractionalWJ;
+ fractionalNWIJ = (potentialNWIJ > 0.) ? fractionalNWI : fractionalNWJ;
+ fractionalWIK = (potentialWIK > 0.) ? fractionalWI : fractionalWK;
+ fractionalNWIK = (potentialNWIK > 0.) ? fractionalNWI : fractionalNWK;
-private:
- static const int velocityType_ = GET_PROP_VALUE(TypeTag, VelocityFormulation); //!< gives kind of velocity used (\f$ 0 = v_w\f$, \f$ 1 = v_n\f$, \f$ 2 = v_t\f$)
-};
+ fractionalWIJ = (potentialWIJ == 0.) ? fMeanWIJ : fractionalWIJ;
+ fractionalNWIJ = (potentialNWIJ == 0.) ? fMeanNWIJ : fractionalNWIJ;
+ fractionalWIK = (potentialWIK == 0.) ? fMeanWIK : fractionalWIK;
+ fractionalNWIK = (potentialNWIK == 0.) ? fMeanNWIK : fractionalNWIK;
-template<class TypeTag>
-void FVVelocity2Padaptive<TypeTag>::calculateVelocity()
-{
- //reset velocity
- this->problem().variables().velocity() = Dune::FieldVector<Scalar, dimWorld>(0.0);
- this->problem().variables().velocitySecondPhase() = Dune::FieldVector<Scalar, dimWorld>(0.0);
+ //calculate velocities and the gravity term
+ Dune::FieldVector < Scalar, dimWorld > velocityW(unitOuterNormal);
+ Dune::FieldVector < Scalar, dimWorld > velocityNW(unitOuterNormal);
+ Dune::FieldVector < Scalar, dimWorld > gravityTermW(unitOuterNormal);
+ Dune::FieldVector < Scalar, dimWorld > gravityTermNW(unitOuterNormal);
- BoundaryTypes bcType;
+ gravityTermW *= lambdaWIJ * kMean * ng;
+ gravityTermW *= densityWIJ - (lJ / l) * (kI + kK) / kI * (densityWIK - densityWIJ) / 2;
+ gravityTermNW *= lambdaNWIJ * kMean * ng;
+ gravityTermNW *= densityNWIJ - (lJ / l) * (kI + kK) / kI * (densityNWIK - densityNWIJ) / 2;
- // compute update vector
- ElementIterator eItEnd = this->problem().gridView().template end<0>();
- for (ElementIterator eIt = this->problem().gridView().template begin<0>(); eIt != eItEnd; ++eIt)
- {
-#if HAVE_MPI
- if (eIt->partitionType() == Dune::GhostEntity || eIt->partitionType() == Dune::OverlapEntity)
+ switch (pressureType_)
{
- continue;
- }
-#endif
-
- // cell index
- int globalIdxI = this->problem().variables().index(*eIt);
-
- GlobalPosition globalPos = eIt->geometry().center();
-
- Scalar pressI = this->problem().variables().pressure()[globalIdxI];
- Scalar pcI = this->problem().variables().capillaryPressure(globalIdxI);
- Scalar lambdaWI = this->problem().variables().mobilityWetting(globalIdxI);
- Scalar lambdaNWI = this->problem().variables().mobilityNonwetting(globalIdxI);
- Scalar fractionalWI = this->problem().variables().fracFlowFuncWetting(globalIdxI);
- Scalar fractionalNWI = this->problem().variables().fracFlowFuncNonwetting(globalIdxI);
- Scalar densityWI = this->problem().variables().densityWetting(globalIdxI);
- Scalar densityNWI = this->problem().variables().densityNonwetting(globalIdxI);
-
- // run through all intersections with neighbors and boundary
- IntersectionIterator isItEnd = this->problem().gridView().iend(*eIt);
- for (IntersectionIterator
- isIt = this->problem().gridView().ibegin(*eIt); isIt
- !=isItEnd; ++isIt)
+ case pw:
{
- // local number of facet
- int isIndex = isIt->indexInInside();
+ velocityW *= lambdaWIJ * kMean * (cellDataI.pressure(wPhaseIdx) - (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx)) / 2.0) / l;
+ velocityNW *= lambdaNWIJ * kMean * (cellDataI.pressure(nPhaseIdx) - (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)) / 2.0) / l;
- Dune::FieldVector<Scalar,dimWorld> unitOuterNormal = isIt->centerUnitOuterNormal();
+ velocityW += gravityTermW;
+ velocityNW += gravityTermNW;
+ break;
+ }
+ case pn:
+ {
+ velocityNW *= lambdaNWIJ * kMean * (cellDataI.pressure(nPhaseIdx) - (cellDataJ.pressure(nPhaseIdx) + cellDataK.pressure(nPhaseIdx)) / 2.0) / l;
+ velocityW *= lambdaWIJ * kMean * (cellDataI.pressure(wPhaseIdx) - (cellDataJ.pressure(wPhaseIdx) + cellDataK.pressure(wPhaseIdx)) / 2.0) / l;
- // handle interior face
- if (isIt->neighbor())
- {
- // access neighbor
- ElementPointer neighborPointer = isIt->outside();
- int globalIdxJ = this->problem().variables().index(*neighborPointer);
+ velocityW += gravityTermW;
+ velocityNW += gravityTermNW;
+ break;
+ }
+ case pglobal:
+ {
+ Scalar pressJK = (cellDataJ.globalPressure() + cellDataK.globalPressure()) / 2;
- if (neighborPointer->level()==eIt.level())
- {
- // neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
-
- // distance vector between barycenters
- Dune::FieldVector<Scalar,dimWorld> distVec = globalPosNeighbor - globalPos;
-
- // compute distance between cell centers
- Scalar dist = distVec * unitOuterNormal;
-// Scalar dist = distVec.two_norm();
-
- // compute vectorized permeabilities
- FieldMatrix meanPermeability(0);
-
- this->problem().spatialParameters().meanK(meanPermeability,
- this->problem().spatialParameters().intrinsicPermeability(*eIt),
- this->problem().spatialParameters().intrinsicPermeability(*neighborPointer));
-
- Dune::FieldVector<Scalar, dim> permeability(0);
- meanPermeability.mv(unitOuterNormal, permeability);
-
- Scalar pressJ = this->problem().variables().pressure()[globalIdxJ];
- Scalar pcJ = this->problem().variables().capillaryPressure(globalIdxJ);
- Scalar lambdaWJ = this->problem().variables().mobilityWetting(globalIdxJ);
- Scalar lambdaNWJ = this->problem().variables().mobilityNonwetting(globalIdxJ);
- Scalar fractionalWJ = this->problem().variables().fracFlowFuncWetting(globalIdxJ);
- Scalar fractionalNWJ = this->problem().variables().fracFlowFuncNonwetting(globalIdxJ);
- Scalar densityWJ = this->problem().variables().densityWetting(globalIdxJ);
- Scalar densityNWJ = this->problem().variables().densityNonwetting(globalIdxJ);
-
- //calculate potential gradients
- Scalar potentialW = 0;
- Scalar potentialNW = 0;
-
- potentialW = this->problem().variables().potentialWetting(globalIdxI, isIndex);
- potentialNW = this->problem().variables().potentialNonwetting(globalIdxI, isIndex);
-
- Scalar densityW = (potentialW> 0.) ? densityWI : densityWJ;
- Scalar densityNW = (potentialNW> 0.) ? densityNWI : densityNWJ;
-
- densityW = (potentialW == 0.) ? 0.5 * (densityWI + densityWJ) : densityW;
- densityNW = (potentialNW == 0.) ? 0.5 * (densityNWI + densityNWJ) : densityNW;
-
- switch (this->pressureType)
- {
- case pw:
- {
- potentialW = (pressI - pressJ);
- potentialNW = (pressI - pressJ+ pcI - pcJ);
- break;
- }
- case pn:
- {
- potentialW = (pressI - pressJ - pcI + pcJ);
- potentialNW = (pressI - pressJ);
- break;
- }
- case pglobal:
- {
- potentialW = (pressI - pressJ - 0.5 * (fractionalNWI+fractionalNWJ)*(pcI - pcJ));
- potentialNW = (pressI - pressJ + 0.5 * (fractionalWI+fractionalWJ)*(pcI - pcJ));
- break;
- }
- }
-
- potentialW += densityW * (distVec * this->gravity);//delta z/delta x in unitOuterNormal[z]
- potentialNW += densityNW * (distVec * this->gravity);
-
- //store potentials for further calculations (velocity, saturation, ...)
- this->problem().variables().potentialWetting(globalIdxI, isIndex) = potentialW;
- this->problem().variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
-
- //do the upwinding of the mobility depending on the phase potentials
- Scalar lambdaW = (potentialW > 0.) ? lambdaWI : lambdaWJ;
- lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaW;
- Scalar lambdaNW = (potentialNW > 0.) ? lambdaNWI : lambdaNWJ;
- lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNW;
- densityW = (potentialW> 0.) ? densityWI : densityWJ;
- densityW = (potentialW == 0.) ? 0.5 * (densityWI + densityWJ) : densityW;
- densityNW = (potentialNW> 0.) ? densityNWI : densityNWJ;
- densityNW = (potentialNW == 0.) ? 0.5 * (densityNWI + densityNWJ) : densityNW;
-
- //calculate the gravity term
- Dune::FieldVector<Scalar,dimWorld> velocityW(permeability);
- Dune::FieldVector<Scalar,dimWorld> velocityNW(permeability);
- Dune::FieldVector<Scalar,dimWorld> gravityTermW(unitOuterNormal);
- Dune::FieldVector<Scalar,dimWorld> gravityTermNW(unitOuterNormal);
-
- gravityTermW *= (this->gravity*permeability)*(lambdaW * densityW);
- gravityTermNW *= (this->gravity*permeability)*(lambdaNW * densityNW);
-
- //calculate velocity depending on the pressure used -> use pc = pn - pw
- switch (this->pressureType)
- {
- case pw:
- {
- velocityW *= lambdaW * (pressI - pressJ)/dist;
- velocityNW *= lambdaNW * (pressI - pressJ)/dist + 0.5 * (lambdaNWI + lambdaNWJ) * (pcI - pcJ) / dist;
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- case pn:
- {
- velocityW *= lambdaW * (pressI - pressJ)/dist - 0.5 * (lambdaWI + lambdaWJ) * (pcI - pcJ) / dist;
- velocityNW *= lambdaNW * (pressI - pressJ) / dist;
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- case pglobal:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = permeability;
- this->problem().variables().velocity()[globalIdxI][isIndex]*= (lambdaW + lambdaNW)* (pressI - pressJ )/ dist;
- this->problem().variables().velocity()[globalIdxI][isIndex] += gravityTermW;
- this->problem().variables().velocity()[globalIdxI][isIndex] += gravityTermNW;
- break;
- }
- }
-
- //store velocities
- switch (velocityType_)
- {
- case vw:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
- break;
- }
- case vn:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
- break;
- }
- case vt:
- {
- switch (this->pressureType)
- {
- case pw:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
- break;
- }
- case pn:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
- break;
- }
- }
- break;
- }
- } //end of switch (velocityType_)
- } //End of case "same level"
-
- if (neighborPointer->level()==eIt.level()+1)
- {
- int globalIdxK = 0;
- ElementPointer thirdCellPointer = isIt->outside();
- bool foundK=false;
- bool foundIJ=false;
- // We are looking for two things:
- // IsIndexJ, the index of the interface from the neighbor-cell point of view
- // GlobalIdxK, the index of the third cell
- // for efficienty this is done in one IntersectionIterator-Loop
-
- // Intersectioniterator around cell J
- IntersectionIterator isItEndJ = this->problem().gridView().iend(*neighborPointer);
- for (IntersectionIterator isItJ = this->problem().gridView().ibegin(*neighborPointer); isItJ != isItEndJ; ++isItJ)
- {
- if (isItJ->neighbor())
- {
- ElementPointer neighborPointer2 = isItJ->outside();
-
- // Neighbor of neighbor is Cell I?
- if (this->problem().variables().index(*neighborPointer2)==globalIdxI)
- {
- foundIJ=true;
- }
- else
- {
- if (neighborPointer2->level()==eIt.level()+1)
- {
- // To verify, we found the correct Cell K, we check
- // - is level(K)=level(J)?
- // - is distance(IJ)=distance(IK)?
- // - K is neighbor of J.
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
- const GlobalPosition& globalPosThirdCell = neighborPointer2->geometry().center();
- Dune::FieldVector<Scalar, dimWorld> distVecIJ = globalPosNeighbor - globalPos;
- Dune::FieldVector<Scalar, dimWorld> distVecIK = globalPosThirdCell - globalPos;
- if (((distVecIK.two_norm()-distVecIJ.two_norm()))/distVecIJ.two_norm() < 0.001)
- {
- globalIdxK= this->problem().variables().index(*neighborPointer2);
- thirdCellPointer = neighborPointer2;
- foundK=true;
- }
-
- }
- }
- }
- if (foundIJ && foundK) break;
- }
-
- int isIndexJ = isIt->indexInOutside();
-
- // neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
- const GlobalPosition& globalPosInterface = isIt->geometry().center();
-
- Dune::FieldVector<Scalar,dimWorld> distVec = globalPosInterface - globalPos;
- Scalar lI= distVec*unitOuterNormal;
- distVec = globalPosNeighbor - globalPosInterface;
- Scalar lJ= distVec*unitOuterNormal;
- Scalar l=lI+lJ;
-
- FieldMatrix permeabilityI(0);
- FieldMatrix permeabilityJ(0);
- FieldMatrix permeabilityK(0);
-
- this->problem().spatialParameters().meanK(permeabilityI,
- this->problem().spatialParameters().intrinsicPermeability(*eIt));
- this->problem().spatialParameters().meanK(permeabilityJ,
- this->problem().spatialParameters().intrinsicPermeability(*neighborPointer));
- this->problem().spatialParameters().meanK(permeabilityK,
- this->problem().spatialParameters().intrinsicPermeability(*thirdCellPointer));
-
- // Calculate permeablity component normal to interface
- Scalar kI, kJ, kK, ng, kMean;//, gI, gJ, gK;
- Dune::FieldVector<Scalar, dim> permI(0);
- Dune::FieldVector<Scalar, dim> permJ(0);
- Dune::FieldVector<Scalar, dim> permK(0);
-
- permeabilityI.mv(unitOuterNormal, permI);
- permeabilityJ.mv(unitOuterNormal, permJ);
- permeabilityK.mv(unitOuterNormal, permK);
-
- // kI,kJ,kK = (n^T)Kn
- kI=unitOuterNormal*permI;
- kJ=unitOuterNormal*permJ;
- kK=unitOuterNormal*permK;
- kMean=kI*kJ*kK*l/(kJ*kK*lI+kI*(kJ+kK)/2*lJ);
-
- ng=this->gravity*unitOuterNormal;
-
- // find secondary variables in cell J and K
- Scalar pressJ = this->problem().variables().pressure()[globalIdxJ];
- Scalar pcJ = this->problem().variables().capillaryPressure(globalIdxJ);
- Scalar lambdaWJ = this->problem().variables().mobilityWetting(globalIdxJ);
- Scalar lambdaNWJ = this->problem().variables().mobilityNonwetting(globalIdxJ);
- Scalar densityWJ = this->problem().variables().densityWetting(globalIdxJ);
- Scalar densityNWJ = this->problem().variables().densityNonwetting(globalIdxJ);
- Scalar fractionalWJ = this->problem().variables().fracFlowFuncWetting(globalIdxJ);
- Scalar fractionalNWJ = this->problem().variables().fracFlowFuncNonwetting(globalIdxJ);
-
- Scalar pressK = this->problem().variables().pressure()[globalIdxK];
- Scalar pcK = this->problem().variables().capillaryPressure(globalIdxK);
- Scalar densityWK = this->problem().variables().densityWetting(globalIdxK);
- Scalar densityNWK = this->problem().variables().densityNonwetting(globalIdxK);
- Scalar fractionalWK = this->problem().variables().fracFlowFuncWetting(globalIdxK);
- Scalar fractionalNWK = this->problem().variables().fracFlowFuncNonwetting(globalIdxK);
-
- Scalar pressJK=(pressJ+pressK)/2;
- Scalar pcJK=(pcJ+pcK)/2;
-
-
- // calculate potential gradients
- // reuse potentials from fvpressure2padaptive
-
- Scalar potentialWIJ = this->problem().variables().potentialWetting(globalIdxI, isIndex);
- Scalar potentialNWIJ = this->problem().variables().potentialNonwetting(globalIdxI, isIndex);
- // We dont know the insideIndex of interface IK
- Scalar potentialWIK = potentialWIJ;
- Scalar potentialNWIK = potentialNWIJ;
- // preliminary potential. The "real" ones are found below
-
- // Comment: reversed weighting is plausible, too (swap lJ and lI)
- Scalar rhoMeanWIJ = (lJ*densityWI + lI*densityWJ)/l;
- Scalar rhoMeanNWIJ = (lJ*densityNWI + lI*densityNWJ)/l;
- Scalar rhoMeanWIK = (lJ*densityWI + lI*densityWK)/l;
- Scalar rhoMeanNWIK = (lJ*densityNWI + lI*densityNWK)/l;
-
- Scalar fMeanWIJ = (lJ*fractionalWI+lI*fractionalWJ)/l;
- Scalar fMeanNWIJ = (lJ*fractionalNWI+lI*fractionalNWJ)/l;
- Scalar fMeanWIK = (lJ*fractionalWI+lI*fractionalWK)/l;
- Scalar fMeanNWIK = (lJ*fractionalNWI+lI*fractionalNWK)/l;
-
- // Upwinding for finding the upwinding direction
- Scalar densityWIJ = (potentialWIJ> 0.) ? densityWI : densityWJ;
- Scalar densityNWIJ = (potentialNWIJ> 0.) ? densityNWI : densityNWJ;
- Scalar densityWIK = (potentialWIK> 0.) ? densityWI : densityWK;
- Scalar densityNWIK = (potentialNWIK> 0.) ? densityNWI : densityNWK;
-
- densityWIJ = (potentialWIJ == 0.) ? rhoMeanWIJ : densityWIJ;
- densityNWIJ = (potentialNWIJ == 0.) ? rhoMeanNWIJ : densityNWIJ;
- densityWIK = (potentialWIK == 0.) ? rhoMeanWIK : densityWIK;
- densityNWIK = (potentialNWIK == 0.) ? rhoMeanNWIK : densityNWIK;
-
- Scalar fractionalWIJ = (potentialWIJ> 0.) ? fractionalWI : fractionalWJ;
- Scalar fractionalNWIJ = (potentialNWIJ> 0.) ? fractionalNWI : fractionalNWJ;
- Scalar fractionalWIK = (potentialWIK> 0.) ? fractionalWI : fractionalWK;
- Scalar fractionalNWIK = (potentialNWIK> 0.) ? fractionalNWI : fractionalNWK;
-
- fractionalWIJ = (potentialWIJ == 0.) ? fMeanWIJ : fractionalWIJ;
- fractionalNWIJ = (potentialNWIJ == 0.) ? fMeanNWIJ : fractionalNWIJ;
- fractionalWIK = (potentialWIK == 0.) ? fMeanWIK : fractionalWIK;
- fractionalNWIK = (potentialNWIK == 0.) ? fMeanNWIK : fractionalNWIK;
-
- switch (this->pressureType)
- {
- case pw:
- {
- potentialWIJ = (pressI-pressJK)/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI+pcI-(pressJK+pcJK))/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-pressJK)/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI+pcI-(pressJK+pcJK))/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- case pn:
- {
- potentialWIJ = (pressI-pcI-(pressJK-pcJK))/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI-pressJK)/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-pcI-(pressJK-pcJK))/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI-pressJK)/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- case pglobal:
- {
- potentialWIJ = (pressI-fractionalNWIJ*pcI-(pressJK-fractionalNWIJ*pcJK))/l+
- (densityWIJ-lJ/l*(kI+kK)/kI*(densityWIK-densityWIJ)/2)*ng;
- potentialNWIJ = (pressI+fractionalWIJ*pcI-(pressJK+fractionalWIJ*pcJK))/l+
- (densityNWIJ-lJ/l*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2)*ng;
- potentialWIK = (pressI-fractionalNWIK*pcI-(pressJK-fractionalNWIK*pcJK))/l+
- (densityWIK-lJ/l*(kI+kK)/kI*(densityWIJ-densityWIK)/2)*ng;
- potentialNWIK = (pressI+fractionalWIK*pcI-(pressJK+fractionalWIK*pcJK))/l+
- (densityNWIK-lJ/l*(kI+kK)/kI*(densityNWIJ-densityNWIK)/2)*ng;
- break;
- }
- }
-
- //store potentials for further calculations (velocity, saturation, ...)
- // these quantities only have correct sign (needed for upwinding)
- // potentials are defined slightly different for adaptive scheme
- this->problem().variables().potentialWetting(globalIdxI, isIndex) = potentialWIJ;
- this->problem().variables().potentialNonwetting(globalIdxI, isIndex) = potentialNWIJ;
- this->problem().variables().potentialWetting(globalIdxJ, isIndexJ) = -potentialWIJ;
- this->problem().variables().potentialNonwetting(globalIdxJ, isIndexJ) = -potentialNWIJ;
-
-
- //do the upwinding of the mobility depending on the phase potentials
- Scalar lambdaWIJ = (potentialWIJ > 0.) ? lambdaWI : lambdaWJ;
- lambdaWIJ = (potentialWIJ == 0) ? 0.5 * (lambdaWI + lambdaWJ) : lambdaWIJ;
- Scalar lambdaNWIJ = (potentialNWIJ > 0.) ? lambdaNWI : lambdaNWJ;
- lambdaNWIJ = (potentialNWIJ == 0) ? 0.5 * (lambdaNWI + lambdaNWJ) : lambdaNWIJ;
-
- densityWIJ = (potentialWIJ> 0.) ? densityWI : densityWJ;
- densityNWIJ = (potentialNWIJ> 0.) ? densityNWI : densityNWJ;
- densityWIK = (potentialWIK> 0.) ? densityWI : densityWK;
- densityNWIK = (potentialNWIK> 0.) ? densityNWI : densityNWK;
-
- densityWIJ = (potentialWIJ == 0.) ? rhoMeanWIJ : densityWIJ;
- densityNWIJ = (potentialNWIJ == 0.) ? rhoMeanNWIJ : densityNWIJ;
- densityWIK = (potentialWIK == 0.) ? rhoMeanWIK : densityWIK;
- densityNWIK = (potentialNWIK == 0.) ? rhoMeanNWIK : densityNWIK;
-
- fractionalWIJ = (potentialWIJ> 0.) ? fractionalWI : fractionalWJ;
- fractionalNWIJ = (potentialNWIJ> 0.) ? fractionalNWI : fractionalNWJ;
- fractionalWIK = (potentialWIK> 0.) ? fractionalWI : fractionalWK;
- fractionalNWIK = (potentialNWIK> 0.) ? fractionalNWI : fractionalNWK;
-
- fractionalWIJ = (potentialWIJ == 0.) ? fMeanWIJ : fractionalWIJ;
- fractionalNWIJ = (potentialNWIJ == 0.) ? fMeanNWIJ : fractionalNWIJ;
- fractionalWIK = (potentialWIK == 0.) ? fMeanWIK : fractionalWIK;
- fractionalNWIK = (potentialNWIK == 0.) ? fMeanNWIK : fractionalNWIK;
-
- //calculate velocities and the gravity term
- Dune::FieldVector<Scalar,dimWorld> velocityW(unitOuterNormal);
- Dune::FieldVector<Scalar,dimWorld> velocityNW(unitOuterNormal);
- Dune::FieldVector<Scalar,dimWorld> gravityTermW(unitOuterNormal);
- Dune::FieldVector<Scalar,dimWorld> gravityTermNW(unitOuterNormal);
-
- gravityTermW *= lambdaWIJ*kMean*ng;
- gravityTermW *= densityWIJ-(lJ/l)*(kI+kK)/kI*(densityWIK-densityWIJ)/2;
- gravityTermNW *= lambdaNWIJ*kMean*ng;
- gravityTermNW *= densityNWIJ-(lJ/l)*(kI+kK)/kI*(densityNWIK-densityNWIJ)/2;
-
- switch (this->pressureType)
- {
- case pw:
- {
- velocityW *= lambdaWIJ*kMean*(pressI-pressJK)/l;
- velocityNW *= lambdaNWIJ*kMean*(pressI+pcI-(pressJK+pcJK))/l;
-
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- case pn:
- {
- velocityNW *= lambdaNWIJ*kMean*(pressI-pressJK)/l;
- velocityW *= lambdaWIJ*kMean*(pressI-pcI-(pressJK-pcJK))/l;
-
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- case pglobal:
- {
- velocityW *= lambdaWIJ*kMean*(pressI-fractionalNWIJ*pcI-(pressJK-fractionalNWIJ*pcJK))/l;
- velocityNW *= lambdaNWIJ*kMean*(pressI+fractionalWIJ*pcI-(pressJK+fractionalWIJ*pcJK))/l;
-
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- }
-
- switch (velocityType_)
- {
- case vw:
- {
- Dune::FieldVector<Scalar,dimWorld> vel(velocityW);
- vel*=0.5;//0.5 -> only one hanging node per face!
- this->problem().variables().velocity()[globalIdxI][isIndex] += vel;
- vel = velocityNW;
- vel*=0.5;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += vel;
- this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityW;
- this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityNW;
- break;
- }
- case vn:
- {
- Dune::FieldVector<Scalar,dimWorld> vel(velocityNW);
- vel*=0.5;
- this->problem().variables().velocity()[globalIdxI][isIndex] += (vel);
- vel = velocityW;
- vel*=0.5;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += (vel);
- this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityW;
- break;
- }
- case vt:
- {
- switch (this->pressureType)
- {
- case pw:
- {
- Dune::FieldVector<Scalar,dimWorld> vel(velocityW);
- vel+=velocityNW;
- this->problem().variables().velocity()[globalIdxI][isIndex] += (vel *=0.5);
- vel = velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += (vel*=0.5);
- this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityNW;
-
- break;
- }
- case pn:
- {
- Dune::FieldVector<Scalar,dimWorld> vel(velocityW);
- vel+=velocityNW;
- this->problem().variables().velocity()[globalIdxI][isIndex] += (vel *=0.5);
- vel = velocityW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] += (vel*=0.5);
- this->problem().variables().velocity()[globalIdxJ][isIndexJ] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxJ][isIndexJ] = velocityW;
- break;
- }
- }
- break;
- }
- }
- }
- }//end intersection with neighbor
+ velocityW *= lambdaWIJ * kMean * (cellDataI.globalPressure() - pressJK) / l;
+ velocityW += gravityTermW;
+ velocityW += gravityTermNW;
+ velocityNW = 0;
+ break;
+ }
+ }
- // handle boundary face
- if (isIt->boundary())
- {
- //get boundary type
- this->problem().boundaryTypes(bcType, *isIt);
- PrimaryVariables boundValues(0.0);
+ cellDataJ.fluxData().setVelocity(wPhaseIdx, isIndexJ, velocityW);
+ cellDataJ.fluxData().setVelocity(nPhaseIdx, isIndexJ, velocityNW);
+ cellDataJ.fluxData().setVelocityMarker(isIndexJ);
- if (bcType.isDirichlet(eqIdxPress))
- {
- this->problem().dirichlet(boundValues, *isIt);
-
- // center of face in global coordinates
- GlobalPosition globalPosFace = isIt->geometry().center();
-
- // distance vector between barycenters
- Dune::FieldVector<Scalar,dimWorld> distVec = globalPosFace - globalPos;
-
- // compute distance between cell centers
- Scalar dist = distVec * unitOuterNormal;
-// Scalar dist = distVec.two_norm();
-
- //permeability vector at boundary
- // compute vectorized permeabilities
- FieldMatrix meanPermeability(0);
-
- this->problem().spatialParameters().meanK(meanPermeability,
- this->problem().spatialParameters().intrinsicPermeability(*eIt));
-
- Dune::FieldVector<Scalar, dim> permeability(0);
- meanPermeability.mv(unitOuterNormal, permeability);
-
- Scalar satBound = 0;
- if (bcType.isDirichlet(eqIdxSat))
- {
- satBound = boundValues[saturationIdx];
- }
- else
- {
- satBound = this->problem().variables().saturation()[globalIdxI];
- }
-
- //determine phase saturations from primary saturation variable
- Scalar satW;
- //Scalar satNW;
- switch (this->saturationType)
- {
- case Sw:
- {
- satW = satBound;
- //satNW = 1-satBound;
- break;
- }
- case Sn:
- {
- satW = 1-satBound;
- //satNW = satBound;
- break;
- }
- default:
- {
- DUNE_THROW(Dune::RangeError, "saturation type not implemented");
- break;
- }
- }
- Scalar pressBound = boundValues[pressureIdx];
- Scalar pcBound = MaterialLaw::pC(this->problem().spatialParameters().materialLawParams(*eIt), satW);
-
- //determine phase pressures from primary pressure variable
- Scalar pressW=0;
- Scalar pressNW=0;
- switch (this->pressureType)
- {
- case pw:
- {
- pressW = pressBound;
- pressNW = pressBound + pcBound;
- break;
- }
- case pn:
- {
- pressW = pressBound - pcBound;
- pressNW = pressBound;
- break;
- }
- }
-
- //get temperature at current position
- Scalar temperature = this->problem().temperature(*eIt);
-
- Scalar densityWBound = 0;
- Scalar densityNWBound = 0;
- Scalar lambdaWBound = 0;
- Scalar lambdaNWBound = 0;
- if (this->compressibility)
- {
- FluidState fluidState;
- fluidState.update(satW, pressW, pressNW, temperature);
- densityWBound = FluidSystem::phaseDensity(wPhaseIdx, temperature, pressW, fluidState);
- densityNWBound = FluidSystem::phaseDensity(nPhaseIdx, temperature, pressNW, fluidState);
- Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx, temperature, pressW, fluidState);
- Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx, temperature, pressNW, fluidState);
- lambdaWBound = MaterialLaw::krw(this->problem().spatialParameters().materialLawParams(*eIt), satW) / viscosityWBound * densityWBound;
- lambdaNWBound = MaterialLaw::krn(this->problem().spatialParameters().materialLawParams(*eIt), satW) / viscosityNWBound * densityNWBound;
- }
- else
- {
- Scalar referencePressure = this->problem().referencePressure(*eIt);
- FluidState fluidState;
- fluidState.update(satW, referencePressure, referencePressure, temperature);
- densityWBound = FluidSystem::phaseDensity(wPhaseIdx, temperature, referencePressure, fluidState);
- densityNWBound = FluidSystem::phaseDensity(nPhaseIdx, temperature, referencePressure, fluidState);
- Scalar viscosityWBound = FluidSystem::phaseViscosity(wPhaseIdx, temperature, referencePressure, fluidState);
- Scalar viscosityNWBound = FluidSystem::phaseViscosity(nPhaseIdx, temperature, referencePressure, fluidState);
- lambdaWBound = MaterialLaw::krw(this->problem().spatialParameters().materialLawParams(*eIt), satW) / viscosityWBound;
- lambdaNWBound = MaterialLaw::krn(this->problem().spatialParameters().materialLawParams(*eIt), satW) / viscosityNWBound;
- }
-
- Scalar potentialW = 0;
- Scalar potentialNW = 0;
-
- potentialW = this->problem().variables().potentialWetting(globalIdxI, isIndex);
- potentialNW = this->problem().variables().potentialNonwetting(globalIdxI, isIndex);
-
- Scalar densityW = (potentialW> 0.) ? densityWI : densityWBound;
- Scalar densityNW = (potentialNW> 0.) ? densityNWI : densityNWBound;
-
- densityW = (potentialW == 0.) ? 0.5 * (densityWI + densityWBound) : densityW;
- densityNW = (potentialNW == 0.) ? 0.5 * (densityNWI + densityNWBound) : densityNW;
-
- //calculate potential gradient
- switch (this->pressureType)
- {
- case pw:
- {
- potentialW = (pressI - pressBound);
- potentialNW = (pressI + pcI - pressBound - pcBound);
- break;
- }
- case pn:
- {
- potentialW = (pressI - pcI - pressBound + pcBound);
- potentialNW = (pressI - pressBound);
- break;
- }
- case pglobal:
- {
- potentialW = (pressI - pressBound - fractionalNWI * (pcI - pcBound));
- potentialNW = (pressI - pressBound + fractionalWI * (pcI - pcBound));
- break;
- }
- }
-
- potentialW += densityW * (distVec * this->gravity);
- potentialNW += densityNW * (distVec * this->gravity);
-
- //store potential gradients for further calculations
- this->problem().variables().potentialWetting(globalIdxI, isIndex) = potentialW;
- this->problem().variables().potentialNonwetting(globalIdxI, isIndex) = potentialNW;
-
- //do the upwinding of the mobility depending on the phase potentials
- Scalar lambdaW = (potentialW > 0.) ? lambdaWI : lambdaWBound;
- lambdaW = (potentialW == 0) ? 0.5 * (lambdaWI + lambdaWBound) : lambdaW;
- Scalar lambdaNW = (potentialNW > 0.) ? lambdaNWI : lambdaNWBound;
- lambdaNW = (potentialNW == 0) ? 0.5 * (lambdaNWI + lambdaNWBound) : lambdaNW;
- densityW = (potentialW> 0.) ? densityWI : densityWBound;
- densityW = (potentialW == 0.) ? 0.5 * (densityWI + densityWBound) : densityW;
- densityNW = (potentialNW> 0.) ? densityNWI : densityNWBound;
- densityNW = (potentialNW == 0.) ? 0.5 * (densityNWI + densityNWBound) : densityNW;
-
- //calculate the gravity term
- Dune::FieldVector<Scalar,dimWorld> velocityW(permeability);
- Dune::FieldVector<Scalar,dimWorld> velocityNW(permeability);
- Dune::FieldVector<Scalar,dimWorld> gravityTermW(unitOuterNormal);
- Dune::FieldVector<Scalar,dimWorld> gravityTermNW(unitOuterNormal);
-
- gravityTermW *= (this->gravity*permeability)*(lambdaW * densityW);
- gravityTermNW *= (this->gravity*permeability)*(lambdaNW * densityNW);
-
- //calculate velocity depending on the pressure used -> use pc = pn - pw
- switch (this->pressureType)
- {
- case pw:
- {
- velocityW *= lambdaW * (pressI - pressBound)/dist;
- velocityNW *= lambdaNW * (pressI - pressBound)/dist + 0.5 * (lambdaNWI + lambdaNWBound) * (pcI - pcBound) / dist;
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- case pn:
- {
- velocityW *= lambdaW * (pressI - pressBound)/dist - 0.5 * (lambdaWI + lambdaWBound) * (pcI - pcBound) / dist;
- velocityNW *= lambdaNW * (pressI - pressBound) / dist;
- velocityW += gravityTermW;
- velocityNW += gravityTermNW;
- break;
- }
- case pglobal:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = permeability;
- this->problem().variables().velocity()[globalIdxI][isIndex] *= (lambdaW + lambdaNW)* (pressI - pressBound )/ dist;
- this->problem().variables().velocity()[globalIdxI][isIndex] += gravityTermW;
- this->problem().variables().velocity()[globalIdxI][isIndex] += gravityTermNW;
- break;
- }
- }
-
- //store velocities
- switch (velocityType_)
- {
- case vw:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
- break;
- }
- case vn:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
- break;
- }
- case vt:
- {
- switch (this->pressureType)
- {
- case pw:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
- break;
- }
- case pn:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
- break;
- }
- }
- break;
- }
- }
- }//end dirichlet boundary
-
- else if (bcType.isNeumann(eqIdxPress))
- {
- this->problem().neumann(boundValues, *isIt);
-
- Dune::FieldVector<Scalar,dimWorld> velocityW(unitOuterNormal);
- Dune::FieldVector<Scalar,dimWorld> velocityNW(unitOuterNormal);
-
- velocityW *= boundValues[wPhaseIdx];
- velocityNW *= boundValues[nPhaseIdx];
-
- if (!this->compressibility)
- {
- velocityW /= densityWI;
- velocityNW /= densityNWI;
- }
-
- switch (velocityType_)
- {
- case vw:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
- break;
- }
- case vn:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
- break;
- }
- case vt:
- {
- switch (this->pressureType)
- {
- case pw:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityNW;
- break;
- }
- case pn:
- {
- this->problem().variables().velocity()[globalIdxI][isIndex] = velocityW + velocityNW;
- this->problem().variables().velocitySecondPhase()[globalIdxI][isIndex] = velocityW;
- break;
- }
- }
- break;
- }
- }
- }//end neumann boundary
- else
- {
- DUNE_THROW(Dune::NotImplemented, "No valid boundary condition type defined for pressure equation!");
- }
- }//end boundary
- }// end all intersections
- }// end grid traversal
-// printvector(std::cout, this->problem().variables().velocity(), "velocity", "row", 4, 1, 3);
+ velocityW *= 0.5;
+ velocityNW *= 0.5;
+ cellDataI.fluxData().setVelocity(wPhaseIdx, isIndexI, velocityW);
+ cellDataI.fluxData().setVelocity(nPhaseIdx, isIndexI, velocityNW);
+ cellDataI.fluxData().setVelocityMarker(isIndexI);
+ }
return;
}
}