New model 2pncmin moved from dumux-devel to dumux-stable. The changes are...

New model 2pncmin moved from dumux-devel to dumux-stable. The changes are reviewed by J.Hommel. The model was directly or indirectly used in publications (Jambhekar et. al 2015 and Hommel. et.al 2015) 

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@15137 2fb0f335-1f38-0410-981e-8018bf24f1b0

dumux/implicit/2pncmin/2pncminfluxvariables.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   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 data which is required to calculate
+ *        all fluxes of components over a face of a finite volume for
+ *        the two-phase two-component mineralization model fully implicit model.
+ */
+#ifndef DUMUX_2PNCMIN_FLUX_VARIABLES_HH
+#define DUMUX_2PNCMIN_FLUX_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+#include <dumux/common/spline.hh>
+#include <dumux/implicit/2pnc/2pncfluxvariables.hh>
+#include "2pncminproperties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPNCMinModel
+ * \ingroup ImplicitFluxVariables
+ * \brief Contains the data which is required to calculate
+ *        all fluxes of components over a face of a finite volume for
+ *        the two-phase n-component mineralization fully implicit model.
+ *
+ * This means pressure and concentration gradients, phase densities at
+ * the integration point, etc.
+ */
+
+template <class TypeTag>
+class TwoPNCMinFluxVariables : public TwoPNCFluxVariables<TypeTag>
+{
+    typedef TwoPNCFluxVariables<TypeTag> ParentType;
+    typedef TwoPNCMinFluxVariables<TypeTag> ThisType;
+    
+    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 GridView::ctype CoordScalar;
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+
+    enum {
+        dim = GridView::dimension,
+        dimWorld = GridView::dimensionworld,
+        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+        numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+    };
+
+    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+    typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename FVElementGeometry::SubControlVolume SCV;
+    typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
+
+    typedef Dune::FieldVector<CoordScalar, dimWorld> DimVector;
+    typedef Dune::FieldMatrix<CoordScalar, dim, dim> DimMatrix;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    enum {
+        wPhaseIdx = FluidSystem::wPhaseIdx,
+        nPhaseIdx = FluidSystem::nPhaseIdx,
+        wCompIdx  = FluidSystem::wCompIdx,
+    };
+
+public:
+    /*!
+     * \brief The constructor
+     *
+     * \param problem The problem
+     * \param element The finite element
+     * \param fvGeometry The finite-volume geometry in the fully implicit scheme
+     * \param fIdx The local index of the sub-control-volume face
+     * \param elemVolVars The volume variables of the current element
+     * \param onBoundary Evaluate flux at inner sub-control-volume face or on a boundary face
+     */
+    TwoPNCMinFluxVariables(const Problem &problem,
+                     const Element &element,
+                     const FVElementGeometry &fvGeometry,
+                     const int fIdx,
+                     const ElementVolumeVariables &elemVolVars,
+                     const bool onBoundary = false)
+    : ParentType(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary) 
+    {
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+            this->density_[phaseIdx] = Scalar(0);
+            this->molarDensity_[phaseIdx] = Scalar(0);
+            this->potentialGrad_[phaseIdx] = Scalar(0);
+            for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+            {
+            	this->massFractionGrad_[phaseIdx][compIdx] = Scalar(0);
+            	this->moleFractionGrad_[phaseIdx][compIdx] = Scalar(0);
+            }
+        }
+        this->calculateGradients_(problem, element, elemVolVars);
+        this->calculateVelocities_(problem, element, elemVolVars);
+        this->calculateporousDiffCoeff_(problem, element, elemVolVars);
+    };
+
+protected:    
+    void calculateVelocities_(const Problem &problem,
+                              const Element &element,
+                              const ElementVolumeVariables &elemVolVars)
+    {
+        const SpatialParams &spatialParams = problem.spatialParams();
+        // multiply the pressure potential with the intrinsic permeability
+        DimMatrix K(0.0);
+
+        for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
+        {
+            const VolumeVariables &volVarsI = elemVolVars[this->face().i];
+            const VolumeVariables &volVarsJ = elemVolVars[this->face().j];
+
+            auto K_i = spatialParams.intrinsicPermeability(element,this->fvGeometry_,this->face().i);
+            K_i *= volVarsI.permFactor();
+
+            auto K_j = spatialParams.intrinsicPermeability(element,this->fvGeometry_,this->face().j);
+            K_j *= volVarsJ.permFactor();
+
+            spatialParams.meanK(K,K_i,K_j);
+
+            K.mv(this->potentialGrad_[phaseIdx], this->Kmvp_[phaseIdx]);
+            this->KmvpNormal_[phaseIdx] = - (this->Kmvp_[phaseIdx] * this->face().normal);
+        }
+
+        // set the upstream and downstream vertices
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+        {
+            this->upstreamIdx_[phaseIdx] = this->face().i;
+            this->downstreamIdx_[phaseIdx] = this->face().j;
+
+            if (this->KmvpNormal_[phaseIdx] < 0) {
+                std::swap(this->upstreamIdx_[phaseIdx],
+                          this->downstreamIdx_[phaseIdx]);
+            }
+        }
+    }
+};
+
+} // end namespace
+
+#endif

dumux/implicit/2pncmin/2pncminindices.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+
+/*!
+ * \file
+ * \brief Defines the indices required for the two-phase n-component mineralization 
+ *        fully implicit model.
+ */
+#ifndef DUMUX_2PNCMIN_INDICES_HH
+#define DUMUX_2PNCMIN_INDICES_HH
+
+#include <dumux/implicit/2pnc/2pncindices.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPNCMinModel
+ * \ingroup ImplicitIndices
+ * \brief The indices for the isothermal two-phase n-component mineralization model.
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int PVOffset = 0>
+	class TwoPNCMinIndices: public TwoPNCIndices<TypeTag, PVOffset>
+{ 
+};
+
+// \}
+
+}
+
+#endif

