diff --git a/dumux/decoupled/2p2c/2p2cproperties.hh b/dumux/decoupled/2p2c/2p2cproperties.hh
index ff40676991d13524bb966e9a5565ae0f937e7124..851ed2018e0d23d74f83b89e12abe5f014269f17 100644
--- a/dumux/decoupled/2p2c/2p2cproperties.hh
+++ b/dumux/decoupled/2p2c/2p2cproperties.hh
@@ -75,7 +75,7 @@ NEW_PROP_TAG( VelocityFormulation); //!< The formulation of the model
NEW_PROP_TAG( EnableCompressibility);// ! Returns whether compressibility is allowed
NEW_PROP_TAG( EnableCapillarity); //!< Returns whether capillarity is regarded
NEW_PROP_TAG( BoundaryMobility );
-NEW_PROP_TAG( NumDensityTransport );
+NEW_PROP_TAG( RestrictFluxInTransport ); //!< Restrict flux if direction reverses after pressure equation
NEW_PROP_TAG( ErrorTermFactor );
NEW_PROP_TAG( ErrorTermLowerBound );
NEW_PROP_TAG( ErrorTermUpperBound );
@@ -149,8 +149,8 @@ SET_BOOL_PROP(DecoupledTwoPTwoC, EnableCompressibility, true);
//! Faces are regarded from both sides
SET_BOOL_PROP(DecoupledTwoPTwoC, VisitFacesOnlyOnce, false);
SET_BOOL_PROP(DecoupledTwoPTwoC, EnableCapillarity, false);
-//! Disable transport of numerical density (e.g. inclusion of error in transport)
-SET_BOOL_PROP(DecoupledTwoPTwoC, NumDensityTransport, false);
+//!< Restrict (no upwind) flux in transport step if direction reverses after pressure equation
+SET_BOOL_PROP(DecoupledTwoPTwoC, RestrictFluxInTransport, true);
SET_PROP(DecoupledTwoPTwoC, BoundaryMobility)
{
diff --git a/dumux/decoupled/2p2c/cellData2p2c.hh b/dumux/decoupled/2p2c/cellData2p2c.hh
index 2e069cb336910ed35084c8439919a4486519a580..f9f721dce19153abd7f3a403287b022d5a5ebaf6 100644
--- a/dumux/decoupled/2p2c/cellData2p2c.hh
+++ b/dumux/decoupled/2p2c/cellData2p2c.hh
@@ -24,7 +24,7 @@
#include "2p2cproperties.hh"
#include <dumux/decoupled/2p2c/dec2p2cfluidstate.hh>
-//#include "fluxData2p2c.hh"
+#include "fluxData2p2c.hh"
/**
* @file
@@ -51,7 +51,7 @@ class CellData2P2C
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-// typedef FluxData2P2C<TypeTag> FluxData;
+ typedef FluxData2P2C<TypeTag> FluxData;
typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
enum
@@ -76,7 +76,6 @@ protected:
Scalar mobility_[numPhases];
//numerical quantities
- Scalar numericalDensity_[numPhases];
Scalar volumeError_;
Scalar errorCorrection_;
Scalar dv_dp_;
@@ -87,7 +86,7 @@ protected:
Scalar perimeter_;
FluidState* fluidState_;
-// FluxData fluxData_;
+ FluxData fluxData_;
public:
//! Constructor for a local CellData object
@@ -96,8 +95,6 @@ public:
{
for (int i = 0; i < numPhases;i++)
{
- mobility_[i] = 0.0;
- numericalDensity_[i] = 0.0;
mobility_[i] = 0.0;
dv_[i] = 0.0;
}
@@ -109,15 +106,15 @@ public:
perimeter_ = 0.;
}
-// FluxData& setFluxData()
-// {
-// return fluxData_;
-// }
-//
-// const FluxData& fluxData() const
-// {
-// return fluxData_;
-// }
+ FluxData& fluxData()
+ {
+ return fluxData_;
+ }
+
+ const FluxData& fluxData() const
+ {
+ return fluxData_;
+ }
/*! \name Acess to primary variables */
@@ -189,22 +186,6 @@ public:
mobility_[phaseIdx]=value;
}
- //! Return numerical density \f$\mathrm{[kg/m^3]}\f$.
- /** Returns the real fluid density if the whole pore volume
- * was filled. This quantity equals the real density scaled with
- * the volume error.
- * @param phaseIdx index of the Phase
- */
- Scalar& numericalDensity(int phaseIdx)
- {
- return numericalDensity_[phaseIdx];
- }
- //! Return numerical density
- const Scalar& numericalDensity(int phaseIdx) const
- {
- return numericalDensity_[phaseIdx];
- }
-
//! Return the volume error [-].
/** This quantity stands for the deviation of real fluid volume
* to available pore space.
@@ -377,6 +358,23 @@ public:
volumeDerivativesAvailable_ = false;
}
+ //! Indicates if current cell is the upwind cell for a given interface
+ /* @param indexInInside Local face index seen from current cell
+ * @param phaseIdx The index of the phase
+ */
+ const bool& isUpwindCell(int indexInInside, int phaseIdx) const
+ {
+ return fluxData_.isUpwindCell(indexInInside, phaseIdx);
+ }
+ //! Specifies if current cell is the upwind cell for a given interface
+ /* @param indexInInside Local face index seen from current cell
+ * @param phaseIdx The index of the phase
+ * @param value Value: true (->outflow) or false (-> inflow)
+ */
+ void setUpwindCell(int indexInInside, int phaseIdx, bool value)
+ {
+ fluxData_.setUpwindCell(indexInInside, phaseIdx, value);
+ }
};
}
diff --git a/dumux/decoupled/2p2c/fluxData2p2c.hh b/dumux/decoupled/2p2c/fluxData2p2c.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f2b893d4d04484f9c9d4d359145c4a0a8e30a9bf
--- /dev/null
+++ b/dumux/decoupled/2p2c/fluxData2p2c.hh
@@ -0,0 +1,149 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * Copyright (C) 2011 by Markus Wolff *
+ * Institute for Modelling Hydraulic and Environmental Systems *
+ * University of Stuttgart, Germany *
+ * email: <givenname>.<name>@iws.uni-stuttgart.de *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 2 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+#ifndef DUMUX_FLUXDATA2P2C_HH
+#define DUMUX_FLUXDATA2P2C_HH
+
+/**
+ * @file
+ * @brief Class including the variables and data of discretized data of the constitutive relations
+ * @author Markus Wolff
+ */
+
+namespace Dumux
+{
+/*!
+ * \ingroup IMPES
+ */
+//! Class including the variables and data of discretized data of the constitutive relations.
+/*! The variables of two-phase flow, which are one pressure and one saturation are stored in this class.
+ * Additionally, a velocity needed in the transport part of the decoupled two-phase flow is stored, as well as discretized data of constitutive relationships like
+ * mobilities, fractional flow functions and capillary pressure. Thus, they have to be callculated just once in every time step or every iteration step.
+ *
+ * @tparam TypeTag The Type Tag
+ 1*/
+template<class TypeTag>
+class FluxData2P2C
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Indices)) Indices;
+
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
+ };
+
+ enum
+ {
+ numEquations = GET_PROP_VALUE(TypeTag, PTAG(NumEq))
+ };
+
+ typename Dune::BlockVector<typename Dune::FieldVector<bool, numEquations>> isUpwindCell_;
+
+public:
+
+ //! Constructs a VariableClass object
+ /**
+ * @param gridView a DUNE gridview object corresponding to diffusion and transport equation
+ */
+
+ FluxData2P2C()
+ {
+ isUpwindCell_.resize(2 * dim);
+ for (int i = 0; i<2*dim; i++)
+ {
+ isUpwindCell_[i] = false;
+ }
+ }
+ void resize(int size)
+ {
+ isUpwindCell_.resize(size);
+ }
+ /** functions returning upwind information **/
+ const bool& isUpwindCell(int indexInInside, int equationIdx) const
+ {
+ return isUpwindCell_[indexInInside][equationIdx];
+ }
+
+ void setUpwindCell(int indexInInside, int equationIdx, bool value)
+ {
+ isUpwindCell_[indexInInside][equationIdx] = value;
+ }
+
+};
+}
+#endif
+
+
+/*** in transport module:
+ * // upwind mobility
+ double lambdaW, lambdaNW;
+ if (potentialW >= 0.)
+ {
+ lambdaW = cellDataI.mobility(wPhaseIdx);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiWEqIdx, true);
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, false);
+ }
+ else
+ {
+ lambdaW = cellDataJ.mobility(wPhaseIdx);
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiWEqIdx, true);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiWEqIdx, false);
+ }
+
+ if (potentialNW >= 0.)
+ {
+ lambdaNW = cellDataI.mobility(nPhaseIdx);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiNEqIdx, true);
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, false);
+ }
+ else
+ {
+ lambdaW = cellDataJ.mobility(nPhaseIdx);
+ cellDataJ.setUpwindCell(intersection.indexInOutside(), contiNEqIdx, true);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiNEqIdx, false);
+ }
+
+
+ in cellData:
+
+ //! Acess to flux data
+ bool& isUpwindCell(int indexInInside, int equationIdx) const
+ {
+ return fluxData_.isUpwindCell(indexInInside, equationIdx);
+ }
+
+ void setUpwindCell(int indexInInside, int equationIdx, bool value)
+ {
+ fluxData_.setUpwindCell(indexInInside, equationIdx, value);
+ }
+ */
+
diff --git a/dumux/decoupled/2p2c/fvpressure2p2c.hh b/dumux/decoupled/2p2c/fvpressure2p2c.hh
index 06adf310c2faed323e39dd928d2e0a8c64f8e930..da0330756c9498627d66d4a2906834bfce15948e 100644
--- a/dumux/decoupled/2p2c/fvpressure2p2c.hh
+++ b/dumux/decoupled/2p2c/fvpressure2p2c.hh
@@ -71,6 +71,8 @@ namespace Dumux
template<class TypeTag> class FVPressure2P2C
: public FVPressureCompositional<TypeTag>
{
+ //the model implementation
+ typedef typename GET_PROP_TYPE(TypeTag, PressureModel) Implementation;
typedef FVPressure<TypeTag> BaseType;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
@@ -188,6 +190,14 @@ protected:
Scalar ErrorTermLowerBound_; //!< Handling of error term: lower bound for error dampening
Scalar ErrorTermUpperBound_; //!< Handling of error term: upper bound for error dampening
static constexpr int pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
+private:
+ //! Returns the implementation of the problem (i.e. static polymorphism)
+ Implementation &asImp_()
+ { return *static_cast<Implementation *>(this);}
+
+ //! \copydoc Dumux::IMPETProblem::asImp_()
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation *>(this);}
};
@@ -313,6 +323,8 @@ void FVPressure2P2C<TypeTag>::getStorage(Dune::FieldVector<Scalar, 2>& storageEn
fac * (1 + x_mi - hifac*x_mi/(1-x_mi) + (hifac/(1-x_mi)-1)*erri/maxError)
* cellDataI.volumeError() * volume;
}
+ else
+ problem().variables().cellData(globalIdxI).errorCorrection() = 0.;
return;
}
@@ -442,7 +454,6 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
Scalar graddv_dC2 = (cellDataJ.dv(nPhaseIdx)
- cellDataI.dv(nPhaseIdx)) / dist;
-
// potentialW = problem().variables().potentialWetting(globalIdxI, isIndex);
// potentialNW = problem().variables().potentialNonwetting(globalIdxI, isIndex);
//
@@ -455,8 +466,8 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
Scalar densityW = rhoMeanW;
Scalar densityNW = rhoMeanNW;
- potentialW = (unitOuterNormal * unitDistVec) *(cellDataI.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx))/dist;
- potentialNW = (unitOuterNormal * unitDistVec) *(cellDataI.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx))/dist;
+ potentialW = (cellDataI.pressure(wPhaseIdx) - cellDataJ.pressure(wPhaseIdx))/dist;
+ potentialNW = (cellDataI.pressure(nPhaseIdx) - cellDataJ.pressure(nPhaseIdx))/dist;
potentialW += densityW * (unitDistVec * gravity_);
potentialNW += densityNW * (unitDistVec * gravity_);
@@ -510,10 +521,10 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
}
//calculate current matrix entry
- entries[0] = faceArea * (lambdaW * dV_w + lambdaN * dV_n);
+ entries[matrix] = faceArea * (lambdaW * dV_w + lambdaN * dV_n)* (unitOuterNormal * unitDistVec);
if(enableVolumeIntegral)
- entries[0] -= volume * faceArea / perimeter * (lambdaW * gV_w + lambdaN * gV_n); // = boundary integral - area integral
- entries[0] *= fabs((permeability*unitOuterNormal)/(dist));
+ entries[matrix] -= volume * faceArea / perimeter * (lambdaW * gV_w + lambdaN * gV_n); // = boundary integral - area integral
+ entries[matrix] *= fabs((permeability*unitDistVec)/(dist));
//calculate right hand side
entries[rhs] = faceArea * (unitOuterNormal * unitDistVec) * (densityW * lambdaW * dV_w + densityNW * lambdaN * dV_n);
@@ -531,7 +542,7 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
case pw:
{
//add capillary pressure term to right hand side
- entries[rhs] += lambdaN * dV_n * (permeability * pCGradient) * faceArea;
+ entries[rhs] += lambdaN * dV_n * (permeability * pCGradient) * faceArea* (unitOuterNormal * unitDistVec);
if(enableVolumeIntegral)
entries[rhs]-= lambdaN * gV_n * (permeability * pCGradient) * volume * faceArea / perimeter;
break;
@@ -539,7 +550,7 @@ void FVPressure2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries,
case pn:
{
//add capillary pressure term to right hand side
- entries[rhs] -= lambdaW * dV_w * (permeability * pCGradient) * faceArea;
+ entries[rhs] -= lambdaW * dV_w * (permeability * pCGradient) * faceArea* (unitOuterNormal * unitDistVec);
if(enableVolumeIntegral)
entries[rhs]+= lambdaW * gV_w * (permeability * pCGradient) * volume * faceArea / perimeter;
break;
@@ -826,7 +837,7 @@ void FVPressure2P2C<TypeTag>::updateMaterialLaws()
CellData& cellData = problem().variables().cellData(globalIdx);
- this->updateMaterialLawsInElement(*eIt);
+ asImp_().updateMaterialLawsInElement(*eIt);
maxError = std::max(maxError, fabs(cellData.volumeError()));
}
@@ -973,8 +984,8 @@ void FVPressure2P2C<TypeTag>::updateMaterialLawsInElement(const Element& element
// determine volume mismatch between actual fluid volume and pore volume
Scalar sumConc = (cellData.totalConcentration(wCompIdx)
+ cellData.totalConcentration(nCompIdx));
- Scalar massw = cellData.numericalDensity(wPhaseIdx) = sumConc * fluidState.phaseMassFraction(wPhaseIdx);
- Scalar massn = cellData.numericalDensity(nPhaseIdx) = sumConc * fluidState.phaseMassFraction(nPhaseIdx);
+ Scalar massw = sumConc * fluidState.phaseMassFraction(wPhaseIdx);
+ Scalar massn = sumConc * fluidState.phaseMassFraction(nPhaseIdx);
if ((cellData.density(wPhaseIdx)*cellData.density(nPhaseIdx)) == 0)
DUNE_THROW(Dune::MathError, "Decoupled2p2c::postProcessUpdate: try to divide by 0 density");
diff --git a/dumux/decoupled/2p2c/fvtransport2p2c.hh b/dumux/decoupled/2p2c/fvtransport2p2c.hh
index a36c9de80428076da92369000b3282d2749f62e0..e74bf048b9b5eba3c42813f54fd9f57c6dc02b80 100644
--- a/dumux/decoupled/2p2c/fvtransport2p2c.hh
+++ b/dumux/decoupled/2p2c/fvtransport2p2c.hh
@@ -109,9 +109,11 @@ public:
void updateTransportedQuantity(TransportSolutionType& updateVector);
+ // Function which calculates the flux update
void getFlux(Dune::FieldVector<Scalar, 2>&, Dune::FieldVector<Scalar, 2>&,
- const Intersection&, const CellData&);
+ const Intersection&, CellData&);
+ // Function which calculates the boundary flux update
void getFluxOnBoundary(Dune::FieldVector<Scalar, 2>&, Dune::FieldVector<Scalar, 2>&,
const Intersection&, const CellData&);
@@ -192,6 +194,7 @@ public:
totalConcentration_[wCompIdx].resize(problem.gridView().size(0));
totalConcentration_[nCompIdx].resize(problem.gridView().size(0));
+ restrictFluxInTransport_ = GET_PARAM(TypeTag,bool, RestrictFluxInTransport);
}
virtual ~FVTransport2P2C()
@@ -201,8 +204,10 @@ public:
protected:
TransportSolutionType totalConcentration_;
Problem& problem_;
+ bool impet_;
static const int pressureType = GET_PROP_VALUE(TypeTag, PressureFormulation); //!< gives kind of pressure used (\f$ 0 = p_w \f$, \f$ 1 = p_n \f$, \f$ 2 = p_{global} \f$)
+ bool restrictFluxInTransport_;
bool switchNormals;
};
@@ -230,6 +235,8 @@ void FVTransport2P2C<TypeTag>::update(const Scalar t, Scalar& dt,
{
// initialize dt very large
dt = 1E100;
+ // store if we do update Estimate for flux functions
+ impet_ = impet;
// resize update vector and set to zero
updateVec.resize(GET_PROP_VALUE(TypeTag, NumComponents));
@@ -325,7 +332,8 @@ void FVTransport2P2C<TypeTag>::updateTransportedQuantity(TransportSolutionType&
\sum_{\alpha} \varrho_{\alpha} \lambda_{\alpha} \sum_{\kappa} X^{\kappa}_{\alpha}
\left( \frac{p_{\alpha,j}^t - p^{t}_{\alpha,i}}{\Delta x} + \varrho_{\alpha} \mathbf{g}\right) \f]
* Here, \f$ \mathbf{u} \f$ is the normalized vector connecting the cell centers, and \f$ \mathbf{n}_{\gamma_{ij}} \f$
- * represents the normal of the face \f$ \gamma{ij} \f$.
+ * represents the normal of the face \f$ \gamma{ij} \f$. Due to the nature of the Primay Variable, the (volume-)specific
+ * total mass concentration, this represents a mass flux per cell volume.
* \param fluxEntries The flux entries, mass influx from cell \f$j\f$ to \f$i\f$.
* \param timestepFlux flow velocities for timestep estimation
* \param intersection The intersection
@@ -335,8 +343,10 @@ template<class TypeTag>
void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries,
Dune::FieldVector<Scalar, 2>& timestepFlux,
const Intersection& intersection,
- const CellData& cellDataI)
+ CellData& cellDataI)
{
+ fluxEntries = 0.;
+ timestepFlux = 0.;
// cell information
ElementPointer elementI= intersection.inside();
int globalIdxI = problem().variables().index(*elementI);
@@ -360,17 +370,9 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
, 1e-2);
Scalar densityWI (0.), densityNWI(0.);
- if (GET_PROP_VALUE(TypeTag, NumDensityTransport))
- {
- densityWI= cellDataI.numericalDensity(wPhaseIdx);
- densityNWI = cellDataI.numericalDensity(nPhaseIdx);
+ densityWI= cellDataI.density(wPhaseIdx);
+ densityNWI = cellDataI.density(nPhaseIdx);
- }
- else
- {
- densityWI= cellDataI.density(wPhaseIdx);
- densityNWI = cellDataI.density(nPhaseIdx);
- }
// face properties
GlobalPosition unitOuterNormal = intersection.centerUnitOuterNormal();
if (switchNormals)
@@ -401,16 +403,8 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
// phase densities in neighbor
Scalar densityWJ (0.), densityNWJ(0.);
- if (GET_PROP_VALUE(TypeTag, NumDensityTransport))
- {
- densityWJ = cellDataJ.numericalDensity(wPhaseIdx);
- densityNWJ = cellDataJ.numericalDensity(nPhaseIdx);
- }
- else
- {
- densityWJ = cellDataJ.density(wPhaseIdx);
- densityNWJ = cellDataJ.density(nPhaseIdx);
- }
+ densityWJ = cellDataJ.density(wPhaseIdx);
+ densityNWJ = cellDataJ.density(nPhaseIdx);
// average phase densities with central weighting
double densityW_mean = (densityWI + densityWJ) * 0.5;
@@ -440,7 +434,7 @@ void FVTransport2P2C<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& fluxEntries
{
potentialW = (K * unitOuterNormal) * (pressI - pressJ - pcI + pcJ) / (dist);
potentialNW = (K * unitOuterNormal) * (pressI - pressJ) / (dist);
-break;
+ break;
}
}
// add gravity term
@@ -448,16 +442,110 @@ break;
potentialW += (K * gravity_) * (unitOuterNormal * unitDistVec) * densityW_mean;
// upwind mobility
- double lambdaW, lambdaNW;
- if (potentialW >= 0.)
- lambdaW = cellDataI.mobility(wPhaseIdx);
- else
- lambdaW = cellDataJ.mobility(wPhaseIdx);
+ double lambdaW(0.), lambdaNW(0.);
+ if(!impet_ or !restrictFluxInTransport_) // perform a simple uwpind scheme
+ {
+ if (potentialW >= 0.)
+ {
+ lambdaW = cellDataI.mobility(wPhaseIdx);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiWEqIdx, true);
- if (potentialNW >= 0.)
- lambdaNW = cellDataI.mobility(nPhaseIdx);
- else
- lambdaNW = cellDataJ.mobility(nPhaseIdx);
+ }
+ else
+ {
+ lambdaW = cellDataJ.mobility(wPhaseIdx);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiWEqIdx, false);
+ }
+
+ if (potentialNW >= 0.)
+ {
+ lambdaNW = cellDataI.mobility(nPhaseIdx);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiNEqIdx, true);
+ }
+ else
+ {
+ lambdaNW = cellDataJ.mobility(nPhaseIdx);
+ cellDataI.setUpwindCell(intersection.indexInInside(), contiNEqIdx, false);
+ }
+ }
+ else // if new potentials indicate same flow direction as in P.E., perform upwind
+ {
+ if (potentialW >= 0. && cellDataI.isUpwindCell(intersection.indexInInside(), contiWEqIdx))
+ lambdaW = cellDataI.mobility(wPhaseIdx);
+ else if (potentialW < 0. && !cellDataI.isUpwindCell(intersection.indexInInside(), contiWEqIdx))
+ lambdaW = cellDataJ.mobility(wPhaseIdx);
+ else // potential of w-phase does not coincide with that of P.E.
+ {
+ bool isUpwindCell = cellDataI.isUpwindCell(intersection.indexInInside(), contiWEqIdx);
+ //check if harmonic weithing is necessary
+ if (!isUpwindCell && !(cellDataI.mobility(wPhaseIdx) != 0. && cellDataJ.mobility(wPhaseIdx) == 0.)) // check if outflow induce neglected phase flux
+ lambdaW = cellDataI.mobility(wPhaseIdx);
+ else if (isUpwindCell && !(cellDataJ.mobility(wPhaseIdx) != 0. && cellDataI.mobility(wPhaseIdx) == 0.)) // check if inflow induce neglected phase flux
+ lambdaW = cellDataJ.mobility(wPhaseIdx);
+ else
+ {
+ //a) perform harmonic averageing
+ fluxEntries[wCompIdx] -= potentialW * faceArea / volume
+ * harmonicMean(cellDataI.massFraction(wPhaseIdx, wCompIdx) * cellDataI.mobility(wPhaseIdx) * cellDataI.density(wPhaseIdx),
+ cellDataJ.massFraction(wPhaseIdx, wCompIdx) * cellDataJ.mobility(wPhaseIdx) * cellDataJ.density(wPhaseIdx));
+ fluxEntries[nCompIdx] -= potentialW * faceArea / volume
+ * harmonicMean(cellDataI.massFraction(wPhaseIdx, nCompIdx) * cellDataI.mobility(wPhaseIdx) * cellDataI.density(wPhaseIdx),
+ cellDataJ.massFraction(wPhaseIdx, nCompIdx) * cellDataJ.mobility(wPhaseIdx) * cellDataJ.density(wPhaseIdx));
+ // b) timestep control
+ // for timestep control : influx
+ timestepFlux[0] += std::max(0.,
+ - potentialW * faceArea / volume * harmonicMean(cellDataI.density(wPhaseIdx),cellDataJ.density(wPhaseIdx)));
+ // outflux
+ timestepFlux[1] += std::max(0.,
+ potentialW * faceArea / volume * harmonicMean(cellDataI.density(wPhaseIdx),cellDataJ.density(wPhaseIdx)));
+
+ //c) stop further standard calculations
+ potentialW = 0;
+
+ //d) output (only for one side)
+ if(potentialW >= 0.)
+ Dune::dinfo << "harmonicMean flux of phase w used from cell" << globalIdxI<< " into " << globalIdxJ <<"\n";
+ }
+ }
+
+ if (potentialNW >= 0. && cellDataI.isUpwindCell(intersection.indexInInside(), contiNEqIdx))
+ lambdaNW = cellDataI.mobility(nPhaseIdx);
+ else if (potentialNW < 0. && !cellDataI.isUpwindCell(intersection.indexInInside(), contiNEqIdx))
+ lambdaNW = cellDataJ.mobility(nPhaseIdx);
+ else // potential of n-phase does not coincide with that of P.E.
+ {
+ bool isUpwindCell = cellDataI.isUpwindCell(intersection.indexInInside(), contiNEqIdx);
+ //check if harmonic weithing is necessary
+ if (!isUpwindCell && !(cellDataI.mobility(nPhaseIdx) != 0. && cellDataJ.mobility(nPhaseIdx) == 0.)) // check if outflow induce neglected phase flux
+ lambdaNW = cellDataI.mobility(nPhaseIdx);
+ else if (isUpwindCell && !(cellDataJ.mobility(nPhaseIdx) != 0. && cellDataI.mobility(nPhaseIdx) == 0.)) // check if inflow induce neglected phase flux
+ lambdaNW = cellDataJ.mobility(nPhaseIdx);
+ else
+ {
+ //a) perform harmonic averageing
+ fluxEntries[wCompIdx] -= potentialNW * faceArea / volume
+ * harmonicMean(cellDataI.massFraction(nPhaseIdx, wCompIdx) * cellDataI.mobility(nPhaseIdx) * cellDataI.density(nPhaseIdx),
+ cellDataJ.massFraction(nPhaseIdx, wCompIdx) * cellDataJ.mobility(nPhaseIdx) * cellDataJ.density(nPhaseIdx));
+ fluxEntries[nCompIdx] -= potentialNW * faceArea / volume
+ * harmonicMean(cellDataI.massFraction(nPhaseIdx, nCompIdx) * cellDataI.mobility(nPhaseIdx) * cellDataI.density(nPhaseIdx),
+ cellDataJ.massFraction(nPhaseIdx, nCompIdx) * cellDataJ.mobility(nPhaseIdx) * cellDataJ.density(nPhaseIdx));
+ // b) timestep control
+ // for timestep control : influx
+ timestepFlux[0] += std::max(0.,
+ - potentialNW * faceArea / volume * harmonicMean(cellDataI.density(nPhaseIdx),cellDataJ.density(nPhaseIdx)));
+ // outflux
+ timestepFlux[1] += std::max(0.,
+ potentialNW * faceArea / volume * harmonicMean(cellDataI.density(nPhaseIdx),cellDataJ.density(nPhaseIdx)));
+
+ //c) stop further standard calculations
+ potentialNW = 0;
+
+ //d) output (only for one side)
+ if(potentialNW >= 0.)
+ Dune::dinfo << "harmonicMean flux of phase n used from cell" << globalIdxI<< " into " << globalIdxJ <<"\n";
+ }
+ }
+ }
// calculate and standardized velocity
double velocityJIw = std::max((-lambdaW * potentialW) * faceArea / volume, 0.0);
@@ -465,25 +553,26 @@ break;
double velocityJIn = std::max((-lambdaNW * potentialNW) * faceArea / volume, 0.0);
double velocityIJn = std::max(( lambdaNW * potentialNW) * faceArea / volume, 0.0);
- // for timestep control
- timestepFlux[0] = velocityJIw + velocityJIn;
+ // for timestep control : influx
+ timestepFlux[0] += velocityJIw + velocityJIn;
double foutw = velocityIJw/SwmobI;
double foutn = velocityIJn/SnmobI;
if (std::isnan(foutw) || std::isinf(foutw) || foutw < 0) foutw = 0;
if (std::isnan(foutn) || std::isinf(foutn) || foutn < 0) foutn = 0;
- timestepFlux[1] = foutw + foutn;
+ timestepFlux[1] += foutw + foutn;
- fluxEntries[wCompIdx] =
+ fluxEntries[wCompIdx] +=
velocityJIw * cellDataJ.massFraction(wPhaseIdx, wCompIdx) * densityWJ
- velocityIJw * cellDataI.massFraction(wPhaseIdx, wCompIdx) * densityWI
+ velocityJIn * cellDataJ.massFraction(nPhaseIdx, wCompIdx) * densityNWJ
- velocityIJn * cellDataI.massFraction(nPhaseIdx, wCompIdx) * densityNWI;
- fluxEntries[nCompIdx] =
+ fluxEntries[nCompIdx] +=
velocityJIw * cellDataJ.massFraction(wPhaseIdx, nCompIdx) * densityWJ
- velocityIJw * cellDataI.massFraction(wPhaseIdx, nCompIdx) * densityWI
+ velocityJIn * cellDataJ.massFraction(nPhaseIdx, nCompIdx) * densityNWJ
- velocityIJn * cellDataI.massFraction(nPhaseIdx, nCompIdx) * densityNWI;
+
return;
}
//! Get flux on Boundary
@@ -528,17 +617,8 @@ void FVTransport2P2C<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& f
, 1e-2);
Scalar densityWI (0.), densityNWI(0.);
- if (GET_PROP_VALUE(TypeTag, NumDensityTransport))
- {
- densityWI= cellDataI.numericalDensity(wPhaseIdx);
- densityNWI = cellDataI.numericalDensity(nPhaseIdx);
-
- }
- else
- {
- densityWI= cellDataI.density(wPhaseIdx);
- densityNWI = cellDataI.density(nPhaseIdx);
- }
+ densityWI= cellDataI.density(wPhaseIdx);
+ densityNWI = cellDataI.density(nPhaseIdx);
// face properties
GlobalPosition unitOuterNormal = intersection.centerUnitOuterNormal();