[1pncmin] add model to next

not yet working

dumux/porousmediumflow/1pncmin/CMakeLists.txt

+add_subdirectory("implicit")

dumux/porousmediumflow/1pncmin/implicit/CMakeLists.txt

+
+#install headers
+install(FILES
+indices.hh
+localresidual.hh
+model.hh
+properties.hh
+propertydefaults.hh
+volumevariables.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/1pncmin/implicit)

dumux/porousmediumflow/1pncmin/implicit/indices.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 one-phase n-component mineralization
+ *        fully implicit model.
+ */
+#ifndef DUMUX_1PNCMIN_INDICES_HH
+#define DUMUX_1PNCMIN_INDICES_HH
+
+#include <dumux/porousmediumflow/1pnc/implicit/indices.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup OnePNCMinModel
+ * \ingroup ImplicitIndices
+ * \brief The indices for the isothermal one-phase n-component mineralization model.
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int PVOffset = 0>
+    class OnePNCMinIndices: public OnePNCIndices<TypeTag, PVOffset>
+{
+};
+
+// \}
+
+}
+
+#endif

dumux/porousmediumflow/1pncmin/implicit/localresidual.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 one-phase n-component mineralisation box model.
+ */
+
+#ifndef DUMUX_1PNCMIN_LOCAL_RESIDUAL_BASE_HH
+#define DUMUX_1PNCMIN_LOCAL_RESIDUAL_BASE_HH
+
+#include "properties.hh"
+#include <dumux/porousmediumflow/compositional/localresidual.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup OnePNCMinModel
+ * \ingroup ImplicitLocalResidual
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ *        using the one-phase n-component mineralization fully implicit model.
+ *
+ * This class is used to fill the gaps in ImplicitLocalResidual for the one-phase n-component flow.
+ */
+template<class TypeTag>
+class OnePNCMinLocalResidual: public OnePNCLocalResidual<TypeTag>
+{
+protected:
+//     typedef OnePNCLocalResidual<TypeTag> ParentType;
+    using ParentType = CompositionalLocalResidual<TypeTag>;
+    using ThisType = TwoPNCMinLocalResidual<TypeTag>;
+    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables)
+//     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, BoundaryTypes) BoundaryTypes;
+//     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),
+
+        pressureIdx = Indices::pressureIdx,
+        firstMoleFracIdx = Indices::firstMoleFracIdx,
+
+        phaseIdx = Indices::phaseIdx,
+        sPhaseIdx = 1, // don't use the sPhaseIdx of the fluidsystem
+
+        conti0EqIdx = Indices::conti0EqIdx,
+    };
+
+//     typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+//     typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
+//     typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
+//     typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+public:
+    /*!
+     * \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 1pnc model, here, the storage of solid phases is included too.
+     *
+     *  \param scv the SCV (sub-control-volume)
+     *  \param volVars The volume variables of the right time step
+     */
+    PrimaryVariables computeStorage(const SubControlVolume& scv,
+                                    const VolumeVariables& volVars) const
+    {
+        //call parenttype function
+        auto storage = ParentType::computeStorage(scv, volVars);
+
+        // Compute storage term of all solid (precipitated) phases (excluding the non-reactive matrix)
+         for (int Idx = sPhaseIdx; Idx < numPhases + numSPhases; ++Idx)
+        {
+            auto eqIdx = conti0EqIdx + numComponents-numPhases + Idx;
+            storage[eqIdx] += volVars.precipitateVolumeFraction(Idx)*volVars.molarDensity(Idx);
+        }
+
+        return storage;
+    }
+};
+
+} // end namespace Dumux
+
+#endif

dumux/porousmediumflow/1pncmin/implicit/properties.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 OnePNCMinModel
+ *
+ * \file
+ *
+ * \brief Defines the properties required for the one-phase n-component mineralization
+ *        fully implicit model.
+ */
+#ifndef DUMUX_1PNCMIN_PROPERTIES_HH
+#define DUMUX_1PNCMIN_PROPERTIES_HH
+
+#include <dumux/porousmediumflow/1pnc/implicit/properties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the isothermal two phase n component mineralisation problems
+NEW_TYPE_TAG(OnePNCMin, INHERITS_FROM(OnePNC));
+NEW_TYPE_TAG(BoxOnePNCMin, INHERITS_FROM(BoxModel, OnePNCMin));
+NEW_TYPE_TAG(CCOnePNCMin, INHERITS_FROM(CCModel, OnePNCMin));
+
+//! The type tags for the corresponding non-isothermal problems
+NEW_TYPE_TAG(OnePNCMinNI, INHERITS_FROM(OnePNCMin, NonIsothermal));
+NEW_TYPE_TAG(BoxOnePNCMinNI, INHERITS_FROM(BoxModel, OnePNCMinNI));
+NEW_TYPE_TAG(CCOnePNCMinNI, INHERITS_FROM(CCModel, OnePNCMinNI));
+
+//////////////////////////////////////////////////////////////////
+// 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(OnePNCMinIndices); //!< Enumerations for the 2pncMin models
+NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
+
+}
+}
+
+#endif

