From 3769ec5ef5b8236536d0cde8616179b49a2085a9 Mon Sep 17 00:00:00 2001
From: Holger Class <holger.class@iws.uni-stuttgart.de>
Date: Mon, 26 Aug 2013 11:27:47 +0000
Subject: [PATCH] the 3p module is new to Dumux
git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@11262 2fb0f335-1f38-0410-981e-8018bf24f1b0
---
dumux/implicit/3p/3pindices.hh | 69 ++++++
dumux/implicit/3p/3plocalresidual.hh | 207 ++++++++++++++++
dumux/implicit/3p/3pmodel.hh | 275 ++++++++++++++++++++++
dumux/implicit/3p/3pproperties.hh | 69 ++++++
dumux/implicit/3p/3ppropertydefaults.hh | 122 ++++++++++
dumux/implicit/3p/3pvolumevariables.hh | 299 ++++++++++++++++++++++++
dumux/implicit/3p/Makefile.am | 4 +
7 files changed, 1045 insertions(+)
create mode 100644 dumux/implicit/3p/3pindices.hh
create mode 100644 dumux/implicit/3p/3plocalresidual.hh
create mode 100644 dumux/implicit/3p/3pmodel.hh
create mode 100644 dumux/implicit/3p/3pproperties.hh
create mode 100644 dumux/implicit/3p/3ppropertydefaults.hh
create mode 100644 dumux/implicit/3p/3pvolumevariables.hh
create mode 100644 dumux/implicit/3p/Makefile.am
diff --git a/dumux/implicit/3p/3pindices.hh b/dumux/implicit/3p/3pindices.hh
new file mode 100644
index 0000000000..3b052b84e0
--- /dev/null
+++ b/dumux/implicit/3p/3pindices.hh
@@ -0,0 +1,69 @@
+// -*- 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 3p model.
+ */
+#ifndef DUMUX_3P_INDICES_HH
+#define DUMUX_3P_INDICES_HH
+
+#include "3pproperties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ThreePModel
+ * \ingroup ImplicitIndices
+ * \brief The indices for the isothermal 3p model.
+ *
+ * \tparam formulation The formulation, only pgSwSn
+ * \tparam PVOffset The first index in a primary variable vector.
+ */
+template <class TypeTag, int PVOffset = 0>
+class ThreePIndices
+{
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+public:
+ // Phase indices
+ static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< index of the wetting liquid phase
+ static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< index of the nonwetting liquid phase
+ static const int gPhaseIdx = FluidSystem::gPhaseIdx; //!< index of the gas phase
+
+
+ // Primary variable indices
+ static const int pressureIdx = PVOffset + 0; //!< Index for gas phase pressure in a solution vector
+ static const int swIdx = PVOffset + 1; //!< Index of water (more wetting than the other liquid) saturation
+ static const int snIdx = PVOffset + 2; //!< Index of (e.g.) NAPL saturation
+
+
+ // equation indices
+ static const int conti0EqIdx = PVOffset + wPhaseIdx; //!< Index of the mass conservation equation for the water component
+ static const int conti1EqIdx = conti0EqIdx + nPhaseIdx; //!< Index of the mass conservation equation for the contaminant component
+ static const int conti2EqIdx = conti0EqIdx + gPhaseIdx; //!< Index of the mass conservation equation for the air component
+
+ static const int contiWEqIdx = conti0EqIdx + wPhaseIdx; //!< index of the mass conservation equation for the water component
+ static const int contiNEqIdx = conti0EqIdx + nPhaseIdx; //!< index of the mass conservation equation for the contaminant component
+ static const int contiGEqIdx = conti0EqIdx + gPhaseIdx; //!< index of the mass conservation equation for the air component
+};
+
+}
+
+#endif
diff --git a/dumux/implicit/3p/3plocalresidual.hh b/dumux/implicit/3p/3plocalresidual.hh
new file mode 100644
index 0000000000..bb6d42324d
--- /dev/null
+++ b/dumux/implicit/3p/3plocalresidual.hh
@@ -0,0 +1,207 @@
+// -*- 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 three-phase fully implicit model.
+ */
+#ifndef DUMUX_3P_LOCAL_RESIDUAL_HH
+#define DUMUX_3P_LOCAL_RESIDUAL_HH
+
+#include "3pproperties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup ThreePModel
+ * \ingroup ImplicitLocalResidual
+ * \brief Element-wise calculation of the Jacobian matrix for problems
+ * using the three-phase fully implicit model.
+ *
+ * This class is used to fill the gaps in BoxLocalResidual for the 3P flow.
+ */
+template<class TypeTag>
+class ThreePLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
+{
+protected:
+ typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
+
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum {
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+
+ conti0EqIdx = Indices::conti0EqIdx,//!< Index of the mass conservation equation for the water component
+ conti1EqIdx = Indices::conti1EqIdx,//!< Index of the mass conservation equation for the contaminant component
+ conti2EqIdx = Indices::conti2EqIdx,//!< Index of the mass conservation equation for the gas component
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+ };
+
+ typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
+
+public:
+ /*!
+ * \brief Evaluate the amount all conservation quantities
+ * (e.g. phase mass) within a sub-control volume.
+ *
+ * The result should be averaged over the volume (e.g. phase mass
+ * inside a sub control volume divided by the volume)
+ *
+ * \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, const 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 components within all phases
+ storage = 0;
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ storage[conti0EqIdx + phaseIdx] +=
+ volVars.porosity()
+ * volVars.saturation(phaseIdx)
+ * volVars.density(phaseIdx);
+ }
+ }
+
+ /*!
+ * \brief Evaluates the total flux of all conservation quantities
+ * over a face of a sub-control volume.
+ *
+ * \param flux The flux over the SCV (sub-control-volume) face for each component
+ * \param faceIdx The index of the SCV face
+ * \param onBoundary A boolean variable to specify whether the flux variables
+ * are calculated for interior SCV faces or boundary faces, default=false
+ */
+ void computeFlux(PrimaryVariables &flux, const int faceIdx, const bool onBoundary=false) const
+ {
+ FluxVariables fluxVars(this->problem_(),
+ this->element_(),
+ this->fvGeometry_(),
+ faceIdx,
+ this->curVolVars_(),
+ onBoundary);
+
+ flux = 0;
+ asImp_()->computeAdvectiveFlux(flux, fluxVars);
+ }
+
+ /*!
+ * \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
+ {
+ Scalar massUpwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
+
+ ////////
+ // advective fluxes of all components in all phases
+ ////////
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ // data attached to upstream and the downstream vertices
+ // of the current phase
+ const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
+ const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
+
+ // add advective flux of current phase
+ // if alpha > 0 und alpha < 1 then both upstream and downstream
+ // nodes need their contribution
+ // if alpha == 1 (which is mostly the case) then, the downstream
+ // node is not evaluated
+ int eqIdx = conti0EqIdx + phaseIdx;
+ flux[eqIdx] += fluxVars.volumeFlux(phaseIdx)
+ * (massUpwindWeight
+ * up.fluidState().density(phaseIdx)
+ +
+ (1.0 - massUpwindWeight)
+ * dn.fluidState().density(phaseIdx));
+ }
+ }
+
+ /*!
+ * \brief Adds the diffusive flux to the flux vector over
+ * the face of a sub-control volume.
+ *
+ * \param flux The diffusive flux over the sub-control-volume face for each phase
+ * \param fluxVars The flux variables at the current SCV
+ *
+ * It doesn't do anything in three-phase model but may be used by the
+ * non-isothermal three-phase models to calculate diffusive heat
+ * fluxes
+ */
+ void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
+ {
+ // diffusive fluxes
+ flux += 0.0;
+ }
+
+ /*!
+ * \brief Calculate the source term of the equation
+ *
+ * \param source The source/sink in the SCV for each phase
+ * \param scvIdx The index of the SCV
+ */
+ void computeSource(PrimaryVariables &source, const int scvIdx)
+ {
+ this->problem_().solDependentSource(source,
+ this->element_(),
+ this->fvGeometry_(),
+ scvIdx,
+ this->curVolVars_());
+ }
+
+protected:
+ Implementation *asImp_()
+ {
+ return static_cast<Implementation *> (this);
+ }
+
+ const Implementation *asImp_() const
+ {
+ return static_cast<const Implementation *> (this);
+ }
+};
+
+} // end namepace
+
+#endif
diff --git a/dumux/implicit/3p/3pmodel.hh b/dumux/implicit/3p/3pmodel.hh
new file mode 100644
index 0000000000..6a40722ad7
--- /dev/null
+++ b/dumux/implicit/3p/3pmodel.hh
@@ -0,0 +1,275 @@
+// -*- 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 scheme to the three-phase flow model.
+ *
+ * The model is designed for simulating three fluid phases with water, gas, and
+ * a liquid contaminant (NAPL - non-aqueous phase liquid)
+ */
+#ifndef DUMUX_3P_MODEL_HH
+#define DUMUX_3P_MODEL_HH
+
+#include <dumux/implicit/common/implicitvelocityoutput.hh>
+#include "3pproperties.hh"
+
+namespace Dumux
+{
+/*!
+ * \ingroup ThreePModel
+ * \brief Adaption of the fully implicit scheme to the three-phase flow model.
+ *
+ * This model implements three-phase flow of three fluid phases
+ * \f$\alpha \in \{ water, gas, NAPL \}\f$
+ * The standard multiphase Darcy
+ * approach is used as the equation for the conservation of momentum.
+ *
+ * By inserting this into the equations for the conservation of the
+ * components, the well-known multiphase flow equation is obtained.
+ *
+ * All equations are discretized using a vertex-centered finite volume (box)
+ * or cell-centered finite volume scheme as spatial
+ * and the implicit Euler method as time discretization.
+ *
+ * The model uses commonly applied auxiliary conditions like
+ * \f$S_w + S_n + S_g = 1\f$ for the saturations.
+ * Furthermore, the phase pressures are related to each other via
+ * capillary pressures between the fluid phases, which are functions of
+ * the saturation, e.g. according to the approach of Parker et al.
+ *
+ * The used primary variables are gas phase pressure, water saturation and
+ * NAPL saturation.
+ */
+template<class TypeTag>
+class ThreePModel: public GET_PROP_TYPE(TypeTag, BaseModel)
+{
+ 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, SolutionVector) SolutionVector;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+ enum {
+ dim = GridView::dimension,
+ dimWorld = GridView::dimensionworld,
+
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+
+ pressureIdx = Indices::pressureIdx,
+ swIdx = Indices::swIdx,
+ snIdx = Indices::snIdx,
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+ };
+
+
+ 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;
+ 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);
+
+ if (isBox)
+ {
+ VertexIterator vIt = this->gridView_().template begin<dim> ();
+ const VertexIterator &vEndIt = this->gridView_().template end<dim> ();
+ for (; vIt != vEndIt; ++vIt)
+ {
+ int globalIdx = this->dofMapper().map(*vIt);
+ const GlobalPosition &globalPos = vIt->geometry().corner(0);
+ }
+ }
+ else
+ {
+ ElementIterator eIt = this->gridView_().template begin<0> ();
+ const ElementIterator &eEndIt = this->gridView_().template end<0> ();
+ for (; eIt != eEndIt; ++eIt)
+ {
+ int globalIdx = this->dofMapper().map(*eIt);
+ const GlobalPosition &globalPos = eIt->geometry().center();
+ }
+ }
+ }
+
+ /*!
+ * \brief Called by the update() method if applying the newton
+ * method was unsuccessful.
+ */
+ void updateFailed()
+ {
+ ParentType::updateFailed();
+ };
+
+ /*!
+ * \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();
+ }
+
+ /*!
+ * \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<double, 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 *saturation[numPhases];
+ ScalarField *pressure[numPhases];
+ ScalarField *density[numPhases];
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ saturation[phaseIdx] = writer.allocateManagedBuffer(numDofs);
+ pressure[phaseIdx] = writer.allocateManagedBuffer(numDofs);
+ density[phaseIdx] = writer.allocateManagedBuffer(numDofs);
+ }
+
+ ScalarField *temperature = writer.allocateManagedBuffer (numDofs);
+ ScalarField *poro = writer.allocateManagedBuffer(numDofs);
+ ScalarField *perm = writer.allocateManagedBuffer(numDofs);
+ VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
+ VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
+ VectorField *velocityG = 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);
+ (*velocityG)[i] = Scalar(0);
+ }
+ }
+
+ unsigned numElements = this->gridView_().size(0);
+ ScalarField *rank = writer.allocateManagedBuffer (numElements);
+
+ 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();
+
+ FVElementGeometry fvGeometry;
+ fvGeometry.update(this->gridView_(), *eIt);
+
+
+ ElementVolumeVariables elemVolVars;
+ elemVolVars.update(this->problem_(),
+ *eIt,
+ fvGeometry,
+ false /* oldSol? */);
+
+ for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+ {
+ int globalIdx = this->dofMapper().map(*eIt, scvIdx, dofCodim);
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
+ (*saturation[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().saturation(phaseIdx);
+ (*pressure[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().pressure(phaseIdx);
+ (*density[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().density(phaseIdx);
+ }
+
+ (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
+ (*perm)[globalIdx] = elemVolVars[scvIdx].permeability();
+ (*temperature)[globalIdx] = elemVolVars[scvIdx].temperature();
+ }
+
+ // velocity output
+ velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, *eIt, wPhaseIdx);
+ velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, nPhaseIdx);
+ velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, *eIt, gPhaseIdx);
+
+ }
+
+ writer.attachDofData(*saturation[wPhaseIdx], "sw", isBox);
+ writer.attachDofData(*saturation[nPhaseIdx], "sn", isBox);
+ writer.attachDofData(*saturation[gPhaseIdx], "sg", isBox);
+ writer.attachDofData(*pressure[wPhaseIdx], "pw", isBox);
+ writer.attachDofData(*pressure[nPhaseIdx], "pn", isBox);
+ writer.attachDofData(*pressure[gPhaseIdx], "pg", isBox);
+ writer.attachDofData(*density[wPhaseIdx], "rhow", isBox);
+ writer.attachDofData(*density[nPhaseIdx], "rhon", isBox);
+ writer.attachDofData(*density[gPhaseIdx], "rhog", isBox);
+
+ writer.attachDofData(*poro, "porosity", isBox);
+ writer.attachDofData(*perm, "permeability", isBox);
+ writer.attachDofData(*temperature, "temperature", isBox);
+
+ if (velocityOutput.enableOutput()) // check if velocity output is demanded
+ {
+ writer.attachDofData(*velocityW, "velocityW", isBox, dim);
+ writer.attachDofData(*velocityN, "velocityN", isBox, dim);
+ writer.attachDofData(*velocityG, "velocityG", isBox, dim);
+ }
+
+ writer.attachCellData(*rank, "process rank");
+ }
+
+};
+
+}
+
+#include "3ppropertydefaults.hh"
+
+#endif
diff --git a/dumux/implicit/3p/3pproperties.hh b/dumux/implicit/3p/3pproperties.hh
new file mode 100644
index 0000000000..f6b3949618
--- /dev/null
+++ b/dumux/implicit/3p/3pproperties.hh
@@ -0,0 +1,69 @@
+// -*- 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 ThreePModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines the properties required for the 3p fully implicit model.
+ */
+#ifndef DUMUX_3P_PROPERTIES_HH
+#define DUMUX_3P_PROPERTIES_HH
+
+#include <dumux/implicit/box/boxproperties.hh>
+#include <dumux/implicit/cellcentered/ccproperties.hh>
+
+namespace Dumux
+{
+
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! The type tags for the implicit three-phase problems
+NEW_TYPE_TAG(ThreeP);
+NEW_TYPE_TAG(BoxThreeP, INHERITS_FROM(BoxModel, ThreeP));
+NEW_TYPE_TAG(CCThreeP, INHERITS_FROM(CCModel, ThreeP));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
+NEW_PROP_TAG(Indices); //!< Enumerations for 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(FluidState); //!< the phases' state
+
+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(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
+NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the upwind parameter for the mobility
+NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
+NEW_PROP_TAG(BaseFluxVariables); //! The base flux variables
+NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
+NEW_PROP_TAG(VtkAddVelocity); //!< Returns whether velocity vectors are written into the vtk output
+}
+}
+
+#endif
diff --git a/dumux/implicit/3p/3ppropertydefaults.hh b/dumux/implicit/3p/3ppropertydefaults.hh
new file mode 100644
index 0000000000..cd31aa893b
--- /dev/null
+++ b/dumux/implicit/3p/3ppropertydefaults.hh
@@ -0,0 +1,122 @@
+// -*- 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 ThreePModel
+ */
+/*!
+ * \file
+ *
+ * \brief Defines default values for most properties required by the
+ * 3p fully implicit model.
+ */
+#ifndef DUMUX_3P_PROPERTY_DEFAULTS_HH
+#define DUMUX_3P_PROPERTY_DEFAULTS_HH
+
+#include "3pindices.hh"
+
+#include "3pmodel.hh"
+#include "3pindices.hh"
+#include "3pvolumevariables.hh"
+#include "3pproperties.hh"
+#include "3plocalresidual.hh"
+
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
+#include <dumux/material/spatialparams/implicitspatialparams.hh>
+#include <dumux/material/fluidstates/immisciblefluidstate.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
+
+
+namespace Dumux
+{
+
+namespace Properties {
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+
+/*!
+ * \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(ThreeP, NumPhases)
+{
+ private:
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+ public:
+ static const int value = FluidSystem::numPhases;
+ static_assert(value == 3,
+ "Only fluid systems with 3 phases are supported by the 3p model!");
+};
+
+SET_INT_PROP(ThreeP, NumEq, 3); //!< set the number of equations to 2
+
+/*!
+ * \brief Set the property for the material parameters by extracting
+ * it from the material law.
+ */
+SET_TYPE_PROP(ThreeP, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
+
+//! The local residual function of the conservation equations
+SET_TYPE_PROP(ThreeP, LocalResidual, ThreePLocalResidual<TypeTag>);
+
+//! the Model property
+SET_TYPE_PROP(ThreeP, Model, ThreePModel<TypeTag>);
+
+//! the VolumeVariables property
+SET_TYPE_PROP(ThreeP, VolumeVariables, ThreePVolumeVariables<TypeTag>);
+
+//! the FluxVariables property
+SET_TYPE_PROP(ThreeP, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
+
+//! the upwind factor for the mobility.
+SET_SCALAR_PROP(ThreeP, ImplicitMassUpwindWeight, 1.0);
+
+//! set default mobility upwind weight to 1.0, i.e. fully upwind
+SET_SCALAR_PROP(ThreeP, ImplicitMobilityUpwindWeight, 1.0);
+
+//! The indices required by the isothermal 3p model
+SET_TYPE_PROP(ThreeP, Indices, ThreePIndices<TypeTag, /*PVOffset=*/0>);
+
+//! The spatial parameters to be employed.
+//! Use ImplicitSpatialParams by default.
+SET_TYPE_PROP(ThreeP, SpatialParams, ImplicitSpatialParams<TypeTag>);
+
+SET_PROP(ThreeP, FluidState)
+{
+private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+public:
+ typedef ImmiscibleFluidState<Scalar, FluidSystem> type;
+};
+
+// disable velocity output by default
+SET_BOOL_PROP(ThreeP, VtkAddVelocity, false);
+
+// enable gravity by default
+SET_BOOL_PROP(ThreeP, ProblemEnableGravity, true);
+
+}
+
+}
+
+#endif
diff --git a/dumux/implicit/3p/3pvolumevariables.hh b/dumux/implicit/3p/3pvolumevariables.hh
new file mode 100644
index 0000000000..6f4f0a7ba8
--- /dev/null
+++ b/dumux/implicit/3p/3pvolumevariables.hh
@@ -0,0 +1,299 @@
+// -*- 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 three-phase model.
+ */
+#ifndef DUMUX_3P_VOLUME_VARIABLES_HH
+#define DUMUX_3P_VOLUME_VARIABLES_HH
+
+#include <dumux/implicit/common/implicitmodel.hh>
+#include <dumux/common/math.hh>
+
+#include <dune/common/collectivecommunication.hh>
+#include <vector>
+#include <iostream>
+
+#include "3pproperties.hh"
+
+#include <dumux/material/constants.hh>
+#include <dumux/material/fluidstates/immisciblefluidstate.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ThreePModel
+ * \brief Contains the quantities which are are constant within a
+ * finite volume in the two-phase, two-component model.
+ */
+template <class TypeTag>
+class ThreePVolumeVariables : 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, 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,
+
+ numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+
+ swIdx = Indices::swIdx,
+ snIdx = Indices::snIdx,
+ pressureIdx = Indices::pressureIdx
+ };
+
+ static const Scalar R; // universal gas constant
+
+ typedef typename GridView::template Codim<0>::Entity Element;
+
+ enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
+
+public:
+ //! The type of the object returned by the fluidState() method
+ typedef Dumux::ImmiscibleFluidState<Scalar, FluidSystem> FluidState;
+
+
+ /*!
+ * \copydoc ImplicitVolumeVariables::update
+ */
+ void update(const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx,
+ bool isOldSol)
+ {
+ ParentType::update(priVars,
+ problem,
+ element,
+ fvGeometry,
+ scvIdx,
+ isOldSol);
+
+ // capillary pressure parameters
+ const MaterialLawParams &materialParams =
+ problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
+
+ int globalIdx = problem.model().dofMapper().map(element, scvIdx, dofCodim);
+
+ Scalar temp = Implementation::temperature_(priVars, problem, element, fvGeometry, scvIdx);
+ fluidState_.setTemperature(temp);
+
+ sw_ = priVars[swIdx];
+ sn_ = priVars[snIdx];
+ sg_ = 1. - sw_ - sn_;
+
+ Valgrind::CheckDefined(sg_);
+
+ fluidState_.setSaturation(wPhaseIdx, sw_);
+ fluidState_.setSaturation(gPhaseIdx, sg_);
+ fluidState_.setSaturation(nPhaseIdx, sn_);
+
+ /* now the pressures */
+ pg_ = priVars[pressureIdx];
+
+ // calculate capillary pressures
+ Scalar pcgw = MaterialLaw::pcgw(materialParams, sw_);
+ Scalar pcnw = MaterialLaw::pcnw(materialParams, sw_);
+ Scalar pcgn = MaterialLaw::pcgn(materialParams, sw_ + sn_);
+
+ Scalar pcAlpha = MaterialLaw::pcAlpha(materialParams, sn_);
+ Scalar pcNW1 = 0.0; // TODO: this should be possible to assign in the problem file
+
+ pn_ = pg_- pcAlpha * pcgn - (1.-pcAlpha)*(pcgw - pcNW1);
+ pw_ = pn_ - pcAlpha * pcnw - (1.-pcAlpha)*pcNW1;
+
+ fluidState_.setPressure(wPhaseIdx, pw_);
+ fluidState_.setPressure(gPhaseIdx, pg_);
+ fluidState_.setPressure(nPhaseIdx, pn_);
+
+ Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
+ Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
+ Scalar rhoN = FluidSystem::density(fluidState_, nPhaseIdx);
+
+ fluidState_.setDensity(wPhaseIdx, rhoW);
+ fluidState_.setDensity(gPhaseIdx, rhoG);
+ fluidState_.setDensity(nPhaseIdx, rhoN);
+
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+ // Mobilities
+ const Scalar mu =
+ FluidSystem::viscosity(fluidState_,
+ phaseIdx);
+ fluidState_.setViscosity(phaseIdx,mu);
+
+ Scalar kr;
+ kr = MaterialLaw::kr(materialParams, phaseIdx,
+ fluidState_.saturation(wPhaseIdx),
+ fluidState_.saturation(nPhaseIdx),
+ fluidState_.saturation(gPhaseIdx));
+ mobility_[phaseIdx] = kr / mu;
+ Valgrind::CheckDefined(mobility_[phaseIdx]);
+ }
+
+ // porosity
+ porosity_ = problem.spatialParams().porosity(element,
+ fvGeometry,
+ scvIdx);
+ Valgrind::CheckDefined(porosity_);
+
+ // permeability
+ permeability_ = problem.spatialParams().intrinsicPermeability(element,
+ fvGeometry,
+ scvIdx);
+ Valgrind::CheckDefined(permeability_);
+
+ // energy related quantities not contained in the fluid state
+ asImp_().updateEnergy_(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
+ }
+
+ /*!
+ * \brief Returns the phase state for the control-volume.
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+
+ /*!
+ * \brief Returns the effective saturation of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar saturation(const 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(const int phaseIdx) const
+ { return fluidState_.density(phaseIdx); }
+
+ /*!
+ * \brief Returns the effective pressure of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar pressure(const int phaseIdx) const
+ { return fluidState_.pressure(phaseIdx); }
+
+ /*!
+ * \brief Returns temperature inside the sub-control volume.
+ *
+ * Note that we assume thermodynamic equilibrium, i.e. the
+ * temperatures of the rock matrix and of all fluid phases are
+ * identical.
+ */
+ Scalar temperature() const
+ { return fluidState_.temperature(/*phaseIdx=*/0); }
+
+ /*!
+ * \brief Returns the effective mobility of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The phase index
+ */
+ Scalar mobility(const int phaseIdx) const
+ {
+ return mobility_[phaseIdx];
+ }
+
+ /*!
+ * \brief Returns the effective capillary pressure within the control volume.
+ */
+ Scalar capillaryPressure() const
+ { return fluidState_.capillaryPressure(); }
+
+ /*!
+ * \brief Returns the average porosity within the control volume.
+ */
+ Scalar porosity() const
+ { return porosity_; }
+
+ /*!
+ * \brief Returns the permeability within the control volume.
+ */
+ Scalar permeability() const
+ { return permeability_; }
+
+protected:
+
+ static Scalar temperature_(const PrimaryVariables &priVars,
+ const Problem &problem,
+ const Element &element,
+ const FVElementGeometry &fvGeometry,
+ const int scvIdx)
+ {
+ return problem.temperatureAtPos(fvGeometry.subContVol[scvIdx].global);
+ }
+
+ /*!
+ * \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 sw_, sg_, sn_, pg_, pw_, pn_;
+
+ Scalar porosity_; //!< Effective porosity within the control volume
+ Scalar permeability_; //!< Effective porosity within the control volume
+ Scalar mobility_[numPhases]; //!< Effective mobility within the control volume
+ Scalar bulkDensTimesAdsorpCoeff_; //!< the basis for calculating adsorbed NAPL
+ FluidState fluidState_;
+
+private:
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+};
+
+template <class TypeTag>
+const typename ThreePVolumeVariables<TypeTag>::Scalar ThreePVolumeVariables<TypeTag>::R = Constants<typename GET_PROP_TYPE(TypeTag, Scalar)>::R;
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/3p/Makefile.am b/dumux/implicit/3p/Makefile.am
new file mode 100644
index 0000000000..f50519766e
--- /dev/null
+++ b/dumux/implicit/3p/Makefile.am
@@ -0,0 +1,4 @@
+3pdir = $(includedir)/dumux/implicit/3p
+3p_HEADERS = *.hh
+
+include $(top_srcdir)/am/global-rules
--
GitLab