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

changed finite volume discretization of decoupled 2p models to new

structure



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7412 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 38216873
......@@ -37,7 +37,7 @@
//Dumux-includes
#include <dumux/decoupled/common/impetproperties.hh>
#include <dumux/decoupled/2p/transport/transportproperties.hh>
#include <dumux/decoupled/2p/transport/transportproperties2p.hh>
#include "2pindices.hh"
namespace Dumux
......@@ -69,7 +69,7 @@ namespace Properties
//////////////////////////////////////////////////////////////////
//! The type tag for the two-phase problems
NEW_TYPE_TAG(DecoupledTwoP, INHERITS_FROM(IMPET, Transport))
NEW_TYPE_TAG(DecoupledTwoP, INHERITS_FROM(IMPET, TransportTwoP))
;
//////////////////////////////////////////////////////////////////
......
......@@ -87,7 +87,7 @@ public:
*/
CellData2P() :
saturation_({0.0, 0.0}), pressure_({0.0, 0.0}), capillaryPressure_(0), mobility_({0.0, 0.0}), fracFlowFunc_({0.0, 0.0})
saturation_({0.0, 0.0}), pressure_({0.0, 0.0}), capillaryPressure_(0), mobility_({0.0, 0.0}), fracFlowFunc_({0.0, 0.0}), update_(0)
{
}
......@@ -121,27 +121,42 @@ public:
return pressure_[phaseIdx];
}
Scalar pressure(int phaseIdx) const
{
return pressure_[phaseIdx];
}
void setPressure(int phaseIdx, Scalar press)
{
pressure_[phaseIdx] = press;
}
Scalar& globalPressure()
Scalar globalPressure()
{
return pressure_[wPhaseIdx];
}
const Scalar& globalPressure() const
Scalar globalPressure() const
{
return pressure_[wPhaseIdx];
}
void setGlobalPressure(Scalar press)
{
pressure_[wPhaseIdx] = press;
}
//! Return saturation vector
Scalar saturation(int phaseIdx)
{
return saturation_[phaseIdx];
}
Scalar saturation(int phaseIdx) const
{
return saturation_[phaseIdx];
}
void setSaturation(int phaseIdx, Scalar sat)
{
saturation_[phaseIdx] = sat;
......@@ -152,33 +167,48 @@ public:
//////////////////////////////////////////////////////////////
//! Return phase mobilities
Scalar& mobility(int phaseIdx)
Scalar mobility(int phaseIdx)
{
return mobility_[phaseIdx];
}
const Scalar& mobility(int phaseIdx) const
Scalar mobility(int phaseIdx) const
{
return mobility_[phaseIdx];
}
void setMobility(int phaseIdx, Scalar mobility)
{
mobility_[phaseIdx] = mobility;
}
//! Return phase fractional flow functions
Scalar& fracFlowFunc(int phaseIdx)
Scalar fracFlowFunc(int phaseIdx)
{
return fracFlowFunc_[phaseIdx];
}
const Scalar& fracFlowFunc(int phaseIdx) const
Scalar fracFlowFunc(int phaseIdx) const
{
return fracFlowFunc_[phaseIdx];
}
void setFracFlowFunc(int phaseIdx, Scalar fracFlowFunc)
{
fracFlowFunc_[phaseIdx] = fracFlowFunc;
}
//! Return capillary pressure vector
Scalar capillaryPressure()
{
return capillaryPressure_;
}
Scalar capillaryPressure() const
{
return capillaryPressure_;
}
void setCapillaryPressure(Scalar pc)
{
capillaryPressure_ = pc;
......@@ -190,12 +220,12 @@ public:
}
//! Return density vector
Scalar& volumeCorrection()
Scalar volumeCorrection()
{
return update_;
}
const Scalar& volumeCorrection() const
Scalar volumeCorrection() const
{
return update_;
}
......@@ -210,7 +240,6 @@ public:
DUNE_THROW(Dune::NotImplemented,"density function not implemented in cell data of incompressible models!");
}
//! Return density vector
Scalar viscosity(int phaseIdx)
{
......@@ -262,7 +291,7 @@ public:
*/
CellData2P() :
mobility_({0.0, 0.0}), fracFlowFunc_({0.0, 0.0})
mobility_({0.0, 0.0}), fracFlowFunc_({0.0, 0.0}), update_(0)
{
}
......@@ -303,15 +332,20 @@ public:
void setPressure(int phaseIdx, Scalar press)
{
DUNE_THROW(Dune::NotImplemented,"setPressure function defined for compressible models!");
fluidState_.setPressure(phaseIdx, press);
}
Scalar globalPressure()
{
DUNE_THROW(Dune::NotImplemented,"no global pressure defined for compressible models!");
}
Scalar& globalPressure()
Scalar globalPressure() const
{
DUNE_THROW(Dune::NotImplemented,"no global pressure defined for compressible models!");
}
const Scalar& globalPressure() const
void setGlobalPressure(Scalar press)
{
DUNE_THROW(Dune::NotImplemented,"no global pressure defined for compressible models!");
}
......@@ -330,7 +364,7 @@ public:
void setSaturation(int phaseIdx, Scalar sat)
{
DUNE_THROW(Dune::NotImplemented,"setSaturation function defined for compressible models!");
fluidState_.setSaturation(phaseIdx, sat);
}
//////////////////////////////////////////////////////////////
......@@ -338,27 +372,37 @@ public:
//////////////////////////////////////////////////////////////
//! Return phase mobilities
Scalar& mobility(int phaseIdx)
Scalar mobility(int phaseIdx)
{
return mobility_[phaseIdx];
}
const Scalar& mobility(int phaseIdx) const
Scalar mobility(int phaseIdx) const
{
return mobility_[phaseIdx];
}
void setMobility(int phaseIdx, Scalar mobility)
{
mobility_[phaseIdx] = mobility;
}
//! Return phase fractional flow functions
Scalar& fracFlowFunc(int phaseIdx)
Scalar fracFlowFunc(int phaseIdx)
{
return fracFlowFunc_[phaseIdx];
}
const Scalar& fracFlowFunc(int phaseIdx) const
Scalar fracFlowFunc(int phaseIdx) const
{
return fracFlowFunc_[phaseIdx];
}
void setFracFlowFunc(int phaseIdx, Scalar fracFlowFunc)
{
fracFlowFunc_[phaseIdx] = fracFlowFunc;
}
//! Return capillary pressure vector
Scalar capillaryPressure()
{
......@@ -385,6 +429,10 @@ public:
return fluidState_.density(phaseIdx);
}
void setDensity(int phaseIdx, Scalar density)
{
fluidState_.setDensity(phaseIdx, density);
}
//! Return density vector
Scalar viscosity(int phaseIdx)
......@@ -397,21 +445,27 @@ public:
return fluidState_.viscosity(phaseIdx);
}
void setViscosity(int phaseIdx, Scalar viscosity)
{
fluidState_.setViscosity(phaseIdx, viscosity);
}
void setUpdate(Scalar update)
{
update_ = update;
}
//! Return density vector
Scalar& volumeCorrection()
Scalar volumeCorrection()
{
return update_;
}
const Scalar& volumeCorrection() const
Scalar volumeCorrection() const
{
return update_;
}
};
}
#endif
......@@ -158,15 +158,6 @@ public:
return;
}
void serialize()
{
return;
}
void deserialize(double t)
{
return;
}
/*!
* \brief Returns the temperature within the domain.
*
......
......@@ -71,6 +71,7 @@ private:
VelocityVector velocity_[numPhases];
Scalar potential_[2 * dim][numPhases];
bool velocityMarker_[2 * dim];
public:
......@@ -89,6 +90,7 @@ public:
potential_[face][phase] = 0.0;
}
velocityMarker_[face] = false;
}
}
......@@ -97,25 +99,19 @@ public:
////////////////////////////////////////////////////////////
//! Return the velocity
VelocityVector& velocity(int phaseIdx)
const FieldVector& velocity(int phaseIdx, int indexInInside)
{
return velocity_[phaseIdx];
}
const VelocityVector& velocity(int phaseIdx) const
{
return velocity_[phaseIdx];
return velocity_[phaseIdx][indexInInside];
}
//! Return the velocity
FieldVector& velocity(int phaseIdx, int indexInInside)
const FieldVector& velocity(int phaseIdx, int indexInInside) const
{
return velocity_[phaseIdx][indexInInside];
}
const FieldVector& velocity(int phaseIdx, int indexInInside) const
void setVelocity(int phaseIdx, int indexInInside, FieldVector& velocity)
{
return velocity_[phaseIdx][indexInInside];
velocity_[phaseIdx][indexInInside] = velocity;
}
//! Return the velocity
......@@ -131,6 +127,22 @@ public:
+ velocity_[nPhaseIdx][indexInInside];
}
void setVelocityMarker(int indexInInside)
{
velocityMarker_[indexInInside] = true;
}
bool haveVelocity(int indexInInside)
{
return velocityMarker_[indexInInside];
}
void resetVelocityMarker()
{
for (int i = 0; i < 2*dim; i++)
velocityMarker_[i] = false;
}
bool isUpwindCell(int phaseIdx, int indexInInside)
{
return (potential_[indexInInside][phaseIdx] >= 0.);
......@@ -156,11 +168,6 @@ public:
potential_[indexInInside][phaseIdx] = pot;
}
void setUpwindCell(int phaseIdx, int indexInInside, Scalar pot) const
{
potential_[indexInInside][phaseIdx] = pot;
}
};
}
#endif
......@@ -73,9 +73,8 @@ private:
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::Intersection Intersection;
typedef Dune::FieldVector<Scalar, dim> FieldVector;
typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
......@@ -89,26 +88,20 @@ public:
* @param[in] pcGradient gradient of capillary pressure between element I and J
* \return capillary pressure term of the saturation equation
*/
FieldVector operator() (const Element& element, const int indexInInside, Scalar satI, Scalar satJ, const FieldVector& pcGradient) const
void getFlux (FieldVector& flux, const Intersection& intersection, Scalar satI, Scalar satJ, const FieldVector& pcGradient) const
{
ElementPointer element = intersection.inside();
// get global coordinate of cell center
const GlobalPosition& globalPos = element.geometry().center();
const GlobalPosition& globalPos = element->geometry().center();
IntersectionIterator isItEnd = problem_.gridView().iend(element);
IntersectionIterator isIt = problem_.gridView().ibegin(element);
for (; isIt != isItEnd; ++isIt)
{
if(isIt->indexInInside() == indexInInside)
break;
}
int globalIdxI = problem_.variables().index(element);
int globalIdxI = problem_.variables().index(*element);
CellData& CellDataI = problem_.variables().cellData(globalIdxI);
// get geometry type of face
//Dune::GeometryType faceGT = isIt->geometryInInside().type();
Scalar temperature = problem_.temperature(element);
Scalar referencePressure = problem_.referencePressure(element);
Scalar temperature = problem_.temperature(*element);
Scalar referencePressure = problem_.referencePressure(*element);
//get lambda_bar = lambda_n*f_w
Scalar mobBar = 0;
......@@ -126,18 +119,18 @@ public:
fluidState.setPressure(wPhaseIdx, referencePressure);
fluidState.setPressure(nPhaseIdx, referencePressure);
fluidState.setTemperature(temperature);
mobilityWI = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(element), satI);
mobilityWI = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*element), satI);
mobilityWI /= FluidSystem::viscosity(fluidState, wPhaseIdx);
mobilityNWI = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(element), satI);
mobilityNWI = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*element), satI);
mobilityNWI /= FluidSystem::viscosity(fluidState, nPhaseIdx);
}
FieldMatrix meanPermeability(0);
if (isIt->neighbor())
if (intersection.neighbor())
{
// access neighbor
ElementPointer neighborPointer = isIt->outside();
ElementPointer neighborPointer = intersection.outside();
int globalIdxJ = problem_.variables().index(*neighborPointer);
CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
......@@ -156,7 +149,7 @@ public:
// get permeability
problem_.spatialParameters().meanK(meanPermeability,
problem_.spatialParameters().intrinsicPermeability(element),
problem_.spatialParameters().intrinsicPermeability(*element),
problem_.spatialParameters().intrinsicPermeability(*neighborPointer));
......@@ -188,7 +181,7 @@ public:
{
// get permeability
problem_.spatialParameters().meanK(meanPermeability,
problem_.spatialParameters().intrinsicPermeability(element));
problem_.spatialParameters().intrinsicPermeability(*element));
Scalar mobilityWJ = 0;
Scalar mobilityNWJ = 0;
......@@ -198,9 +191,9 @@ public:
fluidState.setPressure(wPhaseIdx, referencePressure);
fluidState.setPressure(nPhaseIdx, referencePressure);
fluidState.setTemperature(temperature);
mobilityWJ = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(element), satJ);
mobilityWJ = MaterialLaw::krw(problem_.spatialParameters().materialLawParams(*element), satJ);
mobilityWJ /= FluidSystem::viscosity(fluidState, wPhaseIdx);
mobilityNWJ = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(element), satJ);
mobilityNWJ = MaterialLaw::krn(problem_.spatialParameters().materialLawParams(*element), satJ);
mobilityNWJ /= FluidSystem::viscosity(fluidState, nPhaseIdx);
Scalar mobWMean = 0.5 * (mobilityWI + mobilityWJ);
......@@ -210,21 +203,18 @@ public:
}
// set result to K*grad(pc)
FieldVector result(0);
meanPermeability.umv(pcGradient, result);
meanPermeability.mv(pcGradient, flux);
// set result to f_w*lambda_n*K*grad(pc)
result *= mobBar;
return result;
flux *= mobBar;
}
/*! @brief Constructs a CapillaryDiffusion object
* @param problem an object of class Dumux::TransportProblem or derived
* @param preComput if preCompute = true previous calculated mobilities are taken, if preCompute = false new mobilities will be computed (for implicit Scheme)
*/
CapillaryDiffusion (Problem& problem, const bool preComput = true)
: DiffusivePart<TypeTag>(problem), problem_(problem), preComput_(preComput)
CapillaryDiffusion (Problem& problem)
: DiffusivePart<TypeTag>(problem), problem_(problem), preComput_(GET_PROP_VALUE(TypeTag, PrecomputedConstRels))
{}
private:
......
......@@ -49,7 +49,8 @@ private:
enum{dim = GridView::dimension, dimWorld = GridView::dimensionworld};
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef typename GridView::Intersection Intersection;
typedef Dune::FieldVector<Scalar, dimWorld> FieldVector;
public:
//! Returns convective term
......@@ -59,10 +60,9 @@ public:
* @param[in] sat saturation of current element
* \return convective term of an advection-diffusion equation
*/
Scalar operator() (const Element& element, const int indexInInside, const Scalar sat) const
Scalar getFlux(const Intersection& intersection, const Scalar sat) const
{
Scalar trivial(0);
return trivial;
return 0.0;
}
//! Returns convective term
/*! Returns convective term for current element face
......@@ -72,11 +72,8 @@ public:
* @param[in] satJ saturation of neighbor element
* \return convective term of an advection-diffusion equation
*/
GlobalPosition operator() (const Element& element, const int indexInInside, const Scalar satI, const Scalar satJ) const
{
GlobalPosition trivial(0);
return trivial;
}
void getFlux(FieldVector& flux, const Intersection& intersection, const Scalar satI, const Scalar satJ) const
{}
//! The constructor
/*
......
......@@ -47,6 +47,7 @@ private:
enum{dim = GridView::dimension};
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::Intersection Intersection;
typedef Dune::FieldVector<Scalar, dim> FieldVector;
public:
......@@ -59,11 +60,8 @@ public:
* @param[in] pcGradient gradient of capillary pressure between element I and J
* \return diffusive term of an advection-diffusion equation
*/
FieldVector operator() (const Element& element, const int indexInInside, Scalar satI, Scalar satJ, const FieldVector& pcGradient) const
{
FieldVector trivial(0);
return trivial;
}
void getFlux(FieldVector& flux, const Intersection& intersection, Scalar satI, Scalar satJ, const FieldVector& pcGradient) const
{}
//! Returns diffusive term
/*! Returns diffusive term for current element face
......@@ -74,12 +72,9 @@ public:
* @param[in] time time
* \return diffusive term of an advection-diffusion equation
*/
FieldVector operator() (const Element& element, const int indexInInside,
void getFlux(FieldVector& flux, const Intersection& intersection,
const Scalar satIntersection, const FieldVector& satGradient, const Scalar time) const
{
FieldVector trivial(0);
return trivial;
}
{}
//! Returns diffusive term
/*! Returns diffusive term for current element face
......@@ -92,13 +87,10 @@ public:
* @param[in] satJ saturation of neighbor element
* \return diffusive term of an advection-diffusion equation
*/
FieldVector operator() (const Element& element, const int indexInInside,
void getFlux(FieldVector& flux, const Intersection& intersection,
const Scalar satIntersection, const FieldVector& satGradient, const Scalar time,
Scalar satI, Scalar satJ) const
{
FieldVector trivial(0);
return trivial;
}
{}
//! The constructor
......
......@@ -97,7 +97,18 @@ public:
* \param element element on which the CFL-criterion is evaluated
* \return fluxFunction for the calculation of the CFL-time step (\f$ 1/F_i\f$)
*/
Scalar getCFLFluxFunction(const GlobalPosition& globalPos, const Element& element)
Scalar getCFLFluxFunction(const Element& element)
{
return 0.0;
}
//! adds a flux to the cfl-criterion evaluation
/*!
* \param globalPos global position
* \param element element on which the CFL-criterion is evaluated
* \return fluxFunction for the calculation of the CFL-time step (\f$ 1/F_i\f$)
*/
Scalar getDt(const Element& element)
{
return 0.0;
}
......
......@@ -426,9 +426,9 @@ public:
}
}
Scalar getCflFluxFunction(const GlobalPosition& globalPos, const Element& element)