dumux/implicit/2pncmin/2pncminlocalresidual.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   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 n-component mineralisation box model.
+ */
+
+#ifndef DUMUX_2PNCMIN_LOCAL_RESIDUAL_BASE_HH
+#define DUMUX_2PNCMIN_LOCAL_RESIDUAL_BASE_HH
+
+#include "2pncminproperties.hh"
+#include <dumux/implicit/2pnc/2pnclocalresidual.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPNCMinModel
+ * \ingroup ImplicitLocalResidual
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ *        using the two-phase n-component mineralization fully implicit box model.
+ *
+ * This class is used to fill the gaps in ImplicitLocalResidual for the two-phase n-component flow.
+ */
+template<class TypeTag>
+class TwoPNCMinLocalResidual: public TwoPNCLocalResidual<TypeTag>
+{
+protected:
+	typedef TwoPNCLocalResidual<TypeTag> ParentType;
+    typedef TwoPNCMinLocalResidual<TypeTag> ThisType; 
+	typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementSolutionVector) ElementSolutionVector;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
+
+
+    enum
+    {
+        numEq = GET_PROP_VALUE(TypeTag, NumEq),
+        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+        numSPhases = GET_PROP_VALUE(TypeTag, NumSPhases),
+        numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+
+        replaceCompEqIdx = GET_PROP_VALUE(TypeTag, ReplaceCompEqIdx),
+
+        pressureIdx = Indices::pressureIdx,
+        switchIdx = Indices::switchIdx,
+
+        wPhaseIdx = FluidSystem::wPhaseIdx,
+        nPhaseIdx = FluidSystem::nPhaseIdx,
+
+        wCompIdx = FluidSystem::wCompIdx,
+        nCompIdx = FluidSystem::nCompIdx,
+
+        conti0EqIdx = Indices::conti0EqIdx,
+
+        wPhaseOnly = Indices::wPhaseOnly,
+        nPhaseOnly = Indices::nPhaseOnly,
+        bothPhases = Indices::bothPhases,
+
+        plSg = TwoPNCFormulation::plSg,
+        pgSl = TwoPNCFormulation::pgSl,
+        formulation = GET_PROP_VALUE(TypeTag, Formulation)
+    };
+
+    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+    typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+public:
+    /*!
+     * \brief Constructor. Sets the upwind weight.
+     */
+    TwoPNCMinLocalResidual()
+    {
+        // 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
+        this->massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
+    };
+
+    /*!
+     * \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).
+     * In contrast to the 2pnc model, here, the storage of solid phases is included too.
+     *
+     *  \param storage 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 &storage, int scvIdx, bool usePrevSol) const
+  {
+      //call parenttype function
+      ParentType::computeStorage(storage, scvIdx, usePrevSol);
+      
+      const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_()
+      : this->curVolVars_();
+    const VolumeVariables &volVars = elemVolVars[scvIdx];
+
+    // Compute storage term of all solid (precipitated) phases (excluding the non-reactive matrix)
+    for (int phaseIdx = numPhases; phaseIdx < numPhases + numSPhases; ++phaseIdx)
+    {
+      int eqIdx = conti0EqIdx + numComponents-numPhases + phaseIdx;
+      storage[eqIdx] += volVars.precipitateVolumeFraction(phaseIdx)*volVars.molarDensity(phaseIdx);
+    }
+
+      Valgrind::CheckDefined(storage);
+  }
+};
+} // end namespace
+
+#endif

dumux/implicit/2pncmin/2pncminmodel.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   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 fully implicit box scheme to the two-phase n-component flow model.
+*/
+
+#ifndef DUMUX_2PNCMIN_MODEL_HH
+#define DUMUX_2PNCMIN_MODEL_HH
+
+#include "2pncminproperties.hh"
+#include "2pncminindices.hh"
+
+#include <dumux/material/constants.hh>
+ #include <dumux/implicit/2pnc/2pncmodel.hh>
+#include "2pncminlocalresidual.hh"
+#include <dumux/implicit/common/implicitvelocityoutput.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPNCMinModel
+ * \brief Adaption of the fully implicit scheme to the
+ *        two-phase n-component fully implicit model.
+ *
+ * This model implements two-phase n-component flow of two compressible and
+ * partially miscible fluids \f$\alpha \in \{ w, n \}\f$ composed of the n components
+ * \f$\kappa \in \{ w, a,\cdots \}\f$. The standard multiphase Darcy
+ * approach is used as the equation for the conservation of momentum:
+ * \f[
+ v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mbox{\bf K}
+ \left(\text{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right)
+ * \f]
+ *
+ * By inserting this into the equations for the conservation of the
+ * components, one gets one transport equation for each component
+ * \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\}
+ \nonumber \\ \nonumber \\
+    &-& \sum_\alpha \text{div} \left\{{\bf D_{\alpha, pm}^\kappa} \varrho_{\alpha} \text{grad}\, X^\kappa_{\alpha} \right\}
+ - \sum_\alpha q_\alpha^\kappa = 0 \qquad \kappa \in \{w, a,\cdots \} \, ,
+ \alpha \in \{w, g\}
+ \f}
+ *
+ * All equations are discretized using a vertex-centered finite volume (box)
+ * or cell-centered finite volume scheme (this is not done for 2pnc approach yet, however possible) as  
+ * spatial and the implicit Euler method as time discretization.
+ *
+ * By using constitutive relations for the capillary pressure \f$p_c =
+ * p_n - p_w\f$ and relative permeability \f$k_{r\alpha}\f$ and taking
+ * advantage of the fact that \f$S_w + S_n = 1\f$ and \f$X^\kappa_w + X^\kappa_n = 1\f$, the number of
+ * unknowns can be reduced to number of components.
+ *
+ * The used primary variables are, like in the two-phase model, either \f$p_w\f$ and \f$S_n\f$
+ * or \f$p_n\f$ and \f$S_w\f$. The formulation which ought to be used can be
+ * specified by setting the <tt>Formulation</tt> property to either
+ * TwoPTwoCIndices::pWsN or TwoPTwoCIndices::pNsW. By
+ * default, the model uses \f$p_w\f$ and \f$S_n\f$.
+ *
+ * Moreover, the second primary variable depends on the phase state, since a
+ * primary variable switch is included. The phase state is stored for all nodes
+ * of the system. The model is uses mole fractions.
+ *Following cases can be distinguished:
+ * <ul>
+ *  <li> Both phases are present: The saturation is used (either \f$S_n\f$ or \f$S_w\f$, dependent on the chosen <tt>Formulation</tt>),
+ *      as long as \f$ 0 < S_\alpha < 1\f$</li>.
+ *  <li> Only wetting phase is present: The mass fraction of, e.g., air in the wetting phase \f$X^a_w\f$ is used,
+ *      as long as the maximum mass fraction is not exceeded (\f$X^a_w<X^a_{w,max}\f$)</li>
+ *  <li> Only non-wetting phase is present: The mass fraction of, e.g., water in the non-wetting phase, \f$X^w_n\f$, is used,
+ *      as long as the maximum mass fraction is not exceeded (\f$X^w_n<X^w_{n,max}\f$)</li>
+ * </ul>
+ */
+
+template<class TypeTag>
+class TwoPNCMinModel: public TwoPNCModel<TypeTag>
+{
+    typedef TwoPNCMinModel<TypeTag> ThisType;
+    typedef TwoPNCModel<TypeTag> ParentType;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+    typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, VertexMapper) VertexMapper;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementMapper) ElementMapper;
+    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    typedef Dumux::Constants<Scalar> Constant;
+
+    enum {
+        dim = GridView::dimension,
+        dimWorld = GridView::dimensionworld,
+
+        numEq = GET_PROP_VALUE(TypeTag, NumEq),
+        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+        numSPhases = GET_PROP_VALUE(TypeTag, NumSPhases),
+        numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+        numSecComponents = GET_PROP_VALUE(TypeTag, NumSecComponents),
+        numMajorComponents = GET_PROP_VALUE(TypeTag, NumMajorComponents),
+
+        pressureIdx = Indices::pressureIdx,
+        switchIdx = Indices::switchIdx,
+
+        wPhaseIdx = Indices::wPhaseIdx,
+        nPhaseIdx = Indices::nPhaseIdx,
+
+        wCompIdx = FluidSystem::wCompIdx,
+        nCompIdx = FluidSystem::nCompIdx,
+
+        wPhaseOnly = Indices::wPhaseOnly,
+        nPhaseOnly = Indices::nPhaseOnly,
+        bothPhases = Indices::bothPhases,
+
+        plSg = TwoPNCFormulation::plSg,
+        pgSl = TwoPNCFormulation::pgSl,
+        formulation = GET_PROP_VALUE(TypeTag, Formulation),
+		useElectrochem = GET_PROP_VALUE(TypeTag, useElectrochem)
+    };
+
+    typedef typename GridView::template Codim<dim>::Entity Vertex;
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+    typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
+
+    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+    typedef typename GridView::ctype CoordScalar;
+    typedef Dune::FieldMatrix<CoordScalar, dimWorld, dimWorld> Tensor;
+    typedef Dune::FieldVector<Scalar, numPhases> PhasesVector;
+
+    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+        enum { dofCodim = isBox ? dim : 0 };
+
+public:
+
+    /*!
+     * \brief Append all quantities of interest which can be derived
+     *        from the solution of the current time step to the VTK
+     *        writer.
+     *
+     * \param sol The solution vector
+     * \param writer The writer for multi-file VTK datasets
+     */
+    template<class MultiWriter>
+    //additional output of the permeability and the precipitate volume fractions
+    void addOutputVtkFields(const SolutionVector &sol,
+                            MultiWriter &writer)
+    {
+        typedef Dune::BlockVector<Dune::FieldVector<Scalar, 1> > ScalarField;
+        typedef Dune::BlockVector<Dune::FieldVector<double, dim> > VectorField;
+
+        // get the number of degrees of freedom
+        unsigned numDofs = this->numDofs();
+
+        // create the required scalar fields
+
+        ScalarField *Sg		   = writer.allocateManagedBuffer (numDofs);
+        ScalarField *Sl 	   = writer.allocateManagedBuffer (numDofs);
+        ScalarField *pg		   = writer.allocateManagedBuffer (numDofs);
+        ScalarField *pl		   = writer.allocateManagedBuffer (numDofs);
+        ScalarField *pc 	   = writer.allocateManagedBuffer (numDofs);
+        ScalarField *rhoL	   = writer.allocateManagedBuffer (numDofs);
+        ScalarField *rhoG	   = writer.allocateManagedBuffer (numDofs);
+        ScalarField *mobL	   = writer.allocateManagedBuffer (numDofs);
+        ScalarField *mobG 	   = writer.allocateManagedBuffer (numDofs);
+        ScalarField *phasePresence = writer.allocateManagedBuffer (numDofs);
+        ScalarField *temperature   = writer.allocateManagedBuffer (numDofs);
+        ScalarField *poro          = writer.allocateManagedBuffer (numDofs);
+        ScalarField *boxVolume     = writer.allocateManagedBuffer (numDofs);
+//      ScalarField *potential     = writer.allocateManagedBuffer (numDofs);
+        ScalarField *cellNum 	   = writer.allocateManagedBuffer (numDofs);
+        ScalarField *permFactor    = writer.allocateManagedBuffer (numDofs);
+        ScalarField *precipitateVolumeFraction[numSPhases] ;
+//      ScalarField *saturationIdx[numSPhases];
+
+	for (int i = 0; i < numSPhases; ++i)
+        {
+        	precipitateVolumeFraction[i]= writer.allocateManagedBuffer (numDofs);
+        }
+
+        ScalarField *massFraction[numPhases][numComponents];
+        for (int i = 0; i < numPhases; ++i)
+            for (int j = 0; j < numComponents; ++j)
+                massFraction[i][j] = writer.allocateManagedBuffer(numDofs);
+
+        ScalarField *molarity[numComponents];
+        for (int j = 0; j < numComponents ; ++j)
+            molarity[j] = writer.allocateManagedBuffer(numDofs);
+
+        ScalarField *Perm[dim];
+        for (int j = 0; j < dim; ++j) //Permeability only in main directions xx and yy
+            Perm[j] = writer.allocateManagedBuffer(numDofs);
+
+        *boxVolume = 0;
+
+        VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
+        VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
+        ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
+
+        if (velocityOutput.enableOutput()) // check if velocity output is demanded
+              {
+                  // initialize velocity fields
+                  for (unsigned int i = 0; i < numDofs; ++i)
+                  {
+                      (*velocityN)[i] = Scalar(0);
+                      (*velocityW)[i] = Scalar(0);
+                      (*cellNum)[i] = Scalar(0.0);
+                  }
+              }
+
+        unsigned numElements = this->gridView_().size(0);
+        ScalarField *rank =
+                writer.allocateManagedBuffer (numElements);
+
+        FVElementGeometry fvGeometry;
+        VolumeVariables volVars;
+        ElementVolumeVariables elemVolVars;
+
+        ElementIterator elemIt = this->gridView_().template begin<0>();
+        ElementIterator elemEndIt = this->gridView_().template end<0>();
+        for (; elemIt != elemEndIt; ++elemIt)
+        {
+            int idx = this->problem_().elementMapper().map(*elemIt);
+            (*rank)[idx] = this->gridView_().comm().rank();
+            fvGeometry.update(this->gridView_(), *elemIt);
+            
+            elemVolVars.update(this->problem_(),
+                    *elemIt,
+                    fvGeometry,
+                    false /* oldSol? */);
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
+            int numVerts = elemIt->subEntities(dim);
+#else
+            int numVerts = elemIt->template count<dim> ();
+#endif   
+            for (int i = 0; i < numVerts; ++i)
+            {
+                int globalIdx = this->vertexMapper().map(*elemIt, i, dim);
+                volVars.update(sol[globalIdx],
+                               this->problem_(),
+                               *elemIt,
+                               fvGeometry,
+                               i,
+                               false);
+
+                (*Sg)[globalIdx]     		 = volVars.saturation(nPhaseIdx);
+                (*Sl)[globalIdx]     		 = volVars.saturation(wPhaseIdx);
+                (*pg)[globalIdx]  	 		 = volVars.pressure(nPhaseIdx);
+                (*pl)[globalIdx]   	 		 = volVars.pressure(wPhaseIdx);
+                (*pc)[globalIdx]	 		 = volVars.capillaryPressure();
+                (*rhoL)[globalIdx]	 	 	 = volVars.fluidState().density(wPhaseIdx);
+                (*rhoG)[globalIdx]	  		 = volVars.fluidState().density(nPhaseIdx);
+                (*mobL)[globalIdx] 	 	     = volVars.mobility(wPhaseIdx);
+                (*mobG)[globalIdx] 	 	     = volVars.mobility(nPhaseIdx);
+                (*boxVolume)[globalIdx]		+= fvGeometry.subContVol[i].volume;
+                (*poro)[globalIdx] 		     = volVars.porosity();
+
+                for (int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
+                {
+                    (*precipitateVolumeFraction[sPhaseIdx])[globalIdx]   = volVars.precipitateVolumeFraction(sPhaseIdx + numPhases);
+                }
+                (*temperature)[globalIdx] 	 = volVars.temperature();
+                (*permFactor)[globalIdx] 	 = volVars.permFactor();
+                (*phasePresence)[globalIdx]  = this->staticDat_[globalIdx].phasePresence;
+//                 (*potential)[globalIdx]      = (volVars.pressure(wPhaseIdx)-1e5)
+//                         /volVars.fluidState().density(wPhaseIdx)
+//                         /9.81
+//                         - (zmax - globalPos[dim-1]);
+// 
+// 
+//                 for (int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
+//                 {
+//                     (*saturationIdx[sPhaseIdx])[globalIdx]   = volVars.saturationIdx(sPhaseIdx + numPhases);
+//                 }
+
+
+                for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+                    for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                    {
+                        (*massFraction[phaseIdx][compIdx])[globalIdx]= volVars.fluidState().massFraction(phaseIdx,compIdx);
+
+                        Valgrind::CheckDefined((*massFraction[phaseIdx][compIdx])[globalIdx]);
+
+                    }
+                for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+                    (*molarity[compIdx])[globalIdx] = (volVars.fluidState().molarity(wPhaseIdx, compIdx));
+
+                Tensor K = this->perm_(this->problem_().spatialParams().intrinsicPermeability(*elemIt, fvGeometry, i));
+
+                for (int j = 0; j<dim; ++j)
+                    (*Perm[j])[globalIdx] = K[j][j] * volVars.permFactor();
+            };
+
+            // velocity output
+            if(velocityOutput.enableOutput()){
+                velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *elemIt, wPhaseIdx);
+                velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *elemIt, nPhaseIdx);
+            }
+        }  // loop over element
+
+        writer.attachVertexData(*Sg, "Sg");
+        writer.attachVertexData(*Sl, "Sl");
+        writer.attachVertexData(*pg, "pg");
+        writer.attachVertexData(*pl, "pl");
+        writer.attachVertexData(*pc, "pc");
+        writer.attachVertexData(*rhoL, "rhoL");
+        writer.attachVertexData(*rhoG, "rhoG");
+        writer.attachVertexData(*mobL, "mobL");
+        writer.attachVertexData(*mobG, "mobG");
+        writer.attachVertexData(*poro, "porosity");
+        writer.attachVertexData(*permFactor, "permFactor");
+        writer.attachVertexData(*temperature, "temperature");
+        writer.attachVertexData(*phasePresence, "phase presence");
+//         writer.attachVertexData(*potential, "potential");
+        writer.attachVertexData(*boxVolume, "boxVolume");
+        
+    
+        for (int i = 0; i < numSPhases; ++i)
+        {
+            std::ostringstream oss;
+            oss << "precipitateVolumeFraction_"
+                << FluidSystem::phaseName(numPhases + i);
+            writer.attachDofData(*precipitateVolumeFraction[i], oss.str().c_str(), isBox);
+        }
+
+//         for (int i = 0; i < numSPhases; ++i)
+//         {
+//             std::ostringstream oss;
+//             oss << "saturationIndex_"
+//                 << FluidSystem::phaseName(numPhases + i);
+//             writer.attachDofData(*saturationIdx[i], oss.str().c_str(), isBox);
+//         }
+
+        writer.attachVertexData(*Perm[0], "Kxx");
+        if (dim >= 2)
+            writer.attachVertexData(*Perm[1], "Kyy");
+        if (dim == 3)
+            writer.attachVertexData(*Perm[2], "Kzz");
+
+        for (int i = 0; i < numPhases; ++i)
+        {
+            for (int j = 0; j < numComponents; ++j)
+            {
+                std::ostringstream oss;
+                oss << "X^"
+                << FluidSystem::phaseName(i)
+                << "_"
+                << FluidSystem::componentName(j);
+                writer.attachVertexData(*massFraction[i][j], oss.str().c_str());
+            }
+        }
+
+        for (int j = 0; j < numComponents; ++j)
+        {
+            std::ostringstream oss;
+            oss << "m^w_"
+                << FluidSystem::componentName(j);
+            writer.attachVertexData(*molarity[j], oss.str().c_str());
+        }
+
+        if (velocityOutput.enableOutput()) // check if velocity output is demanded
+        {
+            writer.attachDofData(*velocityW,  "velocityW", isBox, dim);
+            writer.attachDofData(*velocityN,  "velocityN", isBox, dim);
+        }
+
+        writer.attachCellData(*rank, "process rank");
+    }
+
+    /*!
+     * \brief Update the static data of all vertices in the grid.
+     *
+     * \param curGlobalSol The current global solution
+     * \param oldGlobalSol The previous global solution
+     */
+    void updateStaticData(SolutionVector &curGlobalSol,
+                          const SolutionVector &oldGlobalSol)
+    {
+        bool wasSwitched = false;
+
+        for (unsigned i = 0; i < this->staticDat_.size(); ++i)
+            this->staticDat_[i].visited = false;
+
+        FVElementGeometry fvGeometry;
+        static VolumeVariables volVars;
+        ElementIterator it = this->gridView_().template begin<0> ();
+        const ElementIterator &endit = this->gridView_().template end<0> ();
+        for (; it != endit; ++it)
+        {
+            fvGeometry.update(this->gridView_(), *it);
+            for (int i = 0; i < fvGeometry.numScv; ++i)
+            {
+                int globalIdx = this->vertexMapper().map(*it, i, dim);
+
+                if (this->staticDat_[globalIdx].visited)
+                    continue;
+
+                this->staticDat_[globalIdx].visited = true;
+                volVars.update(curGlobalSol[globalIdx],
+                               this->problem_(),
+                               *it,
+                               fvGeometry,
+                               i,
+                               false);
+                const GlobalPosition &global = it->geometry().corner(i);
+                if (primaryVarSwitch_(curGlobalSol,
+                                      volVars,
+                                      globalIdx,
+                                      global))
+                { wasSwitched = true;
+                std::cout<<"Switch works :) "<<std::endl;
+                }
+            }
+        }
+
+        // make sure that if there was a variable switch in an
+        // other partition we will also set the switch flag
+        // for our partition.
+        if (this->gridView_().comm().size() > 1)
+        wasSwitched = this->gridView_().comm().max(wasSwitched);
+
+        setSwitched_(wasSwitched);
+    } 
+protected:
+	 
+    /*!
+     * \brief Set whether there was a primary variable switch after in
+     *        the last timestep.
+     */
+    void setSwitched_(bool yesno)
+    {
+        switchFlag_ = yesno;
+    }
+    
+    /*!
+     * \copydoc 2pnc::primaryVarSwitch_
+     */
+    bool primaryVarSwitch_(SolutionVector &globalSol,
+                           const VolumeVariables &volVars, int globalIdx,
+                           const GlobalPosition &globalPos)
+    {
+
+            // evaluate primary variable switch
+            bool wouldSwitch = false;
+            int phasePresence = this->staticDat_[globalIdx].phasePresence;
+            int newPhasePresence = phasePresence;
+
+            //check if a primary variable switch is necessary
+            if (phasePresence == bothPhases)
+            {
+                Scalar Smin = 0.0; //saturation threshold
+                if (this->staticDat_[globalIdx].wasSwitched)
+                    Smin = -0.01;
+
+                //if saturation of liquid phase is smaller 0 switch
+                if (volVars.saturation(wPhaseIdx) <= Smin)
+                {
+                    wouldSwitch = true;
+                    //liquid phase has to disappear
+                    std::cout << "Liquid Phase disappears at vertex " << globalIdx
+                                << ", coordinated: " << globalPos << ", Sl: "
+                                << volVars.saturation(wPhaseIdx) << std::endl;
+                    newPhasePresence = nPhaseOnly;
+
+                    //switch not depending on formulation
+                    //switch "Sl" to "xgH20"
+                    globalSol[globalIdx][switchIdx]
+                            = volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx /*H2O*/);
+                    //Here unlike 2pnc model we do not switch all components to to mole fraction in gas phase
+                }
+                //if saturation of gas phase is smaller than 0 switch
+                else if (volVars.saturation(nPhaseIdx) <= Smin)
+                {
+                    wouldSwitch = true;
+                    //gas phase has to disappear
+                    std::cout << "Gas Phase disappears at vertex " << globalIdx
+                                << ", coordinated: " << globalPos << ", Sg: "
+                                << volVars.saturation(nPhaseIdx) << std::endl;
+                    newPhasePresence = wPhaseOnly;
+
+                    //switch "Sl" to "xlN2"
+                    globalSol[globalIdx][switchIdx]
+                            = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx /*N2*/);
+                }
+            }
+            else if (phasePresence == nPhaseOnly)
+            {
+            Scalar sumxl = 0;
+            //Calculate sum of mole fractions (water and air) in the hypothetical liquid phase
+            //WARNING: Here numComponents is replaced by numMajorComponents as the solutes 
+            //are only present in the liquid phase and cannot condense as the liquid (water).
+            for (int compIdx = 0; compIdx < numMajorComponents; compIdx++)
+                {
+                    sumxl += volVars.fluidState().moleFraction(wPhaseIdx, compIdx);
+                }
+                    Scalar xlmax = 1.0;
+                    if (sumxl > xlmax)
+                wouldSwitch = true;
+                    if (this->staticDat_[globalIdx].wasSwitched)
+                xlmax *=1.02;
+
+            //if the sum of the mole fractions would be larger than
+            //1, wetting phase appears
+                    if (sumxl/*sum of mole fractions*/ > xlmax/*1*/)
+                    {
+                        // liquid phase appears
+                        std::cout << "Liquid Phase appears at vertex " << globalIdx
+                                << ", coordinated: " << globalPos << ", sumxl: "
+                                << sumxl << std::endl;
+                        newPhasePresence = bothPhases;
+                        if (formulation == pgSl)
+                            globalSol[globalIdx][switchIdx] = 0.0;
+                        else if (formulation == plSg)
+                            globalSol[globalIdx][switchIdx] = 1.0;
+                    //Here unlike 2pnc model we do not switch all components to to mole fraction in gas phase
+                    }
+            }
+            else if (phasePresence == wPhaseOnly)
+            {
+                Scalar xgmax = 1;
+                Scalar sumxg = 0;
+                //Calculate sum of mole fractions in the hypothetical liquid phase
+                for (int compIdx = 0; compIdx < numMajorComponents; compIdx++)
+                {
+                    sumxg += volVars.fluidState().moleFraction(nPhaseIdx, compIdx);
+                }
+                if (sumxg > xgmax)
+                    wouldSwitch = true;
+                if (this->staticDat_[globalIdx].wasSwitched)
+                    xgmax *=1.02;
+                //liquid phase appears if sum is larger than one
+                if (sumxg > xgmax)
+                {
+                    std::cout << "Gas Phase appears at vertex " << globalIdx
+                            << ", coordinated: " << globalPos << ", sumxg: "
+                            << sumxg << std::endl;
+                    newPhasePresence = bothPhases;
+                    //saturation of the liquid phase set to 0.9999 (if formulation pgSl and vice versa)
+                    if (formulation == pgSl)
+                        globalSol[globalIdx][switchIdx] = 0.999;
+                    else if (formulation == plSg)
+                        globalSol[globalIdx][switchIdx] = 0.001;
+
+                }
+            }
+            this->staticDat_[globalIdx].phasePresence = newPhasePresence;
+            this->staticDat_[globalIdx].wasSwitched = wouldSwitch;
+            return phasePresence != newPhasePresence;
+        }
+        // parameters given in constructor
+        bool switchFlag_;
+};
+
+}
+
+#include "2pncminpropertydefaults.hh"
+
+#endif