dumux/porousmediumflow/1pncmin/implicit/propertydefaults.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 OnePNCMinModel
+ * \file
+ *
+ * \brief Defines default values for most properties required by the
+ *        two-phase n-component mineralization fully implicit model.
+ */
+#ifndef DUMUX_1PNCMIN_PROPERTY_DEFAULTS_HH
+#define DUMUX_1PNCMIN_PROPERTY_DEFAULTS_HH
+
+#include "model.hh"
+#include "indices.hh"
+#include "volumevariables.hh"
+#include "properties.hh"
+
+#include <dumux/porousmediumflow/1pnc/implicit/newtoncontroller.hh>
+#include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
+#include <dumux/material/spatialparams/implicit.hh>
+#include <dumux/material/fluidmatrixinteractions/1p/thermalconductivityaverage.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(OnePNCMin, NumSComponents)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+public:
+    static const int value = FluidSystem::numSComponents;
+};
+/*!
+ * \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(OnePNCMin, 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(OnePNCMin, NumEq)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+//     static const int value = FluidSystem::numComponents + FluidSystem::numSPhases;
+    static const int value = FluidSystem::numComponents + FluidSystem::numSComponents; //steamaircao2h2 has 2 components in the fluidphase
+};
+
+/*!
+ * \brief The fluid state which is used by the volume variables to
+ *        store the thermodynamic state. This should be chosen
+ *        appropriately for the model ((non-)isothermal, equilibrium, ...).
+ *        This can be done in the problem.
+ */
+SET_PROP(OnePNCMin, FluidState){
+    private:
+        typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+        typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    public:
+        typedef CompositionalFluidState<Scalar, FluidSystem> type;
+};
+
+//! Use the 2pncmin local residual operator
+SET_TYPE_PROP(OnePNCMin,
+              LocalResidual,
+              OnePNCMinLocalResidual<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(OnePNCMin, Model, OnePNCMinModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(OnePNCMin, VolumeVariables, OnePNCMinVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+// SET_TYPE_PROP(OnePNCMin, FluxVariables, OnePNCMinFluxVariables<TypeTag>);
+
+//! The indices required by the isothermal 2pNcMin model
+SET_TYPE_PROP(OnePNCMin, Indices, OnePNCMinIndices <TypeTag, /*PVOffset=*/0>);
+
+//! 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(OnePNCMin, SpatialParamsForchCoeff, 0.55);
+
+//set isothermal VolumeVariables
+SET_TYPE_PROP(OnePNCMin, IsothermalVolumeVariables, OnePNCMinVolumeVariables<TypeTag>);
+
+//! Somerton is used as default model to compute the effective thermal heat conductivity
+SET_PROP(OnePNCMinNI, ThermalConductivityModel)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+public:
+    typedef ThermalConductivityAverage<Scalar> type;
+};
+
+
+SET_BOOL_PROP(OnePNCMinNI, NiOutputLevel, 0);
+
+//////////////////////////////////////////////////////////////////
+// Property values for isothermal model required for the general non-isothermal model
+//////////////////////////////////////////////////////////////////
+
+// set isothermal Model
+SET_TYPE_PROP(OnePNCMinNI, IsothermalModel, OnePNCMinModel<TypeTag>);
+
+// set isothermal FluxVariables
+//SET_TYPE_PROP(OnePNCMinNI, IsothermalFluxVariables, OnePNCMinFluxVariables<TypeTag>);
+
+//set isothermal VolumeVariables
+SET_TYPE_PROP(OnePNCMinNI, IsothermalVolumeVariables, OnePNCMinVolumeVariables<TypeTag>);
+
+//set isothermal LocalResidual
+SET_TYPE_PROP(OnePNCMinNI, IsothermalLocalResidual, OnePNCMinLocalResidual<TypeTag>);
+
+//set isothermal Indices
+SET_TYPE_PROP(OnePNCMinNI, IsothermalIndices, OnePNCMinIndices<TypeTag, /*PVOffset=*/0>);
+
+//set isothermal NumEq
+SET_PROP(OnePNCMinNI, IsothermalNumEq)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+    static const int value = FluidSystem::numComponents +FluidSystem::numSComponents;// in NonIsothermal 1 is substracted by default
+};
+}
+}
+
+#endif

