diff --git a/dumux/boxmodels/1p2c/1p2cboundaryvariables.hh b/dumux/boxmodels/1p2c/1p2cboundaryvariables.hh
deleted file mode 100644
index eb7c4c75ed66dfb16e3a8ffa3e42e7f3a48b244e..0000000000000000000000000000000000000000
--- a/dumux/boxmodels/1p2c/1p2cboundaryvariables.hh
+++ /dev/null
@@ -1,324 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- * 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 the boundary of a finite volume.
- *
- * This means pressure and temperature gradients, phase densities at
- * the integration point of the boundary, etc.
- */
-#ifndef DUMUX_1P2C_BOUNDARY_VARIABLES_HH
-#define DUMUX_1P2C_BOUNDARY_VARIABLES_HH
-
-#include <dumux/common/math.hh>
-#include "1p2cproperties.hh"
-
-namespace Dumux
-{
-
-/*!
- * \ingroup OnePTwoCModel
- * \brief This template class contains the data which is required to
- * calculate the fluxes of the fluid phases over the boundary of a
- * finite volume for the 1p2c model.
- *
- * This means pressure and velocity gradients, phase density and viscosity at
- * the integration point of the boundary, etc.
- */
-template <class TypeTag>
-class OnePTwoCBoundaryVariables
-{
- 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, OnePTwoCIndices) Indices;
- enum {
- phaseIdx = Indices::phaseIdx,
- comp0Idx = Indices::comp0Idx,
- comp1Idx = Indices::comp1Idx
- };
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
- enum { dim = GridView::dimension };
-
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dune::FieldVector<Scalar, dim> Vector;
- typedef Dune::FieldMatrix<Scalar, dim, dim> Tensor;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename FVElementGeometry::BoundaryFace BoundaryFace;
-
-public:
- OnePTwoCBoundaryVariables(const Problem &problem,
- const Element &element,
- const FVElementGeometry &elemGeom,
- int boundaryFaceIdx,
- const ElementVolumeVariables &elemDat,
- int scvIdx)
- : fvElemGeom_(elemGeom), scvIdx_(scvIdx)
- {
- boundaryFace_ = &fvElemGeom_.boundaryFace[boundaryFaceIdx];
-
- pressureAtBIP_ = Scalar(0);
- massFractionAtBIP_ = Scalar(0);
- densityAtIP_ = Scalar(0);
- viscosityAtIP_ = Scalar(0);
- molarDensityAtIP_ = Scalar(0);
- potentialGrad_ = Scalar(0);
- concentrationGrad_ = Scalar(0);
- moleFracGrad_ = Scalar(0);
- massFracGrad_ = Scalar(0);
- K_= Scalar(0);
-
- calculateBoundaryValues_(problem, element, elemDat);
- };
-
- /*!
- * \brief Return the pressure potential multiplied with the
- * intrinsic permeability which goes from vertex i to
- * vertex j.
- *
- * Note that the length of the face's normal is the area of the
- * phase, so this is not the actual velocity by the integral of
- * the velocity over the face's area. Also note that the phase
- * mobility is not yet included here since this would require a
- * decision on the upwinding approach (which is done in the
- * actual model).
- */
- Scalar KmvpNormal() const
- { return KmvpNormal_; }
-
- /*!
- * \brief The binary diffusion coefficient for each fluid phase in the porous medium.
- */
- Scalar porousDiffCoeff() const
- { return porousDiffCoeff_; };
-
- /*!
- * \brief Return pressure \f$\mathrm{[Pa]}\f$ of a phase at the integration
- * point.
- */
- Scalar pressureAtBIP() const
- { return pressureAtBIP_; }
-
- /*!
- * \brief Return massFraction \f$\mathrm{[-]}\f$ of component 1 at the integration
- * point.
- */
- Scalar massFractionAtBIP() const
- { return massFractionAtBIP_; }
-
- /*!
- * \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration
- * point.
- */
- Scalar densityAtIP() const
- { return densityAtIP_; }
-
- /*!
- * \brief Return viscosity \f$\mathrm{[Pa s]}\f$ of a phase at the integration
- * point.
- */
- Scalar viscosityAtIP() const
- { return viscosityAtIP_; }
-
- /*!
- * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ of a phase at the integration
- * point.
- */
- Scalar molarDensityAtIP() const
- { return molarDensityAtIP_; }
-
- /*!
- * \brief The concentration gradient of a component in a phase.
- *
- * \param compIdx The index of the considered component
- */
- const Vector &concentrationGrad(int compIdx) const
- {
- if (compIdx != 1)
- { DUNE_THROW(Dune::InvalidStateException,
- "The 1p2c model is supposed to need "
- "only the concentration gradient of "
- "the second component!"); }
- return concentrationGrad_;
- };
-
- /*!
- * \brief The molar concentration gradient of a component in a phase.
- *
- * \param compIdx The index of the considered component
- */
- const Vector &moleFracGrad(int compIdx) const
- {
- if (compIdx != 1)
- { DUNE_THROW(Dune::InvalidStateException,
- "The 1p2c model is supposed to need "
- "only the concentration gradient of "
- "the second component!"); }
- return moleFracGrad_;
- };
-
- /*!
- * \brief The mass fraction gradient of a component in a phase.
- *
- * \param compIdx The index of the considered component
- */
- const Vector &massFracGrad(int compIdx) const
- {
- if (compIdx != 1)
- { DUNE_THROW(Dune::InvalidStateException,
- "The 1p2c model is supposed to need "
- "only the concentration gradient of "
- "the second component!"); }
- return massFracGrad_;
- };
-
- const FVElementGeometry &fvElemGeom() const
- { return fvElemGeom_; }
-
- const BoundaryFace& boundaryFace() const
- { return *boundaryFace_; }
-
-protected:
-
- /*!
- * \brief Transforms scalar permeability into tensor.
- *
- * \param k scalar permeability
- *
- */
- void calculateK_(Scalar k)
- {
- K_[0][0] = K_[1][1] = k;
- }
-
- void calculateK_(Tensor K)
- {
- K_ = K;
- }
- /*!
- * \brief Calculation of the pressure and concentration gradients
- *
- * \param problem The considered problem file
- * \param element The considered element of the grid
- * \param elemDat The parameters stored in the considered element
- */
- void calculateBoundaryValues_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemDat)
- {
- Vector tmp(0.0);
-
- // calculate gradients and secondary variables at IPs of the boundary
- for (int idx = 0;
- idx < fvElemGeom_.numVertices;
- idx++) // loop over adjacent vertices
- {
- // FE gradient at vertex idx
- const Vector& feGrad = boundaryFace_->grad[idx];
-
- pressureAtBIP_ += elemDat[idx].pressure() *
- boundaryFace_->shapeValue[idx];
-
- // compute sum of pressure gradients for each phase
- // the pressure gradient
- tmp = feGrad;
- tmp *= elemDat[idx].pressure();
- potentialGrad_ += tmp;
-
- massFractionAtBIP_ +=
- elemDat[idx].massFraction(comp1Idx)
- * boundaryFace_->shapeValue[idx];
- // the concentration gradient of the non-wetting
- // component in the wetting phase
- tmp = feGrad;
- tmp *= elemDat[idx].fluidState().molarity(phaseIdx, comp1Idx);
- concentrationGrad_ += tmp;
-
- tmp = feGrad;
- tmp *= elemDat[idx].fluidState().moleFraction(phaseIdx, comp1Idx);
- moleFracGrad_ += tmp;
-
- tmp = feGrad;
- tmp *= elemDat[idx].fluidState().massFraction(phaseIdx, comp1Idx);
- massFracGrad_ += tmp;
-
-
- densityAtIP_ += elemDat[idx].density()*boundaryFace_->shapeValue[idx];
- viscosityAtIP_ += elemDat[idx].viscosity()*boundaryFace_->shapeValue[idx];
- molarDensityAtIP_ += elemDat[idx].molarDensity()*boundaryFace_->shapeValue[idx];
- }
-
- tmp = problem.gravity();
- tmp *= densityAtIP_;
-
- potentialGrad_ -= tmp;
-
- calculateK_(problem.spatialParameters().intrinsicPermeability(element, fvElemGeom_, scvIdx_));
- Vector Kmvp;
- K_.mv(potentialGrad_, Kmvp);
- KmvpNormal_ = -(Kmvp*boundaryFace_->normal);
-
- const VolumeVariables &vertDat = elemDat[scvIdx_];
-
- Scalar tau = problem.spatialParameters().tortuosity(element, fvElemGeom_, scvIdx_);
-
- porousDiffCoeff_ = vertDat.porosity()*tau*vertDat.diffCoeff();
- }
-
- const FVElementGeometry &fvElemGeom_;
- const BoundaryFace *boundaryFace_;
-
- // gradients
- Vector potentialGrad_;
- Vector concentrationGrad_;
- Vector moleFracGrad_;
- Vector massFracGrad_;
-
- // quanitities at the integration point
- Scalar pressureAtBIP_;
- Scalar massFractionAtBIP_;
- Scalar densityAtIP_;
- Scalar viscosityAtIP_;
- Scalar molarDensityAtIP_;
-
- //intrinsic permeability tensor
- Tensor K_;
- // intrinsic permeability times pressure potential gradient
- // projected on the face normal
- Scalar KmvpNormal_;
-
- // the diffusion coefficient for the porous medium
- Scalar porousDiffCoeff_;
-
- int scvIdx_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/boxmodels/1p2c/1p2cfluxvariables.hh b/dumux/boxmodels/1p2c/1p2cfluxvariables.hh
index 6da4fbb22370f5f98ad6ff1aa73b1c298493079e..182b4d8a0f34fe02246b4b2ba8174306f9eeacad 100644
--- a/dumux/boxmodels/1p2c/1p2cfluxvariables.hh
+++ b/dumux/boxmodels/1p2c/1p2cfluxvariables.hh
@@ -96,12 +96,11 @@ public:
OnePTwoCFluxVariables(const Problem &problem,
const Element &element,
const FVElementGeometry &elemGeom,
- int scvfIdx,
- const ElementVolumeVariables &elemDat)
- : fvElemGeom_(elemGeom)
+ int faceIdx,
+ const ElementVolumeVariables &elemDat,
+ bool onBoundary = false)
+ : fvGeom_(elemGeom), faceIdx_(faceIdx), onBoundary_(onBoundary)
{
- scvfIdx_ = scvfIdx;
-
viscosityAtIP_ = Scalar(0);
molarDensityAtIP_ = Scalar(0);
densityAtIP_ = Scalar(0);
@@ -141,10 +140,16 @@ public:
{ return Kmvp_; }
/*!
- * \brief Return the subcontrol volume face.
+ * \brief The face of the current sub-control volume. This may be either
+ * an inner sub-control-volume face or a SCV face on the boundary.
*/
- const SCVFace &face() const
- { return fvElemGeom_.subContVolFace[scvfIdx_]; }
+ const SCVFace &face() const
+ {
+ if (onBoundary_)
+ return fvGeom_.boundaryFace[faceIdx_];
+ else
+ return fvGeom_.subContVolFace[faceIdx_];
+ }
/*!
* \brief Return the intrinsic permeability tensor.
@@ -298,7 +303,7 @@ protected:
// use finite-element gradients
tmp = 0.0;
for (int idx = 0;
- idx < fvElemGeom_.numVertices;
+ idx < fvGeom_.numVertices;
idx++) // loop over adjacent vertices
{
// FE gradient at vertex idx
@@ -362,8 +367,8 @@ protected:
// estimate the gravitational acceleration at a given SCV face
// using the arithmetic mean
- Vector f(problem.boxGravity(element, fvElemGeom_, face().i));
- f += problem.boxGravity(element, fvElemGeom_, face().j);
+ Vector f(problem.boxGravity(element, fvGeom_, face().i));
+ f += problem.boxGravity(element, fvGeom_, face().j);
f /= 2;
// make it a force
@@ -390,10 +395,10 @@ protected:
const SpatialParameters &sp = problem.spatialParameters();
sp.meanK(K_,
sp.intrinsicPermeability(element,
- fvElemGeom_,
+ fvGeom_,
face().i),
sp.intrinsicPermeability(element,
- fvElemGeom_,
+ fvGeom_,
face().j));
}
@@ -493,8 +498,9 @@ protected:
dispersionTensor_[i][i] += vNorm*dispersivity[1];
}
- const FVElementGeometry &fvElemGeom_;
- int scvfIdx_;
+ const FVElementGeometry &fvGeom_;
+ const int faceIdx_;
+ const bool onBoundary_;
//! pressure potential gradient
Vector potentialGrad_;
diff --git a/dumux/boxmodels/1p2c/1p2clocalresidual.hh b/dumux/boxmodels/1p2c/1p2clocalresidual.hh
index 85599a2bcdebc9236989f757e931df3498a20546..86650517636702860f1599709072a0f894c14d62 100644
--- a/dumux/boxmodels/1p2c/1p2clocalresidual.hh
+++ b/dumux/boxmodels/1p2c/1p2clocalresidual.hh
@@ -38,7 +38,6 @@
#include <dumux/boxmodels/1p2c/1p2cproperties.hh>
#include <dumux/boxmodels/1p2c/1p2cvolumevariables.hh>
#include <dumux/boxmodels/1p2c/1p2cfluxvariables.hh>
-#include <dumux/boxmodels/1p2c/1p2cboundaryvariables.hh>
#include <dune/common/collectivecommunication.hh>
#include <vector>
@@ -67,7 +66,6 @@ protected:
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryVariables) BoundaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
@@ -157,14 +155,15 @@ public:
* \param flux The flux over the SCV (sub-control-volume) face for each component
* \param faceId The index of the considered face of the sub control volume
*/
- void computeFlux(PrimaryVariables &flux, int faceId) const
+ void computeFlux(PrimaryVariables &flux, int faceId, bool onBoundary=false) const
{
flux = 0;
FluxVariables fluxVars(this->problem_(),
this->elem_(),
this->fvElemGeom_(),
faceId,
- this->curVolVars_());
+ this->curVolVars_(),
+ onBoundary);
asImp_()->computeAdvectiveFlux(flux, fluxVars);
asImp_()->computeDiffusiveFlux(flux, fluxVars);
@@ -286,21 +285,13 @@ public:
int scvIdx,
int boundaryFaceIdx)
{
- // temporary vector to store the neumann boundary fluxes
const BoundaryTypes &bcTypes = this->bcTypes_(scvIdx);
- // deal with neumann boundaries
+ // deal with outflow boundaries
if (bcTypes.hasOutflow())
{
- const BoundaryVariables boundaryVars(this->problem_(),
- this->elem_(),
- this->fvElemGeom_(),
- boundaryFaceIdx,
- this->curVolVars_(),
- scvIdx);
-
//calculate outflow fluxes
PrimaryVariables values(0.0);
- asImp_()->computeOutflowValues_(values, boundaryVars, scvIdx, boundaryFaceIdx);
+ asImp_()->computeFlux(values, boundaryFaceIdx, true);
Valgrind::CheckDefined(values);
for (int equationIdx = 0; equationIdx < numEq; ++equationIdx)
@@ -313,56 +304,6 @@ public:
}
}
- /*!
- * \brief Compute the fluxes at the outflow boundaries
- */
- void computeOutflowValues_(PrimaryVariables &values,
- const BoundaryVariables &boundaryVars,
- const int scvIdx,
- const int boundaryFaceIdx)
-
- {
- const VolumeVariables& vertVars = this->curVolVars_()[scvIdx];
-
- // mass balance
- if(!useMoles) //use massfractions
- {
- values[contiEqIdx] += boundaryVars.KmvpNormal()*vertVars.density()/vertVars.viscosity();
- }
- else //use molefractions
- {
- values[contiEqIdx] += boundaryVars.KmvpNormal()*vertVars.molarDensity()/vertVars.viscosity();
- }
-
- // component transport
- if(!useMoles)//use massfractions
- {
- // advective flux
- values[transEqIdx]+=
- boundaryVars.KmvpNormal()*vertVars.density()/vertVars.viscosity()
- *vertVars.fluidState().massFraction(phaseIdx, comp1Idx);
-
- // diffusive flux of comp1 component in phase0
- Scalar tmp = -(boundaryVars.massFracGrad(comp1Idx)*boundaryVars.boundaryFace().normal);
-
- tmp *= boundaryVars.porousDiffCoeff()*boundaryVars.densityAtIP();
- values[transEqIdx] += tmp;//* FluidSystem::molarMass(comp1Idx);
- }
- else //use molefractions
- {
- // advective flux
- values[transEqIdx]+=
- boundaryVars.KmvpNormal()*vertVars.molarDensity()/vertVars.viscosity()
- *vertVars.fluidState().moleFraction(phaseIdx, comp1Idx);
-
- // diffusive flux of comp1 component in phase0
- Scalar tmp = -(boundaryVars.moleFracGrad(comp1Idx)*boundaryVars.boundaryFace().normal);
-
- tmp *= boundaryVars.porousDiffCoeff()*boundaryVars.molarDensityAtIP();
- values[transEqIdx] += tmp;
- }
- }
-
Implementation *asImp_()
{ return static_cast<Implementation *> (this); }
const Implementation *asImp_() const
diff --git a/dumux/boxmodels/1p2c/1p2cpropertydefaults.hh b/dumux/boxmodels/1p2c/1p2cpropertydefaults.hh
index ce1baeae082107dcab4bcbe66744fea820a99ddf..b073b159e5e796cc065db6081416574387d449aa 100644
--- a/dumux/boxmodels/1p2c/1p2cpropertydefaults.hh
+++ b/dumux/boxmodels/1p2c/1p2cpropertydefaults.hh
@@ -41,7 +41,6 @@
#include "1p2cvolumevariables.hh"
#include "1p2cfluxvariables.hh"
#include "1p2cindices.hh"
-#include "1p2cboundaryvariables.hh"
namespace Dumux
{
@@ -69,9 +68,6 @@ SET_TYPE_PROP(BoxOnePTwoC, VolumeVariables, OnePTwoCVolumeVariables<TypeTag>);
//! define the FluxVariables
SET_TYPE_PROP(BoxOnePTwoC, FluxVariables, OnePTwoCFluxVariables<TypeTag>);
-// the BoundaryVariables property
-SET_TYPE_PROP(BoxOnePTwoC, BoundaryVariables, OnePTwoCBoundaryVariables<TypeTag>);
-
//! set default upwind weight to 1.0, i.e. fully upwind
SET_SCALAR_PROP(BoxOnePTwoC, UpwindWeight, 1.0);
diff --git a/dumux/boxmodels/2p2c/2p2clocalresidual.hh b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
index 8a23eabb40d6a5f7676ce6b6262abbde4fe1e2c9..883fe53f9a90ab73bc583b45e770a49beba276b2 100644
--- a/dumux/boxmodels/2p2c/2p2clocalresidual.hh
+++ b/dumux/boxmodels/2p2c/2p2clocalresidual.hh
@@ -213,6 +213,7 @@ class TwoPTwoCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
*
* \param flux The flux over the SCV (sub-control-volume) face for each component
* \param faceIdx The index of the SCV face
+ * \param onBoundary Evaluate flux at inner SCV face or on a boundary face
*/
void computeFlux(PrimaryVariables &flux, int faceIdx, bool onBoundary=false) const
{