dumux/implicit/2pncmin/2pncminproperties.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   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 Properties
+ * \ingroup ImplicitProperties
+ * \ingroup TwoPNCMinModel
+ *
+ * \file
+ *
+ * \brief Defines the properties required for the two-phase n-component mineralization
+ *        fully implicit model.
+ */
+#ifndef DUMUX_2PNCMIN_PROPERTIES_HH
+#define DUMUX_2PNCMIN_PROPERTIES_HH
+
+#include <dumux/implicit/2pnc/2pncproperties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the isothermal two phase n component mineralisation problems
+NEW_TYPE_TAG(BoxTwoPNCMin, INHERITS_FROM(BoxTwoPNC));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumSPhases); //!< Number of solid phases in the system
+NEW_PROP_TAG(NumFSPhases); //!< Number of fluid and solid phases in the system
+NEW_PROP_TAG(NumSComponents); //!< Number of solid components in the system
+NEW_PROP_TAG(NumPSComponents); //!< Number of fluid and solid components in the system
+NEW_PROP_TAG(NumTraceComponents); //!< Number of trace fluid components which are not considered in the calculation of the phase density
+NEW_PROP_TAG(NumSecComponents); //!< Number of secondary components which are not primary variables
+NEW_PROP_TAG(TwoPNCMinIndices); //!< Enumerations for the 2pncMin models
+NEW_PROP_TAG(useSalinity); //!< Determines if salinity is used
+NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
+
+}
+}
+
+#endif

