diff --git a/dumux/boxmodels/2p2c/holle3p3cindices.hh b/dumux/boxmodels/2p2c/holle3p3cindices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5f674d2a3920c25bd0194304b8f9b0fe08b6e7c3
--- /dev/null
+++ b/dumux/boxmodels/2p2c/holle3p3cindices.hh
@@ -0,0 +1,98 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * 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 indices required for the holle3p3c BOX model.
+ */
+#ifndef DUMUX_HOLLE3P3C_INDICES_HH
+#define DUMUX_HOLLE3P3C_INDICES_HH
+
+namespace Dumux
+{
+/*!
+ * \ingroup HolleThreePThreeCModel
+ */
+// \{
+
+/*!
+ * \brief Enumerates the formulations which the holle3p3c model accepts.
+ */
+struct HolleThreePThreeCFormulation
+{
+ enum {
+ pgSwSn,
+ };
+};
+
+/*!
+ * \brief The indices for the isothermal HolleThreePThreeC model.
+ *
+ * \tparam formulation The formulation, only pgSwSn
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag,
+ int formulation = HolleThreePThreeCFormulation::pgSwSn,
+ int PVOffset = 0>
+class HolleThreePThreeCIndices
+{
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ // Phase indices
+ static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< Index of the water phase
+ static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< Index of the NAPL phase
+ static const int gPhaseIdx = FluidSystem::gPhaseIdx; //!< Index of the gas phase
+
+ // Component indices
+ static const int wCompIdx = 0; //!< Index of the water component
+ static const int cCompIdx = 1; //!< Index of the NAPL component
+ static const int aCompIdx = 2; //!< Index of the air component
+
+ // present phases (-> 'pseudo' primary variable)
+ static const int threePhases = 1; //!< All three phases are present
+ static const int wPhaseOnly = 2; //!< Only the water phase is present
+ static const int gnPhaseOnly = 3; //!< Only gas and NAPL phases are present
+ static const int wnPhaseOnly = 4; //!< Only water and NAPL phases are present
+ static const int gPhaseOnly = 5; //!< Only gas phase is present
+ static const int wgPhaseOnly = 6; //!< Only water and gas phases are present
+
+ // Primary variable indices
+ static const int pressureIdx = PVOffset + 0; //!< Index for gas phase pressure in a solution vector
+ static const int switch1Idx = PVOffset + 1; //!< Index 1 of saturation or mole fraction
+ static const int switch2Idx = PVOffset + 2; //!< Index 2 of saturation or mole fraction
+
+ static const int pgIdx = pressureIdx; //!< Index for gas phase pressure in a solution vector
+ static const int SOrX1Idx = switch1Idx; //!< Index of the either the saturation of the gas phase or the mass fraction secondary component if a phase is not present
+ static const int SOrX2Idx = switch2Idx; //!< Index of the either the saturation of the gas phase or the mass fraction secondary component if a phase is not present
+
+ // equation indices
+ static const int contiWEqIdx = PVOffset + wCompIdx; //!< Index of the mass conservation equation for the water component
+ static const int contiCEqIdx = PVOffset + cCompIdx; //!< Index of the mass conservation equation for the contaminant component
+ static const int contiAEqIdx = PVOffset + aCompIdx; //!< Index of the mass conservation equation for the air component
+};
+
+// \}
+
+}
+
+#endif
diff --git a/dumux/boxmodels/2p2c/holle3p3clocalresidual.hh b/dumux/boxmodels/2p2c/holle3p3clocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..86090baf9ed7ff0580fbaa9633ade5f55dbedf1b
--- /dev/null
+++ b/dumux/boxmodels/2p2c/holle3p3clocalresidual.hh
@@ -0,0 +1,336 @@
+/*****************************************************************************
+ * Copyright (C) 2011- by Holger Class *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-2009 by Bernd Flemisch *
+ * Copyright (C) 2009-2010 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 Element-wise calculation of the Jacobian matrix for problems
+ * using the two-phase two-component box model.
+ */
+
+#ifndef DUMUX_NEW_HOLLE3P3C_LOCAL_RESIDUAL_BASE_HH
+#define DUMUX_NEW_HOLLE3P3C_LOCAL_RESIDUAL_BASE_HH
+
+#include <dumux/boxmodels/common/boxmodel.hh>
+#include <dumux/common/math.hh>
+
+#include "holle3p3cproperties.hh"
+
+#include "holle3p3cvolumevariables.hh"
+
+#include "holle3p3cfluxvariables.hh"
+
+#include "holle3p3cnewtoncontroller.hh"
+
+#include <iostream>
+#include <vector>
+
+//#define VELOCITY_OUTPUT 1 // uncomment this line if an output of the velocity is needed
+
+namespace Dumux
+{
+/*!
+ * \ingroup HolleThreePThreeCModel
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the two-phase two-component box model.
+ *
+ * This class is used to fill the gaps in BoxLocalResidual for the 3P3C flow.
+ */
+template<class TypeTag>
+class HolleThreePThreeCLocalResidual: public BoxLocalResidual<TypeTag>
+{
+protected:
+ typedef HolleThreePThreeCLocalResidual<TypeTag> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalResidual)) Implementation;
+ typedef BoxLocalResidual<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementSolutionVector)) ElementSolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(BoundaryTypes)) BoundaryTypes;
+
+ typedef HolleThreePThreeCFluidState<TypeTag> FluidState;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(HolleThreePThreeCIndices)) Indices;
+
+ enum
+ {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)),
+ numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+ numComponents = GET_PROP_VALUE(TypeTag, PTAG(NumComponents)),
+
+ pressureIdx = Indices::pressureIdx,
+ switch1Idx = Indices::switch1Idx,
+ switch2Idx = Indices::switch2Idx,
+
+ contiWEqIdx = Indices::contiWEqIdx,
+ contiCEqIdx = Indices::contiCEqIdx,
+ contiAEqIdx = Indices::contiAEqIdx,
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ wCompIdx = Indices::wCompIdx,
+ cCompIdx = Indices::cCompIdx,
+ aCompIdx = Indices::aCompIdx,
+
+ lPhaseOnly = Indices::lPhaseOnly,
+ gPhaseOnly = Indices::gPhaseOnly,
+ bothPhases = Indices::bothPhases,
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementVolumeVariables)) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementBoundaryTypes)) ElementBoundaryTypes;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluxVariables)) FluxVariables;
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::ctype CoordScalar;
+
+ typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
+ typedef Dune::FieldVector<Scalar, dim> LocalPosition;
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+ typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
+
+ static const Scalar mobilityUpwindAlpha =
+ GET_PROP_VALUE(TypeTag, PTAG(MobilityUpwindAlpha));
+
+public:
+ /*!
+ * \brief Evaluate the storage term of the current solution in a
+ * single phase.
+ *
+ * \param element The element
+ * \param phaseIdx The index of the fluid phase
+ */
+ void evalPhaseStorage(const Element &element, int phaseIdx)
+ {
+ FVElementGeometry fvGeom;
+ fvGeom.update(this->gridView_(), element);
+ ElementBoundaryTypes bcTypes;
+ bcTypes.update(this->problem_(), element, fvGeom);
+ ElementVolumeVariables volVars;
+ volVars.update(this->problem_(), element, fvGeom, false);
+
+ this->residual_.resize(fvGeom.numVertices);
+ this->residual_ = 0;
+
+ this->elemPtr_ = &element;
+ this->fvElemGeomPtr_ = &fvGeom;
+ this->bcTypesPtr_ = &bcTypes;
+ this->prevVolVarsPtr_ = 0;
+ this->curVolVarsPtr_ = &volVars;
+ evalPhaseStorage_(phaseIdx);
+ }
+
+ /*!
+ * \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 &volVars = elemVolVars[scvIdx];
+
+HIERGEHTSWEITER
+
+ // compute storage term of all components within all phases
+ result = 0;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
+ result[eqIdx] += volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx)
+ * volVars.fluidState().massFrac(phaseIdx, compIdx);
+ }
+ }
+ result *= volVars.porosity();
+ }
+
+ /*!
+ * \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 for each component
+ * \param faceIdx The index of the SCV face
+ */
+ void computeFlux(PrimaryVariables &flux, int faceIdx) const
+ {
+ FluxVariables vars(this->problem_(),
+ this->elem_(),
+ this->fvElemGeom_(),
+ faceIdx,
+ this->curVolVars_());
+
+ flux = 0;
+ asImp_()->computeAdvectiveFlux(flux, vars);
+ asImp_()->computeDiffusiveFlux(flux, vars);
+ }
+
+ /*!
+ * \brief Evaluates the advective mass flux of all components 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 &vars) const
+ {
+ ////////
+ // advective fluxes of all components in all phases
+ ////////
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ // data attached to upstream and the downstream vertices
+ // of the current phase
+ const VolumeVariables &up =
+ this->curVolVars_(vars.upstreamIdx(phaseIdx));
+ const VolumeVariables &dn =
+ this->curVolVars_(vars.downstreamIdx(phaseIdx));
+
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
+ // add advective flux of current component in current
+ // phase
+ if (mobilityUpwindAlpha > 0.0)
+ // upstream vertex
+ flux[eqIdx] += vars.KmvpNormal(phaseIdx)
+ * mobilityUpwindAlpha * (up.density(phaseIdx)
+ * up.mobility(phaseIdx) * up.fluidState().massFrac(
+ phaseIdx, compIdx));
+ if (mobilityUpwindAlpha < 1.0)
+ // downstream vertex
+ flux[eqIdx] += vars.KmvpNormal(phaseIdx) * (1
+ - mobilityUpwindAlpha) * (dn.density(phaseIdx)
+ * dn.mobility(phaseIdx) * dn.fluidState().massFrac(
+ phaseIdx, compIdx));
+ }
+ }
+ }
+
+ /*!
+ * \brief Adds the diffusive mass flux of all components over
+ * a face of a subcontrol volume.
+ *
+ * \param flux The diffusive flux over the sub-control-volume face for each component
+ * \param vars The flux variables at the current SCV
+ */
+ void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &vars) const
+ {
+ // add diffusive flux of gas component in liquid phase
+ Scalar tmp =
+ - vars.porousDiffCoeff(lPhaseIdx) *
+ vars.molarDensityAtIP(lPhaseIdx) *
+ (vars.molarConcGrad(lPhaseIdx) * vars.face().normal);
+ flux[contiGEqIdx] += tmp * FluidSystem::molarMass(gCompIdx);
+ flux[contiLEqIdx] -= tmp * FluidSystem::molarMass(lCompIdx);
+
+ // add diffusive flux of liquid component in gas phase
+ tmp =
+ - vars.porousDiffCoeff(gPhaseIdx) *
+ vars.molarDensityAtIP(gPhaseIdx) *
+ (vars.molarConcGrad(gPhaseIdx) * vars.face().normal);
+ flux[contiLEqIdx] += tmp * FluidSystem::molarMass(lCompIdx);
+ flux[contiGEqIdx] -= tmp * FluidSystem::molarMass(gCompIdx);
+ }
+
+ /*!
+ * \brief Calculate the source term of the equation
+ *
+ * \param q The source/sink in the SCV for each component
+ * \param localVertexIdx The index of the SCV
+ */
+ void computeSource(PrimaryVariables &q, int localVertexIdx)
+ {
+ this->problem_().boxSDSource(q,
+ this->elem_(),
+ this->fvElemGeom_(),
+ localVertexIdx,
+ this->curVolVars_());
+ }
+
+protected:
+ void evalPhaseStorage_(int phaseIdx)
+ {
+ // evaluate the storage terms of a single phase
+ for (int i=0; i < this->fvElemGeom_().numVertices; i++) {
+ PrimaryVariables &result = this->residual_[i];
+ const ElementVolumeVariables &elemVolVars = this->curVolVars_();
+ const VolumeVariables &volVars = elemVolVars[i];
+
+ // compute storage term of all components within all phases
+ result = 0;
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ int eqIdx = (compIdx == lCompIdx) ? contiLEqIdx : contiGEqIdx;
+ result[eqIdx] += volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx)
+ * volVars.fluidState().massFrac(phaseIdx, compIdx);
+ }
+
+ result *= volVars.porosity();
+ result *= this->fvElemGeom_().subContVol[i].volume;
+ }
+ }
+
+
+ Implementation *asImp_()
+ {
+ return static_cast<Implementation *> (this);
+ }
+ const Implementation *asImp_() const
+ {
+ return static_cast<const Implementation *> (this);
+ }
+
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/boxmodels/2p2c/holle3p3cnewtoncontroller.hh b/dumux/boxmodels/2p2c/holle3p3cnewtoncontroller.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1453d0d2598d4c3b6c0bcddd160cb447e7a2887f
--- /dev/null
+++ b/dumux/boxmodels/2p2c/holle3p3cnewtoncontroller.hh
@@ -0,0 +1,181 @@
+/*****************************************************************************
+ * Copyright (C) 2008 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 holle3p3c 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_HOLLE3P3C_NEWTON_CONTROLLER_HH
+#define DUMUX_HOLLE3P3C_NEWTON_CONTROLLER_HH
+
+#include "holle3p3cproperties.hh"
+
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+namespace Dumux {
+/*!
+ * \ingroup HolleThreePThreeCModel
+ * \brief A holle3p3c 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 HolleThreePThreeCNewtonController : public NewtonController<TypeTag>
+{
+ typedef NewtonController<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonController)) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Model;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(NewtonMethod)) NewtonMethod;
+
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(HolleThreePThreeCIndices)) Indices;
+
+ enum {
+ pressureIdx = Indices::pressureIdx,
+ switch1Idx = Indices::switch1Idx
+ switch2Idx = Indices::switch2Idx
+ };
+
+ enum { enablePartialReassemble = GET_PROP_VALUE(TypeTag, PTAG(EnablePartialReassemble)) };
+
+public:
+ HolleThreePThreeCNewtonController()
+ {
+ this->setRelTolerance(1e-7);
+ this->setTargetSteps(9);
+ this->setMaxSteps(18);
+ };
+
+
+ /*!
+ * \brief
+ * Suggest a new time step size based either on the number of newton
+ * iterations required or on the variable switch
+ *
+ * \param u The current global solution vector
+ * \param uLastIter The previous global solution vector
+ *
+ */
+ void newtonEndStep(SolutionVector &uCurrentIter,
+ const SolutionVector &uLastIter)
+ {
+ // call the method of the base class
+ this->method().model().updateStaticData(uCurrentIter, uLastIter);
+ ParentType::newtonEndStep(uCurrentIter, uLastIter);
+ }
+
+ /*!
+ * \brief Update the current solution function with a delta vector.
+ *
+ * The error estimates required for the newtonConverged() and
+ * newtonProceed() methods should be updated here.
+ *
+ * Different update strategies, such as line search and chopped
+ * updates can be implemented. The default behaviour is just to
+ * subtract deltaU from uLastIter.
+ *
+ * \param uCurrentIter The solution after the current Newton iteration
+ * \param uLastIter The solution after the last Newton iteration
+ * \param deltaU The delta as calculated from solving the linear
+ * system of equations. This parameter also stores
+ * the updated solution.
+ */
+ void newtonUpdate(SolutionVector &uCurrentIter,
+ const SolutionVector &uLastIter,
+ const SolutionVector &deltaU)
+ {
+ this->writeConvergence_(uLastIter, deltaU);
+
+ this->newtonUpdateRelError(uLastIter, deltaU);
+
+ // compute the vertex and element colors for partial
+ // reassembly
+ if (enablePartialReassemble) {
+ Scalar reassembleTol = Dumux::geometricMean(this->error_, 0.1*this->tolerance_);
+ reassembleTol = std::max(reassembleTol, 0.1*this->tolerance_);
+ this->model_().jacobianAssembler().updateDiscrepancy(uLastIter, deltaU);
+ this->model_().jacobianAssembler().computeColors(reassembleTol);
+ }
+
+ if (GET_PROP_VALUE(TypeTag, PTAG(NewtonUseLineSearch)))
+ lineSearchUpdate_(uCurrentIter, uLastIter, deltaU);
+ else {
+ uCurrentIter = uLastIter;
+ uCurrentIter -= deltaU;
+ }
+ }
+
+
+ /*!
+ * \brief
+ * Returns true if the current solution can be considered to
+ * be accurate enough
+ */
+ bool newtonConverged()
+ {
+ if (this->method().model().switched())
+ return false;
+
+ return ParentType::newtonConverged();
+ };
+
+private:
+ void lineSearchUpdate_(SolutionVector &uCurrentIter,
+ const SolutionVector &uLastIter,
+ const SolutionVector &deltaU)
+ {
+ Scalar lambda = 1.0;
+ Scalar globDef;
+
+ // calculate the residual of the current solution
+ SolutionVector tmp(uLastIter);
+ Scalar oldGlobDef = this->method().model().globalResidual(tmp, uLastIter);
+
+ while (true) {
+ uCurrentIter = deltaU;
+ uCurrentIter *= -lambda;
+ uCurrentIter += uLastIter;
+
+ // calculate the residual of the current solution
+ globDef = this->method().model().globalResidual(tmp, uCurrentIter);
+
+ if (globDef < oldGlobDef || lambda <= 1.0/8) {
+ this->endIterMsg() << ", defect " << oldGlobDef << "->" << globDef << "@lambda=" << lambda;
+ return;
+ }
+
+ // try with a smaller update
+ lambda /= 2;
+ }
+ };
+
+};
+}
+
+#endif
diff --git a/dumux/boxmodels/2p2c/holle3p3cproblem.hh b/dumux/boxmodels/2p2c/holle3p3cproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f4a0f34bd59067d3a271b1c993df5b95ad47ce70
--- /dev/null
+++ b/dumux/boxmodels/2p2c/holle3p3cproblem.hh
@@ -0,0 +1,128 @@
+/*****************************************************************************
+ * 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 two-phase,
+ * two-component box model
+ */
+#ifndef DUMUX_HOLLE3p3C_PROBLEM_HH
+#define DUMUX_HOLLE3p3C_PROBLEM_HH
+
+#include <dumux/boxmodels/common/boxproblem.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup HolleThreePThreeCModel
+ * \brief Base class for all problems which use the two-phase, two-component box model
+ *
+ */
+template<class TypeTag>
+class HolleThreePThreeCProblem : public BoxProblem<TypeTag>
+{
+ typedef BoxProblem<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Model)) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GridView::Grid Grid;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(TimeManager)) TimeManager;
+
+ // material properties
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+
+ enum {
+ dim = Grid::dimension,
+ dimWorld = Grid::dimensionworld
+ };
+
+ typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+ /*!
+ * \brief The constructor
+ *
+ * \param timeManager The time manager
+ * \param gridView The grid view
+ */
+ HolleThreePThreeCProblem(TimeManager &timeManager, const GridView &gridView)
+ : ParentType(timeManager, gridView),
+ gravity_(0),
+ spatialParams_(gridView)
+ {
+ if (GET_PROP_VALUE(TypeTag, PTAG(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
+ { DUNE_THROW(Dune::NotImplemented, "The Problem must implement a temperature() method for isothermal problems!"); };
+
+ /*!
+ * \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 Returns the spatial parameters object.
+ */
+ SpatialParameters &spatialParameters()
+ { return spatialParams_; }
+
+ /*!
+ * \brief Returns the spatial parameters object.
+ */
+ const SpatialParameters &spatialParameters() const
+ { return spatialParams_; }
+
+ // \}
+
+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_;
+
+ // spatial parameters
+ SpatialParameters spatialParams_;
+};
+
+}
+
+#endif
diff --git a/dumux/boxmodels/2p2c/holle3p3cproperties.hh b/dumux/boxmodels/2p2c/holle3p3cproperties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..739206fffa1ba28cab93cf9eee164d3108934e20
--- /dev/null
+++ b/dumux/boxmodels/2p2c/holle3p3cproperties.hh
@@ -0,0 +1,68 @@
+/*****************************************************************************
+ * Copyright (C) 2011- by Holger Class *
+ * Copyright (C) 2008-2010 by Andreas Lauser *
+ * Copyright (C) 2008-2009 by Melanie Darcis *
+ * Copyright (C) 2008-2009 by Klaus Mosthaf *
+ * Copyright (C) 2008-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/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup HolleThreePThreeCModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines the properties required for the holle3p3c BOX model.
+ */
+#ifndef DUMUX_HOLLE3P3C_PROPERTIES_HH
+#define DUMUX_HOLLE3P3C_PROPERTIES_HH
+
+#include <dumux/boxmodels/common/boxproperties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the isothermal single phase problems
+NEW_TYPE_TAG(BoxHolleThreePThreeC, INHERITS_FROM(BoxModel));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
+NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system
+NEW_PROP_TAG(HolleThreePThreeCIndices); //!< Enumerations for the holle3p3c models
+NEW_PROP_TAG(Formulation); //!< The formulation of the model
+NEW_PROP_TAG(SpatialParameters); //!< The type of the spatial parameters
+NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
+
+NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters)
+NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extracted from the spatial parameters)
+
+NEW_PROP_TAG(EnableGravity); //!< Returns whether gravity is considered in the problem
+NEW_PROP_TAG(MobilityUpwindAlpha); //!< The value of the upwind parameter for the mobility
+}
+}
+
+#endif
diff --git a/dumux/boxmodels/2p2c/holle3p3cpropertydefaults.hh b/dumux/boxmodels/2p2c/holle3p3cpropertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..5a8d4ffedc0e72d1442b92d098c2cfe6d145417a
--- /dev/null
+++ b/dumux/boxmodels/2p2c/holle3p3cpropertydefaults.hh
@@ -0,0 +1,165 @@
+/*****************************************************************************
+ * Copyright (C) 2011 by Holger Class *
+ * 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/>. *
+ *****************************************************************************/
+/*!
+ * \ingroup HolleThreePThreeCModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines default values for most properties required by the
+ * holle3p3c box model.
+ */
+#ifndef DUMUX_HOLLE2P2C_PROPERTY_DEFAULTS_HH
+#define DUMUX_HOLLE2P2C_PROPERTY_DEFAULTS_HH
+
+#include "holle3p3cindices.hh"
+
+#include "holle3p3cmodel.hh"
+#include "holle3p3cproblem.hh"
+#include "holle3p3cindices.hh"
+#include "holle3p3cfluxvariables.hh"
+#include "holle3p3cvolumevariables.hh"
+#include "holle3p3cfluidstate.hh"
+#include "holle3p3cproperties.hh"
+#include "holle3p3cnewtoncontroller.hh"
+// #include "holle3p3cboundaryvariables.hh"
+
+namespace Dumux
+{
+
+namespace Properties {
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+
+/*!
+ * \brief Set the property for the number of components.
+ *
+ * We just forward the number from the fluid system and use an static
+ * assert to make sure it is 2.
+ */
+SET_PROP(BoxHolleThreePThreeC, NumComponents)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numComponents;
+
+ static_assert(value == 3,
+ "Only fluid systems with 3 components are supported by the 3p3c model!");
+};
+
+/*!
+ * \brief Set the property for the number of fluid phases.
+ *
+ * We just forward the number from the fluid system and use an static
+ * assert to make sure it is 2.
+ */
+SET_PROP(BoxHolleThreePThreeC, NumPhases)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numPhases;
+ static_assert(value == 3,
+ "Only fluid systems with 3 phases are supported by the 3p3c model!");
+};
+
+SET_INT_PROP(BoxHolleThreePThreeC, NumEq, 3); //!< set the number of equations to 2
+
+//! Set the default formulation to pg-Sw-Sn
+SET_INT_PROP(BoxHolleThreePThreeC,
+ Formulation,
+ HolleThreePThreeCFormulation::pgSwSn);
+
+/*!
+ * \brief Set the property for the material law by retrieving it from
+ * the spatial parameters.
+ */
+SET_PROP(BoxHolleThreePThreeC, MaterialLaw)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SpatialParameters)) SpatialParameters;
+
+public:
+ typedef typename SpatialParameters::MaterialLaw type;
+};
+
+/*!
+ * \brief Set the property for the material parameters by extracting
+ * it from the material law.
+ */
+SET_PROP(BoxHolleThreePThreeC, MaterialLawParams)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
+
+public:
+ typedef typename MaterialLaw::Params type;
+};
+
+//! Use the holle3p3c local jacobian operator for the holle3p3c model
+SET_TYPE_PROP(BoxHolleThreePThreeC,
+ LocalResidual,
+ HolleThreePThreeCLocalResidual<TypeTag>);
+
+//! Use the holle3p3c specific newton controller for the holle3p3c model
+SET_TYPE_PROP(BoxHolleThreePThreeC, NewtonController, HolleThreePThreeCNewtonController<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(BoxHolleThreePThreeC, Model, HolleThreePThreeCModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxHolleThreePThreeC, VolumeVariables, HolleThreePThreeCVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxHolleThreePThreeC, FluxVariables, HolleThreePThreeCFluxVariables<TypeTag>);
+
+//! the BoundaryVariables property
+// SET_TYPE_PROP(BoxHolleThreePThreeC, BoundaryVariables, HolleThreePThreeCBoundaryVariables<TypeTag>);
+
+//! the upwind factor for the mobility.
+SET_SCALAR_PROP(BoxHolleThreePThreeC, MobilityUpwindAlpha, 1.0);
+
+//! The indices required by the isothermal holle3p3c model
+SET_PROP(BoxHolleThreePThreeC,
+ HolleThreePThreeCIndices)
+{ private:
+ enum { Formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation)) };
+public:
+ typedef HolleThreePThreeCIndices<TypeTag, Formulation, 0> type;
+};
+
+// enable jacobian matrix recycling by default
+SET_BOOL_PROP(BoxHolleThreePThreeC, EnableJacobianRecycling, true);
+// enable partial reassembling by default
+SET_BOOL_PROP(BoxHolleThreePThreeC, EnablePartialReassemble, true);
+// enable time-step ramp up by default
+SET_BOOL_PROP(BoxHolleThreePThreeC, EnableTimeStepRampUp, false);
+
+//
+}
+
+}
+
+#endif
diff --git a/dumux/boxmodels/2p2c/holle3p3cvolumevariables.hh b/dumux/boxmodels/2p2c/holle3p3cvolumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7d87e3b2c8bb20b3e4f0c9ea4f558e466de3598a
--- /dev/null
+++ b/dumux/boxmodels/2p2c/holle3p3cvolumevariables.hh
@@ -0,0 +1,302 @@
+/*****************************************************************************
+ * Copyright (C) 2011- by Holger Class *
+ * Copyright (C) 2008,2009 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 Contains the quantities which are constant within a
+ * finite volume in the two-phase, two-component model.
+ */
+#ifndef DUMUX_HOLLE3P3C_VOLUME_VARIABLES_HH
+#define DUMUX_HOLLE3P3C_VOLUME_VARIABLES_HH
+
+#include <dumux/boxmodels/common/boxmodel.hh>
+#include <dumux/common/math.hh>
+
+#include <dune/common/collectivecommunication.hh>
+#include <vector>
+#include <iostream>
+
+#include "holle3p3cproperties.hh"
+#include "holle3p3cfluidstate.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup HolleThreePThreeCModel
+ * \brief Contains the quantities which are are constant within a
+ * finite volume in the two-phase, two-component model.
+ */
+template <class TypeTag>
+class HolleThreePThreeCVolumeVariables : public BoxVolumeVariables<TypeTag>
+{
+ typedef BoxVolumeVariables<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(VolumeVariables)) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(PrimaryVariables)) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLawParams)) MaterialLawParams;
+
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(HolleThreePThreeCIndices)) Indices;
+ enum {
+ dim = GridView::dimension,
+
+ numPhases = GET_PROP_VALUE(TypeTag, PTAG(NumPhases)),
+ formulation = GET_PROP_VALUE(TypeTag, PTAG(Formulation)),
+
+ wCompIdx = Indices::wCompIdx,
+ aCompIdx = Indices::aCompIdx,
+ cCompIdx = Indices::cCompIdx,
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx
+ nPhaseIdx = Indices::nPhaseIdx
+ };
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef HolleThreePThreeCFluidState<TypeTag> FluidState;
+
+public:
+ /*!
+ * \brief Update all quantities for a given control volume.
+ *
+ * \param priVars The primary variables
+ * \param problem The problem
+ * \param element The element
+ * \param elemGeom The finite-volume geometry in the box scheme
+ * \param scvIdx The local index of the SCV (sub-control volume)
+ * \param isOldSol Evaluate function with solution of current or previous time step
+ */
+ void update(const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ bool isOldSol)
+ {
+ ParentType::update(priVars,
+ problem,
+ element,
+ elemGeom,
+ scvIdx,
+ isOldSol);
+
+ asImp().updateTemperature_(priVars,
+ element,
+ elemGeom,
+ scvIdx,
+ problem);
+
+ // capillary pressure parameters
+ const MaterialLawParams &materialParams =
+ problem.spatialParameters().materialLawParams(element, elemGeom, scvIdx);
+
+ int globalVertIdx = problem.model().dofMapper().map(element, scvIdx, dim);
+ int phasePresence = problem.model().phasePresence(globalVertIdx, isOldSol);
+
+ // calculate phase state
+ fluidState_.update(priVars, materialParams, temperature(), phasePresence);
+ Valgrind::CheckDefined(fluidState_);
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ // Mobilities
+ const Scalar mu =
+ FluidSystem::phaseViscosity(phaseIdx,
+ fluidState().temperature(),
+ fluidState().phasePressure(phaseIdx),
+ fluidState());
+ Scalar kr;
+ if (phaseIdx == wPhaseIdx)
+ kr = MaterialLaw::krw(materialParams, fluidState.saturation());
+ else if (phaseIdx == nPhaseIdx)
+ kr = MaterialLaw::krn(materialParams, fluidState.saturation());
+ else // (phaseIdx == gPhaseIdx)
+ kr = MaterialLaw::krg(materialParams, fluidState.saturation());
+ mobility_[phaseIdx] = kr / mu;
+ Valgrind::CheckDefined(mobility_[phaseIdx]);
+ }
+
+/* ACHTUNG Die Umrechnung in einen effektiven Diffusionskoeffizienten im Dreikomponenten-System erfolgt an anderer Stelle im fluidSystem */
+
+ // diffusivity coefficents
+ diffCoeff_[gPhaseIdx][wCompIdx] =
+ FluidSystem::diffCoeff(gphaseIdx,
+ wCompIdx,
+ fluidState_.temperature(),
+ fluidState_.phasePressure(gphaseIdx),
+ fluidState_);
+ diffCoeff_[gPhaseIdx][cCompIdx] =
+ FluidSystem::diffCoeff(gphaseIdx,
+ cCompIdx,
+ fluidState_.temperature(),
+ fluidState_.phasePressure(gphaseIdx),
+ fluidState_);
+ diffCoeff_[gPhaseIdx][aCompIdx] = 0.0; // dummy, should not be used !
+
+ /* no diffusion in water phase considered at the moment */
+ diffCoeff_[wPhaseIdx][aCompIdx] =
+ FluidSystem::diffCoeff(wphaseIdx,
+ aCompIdx,
+ fluidState_.temperature(),
+ fluidState_.phasePressure(wphaseIdx),
+ fluidState_);
+ diffCoeff_[gPhaseIdx][cCompIdx] =
+ FluidSystem::diffCoeff(wphaseIdx,
+ cCompIdx,
+ fluidState_.temperature(),
+ fluidState_.phasePressure(wphaseIdx),
+ fluidState_);
+ diffCoeff_[wPhaseIdx][wCompIdx] = 0.0; // dummy, should not be used !
+
+ /* no diffusion in NAPL phase considered at the moment */
+ diffCoeff_[nPhaseIdx][cCompIdx] = 0.0;
+ diffCoeff_[nPhaseIdx][wCompIdx] = 0.0;
+ diffCoeff_[nPhaseIdx][aCompIdx] = 0.0;
+
+
+ Valgrind::CheckDefined(diffCoeff_);
+
+ // porosity
+ porosity_ = problem.spatialParameters().porosity(element,
+ elemGeom,
+ scvIdx);
+ Valgrind::CheckDefined(porosity_);
+ }
+
+ /*!
+ * \brief Returns the phase state for the control-volume.
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+
+ /*!
+ * \brief Returns the effective saturation of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar saturation(int phaseIdx) const
+ { return fluidState_.saturation(phaseIdx); }
+
+ /*!
+ * \brief Returns the mass density of a given phase within the
+ * control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar density(int phaseIdx) const
+ { return fluidState_.density(phaseIdx); }
+
+ /*!
+ * \brief Returns the molar density of a given phase within the
+ * control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar molarDensity(int phaseIdx) const
+ { return fluidState_.density(phaseIdx) / fluidState_.averageMolarMass(phaseIdx); }
+
+ /*!
+ * \brief Returns the effective pressure of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar pressure(int phaseIdx) const
+ { return fluidState_.phasePressure(phaseIdx); }
+
+ /*!
+ * \brief Returns temperature inside the sub-control volume.
+ *
+ * Note that we assume thermodynamic equilibrium, i.e. the
+ * temperatures of the rock matrix and of all fluid phases are
+ * identical.
+ */
+ Scalar temperature() const
+ { return temperature_; }
+
+ /*!
+ * \brief Returns the effective mobility of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar mobility(int phaseIdx) const
+ {
+ return mobility_[phaseIdx];
+ }
+
+ /*!
+ * \brief Returns the effective capillary pressure within the control volume.
+ */
+ Scalar capillaryPressure() const
+ { return fluidState_.capillaryPressure(); }
+
+ /*!
+ * \brief Returns the average porosity within the control volume.
+ */
+ Scalar porosity() const
+ { return porosity_; }
+
+ /*!
+ * \brief Returns the diffusivity coefficient matrix
+ */
+ Dune::FieldMatrix<Scalar, numPhases, 3> diffCoeff(int phaseIdx) const
+ { return diffCoeff_; }
+
+
+protected:
+
+ void updateTemperature_(const PrimaryVariables &priVars,
+ const Element &element,
+ const FVElementGeometry &elemGeom,
+ int scvIdx,
+ const Problem &problem)
+ {
+ temperature_ = problem.temperature(element, elemGeom, scvIdx);
+ }
+
+ Scalar temperature_; //!< Temperature within the control volume
+ Scalar porosity_; //!< Effective porosity within the control volume
+ Scalar mobility_[numPhases]; //!< Effective mobility within the control volume
+ /* We need a tensor here !! */
+ //!< Binary diffusion coefficients of the 3 components in the phases
+ Dune::FieldMatrix<Scalar, numPhases, 3> diffCoeff_;
+ FluidState fluidState_;
+
+private:
+ Implementation &asImp()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp() const
+ { return *static_cast<const Implementation*>(this); }
+};
+
+} // end namepace
+
+#endif