dumux/porousmediumflow/1pncmin/implicit/volumevariables.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_1PNCMIN_VOLUME_VARIABLES_HH
+#define DUMUX_1PNCMIN_VOLUME_VARIABLES_HH
+
+
+#include <dumux/common/math.hh>
+#include <dumux/implicit/model.hh>
+#include <dumux/material/fluidstates/compositional.hh>
+#include <dumux/porousmediumflow/1pnc/implicit/volumevariables.hh>
+
+#include "properties.hh"
+#include "indices.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup OnePNCMinModel
+ * \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 OnePNCMinVolumeVariables : public OnePNCVolumeVariables<TypeTag>
+{
+    // base type is used for energy related quantites
+    using BaseType = ImplicitVolumeVariables<TypeTag>;
+
+    using ParentType = OnePNCVolumeVariables<TypeTag>;
+    using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+//     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+
+    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),
+        numSComponents = GET_PROP_VALUE(TypeTag, NumSComponents), //if there is more than 1 component in the solid phase
+
+        phaseCompIdx = Indices::phaseCompIdx,
+
+        // phase indices
+        phaseIdx = FluidSystem::gPhaseIdx,
+        cPhaseIdx = Indices::phaseIdx +1,
+        hPhaseIdx = Indices::phaseIdx +2,
+
+        // primary variable indices
+        pressureIdx = Indices::pressureIdx,
+        firstMoleFracIdx = Indices::firstMoleFracIdx,
+
+    };
+
+    using Element = typename GridView::template Codim<0>::Entity;
+    using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
+    using CoordScalar = typename Grid::ctype;
+
+    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+    enum { dofCodim = isBox ? dim : 0 };
+
+public:
+
+    using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+    /*!
+     * \copydoc ImplicitVolumeVariables::update
+     */
+    void update(const ElementSolutionVector &elemSol,
+                const Problem &problem,
+                const Element &element,
+                const SubControlVolume& scv)
+    {
+        ParentType::update(elemSol, problem, element, scv);
+
+//         ParentType::update(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
+//         completeFluidState(priVars, problem, element, fvGeometry, scvIdx, this->fluidState_, isOldSol);
+
+        /////////////
+        // calculate the remaining quantities
+        /////////////
+
+        auto&& priVars = isBox ? elemSol[scv.index()] : elemSol[0];
+
+        // porosity evaluation
+        initialPorosity_ = problem.spatialParams().porosity(element, scv);
+        minimumPorosity_ = problem.spatialParams().porosityMin(element, scv);
+        maximumPorosity_ = problem.spatialParams().porosityMax(element, scv);
+
+        sumPrecipitates_ = 0.0;
+        for(int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
+        {
+           precipitateVolumeFraction_[sPhaseIdx] = priVars[numComponents + sPhaseIdx];
+           sumPrecipitates_+= precipitateVolumeFraction_[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_ = 1 - sumPrecipitates_;
+
+        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_(elemSol, problem, element, scv);
+    }
+
+   /*!
+    * \copydoc ImplicitModel::completeFluidState
+    * \param isOldSol Specifies whether this is the previous solution or the current one
+    */
+    static void completeFluidState(const ElementSolutionVector& elemSol,
+                                   const Problem& problem,
+                                   const Element& element,
+                                   const SubControlVolume& scv,
+                                   FluidState& fluidState)
+
+    {
+//         Scalar t = ParentType::temperature(elemSol, problem, element, scv); //old
+        Scalar t = BaseType::temperature(elemSol, problem, element, scv);
+        fluidState.setTemperature(t);
+
+        /////////////
+        // set the saturations
+        /////////////
+
+        fluidState.setSaturation(phaseIdx, 1.0 );
+
+        /////////////
+        // set the pressures of the fluid phase
+        /////////////
+        const auto& priVars = ParentType::extractDofPriVars(elemSol, scv);
+        fluidState.setPressure(phaseIdx, priVars[pressureIdx]);
+
+        /////////////
+        // calculate the phase compositions
+        /////////////
+
+        typename FluidSystem::ParameterCache paramCache;
+
+        Dune::FieldVector<Scalar, numComponents> moleFrac;
+
+        Scalar sumMoleFracNotWater = 0;
+        for (int compIdx=firstMoleFracIdx; compIdx<numComponents; ++compIdx)
+        {
+            moleFrac[compIdx] = priVars[compIdx];
+            sumMoleFracNotWater+=moleFrac[compIdx];
+        }
+
+        moleFrac[0] = 1 -sumMoleFracNotWater;
+
+//         //mole fractions for the solid phase
+//         moleFrac[firstMoleFracIdx+1]= priVars[firstMoleFracIdx+1];
+//         moleFrac[firstMoleFracIdx +2] = 1- moleFrac[firstMoleFracIdx+1];
+
+        // convert mass to mole fractions and set the fluid state
+        for (int compIdx=0; compIdx<numComponents; ++compIdx)
+        {
+            fluidState.setMoleFraction(phaseIdx, compIdx, moleFrac[compIdx]);
+        }
+
+//         // set mole fractions for the solid phase
+//         fluidState.setMoleFraction(cPhaseIdx, firstMoleFracIdx+1, moleFrac[firstMoleFracIdx+1]);
+//         fluidState.setMoleFraction(hPhaseIdx, firstMoleFracIdx+2, moleFrac[firstMoleFracIdx+2]);
+
+        paramCache.updateAll(fluidState);
+
+//         Scalar h = Implementation::enthalpy_(fluidState, paramCache, phaseIdx);
+        Scalar h = BaseType::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 permeabilityFactor() const
+    { return permeabilityFactor_; }
+
+    /*!
+     * \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);
+        else if (phaseIdx > numPhases)
+            return FluidSystem::precipitateDensity(phaseIdx);
+        else
+            DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+    }
+    /*!
+     * \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 <  1)
+            return this->fluidState_.molarDensity(phaseIdx);
+        else if (phaseIdx >= 1){
+            /*Attention: sPhaseIdx of the fluidsystem and the model can be different.*/
+//             std::cout << "FluidSystem::precipitateMolarDensity("<<phaseIdx<<") = " << FluidSystem::precipitateMolarDensity(phaseIdx) << "\n";
+//             for (int phaseIdx=1; phaseIdx<numSPhases; ++phaseIdx)
+            return FluidSystem::precipitateMolarDensity(phaseIdx);
+        }
+        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
+     * \f$\mathrm{molality}=\frac{n_\mathrm{component}}{m_\mathrm{solvent}}
+     * =\frac{n_\mathrm{component}}{n_\mathrm{solvent}*M_\mathrm{solvent}}\f$
+     * compIdx of the main component (solvent) in the
+     * phase is equal to the phaseIdx
+     */
+    Scalar molality(int phaseIdx, int compIdx) const // [moles/Kg]
+    { return this->fluidState_.moleFraction(phaseIdx, compIdx)
+                  /(this->fluidState_.moleFraction(phaseIdx, phaseIdx)
+                  * FluidSystem::molarMass(phaseIdx));}
+
+   /*!
+    * Circumvents the inheritance architecture of the nonisothermal model
+    */
+    static Scalar callProtectedTemperature(const ElementSolutionVector& elemSol,
+                                        const Problem& problem,
+                                        const Element& element,
+                                        const SubControlVolume& scv)
+    {
+//          return Implementation::temperature_(priVars, problem,element, fvGeometry, scvIdx);
+        return BaseType::temperature(elemSol, problem, element, scv);
+    }
+
+    /*!
+    * Circumvents the inheritance architecture of the ninisothermal model
+    */
+    void callProtectedUpdateEnergy(const ElementSolutionVector& elemSol,
+                                        const Problem& problem,
+                                        const Element& element,
+                                        const SubControlVolume& scv)
+    {
+        asImp_().updateEnergy_(elemSol, problem, element, scv);
+    };
+
+protected:
+    friend class OnePNCVolumeVariables<TypeTag>;
+    static Scalar temperature_(const ElementSolutionVector& elemSol,
+                                        const Problem& problem,
+                                        const Element& element,
+                                        const SubControlVolume& scv)
+    {
+        return problem.temperatureAtPos(scv);
+    }
+
+    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 ElementSolutionVector& elemSol,
+                       const Problem& problem,
+                       const Element& element,
+                       const SubControlVolume& scv)
+    { };
+
+    Scalar precipitateVolumeFraction_[numSPhases];
+    Scalar permeabilityFactor_;
+    Scalar initialPorosity_;
+    Scalar initialPermeability_;
+    Scalar minimumPorosity_;
+    Scalar maximumPorosity_;
+    Scalar sumPrecipitates_;
+
+
+private:
+    Implementation &asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    const Implementation &asImp_() const
+    { return *static_cast<const Implementation*>(this); }
+};
+
+} // end namespace
+
+#endif

dumux/porousmediumflow/CMakeLists.txt

 add_subdirectory("1p")
 add_subdirectory("1p2c")
 add_subdirectory("1pnc")
+add_subdirectory("1pncmin")
 add_subdirectory("2p")
 add_subdirectory("2p1c")
 add_subdirectory("2p2c")