dumux/implicit/2pncmin/2pncminpropertydefaults.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   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 Properties
+ * \ingroup ImplicitProperties
+ * \ingroup TwoPNCMinModel
+ * \file
+ *
+ * \brief Defines default values for most properties required by the
+ *        two-phase n-component mineralization fully implicit model.
+ */
+#ifndef DUMUX_2PNCMIN_PROPERTY_DEFAULTS_HH
+#define DUMUX_2PNCMIN_PROPERTY_DEFAULTS_HH
+
+#include "2pncminindices.hh"
+#include "2pncminmodel.hh"
+#include "2pncminindices.hh"
+#include "2pncminfluxvariables.hh"
+#include "2pncminvolumevariables.hh"
+#include "2pncminproperties.hh"
+
+#include <dumux/implicit/2pnc/2pncnewtoncontroller.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
+#include <dumux/material/spatialparams/implicitspatialparams.hh>
+
+namespace Dumux
+{
+
+namespace Properties {
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+
+/*!
+ * \brief Set the property for the number of secondary components.
+ * Secondary components are components calculated from 
+ * primary components by equilibrium relations and 
+ * do not have mass balance equation on their own.
+ * These components are important in the context of bio-mineralization applications.
+ * We just forward the number from the fluid system
+ *
+ */
+SET_PROP(BoxTwoPNCMin, NumSecComponents)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+    static const int value = FluidSystem::numSecComponents;
+
+};
+/*!
+ * \brief Set the property for the number of solid phases, excluding the non-reactive matrix.
+ *
+ * We just forward the number from the fluid system
+ *
+ */
+SET_PROP(BoxTwoPNCMin, NumSPhases)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+    static const int value = FluidSystem::numSPhases;
+};
+
+/*!
+ * \brief Set the property for the number of equations. 
+ * For each component and each precipitated mineral/solid phase one equation has to
+ * be solved.
+ */
+SET_PROP(BoxTwoPNCMin, NumEq)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+    static const int value = FluidSystem::numComponents + FluidSystem::numSPhases;
+};
+
+//! Use the 2pncmin local residual operator
+SET_TYPE_PROP(BoxTwoPNCMin,
+              LocalResidual,
+              TwoPNCMinLocalResidual<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(BoxTwoPNCMin, Model, TwoPNCMinModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxTwoPNCMin, VolumeVariables, TwoPNCMinVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxTwoPNCMin, FluxVariables, TwoPNCMinFluxVariables<TypeTag>);
+
+//! The indices required by the isothermal 2pNcMin model
+SET_TYPE_PROP(BoxTwoPNCMin, Indices, TwoPNCMinIndices <TypeTag, /*PVOffset=*/0>);
+
+//! disable useSalinity for the calculation of osmotic pressure by default
+SET_BOOL_PROP(BoxTwoPNCMin, useSalinity, false);
+
+
+//! default value for the forchheimer coefficient
+// Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90.
+//        Actually the Forchheimer coefficient is also a function of the dimensions of the
+//        porous medium. Taking it as a constant is only a first approximation
+//        (Nield, Bejan, Convection in porous media, 2006, p. 10)
+SET_SCALAR_PROP(BoxModel, SpatialParamsForchCoeff, 0.55);
+
+}
+
+}
+
+#endif

