diff --git a/dumux/decoupled/1p/1pindices.hh b/dumux/decoupled/1p/1pindices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7541b7229b6088a453bfb7a543daf8c00d19d1fa
--- /dev/null
+++ b/dumux/decoupled/1p/1pindices.hh
@@ -0,0 +1,51 @@
+// -*- 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
diff --git a/dumux/decoupled/1p/1pproperties.hh b/dumux/decoupled/1p/1pproperties.hh
index 8bc1411ff5c4414582a0474f46c30eb357744ddd..8a506b8f2dae97a7b37c5da74a4fd780f85f5df2 100644
--- a/dumux/decoupled/1p/1pproperties.hh
+++ b/dumux/decoupled/1p/1pproperties.hh
@@ -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:
diff --git a/dumux/decoupled/1p/cellData1p.hh b/dumux/decoupled/1p/cellData1p.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9b02de0f9fcf38641c50e70ccd98ed7be2a40081
--- /dev/null
+++ b/dumux/decoupled/1p/cellData1p.hh
@@ -0,0 +1,99 @@
+/*****************************************************************************
+ * 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
diff --git a/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh b/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh
index 0c948880e1c3aae378f8b75e9728b9c3547ab915..743915fe53c9404ae8f4f6c876015fe793c23c3f 100644
--- a/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh
+++ b/dumux/decoupled/1p/diffusion/diffusionproblem1p.hh
@@ -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
{
diff --git a/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh b/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh
index e43e2928a06b6528ea1b2d5f8c63919715b8cf43..9358efcadd68728db70aa5c32c065a76fad8a5ac 100644
--- a/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh
+++ b/dumux/decoupled/1p/diffusion/fv/fvpressure1p.hh
@@ -25,6 +25,7 @@
// dumux environment
+#include <dumux/decoupled/common/fv/fvpressure.hh>
#include <dumux/decoupled/1p/1pproperties.hh>
/**
@@ -51,12 +52,15 @@ namespace Dumux
* @tparam TypeTag The Type Tag
*
*/
-template<class TypeTag> class FVPressure1P
+template<class TypeTag> class FVPressure1P: public FVPressure<TypeTag>
{
+ typedef FVPressure<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;
@@ -65,6 +69,8 @@ template<class TypeTag> class FVPressure1P
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 SolutionTypes::ScalarSolution ScalarSolutionType;
enum
{
@@ -73,7 +79,12 @@ template<class TypeTag> class FVPressure1P
enum
{
- pressEqIdx = 0 // only one equation!
+ pressEqIdx = Indices::pressEqIdx // only one equation!
+ };
+
+ enum
+ {
+ rhs = ParentType::rhs, matrix = ParentType::matrix,
};
typedef typename GridView::Traits::template Codim<0>::Entity Element;
@@ -81,35 +92,25 @@ template<class TypeTag> class FVPressure1P
typedef typename GridView::Grid Grid;
typedef typename GridView::template Codim<0>::EntityPointer ElementPointer;
typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GridView::Intersection Intersection;
typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
typedef Dune::FieldMatrix<Scalar, dim, dim> FieldMatrix;
- typedef typename GET_PROP_TYPE(TypeTag, PressureCoefficientMatrix) Matrix;
- typedef typename GET_PROP_TYPE(TypeTag, PressureRHSVector) Vector;
- //initializes the matrix to store the system of equations
- void initializeMatrix();
-
- //function which assembles the system of equations to be solved
- void assemble(bool first);
-
- //solves the system of equations to get the spatial distribution of the pressure
- void solve();
+public:
+ void getSource(Dune::FieldVector<Scalar, 2>&, const Element&, const CellData&, const bool);
-protected:
- //! Returns reference to the instance of the problem definition
- Problem& problem()
- {
- return problem_;
- }
- //! Returns reference to the instance of the problem definition
- const Problem& problem() const
+ void getStorage(Dune::FieldVector<Scalar, 2>& entries, const Element& element, const CellData& cellData, const bool first)
{
- return problem_;
+ entries = 0;
}
-public:
+ void getFlux(Dune::FieldVector<Scalar, 2>&, const Intersection&, const CellData&, const bool);
+
+ void getFluxOnBoundary(Dune::FieldVector<Scalar, 2>&,
+ const Intersection&, const CellData&, const bool);
+
//! Initializes the problem
/*!
* @param solveTwice repeats the pressure calculation step
@@ -121,45 +122,39 @@ public:
void initialize(bool solveTwice = true)
{
- assemble(true);
- solve();
+ ParentType::initialize();
+ this->assemble(true);
+ this->solve();
if (solveTwice)
{
- assemble(false);
- solve();
+ this->assemble(false);
+ this-> solve();
}
+ storePressureSolution();
return;
}
- //! Calculates the pressure.
- /*!
- * @param solveTwice without any function here!
- *
- * Calculates the pressure \f$p\f$ as solution of the boundary value
- * \f[ \text{div}\, \boldsymbol{v} = q, \f]
- * subject to appropriate boundary conditions.
- */
- void pressure(bool solveTwice = true)
+ void update()
{
- assemble(false);
- solve();
-
- return;
+ ParentType::update();
+ storePressureSolution();
}
- // serialization methods
- //! Function needed for restart option.
- template<class Restarter>
- void serialize(Restarter &res)
+ void storePressureSolution()
{
- return;
+ int size = problem_.gridView().size(0);
+ for (int i = 0; i < size; i++)
+ {
+ CellData& cellData = problem_.variables().cellData(i);
+ storePressureSolution(i, cellData);
+ }
}
- //! Function needed for restart option.
- template<class Restarter>
- void deserialize(Restarter &res)
+ void storePressureSolution(int globalIdx, CellData& cellData)
{
- return;
+ Scalar press = this->pressure()[globalIdx];
+
+ cellData.setPressure(press);
}
//! \brief Writes data files
@@ -167,321 +162,171 @@ public:
template<class MultiWriter>
void addOutputVtkFields(MultiWriter &writer)
{
- typename Variables::ScalarSolutionType *pressure = writer.allocateManagedBuffer (
+ ScalarSolutionType *pressure = writer.allocateManagedBuffer (
problem_.gridView().size(0));
- *pressure = problem_.variables().pressure();
+ *pressure = this->pressure();
writer.attachCellData(*pressure, "pressure");
return;
}
- void setPressureHard(Scalar pressure, int globalIdx)
- {
- setPressHard_ = true;
- pressHard_ = pressure;
- idxPressHard_ = globalIdx;
- }
-
- void unsetPressureHard(int globalIdx)
- {
- setPressHard_ = false;
- pressHard_ = 0.0;
- idxPressHard_ = 0.0;
- }
-
//! Constructs a FVPressure1P object
/**
* \param problem a problem class object
*/
FVPressure1P(Problem& problem) :
- problem_(problem), A_(problem.variables().gridSize(), problem.variables().gridSize(), (2 * dim + 1)
- * problem.variables().gridSize(), Matrix::random), f_(problem.variables().gridSize()),
- pressHard_(0),
- idxPressHard_(0),
- setPressHard_(false),
- gravity(
- problem.gravity())
+ ParentType(problem), problem_(problem),
+ gravity_(
+ problem.gravity())
{
- initializeMatrix();
+ 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);
}
private:
Problem& problem_;
- Matrix A_;
- Dune::BlockVector<Dune::FieldVector<Scalar, 1> > f_;
- Scalar pressHard_;
- Scalar idxPressHard_;
- bool setPressHard_;
-protected:
- const GlobalPosition& gravity; //!< vector including the gravity constant
+ const GlobalPosition& gravity_; //!< vector including the gravity constant
+ Scalar density_;
+ Scalar viscosity_;
};
-//!initializes the matrix to store the system of equations
+//!function which calculates the source entry
template<class TypeTag>
-void FVPressure1P<TypeTag>::initializeMatrix()
+void FVPressure1P<TypeTag>::getSource(Dune::FieldVector<Scalar, 2>& entries, const Element& element
+ , const CellData& cellData, const bool first)
{
- // determine matrix row sizes
- 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);
-
- // initialize row size
- int rowSize = 1;
-
- // run through all intersections with neighbors
- IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
- for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
- if (isIt->neighbor())
- rowSize++;
- A_.setrowsize(globalIdxI, rowSize);
- }
- A_.endrowsizes();
+ // cell volume, assume linear map here
+ Scalar volume = element.geometry().volume();
- // determine position of matrix entries
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
- {
- // cell index
- int globalIdxI = problem_.variables().index(*eIt);
-
- // add diagonal index
- A_.addindex(globalIdxI, globalIdxI);
-
- // run through all intersections with neighbors
- IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
- for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt); isIt != isItEnd; ++isIt)
- if (isIt->neighbor())
- {
- // access neighbor
- ElementPointer outside = isIt->outside();
- int globalIdxJ = problem_.variables().index(*outside);
-
- // add off diagonal index
- A_.addindex(globalIdxI, globalIdxJ);
- }
- }
- A_.endindices();
+ // get sources from problem
+ PrimaryVariables sourcePhase(0.0);
+ problem_.source(sourcePhase, element);
+ sourcePhase /= density_;
+
+ entries[rhs] = volume * sourcePhase;
return;
}
-//!function which assembles the system of equations to be solved
+//!function which calculates internal flux entries
template<class TypeTag>
-void FVPressure1P<TypeTag>::assemble(bool first)
+void FVPressure1P<TypeTag>::getFlux(Dune::FieldVector<Scalar, 2>& entries, const Intersection& intersection
+ , const CellData& cellDataI, const bool first)
{
- // initialization: set matrix A_ to zero
- A_ = 0;
- f_ = 0;
+ ElementPointer elementI = intersection.inside();
+ ElementPointer elementJ = intersection.outside();
- BoundaryTypes bcType;
+ // get global coordinates of cell centers
+ const GlobalPosition& globalPosI = elementI->geometry().center();
+ const GlobalPosition& globalPosJ = elementJ->geometry().center();
- ElementIterator eItEnd = problem_.gridView().template end<0> ();
- for (ElementIterator eIt = problem_.gridView().template begin<0> (); eIt != eItEnd; ++eIt)
- {
- // get global coordinate of cell center
- const GlobalPosition& globalPos = eIt->geometry().center();
+ //get face normal
+ const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
- // cell index
- int globalIdxI = problem_.variables().index(*eIt);
+ // get face area
+ Scalar faceArea = intersection.geometry().volume();
- // cell volume, assume linear map here
- Scalar volume = eIt->geometry().volume();
+ // distance vector between barycenters
+ GlobalPosition distVec = globalPosJ - globalPosI;
- Scalar temperatureI = problem_.temperature(*eIt);
- Scalar referencePressI = problem_.referencePressure(*eIt);
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
- Scalar densityI = Fluid::density(temperatureI, referencePressI);
- Scalar viscosityI = Fluid::viscosity(temperatureI, referencePressI);
+ // compute vectorized permeabilities
+ FieldMatrix meanPermeability(0);
- // set right side to zero
- PrimaryVariables source(0.0);
- problem_.source(source, *eIt);
- source /= densityI;
+ problem_.spatialParameters().meanK(meanPermeability, problem_.spatialParameters().intrinsicPermeability(*elementI),
+ problem_.spatialParameters().intrinsicPermeability(*elementJ));
- f_[globalIdxI] = source *= volume;
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ meanPermeability.mv(unitOuterNormal, permeability);
- int isIndex = 0;
- IntersectionIterator isItEnd = problem_.gridView().iend(*eIt);
- for (IntersectionIterator isIt = problem_.gridView().ibegin(*eIt);
- isIt != isItEnd;
- ++isIt, ++isIndex)
- {
- // get normal vector
- Dune::FieldVector < Scalar, dimWorld > unitOuterNormal = isIt->centerUnitOuterNormal();
+ permeability/=viscosity_;
- // get face volume
- Scalar faceArea = isIt->geometry().volume();
+ //calculate current matrix entry
+ entries[matrix] = ((permeability * unitOuterNormal) / dist) * faceArea;
- // handle interior face
- if (isIt->neighbor())
- {
- // access neighbor
- ElementPointer neighborPointer = isIt->outside();
- int globalIdxJ = problem_.variables().index(*neighborPointer);
+ //calculate right hand side
+ entries[rhs] = density_ * (permeability * gravity_) * faceArea;
- // neighbor cell center in global coordinates
- const GlobalPosition& globalPosNeighbor = neighborPointer->geometry().center();
-
- // distance vector between barycenters
- Dune::FieldVector < Scalar, dimWorld > distVec = globalPosNeighbor - globalPos;
-
- // compute distance between cell centers
- Scalar dist = distVec.two_norm();
-
- // compute vectorized permeabilities
- FieldMatrix meanPermeability(0);
-
- problem_.spatialParameters().meanK(meanPermeability,
- problem_.spatialParameters().intrinsicPermeability(*eIt),
- problem_.spatialParameters().intrinsicPermeability(*neighborPointer));
-
- Dune::FieldVector<Scalar, dim> permeability(0);
- meanPermeability.mv(unitOuterNormal, permeability);
-
- permeability/=viscosityI;
-
- Scalar temperatureJ = problem_.temperature(*neighborPointer);
- Scalar referencePressJ = problem_.referencePressure(*neighborPointer);
-
- Scalar densityJ = Fluid::density(temperatureJ, referencePressJ);
-
- Scalar rhoMean = 0.5 * (densityI + densityJ);
-
- // update diagonal entry
- Scalar entry;
-
- //calculate potential gradients
- Scalar potential = 0;
-
- Scalar density = 0;
-
- //if we are at the very first iteration we can't calculate phase potentials
- if (!first)
- {
- potential = problem_.variables().potential(globalIdxI, isIndex);
-
- density = (potential > 0.) ? densityI : densityJ;
-
- density = (potential == 0.) ? rhoMean : density;
-
- potential = (problem_.variables().pressure()[globalIdxI]
- - problem_.variables().pressure()[globalIdxJ]) / dist;
-
- potential += density * (unitOuterNormal * gravity);
-
- //store potentials for further calculations (velocity, saturation, ...)
- problem_.variables().potential(globalIdxI, isIndex) = potential;
- }
-
- //do the upwinding depending on the potentials
-
- density = (potential > 0.) ? densityI : densityJ;
-
- density = (potential == 0) ? rhoMean : density;
-
- //calculate current matrix entry
- entry = ((permeability * unitOuterNormal) / dist) * faceArea;
-
- //calculate right hand side
- Scalar rightEntry = density * (permeability * gravity) * faceArea;
-
- //set right hand side
- f_[globalIdxI] -= rightEntry;
-
- // set diagonal entry
- A_[globalIdxI][globalIdxI] += entry;
-
- // set off-diagonal entry
- A_[globalIdxI][globalIdxJ] = -entry;
- }
-
- // boundary face
+ return;
+}
- else if (isIt->boundary())
- {
- // center of face in global coordinates
- const GlobalPosition& globalPosFace = isIt->geometry().center();
+//!function which calculates internal flux entries
+template<class TypeTag>
+void FVPressure1P<TypeTag>::getFluxOnBoundary(Dune::FieldVector<Scalar, 2>& entries,
+const Intersection& intersection, const CellData& cellData, const bool first)
+{
+ ElementPointer element = intersection.inside();
- //get boundary condition for boundary face center
- problem_.boundaryTypes(bcType, *isIt);
- PrimaryVariables boundValues(0.0);
+ // get global coordinates of cell centers
+ const GlobalPosition& globalPosI = element->geometry().center();
- if (bcType.isDirichlet(pressEqIdx))
- {
- problem_.dirichlet(boundValues, *isIt);
+ // center of face in global coordinates
+ const GlobalPosition& globalPosJ = intersection.geometry().center();
- Dune::FieldVector < Scalar, dimWorld > distVec(globalPosFace - globalPos);
- Scalar dist = distVec.two_norm();
+ //get face normal
+ const Dune::FieldVector<Scalar, dim>& unitOuterNormal = intersection.centerUnitOuterNormal();
- //permeability vector at boundary
- // compute vectorized permeabilities
- FieldMatrix meanPermeability(0);
+ // get face area
+ Scalar faceArea = intersection.geometry().volume();
- problem_.spatialParameters().meanK(meanPermeability,
- problem_.spatialParameters().intrinsicPermeability(*eIt));
+ // distance vector between barycenters
+ GlobalPosition distVec = globalPosJ - globalPosI;
- //permeability vector at boundary
- Dune::FieldVector < Scalar, dim > permeability(0);
- meanPermeability.mv(unitOuterNormal, permeability);
+ // compute distance between cell centers
+ Scalar dist = distVec.two_norm();
- permeability/= viscosityI;
+ BoundaryTypes bcType;
+ problem_.boundaryTypes(bcType, intersection);
+ PrimaryVariables boundValues(0.0);
- //get dirichlet pressure boundary condition
- Scalar pressBound = boundValues;
+ if (bcType.isDirichlet(pressEqIdx))
+ {
+ problem_.dirichlet(boundValues, intersection);
- //calculate current matrix entry
- Scalar entry = ((permeability * unitOuterNormal) / dist) * faceArea;
+ //permeability vector at boundary
+ // compute vectorized permeabilities
+ FieldMatrix meanPermeability(0);
- //calculate right hand side
- Scalar rightEntry = densityI * (permeability * gravity) * faceArea;
+ problem_.spatialParameters().meanK(meanPermeability,
+ problem_.spatialParameters().intrinsicPermeability(*element));
- // set diagonal entry and right hand side entry
- A_[globalIdxI][globalIdxI] += entry;
- f_[globalIdxI] += entry * pressBound;
- f_[globalIdxI] -= rightEntry;
- }
- //set neumann boundary condition
+ Dune::FieldVector<Scalar, dim> permeability(0);
+ meanPermeability.mv(unitOuterNormal, permeability);
- else if (bcType.isNeumann(pressEqIdx))
- {
- problem_.neumann(boundValues, *isIt);
- Scalar J = boundValues /= densityI;
+ permeability/= viscosity_;
- f_[globalIdxI] -= J * faceArea;
- }
- }
- } // end all intersections
- } // end grid traversal
- return;
-}
+ //get dirichlet pressure boundary condition
+ Scalar pressBound = boundValues;
-//!solves the system of equations to get the spatial distribution of the pressure
-template<class TypeTag>
-void FVPressure1P<TypeTag>::solve()
-{
- typedef typename GET_PROP_TYPE(TypeTag, LinearSolver) Solver;
+ //calculate current matrix entry
+ entries[matrix] = ((permeability * unitOuterNormal) / dist) * faceArea;
+ entries[rhs] = entries[matrix] * pressBound;
- int verboseLevelSolver = GET_PARAM(TypeTag, int, LinearSolver, Verbosity);
+ //calculate right hand side
+ entries[rhs] -= density_ * (permeability * gravity_) * faceArea;
- if (verboseLevelSolver)
- std::cout << "FVPressure1P: solve for pressure" << std::endl;
+ }
+ //set neumann boundary condition
+ else if (bcType.isNeumann(pressEqIdx))
+ {
+ problem_.neumann(boundValues, intersection);
+ Scalar J = boundValues /= density_;
- if (setPressHard_)
+ entries[rhs] = -(J * faceArea);
+ }
+ else
{
- A_[idxPressHard_] = 0;
- A_[idxPressHard_][idxPressHard_] = 1;
- f_[idxPressHard_] = pressHard_;
+ DUNE_THROW(Dune::NotImplemented, "No valid boundary condition type defined for pressure equation!");
}
- Solver solver(problem_);
- solver.solve(A_, problem_.variables().pressure(), f_);
- // printmatrix(std::cout, A_, "global stiffness matrix", "row", 11, 3);
- // printvector(std::cout, f_, "right hand side", "row", 200, 1, 3);
- // printvector(std::cout, (problem_.variables().pressure()), "pressure", "row", 200, 1, 3);
-
return;
}
diff --git a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
index d0103bda52bb0f6710fe3bedb0b337e63ca81b61..3f89aaa9be268a72476bade22984def71355f48e 100644
--- a/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
+++ b/dumux/decoupled/1p/diffusion/fv/fvvelocity1p.hh
@@ -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;
}
}
diff --git a/dumux/decoupled/1p/fluxData1p.hh b/dumux/decoupled/1p/fluxData1p.hh
new file mode 100644
index 0000000000000000000000000000000000000000..14d77a07c1a28bd618409a6319b29aa5b383c9f3
--- /dev/null
+++ b/dumux/decoupled/1p/fluxData1p.hh
@@ -0,0 +1,144 @@
+/*****************************************************************************
+ * 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_FLUXDATA1P_HH
+#define DUMUX_FLUXDATA1P_HH
+
+#include "1pproperties.hh"
+
+/**
+ * @file
+ * @brief Class including the variables and data of discretized data of the constitutive relations
+ * @author Markus Wolff
+ */
+
+namespace Dumux
+{
+/*!
+ * \ingroup IMPES
+ */
+//! Class including the variables and data of discretized data of the constitutive relations.
+/*! 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 FluxData1P
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GridView::Traits::template Codim<0>::Entity Element;
+
+ enum
+ {
+ dim = GridView::dimension, dimWorld = GridView::dimensionworld
+ };
+
+ typedef Dune::FieldVector<Scalar, dim> FieldVector;
+ typedef Dune::FieldVector<FieldVector, 2 * dim> VelocityVector;
+
+ VelocityVector velocity_;
+ Scalar potential_[2 * dim];
+ bool velocityMarker_[2 * dim];
+
+public:
+
+ //! Constructs a FluxData object
+ /**
+ *
+ */
+
+ FluxData1P()
+ {
+ for (int face = 0; face < 2*dim; face++)
+ {
+ velocity_[face] = FieldVector(0.0);
+ potential_[face] = 0.0;
+ velocityMarker_[face] = false;
+ }
+ }
+
+ ////////////////////////////////////////////////////////////
+ // functions returning the vectors of the primary variables
+ ////////////////////////////////////////////////////////////
+
+ //! Return the velocity
+ const FieldVector& velocity(int indexInInside)
+ {
+ return velocity_[indexInInside];
+ }
+
+ const FieldVector& velocity(int indexInInside) const
+ {
+ return velocity_[indexInInside];
+ }
+
+ void setVelocity(int indexInInside, FieldVector& velocity)
+ {
+ velocity_[indexInInside] = velocity;
+ }
+
+ void setVelocityMarker(int indexInInside)
+ {
+ velocityMarker_[indexInInside] = true;
+ }
+
+ bool haveVelocity(int indexInInside)
+ {
+ return velocityMarker_[indexInInside];
+ }
+
+ void resetVelocityMarker()
+ {
+ for (int i = 0; i < 2*dim; i++)
+ velocityMarker_[i] = false;
+ }
+
+ bool isUpwindCell(int indexInInside)
+ {
+ return (potential_[indexInInside] >= 0.);
+ }
+
+ bool isUpwindCell(int indexInInside) const
+ {
+ return (potential_[indexInInside] >= 0.);
+ }
+
+ Scalar potential(int indexInInside)
+ {
+ return potential_[indexInInside];
+ }
+
+ Scalar potential(int indexInInside) const
+ {
+ return potential_[indexInInside];
+ }
+
+ void setPotential(int indexInInside, Scalar pot)
+ {
+ potential_[indexInInside] = pot;
+ }
+
+};
+}
+#endif