diff --git a/dumux/freeflow/stokes/stokesfluxvariables.hh b/dumux/freeflow/stokes/stokesfluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d9ec717ab4fba415fef7168739858d89ab699782
--- /dev/null
+++ b/dumux/freeflow/stokes/stokesfluxvariables.hh
@@ -0,0 +1,280 @@
+/*****************************************************************************
+ * 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
diff --git a/dumux/freeflow/stokes/stokesindices.hh b/dumux/freeflow/stokes/stokesindices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f50133a1fde0ab9978df913ab857833643e4e120
--- /dev/null
+++ b/dumux/freeflow/stokes/stokesindices.hh
@@ -0,0 +1,67 @@
+/*****************************************************************************
+ * 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
diff --git a/dumux/freeflow/stokes/stokeslocaljacobian.hh b/dumux/freeflow/stokes/stokeslocaljacobian.hh
new file mode 100644
index 0000000000000000000000000000000000000000..39bb433c3e937d8593743e900ab34f4bdd5860b0
--- /dev/null
+++ b/dumux/freeflow/stokes/stokeslocaljacobian.hh
@@ -0,0 +1,82 @@
+/*****************************************************************************
+ * 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
diff --git a/dumux/freeflow/stokes/stokeslocalresidual.hh b/dumux/freeflow/stokes/stokeslocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..65c4855cba699e35fbf3aa34ad0ab7b0cd324cae
--- /dev/null
+++ b/dumux/freeflow/stokes/stokeslocalresidual.hh
@@ -0,0 +1,638 @@
+/*****************************************************************************
+ * Copyright (C) 2009-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 Element-wise calculation of the Jacobian matrix for problems
+ * using the stokes box model.
+ */
+
+#ifndef DUMUX_STOKES_LOCAL_RESIDUAL_BASE_HH
+#define DUMUX_STOKES_LOCAL_RESIDUAL_BASE_HH
+
+#include <dumux/boxmodels/common/boxmodel.hh>
+//#include <dumux/boxmodels/common/boxcouplinglocalresidual.hh>
+
+#include "stokesproperties.hh"
+#include "stokesvolumevariables.hh"
+#include "stokesfluxvariables.hh"
+
+#include <dune/common/collectivecommunication.hh>
+#include <vector>
+#include <iostream>
+
+#include<dune/grid/common/grid.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup StokesModel
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the stokes box model.
+ *
+ * This class is also used for the non-isothermal and the stokes transport
+ * model, which means that it uses static polymorphism.
+ */
+template<class TypeTag>
+class StokesLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
+{
+protected:
+ typedef typename GET_PROP_TYPE(TypeTag, BaseLocalResidual) ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, StokesIndices) Indices;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ enum {
+ dim = GridView::dimension,
+ numEq = GET_PROP_VALUE(TypeTag, NumEq)
+ };
+ enum {
+ massBalanceIdx = Indices::massBalanceIdx, //!< Index of the mass balance
+
+ momentumXIdx = Indices::momentumXIdx, //!< Index of the x-component of the momentum balance
+ momentumYIdx = Indices::momentumYIdx, //!< Index of the y-component of the momentum balance
+ momentumZIdx = Indices::momentumZIdx, //!< Index of the z-component of the momentum balance
+ lastMomentumIdx = Indices::lastMomentumIdx //!< Index of the last component of the momentum balance
+ };
+ enum {
+ dimXIdx = Indices::dimXIdx, //!< Index for the first component of a vector
+ dimYIdx = Indices::dimYIdx, //!< Index for the second component of a vector
+ dimZIdx = Indices::dimZIdx //!< Index for the third component of a vector
+ };
+
+ typedef typename GridView::ctype CoordScalar;
+ typedef Dune::GenericReferenceElements<Scalar, dim> ReferenceElements;
+ typedef Dune::GenericReferenceElement<Scalar, dim> ReferenceElement;
+
+ typedef Dune::FieldVector<CoordScalar, dim> GlobalPosition;
+ typedef Dune::FieldVector<Scalar, dim> ScalarGradient;
+ typedef Dune::FieldVector<Scalar, dim> FieldVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ typedef typename GridView::Intersection Intersection;
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+ typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+
+public:
+ /*!
+ * \brief Constructor. Sets the upwind weight.
+ */
+ StokesLocalResidual()
+ {
+ // retrieve the upwind weight for the mass conservation equations. Use the value
+ // specified via the property system as default, and overwrite
+ // it by the run-time parameter from the Dune::ParameterTree
+ massUpwindWeight_ = GET_PARAM(TypeTag, Scalar, MassUpwindWeight);
+ stabilizationAlpha_ = GET_PARAM(TypeTag, Scalar, StabilizationAlpha);
+ stabilizationBeta_ = GET_PARAM(TypeTag, Scalar, StabilizationBeta);
+ };
+
+ /*!
+ * \brief Evaluate the amount all conservation quantities
+ * (e.g. phase mass) within a sub-control volume.
+ *
+ * The result should be averaged over the volume (e.g. phase mass
+ * inside a sub control volume divided by the volume)
+ *
+ * \param result The mass of the component within the sub-control volume
+ * \param scvIdx The SCV (sub-control-volume) index
+ * \param usePrevSol Evaluate function with solution of current or previous time step
+ */
+ void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
+ {
+ // if flag usePrevSol is set, the solution from the previous
+ // time step is used, otherwise the current solution is
+ // used. The secondary variables are used accordingly. This
+ // is required to compute the derivative of the storage term
+ // using the implicit Euler method.
+ const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_()
+ : this->curVolVars_();
+ const VolumeVariables &vertexData = elemVolVars[scvIdx];
+
+ result = 0.0;
+
+ // mass balance
+ result[massBalanceIdx] = vertexData.density();
+
+ // momentum balance
+ for (int momentumIdx = momentumXIdx; momentumIdx <= lastMomentumIdx; ++momentumIdx)
+ result[momentumIdx] = vertexData.density()
+ * vertexData.velocity()[momentumIdx-momentumXIdx];
+ }
+
+ /*!
+ * \brief Evaluates the total flux of all conservation quantities
+ * over a face of a sub-control volume.
+ *
+ * \param flux The flux over the SCV (sub-control-volume) face
+ * \param faceIdx The index of the SCV face
+ */
+ void computeFlux(PrimaryVariables &flux, int faceIdx, bool onBoundary=false) const
+ {
+ const FluxVariables fluxVars(this->problem_(),
+ this->elem_(),
+ this->fvElemGeom_(),
+ faceIdx,
+ this->curVolVars_(),
+ onBoundary);
+ flux = 0.0;
+
+ asImp_()->computeAdvectiveFlux(flux, fluxVars);
+ Valgrind::CheckDefined(flux);
+ asImp_()->computeDiffusiveFlux(flux, fluxVars);
+ Valgrind::CheckDefined(flux);
+ }
+
+ /*!
+ * \brief Evaluates the advective fluxes over
+ * a face of a subcontrol volume.
+ *
+ * \param flux The advective flux over the sub-control-volume face for each component
+ * \param vars The flux variables at the current SCV
+ */
+ void computeAdvectiveFlux(PrimaryVariables &flux,
+ const FluxVariables &fluxVars) const
+ {
+ // if the momentum balance has a dirichlet b.c., the mass balance
+ // is replaced, thus we do not need to calculate outflow fluxes here
+ if (fluxVars.onBoundary())
+ if (momentumBalanceDirichlet_(this->bcTypes_(fluxVars.upstreamIdx())))
+ return;
+
+ // data attached to upstream and the downstream vertices
+ const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx());
+ const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx());
+
+ // mass balance with upwinded density
+ FieldVector massBalanceResidual = fluxVars.velocityAtIP();
+ if (massUpwindWeight_ == 1.0) // fully upwind
+ massBalanceResidual *= up.density();
+ else
+ massBalanceResidual *= massUpwindWeight_ * up.density() +
+ (1.-massUpwindWeight_) * dn.density();
+
+ if (!fluxVars.onBoundary())
+ {
+ // stabilization of the mass balance
+ // with 0.5*alpha*(V_i + V_j)*grad P
+ FieldVector stabilizationTerm = fluxVars.pressureGradAtIP();
+ stabilizationTerm *= stabilizationAlpha_*
+ fluxVars.averageSCVVolume();
+ massBalanceResidual += stabilizationTerm;
+ }
+
+ flux[massBalanceIdx] +=
+ massBalanceResidual*fluxVars.face().normal;
+
+ // momentum balance - pressure is evaluated as volume term
+ // at the center of the SCV in computeSource
+ // viscosity is upwinded
+
+ // compute symmetrized gradient for the momentum flux:
+ // mu (grad v + (grad v)^t)
+ Dune::FieldMatrix<Scalar, dim, dim> symmVelGrad = fluxVars.velocityGradAtIP();
+ for (int i=0; i<dim; ++i)
+ for (int j=0; j<dim; ++j)
+ symmVelGrad[i][j] += fluxVars.velocityGradAtIP()[j][i];
+
+ FieldVector velGradComp(0.);
+ for (int velIdx = 0; velIdx < dim; ++velIdx)
+ {
+ velGradComp = symmVelGrad[velIdx];
+
+ // TODO: dilatation term has to be accounted for in outflow, coupling, neumann
+// velGradComp[velIdx] += 2./3*fluxVars.velocityDivAtIP;
+
+ if (massUpwindWeight_ == 1.0) // fully upwind
+ velGradComp *= up.viscosity();
+ else
+ velGradComp *= massUpwindWeight_ * up.viscosity() +
+ (1.-massUpwindWeight_) * dn.viscosity();
+
+ flux[momentumXIdx + velIdx] -=
+ velGradComp*fluxVars.face().normal;
+
+// gravity is accounted for in computeSource; alternatively:
+// Scalar gravityTerm = fluxVars.densityAtIP *
+// this->problem_().gravity()[dim-1] *
+// fluxVars.face().ipGlobal[dim-1]*
+// fluxVars.face().normal[velIdx];
+// flux[momentumXIdx + velIdx] -=
+// gravityTerm;
+
+ }
+ }
+
+ /*!
+ * \brief Adds the diffusive flux to the flux vector over
+ * the face of a sub-control volume.
+ *
+ * It doesn't do anything in the Stokes model but is used by the
+ * transport and non-isothermal models to calculate diffusive and
+ * conductive fluxes
+ *
+ * \param flux The diffusive flux over the sub-control-volume face for each component
+ * \param vars The flux variables at the current sub control volume face
+ */
+ void computeDiffusiveFlux(PrimaryVariables &flux,
+ const FluxVariables &fluxData) const
+ { }
+
+ /*!
+ * \brief Calculate the source term of the equation,
+ * compute the pressure gradient at the center of a SCV
+ * and evaluate the gravity term
+ *
+ * \param q The source/sink in the sub control volume for each component
+ * \param localVertexIdx The index of the sub control volume
+ */
+ void computeSource(PrimaryVariables &q, int localVertexIdx)
+ {
+ const ElementVolumeVariables &elemVolVars = this->curVolVars_();
+ const VolumeVariables &vertexData = elemVolVars[localVertexIdx];
+
+ // retrieve the source term intrinsic to the problem
+ this->problem_().source(q,
+ this->elem_(),
+ this->fvElemGeom_(),
+ localVertexIdx);
+
+ // ATTENTION: The source term of the mass balance has to be chosen as
+ // div (q_momentum) in the problem file
+ const Scalar alphaH2 = stabilizationAlpha_*
+ this->fvElemGeom_().subContVol[localVertexIdx].volume;
+ q[massBalanceIdx] *= alphaH2; // stabilization of the source term
+
+ // pressure gradient at the center of the SCV,
+ // the pressure is discretized as volume term,
+ // while -mu grad v is calculated in computeFlux
+ ScalarGradient pressureGradAtSCVCenter(0.0);
+ ScalarGradient grad(0.0);
+
+ for (int vertexIdx = 0; vertexIdx < this->fvElemGeom_().numVertices; vertexIdx++)
+ {
+ grad = this->fvElemGeom_().subContVol[localVertexIdx].gradCenter[vertexIdx];
+ Valgrind::CheckDefined(grad);
+ grad *= elemVolVars[vertexIdx].pressure();
+
+ pressureGradAtSCVCenter += grad;
+ }
+
+ // add the component of the pressure gradient to the respective part
+ // of the momentum equation and take the gravity term into account
+ // signs are inverted, since q is subtracted
+ for (int dimIdx=0; dimIdx<dim; ++dimIdx)
+ {
+ q[momentumXIdx + dimIdx] -= pressureGradAtSCVCenter[dimIdx];
+ q[momentumXIdx + dimIdx] += vertexData.density()*this->problem_().gravity()[dimIdx];
+ }
+ }
+
+ // the stokes model needs a modified treatment of the BCs
+ void evalBoundary_()
+ {
+ assert(this->residual_.size() == this->fvElemGeom_().numVertices);
+ const ReferenceElement &refElem = ReferenceElements::general(this->elem_().geometry().type());
+
+ // loop over vertices of the element
+ for (int vertexIdx = 0; vertexIdx < this->fvElemGeom_().numVertices; vertexIdx++)
+ {
+ // consider only SCVs on the boundary
+ if (this->fvElemGeom_().subContVol[vertexIdx].inner)
+ continue;
+
+ // important at corners of the grid
+ FieldVector momentumResidual(0.0);
+ FieldVector averagedNormal(0.0);
+ int numberOfOuterFaces = 0;
+ // evaluate boundary conditions for the intersections of
+ // the current element
+ const BoundaryTypes &bcTypes = this->bcTypes_(vertexIdx);
+ IntersectionIterator isIt = this->gridView_().ibegin(this->elem_());
+ const IntersectionIterator &endIt = this->gridView_().iend(this->elem_());
+ for (; isIt != endIt; ++isIt)
+ {
+ // handle only intersections on the boundary
+ if (!isIt->boundary())
+ continue;
+
+ // assemble the boundary for all vertices of the current face
+ const int faceIdx = isIt->indexInInside();
+ const int numFaceVertices = refElem.size(faceIdx, 1, dim);
+
+ // loop over the single vertices on the current face
+ for (int faceVertIdx = 0; faceVertIdx < numFaceVertices; ++faceVertIdx)
+ {
+ // only evaluate, if we consider the same face vertex as in the outer
+ // loop over the element vertices
+ if (refElem.subEntity(faceIdx, 1, faceVertIdx, dim)
+ != vertexIdx)
+ continue;
+
+ const int boundaryFaceIdx = this->fvElemGeom_().boundaryFaceIndex(faceIdx, faceVertIdx);
+ const FluxVariables boundaryVars(this->problem_(),
+ this->elem_(),
+ this->fvElemGeom_(),
+ boundaryFaceIdx,
+ this->curVolVars_(),
+ true);
+
+ // the computed residual of the momentum equations is stored
+ // into momentumResidual for the replacement of the mass balance
+ // in case of Dirichlet conditions for the momentum balance;
+ // the fluxes at the boundary are added in the second step
+ if (momentumBalanceDirichlet_(bcTypes))
+ {
+ FieldVector muGradVelNormal(0.);
+ const FieldVector &boundaryFaceNormal =
+ boundaryVars.face().normal;
+
+ boundaryVars.velocityGradAtIP().umv(boundaryFaceNormal, muGradVelNormal);
+ muGradVelNormal *= boundaryVars.viscosityAtIP();
+
+ for (int i=0; i < this->residual_.size(); i++)
+ Valgrind::CheckDefined(this->residual_[i]);
+ for (int dimIdx=0; dimIdx < dim; ++dimIdx)
+ momentumResidual[dimIdx] = this->residual_[vertexIdx][momentumXIdx+dimIdx];
+
+ //Sign is right!!!: boundary flux: -mu grad v n
+ //but to compensate outernormal -> residual - (-mu grad v n)
+ momentumResidual += muGradVelNormal;
+ averagedNormal += boundaryFaceNormal;
+ }
+
+ // evaluate fluxes at a single boundary segment
+ asImp_()->evalNeumannSegment_(isIt, vertexIdx, boundaryFaceIdx, boundaryVars);
+ asImp_()->evalOutflowSegment_(isIt, vertexIdx, boundaryFaceIdx, boundaryVars);
+
+ // count the number of outer faces to determine, if we are on
+ // a corner point and if an interpolation should be done
+ numberOfOuterFaces++;
+ } // end loop over face vertices
+ } // end loop over intersections
+
+ // replace defect at the corner points of the grid
+ // by the interpolation of the primary variables
+ if(!bcTypes.isDirichlet(massBalanceIdx))
+ {
+ if (momentumBalanceDirichlet_(bcTypes))
+ replaceMassbalanceResidual_(momentumResidual, averagedNormal, vertexIdx);
+ else // de-stabilize (remove alpha*grad p - alpha div f
+ // from computeFlux on the boundary)
+ removeStabilizationAtBoundary_(vertexIdx);
+ }
+ if (numberOfOuterFaces == 2)
+ interpolateCornerPoints_(bcTypes, vertexIdx);
+ } // end loop over element vertices
+
+ // evaluate the dirichlet conditions of the element
+ if (this->bcTypes_().hasDirichlet())
+ asImp_()->evalDirichlet_();
+ }
+
+protected:
+ /*!
+ * \brief Add Neumann boundary conditions for a single sub-control
+ * volume face to the local residual.
+ */
+ void evalNeumannSegment_(const IntersectionIterator &isIt,
+ const int scvIdx,
+ const int boundaryFaceIdx,
+ const FluxVariables &boundaryVars)
+ {
+ const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
+
+ // deal with neumann boundaries
+ if (bcTypes.hasNeumann())
+ {
+ // call evalNeumannSegment_() of the base class first
+ ParentType::evalNeumannSegment_(isIt, scvIdx, boundaryFaceIdx);
+
+ // temporary vector to store the neumann boundary fluxes
+ PrimaryVariables values(0.0);
+ if (momentumBalanceHasNeumann_(bcTypes))
+ {
+ // Neumann BC of momentum equation needs special treatment
+ // mathematically Neumann BC: p n - mu grad v n = q
+ // boundary terms: -mu grad v n
+ // implement q * A (from evalBoundarySegment) - p n(unity) A
+ FieldVector pressureCorrection(boundaryVars.face().normal);
+ pressureCorrection *= this->curVolVars_(scvIdx).pressure();
+ for (int momentumIdx = momentumXIdx; momentumIdx <= lastMomentumIdx; momentumIdx++)
+ if(bcTypes.isNeumann(momentumIdx))
+ this->residual_[scvIdx][momentumIdx] += pressureCorrection[momentumIdx];
+
+ // beta-stabilization at the boundary
+ // in case of neumann conditions for the momentum equation;
+ // calculate mu grad v t t
+ // center in the face of the reference element
+ FieldVector tangent;
+ const FieldVector& elementUnitNormal = isIt->centerUnitOuterNormal();
+ tangent[0] = elementUnitNormal[1]; //TODO: 3D
+ tangent[1] = -elementUnitNormal[0];
+ FieldVector tangentialVelGrad;
+ boundaryVars.velocityGradAtIP().mv(tangent, tangentialVelGrad);
+ tangentialVelGrad *= boundaryVars.viscosityAtIP();
+
+ this->residual_[scvIdx][massBalanceIdx] -= stabilizationBeta_*0.5*
+ this->curVolVars_(scvIdx).pressure();
+ this->residual_[scvIdx][massBalanceIdx] -= stabilizationBeta_*0.5*
+ (tangentialVelGrad*tangent);
+
+ for (int momentumIdx = momentumXIdx; momentumIdx <= lastMomentumIdx; momentumIdx++)
+ this->residual_[scvIdx][massBalanceIdx] -= stabilizationBeta_*0.5
+ * values[momentumIdx]*elementUnitNormal[momentumIdx-momentumXIdx];
+ Valgrind::CheckDefined(this->residual_);
+ }
+ }
+ }
+
+ void evalOutflowSegment_(const IntersectionIterator &isIt,
+ const int scvIdx,
+ const int boundaryFaceIdx,
+ const FluxVariables &boundaryVars)
+ {
+ // temporary vector to store the neumann boundary fluxes
+ const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
+
+ // deal with outflow boundaries
+ if (bcTypes.hasOutflow())
+ {
+ PrimaryVariables values(0.0);
+
+ asImp_()->computeFlux(values, boundaryFaceIdx, /*onBoundary=*/true);
+ Valgrind::CheckDefined(values);
+
+ for (int eqIdx = 0; eqIdx < numEq; ++eqIdx)
+ {
+ if (!bcTypes.isOutflow(eqIdx) )
+ continue;
+ // do not calculate outflow for the mass balance
+ // if the momentum balance is dirichlet -
+ // it is replaced in that case
+ if (eqIdx==massBalanceIdx && momentumBalanceDirichlet_(bcTypes))
+ continue;
+ // deduce outflow
+ this->residual_[scvIdx][eqIdx] += values[eqIdx];
+ }
+
+ // beta-stabilization at the boundary in case of outflow condition
+ // for the momentum balance
+ if(momentumBalanceOutflow_(bcTypes) && stabilizationBeta_ != 0)
+ {
+ // calculate mu grad v t t for beta-stabilization
+ // center in the face of the reference element
+ FieldVector tangent;
+ const FieldVector& elementUnitNormal = isIt->centerUnitOuterNormal();
+ tangent[0] = elementUnitNormal[1];
+ tangent[1] = -elementUnitNormal[0];
+ FieldVector tangentialVelGrad;
+ boundaryVars.velocityGradAtIP().mv(tangent, tangentialVelGrad);
+
+ this->residual_[scvIdx][massBalanceIdx] -= 0.5*stabilizationBeta_
+ * boundaryVars.viscosityAtIP()
+ * (tangentialVelGrad*tangent);
+ }
+ }
+ }
+
+ void removeStabilizationAtBoundary_(const int vertexIdx)
+ {
+ // loop over the edges of the element
+ for (int faceIdx = 0; faceIdx < this->fvElemGeom_().numEdges; faceIdx++)
+ {
+ const FluxVariables fluxVars(this->problem_(),
+ this->elem_(),
+ this->fvElemGeom_(),
+ faceIdx,
+ this->curVolVars_());
+
+ const int i = this->fvElemGeom_().subContVolFace[faceIdx].i;
+ const int j = this->fvElemGeom_().subContVolFace[faceIdx].j;
+
+ if (i != vertexIdx && j != vertexIdx)
+ continue;
+
+ const Scalar alphaH2 = stabilizationAlpha_*
+ fluxVars.averageSCVVolume();
+ Scalar stabilizationTerm = fluxVars.pressureGradAtIP() *
+ this->fvElemGeom_().subContVolFace[faceIdx].normal;
+
+ stabilizationTerm *= alphaH2;
+
+ if (vertexIdx == i)
+ this->residual_[i][massBalanceIdx] += stabilizationTerm;
+ if (vertexIdx == j)
+ this->residual_[j][massBalanceIdx] -= stabilizationTerm;
+ }
+
+ //destabilize source term
+ PrimaryVariables q(0.0);
+ this->problem_().source(q,
+ this->elem_(),
+ this->fvElemGeom_(),
+ vertexIdx);
+ const Scalar alphaH2 = stabilizationAlpha_*this->fvElemGeom_().subContVol[vertexIdx].volume;
+ this->residual_[vertexIdx][massBalanceIdx] += alphaH2*q[massBalanceIdx]*
+ this->fvElemGeom_().subContVol[vertexIdx].volume;
+ }
+
+ void interpolateCornerPoints_(const BoundaryTypes &bcTypes, const int vertexIdx)
+ {
+ // TODO: 3D
+ if (bcTypes.isCouplingInflow(massBalanceIdx) || bcTypes.isCouplingOutflow(massBalanceIdx))
+ {
+ if (vertexIdx == 0 || vertexIdx == 3)
+ this->residual_[vertexIdx][massBalanceIdx] =
+ this->curPrimaryVars_(0)[massBalanceIdx]-this->curPrimaryVars_(3)[massBalanceIdx];
+ if (vertexIdx == 1 || vertexIdx == 2)
+ this->residual_[vertexIdx][massBalanceIdx] =
+ this->curPrimaryVars_(1)[massBalanceIdx]-this->curPrimaryVars_(2)[massBalanceIdx];
+ }
+ else
+ {
+ if (!bcTypes.isDirichlet(massBalanceIdx)) // do nothing in case of dirichlet
+ this->residual_[vertexIdx][massBalanceIdx] =
+ this->curPrimaryVars_(0)[massBalanceIdx]+this->curPrimaryVars_(3)[massBalanceIdx]-
+ this->curPrimaryVars_(1)[massBalanceIdx]-this->curPrimaryVars_(2)[massBalanceIdx];
+ }
+ }
+
+ void replaceMassbalanceResidual_(const FieldVector& momentumResidual,
+ FieldVector& averagedNormal,
+ const int vertexIdx)
+ {
+ assert(averagedNormal.two_norm() != 0.0);
+
+ // divide averagedNormal by its length
+ averagedNormal /= averagedNormal.two_norm();
+ // replace the mass balance by the sum of the momentum balances
+ this->residual_[vertexIdx][massBalanceIdx] = momentumResidual*averagedNormal;
+ }
+
+ // returns true, if all conditions for the momentum balance are dirichlet
+ bool momentumBalanceDirichlet_(const BoundaryTypes& bcTypes) const
+ {
+ for (int momentumIdx=momentumXIdx; momentumIdx<=lastMomentumIdx; ++momentumIdx)
+ if (!bcTypes.isDirichlet(momentumIdx))
+ return false;
+ return true;
+ }
+
+ // returns true, if one condition of the momentum balance is neumann
+ bool momentumBalanceHasNeumann_(const BoundaryTypes& bcTypes) const
+ {
+ for (int momentumIdx=momentumXIdx; momentumIdx<=lastMomentumIdx; ++momentumIdx)
+ if (bcTypes.isNeumann(momentumIdx))
+ return true;
+ return false;
+ }
+
+ // returns true, if all conditions for the momentum balance are outlow
+ bool momentumBalanceOutflow_(const BoundaryTypes& bcTypes) const
+ {
+ for (int momentumIdx=momentumXIdx; momentumIdx<=lastMomentumIdx; ++momentumIdx)
+ if (!bcTypes.isOutflow(momentumIdx))
+ return false;
+ return true;
+ }
+
+ Scalar massUpwindWeight_;
+ Scalar stabilizationAlpha_;
+ Scalar stabilizationBeta_;
+
+ Implementation *asImp_()
+ { return static_cast<Implementation *>(this); }
+ const Implementation *asImp_() const
+ { return static_cast<const Implementation *>(this); }
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/stokes/stokesmodel.hh b/dumux/freeflow/stokes/stokesmodel.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1d2de2eb2ffd917c6f6d0aabf5874e97879ec400
--- /dev/null
+++ b/dumux/freeflow/stokes/stokesmodel.hh
@@ -0,0 +1,281 @@
+/*****************************************************************************
+ * 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)
+ {
+ int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ volVars.update(sol[globalIdx],
+ this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ i,
+ false);
+
+ // int globalIdx = dofMapper.map(*elemIt, i, dim);
+ vol = fvElemGeom.subContVol[i].volume;
+ poro = volVars.porosity;
+ rhoN = volVars.density;
+ pW = volVars.pressure;
+ // Te = asImp_()->temperature((*sol)[globalIdx]);
+
+
+ massNPhase = vol * poro * rhoN;
+
+ // get minimum and maximum values of primary variables
+ minP = std::min(minP, pW);
+ maxP = std::max(maxP, pW);
+ // minTe = std::min(minTe, Te);
+ // maxTe = std::max(maxTe, Te);
+
+ // calculate total mass
+ mass[0] += massNPhase; // mass nonwetting phase
+ }
+ }
+
+ // IF PARALLEL: calculate total mass including all processors
+ // also works for sequential calculation
+ mass = this->problem_.gridView().comm().sum(mass);
+
+ if(this->problem_.gridView().comm().rank() == 0) // IF PARALLEL: only print by processor with rank() == 0
+ {
+ // print minimum and maximum values
+ std::cout << "nonwetting phase saturation: min = "<< minSat
+ << ", max = "<< maxSat << std::endl;
+ std::cout << "wetting phase pressure: min = "<< minP
+ << ", max = "<< maxP << std::endl;
+ // std::cout << "temperature: min = "<< minTe
+ // << ", max = "<< maxTe << std::endl;
+ }
+ }
+
+
+ /*!
+ * \brief Append all quantities of interest which can be derived
+ * from the solution of the current time step to the VTK
+ * writer.
+ */
+ template <class MultiWriter>
+ void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer)
+ {
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > VelocityField;
+
+ // create the required scalar fields
+ unsigned numVertices = this->gridView_().size(dim);
+ ScalarField &pN = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &delP = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &rho = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &mu = *writer.allocateManagedBuffer(numVertices);
+ VelocityField &velocity = *writer.template allocateManagedBuffer<Scalar, dim> (numVertices);
+
+ unsigned numElements = this->gridView_().size(0);
+ ScalarField &rank = *writer.allocateManagedBuffer(numElements);
+
+ FVElementGeometry fvElemGeom;
+ VolumeVariables volVars;
+ ElementBoundaryTypes elemBcTypes;
+
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ ElementIterator endit = this->gridView_().template end<0>();
+ for (; elemIt != endit; ++elemIt)
+ {
+ int idx = this->elementMapper().map(*elemIt);
+ rank[idx] = this->gridView_().comm().rank();
+
+ fvElemGeom.update(this->gridView_(), *elemIt);
+ elemBcTypes.update(this->problem_(), *elemIt, fvElemGeom);
+
+ int numLocalVerts = elemIt->template count<dim>();
+ for (int i = 0; i < numLocalVerts; ++i)
+ {
+ int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ volVars.update(sol[globalIdx],
+ this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ i,
+ false);
+
+ pN[globalIdx] = volVars.pressure();
+ delP[globalIdx] = volVars.pressure() - 1e5;
+ rho[globalIdx] = volVars.density();
+ mu[globalIdx] = volVars.viscosity();
+ velocity[globalIdx] = volVars.velocity();
+ };
+ }
+ writer.attachVertexData(pN, "P");
+ writer.attachVertexData(delP, "delP");
+ writer.attachVertexData(rho, "rho");
+ writer.attachVertexData(mu, "mu");
+ writer.attachVertexData(velocity, "v", dim);
+ }
+};
+}
+
+#endif
diff --git a/dumux/freeflow/stokes/stokesnewtoncontroller.hh b/dumux/freeflow/stokes/stokesnewtoncontroller.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f3bbec33762d8bc48df9a019447368e86ac25589
--- /dev/null
+++ b/dumux/freeflow/stokes/stokesnewtoncontroller.hh
@@ -0,0 +1,53 @@
+/*****************************************************************************
+ * 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 A stokes specific controller for the newton solver.
+ */
+#ifndef DUMUX_STOKES_NEWTON_CONTROLLER_HH
+#define DUMUX_STOKES_NEWTON_CONTROLLER_HH
+
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+namespace Dumux {
+/*!
+ * \brief A Stokes specific controller for the newton solver.
+ */
+template <class TypeTag>
+class StokesNewtonController : public NewtonController<TypeTag>
+{
+ typedef NewtonController<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+public:
+ StokesNewtonController(const Problem &problem)
+ : ParentType(problem)
+ {
+ Dune::FMatrixPrecision<>::set_singular_limit(1e-35);
+
+ this->setRelTolerance(1e-4);
+ this->setTargetSteps(15);
+ this->setMaxSteps(23);
+ };
+};
+}
+
+#endif
diff --git a/dumux/freeflow/stokes/stokesproblem.hh b/dumux/freeflow/stokes/stokesproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3335a5968b675fe85ce04c28ff4ba0870e3cfd75
--- /dev/null
+++ b/dumux/freeflow/stokes/stokesproblem.hh
@@ -0,0 +1,124 @@
+/*****************************************************************************
+ * 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 Base class for all stokes problems which use the box scheme
+ */
+#ifndef DUMUX_STOKES_PROBLEM_HH
+#define DUMUX_STOKES_PROBLEM_HH
+
+#include "dumux/freeflow/stokes/stokesnewtoncontroller.hh"
+
+#include <dumux/boxmodels/common/boxproblem.hh>
+
+
+namespace Dumux
+{
+/*!
+ * \ingroup StokesProblems
+ * \brief Base class for all problems which use the stokes box model
+ *
+ * \todo Please doc me more!
+ */
+template<class TypeTag>
+class StokesProblem : public BoxProblem<TypeTag>
+{
+ typedef BoxProblem<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+ typedef typename GridView::Grid Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::Intersection Intersection;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+ enum {
+ dim = Grid::dimension,
+ dimWorld = Grid::dimensionworld
+ };
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ StokesProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView),
+ gravity_(0)
+ {
+// if (GET_PARAM(TypeTag, bool, EnableGravity))
+// gravity_[dim-1] = -9.81;
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+
+ /*!
+ * \brief Returns the temperature within the domain.
+ *
+ * This method MUST be overwritten by the actual problem.
+ */
+ Scalar temperature() const
+ { return asImp_()->temperature(); };
+
+ /*!
+ * \brief Returns the acceleration due to gravity.
+ *
+ * If the <tt>EnableGravity</tt> property is true, this means
+ * \f$\boldsymbol{g} = ( 0,\dots,\ -9.81)^T \f$, else \f$\boldsymbol{g} = ( 0,\dots, 0)^T \f$
+ */
+ const GlobalPosition &gravity() const
+ { return gravity_; }
+
+ /*!
+ * \brief Evaluate the intrinsic permeability
+ * at the corner of a given element
+ *
+ * \return permeability in x-direction
+ */
+ Scalar permeability(const Element &element,
+ const FVElementGeometry &fvElemGeom,
+ const Intersection &is,
+ int scvIdx,
+ int boundaryFaceIdx) const
+ { DUNE_THROW(Dune::NotImplemented, "permeability()"); }
+
+ // \}
+
+private:
+ //! Returns the implementation of the problem (i.e. static polymorphism)
+ Implementation *asImp_()
+ { return static_cast<Implementation *>(this); }
+
+ //! \copydoc asImp_()
+ const Implementation *asImp_() const
+ { return static_cast<const Implementation *>(this); }
+
+ GlobalPosition gravity_;
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/stokes/stokesproperties.hh b/dumux/freeflow/stokes/stokesproperties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d8a3e45bcb44fe852770ea01bff92d9ab809f6c5
--- /dev/null
+++ b/dumux/freeflow/stokes/stokesproperties.hh
@@ -0,0 +1,197 @@
+/*****************************************************************************
+ * 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 Defines the properties required for the Stokes BOX model.
+ */
+
+#ifndef DUMUX_STOKESPROPERTIES_HH
+#define DUMUX_STOKESPROPERTIES_HH
+
+#include "stokesindices.hh"
+
+#include <dumux/boxmodels/common/boxproperties.hh>
+#include <dumux/freeflow/stokes/stokeslocaljacobian.hh>
+
+#include <dumux/material/fluidsystems/gasphase.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/components/nullcomponent.hh>
+
+#include <dumux/material/fluidsystems/1pfluidsystem.hh>
+#include <dumux/material/fluidstates/immisciblefluidstate.hh>
+
+
+
+namespace Dumux
+{
+/*!
+ * \addtogroup StokesModel
+ */
+// \{
+
+////////////////////////////////
+// forward declarations
+////////////////////////////////
+template<class TypeTag>
+class StokesModel;
+
+template<class TypeTag>
+class StokesLocalResidual;
+
+template <class TypeTag>
+class StokesVolumeVariables;
+
+template <class TypeTag>
+class StokesFluxVariables;
+
+template <class TypeTag>
+class StokesNewtonController;
+
+// \}
+
+////////////////////////////////
+// properties
+////////////////////////////////
+namespace Properties
+{
+
+/*!
+ * \addtogroup StokesModel
+ */
+// \{
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the stokes problems
+NEW_TYPE_TAG(BoxStokes, INHERITS_FROM(BoxModel));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+//NEW_PROP_TAG(Soil); //!< The type of the soil properties object
+NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
+NEW_PROP_TAG(MassUpwindWeight); //!< The value of the upwind parameter for the mobility
+NEW_PROP_TAG(StokesIndices); //!< Enumerations for the Stokes models
+NEW_PROP_TAG(Fluid);
+NEW_PROP_TAG(FluidSystem);
+NEW_PROP_TAG(FluidState);
+NEW_PROP_TAG(StabilizationAlpha);
+NEW_PROP_TAG(StabilizationBeta);
+
+NEW_PROP_TAG(CalculateVariablesAtCenterOfSCV);
+NEW_PROP_TAG(PhaseIndex);
+
+NEW_PROP_TAG(SpatialParameters);
+
+NEW_PROP_TAG(Scaling); //!Defines Scaling of the model
+//////////////////////////////////////////////////////////////////
+// Properties
+//////////////////////////////////////////////////////////////////
+
+//! The local jacobian operator for the stokes box scheme
+SET_TYPE_PROP(BoxStokes, LocalJacobian, Dumux::StokesLocalJacobian<TypeTag>);
+
+SET_PROP(BoxStokes, NumEq) //!< set the number of equations
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+ static const int dim = Grid::dimension;
+ public:
+ static constexpr int value = 1 + dim;
+};
+
+SET_SCALAR_PROP(BoxStokes, Scaling, 1); //!< set scaling to 1 by default
+
+//! Use the Stokes local residual function for the Stokes model
+SET_TYPE_PROP(BoxStokes,
+ LocalResidual,
+ StokesLocalResidual<TypeTag>);
+
+// Use the Stokes specific newton controller for the Stokes model
+SET_PROP(BoxStokes, NewtonController)
+{public:
+ typedef StokesNewtonController<TypeTag> type;
+};
+
+#if HAVE_SUPERLU
+SET_TYPE_PROP(BoxStokes, LinearSolver, SuperLUBackend<TypeTag>);
+#endif
+
+//! the Model property
+SET_TYPE_PROP(BoxStokes, Model, StokesModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxStokes, VolumeVariables, StokesVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxStokes, FluxVariables, StokesFluxVariables<TypeTag>);
+
+//! the upwind factor for the mobility.
+SET_SCALAR_PROP(BoxStokes, MassUpwindWeight, 1.0);
+
+//! The fluid system to use by default
+SET_PROP(BoxStokes, FluidSystem)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, Fluid) Fluid;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public:
+ typedef FluidSystems::OneP<Scalar, Fluid> type;
+};
+
+SET_PROP(BoxStokes, Fluid)
+{ private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+public:
+ typedef Dumux::LiquidPhase<Scalar, Dumux::NullComponent<Scalar> > type;
+};
+
+SET_TYPE_PROP(BoxStokes, StokesIndices, StokesCommonIndices<TypeTag>);
+
+//! Choose the type of the employed fluid state
+SET_PROP(BoxStokes, FluidState)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+public:
+ typedef Dumux::ImmiscibleFluidState<Scalar, FluidSystem> type;
+
+};
+
+SET_INT_PROP(BoxStokes, PhaseIndex, 0);
+
+// \}
+
+//// enable jacobian matrix recycling by default
+//SET_BOOL_PROP(BoxStokes, EnableJacobianRecycling, false);
+//// enable partial reassembling by default
+//SET_BOOL_PROP(BoxStokes, EnablePartialReassemble, false);
+//// set some Newton properties deviating from the default ones
+//
+//SET_SCALAR_PROP(BoxStokes, NewtonRelTolerance, 1e-4);
+//SET_INT_PROP(BoxStokes, NewtonTargetSteps, 10);
+
+}
+}
+
+#endif
diff --git a/dumux/freeflow/stokes/stokesvolumevariables.hh b/dumux/freeflow/stokes/stokesvolumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7ab71e62591f06091103c1dc7bc333fbe36d73fb
--- /dev/null
+++ b/dumux/freeflow/stokes/stokesvolumevariables.hh
@@ -0,0 +1,230 @@
+/*****************************************************************************
+ * 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 Contains the quantities which are constant within a
+ * finite volume in the stokes model.
+ */
+#ifndef DUMUX_STOKES_VOLUME_VARIABLES_HH
+#define DUMUX_STOKES_VOLUME_VARIABLES_HH
+
+#include "stokesproperties.hh"
+#include "dumux/boxmodels/common/boxvolumevariables.hh"
+
+#include <dumux/material/fluidstates/immisciblefluidstate.hh>
+
+//#include <dumux/material/fluidstates/compositionalfluidstate.hh>
+//#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
+//#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup StokesModel
+ * \ingroup BoxVolumeVariables
+ * \brief Contains the quantities which are are constant within a
+ * finite volume in the Stokes model.
+ */
+template <class TypeTag>
+class StokesVolumeVariables : public BoxVolumeVariables<TypeTag>
+{
+ typedef BoxVolumeVariables<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
+ typedef typename GET_PROP_TYPE(TypeTag, StokesIndices) Indices;
+
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+ enum {
+ momentumXIdx = Indices::momentumXIdx,
+ lastMomentumIdx = Indices::lastMomentumIdx,
+ pressureIdx = Indices::pressureIdx
+ };
+ enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex) };
+
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+
+ typedef Dune::FieldVector<Scalar, dim> VelocityVector;
+
+public:
+
+ /*!
+ * \brief Update all quantities for a given control volume.
+ */
+ void update(const PrimaryVariables &primaryVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int vertIdx,
+ bool isOldSol)
+ {
+ ParentType::update(primaryVars,
+ problem,
+ element,
+ elemGeom,
+ vertIdx,
+ isOldSol);
+
+ completeFluidState(primaryVars, problem, element, elemGeom, vertIdx, fluidState_, isOldSol);
+
+ for (int dimIdx=momentumXIdx; dimIdx<=lastMomentumIdx; ++dimIdx)
+ velocity_[dimIdx] = primaryVars[dimIdx];
+ }
+
+ static void completeFluidState(const PrimaryVariables& primaryVars,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& elemGeom,
+ int scvIdx,
+ FluidState& fluidState,
+ bool isOldSol = false)
+ {
+ Scalar temperature = Implementation::temperature_(primaryVars, problem, element,
+ elemGeom, scvIdx);
+ fluidState.setTemperature(temperature);
+ fluidState.setPressure(phaseIdx, primaryVars[pressureIdx]);
+
+ // create NullParameterCache and do dummy update
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updateAll(fluidState);
+
+ fluidState.setDensity(phaseIdx,
+ FluidSystem::density(fluidState,
+ paramCache,
+ phaseIdx));
+ fluidState.setViscosity(phaseIdx,
+ FluidSystem::viscosity(fluidState,
+ paramCache,
+ phaseIdx));
+
+ // compute and set the enthalpy
+ Scalar h = Implementation::enthalpy_(fluidState, paramCache, phaseIdx);
+ fluidState.setEnthalpy(phaseIdx, h);
+
+// int globalVertIdx = problem.model().dofMapper().map(element, scvIdx, dim);
+ }
+
+ /*!
+ * \brief Returns the phase state for the control-volume.
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+ FluidState &fluidState()
+ { return fluidState_; }
+
+ /*!
+ * \brief Returns the mass density of a given phase within the
+ * control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar density() const
+ { return fluidState_.density(phaseIdx); }
+
+ /*!
+ * \brief Returns the effective pressure of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar pressure() const
+ { return fluidState_.pressure(phaseIdx); }
+
+ /*!
+ * \brief Returns temperature inside the sub-control volume.
+ *
+ * Note that we assume thermodynamic equilibrium, i.e. the
+ * temperature of the rock matrix and of all fluid phases are
+ * identical.
+ */
+ Scalar temperature() const
+ { return fluidState_.temperature(phaseIdx); }
+
+ /*!
+ * \brief Returns the viscosity of the fluid in
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar viscosity() const
+ { return fluidState_.viscosity(phaseIdx); }
+
+ /*!
+ * \brief Returns the velocity
+ */
+ const VelocityVector &velocity() const
+ { return velocity_; }
+
+protected:
+ template<class ParameterCache>
+ static Scalar enthalpy_(const FluidState& fluidState,
+ const ParameterCache& paramCache,
+ int phaseIdx)
+ {
+ return 0;
+ }
+
+ static Scalar temperature_(const PrimaryVariables &primaryVars,
+ const Problem& problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx)
+ {
+ return problem.boxTemperature(element, elemGeom, scvIdx);
+ }
+
+// /*!
+// * \brief Called by update() to compute the energy related quantities
+// */
+// void updateEnergy_(const PrimaryVariables &primaryVars,
+// const Problem &problem,
+// const Element &element,
+// const FVElementGeometry &elemGeom,
+// int vertIdx,
+// bool isOldSol)
+// { }
+
+ VelocityVector velocity_;
+ FluidState fluidState_;
+
+private:
+
+ Implementation &asImp()
+ { return *static_cast<Implementation*>(this); }
+ const Implementation &asImp() const
+ { return *static_cast<const Implementation*>(this); }
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/stokes2c/stokes2cfluxvariables.hh b/dumux/freeflow/stokes2c/stokes2cfluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4dda5ed9719449d7477687a146752f6443907d96
--- /dev/null
+++ b/dumux/freeflow/stokes2c/stokes2cfluxvariables.hh
@@ -0,0 +1,136 @@
+/*****************************************************************************
+ * 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 the
+ * component fluxes over a face of a finite volume.
+ *
+ * This means concentration gradients, diffusion coefficients, mass fractions, etc.
+ * at the integration point.
+ */
+#ifndef DUMUX_STOKES2C_FLUX_VARIABLES_HH
+#define DUMUX_STOKES2C_FLUX_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+NEW_PROP_TAG(Stokes2cIndices); //!< Enumerations for the compositional stokes models
+}
+
+/*!
+ * \ingroup Stokes2cModel
+ * \brief This template class contains data which is required to
+ * calculate the component fluxes over a face of a finite
+ * volume for the compositional stokes model.
+ *
+ * This means concentration gradients, diffusion coefficients, mass fractions, etc.
+ * at the integration point.
+ */
+template <class TypeTag>
+class Stokes2cFluxVariables : public StokesFluxVariables<TypeTag>
+{
+ typedef StokesFluxVariables<TypeTag> ParentType;
+ 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;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, Stokes2cIndices) Indices;
+
+ enum { dim = GridView::dimension };
+ enum {
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx
+ };
+ enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex) };
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef Dune::FieldVector<Scalar, dim> ScalarGradient;
+
+ typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+
+public:
+ Stokes2cFluxVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int faceIdx,
+ const ElementVolumeVariables &elemVolVars,
+ bool onBoundary = false)
+ : ParentType(problem, element, elemGeom, faceIdx, elemVolVars, onBoundary)
+ {
+ calculateValues_(problem, element, elemVolVars);
+ }
+
+ Scalar massFractionAtIP() const
+ { return massFractionAtIP_; }
+
+ Scalar diffusionCoeffAtIP() const
+ { return diffusionCoeffAtIP_; }
+
+ const ScalarGradient &massFractionGradAtIP() const
+ { return massFractionGradAtIP_; }
+
+protected:
+ void calculateValues_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemVolVars)
+ {
+ massFractionAtIP_ = Scalar(0);
+ diffusionCoeffAtIP_ = Scalar(0);
+ massFractionGradAtIP_ = Scalar(0);
+
+ // calculate gradients and secondary variables at IPs
+ for (int idx = 0;
+ idx < this->fvGeom_.numVertices;
+ idx++) // loop over vertices of the element
+ {
+ massFractionAtIP_ += elemVolVars[idx].fluidState().massFraction(phaseIdx, lCompIdx) *
+ this->face().shapeValue[idx];
+ diffusionCoeffAtIP_ += elemVolVars[idx].diffusionCoeff() *
+ this->face().shapeValue[idx];
+
+ // the gradient of the mass fraction at the IP
+ for (int dimIdx=0; dimIdx<dim; ++dimIdx)
+ massFractionGradAtIP_ +=
+ this->face().grad[idx][dimIdx]*
+ elemVolVars[idx].fluidState().massFraction(phaseIdx, lCompIdx);
+ };
+
+ Valgrind::CheckDefined(massFractionAtIP_);
+ Valgrind::CheckDefined(diffusionCoeffAtIP_);
+ Valgrind::CheckDefined(massFractionGradAtIP_);
+ }
+
+ Scalar massFractionAtIP_;
+ Scalar diffusionCoeffAtIP_;
+ ScalarGradient massFractionGradAtIP_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/stokes2c/stokes2cindices.hh b/dumux/freeflow/stokes2c/stokes2cindices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7a955a2dc658ee04fc51eaaf051404b013894b29
--- /dev/null
+++ b/dumux/freeflow/stokes2c/stokes2cindices.hh
@@ -0,0 +1,64 @@
+/*****************************************************************************
+ * 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 compositional Stokes BOX model.
+ */
+#ifndef DUMUX_STOKES2C_INDICES_HH
+#define DUMUX_STOKES2C_INDICES_HH
+
+#include <dumux/freeflow/stokes/stokesindices.hh>
+
+namespace Dumux
+{
+// \{
+
+/*!
+ * \brief The common indices for the compositional stokes model.
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int PVOffset = 0>
+struct Stokes2cCommonIndices : public StokesCommonIndices<TypeTag>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+public:
+ // Phase indices
+ static const int lPhaseIdx = FluidSystem::lPhaseIdx; //!< Index of the liquid phase
+ static const int gPhaseIdx = FluidSystem::gPhaseIdx; //!< Index of the gas phase
+
+ // Component indices
+ static const int lCompIdx = 0; //!< Index of the liquid's primary component
+ static const int gCompIdx = 1; //!< Index of the gas' primary component
+
+ // equation and primary variable indices
+ static const int dim = StokesCommonIndices<TypeTag>::dim;
+ static const int transportIdx = PVOffset + dim+1; //! The index for the transport equation.
+ static const int massOrMoleFracIndex = transportIdx; //! The index for the mass or mole fraction in primary variable vectors.
+
+ static const int phaseIdx = gPhaseIdx; //!< Index of the gas phase (required to use the same fluid system in coupled models)
+};
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/stokes2c/stokes2clocalresidual.hh b/dumux/freeflow/stokes2c/stokes2clocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..07928df75e21f0520363fc1a61cb10a53f3df741
--- /dev/null
+++ b/dumux/freeflow/stokes2c/stokes2clocalresidual.hh
@@ -0,0 +1,225 @@
+/*****************************************************************************
+ * 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 Element-wise calculation of the Jacobian matrix for problems
+ * using the compositional stokes box model.
+ *
+ */
+#ifndef DUMUX_STOKES2C_LOCAL_RESIDUAL_HH
+#define DUMUX_STOKES2C_LOCAL_RESIDUAL_HH
+
+#include <dumux/freeflow/stokes/stokeslocalresidual.hh>
+
+#include <dumux/freeflow/stokes2c/stokes2cvolumevariables.hh>
+#include <dumux/freeflow/stokes2c/stokes2cfluxvariables.hh>
+#include <dumux/freeflow/stokes2c/stokes2cproperties.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Stokes2cModel
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the compositional stokes box model. This is derived
+ * from the stokes box model.
+ */
+template<class TypeTag>
+class Stokes2cLocalResidual : public StokesLocalResidual<TypeTag>
+{
+ typedef StokesLocalResidual<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Stokes2cIndices) Indices;
+
+ enum { dim = GridView::dimension };
+ enum {
+ momentumYIdx = Indices::momentumYIdx, //!< Index of the y-component of the momentum balance
+ massBalanceIdx = Indices::massBalanceIdx, //!< Index of the mass balance
+ transportIdx = Indices::transportIdx //!< Index of the transport equation
+ };
+ enum {
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx
+ };
+ enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex)};
+ enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+public:
+ /*!
+ * \brief Constructor. Sets the upwind weight.
+ */
+ Stokes2cLocalResidual()
+ {
+ // retrieve the upwind weight for the mass conservation equations. Use the value
+ // specified via the property system as default, and overwrite
+ // it by the run-time parameter from the Dune::ParameterTree
+ massUpwindWeight_ = GET_PARAM(TypeTag, Scalar, MassUpwindWeight);
+ };
+
+ /*!
+ * \brief Evaluate the amount the additional quantities to the stokes model
+ * (transport equation).
+ *
+ * The result should be averaged over the volume (e.g. phase mass
+ * inside a sub control volume divided by the volume)
+ */
+ void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
+ {
+ // compute the storage term for the transport equation
+ ParentType::computeStorage(result, scvIdx, usePrevSol);
+
+ // if flag usePrevSol is set, the solution from the previous
+ // time step is used, otherwise the current solution is
+ // used. The secondary variables are used accordingly. This
+ // is required to compute the derivative of the storage term
+ // using the implicit euler method.
+ const ElementVolumeVariables &elemDat = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
+ const VolumeVariables &vertexDat = elemDat[scvIdx];
+
+ // compute the storage of the component
+ result[transportIdx] =
+ vertexDat.density() *
+ vertexDat.fluidState().massFraction(phaseIdx, lCompIdx);
+
+ Valgrind::CheckDefined(vertexDat.density());
+ Valgrind::CheckDefined(vertexDat.fluidState().massFraction(phaseIdx, lCompIdx));
+ }
+
+ /*!
+ * \brief Evaluates the advective component (mass) flux
+ * over a face of a subcontrol volume and writes the result in
+ * the flux vector.
+ *
+ * This method is called by compute flux (base class)
+ */
+ void computeAdvectiveFlux(PrimaryVariables &flux,
+ const FluxVariables &fluxVars) const
+ {
+ // call computation of the advective fluxes of the stokes model
+ // (momentum and mass fluxes)
+ ParentType::computeAdvectiveFlux(flux, fluxVars);
+
+ // vertex data of the upstream and the downstream vertices
+ const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx());
+ const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx());
+
+ Scalar tmp = fluxVars.normalVelocityAtIP();
+
+ if (massUpwindWeight_ > 0.0)
+ tmp *= massUpwindWeight_ * // upwind data
+ up.density() * up.fluidState().massFraction(phaseIdx, lCompIdx);
+ if (massUpwindWeight_ < 1.0)
+ tmp += (1 - massUpwindWeight_) * // rest
+ dn.density() * dn.fluidState().massFraction(phaseIdx, lCompIdx);
+
+ flux[transportIdx] += tmp;
+ Valgrind::CheckDefined(flux[transportIdx]);
+ }
+
+ /*!
+ * \brief Adds the diffusive component flux to the flux vector over
+ * the face of a sub-control volume.
+ */
+ void computeDiffusiveFlux(PrimaryVariables &flux,
+ const FluxVariables &fluxVars) const
+ {
+ // diffusive mass flux
+ ParentType::computeDiffusiveFlux(flux, fluxVars);
+
+ // diffusive component flux
+ for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
+ flux[transportIdx] -=
+ fluxVars.massFractionGradAtIP()[dimIdx] *
+ fluxVars.face().normal[dimIdx] *
+ fluxVars.diffusionCoeffAtIP() *
+ fluxVars.densityAtIP();
+ Valgrind::CheckDefined(flux[transportIdx]);
+ }
+
+ // handle boundary conditions for a single sub-control volume face
+ // evaluate one part of the Dirichlet-like conditions for the massfraction
+ // rest is done in local coupling operator
+ void evalCouplingVertex_(const IntersectionIterator &isIt,
+ const int scvIdx,
+ const int boundaryFaceIdx,
+ const FluxVariables& boundaryVars)
+ {
+ ParentType::evalCouplingVertex_(isIt, scvIdx, boundaryFaceIdx, boundaryVars);
+ const VolumeVariables &volVars = this->curVolVars_()[scvIdx];
+
+ // set mass fraction
+ if (this->bcTypes_(scvIdx).isCouplingOutflow(transportIdx))
+ this->residual_[scvIdx][transportIdx] = volVars.fluidState().massFraction(phaseIdx, lCompIdx);
+ Valgrind::CheckDefined(this->residual_[scvIdx][transportIdx]);
+ }
+
+ /*!
+ * \brief Compute the component fluxes at the outflow boundaries. The
+ * rest is done in the local jacobian of the stokes model
+ */
+ void computeOutflowValues_(PrimaryVariables &values,
+ const FluxVariables &boundaryVars,
+ const int scvIdx,
+ const int boundaryFaceIdx)
+ {
+ ParentType::computeOutflowValues_(values,
+ boundaryVars,
+ scvIdx,
+ boundaryFaceIdx);
+
+ // calculate the outflow component flux
+ // density and massfraction are evaluated at the vertex (upwinding),
+ // the velocity at the IP
+ const Scalar advectiveFlux =
+ boundaryVars.velocityAtIP() *
+ boundaryVars.face().normal *
+ this->curVolVars_(scvIdx).density() *
+ this->curVolVars_(scvIdx).fluidState().massFraction(phaseIdx, lCompIdx);
+ // boundaryVars.massFractionAtIP;
+
+ // the diffusive fluxes are evaluated at the IP
+ const Scalar diffusiveFlux =
+ boundaryVars.massFractionGradAtIP() *
+ boundaryVars.face().normal *
+ boundaryVars.diffusionCoeffAtIP() *
+ boundaryVars.densityAtIP();
+
+ Valgrind::CheckDefined(advectiveFlux);
+ Valgrind::CheckDefined(diffusiveFlux);
+
+ values[transportIdx] = advectiveFlux - diffusiveFlux;
+ }
+
+protected:
+ Scalar massUpwindWeight_;
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/stokes2c/stokes2cmodel.hh b/dumux/freeflow/stokes2c/stokes2cmodel.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5f592701184ad0478320bdd7597c160e2550d729
--- /dev/null
+++ b/dumux/freeflow/stokes2c/stokes2cmodel.hh
@@ -0,0 +1,165 @@
+/*****************************************************************************
+ * 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 Adaption of the BOX scheme to the compositional stokes model (with two components).
+ */
+#ifndef DUMUX_STOKES2C_MODEL_HH
+#define DUMUX_STOKES2C_MODEL_HH
+
+#include <dumux/freeflow/stokes/stokesmodel.hh>
+
+#include "stokes2clocalresidual.hh"
+#include "stokes2cproperties.hh"
+#include "stokes2cproblem.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup BoxProblems
+ * \defgroup Stokes2cBoxProblems Compositional stokes problems
+ */
+
+/*!
+ * \ingroup BoxModels
+ * \defgroup Stokes2cModel Compositional box stokes model
+ */
+
+/*!
+ * \ingroup Stokes2cModel
+ * \brief Adaption of the BOX scheme to the compositional stokes model.
+ *
+ * This model implements a non-isothermal flow of a fluid
+ * \f$\alpha \in \{ w, n \}\f$.
+ * Using the standard Stokes approach a mass balance equation is
+ * solved:
+ * \f{eqnarray*}
+ && \phi \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa S_\alpha )}{\partial t}
+ - \sum_\alpha \text{div} \left\{ \varrho_\alpha X_\alpha^\kappa
+ \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
+ (\text{grad} p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\}\\
+ &-& \sum_\alpha \text{div} \left\{{\bf D_{\alpha, pm}^\kappa} \varrho_{\alpha} \text{grad} X^\kappa_{\alpha} \right\}
+ - \sum_\alpha q_\alpha^\kappa = \quad 0 \qquad \kappa \in \{w, a\} \, ,
+ \alpha \in \{w, n\}
+ * \f}
+ \f}
+ *
+ * This is discretized using a fully-coupled vertex
+ * centered finite volume (box) scheme as spatial and
+ * the implicit Euler method as temporal discretization.
+ *
+ */
+template<class TypeTag>
+class Stokes2cModel : public StokesModel<TypeTag>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Stokes2cIndices) Indices;
+
+ enum { dim = GridView::dimension };
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ enum {
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx
+ };
+ enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex) };
+
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+
+ static const Scalar scale_ = GET_PROP_VALUE(TypeTag, Scaling);
+
+public:
+ /*!
+ * \brief Append all quantities of interest which can be derived
+ * from the solution of the current time step to the VTK
+ * writer.
+ */
+ template <class MultiWriter>
+ void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer)
+ {
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > VelocityField;
+
+ // create the required scalar fields
+ unsigned numVertices = this->gridView_().size(dim);
+ ScalarField &pN = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &delP = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &Xw = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &rho = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &mu = *writer.allocateManagedBuffer(numVertices);
+ VelocityField &velocity = *writer.template allocateManagedBuffer<Scalar, dim> (numVertices);
+
+ unsigned numElements = this->gridView_().size(0);
+ ScalarField &rank = *writer.allocateManagedBuffer(numElements);
+
+ FVElementGeometry fvElemGeom;
+ VolumeVariables volVars;
+ ElementBoundaryTypes elemBcTypes;
+
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ ElementIterator endit = this->gridView_().template end<0>();
+ for (; elemIt != endit; ++elemIt)
+ {
+ int idx = this->elementMapper().map(*elemIt);
+ rank[idx] = this->gridView_().comm().rank();
+
+ fvElemGeom.update(this->gridView_(), *elemIt);
+ elemBcTypes.update(this->problem_(), *elemIt, fvElemGeom);
+
+ int numLocalVerts = elemIt->template count<dim>();
+ for (int i = 0; i < numLocalVerts; ++i)
+ {
+ int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+ volVars.update(sol[globalIdx],
+ this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ i,
+ false);
+
+ pN[globalIdx] = volVars.pressure()*scale_;
+ delP[globalIdx] = volVars.pressure()*scale_ - 1e5;
+ Xw[globalIdx] = volVars.fluidState().massFraction(phaseIdx, lCompIdx);
+ rho[globalIdx] = volVars.density()*scale_*scale_*scale_;
+ mu[globalIdx] = volVars.viscosity()*scale_;
+ velocity[globalIdx] = volVars.velocity();
+ velocity[globalIdx] *= 1/scale_;
+ };
+ }
+ writer.attachVertexData(pN, "P");
+ writer.attachVertexData(delP, "delP");
+ writer.attachVertexData(Xw, "X_TRAIL");
+ writer.attachVertexData(rho, "rho");
+ writer.attachVertexData(mu, "mu");
+ writer.attachVertexData(velocity, "v", dim);
+ }
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/stokes2c/stokes2cnewtoncontroller.hh b/dumux/freeflow/stokes2c/stokes2cnewtoncontroller.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d172c8ca7983d1d75176e640c31b1d9399818abe
--- /dev/null
+++ b/dumux/freeflow/stokes2c/stokes2cnewtoncontroller.hh
@@ -0,0 +1,57 @@
+/*****************************************************************************
+ * 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 A stokes2c specific controller for the newton solver.
+ */
+#ifndef DUMUX_STOKES2C_NEWTON_CONTROLLER_HH
+#define DUMUX_STOKES2C_NEWTON_CONTROLLER_HH
+
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+namespace Dumux {
+/*!
+ * \brief A Stokes2c specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+template <class TypeTag>
+class Stokes2cNewtonController : public NewtonController<TypeTag>
+{
+ typedef NewtonController<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+public:
+ Stokes2cNewtonController(const Problem &problem)
+ : ParentType(problem)
+ {
+ Dune::FMatrixPrecision<>::set_singular_limit(1e-35);
+
+ this->setRelTolerance(1e-4);
+ this->setTargetSteps(15);
+ this->setMaxSteps(23);
+ };
+};
+}
+
+#endif
diff --git a/dumux/freeflow/stokes2c/stokes2cproblem.hh b/dumux/freeflow/stokes2c/stokes2cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a61db6683b0ec39a7533915b3f9e2ef77dfe225d
--- /dev/null
+++ b/dumux/freeflow/stokes2c/stokes2cproblem.hh
@@ -0,0 +1,64 @@
+/*****************************************************************************
+ * Copyright (C) 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. *
+ * *
+ * 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 Base class for all problems which use the stokes transport model.
+ */
+#ifndef DUMUX_STOKESTRANSPORT_PROBLEM_HH
+#define DUMUX_STOKESTRANSPORT_PROBLEM_HH
+
+#include "dumux/freeflow/stokes2c/stokes2cnewtoncontroller.hh"
+
+#include <dumux/freeflow/stokes/stokesproblem.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Stokes2cProblems
+ * \brief Base class for all problems which use the
+ * stokes transport box model.
+ *
+ * \todo Please doc me more!
+ */
+template<class TypeTag>
+class Stokes2cProblem : public StokesProblem<TypeTag>
+{
+ typedef StokesProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+public:
+ Stokes2cProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/stokes2c/stokes2cproperties.hh b/dumux/freeflow/stokes2c/stokes2cproperties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4f01f894ae2ba17b2ef7993a356f937f8dc4dfb9
--- /dev/null
+++ b/dumux/freeflow/stokes2c/stokes2cproperties.hh
@@ -0,0 +1,143 @@
+/*****************************************************************************
+ * Copyright (C) 2008 by Klaus Mosthaf, Andreas Lauser, 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 Defines the properties required for the compositional
+ * stokes BOX model.
+ */
+#ifndef DUMUX_STOKES2CPROPERTIES_HH
+#define DUMUX_STOKES2CPROPERTIES_HH
+
+#include "stokes2cindices.hh"
+
+#include <dumux/freeflow/stokes/stokesproperties.hh>
+#include "stokes2cvolumevariables.hh"
+#include "stokes2cfluxvariables.hh"
+
+#include <dumux/material/fluidstates/compositionalfluidstate.hh>
+
+namespace Dumux
+{
+/*!
+ * \addtogroup Stokes2cModel
+ */
+// \{
+////////////////////////////////
+// forward declarations
+////////////////////////////////
+template<class TypeTag>
+class Stokes2cModel;
+
+template<class TypeTag>
+class Stokes2cLocalResidual;
+
+template <class TypeTag>
+class Stokes2cVolumeVariables;
+
+template <class TypeTag>
+class Stokes2cFluxVariables;
+
+template <class TypeTag>
+class Stokes2cNewtonController;
+
+
+////////////////////////////////
+// properties
+////////////////////////////////
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the compositional stokes problems
+NEW_TYPE_TAG(BoxStokes2c, INHERITS_FROM(BoxStokes));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(Stokes2cIndices); //!< Enumerations for the compositional stokes models
+NEW_PROP_TAG(NumComponents); //!< Number of components
+
+//////////////////////////////////////////////////////////////////
+// Properties
+//////////////////////////////////////////////////////////////////
+
+SET_PROP(BoxStokes2c, NumEq) //!< set the number of equations
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+ static const int dim = Grid::dimension;
+ public:
+ static constexpr int value = 2 + dim;
+};
+
+//! Use the stokes2c local jacobian operator for the compositional stokes model
+SET_TYPE_PROP(BoxStokes2c,
+ LocalResidual,
+ Stokes2cLocalResidual<TypeTag>);
+
+//! Use the stokes2c specific newton controller for the compositional stokes model
+SET_PROP(BoxStokes2c, NewtonController)
+{public:
+ typedef Stokes2cNewtonController<TypeTag> type;
+};
+
+//! the Model property
+SET_TYPE_PROP(BoxStokes2c, Model, Stokes2cModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxStokes2c, VolumeVariables, Stokes2cVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxStokes2c, FluxVariables, Stokes2cFluxVariables<TypeTag>);
+
+//! the Indices for the Stokes2c model
+SET_TYPE_PROP(BoxStokes2c,
+ Stokes2cIndices,
+ Stokes2cCommonIndices<TypeTag>);
+
+//! Set the number of components to 2
+SET_INT_PROP(BoxStokes2c, NumComponents, 2);
+
+//! Choose the type of the employed fluid state
+SET_PROP(BoxStokes2c, FluidState)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+public:
+ typedef Dumux::CompositionalFluidState<Scalar, FluidSystem> type;
+
+};
+
+//! Choose the considered phase
+SET_PROP(BoxStokes2c, PhaseIndex)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Stokes2cIndices) Indices;
+public:
+ static constexpr int value = Indices::gPhaseIdx;
+};
+
+}
+
+}
+#endif
diff --git a/dumux/freeflow/stokes2c/stokes2cvolumevariables.hh b/dumux/freeflow/stokes2c/stokes2cvolumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d3c18e0a4d34c9e9c62e65e77a9ec965f2238d13
--- /dev/null
+++ b/dumux/freeflow/stokes2c/stokes2cvolumevariables.hh
@@ -0,0 +1,142 @@
+/*****************************************************************************
+ * 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 *
+ * *
+ * 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 Contains the quantities which are constant within a
+ * finite volume in the compositional stokes model.
+ */
+#ifndef DUMUX_STOKES2C_VOLUME_VARIABLES_HH
+#define DUMUX_STOKES2C_VOLUME_VARIABLES_HH
+
+#include <dumux/freeflow/stokes/stokesvolumevariables.hh>
+#include "stokes2cproperties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Stokes2cModel
+ * \brief Contains the quantities which are are constant within a
+ * finite volume in the compositional Stokes model.
+ */
+template <class TypeTag>
+class Stokes2cVolumeVariables : public StokesVolumeVariables<TypeTag>
+{
+ typedef StokesVolumeVariables<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, Stokes2cIndices) Indices;
+
+ enum { dim = GridView::dimension };
+ enum {
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx
+ };
+ enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+ enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex) };
+
+ enum { transportIdx = Indices::transportIdx };
+
+public:
+ /*!
+ * \brief Update all additional quantities for a given control volume.
+ */
+ void update(const PrimaryVariables &primaryVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int vertIdx,
+ bool isOldSol)
+ {
+ // set the mole fractions first
+ completeFluidState(primaryVars, problem, element, elemGeom, vertIdx, this->fluidState(), isOldSol);
+
+ // update vertex data for the mass and momentum balance
+ ParentType::update(primaryVars,
+ problem,
+ element,
+ elemGeom,
+ vertIdx,
+ isOldSol);
+
+ // Second instance of a parameter cache.
+ // Could be avoided if diffusion coefficients also
+ // became part of the fluid state.
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updateAll(this->fluidState());
+
+ diffCoeff_ = FluidSystem::binaryDiffusionCoefficient(this->fluidState(),
+ paramCache,
+ phaseIdx,
+ lCompIdx,
+ gCompIdx);
+// diffCoeff_ *= GET_RUNTIME_PARAM(TypeTag, Scalar, factor);
+ Valgrind::CheckDefined(diffCoeff_);
+ };
+
+ static void completeFluidState(const PrimaryVariables& primaryVars,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& elemGeom,
+ int scvIdx,
+ FluidState& fluidState,
+ bool isOldSol = false)
+ {
+ Scalar massFraction[numComponents];
+ massFraction[lCompIdx] = primaryVars[transportIdx];
+ massFraction[gCompIdx] = 1 - massFraction[lCompIdx];
+
+ // calculate average molar mass of the gas phase
+ Scalar M1 = FluidSystem::molarMass(lCompIdx);
+ Scalar M2 = FluidSystem::molarMass(gCompIdx);
+ Scalar X2 = massFraction[gCompIdx];
+ Scalar avgMolarMass = M1*M2/(M2 + X2*(M1 - M2));
+
+ // convert mass to mole fractions and set the fluid state
+ fluidState.setMoleFraction(phaseIdx, lCompIdx, massFraction[lCompIdx]*avgMolarMass/M1);
+ fluidState.setMoleFraction(phaseIdx, gCompIdx, massFraction[gCompIdx]*avgMolarMass/M2);
+ }
+
+ /*!
+ * \brief Returns the binary diffusion coefficients for a phase
+ */
+ Scalar diffusionCoeff() const
+ { return diffCoeff_; }
+
+protected:
+// // this method gets called by the parent class. since this method
+// // is protected, we are friends with our parent..
+// friend class StokesVolumeVariables<TypeTag>;
+
+ Scalar diffCoeff_; //!< Binary diffusion coefficient
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh b/dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4ca1382f7bec9f453ecf47ca4a146098090556c4
--- /dev/null
+++ b/dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh
@@ -0,0 +1,116 @@
+/*****************************************************************************
+ * Copyright (C) 2010 by 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 the
+ * energy fluxes over a face of a finite volume.
+ *
+ * This means concentration temperature gradients and heat conductivity at
+ * integration point.
+ */
+#ifndef DUMUX_STOKES2CNI_FLUX_VARIABLES_HH
+#define DUMUX_STOKES2CNI_FLUX_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Stokes2cniModel
+ * \brief This template class contains data which is required to
+ * calculate the energy fluxes over a face of a finite
+ * volume for the non-isothermal compositional stokes model.
+ *
+ * This means temperature gradients, coefficients etc.
+ * at the integration point.
+ */
+template <class TypeTag>
+class Stokes2cniFluxVariables : public Stokes2cFluxVariables<TypeTag>
+{
+ typedef Stokes2cFluxVariables<TypeTag> ParentType;
+ 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, ElementVolumeVariables) ElementVolumeVariables;
+
+ enum { dim = GridView::dimension };
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef Dune::FieldVector<Scalar, dim> ScalarGradient;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+
+public:
+ Stokes2cniFluxVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int faceIdx,
+ const ElementVolumeVariables &elemVolVars,
+ bool onBoundary = false)
+ : ParentType(problem, element, elemGeom, faceIdx, elemVolVars, onBoundary)
+ {
+ calculateValues_(problem, element, elemVolVars);
+ }
+
+ Scalar heatConductivityAtIP() const
+ { return heatConductivityAtIP_; }
+
+ const ScalarGradient &temperatureGradAtIP() const
+ { return temperatureGradAtIP_; }
+
+protected:
+ void calculateValues_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemVolVars)
+ {
+ heatConductivityAtIP_ = Scalar(0);
+ temperatureGradAtIP_ = Scalar(0);
+
+ // calculate gradients and secondary variables at IPs
+ ScalarGradient tmp(0.0);
+ for (int idx = 0;
+ idx < this->fvGeom_.numVertices;
+ idx++) // loop over vertices of the element
+ {
+ heatConductivityAtIP_ += elemVolVars[idx].heatConductivity() *
+ this->face().shapeValue[idx];
+
+ // the gradient of the temperature at the IP
+ for (int dimIdx=0; dimIdx<dim; ++dimIdx)
+ temperatureGradAtIP_ +=
+ this->face().grad[idx][dimIdx]*
+ elemVolVars[idx].temperature();
+ }
+ Valgrind::CheckDefined(heatConductivityAtIP_);
+ Valgrind::CheckDefined(temperatureGradAtIP_);
+ }
+
+ Scalar heatConductivityAtIP_;
+ ScalarGradient temperatureGradAtIP_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/stokes2cni/stokes2cniindices.hh b/dumux/freeflow/stokes2cni/stokes2cniindices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..055970eba6b36d1f6d4f4e5437164db579dcf6f8
--- /dev/null
+++ b/dumux/freeflow/stokes2cni/stokes2cniindices.hh
@@ -0,0 +1,50 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by 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 non-isothermal compositional Stokes BOX model.
+ */
+#ifndef DUMUX_STOKES2CNI_INDICES_HH
+#define DUMUX_STOKES2CNI_INDICES_HH
+
+#include <dumux/freeflow/stokes2c/stokes2cindices.hh>
+
+namespace Dumux
+{
+// \{
+
+/*!
+ * \brief Enumerations for the non-isothermal compositional stokes model
+ */
+template <class TypeTag, int PVOffset=0>
+struct Stokes2cniCommonIndices : public Stokes2cCommonIndices<TypeTag, PVOffset>
+{
+public:
+ // number of dimensions
+ static const int dim = StokesCommonIndices<TypeTag>::dim;
+ static const int energyIdx = PVOffset + dim+2; //! The index for the energy balance equation.
+ static const int temperatureIdx = energyIdx; //! The index for temperature in primary variable vectors.
+};
+} // end namespace
+
+#endif
diff --git a/dumux/freeflow/stokes2cni/stokes2cnilocalresidual.hh b/dumux/freeflow/stokes2cni/stokes2cnilocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..189ba74db64c3c05f115adeb1957df9493a1ea2b
--- /dev/null
+++ b/dumux/freeflow/stokes2cni/stokes2cnilocalresidual.hh
@@ -0,0 +1,173 @@
+/*****************************************************************************
+ * Copyright (C) 2010 by 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 Element-wise calculation of the Jacobian matrix for problems
+ * using the non-isothermal compositional stokes box model.
+ *
+ */
+#ifndef DUMUX_STOKES2CNI_LOCAL_RESIDUAL_HH
+#define DUMUX_STOKES2CNI_LOCAL_RESIDUAL_HH
+
+#include <dumux/freeflow/stokes2c/stokes2clocalresidual.hh>
+
+#include <dumux/freeflow/stokes2cni/stokes2cnivolumevariables.hh>
+#include <dumux/freeflow/stokes2cni/stokes2cnifluxvariables.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Stokes2cniModel
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the non-isothermal compositional stokes box model. This is derived
+ * from the stokes2c box model.
+ */
+template<class TypeTag>
+class Stokes2cniLocalResidual : public Stokes2cLocalResidual<TypeTag>
+{
+ typedef Stokes2cLocalResidual<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Stokes2cniIndices) Indices;
+
+ enum {
+ dim = GridView::dimension,
+ numEq = GET_PROP_VALUE(TypeTag, NumEq)
+ };
+ enum { energyIdx = Indices::energyIdx }; //!< Index of the transport equation
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+ typedef typename GridView::IntersectionIterator IntersectionIterator;
+
+public:
+ /*!
+ * \brief Constructor. Sets the upwind weight.
+ */
+ Stokes2cniLocalResidual()
+ {
+ // retrieve the upwind weight for the mass conservation equations. Use the value
+ // specified via the property system as default, and overwrite
+ // it by the run-time parameter from the Dune::ParameterTree
+ massUpwindWeight_ = GET_PARAM(TypeTag, Scalar, MassUpwindWeight);
+ };
+
+ /*!
+ * \brief Evaluate the amount the additional quantities to the stokes2c model
+ * (energy equation).
+ *
+ * The result should be averaged over the volume (e.g. phase mass
+ * inside a sub control volume divided by the volume)
+ */
+ void computeStorage(PrimaryVariables &result, int scvIdx, bool usePrevSol) const
+ {
+ // compute the storage term for the transport equation
+ ParentType::computeStorage(result, scvIdx, usePrevSol);
+
+ // if flag usePrevSol is set, the solution from the previous
+ // time step is used, otherwise the current solution is
+ // used. The secondary variables are used accordingly. This
+ // is required to compute the derivative of the storage term
+ // using the implicit euler method.
+ const ElementVolumeVariables &elemDat = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
+ const VolumeVariables &vertexDat = elemDat[scvIdx];
+
+ // compute the storage of energy
+ result[energyIdx] =
+ vertexDat.density() *
+ vertexDat.internalEnergy();
+ }
+
+ /*!
+ * \brief Evaluates the convective energy flux
+ * over a face of a subcontrol volume and writes the result in
+ * the flux vector.
+ *
+ * This method is called by compute flux (base class)
+ */
+ void computeAdvectiveFlux(PrimaryVariables &flux,
+ const FluxVariables &fluxVars) const
+ {
+ // call computation of the advective fluxes of the stokes model
+ // (momentum and mass fluxes)
+ ParentType::computeAdvectiveFlux(flux, fluxVars);
+
+ // vertex data of the upstream and the downstream vertices
+ const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx());
+ const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx());
+
+ Scalar tmp = fluxVars.normalVelocityAtIP();
+
+ tmp *= massUpwindWeight_ * // upwind data
+ up.density() * up.enthalpy() +
+ (1 - massUpwindWeight_) * // rest
+ dn.density() * dn.enthalpy();
+
+ flux[energyIdx] += tmp;
+ Valgrind::CheckDefined(flux[energyIdx]);
+ }
+
+ /*!
+ * \brief Adds the conductive energy flux to the flux vector over
+ * the face of a sub-control volume.
+ */
+ void computeDiffusiveFlux(PrimaryVariables &flux,
+ const FluxVariables &fluxVars) const
+ {
+ // diffusive mass flux
+ ParentType::computeDiffusiveFlux(flux, fluxVars);
+
+ // diffusive heat flux
+ for (int dimIdx = 0; dimIdx < dim; ++dimIdx)
+ flux[energyIdx] -=
+ fluxVars.temperatureGradAtIP()[dimIdx] *
+ fluxVars.face().normal[dimIdx] *
+ fluxVars.heatConductivityAtIP();
+ }
+
+ // handle boundary conditions for a single sub-control volume face
+ // evaluate one part of the Dirichlet-like conditions for the temperature
+ // rest is done in local coupling operator
+ void evalCouplingVertex_(const IntersectionIterator &isIt,
+ const int scvIdx,
+ const int boundaryFaceIdx,
+ const FluxVariables& boundaryVars)
+ {
+ ParentType::evalCouplingVertex_(isIt, scvIdx, boundaryFaceIdx, boundaryVars);
+ const VolumeVariables &volVars = this->curVolVars_()[scvIdx];
+
+ if (this->bcTypes_(scvIdx).isCouplingOutflow(energyIdx))
+ this->residual_[scvIdx][energyIdx] = volVars.temperature();
+ }
+
+protected:
+ Scalar massUpwindWeight_;
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/stokes2cni/stokes2cnimodel.hh b/dumux/freeflow/stokes2cni/stokes2cnimodel.hh
new file mode 100644
index 0000000000000000000000000000000000000000..82a476040bf36af0943ac9ffbf9f6c340487cce6
--- /dev/null
+++ b/dumux/freeflow/stokes2cni/stokes2cnimodel.hh
@@ -0,0 +1,174 @@
+/*****************************************************************************
+ * Copyright (C) 2010 by 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 Adaption of the BOX scheme to the non-isothermal
+ * compositional stokes model (with two components).
+ */
+#ifndef DUMUX_STOKES2CNI_MODEL_HH
+#define DUMUX_STOKES2CNI_MODEL_HH
+
+#include <dumux/freeflow/stokes2c/stokes2cmodel.hh>
+
+#include "stokes2cnilocalresidual.hh"
+#include "stokes2cniproperties.hh"
+#include "stokes2cniproblem.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup BoxProblems
+ * \defgroup Stokes2cniBoxProblems Non-isothermal compositional stokes problems
+ */
+
+/*!
+ * \ingroup BoxModels
+ * \defgroup Stokes2cniModel Non-isothermal compositional box stokes model
+ */
+
+/*!
+ * \ingroup Stokes2cniModel
+ * \brief Adaption of the BOX scheme to the non-isothermal compositional stokes model.
+ *
+ * This model implements a non-isothermal flow of a fluid
+ * \f$\alpha \in \{ w, n \}\f$.
+ * Using the standard Stokes approach a mass balance equation is
+ * solved:
+ * \f{eqnarray*}
+ && \phi \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa S_\alpha )}{\partial t}
+ - \sum_\alpha \text{div} \left\{ \varrho_\alpha X_\alpha^\kappa
+ \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
+ (\text{grad} p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\}\\
+ &-& \sum_\alpha \text{div} \left\{{\bf D_{\alpha, pm}^\kappa} \varrho_{\alpha} \text{grad} X^\kappa_{\alpha} \right\}
+ - \sum_\alpha q_\alpha^\kappa = \quad 0 \qquad \kappa \in \{w, a\} \, ,
+ \alpha \in \{w, n\}
+ * \f}
+ \f}
+ *
+ * This is discretized using a fully-coupled vertex
+ * centered finite volume (box) scheme as spatial and
+ * the implicit Euler method as temporal discretization.
+ *
+ */
+template<class TypeTag>
+class Stokes2cniModel : public Stokes2cModel<TypeTag>
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, Stokes2cniIndices) Indices;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld
+ };
+ enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+ enum {
+ lCompIdx = Indices::lCompIdx,
+ gCompIdx = Indices::gCompIdx
+ };
+ enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex) };
+
+ typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+
+public:
+ /*!
+ * \brief Append all quantities of interest which can be derived
+ * from the solution of the current time step to the VTK
+ * writer.
+ */
+ template <class MultiWriter>
+ void addOutputVtkFields(const SolutionVector &sol, MultiWriter &writer)
+ {
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
+ typedef Dune::BlockVector<Dune::FieldVector<Scalar, dim> > VelocityField;
+
+ // create the required scalar fields
+ unsigned numVertices = this->gridView_().size(dim);
+ ScalarField &pN = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &delP = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &Xw = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &T = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &rho = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &mu = *writer.allocateManagedBuffer(numVertices);
+ ScalarField &h = *writer.allocateManagedBuffer(numVertices);
+// ScalarField &D = *writer.allocateManagedBuffer(numVertices);
+ VelocityField &velocity = *writer.template allocateManagedBuffer<Scalar, dim> (numVertices);
+
+ unsigned numElements = this->gridView_().size(0);
+ ScalarField &rank = *writer.allocateManagedBuffer(numElements);
+
+ FVElementGeometry fvElemGeom;
+ VolumeVariables volVars;
+ ElementBoundaryTypes elemBcTypes;
+
+ ElementIterator elemIt = this->gridView_().template begin<0>();
+ ElementIterator endit = this->gridView_().template end<0>();
+ for (; elemIt != endit; ++elemIt)
+ {
+ int idx = this->elementMapper().map(*elemIt);
+ rank[idx] = this->gridView_().comm().rank();
+
+ fvElemGeom.update(this->gridView_(), *elemIt);
+ elemBcTypes.update(this->problem_(), *elemIt, fvElemGeom);
+
+ int numLocalVerts = elemIt->template count<dim>();
+ for (int vertexIdx = 0; vertexIdx < numLocalVerts; ++vertexIdx)
+ {
+ int globalIdx = this->vertexMapper().map(*elemIt, vertexIdx, dim);
+ volVars.update(sol[globalIdx],
+ this->problem_(),
+ *elemIt,
+ fvElemGeom,
+ vertexIdx,
+ false);
+
+ pN [globalIdx] = volVars.pressure();
+ delP[globalIdx] = volVars.pressure() - 1e5;
+ Xw [globalIdx] = volVars.fluidState().massFraction(phaseIdx, lCompIdx);
+ T [globalIdx] = volVars.temperature();
+ rho [globalIdx] = volVars.density();
+ mu [globalIdx] = volVars.viscosity();
+ h [globalIdx] = volVars.enthalpy();
+// D [globalIdx] = volVars.diffusionCoeff();
+ velocity[globalIdx] = volVars.velocity();
+ };
+ }
+ writer.attachVertexData(pN, "pg");
+ writer.attachVertexData(delP, "delP");
+// writer.attachVertexData(D, "Dwg");
+ writer.attachVertexData(Xw, "X_gH2O");
+ writer.attachVertexData(T, "temperature");
+ writer.attachVertexData(rho, "rhoG");
+ writer.attachVertexData(mu, "mu");
+ writer.attachVertexData(h, "h");
+ writer.attachVertexData(velocity, "v", dim);
+ }
+};
+
+}
+#endif
diff --git a/dumux/freeflow/stokes2cni/stokes2cninewtoncontroller.hh b/dumux/freeflow/stokes2cni/stokes2cninewtoncontroller.hh
new file mode 100644
index 0000000000000000000000000000000000000000..31406474d80e874a1e80d72f5612575bc908d277
--- /dev/null
+++ b/dumux/freeflow/stokes2cni/stokes2cninewtoncontroller.hh
@@ -0,0 +1,62 @@
+/*****************************************************************************
+ * Copyright (C) 2008 by Bernd Flemisch, Andreas Lauser *
+ * Copyright (C) 2010 by Klaus Mosthaf *
+ * 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 A stokes2cni specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+#ifndef DUMUX_STOKES2CNI_NEWTON_CONTROLLER_HH
+#define DUMUX_STOKES2CNI_NEWTON_CONTROLLER_HH
+
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+namespace Dumux {
+/*!
+ * \ingroup Stokes2cniModel
+ * \brief A stokes transport specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+template <class TypeTag>
+class Stokes2cniNewtonController : public NewtonController<TypeTag >
+{
+ typedef NewtonController<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+public:
+ Stokes2cniNewtonController(const Problem &problem)
+ : ParentType(problem)
+ {
+ Dune::FMatrixPrecision<>::set_singular_limit(1e-35);
+
+ this->setRelTolerance(1e-6);
+ this->setTargetSteps(9);
+ this->setMaxSteps(18);
+ };
+};
+}
+
+#endif
diff --git a/dumux/freeflow/stokes2cni/stokes2cniproblem.hh b/dumux/freeflow/stokes2cni/stokes2cniproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5dcf9301c44c0320e10bb165068e040ae9e2b180
--- /dev/null
+++ b/dumux/freeflow/stokes2cni/stokes2cniproblem.hh
@@ -0,0 +1,63 @@
+/*****************************************************************************
+ * Copyright (C) 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. *
+ * *
+ * 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 Base class for all problems which use the stokes2cni model.
+ */
+#ifndef DUMUX_STOKES2CNI_PROBLEM_HH
+#define DUMUX_STOKES2CNI_PROBLEM_HH
+
+#include "dumux/freeflow/stokes2cni/stokes2cninewtoncontroller.hh"
+
+#include <dumux/freeflow/stokes/stokesproblem.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup Stokes2cniProblems
+ * \brief Base class for all problems which use the
+ * stokes2cni box model.
+ *
+ * \todo Please doc me more!
+ */
+template<class TypeTag>
+class Stokes2cniProblem : public StokesProblem<TypeTag>
+{
+ typedef StokesProblem<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+
+public:
+ Stokes2cniProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView)
+ {
+ }
+
+ /*!
+ * \name Problem parameters
+ */
+ // \{
+};
+
+}
+
+#endif
diff --git a/dumux/freeflow/stokes2cni/stokes2cniproperties.hh b/dumux/freeflow/stokes2cni/stokes2cniproperties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..24743d089d41831952c25a408a9ed3a4283250d3
--- /dev/null
+++ b/dumux/freeflow/stokes2cni/stokes2cniproperties.hh
@@ -0,0 +1,122 @@
+/*****************************************************************************
+ * Copyright (C) 2008 by Klaus Mosthaf, Andreas Lauser, 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 Defines the properties required for the non-isothermal compositional
+ * stokes BOX model.
+ */
+#ifndef DUMUX_STOKES2CNIPROPERTIES_HH
+#define DUMUX_STOKES2CNIPROPERTIES_HH
+
+#include "stokes2cniindices.hh"
+
+#include <dumux/freeflow/stokes2c/stokes2cproperties.hh>
+#include "stokes2cnivolumevariables.hh"
+#include "stokes2cnifluxvariables.hh"
+
+namespace Dumux
+{
+/*!
+ * \addtogroup Stokes2cniModel
+ */
+// \{
+////////////////////////////////
+// forward declarations
+////////////////////////////////
+template<class TypeTag>
+class Stokes2cniModel;
+
+template<class TypeTag>
+class Stokes2cniLocalResidual;
+
+template <class TypeTag>
+class Stokes2cniVolumeVariables;
+
+template <class TypeTag>
+class Stokes2cniFluxVariables;
+
+template <class TypeTag>
+class Stokes2cniNewtonController;
+
+
+////////////////////////////////
+// properties
+////////////////////////////////
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the compositional stokes problems
+NEW_TYPE_TAG(BoxStokes2cni, INHERITS_FROM(BoxStokes2c));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(Stokes2cniIndices); //!< Enumerations for the compositional stokes models
+NEW_PROP_TAG(NumComponents); //!< Number of components
+
+//////////////////////////////////////////////////////////////////
+// Properties
+//////////////////////////////////////////////////////////////////
+
+SET_PROP(BoxStokes2cni, NumEq) //!< set the number of equations
+{
+ typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+ static const int dim = Grid::dimension;
+ public:
+ static constexpr int value = 3 + dim;
+};
+
+//! Use the stokes2cni local jacobian operator for the compositional stokes model
+SET_TYPE_PROP(BoxStokes2cni,
+ LocalResidual,
+ Stokes2cniLocalResidual<TypeTag>);
+
+//! Use the stokes2cni specific newton controller for the compositional stokes model
+SET_PROP(BoxStokes2cni, NewtonController)
+{public:
+ typedef Stokes2cniNewtonController<TypeTag> type;
+};
+
+//! the Model property
+SET_TYPE_PROP(BoxStokes2cni, Model, Stokes2cniModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxStokes2cni, VolumeVariables, Stokes2cniVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxStokes2cni, FluxVariables, Stokes2cniFluxVariables<TypeTag>);
+
+// the indices for the Stokes2cni model
+SET_TYPE_PROP(BoxStokes2cni,
+ Stokes2cniIndices,
+ Stokes2cniCommonIndices<TypeTag>);
+
+//! Set the number of components to 2
+SET_INT_PROP(BoxStokes2cni, NumComponents, 2);
+}
+
+}
+#endif
diff --git a/dumux/freeflow/stokes2cni/stokes2cnivolumevariables.hh b/dumux/freeflow/stokes2cni/stokes2cnivolumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ab91d00a82f6325dc89ab8d86f0ebbabbd734417
--- /dev/null
+++ b/dumux/freeflow/stokes2cni/stokes2cnivolumevariables.hh
@@ -0,0 +1,139 @@
+/*****************************************************************************
+ * Copyright (C) 2010 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch, Andreas Lauser *
+ * Institute of Hydraulic Engineering *
+ * University of Stuttgart, Germany *
+ * *
+ * 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 Contains the quantities which are constant within a
+ * finite volume in the non-isothermal compositional stokes model.
+ */
+#ifndef DUMUX_STOKES2CNI_VOLUME_VARIABLES_HH
+#define DUMUX_STOKES2CNI_VOLUME_VARIABLES_HH
+
+#include <dumux/freeflow/stokes2c/stokes2cvolumevariables.hh>
+#include "stokes2cniproperties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup Stokes2cniModel
+ * \brief Contains the quantities which are are constant within a
+ * finite volume in the non-isothermal compositional stokes
+ * model.
+ */
+template <class TypeTag>
+class Stokes2cniVolumeVariables : public Stokes2cVolumeVariables<TypeTag>
+{
+ typedef Stokes2cVolumeVariables<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Stokes2cniIndices) Indices;
+ enum { dim = GridView::dimension };
+ enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
+ enum { phaseIdx = GET_PROP_VALUE(TypeTag, PhaseIndex) };
+ enum { energyIdx = Indices::energyIdx };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
+
+public:
+ /*!
+ * \brief Update all additional quantities for a given control volume.
+ */
+ void update(const PrimaryVariables &primaryVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int vertIdx,
+ bool isOldSol)
+ {
+ // the internal energies and the enthalpies
+ heatConductivity_ = 0.0257; //TODO: value (Source: www.engineeringtoolbox.com/air-properties-d_156.html)
+
+ // vertex update data for the mass balance
+ ParentType::update(primaryVars,
+ problem,
+ element,
+ elemGeom,
+ vertIdx,
+ isOldSol);
+ };
+
+ /*!
+ * \brief Returns the total internal energy of a phase in the
+ * sub-control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar internalEnergy() const
+ { return this->fluidState_.internalEnergy(phaseIdx); };
+
+ /*!
+ * \brief Returns the total enthalpy of a phase in the sub-control
+ * volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar enthalpy() const
+ { return this->fluidState_.enthalpy(phaseIdx); };
+
+ /*!
+ * \brief Returns the total heat capacity \f$\mathrm{[J/(K*m^3]}\f$ of the rock matrix in
+ * the sub-control volume.
+ */
+ Scalar heatConductivity() const
+ { return heatConductivity_; };
+
+
+protected:
+ // this method gets called by the parent class. since this method
+ // is protected, we are friends with our parent..
+ friend class StokesVolumeVariables<TypeTag>;
+
+ static Scalar temperature_(const PrimaryVariables &priVars,
+ const Problem& problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx)
+ {
+ return priVars[energyIdx];
+ }
+
+ template<class ParameterCache>
+ static Scalar enthalpy_(const FluidState& fluidState,
+ const ParameterCache& paramCache,
+ int phaseIdx)
+ {
+ return FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
+ }
+
+ Scalar heatConductivity_;
+};
+
+} // end namepace
+
+#endif