diff --git a/dumux/implicit/2pnc/2pncfluxvariables.hh b/dumux/implicit/2pnc/2pncfluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..165dc8857d7537c259d4416e693764efe256e86e
--- /dev/null
+++ b/dumux/implicit/2pnc/2pncfluxvariables.hh
@@ -0,0 +1,415 @@
+// -*- 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 model fully implicit model.
+ */
+#ifndef DUMUX_2PNC_FLUX_VARIABLES_HH
+#define DUMUX_2PNC_FLUX_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+#include <dumux/common/spline.hh>
+
+#include "2pncproperties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPNCModel
+ * \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 fully implicit model.
+ *
+ * This means pressure and concentration gradients, phase densities at
+ * the integration point, etc.
+ */
+
+template <class TypeTag>
+class TwoPNCFluxVariables : public GET_PROP_TYPE(TypeTag, BaseFluxVariables)
+{
+ typedef typename GET_PROP_TYPE(TypeTag, BaseFluxVariables) BaseFluxVariables;
+ 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
+ */
+ TwoPNCFluxVariables(const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int fIdx,
+ const ElementVolumeVariables &elemVolVars,
+ const bool onBoundary = false)
+ : BaseFluxVariables(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary)
+ {
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ density_[phaseIdx] = Scalar(0);
+ molarDensity_[phaseIdx] = Scalar(0);
+ potentialGrad_[phaseIdx] = Scalar(0);
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ massFractionGrad_[phaseIdx][compIdx] = Scalar(0);
+ moleFractionGrad_[phaseIdx][compIdx] = Scalar(0);
+ }
+ }
+ calculateGradients_(problem, element, elemVolVars);
+ calculateVelocities_(problem, element, elemVolVars);
+ calculateporousDiffCoeff_(problem, element, elemVolVars);
+ };
+
+protected:
+ void calculateGradients_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemVolVars)
+ {
+ // calculate gradients
+ DimVector tmp(0.0);
+ for (int idx = 0;
+ idx < this->fvGeometry_.numScv;
+ idx++) // loop over adjacent vertices
+ {
+ // FE gradient at vertex idx
+ const DimVector &feGrad = face().grad[idx];
+
+ // compute sum of pressure gradients for each phase
+ for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
+ {
+ // the pressure gradient
+ tmp = feGrad;
+ tmp *= elemVolVars[idx].pressure(phaseIdx); //FE grad times phase pressure
+ potentialGrad_[phaseIdx] += tmp;
+ }
+
+ // the concentration gradient of the non-wetting
+ // component in the wetting phase
+
+ for(int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ for(int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ if(compIdx != phaseIdx) //No grad is needed for this case
+ {
+ tmp = feGrad;
+ tmp *= elemVolVars[idx].fluidState().massFraction(phaseIdx, compIdx);
+ massFractionGrad_[phaseIdx][compIdx] += tmp;
+
+ tmp = feGrad;
+ tmp *= elemVolVars[idx].fluidState().moleFraction(phaseIdx, compIdx);
+ moleFractionGrad_[phaseIdx][compIdx] += tmp;
+ }
+ }
+ }
+ }
+
+ // correct the pressure gradients by the hydrostatic
+ // pressure due to gravity
+ for (int phaseIdx=0; phaseIdx < numPhases; phaseIdx++)
+ {
+ int i = face().i;
+ int j = face().j;
+ Scalar fI = rhoFactor_(phaseIdx, i, elemVolVars);
+ Scalar fJ = rhoFactor_(phaseIdx, j, elemVolVars);
+ if (fI + fJ <= 0)
+ fI = fJ = 0.5; // doesn't matter because no phase is
+ // present in both cells!
+ density_[phaseIdx] =
+ (fI*elemVolVars[i].density(phaseIdx) +
+ fJ*elemVolVars[j].density(phaseIdx))
+ /
+ (fI + fJ);
+ // phase density
+ molarDensity_[phaseIdx]
+ =
+ (fI*elemVolVars[i].molarDensity(phaseIdx) +
+ fJ*elemVolVars[j].molarDensity(phaseIdx))
+ /
+ (fI + fJ); //arithmetic averaging
+
+ tmp = problem.gravity();
+ tmp *= density_[phaseIdx];
+
+ potentialGrad_[phaseIdx] -= tmp;
+ }
+ }
+
+ Scalar rhoFactor_(int phaseIdx, int scvIdx, const ElementVolumeVariables &vDat)
+ {
+
+ static const Scalar eps = 1e-2;
+ const Scalar sat = vDat[scvIdx].density(phaseIdx);
+ if (sat > eps)
+ return 0.5;
+ if (sat <= 0)
+ return 0;
+
+ static const Dumux::Spline<Scalar> sp(0, eps, // x0, x1
+ 0, 0.5, // y0, y1
+ 0, 0); // m0, m1
+ return sp.eval(sat);
+ }
+
+ 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++)
+ {
+ auto K_i = spatialParams.intrinsicPermeability(element,this->fvGeometry_,face().i);
+ //K_i *= volVarsI.permFactor();
+
+ auto K_j = spatialParams.intrinsicPermeability(element,this->fvGeometry_,face().j);
+ //K_j *= volVarsJ.permFactor();
+
+ spatialParams.meanK(K,K_i,K_j);
+
+ K.mv(potentialGrad_[phaseIdx], Kmvp_[phaseIdx]);
+ KmvpNormal_[phaseIdx] = - (Kmvp_[phaseIdx] * face().normal);
+ }
+
+ // set the upstream and downstream vertices
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ upstreamIdx_[phaseIdx] = face().i;
+ downstreamIdx_[phaseIdx] = face().j;
+
+ if (KmvpNormal_[phaseIdx] < 0) {
+ std::swap(upstreamIdx_[phaseIdx],
+ downstreamIdx_[phaseIdx]);
+ }
+ }
+ }
+
+ void calculateporousDiffCoeff_(const Problem &problem,
+ const Element &element,
+ const ElementVolumeVariables &elemVolVars)
+ {
+ const VolumeVariables &volVarsI = elemVolVars[face().i];
+ const VolumeVariables &volVarsJ = elemVolVars[face().j];
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ /* If there is no phase saturation on either side of the face
+ * no diffusion takes place */
+
+ if (volVarsI.saturation(phaseIdx) <= 0 ||
+ volVarsJ.saturation(phaseIdx) <= 0)
+ {
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ porousDiffCoeff_[phaseIdx][compIdx] = 0.0;
+ }
+ }
+
+ else
+ {
+ // calculate tortuosity at the nodes i and j needed
+ // for porous media diffusion coefficient
+ Scalar tauI = 1.0/(volVarsI.porosity() * volVarsI.porosity()) *
+ pow(volVarsI.porosity() * volVarsI.saturation(phaseIdx), 7.0/3);
+
+ Scalar tauJ = 1.0/(volVarsJ.porosity() * volVarsJ.porosity()) *
+ pow(volVarsJ.porosity() * volVarsJ.saturation(phaseIdx), 7.0/3);
+ // Diffusion coefficient in the porous medium
+
+ // -> harmonic mean
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ if(phaseIdx==compIdx)
+ porousDiffCoeff_[phaseIdx][compIdx] = 0.0;
+ else
+ {
+ porousDiffCoeff_[phaseIdx][compIdx] = harmonicMean(volVarsI.porosity() * volVarsI.saturation(phaseIdx) * tauI * volVarsI.diffCoeff(phaseIdx, compIdx),
+ volVarsJ.porosity() * volVarsJ.saturation(phaseIdx) * tauJ * volVarsJ.diffCoeff(phaseIdx, compIdx));
+ }
+ }
+ }
+ }
+ }
+
+public:
+ /*!
+ * \brief Return the pressure potential multiplied with the
+ * intrinsic permeability which goes from vertex i to
+ * vertex j.
+ *
+ * Note that the length of the face's normal is the area of the
+ * face, so this is not the actual velocity by the integral of
+ * the velocity over the face's area. Also note that the phase
+ * mobility is not yet included here since this would require a
+ * decision on the upwinding approach (which is done in the
+ * model and/or local residual file).
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar KmvpNormal(int phaseIdx) const
+ { return KmvpNormal_[phaseIdx]; }
+
+ /*!
+ * \brief Return the pressure potential multiplied with the
+ * intrinsic permeability as vector (for velocity output)
+ *
+ * \param phaseIdx The phase index
+ */
+ DimVector Kmvp(int phaseIdx) const
+ { return Kmvp_[phaseIdx]; }
+
+ /*!
+ * \brief Return the local index of the upstream control volume
+ * for a given phase.
+ *
+ * \param phaseIdx The phase index
+ */
+ int upstreamIdx(int phaseIdx) const
+ { return upstreamIdx_[phaseIdx]; }
+
+ /*!
+ * \brief Return the local index of the downstream control volume
+ * for a given phase.
+ *
+ * \param phaseIdx The phase index
+ */
+ int downstreamIdx(int phaseIdx) const
+ { return downstreamIdx_[phaseIdx]; }
+
+ /*!
+ * \brief The binary diffusion coefficient for each fluid phase.
+ *
+ * \param phaseIdx The phase index
+ * \param compIdx The component index
+ */
+ Scalar porousDiffCoeff(int phaseIdx, int compIdx) const
+ { return porousDiffCoeff_[phaseIdx][compIdx];}
+
+ /*!
+ * \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase at the integration
+ * point.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar density(int phaseIdx) const
+ { return density_[phaseIdx]; }
+
+ /*!
+ * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ of a phase at the integration
+ * point.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar molarDensity(int phaseIdx) const
+ { return molarDensity_[phaseIdx]; }
+
+ /*!
+ * \brief The concentration gradient of a component in a phase.
+ *
+ * \param phaseIdx The phase index
+ * \param compIdx The component index
+ */
+ const DimVector &massFractionGrad(int phaseIdx, int compIdx) const
+ { return massFractionGrad_[phaseIdx][compIdx]; }
+
+ /*!
+ * \brief The molar concentration gradient of a component in a phase.
+ *
+ * \param phaseIdx The phase index
+ * \param compIdx The component index
+ */
+ const DimVector &moleFractionGrad(int phaseIdx, int compIdx) const
+ { return moleFractionGrad_[phaseIdx][compIdx]; }
+
+ const SCVFace &face() const
+ {
+ if (this->onBoundary_)
+ return this->fvGeometry_.boundaryFace[this->faceIdx_];
+ else
+ return this->fvGeometry_.subContVolFace[this->faceIdx_];
+ }
+
+protected:
+
+ // gradients
+ DimVector potentialGrad_[numPhases];
+ DimVector massFractionGrad_[numPhases][numComponents];
+ DimVector moleFractionGrad_[numPhases][numComponents];
+
+ // density of each face at the integration point
+ Scalar density_[numPhases], molarDensity_[numPhases];
+
+ // intrinsic permeability times pressure potential gradient
+ DimVector Kmvp_[numPhases];
+ // projected on the face normal
+ Scalar KmvpNormal_[numPhases];
+
+ // local index of the upwind vertex for each phase
+ int upstreamIdx_[numPhases];
+ // local index of the downwind vertex for each phase
+ int downstreamIdx_[numPhases];
+
+ // the diffusion coefficient for the porous medium
+ Dune::FieldMatrix<Scalar, numPhases, numComponents> porousDiffCoeff_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/2pnc/2pncindices.hh b/dumux/implicit/2pnc/2pncindices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..671917a14488bd82a5f8d39c79cdd28ba09a5364
--- /dev/null
+++ b/dumux/implicit/2pnc/2pncindices.hh
@@ -0,0 +1,80 @@
+// -*- 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
+ * fully implicit model.
+ */
+#ifndef DUMUX_2PNC_INDICES_HH
+#define DUMUX_2PNC_INDICES_HH
+#include "2pncproperties.hh"
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPNCModel
+ * \ingroup ImplicitIndices
+ * \brief Enumerates the formulations which the two-phase n-component model accepts.
+ *
+ */
+struct TwoPNCFormulation//TODO: This might need to be change similar to 2p2c indices
+{
+ enum {
+ plSg,
+ pgSl,
+ pnSw = pgSl,
+ pwSn = plSg
+ };
+};
+
+/*!
+ * \ingroup TwoPNCModel
+ * \ingroup ImplicitIndices
+ * \brief The indices for the isothermal two-phase n-component model.
+ *
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int PVOffset = 0>
+class TwoPNCIndices
+{
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+public:
+ // Phase indices
+ static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< Index of the wetting phase
+ static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< Index of the non-wetting phase
+ // present phases (-> 'pseudo' primary variable)
+ static const int wPhaseOnly = 1; //!< Only the non-wetting phase is present
+ static const int nPhaseOnly = 2; //!< Only the wetting phase is present
+ static const int bothPhases = 3; //!< Both phases are present
+
+ // Primary variable indices
+ static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
+ static const int switchIdx = PVOffset + 1; //!< Index of the either the saturation or the mass fraction of the non-wetting/wetting phase
+ // equation indices
+ static const int conti0EqIdx = PVOffset + 0; //!< Reference index for mass conservation equations.
+ static const int contiWEqIdx = conti0EqIdx + FluidSystem::wCompIdx; //!< Index of the mass conservation equation for the wetting phase major component
+ static const int contiNEqIdx = conti0EqIdx + FluidSystem::nCompIdx; //!< Index of the mass conservation equation for the non-wetting phase major component
+};
+
+// \}
+
+}
+
+#endif
diff --git a/dumux/implicit/2pnc/2pnclocalresidual.hh b/dumux/implicit/2pnc/2pnclocalresidual.hh
new file mode 100644
index 0000000000000000000000000000000000000000..43755d5731123449c1672ac274b78bff4002a172
--- /dev/null
+++ b/dumux/implicit/2pnc/2pnclocalresidual.hh
@@ -0,0 +1,370 @@
+// -*- 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_2PNC_LOCAL_RESIDUAL_BASE_HH
+#define DUMUX_2PNC_LOCAL_RESIDUAL_BASE_HH
+
+#include "2pncproperties.hh"
+#include "2pncvolumevariables.hh"
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+#include <iostream>
+#include <vector>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPNCModel
+ * \ingroup ImplicitLocalResidual
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the two-phase n-component 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 TwoPNCLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
+{
+protected:
+ typedef TwoPNCLocalResidual<TypeTag> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
+ typedef BoxLocalResidual<TypeTag> ParentType;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+
+ 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 CompositionalFluidState<Scalar, FluidSystem> FluidState;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum
+ {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ numEq = GET_PROP_VALUE(TypeTag, NumEq),
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ 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;
+
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ typedef typename GridView::ctype CoordScalar;
+
+
+public:
+ /*!
+ * \brief Constructor. Sets the upwind weight.
+ */
+ TwoPNCLocalResidual()
+ {
+ // 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
+ massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
+ };
+
+ /*!
+ * \brief Evaluate the storage term of the current solution in a
+ * single phase.
+ *
+ * \param element The element
+ * \param phaseIdx The index of the fluid phase
+ */
+ void evalPhaseStorage(const Element &element, int phaseIdx)
+ {
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(this->gridView_(), element);
+ ElementBoundaryTypes bcTypes;
+ bcTypes.update(this->problem_(), element, fvGeometry);
+ ElementVolumeVariables volVars;
+ volVars.update(this->problem_(), element, fvGeometry, false);
+
+ this->residual_.resize(fvGeometry.numScv);
+ this->residual_ = 0;
+
+ this->elemPtr_ = &element;
+ this->fvElemGeomPtr_ = &fvGeometry;
+ this->bcTypesPtr_ = &bcTypes;
+ this->prevVolVarsPtr_ = 0;
+ this->curVolVarsPtr_ = &volVars;
+ evalPhaseStorage_(phaseIdx);
+ }
+
+ /*!
+ * \brief Evaluate the amount all conservation quantities
+ * (e.g. phase mass) within a sub-control volume.
+ *
+ * The result should be averaged over the volume (e.g. phase mass
+ * inside a sub control volume divided by the volume)
+ *
+ * \param 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
+ {
+ // if flag usePrevSol is set, the solution from the previous
+ // time step is used, otherwise the current solution is
+ // used. The secondary variables are used accordingly. This
+ // is required to compute the derivative of the storage term
+ // using the implicit euler method.
+ const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_()
+ : this->curVolVars_();
+ const VolumeVariables &volVars = elemVolVars[scvIdx];
+
+ // Compute storage term of all fluid components in the fluid phases
+ storage = 0;
+ for (unsigned int phaseIdx = 0; phaseIdx < numPhases /*+ numSPhases*/; ++phaseIdx)
+ {
+ //if(phaseIdx< numPhases)
+ //{
+ for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx) //H2O, Air, Salt
+ {
+ int eqIdx = conti0EqIdx + compIdx;
+ if (replaceCompEqIdx != eqIdx)
+ {
+ storage[eqIdx] += volVars.molarDensity(phaseIdx)
+ * volVars.saturation(phaseIdx)
+ * volVars.fluidState().moleFraction(phaseIdx, compIdx)
+ * volVars.porosity();
+ }
+ else
+ {
+ storage[replaceCompEqIdx] += volVars.molarDensity(phaseIdx)
+ * volVars.saturation(phaseIdx)
+ * volVars.porosity();
+ }
+ }
+ }
+ Valgrind::CheckDefined(storage);
+ }
+ /*!
+ * \brief Evaluates the total flux of all conservation quantities
+ * over a face of a sub-control volume.
+ *
+ * \param flux The flux over the sub-control-volume face for each component
+ * \param fIdx The index of the sub-control-volume face
+ * \param onBoundary Evaluate flux at inner sub-control-volume face or on a boundary face
+ */
+ void computeFlux(PrimaryVariables &flux, const int fIdx, bool onBoundary = false) const
+ {
+ FluxVariables fluxVars(this->problem_(),
+ this->element_(),
+ this->fvGeometry_(),
+ fIdx,
+ this->curVolVars_(),
+ onBoundary);
+
+ flux = 0;
+ asImp_()->computeAdvectiveFlux(flux, fluxVars);
+ asImp_()->computeDiffusiveFlux(flux, fluxVars);
+ Valgrind::CheckDefined(flux);
+ }
+ /*!
+ * \brief Evaluates the advective mass flux of all components over
+ * a face of a subcontrol volume.
+ *
+ * \param flux The advective flux over the sub-control-volume face for each component
+ * \param fluxVars The flux variables at the current SCV
+ */
+ void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
+ {
+ ////////
+ // advective fluxes of all components in all phases
+ ////////
+ for (unsigned int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ // data attached to upstream and the downstream vertices
+ // of the current phase
+ const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
+ const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
+
+ for (unsigned int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ // add advective flux of current component in current
+ // phase
+ unsigned int eqIdx = conti0EqIdx + compIdx;
+
+ if (eqIdx != replaceCompEqIdx)
+ {
+ // upstream vertex
+ flux[eqIdx] += fluxVars.KmvpNormal(phaseIdx)
+ * (massUpwindWeight_
+ * up.mobility(phaseIdx)
+ * up.molarDensity(phaseIdx)
+ * up.fluidState().moleFraction(phaseIdx, compIdx)
+ +
+ (1.0 - massUpwindWeight_)
+ * dn.mobility(phaseIdx)
+ * dn.molarDensity(phaseIdx)
+ * dn.fluidState().moleFraction(phaseIdx, compIdx));
+
+ Valgrind::CheckDefined(fluxVars.KmvpNormal(phaseIdx));
+ Valgrind::CheckDefined(up.molarDensity(phaseIdx));
+ Valgrind::CheckDefined(dn.molarDensity(phaseIdx));
+ }
+ else
+ {
+ flux[replaceCompEqIdx] += fluxVars.KmvpNormal(phaseIdx)
+ * (massUpwindWeight_
+ * up.molarDensity(phaseIdx)
+ * up.mobility(phaseIdx)
+ +
+ (1.0 - massUpwindWeight_)
+ * dn.molarDensity(phaseIdx)
+ * dn.mobility(phaseIdx));
+
+
+ Valgrind::CheckDefined(fluxVars.KmvpNormal(phaseIdx));
+ Valgrind::CheckDefined(up.molarDensity(phaseIdx));
+ Valgrind::CheckDefined(dn.molarDensity(phaseIdx));
+ }
+ }
+ }
+ }
+
+ /*!
+ * \brief Evaluates the diffusive mass flux of all components over
+ * a face of a sub-control volume.
+ *
+ * \param flux The flux over the sub-control-volume face for each component
+ * \param fluxVars The flux variables at the current sub-control-volume face
+ */
+ void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
+ {
+ //Loop calculates the diffusive flux for every component in a phase. The amount of moles of a component
+ //(eg Air in liquid) in a phase
+ //which is not the main component (eg. H2O in the liquid phase) moved from i to j equals the amount of moles moved
+ //from the main component in a phase (eg. H2O in the liquid phase) from j to i. So two fluxes in each component loop
+ // are calculated in the same phase.
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ //add diffusive fluxes only to the component balances
+ if (replaceCompEqIdx != (conti0EqIdx + compIdx))
+ {
+ Scalar diffCont = - fluxVars.porousDiffCoeff(phaseIdx ,compIdx)
+ * fluxVars.molarDensity(phaseIdx)
+ * (fluxVars.moleFractionGrad(phaseIdx, compIdx)
+ * fluxVars.face().normal);
+ flux[conti0EqIdx + compIdx] += diffCont;
+ flux[conti0EqIdx + phaseIdx] -= diffCont;
+ }
+ }
+ }
+
+ /*!
+ * \brief Evaluates the source term
+ *
+ * \param source The source/sink in the sub-control volume
+ * \param scvIdx The index of the sub-control volume
+ *
+ * Be careful what you use, mole or mass fraction! Think of the units!
+ * If a total mass balance is used, the sum of both components has to be specified as source.
+ */
+ void computeSource(PrimaryVariables &source, int scvIdx)
+ {
+ this->problem_().solDependentSource(source,
+ this->element_(),
+ this->fvGeometry_(),
+ scvIdx,
+ this->curVolVars_());
+
+ Valgrind::CheckDefined(source);
+ }
+
+
+protected:
+
+ void evalPhaseStorage_(int phaseIdx)
+ {
+ // evaluate the storage terms of a single phase
+ for (int i=0; i < this->fvGeometry_().numScv; i++)
+ {
+ PrimaryVariables &result = this->residual_[i];
+ const ElementVolumeVariables &elemVolVars = this->curVolVars_();
+ const VolumeVariables &volVars = elemVolVars[i];
+
+ // compute storage term of all fluid components within all phases
+ result = 0;
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ result[conti0EqIdx + compIdx] += volVars.density(phaseIdx)
+ * volVars.saturation(phaseIdx)
+ * volVars.fluidState().massFraction(phaseIdx, compIdx)
+ * volVars.porosity();
+ }
+ result *= this->fvGeometry_().subContVol[i].volume;
+ }
+ }
+
+ Implementation *asImp_()
+ { return static_cast<Implementation *> (this); }
+ const Implementation *asImp_() const
+ { return static_cast<const Implementation *> (this); }
+
+public:
+ Scalar massUpwindWeight_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/2pnc/2pncmodel.hh b/dumux/implicit/2pnc/2pncmodel.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4d6e779d08aad9dd4af9cd7e5ecebc819daa3b8a
--- /dev/null
+++ b/dumux/implicit/2pnc/2pncmodel.hh
@@ -0,0 +1,819 @@
+// -*- 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_2PNC_MODEL_HH
+#define DUMUX_2PNC_MODEL_HH
+
+#include <dune/common/version.hh>
+// #include <dumux/material/constants.hh>
+
+#include "2pncproperties.hh"
+#include "2pncindices.hh"
+#include "2pnclocalresidual.hh"
+#include <dumux/implicit/common/implicitvelocityoutput.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPNCModel
+ * \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 TwoPNCModel: public GET_PROP_TYPE(TypeTag, BaseModel)
+{
+ typedef TwoPNCModel<TypeTag> ThisType;
+ typedef typename GET_PROP_TYPE(TypeTag, BaseModel) 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),
+ numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
+ 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 CompositionalFluidState<Scalar, FluidSystem> FluidState;
+
+ 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 Initialize the static data with the initial solution.
+ *
+ * \param problem The problem to be solved
+ */
+ void init(Problem &problem)
+ {
+ ParentType::init(problem);
+
+ unsigned numDofs = this->numDofs();
+
+ staticDat_.resize(numDofs);
+
+ setSwitched_(false);
+
+ ElementIterator eIt = this->gridView_().template begin<0>();
+ ElementIterator elemEndIt = this->gridView_().template end<0>();
+ for (; eIt != elemEndIt; ++eIt)
+ {
+ if (!isBox) // i.e. cell-centered discretization
+ {
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
+ int eIdxGlobal = this->dofMapper().index(*eIt);
+#else
+ int eIdxGlobal = this->dofMapper().map(*eIt);
+#endif
+ const GlobalPosition &globalPos = eIt->geometry().center();
+
+ // initialize phase presence
+ staticDat_[eIdxGlobal].phasePresence
+ = this->problem_().initialPhasePresence(*(this->gridView_().template begin<dim>()),
+ eIdxGlobal, globalPos);
+ staticDat_[eIdxGlobal].wasSwitched = false;
+
+ staticDat_[eIdxGlobal].oldPhasePresence
+ = staticDat_[eIdxGlobal].phasePresence;
+ }
+ }
+
+ if (isBox) // i.e. vertex-centered discretization
+ {
+ VertexIterator vIt = this->gridView_().template begin<dim> ();
+ const VertexIterator &vEndIt = this->gridView_().template end<dim> ();
+ for (; vIt != vEndIt; ++vIt)
+ {
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
+ int vIdxGlobal = this->dofMapper().index(*vIt);
+#else
+ int vIdxGlobal = this->dofMapper().map(*vIt);
+#endif
+ const GlobalPosition &globalPos = vIt->geometry().corner(0);
+
+ // initialize phase presence
+ staticDat_[vIdxGlobal].phasePresence
+ = this->problem_().initialPhasePresence(*vIt, vIdxGlobal,
+ globalPos);
+ staticDat_[vIdxGlobal].wasSwitched = false;
+
+ staticDat_[vIdxGlobal].oldPhasePresence
+ = staticDat_[vIdxGlobal].phasePresence;
+ }
+ }
+ }
+
+ /*!
+ * \brief Compute the total storage of all conservation quantities in one phase
+ *
+ * \param storage Contains the storage of each component in one phase
+ * \param phaseIdx The phase index
+ */
+ void globalPhaseStorage(PrimaryVariables &storage, int phaseIdx)
+ {
+ storage = 0;
+
+ ElementIterator eIt = this->gridView_().template begin<0>();
+ const ElementIterator elemEndIt = this->gridView_().template end<0>();
+ for (; eIt != elemEndIt; ++eIt)
+ {
+ if(eIt->partitionType() == Dune::InteriorEntity)
+ {
+ this->localResidual().evalPhaseStorage(*eIt, phaseIdx);
+
+ for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
+ storage += this->localResidual().storageTerm()[i];
+ }
+ }
+
+ this->gridView_().comm().sum(storage);
+ }
+
+ /*!
+ * \brief Called by the update() method if applying the newton
+ * method was unsuccessful.
+ */
+ void updateFailed()
+ {
+ ParentType::updateFailed();
+
+ setSwitched_(false);
+ resetPhasePresence_();
+ }
+
+ /*!
+ * \brief Called by the problem if a time integration was
+ * successful, post processing of the solution is done and the
+ * result has been written to disk.
+ *
+ * This should prepare the model for the next time integration.
+ */
+ void advanceTimeLevel()
+ {
+ ParentType::advanceTimeLevel();
+
+ // update the phase state
+ updateOldPhasePresence_();
+ setSwitched_(false);
+ }
+
+ /*!
+ * \brief Return true if the primary variables were switched for
+ * at least one vertex after the last timestep.
+ */
+ bool switched() const
+ {
+ return switchFlag_;
+ }
+
+ /*!
+ * \brief Returns the phase presence of the current or the old solution of a vertex.
+ *
+ * \param globalVertexIdx The global vertex index
+ * \param oldSol Evaluate function with solution of current or previous time step
+ */
+ int phasePresence(int globalVertexIdx, bool oldSol) const
+ {
+ return oldSol ? staticDat_[globalVertexIdx].oldPhasePresence
+ : staticDat_[globalVertexIdx].phasePresence;
+ }
+
+ /*!
+ * \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>
+ 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 *currentDensity = writer.allocateManagedBuffer (numDofs);
+ ScalarField *reactionSourceH2O = writer.allocateManagedBuffer (numDofs);
+ ScalarField *reactionSourceO2 = writer.allocateManagedBuffer (numDofs);
+
+ 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 *temperature = writer.allocateManagedBuffer (numDofs);
+ ScalarField *poro = writer.allocateManagedBuffer (numDofs);
+ ScalarField *boxVolume = writer.allocateManagedBuffer (numDofs);
+ 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);
+ }
+ }
+
+ ScalarField *moleFraction[numPhases][numComponents];
+ for (int i = 0; i < numPhases; ++i)
+ for (int j = 0; j < numComponents; ++j)
+ moleFraction[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;
+
+ unsigned numElements = this->gridView_().size(0);
+ ScalarField *rank =
+ writer.allocateManagedBuffer (numElements);
+
+ FVElementGeometry fvGeometry;
+ VolumeVariables volVars;
+ ElementVolumeVariables elemVolVars;
+
+ ElementIterator eIt = this->gridView_().template begin<0>();
+ ElementIterator eEndIt = this->gridView_().template end<0>();
+ for (; eIt != eEndIt; ++eIt)
+ {
+ int idx = this->problem_().elementMapper().map(*eIt);
+ (*rank)[idx] = this->gridView_().comm().rank();
+ fvGeometry.update(this->gridView_(), *eIt);
+
+ elemVolVars.update(this->problem_(),
+ *eIt,
+ fvGeometry,
+ false /* oldSol? */);
+
+#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
+ int numVerts = eIt->subEntities(dim);
+#else
+ int numVerts = eIt->template count<dim> ();
+#endif
+ for (int i = 0; i < numVerts; ++i)
+ {
+ int globalIdx = this->vertexMapper().map(*eIt, i, dim);
+ volVars.update(sol[globalIdx],
+ this->problem_(),
+ *eIt,
+ fvGeometry,
+ i,
+ false);
+
+ GlobalPosition globalPos = fvGeometry.subContVol[i].global;
+ //Scalar zmax = this->problem_().bBoxMax()[dim-1];
+ (*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();
+
+ if(useElectrochem){
+
+ //for electrochemistry output we need elemVolVars
+ elemVolVars.update(this->problem_(), *eIt, fvGeometry,/*oldSol=*/false);
+
+ //reactionSource Output
+ PrimaryVariables source;
+ this->problem_().solDependentSource(source, *eIt, fvGeometry, i, elemVolVars);
+
+ (*reactionSourceH2O)[globalIdx] = source[wPhaseIdx];
+ (*reactionSourceO2)[globalIdx] = source[numComponents-1];
+
+ //Current Output
+ (*currentDensity)[globalIdx] = -1.0*source[numComponents-1]*4*Constant::F;
+ //recorrection of the area for output
+ Scalar gridYMin = GET_RUNTIME_PARAM(TypeTag, Scalar, Grid.yMin);
+ Scalar gridYMax = GET_RUNTIME_PARAM(TypeTag, Scalar, Grid.yMax);
+ Scalar nCellsY = GET_RUNTIME_PARAM(TypeTag, Scalar, Grid.CellsY);
+
+ Scalar lengthBox= (gridYMax - gridYMin)/nCellsY;
+
+ (*currentDensity)[globalIdx] = (*currentDensity)[globalIdx]/2*lengthBox/10000; //i in [A/cm^2], cf. paper Oliver
+ }
+
+ (*temperature)[globalIdx] = volVars.temperature();
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ (*moleFraction[phaseIdx][compIdx])[globalIdx]= volVars.fluidState().moleFraction(phaseIdx,compIdx);
+
+ Valgrind::CheckDefined((*moleFraction[phaseIdx][compIdx])[globalIdx]);
+
+ }
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ (*molarity[compIdx])[globalIdx] = (volVars.fluidState().molarity(wPhaseIdx, compIdx));
+
+ Tensor K = perm_(this->problem_().spatialParams().intrinsicPermeability(*eIt, 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, *eIt, wPhaseIdx);
+ velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, 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(*temperature, "temperature");
+ writer.attachVertexData(*boxVolume, "boxVolume");
+
+ if(useElectrochem){
+ writer.attachVertexData(*reactionSourceH2O, "reactionSourceH2O [mol/(sm^2)]");
+ writer.attachVertexData(*reactionSourceO2, "reactionSourceO2 [mol/(sm^2)]");
+ writer.attachVertexData(*currentDensity, "currentDensity [A/cm^2]");
+ }
+
+ 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::componentName(j)
+ << FluidSystem::phaseName(i);
+ writer.attachVertexData(*moleFraction[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 Write the current solution to a restart file.
+ *
+ * \param outStream The output stream of one vertex for the restart file
+ * \param entity The Entity
+ */
+ template<class Entity>
+ void serializeEntity(std::ostream &outStream, const Entity &entity)
+ {
+ // write primary variables
+ ParentType::serializeEntity(outStream, entity);
+
+ int vertIdx = this->dofMapper().map(entity);
+ if (!outStream.good())
+ DUNE_THROW(Dune::IOError, "Could not serialize vertex " << vertIdx);
+
+ outStream << staticDat_[vertIdx].phasePresence << " ";
+ }
+
+ /*!
+ * \brief Reads the current solution for a vertex from a restart
+ * file.
+ *
+ * \param inStream The input stream of one vertex from the restart file
+ * \param entity The Entity
+ */
+ template<class Entity>
+ void deserializeEntity(std::istream &inStream, const Entity &entity)
+ {
+ // read primary variables
+ ParentType::deserializeEntity(inStream, entity);
+
+ // read phase presence
+ int vertIdx = this->dofMapper().map(entity);
+ if (!inStream.good())
+ DUNE_THROW(Dune::IOError,
+ "Could not deserialize vertex " << vertIdx);
+
+ inStream >> staticDat_[vertIdx].phasePresence;
+ staticDat_[vertIdx].oldPhasePresence
+ = staticDat_[vertIdx].phasePresence;
+
+ }
+
+ /*!
+ * \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 < staticDat_.size(); ++i)
+ 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 (staticDat_[globalIdx].visited)
+ continue;
+
+ 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 Data which is attached to each vertex and is not only
+ * stored locally.
+ */
+ struct StaticVars
+ {
+ int phasePresence;
+ bool wasSwitched;
+
+ int oldPhasePresence;
+ bool visited;
+ };
+
+ Tensor perm_(Scalar perm)
+ {
+ Tensor K(0.0);
+
+ for(int i=0; i<dim; i++)
+ K[i][i] = perm;
+
+ return K;
+ }
+
+ Tensor perm_(Tensor perm)
+ {
+ return perm;
+ }
+
+ /*!
+ * \brief Reset the current phase presence of all vertices to the old one.
+ *
+ * This is done after an update failed.
+ */
+ void resetPhasePresence_()
+ {
+ int numDofs = this->gridView_().size(dim);
+ for (int i = 0; i < numDofs; ++i)
+ {
+ staticDat_[i].phasePresence
+ = staticDat_[i].oldPhasePresence;
+ staticDat_[i].wasSwitched = false;
+ }
+ }
+
+ /*!
+ * \brief Set the old phase of all verts state to the current one.
+ */
+ void updateOldPhasePresence_()
+ {
+ int numDofs = this->gridView_().size(dim);
+ for (int i = 0; i < numDofs; ++i)
+ {
+ staticDat_[i].oldPhasePresence
+ = staticDat_[i].phasePresence;
+ staticDat_[i].wasSwitched = false;
+ }
+ }
+
+ /*!
+ * \brief Set whether there was a primary variable switch after in
+ * the last timestep.
+ */
+ void setSwitched_(bool yesno)
+ {
+ switchFlag_ = yesno;
+ }
+
+ // perform variable switch at a vertex; Returns true if a
+ // variable switch was performed.
+ bool primaryVarSwitch_(SolutionVector &globalSol,
+ const VolumeVariables &volVars, int globalIdx,
+ const GlobalPosition &globalPos)
+ {
+
+ // evaluate primary variable switch
+ bool wouldSwitch = false;
+ int phasePresence = staticDat_[globalIdx].phasePresence;
+ int newPhasePresence = phasePresence;
+
+ //check if a primary variable switch is necessary
+ if (phasePresence == bothPhases)
+ {
+ Scalar Smin = 0; //saturation threshold
+ if (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*/);
+
+ //switch all secondary components to mole fraction in gas phase
+ for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
+ globalSol[globalIdx][compIdx] = volVars.fluidState().moleFraction(nPhaseIdx,compIdx);
+ }
+ //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 xlmax = 1;
+ Scalar sumxl = 0;
+ //Calculate sum of mole fractions in the hypothetical liquid phase
+ for (int compIdx = 0; compIdx < numComponents; compIdx++)
+ {
+ sumxl += volVars.fluidState().moleFraction(wPhaseIdx, compIdx);
+ }
+ if (sumxl > xlmax)
+ wouldSwitch = true;
+ if (staticDat_[globalIdx].wasSwitched)
+ xlmax *=1.02;
+ //liquid phase appears if sum is larger than one
+ if (sumxl/*sum of mole fractions*/ > xlmax/*1*/)
+ {
+ std::cout << "Liquid Phase appears at vertex " << globalIdx
+ << ", coordinated: " << globalPos << ", sumxl: "
+ << sumxl << std::endl;
+ newPhasePresence = bothPhases;
+
+ //saturation of the liquid phase set to 0.0001 (if formulation pgSl and vice versa)
+ if (formulation == pgSl)
+ globalSol[globalIdx][switchIdx] = 0.0001;
+ else if (formulation == plSg)
+ globalSol[globalIdx][switchIdx] = 0.9999;
+
+ //switch all secondary components back to liquid mole fraction
+ for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
+ globalSol[globalIdx][compIdx] = volVars.fluidState().moleFraction(wPhaseIdx,compIdx);
+ }
+ }
+ 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 < numComponents; compIdx++)
+ {
+ sumxg += volVars.fluidState().moleFraction(nPhaseIdx, compIdx);
+ }
+ if (sumxg > xgmax)
+ wouldSwitch = true;
+ if (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.9999;
+ else if (formulation == plSg)
+ globalSol[globalIdx][switchIdx] = 0.0001;
+
+ }
+ }
+
+
+ staticDat_[globalIdx].phasePresence = newPhasePresence;
+ staticDat_[globalIdx].wasSwitched = wouldSwitch;
+ return phasePresence != newPhasePresence;
+ }
+
+ // parameters given in constructor
+ std::vector<StaticVars> staticDat_;
+ bool switchFlag_;
+};
+
+}
+
+#include "2pncpropertydefaults.hh"
+
+#endif
diff --git a/dumux/implicit/2pnc/2pncnewtoncontroller.hh b/dumux/implicit/2pnc/2pncnewtoncontroller.hh
new file mode 100644
index 0000000000000000000000000000000000000000..2f970c83b100a4a0c9589bb22b0e7d695183f3e5
--- /dev/null
+++ b/dumux/implicit/2pnc/2pncnewtoncontroller.hh
@@ -0,0 +1,105 @@
+// -*- 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 A two-phase n-component specific controller for the newton solver.
+ */
+#ifndef DUMUX_2PNC_NEWTON_CONTROLLER_HH
+#define DUMUX_2PNC_NEWTON_CONTROLLER_HH
+
+#include "2pncproperties.hh"
+
+#include <dumux/nonlinear/newtoncontroller.hh>
+
+namespace Dumux {
+/*!
+ * \ingroup Newton
+ * \ingroup TwoPNCModel
+ * \brief A two-phase n-component specific controller for the newton solver.
+ *
+ * This controller 'knows' what a 'physically meaningful' solution is
+ * which allows the newton method to abort quicker if the solution is
+ * way out of bounds.
+ */
+template <class TypeTag>
+class TwoPNCNewtonController : public NewtonController<TypeTag>
+{
+ typedef NewtonController<TypeTag> ParentType;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
+
+public:
+ TwoPNCNewtonController(Problem &problem)
+ : ParentType(problem)
+ {};
+
+
+ /*!
+ * \brief
+ * Suggest a new time step size based either on the number of newton
+ * iterations required or on the variable switch
+ *
+ * \param uCurrentIter The current global solution vector
+ * \param uLastIter The previous global solution vector
+ *
+ */
+ void newtonEndStep(SolutionVector &uCurrentIter,
+ const SolutionVector &uLastIter)
+ {
+ int succeeded;
+ try {
+ // call the method of the base class
+ this->method().model().updateStaticData(uCurrentIter, uLastIter);
+ ParentType::newtonEndStep(uCurrentIter, uLastIter);
+
+ succeeded = 1;
+ if (this->gridView_().comm().size() > 1)
+ succeeded = this->gridView_().comm().min(succeeded);
+ }
+ catch (Dumux::NumericalProblem &e)
+ {
+ std::cout << "rank " << this->problem_().gridView().comm().rank()
+ << " caught an exception while updating:" << e.what()
+ << "\n";
+ succeeded = 0;
+ if (this->gridView_().comm().size() > 1)
+ succeeded = this->gridView_().comm().min(succeeded);
+ }
+
+ if (!succeeded) {
+ DUNE_THROW(NumericalProblem,
+ "A process did not succeed in linearizing the system");
+ }
+ }
+
+ /*!
+ * \brief Returns true if the current solution can be considered to
+ * be accurate enough
+ */
+ bool newtonConverged()
+ {
+ if (this->method().model().switched())
+ return false;
+
+ return ParentType::newtonConverged();
+ }
+};
+}
+
+#endif
diff --git a/dumux/implicit/2pnc/2pncproperties.hh b/dumux/implicit/2pnc/2pncproperties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..42c3ed31b5d5b1e2bb3125223eaec9a271f69837
--- /dev/null
+++ b/dumux/implicit/2pnc/2pncproperties.hh
@@ -0,0 +1,73 @@
+// -*- 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 TwoPNCModel
+ *
+ * \file
+ *
+ * \brief Defines the properties required for the two-phase n-component
+ * fully implicit model.
+ */
+#ifndef DUMUX_2PNC_PROPERTIES_HH
+#define DUMUX_2PNC_PROPERTIES_HH
+
+#include <dumux/implicit/box/boxproperties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tag for the isothermal two phase n component mineralisation problems
+NEW_TYPE_TAG(BoxTwoPNC, INHERITS_FROM(BoxModel));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
+NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system
+NEW_PROP_TAG(NumMajorComponents); //!< Number of major fluid components which are considered in the calculation of the phase density
+NEW_PROP_TAG(TwoPNCIndices); //!< Enumerations for the 2pncMin models
+NEW_PROP_TAG(Formulation); //!< The formulation of the model
+NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters
+NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
+NEW_PROP_TAG(Indices); //!< Enumerations for the model
+NEW_PROP_TAG(Chemistry); //!< The chemistry class with which solves equlibrium reactions
+NEW_PROP_TAG(useElectrochem); //!< Determines if Output for Electrochemistry is needed
+
+NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters)
+NEW_PROP_TAG(MaterialLawParams); //!< The parameters of the material law (extracted from the spatial parameters)
+
+NEW_PROP_TAG(ReplaceCompEqIdx); //!< The index of the total mass balance equation, if one component balance is replaced (ReplaceCompEqIdx < NumComponents)
+NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
+NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
+NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the upwind weight for the mass conservation equations
+NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< The value of the upwind parameter for the mobility
+NEW_PROP_TAG(BaseFluxVariables); //! The base flux variables
+}
+}
+
+#endif
diff --git a/dumux/implicit/2pnc/2pncpropertydefaults.hh b/dumux/implicit/2pnc/2pncpropertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..c5620b4d4217b455cfaa3084b728689d2cd7d4a8
--- /dev/null
+++ b/dumux/implicit/2pnc/2pncpropertydefaults.hh
@@ -0,0 +1,173 @@
+// -*- 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 TwoPNCModel
+ * \file
+ *
+ * \brief Defines default values for most properties required by the
+ * two-phase n-component fully implicit model.
+ */
+#ifndef DUMUX_2PNC_PROPERTY_DEFAULTS_HH
+#define DUMUX_2PNC_PROPERTY_DEFAULTS_HH
+
+#include "2pncindices.hh"
+#include "2pncmodel.hh"
+#include "2pncindices.hh"
+#include "2pncfluxvariables.hh"
+#include "2pncvolumevariables.hh"
+#include "2pncproperties.hh"
+#include "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 components.
+ *
+ * We just forward the number from the fluid system
+ *
+ */
+SET_PROP(BoxTwoPNC, NumComponents)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numComponents;
+
+};
+//! Set as default that no component mass balance is replaced by the total mass balance
+SET_PROP(BoxTwoPNC, ReplaceCompEqIdx)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numComponents;
+};
+//! The major components belonging to the existing phases are mentioned here e.g., 2 for water and air being the major component in the liquid and gas phases in a 2 phase system
+SET_PROP(BoxTwoPNC, NumMajorComponents)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numPhases;
+ static_assert(value == 2,
+ "The model is restricted to two-phases, thus number of major components must also be two.");
+};
+
+/*!
+ * \brief Set the property for the number of fluid phases.
+ *
+ * We just forward the number from the fluid system and use an static
+ * assert to make sure it is 2.
+ */
+SET_PROP(BoxTwoPNC, NumPhases)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numPhases;
+ static_assert(value == 2,
+ "Only fluid systems with 2 fluid phases are supported by the 2p-nc model!");
+};
+/*!
+ * \brief Set the property for the number of equations: For each existing component one equation has to be solved.
+ */
+SET_PROP(BoxTwoPNC, NumEq)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(FluidSystem)) FluidSystem;
+
+public:
+ static const int value = FluidSystem::numComponents;
+};
+
+//! Set the default formulation to pl-Sg: This can be over written in the problem.
+SET_INT_PROP(BoxTwoPNC, Formulation, TwoPNCFormulation::plSg);
+
+//! Set the property for the material parameters by extracting it from the material law.
+SET_PROP(BoxTwoPNC, MaterialLawParams)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, PTAG(MaterialLaw)) MaterialLaw;
+
+public:
+ typedef typename MaterialLaw::Params type;
+};
+
+//! Use the 2pnc local residual operator
+SET_TYPE_PROP(BoxTwoPNC,
+ LocalResidual,
+ TwoPNCLocalResidual<TypeTag>);
+
+//! Use the 2pnc newton controller
+SET_TYPE_PROP(BoxTwoPNC, NewtonController, TwoPNCNewtonController<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(BoxTwoPNC, Model, TwoPNCModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(BoxTwoPNC, VolumeVariables, TwoPNCVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(BoxTwoPNC, FluxVariables, TwoPNCFluxVariables<TypeTag>);
+
+//! define the base flux variables to realize Darcy flow
+SET_TYPE_PROP(BoxTwoPNC, BaseFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
+
+//! the upwind weight for the mass conservation equations.
+SET_SCALAR_PROP(BoxTwoPNC, ImplicitMassUpwindWeight, 1.0);
+
+//! Set default mobility upwind weight to 1.0, i.e. fully upwind
+SET_SCALAR_PROP(BoxTwoPNC, ImplicitMobilityUpwindWeight, 1.0);
+
+//! The indices required by the isothermal 2pnc model
+SET_TYPE_PROP(BoxTwoPNC, Indices, TwoPNCIndices <TypeTag, /*PVOffset=*/0>);
+
+//! Use the ImplicitSpatialParams by default
+SET_TYPE_PROP(BoxTwoPNC, SpatialParams, ImplicitSpatialParams<TypeTag>);
+
+//! Enable gravity by default
+SET_BOOL_PROP(BoxTwoPNC, ProblemEnableGravity, true);
+
+//! Disable velocity output by default
+SET_BOOL_PROP(BoxTwoPNC, VtkAddVelocity, false);
+
+//! disable electro-chemistry by default
+SET_BOOL_PROP(BoxTwoPNC, useElectrochem, false);
+
+}
+
+}
+
+#endif
diff --git a/dumux/implicit/2pnc/2pncvolumevariables.hh b/dumux/implicit/2pnc/2pncvolumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9e70e3915dd0781b8d2a56ebf53e4bfdafc9cc33
--- /dev/null
+++ b/dumux/implicit/2pnc/2pncvolumevariables.hh
@@ -0,0 +1,501 @@
+// -*- 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 model.
+ */
+#ifndef DUMUX_2PNC_VOLUME_VARIABLES_HH
+#define DUMUX_2PNC_VOLUME_VARIABLES_HH
+
+#include <iostream>
+#include <vector>
+#include <dune/common/parallel/collectivecommunication.hh>
+
+#include <dumux/implicit/common/implicitmodel.hh>
+#include <dumux/material/fluidstates/compositionalfluidstate.hh>
+#include <dumux/common/math.hh>
+
+#include "2pncproperties.hh"
+#include "2pncindices.hh"
+#include <dumux/material/constraintsolvers/computefromreferencephase2pnc.hh>
+#include <dumux/material/constraintsolvers/miscible2pnccomposition.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup TwoPNCModel
+ * \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 TwoPNCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
+{
+ typedef ImplicitVolumeVariables<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;
+ enum
+ {
+ dim = GridView::dimension,
+ dimWorld=GridView::dimensionworld,
+
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ 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,
+
+ };
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+ 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 &primaryVariables,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ int scvIdx,
+ bool isOldSol)
+ {
+ ParentType::update(primaryVariables,
+ problem,
+ element,
+ fvGeometry,
+ scvIdx,
+ isOldSol);
+
+ completeFluidState(primaryVariables, problem, element, fvGeometry, scvIdx, fluidState_, isOldSol);
+
+ /////////////
+ // calculate the remaining quantities
+ /////////////
+ const MaterialLawParams &materialParams = problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
+
+ // Second instance of a parameter cache.
+ // Could be avoided if diffusion coefficients also
+ // became part of the fluid state.
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updateAll(fluidState_);
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {// relative permeabilities
+ Scalar kr;
+ if (phaseIdx == wPhaseIdx)
+ kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx));
+ else // ATTENTION: krn requires the liquid saturation
+ // as parameter!
+ kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));
+ mobility_[phaseIdx] = kr / fluidState_.viscosity(phaseIdx);
+ Valgrind::CheckDefined(mobility_[phaseIdx]);
+ int compIIdx = phaseIdx;
+ for(int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ int compJIdx = compIdx;
+ // binary diffusion coefficents
+ diffCoeff_[phaseIdx][compIdx] = 0.0;
+ if(compIIdx!= compJIdx)
+ diffCoeff_[phaseIdx][compIdx] = FluidSystem::binaryDiffusionCoefficient(fluidState_,
+ paramCache,
+ phaseIdx,
+ compIIdx,
+ compJIdx);
+ Valgrind::CheckDefined(diffCoeff_[phaseIdx][compIdx]);
+ }
+ }
+
+ // porosity
+ porosity_ = problem.spatialParams().porosity(element,
+ fvGeometry,
+ scvIdx);
+ Valgrind::CheckDefined(porosity_);
+ // energy related quantities not contained in the fluid state
+
+ asImp_().updateEnergy_(primaryVariables, 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& primaryVariables,
+ const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ int scvIdx,
+ FluidState& fluidState,
+ bool isOldSol = false)
+
+ {
+ Scalar t = Implementation::temperature_(primaryVariables, 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 = primaryVariables[switchIdx];
+ else if (formulation == pgSl)
+ Sg = 1.0 - primaryVariables[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, primaryVariables[pressureIdx]);
+ if (primaryVariables[pressureIdx] + pc < 0.0)
+ DUNE_THROW(Dumux::NumericalProblem,"Capillary pressure is too low");
+ fluidState.setPressure(nPhaseIdx, primaryVariables[pressureIdx] + pc);
+ }
+ else if (formulation == pgSl) {
+ fluidState.setPressure(nPhaseIdx, primaryVariables[pressureIdx]);
+ // Here we check for (p_g - pc) in order to ensure that (p_l > 0)
+ if (primaryVariables[pressureIdx] - pc < 0.0)
+ {
+ std::cout<< "p_g: "<< primaryVariables[pressureIdx]<<" Cap_press: "<< pc << std::endl;
+ DUNE_THROW(Dumux::NumericalProblem,"Capillary pressure is too high");
+ }
+ fluidState.setPressure(wPhaseIdx, primaryVariables[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, primaryVariables[compIdx]);
+ }
+
+ Miscible2pNcComposition::solve(fluidState,
+ paramCache,
+ wPhaseIdx, //known phaseIdx
+ /*setViscosity=*/true,
+ /*setInternalEnergy=*/false);
+ }
+ else if (phasePresence == nPhaseOnly){
+
+ Dune::FieldVector<Scalar, numComponents> moleFrac;
+
+
+ moleFrac[wCompIdx] = primaryVariables[switchIdx];
+
+ for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
+ moleFrac[compIdx] = primaryVariables[compIdx];
+
+
+ Scalar sumMoleFracNotGas = 0;
+ for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
+ sumMoleFracNotGas+=moleFrac[compIdx];
+
+ sumMoleFracNotGas += moleFrac[wCompIdx];
+ moleFrac[nCompIdx] = 1 - sumMoleFracNotGas;
+
+
+ // 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] = primaryVariables[compIdx];
+ }
+ moleFrac[nCompIdx] = primaryVariables[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 rho = FluidSystem::density(fluidState, paramCache, phaseIdx);
+ Scalar mu = FluidSystem::viscosity(fluidState, paramCache, phaseIdx);
+
+ fluidState.setDensity(phaseIdx, rho);
+ fluidState.setViscosity(phaseIdx, mu);
+ }
+ }
+
+ /*!
+ * \brief Returns the phase state for the control-volume.
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+
+ /*!
+ * \brief Returns the saturation of a given phase within
+ * the control volume in \f$[-]\f$.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar saturation(int phaseIdx) const
+ { return fluidState_.saturation(phaseIdx); }
+
+ /*!
+ * \brief Returns the mass density of a given phase within the
+ * control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar density(int phaseIdx) const
+ {
+ if (phaseIdx < numPhases)
+ return fluidState_.density(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 < numPhases)
+ return fluidState_.molarDensity(phaseIdx);
+
+ else
+ DUNE_THROW(Dune::InvalidStateException, "Invalid phase index " << phaseIdx);
+ }
+
+ /*!
+ * \brief Returns the effective pressure of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar pressure(int phaseIdx) const
+ {
+ return fluidState_.pressure(phaseIdx);
+ }
+
+ /*!
+ * \brief Returns temperature inside the sub-control volume.
+ *
+ * Note that we assume thermodynamic equilibrium, i.e. the
+ * temperature of the rock matrix and of all fluid phases are
+ * identical.
+ */
+ Scalar temperature() const
+ { return fluidState_.temperature(/*phaseIdx=*/0); }
+
+ /*!
+ * \brief Returns the effective mobility of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar mobility(int phaseIdx) const
+ {
+ return mobility_[phaseIdx];
+ }
+
+ /*!
+ * \brief Returns the effective capillary pressure within the control volume
+ * in \f$[kg/(m*s^2)=N/m^2=Pa]\f$.
+ */
+ Scalar capillaryPressure() const
+ { return fluidState_.pressure(FluidSystem::nPhaseIdx) - fluidState_.pressure(FluidSystem::wPhaseIdx); }
+
+ /*!
+ * \brief Returns the average porosity within the control volume.
+ */
+ Scalar porosity() const
+ { return porosity_; }
+
+
+ /*!
+ * \brief Returns the binary diffusion coefficients for a phase in \f$[m^2/s]\f$.
+ */
+ Scalar diffCoeff(int phaseIdx, int compIdx) const
+ { return diffCoeff_[phaseIdx][compIdx]; }
+
+protected:
+
+ 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 Called by update() to compute the energy related quantities
+ */
+ void updateEnergy_(const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx,
+ bool isOldSol)
+ { };
+
+ Scalar porosity_; //!< Effective porosity within the control volume
+ Scalar mobility_[numPhases]; //!< Effective mobility within the control volume
+ Scalar density_;
+ FluidState fluidState_;
+ Scalar theta_;
+ Scalar InitialPorosity_;
+ Scalar molWtPhase_[numPhases];
+ Dune::FieldMatrix<Scalar, numPhases, numComponents> diffCoeff_;
+
+private:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+
+};
+
+} // end namespace
+
+#endif