Commit f22e76e3 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

richards box model: bring it to stable standards and fix its test

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4182 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent abb4a4c9
......@@ -16,18 +16,20 @@
/*!
* \file
*
* \brief Gives the phase state.
* \brief Calcultes the fluid state from the primary variables in the
* 2p model.
*/
#ifndef DUMUX_2P_PHASE_STATE_HH
#define DUMUX_2P_PHASE_STATE_HH
#ifndef DUMUX_2P_FLUID_STATE_HH
#define DUMUX_2P_FLUID_STATE_HH
#include <dumux/material/fluidstate.hh>
#include <dumux/boxmodels/2p/2pproperties.hh>
#include "2pproperties.hh"
namespace Dumux
{
/*!
* \brief Calcultes the phase state from the primary variables in the
* \brief Calcultes the fluid state from the primary variables in the
* 2p model.
*/
template <class TypeTag>
......@@ -51,7 +53,6 @@ class TwoPFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag, PTAG(Sc
public:
enum { numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)) };
public:
void update(Scalar Sn, Scalar pressW, Scalar pressN, Scalar temperature)
{
Sn_ = Sn;
......
......@@ -83,7 +83,7 @@ public:
}
calculateGradients_(problem, element, elemDat);
calculateVelocities_(problem, element, elemDat);
calculateK_(problem, element, elemDat);
};
public:
......@@ -178,9 +178,9 @@ private:
}
}
void calculateVelocities_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemDat)
void calculateK_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemDat)
{
const SpatialParameters &spatialParams = problem.spatialParameters();
// calculate the intrinsic permeability
......
......@@ -82,9 +82,6 @@ class TwoPModel : public BoxModel<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
......@@ -96,6 +93,8 @@ class TwoPModel : public BoxModel<TypeTag>
wPhaseIdx = Indices::wPhaseIdx,
};
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
public:
/*!
* \brief Returns the relative weight of a primary variable for
......@@ -104,7 +103,7 @@ public:
Scalar primaryVarWeight(int vertIdx, int pvIdx) const
{
if (Indices::pressureIdx == pvIdx)
return 1e-7;
return 1e-5;
return 1;
}
......@@ -176,8 +175,8 @@ public:
writer.addVertexData(Sn, "Sn");
writer.addVertexData(Sw, "Sw");
writer.addVertexData(pW, "pg");
writer.addVertexData(pN, "pn");
writer.addVertexData(pW, "pw");
writer.addVertexData(pC, "pc");
writer.addVertexData(rhoW, "rhoW");
writer.addVertexData(rhoN, "rhoN");
......
......@@ -23,14 +23,13 @@
#ifndef DUMUX_2P_PROPERTY_DEFAULTS_HH
#define DUMUX_2P_PROPERTY_DEFAULTS_HH
#include <dumux/boxmodels/common/boxproperties.hh>
#include "2pmodel.hh"
#include "2pproblem.hh"
#include "2pindices.hh"
#include "2pfluxvariables.hh"
#include "2pvolumevariables.hh"
#include "2pfluidstate.hh"
#include "2pproperties.hh"
#include <dumux/material/fluidsystems/2p_system.hh>
namespace Dumux
......
......@@ -41,25 +41,24 @@ template <class TypeTag>
class TwoPVolumeVariables : public BoxVolumeVariables<TypeTag>
{
typedef BoxVolumeVariables<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
enum {
numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation)),
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
enum {
pwSn = Indices::pwSn,
pnSw = Indices::pnSw,
......@@ -70,21 +69,8 @@ class TwoPVolumeVariables : public BoxVolumeVariables<TypeTag>
nPhaseIdx = Indices::nPhaseIdx
};
typedef typename GET_PROP(TypeTag, PTAG(ReferenceElements)) RefElemProp;
typedef typename RefElemProp::Container ReferenceElements;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GridView::template Codim<0>::Entity Element;
typedef TwoPFluidState<TypeTag> FluidState;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, dim> LocalPosition;
public:
/*!
......
......@@ -38,129 +38,142 @@ class RichardsFluxVariables
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(RichardsIndices)) Indices;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
};
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename FVElementGeometry::SubControlVolume SCV;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(RichardsIndices)) Indices;
typedef Dune::FieldVector<Scalar, dimWorld> Vector;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
typedef typename GridView::template Codim<0>::Entity Element;
public:
RichardsFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
int faceIdx,
const ElementVolumeVariables &elemDat)
: fvElemGeom(elemGeom)
const Element &element,
const FVElementGeometry &elemGeom,
int faceIdx,
const ElementVolumeVariables &elemVolVars)
: fvElemGeom_(elemGeom)
{
face = &fvElemGeom.subContVolFace[faceIdx];
int i = face->i;
int j = face->j;
insideSCV = &fvElemGeom.subContVol[i];
outsideSCV = &fvElemGeom.subContVol[j];
scvfIdx_ = faceIdx;
densityWAtIP = 0;
pressureGrad = Scalar(0);
calculateGradients_(problem, element, elemDat);
calculateVelocities_(problem, element, elemDat);
calculateGradients_(problem, element, elemVolVars);
calculateK_(problem, element, elemVolVars);
};
private:
/*
* \brief Return the intrinsic permeability.
*/
const Tensor &intrinsicPermeability() const
{ return K_; }
/*!
* \brief Return the pressure potential gradient.
*/
const Vector &potentialGradW() const
{ return potentialGrad_; }
/*!
* \brief Given the intrinisc permeability times the pressure
* potential gradient and SCV face normal for a phase,
* return the local index of the downstream control volume
* for a given phase.
*/
int downstreamIdx(Scalar normalFlux) const
{ return (normalFlux >= 0)?face().j:face().i; }
/*!
* \brief Given the intrinisc permeability times the pressure
* potential gradient and SCV face normal for a phase,
* return the local index of the upstream control volume
* for a given phase.
*/
int upstreamIdx(Scalar normalFlux) const
{ return (normalFlux > 0)?face().i:face().j; }
const SCVFace &face() const
{ return fvElemGeom_.subContVolFace[scvfIdx_]; }
protected:
void calculateGradients_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemDat)
const ElementVolumeVariables &elemVolVars)
{
potentialGrad_ = 0.0;
// calculate gradients
GlobalPosition tmp(0.0);
for (int idx = 0;
idx < fvElemGeom.numVertices;
idx < fvElemGeom_.numVertices;
idx++) // loop over adjacent vertices
{
// FE gradient at vertex index
const LocalPosition &feGrad = face->grad[idx];
const Vector &feGrad = face().grad[idx];
// the pressure gradient
tmp = feGrad;
tmp *= elemDat[idx].pW;
pressureGrad += tmp;
// fluid density
densityWAtIP +=
elemDat[idx].densityW*face->shapeValue[idx];
Vector tmp(feGrad);
tmp *= elemVolVars[idx].pressure(wPhaseIdx);
potentialGrad_ += tmp;
}
///////////////
// correct the pressure gradients by the hydrostatic
// pressure due to gravity
tmp = problem.gravity();
tmp *= densityWAtIP;
pressureGrad -= tmp;
///////////////
// calculate the phase density at the integration point. we
// only do this if the wetting phase is present in both cells
Scalar SI = elemVolVars[face().i].saturation(wPhaseIdx);
Scalar SJ = elemVolVars[face().j].saturation(wPhaseIdx);
Scalar rhoI = elemVolVars[face().i].density(wPhaseIdx);
Scalar rhoJ = elemVolVars[face().j].density(wPhaseIdx);
Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
if (fI + fJ == 0)
// doesn't matter because no wetting phase is present in
// both cells!
fI = fJ = 0.5;
Scalar density = (fI*rhoI + fJ*rhoJ)/(fI + fJ);
Vector tmp(problem.gravity());
tmp *= density;
potentialGrad_ -= tmp;
}
void calculateVelocities_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemDat)
void calculateK_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemDat)
{
// calculate the permeability tensor. TODO: this should be
// more flexible
typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
Tensor K;
problem.spatialParameters().meanK(K,
problem.spatialParameters().intrinsicPermeability(element,
fvElemGeom,
face->i),
problem.spatialParameters().intrinsicPermeability(element,
fvElemGeom,
face->j));
// temporary vector for the Darcy velocity
GlobalPosition vDarcy;
K.mv(pressureGrad, vDarcy); // vDarcy = K * grad p
vDarcyNormal = vDarcy*face->normal;
// set the upstream and downstream vertices
upstreamIdx = face->i;
downstreamIdx = face->j;
if (vDarcyNormal > 0)
std::swap(upstreamIdx, downstreamIdx);
const SpatialParameters &sp = problem.spatialParameters();
// calculate the intrinsic permeability
sp.meanK(K_,
sp.intrinsicPermeability(element,
fvElemGeom_,
face().i),
sp.intrinsicPermeability(element,
fvElemGeom_,
face().j));
}
public:
const FVElementGeometry &fvElemGeom;
const SCVFace *face;
const SCV *insideSCV;
const SCV *outsideSCV;
const FVElementGeometry &fvElemGeom_;
int scvfIdx_;
// gradients
GlobalPosition pressureGrad;
// density of the fluid at the integration point
Scalar densityWAtIP;
// darcy velocity in direction of the face normal
Scalar vDarcyNormal;
Vector potentialGrad_;
// local index of the upwind vertex
int upstreamIdx;
// local index of the downwind vertex
int downstreamIdx;
// intrinsic permeability
Tensor K_;
};
} // end namepace
......
// $Id: richardsproperties.hh 3840 2010-07-15 10:14:15Z bernd $
/*****************************************************************************
* Copyright (C) 2009 by Andreas Lauser *
* Institute of Hydraulic Engineering *
* 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, as long as this copyright notice *
* is included in its original form. *
* *
* This program is distributed WITHOUT ANY WARRANTY. *
*****************************************************************************/
/*!
* \file
*
* \brief Indices for the Richards model.
*/
#ifndef DUMUX_RICHARDS_INDICES_HH
#define DUMUX_RICHARDS_INDICES_HH
namespace Dumux
{
/*!
* \addtogroup RichardsModel
*/
// \{
/*!
* \brief Indices for the Richards model.
*/
struct RichardsIndices
{
//////////
// primary variable indices
//////////
//! wetting phase pressure
static const int pwIdx = 0;
//////////
// equation indices
//////////
//! continuity
static const int contiEqIdx = 0;
//////////
// phase indices
//////////
static const int wPhaseIdx = 0;
static const int nPhaseIdx = 1;
};
// \}
} // end namepace
#endif
......@@ -38,40 +38,25 @@ namespace Dumux
template<class TypeTag>
class RichardsLocalResidual : public BoxLocalResidual<TypeTag>
{
typedef RichardsLocalResidual<TypeTag> ThisType;
typedef BoxLocalResidual<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GridView::template Codim<0>::Entity Element;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(RichardsIndices)) Indices;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
pW = Indices::pW,
};
contiEqIdx = Indices::contiEqIdx,
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
typedef Dune::FieldVector<Scalar, dim> LocalPosition;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
};
typedef Dune::FieldVector<Scalar, dimWorld> Vector;
static const Scalar mobilityUpwindAlpha = GET_PROP_VALUE(TypeTag, PTAG(MobilityUpwindAlpha));
public:
/*!
* \brief Evaluate the rate of change of all conservation
......@@ -88,15 +73,16 @@ public:
// used. The secondary variables are used accordingly. This
// is required to compute the derivative of the storage term
// using the implicit euler method.
const ElementVolumeVariables &elemDat = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
const VolumeVariables &vertDat = elemDat[scvIdx];
const VolumeVariables &volVars =
usePrevSol ?
this->prevVolVars_(scvIdx) :
this->curVolVars_(scvIdx);
// partial time derivative of the wetting phase mass
result[pW] =
vertDat.densityW
* vertDat.porosity
* this->prevVolVars_()[scvIdx].dSwdpC // TODO: use derivative for the current solution
* (vertDat.pNreference - vertDat.pW);
result[contiEqIdx] =
volVars.density(wPhaseIdx)
* volVars.saturation(wPhaseIdx)
* volVars.porosity();
}
......@@ -106,25 +92,33 @@ public:
*/
void computeFlux(PrimaryVariables &flux, int faceId) const
{
FluxVariables vars(this->problem_(),
this->elem_(),
this->fvElemGeom_(),
faceId,
this->curVolVars_());
FluxVariables fluxVars(this->problem_(),
this->elem_(),
this->fvElemGeom_(),
faceId,
this->curVolVars_());
// calculate the flux in the normal direction of the
// current sub control volume face
Vector tmpVec;
fluxVars.intrinsicPermeability().mv(fluxVars.potentialGradW(),
tmpVec);
Scalar normalFlux = - (tmpVec*fluxVars.face().normal);
// data attached to upstream and the downstream vertices
const VolumeVariables &up = this->curVolVars_(vars.upstreamIdx);
const VolumeVariables &dn = this->curVolVars_(vars.downstreamIdx);
flux[pW] =
vars.vDarcyNormal*
( mobilityUpwindAlpha*
( up.densityW *
up.mobilityW)
+
(1 - mobilityUpwindAlpha)*
( dn.densityW*
dn.mobilityW));
// of the current phase
const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(normalFlux));
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(normalFlux));
flux[contiEqIdx] =
normalFlux
*
(( mobilityUpwindAlpha)*up.density(wPhaseIdx)*up.mobility(wPhaseIdx)
+
(1 - mobilityUpwindAlpha)*dn.density(wPhaseIdx)*dn.mobility(wPhaseIdx));
// we need the flux from i to j instead of the other way
flux *= -1;
}
/*!
......@@ -137,21 +131,6 @@ public:
this->fvElemGeom_(),
localVertexIdx);
}
/*!
* \brief Return the temperature given the solution vector of a
* finite volume.
*/
template <class PrimaryVariables>
Scalar temperature(const PrimaryVariables &sol)
{ return this->problem_.temperature(); /* constant temperature */ }
private:
ThisType &asImp_()