From e2e8d492ea3f26cb2c2f2fa118e59b8e80107693 Mon Sep 17 00:00:00 2001
From: vishal jambhekar <vishal.jambhekar@iws.uni-stuttgart.de>
Date: Tue, 21 Jul 2015 16:49:57 +0000
Subject: [PATCH] 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
---
.../implicit/2pncmin/2pncminfluxvariables.hh | 162 +++++
dumux/implicit/2pncmin/2pncminindices.hh | 48 ++
.../implicit/2pncmin/2pncminlocalresidual.hh | 140 +++++
dumux/implicit/2pncmin/2pncminmodel.hh | 586 ++++++++++++++++++
dumux/implicit/2pncmin/2pncminproperties.hh | 63 ++
.../2pncmin/2pncminpropertydefaults.hh | 129 ++++
.../2pncmin/2pncminvolumevariables.hh | 554 +++++++++++++++++
7 files changed, 1682 insertions(+)
create mode 100644 dumux/implicit/2pncmin/2pncminfluxvariables.hh
create mode 100644 dumux/implicit/2pncmin/2pncminindices.hh
create mode 100644 dumux/implicit/2pncmin/2pncminlocalresidual.hh
create mode 100644 dumux/implicit/2pncmin/2pncminmodel.hh
create mode 100644 dumux/implicit/2pncmin/2pncminproperties.hh
create mode 100644 dumux/implicit/2pncmin/2pncminpropertydefaults.hh
create mode 100644 dumux/implicit/2pncmin/2pncminvolumevariables.hh
diff --git a/dumux/implicit/2pncmin/2pncminfluxvariables.hh b/dumux/implicit/2pncmin/2pncminfluxvariables.hh
new file mode 100644
index 0000000000..50d5ce61c4
--- /dev/null
+++ b/dumux/implicit/2pncmin/2pncminfluxvariables.hh
@@ -0,0 +1,162 @@
+// -*- 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
diff --git a/dumux/implicit/2pncmin/2pncminindices.hh b/dumux/implicit/2pncmin/2pncminindices.hh
new file mode 100644
index 0000000000..dc1b2e9e0c
--- /dev/null
+++ b/dumux/implicit/2pncmin/2pncminindices.hh
@@ -0,0 +1,48 @@
+// -*- 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
diff --git a/dumux/implicit/2pncmin/2pncminlocalresidual.hh b/dumux/implicit/2pncmin/2pncminlocalresidual.hh
new file mode 100644
index 0000000000..c0fc7f9529
--- /dev/null
+++ b/dumux/implicit/2pncmin/2pncminlocalresidual.hh
@@ -0,0 +1,140 @@
+// -*- 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
diff --git a/dumux/implicit/2pncmin/2pncminmodel.hh b/dumux/implicit/2pncmin/2pncminmodel.hh
new file mode 100644
index 0000000000..54c5d0b3fc
--- /dev/null
+++ b/dumux/implicit/2pncmin/2pncminmodel.hh
@@ -0,0 +1,586 @@
+// -*- 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
diff --git a/dumux/implicit/2pncmin/2pncminproperties.hh b/dumux/implicit/2pncmin/2pncminproperties.hh
new file mode 100644
index 0000000000..18e6c866e6
--- /dev/null
+++ b/dumux/implicit/2pncmin/2pncminproperties.hh
@@ -0,0 +1,63 @@
+// -*- 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
diff --git a/dumux/implicit/2pncmin/2pncminpropertydefaults.hh b/dumux/implicit/2pncmin/2pncminpropertydefaults.hh
new file mode 100644
index 0000000000..516ee45bc6
--- /dev/null
+++ b/dumux/implicit/2pncmin/2pncminpropertydefaults.hh
@@ -0,0 +1,129 @@
+// -*- 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
diff --git a/dumux/implicit/2pncmin/2pncminvolumevariables.hh b/dumux/implicit/2pncmin/2pncminvolumevariables.hh
new file mode 100644
index 0000000000..e3b5883083
--- /dev/null
+++ b/dumux/implicit/2pncmin/2pncminvolumevariables.hh
@@ -0,0 +1,554 @@
+// -*- 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
--
GitLab