Commit 4e59dedf authored by Klaus Mosthaf's avatar Klaus Mosthaf
Browse files

added folder freeflow plus the stokes* header files from the devel

branch


git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7445 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent e57016a2
/*****************************************************************************
* Copyright (C) 2010 by Katherina Baber, Klaus Mosthaf *
* Copyright (C) 2008-2009 by Bernd Flemisch, 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. *
* *
* 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/>. *
*****************************************************************************/
/*!
* \file
*
* \brief This file contains the data which is required to calculate
* all fluxes of the fluid phase over a face of a finite volume.
*
* This means pressure and temperature gradients, phase densities at
* the integration point, etc.
*/
#ifndef DUMUX_STOKES_FLUX_VARIABLES_HH
#define DUMUX_STOKES_FLUX_VARIABLES_HH
#include <dumux/common/math.hh>
namespace Dumux
{
/*!
* \ingroup StokesModel
* \brief This template class contains the data which is required to
* calculate the fluxes of the fluid phases over a face of a
* finite volume for the stokes model.
*
* This means pressure and concentration gradients, phase densities at
* the intergration point, etc.
*/
template <class TypeTag>
class StokesFluxVariables
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
enum { dim = GridView::dimension };
typedef typename GridView::template Codim<0>::Entity Element;
typedef Dune::FieldVector<Scalar, dim> FieldVector;
typedef Dune::FieldVector<Scalar, dim> VelocityVector;
typedef Dune::FieldVector<Scalar, dim> ScalarGradient;
typedef Dune::FieldMatrix<Scalar, dim, dim> VectorGradient;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
public:
StokesFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
int faceIdx,
const ElementVolumeVariables &elemVolVars,
bool onBoundary = false)
: fvGeom_(elemGeom), onBoundary_(onBoundary), faceIdx_(faceIdx)
{
calculateValues_(problem, element, elemVolVars);
determineUpwindDirection_(elemVolVars);
};
protected:
void calculateValues_(const Problem &problem,
const Element &element,
const ElementVolumeVariables &elemVolVars)
{
// calculate gradients and secondary variables at IPs
FieldVector tmp(0.0);
densityAtIP_ = Scalar(0);
viscosityAtIP_ = Scalar(0);
pressureAtIP_ = Scalar(0);
normalVelocityAtIP_ = Scalar(0);
velocityAtIP_ = Scalar(0);
pressureGradAtIP_ = Scalar(0);
velocityGradAtIP_ = Scalar(0);
velocityDivAtIP_ = Scalar(0);
for (int idx = 0;
idx < fvGeom_.numVertices;
idx++) // loop over adjacent vertices
{
// phase density and viscosity at IP
densityAtIP_ += elemVolVars[idx].density() *
face().shapeValue[idx];
viscosityAtIP_ += elemVolVars[idx].viscosity() *
face().shapeValue[idx];
pressureAtIP_ += elemVolVars[idx].pressure() *
face().shapeValue[idx];
// velocity at the IP (fluxes)
VelocityVector velocityTimesShapeValue = elemVolVars[idx].velocity();
velocityTimesShapeValue *= face().shapeValue[idx];
velocityAtIP_ += velocityTimesShapeValue;
// the pressure gradient
tmp = face().grad[idx];
tmp *= elemVolVars[idx].pressure();
pressureGradAtIP_ += tmp;
// take gravity into account
tmp = problem.gravity();
tmp *= densityAtIP_;
// pressure gradient including influence of gravity
pressureGradAtIP_ -= tmp;
// the velocity gradients
for (int dimIdx = 0; dimIdx<dim; ++dimIdx)
{
tmp = face().grad[idx];
tmp *= elemVolVars[idx].velocity()[dimIdx];
velocityGradAtIP_[dimIdx] += tmp;
velocityDivAtIP_ += face().grad[idx][dimIdx]*elemVolVars[idx].velocity()[dimIdx];
}
}
normalVelocityAtIP_ = velocityAtIP_ * face().normal;
Valgrind::CheckDefined(densityAtIP_);
Valgrind::CheckDefined(viscosityAtIP_);
Valgrind::CheckDefined(normalVelocityAtIP_);
Valgrind::CheckDefined(velocityAtIP_);
Valgrind::CheckDefined(pressureGradAtIP_);
Valgrind::CheckDefined(velocityGradAtIP_);
Valgrind::CheckDefined(velocityDivAtIP_);
};
void determineUpwindDirection_(const ElementVolumeVariables &elemVolVars)
{
// set the upstream and downstream vertices
upstreamIdx_ = face().i;
downstreamIdx_ = face().j;
if (normalVelocityAtIP() < 0)
std::swap(upstreamIdx_, downstreamIdx_);
};
public:
/*!
* \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 fvGeom_.boundaryFace[faceIdx_];
else
return fvGeom_.subContVolFace[faceIdx_];
}
/*!
* \brief Return the average volume of the upstream and the downstream sub-control volume;
* this is required for the stabilization
*/
const Scalar averageSCVVolume() const
{
return 0.5*(fvGeom_.subContVol[upstreamIdx_].volume +
fvGeom_.subContVol[downstreamIdx_].volume);
}
/*!
* \brief Return pressure \f$\mathrm{[Pa]}\f$ at the integration
* point.
*/
Scalar pressureAtIP() const
{ return pressureAtIP_; }
/*!
* \brief Return density \f$\mathrm{[kg/m^3]}\f$ at the integration
* point.
*/
Scalar densityAtIP() const
{ return densityAtIP_; }
/*!
* \brief Return viscosity \f$\mathrm{[m^2/s]}\f$ at the integration
* point.
*/
Scalar viscosityAtIP() const
{ return viscosityAtIP_; }
/*!
* \brief Return the normal velocity \f$\mathrm{[m/s]}\f$ at the integration
* point.
*/
Scalar normalVelocityAtIP() const
{ return normalVelocityAtIP_; }
/*!
* \brief Return the pressure gradient at the integration
* point.
*/
const ScalarGradient &pressureGradAtIP() const
{ return pressureGradAtIP_; }
/*!
* \brief Return the velocity at the integration
* point.
*/
const VelocityVector &velocityAtIP() const
{ return velocityAtIP_; }
/*!
* \brief Return the velocity gradient at the integration
* point.
*/
const VectorGradient &velocityGradAtIP() const
{ return velocityGradAtIP_; }
/*!
* \brief Return the divergence of the normal velocity at the integration
* point.
*/
Scalar velocityDivAtIP() const
{ return velocityDivAtIP_; }
/*!
* \brief Return the local index of the upstream control volume
* for a given phase.
*/
int upstreamIdx() const
{ return upstreamIdx_; }
/*!
* \brief Return the local index of the downstream control volume
* for a given phase.
*/
int downstreamIdx() const
{ return downstreamIdx_; }
bool onBoundary() const
{ return onBoundary_; }
protected:
const FVElementGeometry &fvGeom_;
const bool onBoundary_;
// values at the integration point
Scalar densityAtIP_;
Scalar viscosityAtIP_;
Scalar pressureAtIP_;
Scalar normalVelocityAtIP_;
Scalar velocityDivAtIP_;
VelocityVector velocityAtIP_;
// gradients at the IPs
ScalarGradient pressureGradAtIP_;
VectorGradient velocityGradAtIP_;
// local index of the upwind vertex
int upstreamIdx_;
// local index of the downwind vertex
int downstreamIdx_;
// the index of the considered face
int faceIdx_;
};
} // end namepace
#endif
/*****************************************************************************
* Copyright (C) 2011 by Katherina Baber, Klaus Mosthaf *
* Copyright (C) 2008-2009 by Bernd Flemisch, 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. *
* *
* 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/>. *
*****************************************************************************/
/*!
* \file
*
* \brief Defines the indices required for the Stokes BOX model.
*/
#ifndef DUMUX_STOKES_INDICES_HH
#define DUMUX_STOKES_INDICES_HH
namespace Dumux
{
// \{
/*!
* \brief The common indices for the isothermal stokes model.
*
* \tparam PVOffset The first index in a primary variable vector.
*/
template <class TypeTag, int PVOffset = 0>
struct StokesCommonIndices
{
// number of dimensions
typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
static const int dim = Grid::dimension;
// Primary variable indices
static const int momentumXIdx = PVOffset + 0; //!< Index of the x-component of the momentum equation
static const int momentumYIdx = PVOffset + 1; //!< Index of the y-component of the momentum equation
static const int momentumZIdx = PVOffset + 2; //!< Index of the z-component of the momentum equation
static const int lastMomentumIdx = momentumXIdx+dim-1; //!< Index of the last component of the momentum equation
static const int dimXIdx = 0;
static const int dimYIdx = 1;
static const int dimZIdx = 2;
static const int massBalanceIdx = dim; //!< Index of the pressure in a solution vector
static const int pressureIdx = massBalanceIdx;
static const int velocityXIdx = momentumXIdx;
static const int velocityYIdx = momentumYIdx;
static const int velocityZIdx = momentumZIdx;
static const int phaseIdx = 0; //!< Index of the phase (required to use the same fluid system in coupled models)
};
} // end namespace
#endif
/*****************************************************************************
* Copyright (C) 2007 by Peter Bastian *
* Institute of Parallel and Distributed System *
* Department Simulation of Large Systems *
* University of Stuttgart, Germany *
* *
* Copyright (C) 2008-2009 by Klaus Mosthaf *
* Copyright (C) 2007-2009 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. *
* *
* 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/>. *
*****************************************************************************/
/*!
* \file
* \brief Caculates the jacobian of models based on the box scheme element-wise.
*/
#ifndef DUMUX_STOKES_LOCAL_JACOBIAN_HH
#define DUMUX_STOKES_LOCAL_JACOBIAN_HH
#include <dune/istl/matrix.hh>
//#include "boxelementboundarytypes.hh"
namespace Dumux
{
/*!
* \ingroup StokesModel
* \brief Element-wise calculation of the jacobian matrix for models
* based on the box scheme .
*
* \todo Please doc me more!
*/
template<class TypeTag>
class StokesLocalJacobian : public BoxLocalJacobian<TypeTag>
{
private:
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
public:
/*!
* \brief Returns the epsilon value which is added and removed
* from the current solution.
*
* \param elemSol The current solution on the element
* \param scvIdx The local index of the element's vertex for
* which the local derivative ought to be calculated.
* \param pvIdx The index of the primary variable which gets varied
*/
Scalar numericEpsilon_(int scvIdx,
int pvIdx) const
{
Scalar pv = this->curVolVars_[scvIdx].primaryVars()[pvIdx];
if (pvIdx == 0 || pvIdx == 1){
// std::cout << "adjusting eps_ for momentum balance\n";
return 1e-7*(std::abs(pv) + 1);
}
return 1e-9*(std::abs(pv) + 1);
}
};
}
#endif
This diff is collapsed.
/*****************************************************************************
* Copyright (C) 2007 by Peter Bastian *
* Institute of Parallel and Distributed System *
* Department Simulation of Large Systems *
* University of Stuttgart, Germany *
* *
* Copyright (C) 2010 by Katherina Baber, Klaus Mosthaf *
* Copyright (C) 2008-2009 by Bernd Flemisch, 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. *
* *
* 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_STOKES_MODEL_HH
#define DUMUX_STOKES_MODEL_HH
#include "stokeslocalresidual.hh"
#include "stokeslocaljacobian.hh"
#include "stokesproblem.hh"
#include "stokesproperties.hh"
namespace Dumux
{
/*!
* \ingroup BoxProblems
* \defgroup StokesBoxProblems Stokes box problems
*/
/*!
* \ingroup BoxModels
* \defgroup StokesModel Stokes box model
*/
/*!
* \ingroup StokesModel
* \brief Adaption of the BOX scheme to the stokes flow model.
*
* This model implements stokes flow of a single fluid
* discretized by a fully-coupled vertex
* centered finite volume (box) scheme as spatial and
* the implicit Euler method as time discretization.
*/
template<class TypeTag >
class StokesModel : public BoxModel<TypeTag>
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
enum {
dim = GridView::dimension,
dimWorld = GridView::dimensionworld,
numEq = GET_PROP_VALUE(TypeTag, NumEq),
};
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, DofMapper) DofMapper;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
public:
/*!
* \brief Calculate the fluxes across a certain layer in the domain.
* The layer is situated perpendicular to the coordinate axis "coord" and cuts
* the axis at the value "coordValue"
*
*/
void calculateFluxAcrossLayer (const SolutionVector &globalSol, Dune::FieldVector<Scalar, 2> &flux, int coord, Scalar coordVal)
{
GlobalPosition globalI, globalJ;
PrimaryVariables tmpFlux(0.0);
int sign;
FVElementGeometry fvElemGeom;
ElementVolumeVariables elemVolVars;
// Loop over elements
ElementIterator elemIt = this->problem_.gridView().template begin<0>();
ElementIterator endit = this->problem_.gridView().template end<0>();
for (; elemIt != endit; ++elemIt)
{
if (elemIt->partitionType() != Dune::InteriorEntity)
continue;
fvElemGeom.update(this->gridView_(), *elemIt);
elemVolVars.update(this->problem_(), *elemIt, fvElemGeom);
this->localResidual().evalFluxes(*elemIt, elemVolVars);
bool hasLeft = false;
bool hasRight = false;
for (int i = 0; i < fvElemGeom.numVertices; i++) {
const GlobalPosition &global = fvElemGeom.subContVol[i].global;
if (globalI[coord] < coordVal)
hasLeft = true;
else if (globalI[coord] >= coordVal)
hasRight = true;
}
if (!hasLeft || !hasRight)
continue;
for (int i = 0; i < fvElemGeom.numVertices; i++) {
const int &globalIdx =
this->vertexMapper().map(*elemIt, i, dim);
const GlobalPosition &global = fvElemGeom.subContVol[i].global;
if (globalI[coord] < coordVal)
flux += this->localResidual().residual(i);
}
}
flux = this->problem_.gridView().comm().sum(flux);
}
/*!
* \brief Calculate mass in the whole model domain
* and get minimum and maximum values of primary variables
*
*/
void calculateMass(const SolutionVector &sol, Dune::FieldVector<Scalar, 2> &mass)
{
// const DofMapper &dofMapper = this->dofEntityMapper();
VolumeVariables tmp;
Scalar vol, poro, rhoN, rhoW, satN, satW, pW;//, Te;
Scalar massNPhase(0.), massWPhase(0.);
mass = 0;
Scalar minSat = 1e100;
Scalar maxSat = -1e100;
Scalar minP = 1e100;
Scalar maxP = -1e100;
// Scalar minTe = 1e100;
// Scalar maxTe = -1e100;
FVElementGeometry fvElemGeom;
VolumeVariables volVars;
ElementBoundaryTypes elemBcTypes;
// Loop over elements
ElementIterator elemIt = this->problem_.gridView().template begin<0>();
ElementIterator endit = this->problem_.gridView().template end<0>();
for (; elemIt != endit; ++elemIt)
{
if (elemIt->partitionType() != Dune::InteriorEntity)
continue;
fvElemGeom.update(this->gridView_(), *elemIt);
elemBcTypes.update(this->problem_(), *elemIt, fvElemGeom);
int numVerts = elemIt->template count<dim>();
for (int i = 0; i < numVerts; ++i)