Commit 2fed48b4 authored by Markus Wolff's avatar Markus Wolff
Browse files

changed decoupled 1p models to new structure



git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@7431 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 1edee536
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* Copyright (C) 2011 by Markus Wolff *
* 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 two-phase box model.
*/
#ifndef DUMUX_DECOUPLED_1P_INDICES_HH
#define DUMUX_DECOUPLED_1P_INDICES_HH
namespace Dumux
{
/*!
* \ingroup IMPES
*/
// \{
/*!
* \brief The common indices for the 1-p models.
*/
struct DecoupledOnePCommonIndices
{
// Formulations
static const int pressEqIdx = 0;
};
// \}
} // namespace Dumux
#endif
......@@ -33,14 +33,10 @@
#ifndef DUMUX_1PPROPERTIES_HH
#define DUMUX_1PPROPERTIES_HH
//Dune-includes
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
//Dumux-includes
#include <dumux/decoupled/common/decoupledproperties.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include "1pindices.hh"
namespace Dumux
{
......@@ -52,6 +48,9 @@ namespace Dumux
template<class TypeTag>
class VariableClass;
template<class TypeTag>
class CellData1P;
////////////////////////////////
// properties
////////////////////////////////
......@@ -90,8 +89,12 @@ SET_INT_PROP(DecoupledOneP, NumPhases, 1)//!< Single phase system
;
SET_INT_PROP(DecoupledOneP, NumComponents, 1); //!< Each phase consists of 1 pure component
SET_TYPE_PROP(DecoupledOneP, Indices, DecoupledOnePCommonIndices);
SET_TYPE_PROP(DecoupledOneP, Variables, VariableClass<TypeTag>);
SET_TYPE_PROP(DecoupledOneP, CellData, CellData1P<TypeTag>);
SET_PROP(DecoupledOneP, PressureCoefficientMatrix)
{
private:
......
/*****************************************************************************
* Copyright (C) 2011 by Markus Wolff *
* 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_ELEMENTDATA1P_HH
#define DUMUX_ELEMENTDATA1P_HH
#include "1pproperties.hh"
#include "fluxData1p.hh"
/**
* @file
* @brief Class including the variables and data of discretized data of the constitutive relations for one element
* @author Markus Wolff
*/
namespace Dumux
{
/*!
* \ingroup IMPES
*/
//! Class including the variables and data of discretized data of the constitutive relations for one element.
/*! The variables of two-phase flow, which are one pressure and one saturation are stored in this class.
* Additionally, a velocity needed in the transport part of the decoupled two-phase flow is stored, as well as discretized data of constitutive relationships like
* mobilities, fractional flow functions and capillary pressure. Thus, they have to be callculated just once in every time step or every iteration step.
*
* @tparam TypeTag The Type Tag
1*/
template<class TypeTag>
class CellData1P
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef FluxData1P<TypeTag> FluxData;
private:
Scalar pressure_;
FluxData fluxData_;
public:
//! Constructs a VariableClass object
/**
* @param gridView a DUNE gridview object corresponding to diffusion and transport equation
*/
CellData1P() :
pressure_(0.0)
{
}
FluxData& fluxData()
{
return fluxData_;
}
const FluxData& fluxData() const
{
return fluxData_;
}
////////////////////////////////////////////////////////////
// functions returning primary variables
////////////////////////////////////////////////////////////
Scalar pressure()
{
return pressure_;
}
Scalar pressure() const
{
return pressure_;
}
void setPressure(Scalar press)
{
pressure_ = press;
}
};
}
#endif
......@@ -28,9 +28,10 @@
#ifndef DUMUX_DIFFUSIONPROBLEM_1P_HH
#define DUMUX_DIFFUSIONPROBLEM_1P_HH
#include <dumux/decoupled/common/onemodelproblem_old.hh>
#include <dumux/decoupled/common/variableclass_old.hh>
#include <dumux/decoupled/common/onemodelproblem.hh>
#include <dumux/decoupled/common/variableclass.hh>
#include <dumux/decoupled/1p/1pproperties.hh>
#include <dumux/decoupled/1p/cellData1p.hh>
namespace Dumux
{
......
......@@ -28,7 +28,6 @@
* @author Markus Wolff
*/
#include <dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh>
namespace Dumux
{
......@@ -46,27 +45,27 @@ namespace Dumux
*/
template<class TypeTag>
class FVVelocity1P: public FVPressure1P<TypeTag>
class FVVelocity1P
{
typedef FVVelocity1P<TypeTag> ThisType;
typedef FVPressure1P<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
typedef typename GET_PROP_TYPE(TypeTag, Variables) Variables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
typedef typename GET_PROP_TYPE(TypeTag, SpatialParameters) SpatialParameters;
typedef typename GET_PROP_TYPE(TypeTag, Fluid) Fluid;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
typedef typename GET_PROP(TypeTag, SolutionTypes) SolutionTypes;
typedef typename SolutionTypes::PrimaryVariables PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, CellData) CellData;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::Grid Grid;
typedef typename GridView::IndexSet IndexSet;
typedef typename GridView::template Codim<0>::Iterator ElementIterator;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename Grid::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::Intersection Intersection;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
enum
{
......@@ -75,7 +74,7 @@ typedef typename GridView::Traits::template Codim<0>::Entity Element;
enum
{
pressEqIdx = 0 // only one equation!
pressEqIdx = Indices::pressEqIdx // only one equation!
};
typedef Dune::FieldVector<Scalar,dimWorld> GlobalPosition;
......@@ -87,8 +86,15 @@ public:
* \param problem a problem class object
*/
FVVelocity1P(Problem& problem)
: FVPressure1P<TypeTag>(problem), problem_(problem)
{ }
: problem_(problem), gravity_(problem.gravity())
{
const Element& element = *(problem_.gridView().template begin<0> ());
Scalar temperature = problem_.temperature(element);
Scalar referencePress = problem_.referencePressure(element);
density_ = Fluid::density(temperature, referencePress);
viscosity_ = Fluid::viscosity(temperature, referencePress);
}
//! Calculate the velocity.
......@@ -97,25 +103,15 @@ public:
* Given the piecewise constant pressure \f$p\f$,
* this method calculates the velocity field
*/
void calculateVelocity();
void initialize()
{
ParentType::initialize();
void calculateVelocity(const Intersection&, CellData&);
calculateVelocity();
return;
}
void calculateVelocityOnBoundary(const Intersection&, CellData&);
//! \brief Write data files
/* \param name file name */
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
ParentType::addOutputVtkFields(writer);
Dune::BlockVector<Dune::FieldVector<Scalar, dim> > &velocity = *(writer.template allocateManagedBuffer<Scalar, dim> (
problem_.gridView().size(0)));
......@@ -126,6 +122,7 @@ public:
// cell index
int globalIdx = problem_.variables().index(*eIt);
CellData& cellData = problem_.variables().cellData(globalIdx);
Dune::FieldVector<Scalar, 2*dim> flux(0);
// run through all intersections with neighbors and boundary
......@@ -137,7 +134,8 @@ public:
{
int isIndex = isIt->indexInInside();
flux[isIndex] = isIt->geometry().volume() * (isIt->centerUnitOuterNormal() * problem_.variables().velocityElementFace(*eIt, isIndex));
flux[isIndex] = isIt->geometry().volume()
* (isIt->centerUnitOuterNormal() * cellData.fluxData().velocity(isIndex));
}
Dune::FieldVector<Scalar, dim> refVelocity(0);
......@@ -167,173 +165,157 @@ public:
}
private:
Problem &problem_;
const GlobalPosition& gravity_; //!< vector including the gravity constant
Scalar density_;
Scalar viscosity_;
};
template<class TypeTag>
void FVVelocity1P<TypeTag>::calculateVelocity()
void FVVelocity1P<TypeTag>::calculateVelocity(const Intersection& intersection, CellData& cellDataI)
{
BoundaryTypes bcType;
// compute update vector
ElementIterator eItEnd = problem_.gridView().template end<0>();
for (ElementIterator eIt = problem_.gridView().template begin<0>(); eIt != eItEnd; ++eIt)
{
// cell index
int globalIdxI = problem_.variables().index(*eIt);
Scalar pressI = problem_.variables().pressure()[globalIdxI];
Scalar temperatureI = problem_.temperature(*eIt);
Scalar referencePressI = problem_.referencePressure(*eIt);
Scalar densityI = Fluid::density(temperatureI, referencePressI);
Scalar viscosityI = Fluid::viscosity(temperatureI, referencePressI);
ElementPointer elementI = intersection.inside();
ElementPointer elementJ = intersection.outside();
// run through all intersections with neighbors and boundary
IntersectionIterator
isItEnd = problem_.gridView().iend(*eIt);
for (IntersectionIterator
isIt = problem_.gridView().ibegin(*eIt); isIt
!=isItEnd; ++isIt)
{
// local number of facet
int isIndex = isIt->indexInInside();
const GlobalPosition& unitOuterNormal = isIt->centerUnitOuterNormal();
int globalIdxJ = problem_.variables().index(*elementJ);
// handle interior face
if (isIt->neighbor())
{
// access neighbor
ElementPointer neighborPointer = isIt->outside();
int globalIdxJ = problem_.variables().index(*neighborPointer);
CellData& cellDataJ = problem_.variables().cellData(globalIdxJ);
// get global coordinates of cell centers
const GlobalPosition& globalPosI = elementI->geometry().center();
const GlobalPosition& globalPosJ = elementJ->geometry().center();
// cell center in global coordinates
const GlobalPosition& globalPos = eIt->geometry().center();
//get face index
int isIndexI = intersection.indexInInside();
int isIndexJ = intersection.indexInOutside();
// neighbor cell center in global coordinates
const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
//get face normal
const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
// distance vector between barycenters
GlobalPosition distVec = globalPosNeighbor - globalPos;
// distance vector between barycenters
GlobalPosition distVec = globalPosJ - globalPosI;
// compute distance between cell centers
Scalar dist = distVec.two_norm();
// compute distance between cell centers
Scalar dist = distVec.two_norm();
// compute vectorized permeabilities
FieldMatrix meanPermeability(0);
// compute vectorized permeabilities
FieldMatrix meanPermeability(0);
problem_.spatialParameters().meanK(meanPermeability,
problem_.spatialParameters().intrinsicPermeability(*eIt),
problem_.spatialParameters().intrinsicPermeability(*neighborPointer));
problem_.spatialParameters().meanK(meanPermeability, problem_.spatialParameters().intrinsicPermeability(*elementI),
problem_.spatialParameters().intrinsicPermeability(*elementJ));
Dune::FieldVector<Scalar, dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
Dune::FieldVector < Scalar, dim > permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
permeability/=viscosityI;
permeability /= viscosity_;
// std::cout<<"Mean Permeability / Viscosity: "<<meanPermeability<<std::endl;
//calculate potential gradients
Scalar potential = (cellDataI.pressure() - cellDataJ.pressure()) / dist;
Scalar pressJ = problem_.variables().pressure()[globalIdxJ];
potential += density_ * (unitOuterNormal * gravity_);
Scalar temperatureJ = problem_.temperature(*eIt);
Scalar referencePressJ = problem_.referencePressure(*eIt);
//store potentials for further calculations (velocity, saturation, ...)
cellDataI.fluxData().setPotential(isIndexI, potential);
cellDataJ.fluxData().setPotential(isIndexJ, -potential);
Scalar densityJ = Fluid::density(temperatureJ, referencePressJ);
//calculate the gravity term
GlobalPosition velocity(permeability);
velocity *= (cellDataI.pressure() - cellDataJ.pressure()) / dist;
//calculate potential gradients
Scalar potential = problem_.variables().potential(globalIdxI, isIndex);
GlobalPosition gravityTerm(unitOuterNormal);
gravityTerm *= (gravity_ * permeability) * density_;
Scalar density = (potential> 0.) ? densityI : densityJ;
velocity += gravityTerm;
density= (potential == 0.) ? 0.5 * (densityI + densityJ) : density;
//store velocities
cellDataI.fluxData().setVelocity(isIndexI, velocity);
cellDataI.fluxData().setVelocityMarker(isIndexI);
potential = (pressI - pressJ) / dist;
cellDataJ.fluxData().setVelocity(isIndexJ, velocity);
cellDataJ.fluxData().setVelocityMarker(isIndexJ);
return;
}
potential += density * (unitOuterNormal * this->gravity);
template<class TypeTag>
void FVVelocity1P<TypeTag>::calculateVelocityOnBoundary(const Intersection& intersection, CellData& cellData)
{
ElementPointer element = intersection.inside();
//store potentials for further calculations (velocity, saturation, ...)
problem_.variables().potential(globalIdxI, isIndex) = potential;
//get face index
int isIndex = intersection.indexInInside();
//do the upwinding depending on the potentials
density = (potential> 0.) ? densityI : densityJ;
density = (potential == 0.) ? 0.5 * (densityI + densityJ) : density;
//get face normal
const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
//calculate the gravity term
GlobalPosition velocity(permeability);
velocity *= (pressI - pressJ)/dist;
BoundaryTypes bcType;
//get boundary type
problem_.boundaryTypes(bcType, intersection);
PrimaryVariables boundValues(0.0);
GlobalPosition gravityTerm(unitOuterNormal);
gravityTerm *= (this->gravity*permeability)*density;
if (bcType.isDirichlet(pressEqIdx))
{
problem_.dirichlet(boundValues, intersection);
//store velocities
problem_.variables().velocity()[globalIdxI][isIndex] = (velocity + gravityTerm);
// get global coordinates of cell centers
const GlobalPosition& globalPosI = element->geometry().center();
}//end intersection with neighbor
// center of face in global coordinates
const GlobalPosition& globalPosJ = intersection.geometry().center();
// handle boundary face
else if (isIt->boundary())
{
// center of face in global coordinates
const GlobalPosition& globalPosFace = isIt->geometry().center();
// distance vector between barycenters
GlobalPosition distVec = globalPosJ - globalPosI;
//get boundary type
problem_.boundaryTypes(bcType, *isIt);
PrimaryVariables boundValues(0.0);
// compute distance between cell centers
Scalar dist = distVec.two_norm();
if (bcType.isDirichlet(pressEqIdx))
{
problem_.dirichlet(boundValues, *isIt);
//permeability vector at boundary
// compute vectorized permeabilities
FieldMatrix meanPermeability(0);
// cell center in global coordinates
const GlobalPosition& globalPos = eIt->geometry().center();
problem_.spatialParameters().meanK(meanPermeability, problem_.spatialParameters().intrinsicPermeability(*element));
// distance vector between barycenters
GlobalPosition distVec = globalPosFace - globalPos;
//multiply with normal vector at the boundary
Dune::FieldVector < Scalar, dim > permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
permeability /= viscosity_;
// compute distance between cell centers
Scalar dist = distVec.two_norm();
Scalar pressBound = boundValues;
//permeability vector at boundary
// compute vectorized permeabilities
FieldMatrix meanPermeability(0);
//calculate potential gradients
Scalar potential = (cellData.pressure() - pressBound) / dist;
problem_.spatialParameters().meanK(meanPermeability,
problem_.spatialParameters().intrinsicPermeability(*eIt));
potential += density_ * (unitOuterNormal * gravity_);
//multiply with normal vector at the boundary
Dune::FieldVector<Scalar,dim> permeability(0);
meanPermeability.mv(unitOuterNormal, permeability);
permeability/=viscosityI;
//store potentials for further calculations (velocity, saturation, ...)
cellData.fluxData().setPotential(isIndex, potential);
Scalar pressBound = boundValues;
//calculate the gravity term
GlobalPosition velocity(permeability);
velocity *= (cellData.pressure() - pressBound) / dist;
//calculate the gravity term
GlobalPosition velocity(permeability);
velocity *= (pressI - pressBound)/dist;
GlobalPosition gravityTerm(unitOuterNormal);
gravityTerm *= (gravity_ * permeability) * density_;
GlobalPosition gravityTerm(unitOuterNormal);
gravityTerm *= (this->gravity*permeability)*densityI;
velocity += gravityTerm;
problem_.variables().velocity()[globalIdxI][isIndex] = (velocity + gravityTerm);
//store velocities
cellData.fluxData().setVelocity(isIndex, velocity);
cellData.fluxData().setVelocityMarker(isIndex);
}//end dirichlet boundary
} //end dirichlet boundary
else
{
problem_.neumann(boundValues, *isIt);
GlobalPosition velocity(unitOuterNormal);
else
{
problem_.neumann(boundValues, intersection);
GlobalPosition velocity(unitOuterNormal);
velocity *= boundValues[pressEqIdx]/densityI;
velocity *= boundValues[pressEqIdx] / density_;
problem_.variables().velocity()[globalIdxI][isIndex] = velocity;
//store potential gradients for further calculations
cellData.fluxData().setPotential(isIndex, boundValues[pressEqIdx]);
}//end neumann boundary
}//end boundary
}// end all intersections
}// end grid traversal
// printvector(std::cout, problem_.variables().velocity(), "velocity", "row", 4, 1, 3);
//store velocity
cellData.fluxData().setVelocity(isIndex, velocity);
cellData.fluxData().setVelocityMarker(isIndex);
} //end neumann boundary
return;
}