dumux/implicit/2pncmin/2pncminvolumevariables.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   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, n-component mineralization model.
+ */
+#ifndef DUMUX_2PNCMin_VOLUME_VARIABLES_HH
+#define DUMUX_2PNCMin_VOLUME_VARIABLES_HH
+
+#include <dumux/implicit/common/implicitmodel.hh>
+#include <dumux/material/fluidstates/compositionalfluidstate.hh>
+#include <dumux/common/math.hh>
+#include <dune/common/collectivecommunication.hh>
+#include <vector>
+#include <iostream>
+
+#include "2pncminproperties.hh"
+#include "2pncminindices.hh"
+#include <dumux/material/constraintsolvers/computefromreferencephase2pnc.hh>
+#include <dumux/material/constraintsolvers/miscible2pnccomposition.hh>
+#include <dumux/implicit/2pnc/2pncvolumevariables.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPNCMinModel
+ * \ingroup ImplicitVolumeVariables
+ * \brief Contains the quantities which are are constant within a
+ *        finite volume in the two-phase, n-component model.
+ */
+template <class TypeTag>
+class TwoPNCMinVolumeVariables : public TwoPNCVolumeVariables<TypeTag>
+{
+    typedef TwoPNCVolumeVariables<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) Implementation;
+
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    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, MaterialLaw) MaterialLaw;
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+//     typedef typename GET_PROP_TYPE(TypeTag, Chemistry) Chemistry;
+
+    enum {
+        dim = GridView::dimension,
+        dimWorld=GridView::dimensionworld,
+
+        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+        numSPhases =  GET_PROP_VALUE(TypeTag, NumSPhases),
+        numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+        numMajorComponents = GET_PROP_VALUE(TypeTag, NumMajorComponents),
+
+        // formulations
+        formulation = GET_PROP_VALUE(TypeTag, Formulation),
+        plSg = TwoPNCFormulation::plSg,
+        pgSl = TwoPNCFormulation::pgSl,
+
+        // phase indices
+        wPhaseIdx = FluidSystem::wPhaseIdx,
+        nPhaseIdx = FluidSystem::nPhaseIdx,
+
+        // component indices
+        wCompIdx = FluidSystem::wCompIdx,
+        nCompIdx = FluidSystem::nCompIdx,
+
+        // phase presence enums
+        nPhaseOnly = Indices::nPhaseOnly,
+        wPhaseOnly = Indices::wPhaseOnly,
+        bothPhases = Indices::bothPhases,
+
+        // primary variable indices
+        pressureIdx = Indices::pressureIdx,
+        switchIdx = Indices::switchIdx,
+
+        useSalinity = GET_PROP_VALUE(TypeTag, useSalinity)
+    };
+
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+    typedef typename Grid::ctype CoordScalar;
+    typedef Dumux::Miscible2pNcComposition<Scalar, FluidSystem> Miscible2pNcComposition;
+    typedef Dumux::ComputeFromReferencePhase2pNc<Scalar, FluidSystem> ComputeFromReferencePhase2pNc;
+
+    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+    enum { dofCodim = isBox ? dim : 0 };
+public:
+  
+      //! The type of the object returned by the fluidState() method
+      typedef CompositionalFluidState<Scalar, FluidSystem> FluidState;
+    /*!
+     * \copydoc ImplicitVolumeVariables::update
+     */
+    void update(const PrimaryVariables &priVars,
+                const Problem &problem,
+                const Element &element,
+                const FVElementGeometry &fvGeometry,
+                int scvIdx,
+                bool isOldSol)
+    {
+        ParentType::update(priVars,
+                           problem,
+                           element,
+                           fvGeometry,
+                           scvIdx,
+                           isOldSol);
+
+        completeFluidState(priVars, problem, element, fvGeometry, scvIdx, this->fluidState_, isOldSol);
+
+    /////////////
+        // calculate the remaining quantities
+        /////////////
+
+    // porosity evaluation
+    initialPorosity_ = problem.spatialParams().porosity(element, fvGeometry, scvIdx);
+
+    sumPrecipitates_ = 0.0;
+    for(int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
+    {
+       precipitateVolumeFraction_[sPhaseIdx] = priVars[numComponents + sPhaseIdx];
+       sumPrecipitates_+= precipitateVolumeFraction_[sPhaseIdx];
+    }
+    
+//         for(int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
+//     {
+//         Chemistry chemistry; // the non static functions can not be called without abject 
+//         saturationIdx_[sPhaseIdx] = chemistry.omega(sPhaseIdx);
+//     }
+// TODO/FIXME: The salt crust porosity is not clearly defined. However form literature review it is
+//    found that the salt crust have porosity of approx. 10 %. Thus we restrict the decrease in porosity
+//    to this limit. Moreover in the Problem files the precipitation should also be made dependent on local
+//    porosity value, as the porous media media properties change related to salt precipitation will not be
+//    accounted otherwise.
+
+//      this->porosity_ = initialPorosity_ - sumPrecipitates_;
+
+   this->porosity_ = std::max(0.1, std::max(0.0, initialPorosity_ - sumPrecipitates_));
+
+   salinity_= 0.0;
+   moleFractionSalinity_ = 0.0;
+   for (int compIdx = numMajorComponents; compIdx< numComponents; compIdx++)	//sum of the mass fraction of the components
+   {
+       if(this->fluidState_.moleFraction(wPhaseIdx, compIdx)> 0)
+       {
+          salinity_+= this->fluidState_.massFraction(wPhaseIdx, compIdx);
+          moleFractionSalinity_ += this->fluidState_.moleFraction(wPhaseIdx, compIdx);
+       }
+    }
+
+// TODO/FIXME: Different relations for the porosoty-permeability changes are given here. We have to fins a way
+//    so that one can select the relation form the input file.
+
+    // kozeny-Carman relation
+    permeabilityFactor_  =  std::pow(((1-initialPorosity_)/(1-this->porosity_)),2)
+            * std::pow((this->porosity_/initialPorosity_),3);
+
+    // Verma-Pruess relation
+//  permeabilityFactor_  =  100 * std::pow(((this->porosity_/initialPorosity_)-0.9),2);
+
+    // Modified Fair-Hatch relation with final porosity set to 0.2 and E1=1
+//  permeabilityFactor_  =  std::pow((this->porosity_/initialPorosity_),3)
+//         * std::pow((std::pow((1 - initialPorosity_),2/3))+(std::pow((0.2 - initialPorosity_),2/3)),2)
+//         / std::pow((std::pow((1 -this->porosity_),2/3))+(std::pow((0.2 -this->porosity_),2/3)),2);
+
+    //Timur relation with residual water saturation set to 0.001
+//    permeabilityFactor_ =  0.136 * (std::pow(this->porosity_,4.4)) / (2000 * (std::pow(0.001,2)));
+
+    //Timur relation1 with residual water saturation set to 0.001
+//    permeabilityFactor_ =  0.136 * (std::pow(this->porosity_,4.4)) / (200000 * (std::pow(0.001,2)));
+
+
+    //Bern. relation
+   // permeabilityFactor_ = std::pow((this->porosity_/initialPorosity_),8);
+
+    //Tixier relation with residual water saturation set to 0.001
+    //permeabilityFactor_ = (std::pow((250 * (std::pow(this->porosity_,3)) / 0.001),2)) / InitialPermeability_;
+
+    //Coates relation with residual water saturation set to 0.001
+    //permeabilityFactor_ = (std::pow((100 * (std::pow(this->porosity_,2)) * (1-0.001) / 0.001,2))) / InitialPermeability_ ;
+
+
+    // energy related quantities not contained in the fluid state
+    //asImp_().updateEnergy_(priVars, problem,element, fvGeometry, scvIdx, isOldSol);
+    }
+
+   /*!
+    * \copydoc ImplicitModel::completeFluidState
+    * \param isOldSol Specifies whether this is the previous solution or the current one
+    */
+  static void completeFluidState(const PrimaryVariables& priVars,
+                                 const Problem& problem,
+                                 const Element& element,
+                                 const FVElementGeometry& fvGeometry,
+                                 int scvIdx,
+                                 FluidState& fluidState,
+                                 bool isOldSol = false)
+
+    {
+        Scalar t = Implementation::temperature_(priVars, problem, element,fvGeometry, scvIdx);
+        fluidState.setTemperature(t);
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
+        int dofIdxGlobal = problem.model().dofMapper().subIndex(element, scvIdx, dofCodim);
+#else
+        int dofIdxGlobal = problem.model().dofMapper().map(element, scvIdx, dofCodim);
+#endif   
+        int phasePresence = problem.model().phasePresence(dofIdxGlobal, isOldSol);
+
+        /////////////
+        // set the saturations
+        /////////////
+
+        Scalar Sg;
+        if (phasePresence == nPhaseOnly)
+            Sg = 1.0;
+        else if (phasePresence == wPhaseOnly) {
+            Sg = 0.0;
+        }
+        else if (phasePresence == bothPhases) {
+            if (formulation == plSg)
+                Sg = priVars[switchIdx];
+            else if (formulation == pgSl)
+                Sg = 1.0 - priVars[switchIdx];
+            else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
+        }
+        else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
+        fluidState.setSaturation(nPhaseIdx, Sg);
+        fluidState.setSaturation(wPhaseIdx, 1.0 - Sg);
+
+        /////////////
+        // set the pressures of the fluid phases
+        /////////////
+
+        // calculate capillary pressure
+        const MaterialLawParams &materialParams = problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
+        Scalar pc = MaterialLaw::pc(materialParams, 1 - Sg);
+
+        // extract the pressures
+        if (formulation == plSg) {
+            fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]);
+            fluidState.setPressure(nPhaseIdx, priVars[pressureIdx] + pc);
+        }
+        else if (formulation == pgSl) {
+            fluidState.setPressure(nPhaseIdx, priVars[pressureIdx]);
+            fluidState.setPressure(wPhaseIdx, priVars[pressureIdx] - pc);
+        }
+        else DUNE_THROW(Dune::InvalidStateException, "Formulation: " << formulation << " is invalid.");
+
+        /////////////
+        // calculate the phase compositions
+        /////////////
+
+    typename FluidSystem::ParameterCache paramCache;
+
+        // now comes the tricky part: calculate phase composition
+        if (phasePresence == bothPhases) {
+            // both phases are present, phase composition results from
+            // the gas <-> liquid equilibrium. This is
+            // the job of the "MiscibleMultiPhaseComposition"
+            // constraint solver
+
+            // set the known mole fractions in the fluidState so that they
+            // can be used by the Miscible2pNcComposition constraint solver
+            for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
+            {
+            	fluidState.setMoleFraction(wPhaseIdx, compIdx, priVars[compIdx]);
+            }
+
+            Miscible2pNcComposition::solve(fluidState,
+                                            paramCache,
+                                            wPhaseIdx,	//known phaseIdx
+                                            /*setViscosity=*/true,
+                                            /*setInternalEnergy=*/false);
+        }
+        else if (phasePresence == nPhaseOnly){
+
+            Dune::FieldVector<Scalar, numComponents> moleFrac;
+            Dune::FieldVector<Scalar, numComponents> fugCoeffL;
+            Dune::FieldVector<Scalar, numComponents> fugCoeffG;
+
+            for (int compIdx=0; compIdx<numComponents; ++compIdx)
+            {
+                fugCoeffL[compIdx] = FluidSystem::fugacityCoefficient(fluidState,
+                                        paramCache,
+                                        wPhaseIdx,
+                                        compIdx);
+                fugCoeffG[compIdx] = FluidSystem::fugacityCoefficient(fluidState,
+                                        paramCache,
+                                        nPhaseIdx,
+                                        compIdx);
+            }
+            for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
+                moleFrac[compIdx] = (priVars[compIdx]*fugCoeffL[compIdx]*fluidState.pressure(wPhaseIdx))
+                /(fugCoeffG[compIdx]*fluidState.pressure(nPhaseIdx));
+            
+            moleFrac[wCompIdx] =  priVars[switchIdx];
+            Scalar sumMoleFracNotGas = 0;
+            for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
+            {
+            		sumMoleFracNotGas+=moleFrac[compIdx];
+            }
+            sumMoleFracNotGas += moleFrac[wCompIdx];
+            moleFrac[nCompIdx] = 1 - sumMoleFracNotGas;
+
+//          typedef Dune::FieldMatrix<Scalar, numComponents, numComponents> Matrix;
+//          typedef Dune::FieldVector<Scalar, numComponents> Vector;
+
+
+            // Set fluid state mole fractions
+            for (int compIdx=0; compIdx<numComponents; ++compIdx)
+            {
+            	fluidState.setMoleFraction(nPhaseIdx, compIdx, moleFrac[compIdx]);
+            }
+
+            // calculate the composition of the remaining phases (as
+            // well as the densities of all phases). this is the job
+            // of the "ComputeFromReferencePhase2pNc" constraint solver
+            ComputeFromReferencePhase2pNc::solve(fluidState,
+                                             paramCache,
+                                             nPhaseIdx,
+                                             /*setViscosity=*/true,
+                                             /*setInternalEnergy=*/false);
+
+            }
+        else if (phasePresence == wPhaseOnly){
+        
+        // only the liquid phase is present, i.e. liquid phase
+        // composition is stored explicitly.
+        // extract _mass_ fractions in the gas phase
+            Dune::FieldVector<Scalar, numComponents> moleFrac;
+
+            for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
+            {
+            	moleFrac[compIdx] = priVars[compIdx];
+            }
+            moleFrac[nCompIdx] = priVars[switchIdx];
+            Scalar sumMoleFracNotWater = 0;
+            for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
+            {
+            		sumMoleFracNotWater+=moleFrac[compIdx];
+            }
+            sumMoleFracNotWater += moleFrac[nCompIdx];
+            moleFrac[wCompIdx] = 1 -sumMoleFracNotWater;
+
+//             convert mass to mole fractions and set the fluid state
+            for (int compIdx=0; compIdx<numComponents; ++compIdx)
+            {
+                fluidState.setMoleFraction(wPhaseIdx, compIdx, moleFrac[compIdx]);
+            }
+
+//             calculate the composition of the remaining phases (as
+//             well as the densities of all phases). this is the job
+//             of the "ComputeFromReferencePhase2pNc" constraint solver
+            ComputeFromReferencePhase2pNc::solve(fluidState,
+                                             paramCache,
+                                             wPhaseIdx,
+                                             /*setViscosity=*/true,
+                                             /*setInternalEnergy=*/false);
+        }
+        paramCache.updateAll(fluidState);
+        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+        {
+            Scalar h = Implementation::enthalpy_(fluidState, paramCache, phaseIdx);
+            fluidState.setEnthalpy(phaseIdx, h);
+        }
+    }
+    /*!
+     * \brief Returns the volume fraction of the precipitate (solid phase) 
+     * for the given phaseIdx
+     * 
+     * \param phaseIdx the index of the solid phase
+     */
+    Scalar precipitateVolumeFraction(int phaseIdx) const
+    { return precipitateVolumeFraction_[phaseIdx - numPhases]; }
+   
+    /*!
+     * \brief Returns the inital porosity of the 
+     * pure, precipitate-free porous medium
+     */
+    Scalar InitialPorosity() const
+    { return initialPorosity_;}
+    
+    /*!
+     * \brief Returns the inital permeability of the 
+     * pure, precipitate-free porous medium
+     */    
+    Scalar InitialPermeability() const
+    { return InitialPermeability_;}
+    
+    /*!
+     * \brief Returns the factor for the reduction of the initial permeability 
+     * due precipitates in the porous medium
+     */ 	
+    Scalar permFactor() const
+    { return permeabilityFactor_; }
+    
+    /*!
+     * \brief Returns the mole fraction of a component in the phase
+     * 
+     * \param phaseIdx the index of the fluid phase
+     * \param compIdx the index of the component
+     */ 
+    Scalar moleFraction(int phaseIdx, int compIdx) const
+    {
+       return this->fluidState_.moleFraction(phaseIdx, compIdx);
+    }
+    
+    /*!
+     * \brief Returns the mole fraction of the salinity in the liquid phase
+     */
+    Scalar moleFracSalinity() const
+    {
+    	return moleFractionSalinity_;
+    }
+    
+    /*!
+     * \brief Returns the salinity (mass fraction) in the liquid phase
+     */
+    Scalar salinity() const
+    {
+        return salinity_;
+    }
+    
+    /*!
+     * \brief Returns the density of the phase for all fluid and solid phases 
+     * 
+     * \param phaseIdx the index of the fluid phase
+     */ 
+    Scalar density(int phaseIdx) const
+    {
+        if (phaseIdx < numPhases)
+            return this->fluidState_.density(phaseIdx);
+#if SALINIZATION
+        else if (phaseIdx >= numPhases)
+        return FluidSystem::precipitateDensity(phaseIdx);
+#endif
+        else
+            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+    }
+    /*!
+        * \brief Returns the liquid vapor pressure within the control volume.
+        */
+#if SALINIZATION
+    Scalar vaporPressure() const
+    { return FluidSystem::vaporPressure(this->fluidState_.temperature(/*phaseIdx=*/0),this->fluidState_.moleFraction(wPhaseIdx, FluidSystem::NaClIdx)); }
+#endif
+
+    /*!
+     * \brief Returns the mass density of a given phase within the
+     *        control volume.
+     *
+     * \param phaseIdx The phase index
+     */
+    Scalar molarDensity(int phaseIdx) const
+    {
+        if (phaseIdx < numPhases)
+            return this->fluidState_.molarDensity(phaseIdx);
+#if SALINIZATION
+        else if (phaseIdx >= numPhases)
+            return FluidSystem::precipitateMolarDensity(phaseIdx);
+#endif
+        else
+            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+    }
+    
+    /*!
+     * \brief Returns the molality of a component in the phase
+     * 
+     * \param phaseIdx the index of the fluid phase
+     * \param compIdx the index of the component
+     */ 
+     Scalar molality(int phaseIdx, int compIdx) const // [moles/Kg]
+    { return this->fluidState_.moleFraction(phaseIdx, compIdx)/FluidSystem::molarMass(compIdx);}
+
+protected:
+    friend class TwoPNCVolumeVariables<TypeTag>;
+    static Scalar temperature_(const PrimaryVariables &priVars,
+                                const Problem& problem,
+                                const Element &element,
+                                const FVElementGeometry &fvGeometry,
+                                int scvIdx)
+    {
+        return problem.temperatureAtPos(fvGeometry.subContVol[scvIdx].global);
+    }
+    
+    template<class ParameterCache>
+    static Scalar enthalpy_(const FluidState& fluidState,
+                            const ParameterCache& paramCache,
+                            int phaseIdx)
+    {
+        return 0;
+    }
+
+   /*!
+    * \brief Update all quantities for a given control volume.
+    *
+    * \param priVars The solution primary variables
+    * \param problem The problem
+    * \param element The element
+    * \param fvGeometry Evaluate function with solution of current or previous time step
+    * \param scvIdx The local index of the SCV (sub-control volume)
+    * \param isOldSol Evaluate function with solution of current or previous time step
+    */
+        void updateEnergy_(const PrimaryVariables &priVars,
+                           const Problem &problem,
+                           const Element &element,
+                           const FVElementGeometry &fvGeometry,
+                           const int scvIdx,
+                           bool isOldSol)
+        { };
+
+    Scalar precipitateVolumeFraction_[numSPhases];
+//     Scalar saturationIdx_[numSPhases];
+    Scalar permeabilityFactor_;
+    Scalar initialPorosity_;
+    Scalar InitialPermeability_;
+    Scalar sumPrecipitates_;
+    Scalar salinity_;
+    Scalar moleFractionSalinity_;
+
+private:
+    Implementation &asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    const Implementation &asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+};
+
+} // end namespace
+
+#endif