opm-simulators/examples/tutorialproblem_decoupled.hh
Andreas Lauser 15116fb9f8 decoupled tutorial: streamline the code a bit
we do not need to set any default arguments for the type tag in the
problem definition. the type tag which is actually used is chosen in
the main function.
2012-07-12 21:22:53 +02:00

294 lines
12 KiB
C++

// $Id$
/*****************************************************************************
* Copyright (C) 2007-2008 by Klaus Mosthaf *
* Copyright (C) 2007-2008 by Bernd Flemisch *
* Copyright (C) 2008-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 problem for the sequential tutorial
*/
#ifndef DUMUX_TUTORIALPROBLEM_DECOUPLED_HH
#define DUMUX_TUTORIALPROBLEM_DECOUPLED_HH
// the grid includes
#include <dune/grid/sgrid.hh>
// dumux 2p-decoupled environment
#include <dumux/decoupled/2p/impes/impesproblem2p.hh> /*@\label{tutorial-decoupled:parent-problem}@*/
#include <dumux/decoupled/2p/diffusion/fv/fvvelocity2p.hh>
#include <dumux/decoupled/2p/transport/fv/fvsaturation2p.hh>
#include <dumux/decoupled/2p/transport/fv/capillarydiffusion.hh>
// assign parameters dependent on space (e.g. spatial parameters)
#include "tutorialspatialparameters_decoupled.hh" /*@\label{tutorial-decoupled:spatialparameters}@*/
namespace Dumux
{
template<class TypeTag>
class TutorialProblemDecoupled;
//////////
// Specify the properties for the lens problem
//////////
namespace Properties
{
// create a new type tag for the problem
NEW_TYPE_TAG(TutorialProblemDecoupled, INHERITS_FROM(DecoupledTwoP)); /*@\label{tutorial-decoupled:create-type-tag}@*/
// Set the problem property
SET_PROP(TutorialProblemDecoupled, Problem) /*@\label{tutorial-decoupled:set-problem}@*/
{
typedef Dumux::TutorialProblemDecoupled<TypeTag> type;
};
// Set the grid type
SET_PROP(TutorialProblemDecoupled, Grid) /*@\label{tutorial-decoupled:grid-begin}@*/
{
typedef Dune::SGrid<2, 2> type; /*@\label{tutorial-decoupled:set-grid-type}@*/
static type *create() /*@\label{tutorial-decoupled:create-grid-method}@*/
{
typedef typename type::ctype ctype;
Dune::FieldVector<int, 2> cellRes;
Dune::FieldVector<ctype, 2> lowerLeft(0.0);
Dune::FieldVector<ctype, 2> upperRight;
upperRight[0] = 300;
upperRight[1] = 60;
cellRes[0] = 100;
cellRes[1] = 1;
return new Dune::SGrid<2,2>(cellRes,
lowerLeft,
upperRight);
} /*@\label{tutorial-decoupled:grid-end}@*/
};
// Set the wetting phase
SET_PROP(TutorialProblemDecoupled, WettingPhase) /*@\label{tutorial-decoupled:2p-system-start}@*/
{
private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
public:
typedef Dumux::LiquidPhase<Scalar, Dumux::H2O<Scalar> > type; /*@\label{tutorial-decoupled:wettingPhase}@*/
};
// Set the non-wetting phase
SET_PROP(TutorialProblemDecoupled, NonwettingPhase)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
public:
typedef Dumux::LiquidPhase<Scalar, Dumux::Oil<Scalar> > type; /*@\label{tutorial-decoupled:nonwettingPhase}@*/
}; /*@\label{tutorial-decoupled:2p-system-end}@*/
// Set the spatial parameters
SET_PROP(TutorialProblemDecoupled, SpatialParameters) /*@\label{tutorial-decoupled:set-spatialparameters}@*/
{
typedef Dumux::TutorialSpatialParametersDecoupled<TypeTag> type;
};
// Set the model properties
SET_PROP(TutorialProblemDecoupled, TransportModel) /*@\label{tutorial-decoupled:TransportModel}@*/
{
typedef Dumux::FVSaturation2P<TTAG(TutorialProblemDecoupled)> type;
};
SET_PROP(TutorialProblemDecoupled, PressureModel) /*@\label{tutorial-decoupled:PressureModel}@*/
{
typedef Dumux::FVVelocity2P<TTAG(TutorialProblemDecoupled)> type;
};
// model-specific settings
SET_INT_PROP(TutorialProblemDecoupled, VelocityFormulation,
GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices))::velocityW); /*@\label{tutorial-decoupled:velocityFormulation}@*/
SET_TYPE_PROP(TutorialProblemDecoupled, DiffusivePart,
Dumux::CapillaryDiffusion<TypeTag>); /*@\label{tutorial-decoupled:DiffusivePart}@*/
SET_SCALAR_PROP(TutorialProblemDecoupled, CFLFactor, 0.5); /*@\label{tutorial-decoupled:cfl}@*/
// Disable gravity
SET_BOOL_PROP(TutorialProblemDecoupled, EnableGravity, false); /*@\label{tutorial-decoupled:gravity}@*/
} /*@\label{tutorial-decoupled:propertysystem-end}@*/
/*! \ingroup DecoupledProblems
* @brief Problem class for the decoupled tutorial
*/
template<class TypeTag>
class TutorialProblemDecoupled: public IMPESProblem2P<TypeTag, TutorialProblemDecoupled<TypeTag> > /*@\label{tutorial-decoupled:def-problem}@*/
{
typedef TutorialProblemDecoupled<TypeTag> ThisType;
typedef IMPESProblem2P<TypeTag, ThisType> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(TwoPIndices)) Indices;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidState)) FluidState;
enum
{
dim = GridView::dimension, dimWorld = GridView::dimensionworld
};
enum
{
wPhaseIdx = Indices::wPhaseIdx, nPhaseIdx = Indices::nPhaseIdx
};
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GridView::Traits::template Codim<0>::Entity Element;
typedef typename GridView::Intersection Intersection;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldVector<Scalar, dim> LocalPosition;
public:
TutorialProblemDecoupled(const GridView &gridView, const GlobalPosition lowerLeft = GlobalPosition(0.),
const GlobalPosition upperRight = GlobalPosition(0.)) : ParentType(gridView) /*@\label{tutorial-decoupled:constructor-problem}@*/
{ }
//! The problem name.
/*! This is used as a prefix for files generated by the simulation.
*/
const char *name() const /*@\label{tutorial-decoupled:name}@*/
{
return "tutorial_decoupled";
}
//! Returns true if a restart file should be written.
/* The default behaviour is to write no restart file.
*/
bool shouldWriteRestartFile() const /*@\label{tutorial-decoupled:restart}@*/
{
return false;
}
//! Returns true if the current solution should be written to disk (i.e. as a VTK file)
/*! The default behaviour is to write out every the solution for
* very time step. Else, change divisor.
*/
bool shouldWriteOutput() const /*@\label{tutorial-decoupled:output}@*/
{
return this->timeManager().timeStepIndex() > 0 &&
(this->timeManager().timeStepIndex() % 1 == 0);
}
//! Returns the temperature within the domain.
/*! This problem assumes a temperature of 10 degrees Celsius.
*/
Scalar temperature(const GlobalPosition& globalPos, const Element& element) const /*@\label{tutorial-decoupled:temperature}@*/
{
return 273.15 + 10; // -> 10°C
}
//! Returns a constant pressure to enter material laws
/* For incrompressible simulations, a constant pressure is necessary
* to enter the material laws to gain a constant density etc. In the compressible
* case, the pressure is used for the initialization of material laws.
*/
Scalar referencePressure(const GlobalPosition& globalPos, const Element& element) const /*@\label{tutorial-decoupled:refPressure}@*/
{
return 2e5;
}
//! Source of mass \f$ [\frac{kg}{m^3 \cdot s}] \f$
/*! Evaluate the source term for all phases within a given
* volume. The method returns the mass generated (positive) or
* annihilated (negative) per volume unit.
*/
std::vector<Scalar> source(const GlobalPosition& globalPos, const Element& element) /*@\label{tutorial-decoupled:source}@*/
{
return std::vector<Scalar>(2, 0.);
}
//! Type of pressure boundary condition.
/*! Defines the type the boundary condition for the pressure equation,
* either pressure (dirichlet) or flux (neumann).
*/
typename BoundaryConditions::Flags bctypePress(const GlobalPosition& globalPos, const Intersection& intersection) const /*@\label{tutorial-decoupled:bctypePress}@*/
{
if (globalPos[0] < this->bboxMin()[0] + eps_)
return BoundaryConditions::dirichlet;
// all other boundaries
return BoundaryConditions::neumann;
}
//! Type of Transport boundary condition.
/*! Defines the type the boundary condition for the transport equation,
* either saturation (dirichlet) or flux (neumann).
*/
BoundaryConditions::Flags bctypeSat(const GlobalPosition& globalPos, const Intersection& intersection) const /*@\label{tutorial-decoupled:bctypeSat}@*/
{
if (globalPos[0] < this->bboxMin()[0] + eps_)
return Dumux::BoundaryConditions::dirichlet;
else
return Dumux::BoundaryConditions::neumann;
}
//! Value for dirichlet pressure boundary condition \f$ [Pa] \f$.
/*! In case of a dirichlet BC for the pressure equation, the pressure
* have to be defined on boundaries.
*/
Scalar dirichletPress(const GlobalPosition& globalPos, const Intersection& intersection) const /*@\label{tutorial-decoupled:dirichletPress}@*/
{
return 2e5;
}
//! Value for transport dirichlet boundary condition (dimensionless).
/*! In case of a dirichlet BC for the transport equation, a saturation
* has to be defined on boundaries.
*/
Scalar dirichletSat(const GlobalPosition& globalPos, const Intersection& intersection) const /*@\label{tutorial-decoupled:dirichletSat}@*/
{
return 1;
}
//! Value for pressure neumann boundary condition \f$ [\frac{kg}{m^3 \cdot s}] \f$.
/*! In case of a neumann boundary condition, the flux of matter
* is returned as a vector.
*/
std::vector<Scalar> neumannPress(const GlobalPosition& globalPos, const Intersection& intersection) const /*@\label{tutorial-decoupled:neumannPress}@*/
{
std::vector<Scalar> neumannFlux(2,0.0);
if (globalPos[0] > this->bboxMax()[0] - eps_)
{
neumannFlux[nPhaseIdx] = 3e-4;
}
return neumannFlux;
}
//! Value for transport neumann boundary condition \f$ [\frac{kg}{m^3 \cdot s}] \f$.
/*! In case of a neumann boundary condition for the transport equation
* the flux of matter for the primary variable is returned as a scalar.
*/
Scalar neumannSat(const GlobalPosition& globalPos, const Intersection& intersection, Scalar factor) const /*@\label{tutorial-decoupled:neumannSat}@*/
{
return 0;
}
//! Saturation initial condition (dimensionless)
/*! The problem is initialized with the following saturation.
*/
Scalar initSat(const GlobalPosition& globalPos, const Element& element) const /*@\label{tutorial-decoupled:initSat}@*/
{
return 0;
}
private:
GlobalPosition lowerLeft_;
GlobalPosition upperRight_;
static const Scalar eps_ = 1e-6;
};
} //end namespace
#endif