Commit 582da479 authored by Alexander Kissinger's avatar Alexander Kissinger
Browse files

New naming convention for richards box-model


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@8166 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent e7e33255
......@@ -46,7 +46,7 @@ class RichardsFluxVariables
{
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParams;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
......@@ -55,11 +55,11 @@ class RichardsFluxVariables
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GridView::template Codim<0>::Entity Element;
enum { dimWorld = GridView::dimensionworld};
enum { dim = GridView::dimension};
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, dimWorld> Vector;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
typedef Dune::FieldVector<Scalar, dim> DimVector;
typedef Dune::FieldMatrix<Scalar, dim, dim> DimTensor;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
......@@ -80,13 +80,13 @@ public:
*/
RichardsFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &fvElemGeom,
int scvfIdx,
const FVElementGeometry &fvGeometry,
const int faceIdx,
const ElementVolumeVariables &elemVolVars,
const bool onBoundary = false)
: fvElemGeom_(fvElemGeom), onBoundary_(onBoundary)
: fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
{
scvfIdx_ = scvfIdx;
calculateGradients_(problem, element, elemVolVars);
calculateK_(problem, element, elemVolVars);
......@@ -95,13 +95,13 @@ public:
/*!
* \brief Return the intrinsic permeability \f$\mathrm{[m^2]}\f$.
*/
const Tensor &intrinsicPermeability() const
const DimTensor &intrinsicPermeability() const
{ return K_; }
/*!
* \brief Return the pressure potential gradient \f$\mathrm{[Pa/m]}\f$
*/
const Vector &potentialGradW() const
const DimVector &potentialGradW() const
{ return potentialGrad_; }
/*!
......@@ -128,18 +128,22 @@ public:
int upstreamIdx(Scalar normalFlux) const
{ return (normalFlux > 0)?face().i:face().j; }
/*!
* \brief Returns the face of the element's finite volume geometry
* which the flux variables object looks at
* \brief The face of the current sub-control volume. This may be either
* an inner sub-control-volume face or a face on the boundary.
*/
const SCVFace &face() const
{
if (onBoundary_)
return fvElemGeom_.boundaryFace[scvfIdx_];
return fvGeometry_.boundaryFace[faceIdx_];
else
return fvElemGeom_.subContVolFace[scvfIdx_];
return fvGeometry_.subContVolFace[faceIdx_];
}
protected:
void calculateGradients_(const Problem &problem,
const Element &element,
......@@ -148,15 +152,16 @@ protected:
potentialGrad_ = 0.0;
// calculate gradients
for (int idx = 0;
idx < fvElemGeom_.numVertices;
idx < fvGeometry_.numVertices;
idx++) // loop over adjacent vertices
{
// FE gradient at vertex index
const Vector &feGrad = face().grad[idx];
const DimVector &feGrad = face().grad[idx];
// the pressure gradient
Vector tmp(feGrad);
DimVector tmp(feGrad);
tmp *= elemVolVars[idx].pressure(wPhaseIdx);
potentialGrad_ += tmp;
}
......@@ -180,8 +185,8 @@ protected:
// estimate the gravitational acceleration at a given SCV face
// using the arithmetic mean
Vector f(problem.boxGravity(element, fvElemGeom_, face().i));
f += problem.boxGravity(element, fvElemGeom_, face().j);
DimVector f(problem.boxGravity(element, fvGeometry_, face().i));
f += problem.boxGravity(element, fvGeometry_, face().j);
f /= 2;
// make it a force
......@@ -196,26 +201,26 @@ protected:
const Element &element,
const ElementVolumeVariables &elemDat)
{
const SpatialParameters &sp = problem.spatialParameters();
const SpatialParams &spatialParams = problem.spatialParams();
// calculate the intrinsic permeability
sp.meanK(K_,
sp.intrinsicPermeability(element,
fvElemGeom_,
spatialParams.meanK(K_,
spatialParams.intrinsicPermeability(element,
fvGeometry_,
face().i),
sp.intrinsicPermeability(element,
fvElemGeom_,
spatialParams.intrinsicPermeability(element,
fvGeometry_,
face().j));
}
const FVElementGeometry &fvElemGeom_;
int scvfIdx_;
const FVElementGeometry &fvGeometry_;
const int faceIdx_;
const bool onBoundary_;
// gradients
Vector potentialGrad_;
DimVector potentialGrad_;
// intrinsic permeability
Tensor K_;
DimTensor K_;
};
} // end namepace
......
......@@ -55,10 +55,10 @@ class RichardsLocalResidual : public BoxLocalResidual<TypeTag>
};
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
enum { dimWorld = GridView::dimensionworld};
enum { dim = GridView::dimension};
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef Dune::FieldVector<Scalar, dimWorld> Vector;
typedef Dune::FieldVector<Scalar, dim> DimVector;
public:
/*!
......@@ -80,12 +80,12 @@ public:
*
* This function should not include the source and sink terms.
*
* \param result Stores the average mass per unit volume for each phase \f$\mathrm{[kg/m^3]}\f$
* \param storage Stores the average mass per unit volume for each phase \f$\mathrm{[kg/m^3]}\f$
* \param scvIdx The sub control volume index of the current element
* \param usePrevSol Calculate the storage term of the previous solution
* instead of the model's current solution.
*/
void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
void computeStorage(PrimaryVariables &storage, int scvIdx, bool usePrevSol) const
{
// if flag usePrevSol is set, the solution from the previous
// time step is used, otherwise the current solution is
......@@ -98,7 +98,7 @@ public:
this->curVolVars_(scvIdx);
// partial time derivative of the wetting phase mass
result[contiEqIdx] =
storage[contiEqIdx] =
volVars.density(wPhaseIdx)
* volVars.saturation(wPhaseIdx)
* volVars.porosity();
......@@ -117,18 +117,18 @@ public:
* \param onBoundary A boolean variable to specify whether the flux variables
* are calculated for interior SCV faces or boundary faces, default=false
*/
void computeFlux(PrimaryVariables &flux, int scvfIdx, const bool onBoundary=false) const
void computeFlux(PrimaryVariables &flux, int faceIdx, bool onBoundary=false) const
{
FluxVariables fluxVars(this->problem_(),
this->elem_(),
this->fvElemGeom_(),
scvfIdx,
this->element_(),
this->fvGeometry_(),
faceIdx,
this->curVolVars_(),
onBoundary);
// calculate the flux in the normal direction of the
// current sub control volume face
Vector tmpVec;
DimVector tmpVec;
fluxVars.intrinsicPermeability().mv(fluxVars.potentialGradW(),
tmpVec);
Scalar normalFlux = -(tmpVec*fluxVars.face().normal);
......@@ -149,16 +149,16 @@ public:
/*!
* \brief Calculate the source term of the equation
*
* \param q Stores the average source term of all phases inside a
* \param source Stores the average source term of all phases inside a
* sub-control volume of the current element \f$\mathrm{[kg/(m^3 * s)]}\f$
* \param scvIdx The sub control volume index inside the current
* element
*/
void computeSource(PrimaryVariables &q, int scvIdx)
void computeSource(PrimaryVariables &source, int scvIdx)
{
this->problem_().boxSDSource(q,
this->elem_(),
this->fvElemGeom_(),
this->problem_().boxSDSource(source,
this->element_(),
this->fvGeometry_(),
scvIdx,
this->curVolVars_());
}
......
......@@ -46,7 +46,7 @@ class RichardsNewtonController : public NewtonController<TypeTag>
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParams;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
......@@ -90,21 +90,21 @@ public:
return;
// clamp saturation change to at most 20% per iteration
FVElementGeometry fvElemGeom;
FVElementGeometry fvGeomtry;
const GridView &gv = this->problem_().gridView();
ElementIterator eIt = gv.template begin<0>();
const ElementIterator eEndIt = gv.template end<0>();
for (; eIt != eEndIt; ++eIt) {
fvElemGeom.update(gv, *eIt);
for (int i = 0; i < fvElemGeom.numVertices; ++i) {
fvGeomtry.update(gv, *eIt);
for (int i = 0; i < fvGeomtry.numVertices; ++i) {
int globI = this->problem_().vertexMapper().map(*eIt, i, dim);
// calculate the old wetting phase saturation
const SpatialParameters &sp = this->problem_().spatialParameters();
const MaterialLawParams &mp = sp.materialLawParams(*eIt, fvElemGeom, i);
const SpatialParams &spatialParams = this->problem_().spatialParams();
const MaterialLawParams &mp = spatialParams.materialLawParams(*eIt, fvGeomtry, i);
Scalar pcMin = MaterialLaw::pC(mp, 1.0);
Scalar pW = uLastIter[globI][pwIdx];
Scalar pN = std::max(this->problem_().referencePressure(*eIt, fvElemGeom, i),
Scalar pN = std::max(this->problem_().referencePressure(*eIt, fvGeomtry, i),
pW + pcMin);
Scalar pcOld = pN - pW;
Scalar SwOld = std::max<Scalar>(0.0, MaterialLaw::Sw(mp, pcOld));
......
......@@ -84,7 +84,7 @@ public:
* \param scvIdx The local index of the sub control volume inside the element
*/
Scalar referencePressure(const Element &element,
const FVElementGeometry fvGeom,
const FVElementGeometry fvGeometry,
int scvIdx) const
{ DUNE_THROW(Dune::NotImplemented, "referencePressure() method not implemented by the actual problem"); };
// \}
......
......@@ -85,37 +85,37 @@ public:
* \param isOldSol Specifies whether the solution is from
* the previous time step or from the current one
*/
void update(const PrimaryVariables &priVars,
void update(const PrimaryVariables &primaryVariables,
const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
const FVElementGeometry &fvGeometry,
int scvIdx,
bool isOldSol)
{
assert(!FluidSystem::isLiquid(nPhaseIdx));
ParentType::update(priVars,
ParentType::update(primaryVariables,
problem,
element,
elemGeom,
fvGeometry,
scvIdx,
isOldSol);
completeFluidState(priVars, problem, element, elemGeom, scvIdx, fluidState_);
completeFluidState(primaryVariables, problem, element, fvGeometry, scvIdx, fluidState_);
//////////
// specify the other parameters
//////////
const MaterialLawParams &matParams =
problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
relativePermeabilityWetting_ =
MaterialLaw::krw(matParams,
fluidState_.saturation(wPhaseIdx));
porosity_ = problem.spatialParameters().porosity(element, elemGeom, scvIdx);
porosity_ = problem.spatialParams().porosity(element, fvGeometry, scvIdx);
// energy related quantities not contained in the fluid state
asImp_().updateEnergy_(priVars, problem, element, elemGeom, scvIdx, isOldSol);
asImp_().updateEnergy_(primaryVariables, problem, element, fvGeometry, scvIdx, isOldSol);
}
/*!
......@@ -124,19 +124,19 @@ public:
static void completeFluidState(const PrimaryVariables& primaryVariables,
const Problem& problem,
const Element& element,
const FVElementGeometry& elementGeometry,
const FVElementGeometry& fvGeometry,
int scvIdx,
FluidState& fluidState)
{
// temperature
Scalar t = Implementation::temperature_(primaryVariables, problem, element,
elementGeometry, scvIdx);
fvGeometry, scvIdx);
fluidState.setTemperature(t);
// pressures
Scalar pnRef = problem.referencePressure(element, elementGeometry, scvIdx);
Scalar pnRef = problem.referencePressure(element, fvGeometry, scvIdx);
const MaterialLawParams &matParams =
problem.spatialParameters().materialLawParams(element, elementGeometry, scvIdx);
problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
Scalar minPc = MaterialLaw::pC(matParams, 1.0);
fluidState.setPressure(wPhaseIdx, primaryVariables[pwIdx]);
fluidState.setPressure(nPhaseIdx, std::max(pnRef, primaryVariables[pwIdx] + minPc));
......@@ -257,13 +257,13 @@ public:
{ return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx); }
protected:
static Scalar temperature_(const PrimaryVariables &priVars,
static Scalar temperature_(const PrimaryVariables &primaryVariables,
const Problem& problem,
const Element &element,
const FVElementGeometry &elemGeom,
const FVElementGeometry &fvGeometry,
int scvIdx)
{
return problem.boxTemperature(element, elemGeom, scvIdx);
return problem.boxTemperature(element, fvGeometry, scvIdx);
}
/*!
......@@ -272,8 +272,8 @@ protected:
void updateEnergy_(const PrimaryVariables &sol,
const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
int vertIdx,
const FVElementGeometry &fvGeometry,
int scvIdx,
bool isOldSol)
{ }
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment