Commit 0606de03 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

various cleanups

- make all headers in dumux/common and dumux/boxmodels/common self sustained,
  i.e. that they can be included without any preconditions
- remove most irrelevant includes in the common code
- bring the 1p and 1p2c box model up to stable standard (accessor functions in
  the volume and flux variable classes, make the 1p2c model work for
  compressible fluids, various fixes the 1p2c test problem)
- change the isfluid-trail fluid system to standard SI units
- play around with timestep control if timestep ramp-up is enabled
  (problem is still not solved, though)
- split the *properties.hh headers into the declaration of the properties and
  definition of defaults for the common part of the box models as well as for
  the 1p and 1p2c models
- put the indices structure for the 1p and 1p2c models into separate files
- introduce a base class for volume variables of the box models

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4180 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 66b67cdc
......@@ -23,6 +23,8 @@
#ifndef DUMUX_1P_FLUX_VARIABLES_HH
#define DUMUX_1P_FLUX_VARIABLES_HH
#include "1pproperties.hh"
#include <dumux/common/math.hh>
namespace Dumux
......@@ -61,6 +63,9 @@ class OnePFluxVariables
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePIndices)) Indices;
typedef Dune::FieldVector<Scalar, dim> Vector;
typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
public:
OnePFluxVariables(const Problem &problem,
const Element &element,
......@@ -71,24 +76,53 @@ public:
{
scvfIdx_ = faceIdx;
densityAtIP = 0;
viscosityAtIP = 0;
pressureGrad = Scalar(0);
calculateGradients_(problem, element, elemDat);
calculateVelocities_(problem, element, elemDat);
calculateK_(problem, element);
};
const SCVFace &face() const
{ return fvElemGeom_.subContVolFace[scvfIdx_]; }
/*!
* \brief Return the intrinsic permeability.
*/
const Tensor &intrinsicPermeability() const
{ return K_; }
/*!
* \brief Return the pressure potential gradient.
*/
const Vector &potentialGrad() 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 upstream control volume
* for a given phase.
*/
int upstreamIdx(Scalar normalFlux) const
{ return (normalFlux >= 0)?face().i:face().j; }
/*!
* \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; }
private:
void calculateGradients_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemDat)
{
// calculate gradients
GlobalPosition tmp(0.0);
potentialGrad_ = 0.0;
Scalar densityAtIP = 0.0;
// calculate potential gradient
for (int idx = 0;
idx < fvElemGeom_.numVertices;
idx++) // loop over adjacent vertices
......@@ -97,72 +131,50 @@ private:
const LocalPosition &feGrad = face().grad[idx];
// the pressure gradient
tmp = feGrad;
tmp *= elemDat[idx].pressure;
pressureGrad += tmp;
Vector tmp(feGrad);
tmp *= elemDat[idx].pressure();
potentialGrad_ += tmp;
// fluid density
densityAtIP +=
elemDat[idx].density*face().shapeValue[idx];
// fluid viscosity
viscosityAtIP +=
elemDat[idx].viscosity*face().shapeValue[idx];
elemDat[idx].density()*face().shapeValue[idx];
}
// correct the pressure gradients by the hydrostatic
// pressure due to gravity
tmp = problem.gravity();
Vector tmp(problem.gravity());
tmp *= densityAtIP;
pressureGrad -= tmp;
potentialGrad_ -= tmp;
}
void calculateVelocities_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemDat)
void calculateK_(const Problem &problem,
const Element &element)
{
const SpatialParameters &spatialParams = problem.spatialParameters();
typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
Tensor K;
spatialParams.meanK(K,
spatialParams.intrinsicPermeability(element,
fvElemGeom_,
face().i),
spatialParams.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);
spatialParams.meanK(K_,
spatialParams.intrinsicPermeability(element,
fvElemGeom_,
face().i),
spatialParams.intrinsicPermeability(element,
fvElemGeom_,
face().j));
}
public:
protected:
const FVElementGeometry &fvElemGeom_;
int scvfIdx_;
// gradients
GlobalPosition pressureGrad;
// density of the fluid at the integration point
Scalar densityAtIP;
// viscosity of the fluid at the integration point
Scalar viscosityAtIP;
Vector potentialGrad_;
// darcy velocity in direction of the face normal
Scalar vDarcyNormal;
// intrinsic permeability
Tensor K_;
// local index of the upwind vertex
int upstreamIdx;
int upstreamIdx_;
// local index of the downwind vertex
int downstreamIdx;
int downstreamIdx_;
};
} // end namepace
......
// $Id: 1pproperties.hh 3784 2010-06-24 13:43:57Z bernd $
/*****************************************************************************
* Copyright (C) 2008-2010 by Andreas Lauser *
* Copyright (C) 2008 by Bernd Flemisch *
* 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 single phase model.
*/
#ifndef DUMUX_1P_INDICES_HH
#define DUMUX_1P_INDICES_HH
namespace Dumux
{
/*!
* \brief Indices for the single phase model.
*/
struct OnePIndices
{
static const int pressureIdx = 0;
};
} // end namepace
#endif
......@@ -62,13 +62,13 @@ class OnePLocalResidual : public BoxLocalResidual<TypeTag>
pressureIdx = Indices::pressureIdx,
};
static const Scalar upwindWeight = GET_PROP_VALUE(TypeTag, PTAG(UpwindWeight));
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;
typedef Dune::FieldVector<Scalar, dimWorld> Vector;
public:
......@@ -92,7 +92,7 @@ public:
const VolumeVariables &volVars = elemVars[scvIdx];
// partial time derivative of the wetting phase mass
result[pressureIdx] = volVars.density * volVars.porosity;
result[pressureIdx] = volVars.density() * volVars.porosity();
}
......@@ -102,13 +102,25 @@ public:
*/
void computeFlux(PrimaryVariables &flux, int faceId) const
{
FluxVariables vars(this->problem_(),
this->elem_(),
this->fvElemGeom_(),
faceId,
this->curVolVars_());
flux[pressureIdx] = vars.densityAtIP * vars.vDarcyNormal / vars.viscosityAtIP;
FluxVariables fluxVars(this->problem_(),
this->elem_(),
this->fvElemGeom_(),
faceId,
this->curVolVars_());
Vector tmpVec;
fluxVars.intrinsicPermeability().mv(fluxVars.potentialGrad(),
tmpVec);
Scalar normalFlux = - (tmpVec*fluxVars.face().normal);
const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(normalFlux));
const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(normalFlux));
flux[pressureIdx] =
(( upwindWeight)*(up.density()/up.viscosity())
+
(1 - upwindWeight)*(dn.density()/dn.viscosity()))
*
normalFlux;
}
/*!
......
......@@ -118,7 +118,7 @@ public:
i,
false);
(*p)[globalIdx] = volVars.pressure;
(*p)[globalIdx] = volVars.pressure();
};
}
......@@ -128,4 +128,6 @@ public:
};
}
#include "1ppropertydefaults.hh"
#endif
......@@ -30,29 +30,6 @@ namespace Dumux
* \addtogroup OnePBoxModel
*/
// \{
////////////////////////////////
// forward declarations
////////////////////////////////
template<class TypeTag>
class OnePBoxModel;
template<class TypeTag>
class OnePLocalResidual;
template <class TypeTag>
class OnePVolumeVariables;
template <class TypeTag>
class OnePFluxVariables;
/*!
* \brief Indices for the single phase model.
*/
struct OnePIndices
{
static const int pressureIdx = 0;
};
///////////////////////////////////////////////////////////////////////////
// properties for the isothermal single phase model
///////////////////////////////////////////////////////////////////////////
......@@ -74,30 +51,7 @@ NEW_PROP_TAG(OnePIndices); //!< Enumerations for the 1p models
NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters object
NEW_PROP_TAG(Fluid); //!< The fluid for the single-phase problems
NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
//////////////////////////////////////////////////////////////////
// Properties
//////////////////////////////////////////////////////////////////
SET_INT_PROP(BoxOneP, NumEq, 1);
SET_INT_PROP(BoxOneP, NumPhases, 1);
//! The local residual function
SET_TYPE_PROP(BoxOneP,
LocalResidual,
OnePLocalResidual<TypeTag>);
//! the Model property
SET_TYPE_PROP(BoxOneP, Model, OnePBoxModel<TypeTag>);
//! the VolumeVariables property
SET_TYPE_PROP(BoxOneP, VolumeVariables, OnePVolumeVariables<TypeTag>);
//! the FluxVariables property
SET_TYPE_PROP(BoxOneP, FluxVariables, OnePFluxVariables<TypeTag>);
//! The indices required by the isothermal single-phase model
SET_TYPE_PROP(BoxOneP, OnePIndices, OnePIndices);
NEW_PROP_TAG(UpwindWeight); //!< Returns weight of the upwind cell when calculating fluxes
// \}
};
......
// $Id: 1pproperties.hh 3784 2010-06-24 13:43:57Z bernd $
/*****************************************************************************
* Copyright (C) 2008-2009 by Andreas Lauser *
* Copyright (C) 2008 by Bernd Flemisch *
* 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 Defines the properties required for the onephase BOX model.
*/
#ifndef DUMUX_1P_PROPERTY_DEFAULTS_HH
#define DUMUX_1P_PROPERTY_DEFAULTS_HH
#include <dumux/boxmodels/common/boxproperties.hh>
#include "1pmodel.hh"
#include "1plocalresidual.hh"
#include "1pvolumevariables.hh"
#include "1pfluxvariables.hh"
#include "1pindices.hh"
namespace Dumux
{
///////////////////////////////////////////////////////////////////////////
// default property values for the isothermal single phase model
///////////////////////////////////////////////////////////////////////////
namespace Properties {
SET_INT_PROP(BoxOneP, NumEq, 1);
SET_INT_PROP(BoxOneP, NumPhases, 1);
//! The local residual function
SET_TYPE_PROP(BoxOneP,
LocalResidual,
OnePLocalResidual<TypeTag>);
//! the Model property
SET_TYPE_PROP(BoxOneP, Model, OnePBoxModel<TypeTag>);
//! the VolumeVariables property
SET_TYPE_PROP(BoxOneP, VolumeVariables, OnePVolumeVariables<TypeTag>);
//! the FluxVariables property
SET_TYPE_PROP(BoxOneP, FluxVariables, OnePFluxVariables<TypeTag>);
//! The indices required by the isothermal single-phase model
SET_TYPE_PROP(BoxOneP, OnePIndices, OnePIndices);
//! The weight of the upwind control volume when calculating
//! fluxes. Use central differences by default.
SET_SCALAR_PROP(BoxOneP, UpwindWeight, 0.5);
// \}
};
} // end namepace
#endif
......@@ -23,6 +23,7 @@
#define DUMUX_1P_VOLUME_VARIABLES_HH
#include "1pproperties.hh"
#include <dumux/boxmodels/common/boxvolumevariables.hh>
namespace Dumux
{
......@@ -33,8 +34,10 @@ namespace Dumux
* finite volume in the one-phase model.
*/
template <class TypeTag>
class OnePVolumeVariables
class OnePVolumeVariables : public BoxVolumeVariables<TypeTag>
{
typedef BoxVolumeVariables<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
......@@ -71,40 +74,74 @@ public:
int scvIdx,
bool isOldSol)
{
primaryVars_ = priVars;
ParentType::update(priVars, problem, element, elemGeom, scvIdx, isOldSol);
typedef Indices I;
this->asImp_().updateTemperature_(priVars,
element,
elemGeom,
scvIdx,
problem);
Scalar temperature = problem.temperature(element, elemGeom, scvIdx);
pressure = priVars[I::pressureIdx];
density = Fluid::density(temperature, pressure);
viscosity = Fluid::viscosity(temperature, pressure);
pressure_ = priVars[Indices::pressureIdx];
density_ = Fluid::density(temperature(), pressure());
viscosity_ = Fluid::viscosity(temperature(), pressure());
// porosity
porosity = problem.spatialParameters().porosity(element,
elemGeom,
scvIdx);
porosity_ = problem.spatialParameters().porosity(element,
elemGeom,
scvIdx);
};
/*!
* \brief Returns temperature inside the sub-control volume.
*
* Note that we assume thermodynamic equilibrium, i.e. the
* temperature of the rock matrix and of all fluid phases are
* identical.
*/
Scalar temperature() const
{ return temperature_; }
/*!
* \brief Returns the effective pressure of a given phase within
* the control volume.
*/
Scalar pressure() const
{ return pressure_; }
/*!
* \brief Sets the evaluation point used in the by the local jacobian.
* \brief Returns the mass density of a given phase within the
* control volume.
*/
void setEvalPoint(const Implementation *ep)
{ }
Scalar density() const
{ return density_; }
/*!
* \brief Return the vector of primary variables
* \brief Returns the dynamic viscosity of the fluid within the
* control volume.
*/
const PrimaryVariables &primaryVars() const
{ return primaryVars_; }
Scalar viscosity() const
{ return viscosity_; }
Scalar pressure;
Scalar density;
Scalar viscosity;
Scalar porosity;
/*!
* \brief Returns the average porosity within the control volume.
*/
Scalar porosity() const
{ return porosity_; }
protected:
PrimaryVariables primaryVars_;
void updateTemperature_(const PrimaryVariables &priVars,
const Element &element,
const FVElementGeometry &elemGeom,
int scvIdx,
const Problem &problem)
{ temperature_ = problem.temperature(element, elemGeom, scvIdx); }
Scalar temperature_;
Scalar pressure_;
Scalar density_;
Scalar viscosity_;
Scalar porosity_;
};
}
......
......@@ -37,15 +37,20 @@ class OnePTwoCFluidState : public FluidState<typename GET_PROP_TYPE(TypeTag, PTA
OnePTwoCFluidState<TypeTag> >
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(OnePTwoCIndices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
enum {
konti = Indices::konti,
transport = Indices::transport
pressureIdx = Indices::pressureIdx,
x1Idx = Indices::x1Idx,
contiEqIdx = Indices::contiEqIdx,
transEqIdx = Indices::transEqIdx,
phaseIndex = GET_PROP_VALUE(TypeTag, PTAG(PhaseIndex)),
comp1Index = GET_PROP_VALUE(TypeTag, PTAG(Comp1Index)),
comp2Index = GET_PROP_VALUE(TypeTag, PTAG(Comp2Index)),
};
public:
......@@ -62,26 +67,106 @@ public:
temperature_ = temperature;
phasePressure_[0] = primaryVars[konti];
moleFraction_[0][konti] = 1.0 - primaryVars[transport];
moleFraction_[0][transport] = primaryVars[transport];
phasePressure_ = primaryVars[pressureIdx];
x1_ = primaryVars[x1Idx];
meanMolarMass_ =
(1 - x1_)*FluidSystem::molarMass(comp1Index) +
(x1_ )*FluidSystem::molarMass(comp2Index);
density_ = FluidSystem::phaseDensity(phaseIndex, temperature_, phasePressure_, *this);
Valgrind::CheckDefined(x1_);
Valgrind::CheckDefined(phasePressure_);
Valgrind::CheckDefined(density_);
Valgrind::CheckDefined(meanMolarMass_);
Valgrind::CheckDefined(temperature_);
Valgrind::CheckDefined(*this);
}
public:
/*!
* \brief Returns the saturation of a phase.
*/
Scalar saturation(int phaseIdx) const
{ return (phaseIndex == phaseIdx)?1.0:0.0; };
/*!
* \brief Returns the molar fraction of a component in a fluid phase.
*/
Scalar moleFrac(int phaseIdx, int compIdx) const
{
// we are a single phase model!
if (phaseIdx != phaseIndex) return 0.0;
if (compIdx==comp1Index)
return 1-x1_;
else if (compIdx==comp2Index)
return x1_;
return 0.0;
}
/*!
* \brief Returns the total concentration of a phase [mol / m^3].
*
* This is equivalent to the sum of all component concentrations.
*/
Scalar phaseConcentration(int phaseIdx) const