diff --git a/dumux/porousmediumflow/2p1c/implicit/indices.hh b/dumux/porousmediumflow/2p1c/implicit/indices.hh
index 49a07a06e8fea6035c1e9e8d98799e83212b734b..6bd3e53b970df39299c7cdc9e88174b6d1474dd3 100644
--- a/dumux/porousmediumflow/2p1c/implicit/indices.hh
+++ b/dumux/porousmediumflow/2p1c/implicit/indices.hh
@@ -25,11 +25,22 @@
#ifndef DUMUX_2P1C_INDICES_HH
#define DUMUX_2P1C_INDICES_HH
-#include "properties.hh"
-
namespace Dumux
{
+/*!
+ * \ingroup TwoPNCModel
+ * \ingroup ImplicitIndices
+ * \brief Enumerates the formulations which the two-phase n-component model accepts.
+ */
+struct TwoPOneCFormulation
+{
+ enum {
+ pnsw,
+ pwsn
+ };
+};
+
/*!
* \ingroup TwoPOneCModel
* \ingroup ImplicitIndices
@@ -46,12 +57,12 @@ class TwoPOneCIndices
public:
// Phase indices
static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< index of the wetting liquid phase
- static const int gPhaseIdx = FluidSystem::nPhaseIdx; //!< index of the gas phase
+ static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< index of the gas phase
// present phases (-> 'pseudo' primary variable)
static const int twoPhases = 1; //!< All three phases are present
static const int wPhaseOnly = 2; //!< Only the water phase is present
- static const int gPhaseOnly = 3; //!< Only gas phase is present
+ static const int nPhaseOnly = 3; //!< Only gas phase is present
// Primary variable indices
static const int pressureIdx = PVOffset + 0; //!< Index for phase pressure in a solution vector
diff --git a/dumux/porousmediumflow/2p1c/implicit/localresidual.hh b/dumux/porousmediumflow/2p1c/implicit/localresidual.hh
deleted file mode 100644
index 51a772df81669f7b9b2821a65cf12b2256e535fe..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2p1c/implicit/localresidual.hh
+++ /dev/null
@@ -1,272 +0,0 @@
-// -*- 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 one-component fully implicit model.
- *
- * Important note: The 2p1c model requires the use of the non-isothermal extension found in dumux/porousmediumflow/nonisothermal/implicit/
- */
-#ifndef DUMUX_2P1C_LOCAL_RESIDUAL_HH
-#define DUMUX_2P1C_LOCAL_RESIDUAL_HH
-
-#include "properties.hh"
-
-namespace Dumux
-{
-/*!
- * \ingroup TwoPOneCModel
- * \ingroup ImplicitLocalResidual
- * \brief Element-wise calculation of the Jacobian matrix for problems
- * using the two-phase one-component fully implicit model.
- *
- * This class depends on the non-isothermal model.
- *
- */
-template<class TypeTag>
-class TwoPOneCLocalResidual: 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;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes) ElementBoundaryTypes;
-
-
- enum {
- numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
- numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
-
- conti0EqIdx = Indices::conti0EqIdx, //!< Index of the mass conservation equation for the water component
- energyEqIdx = Indices::energyEqIdx, //!< Index of the energy conservation equation
- wPhaseIdx = Indices::wPhaseIdx,
- 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;
-
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- static const bool useBlockingOfSpuriousFlow = GET_PROP_VALUE(TypeTag, UseBlockingOfSpuriousFlow);
-
-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
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- storage[conti0EqIdx] +=
- 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;
- fluxVars.update(this->problem_(),
- this->element_(),
- this->fvGeometry_(),
- faceIdx,
- this->curVolVars_(),
- onBoundary);
-
- flux = 0;
- computeAdvectiveFlux(flux, fluxVars);
- asImp_()->computeDiffusiveFlux(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);
-
- // loop over 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));
-
- Scalar factor = 1.0;
-
- // factor for the blocking of spurious cold water fluxes into the steam zone (Gudbjerg, 2004)
- if(useBlockingOfSpuriousFlow)
- { factor = factor_(up, dn, phaseIdx); }
-
-
- // advective mass flux
- flux[conti0EqIdx] += fluxVars.volumeFlux(phaseIdx)
- * (massUpwindWeight
- * up.fluidState().density(phaseIdx)
- + (1.0 - massUpwindWeight)
- * dn.fluidState().density(phaseIdx))
- * factor ;
-
- // advective heat flux
- flux[energyEqIdx] += fluxVars.volumeFlux(phaseIdx)
- * (massUpwindWeight
- * up.density(phaseIdx)
- * up.enthalpy(phaseIdx)
- + (1.0 - massUpwindWeight)
- * dn.density(phaseIdx)
- * dn.enthalpy(phaseIdx))
- * factor ;
- }
- }
-
- /*!
- * \brief Adds the diffusive mass flux of all components over
- * a face of a subcontrol volume.
- *
- * \param flux The diffusive flux over the sub-control-volume face for each component
- * \param fluxVars The flux variables at the current SCV
- */
- void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
- {}
-
- /*!
- * \brief Returns true if a spurious flow has been detected
- *
- */
- bool spuriousFlowDetected() const
- {
- return spuriousFlowDetected_;
- }
-
- /*!
- * \brief Used to reset the respective flag
- *
- */
- void resetSpuriousFlowDetected()
- {
- spuriousFlowDetected_ = false;
- }
-
-protected:
-
- /*!
- * \brief Calculate the blocking factor which prevents spurious cold water fluxes into the steam zone (Gudbjerg et al., 2005) \cite gudbjerg2004 <BR>
- *
- * \param up The upstream volume variables
- * \param dn The downstream volume variables
- * \param phaseIdx The index of the fluid phase
- */
- Scalar factor_(const VolumeVariables &up, const VolumeVariables &dn,const int phaseIdx) const {
-
- Scalar factor = 1.0;
-
- Scalar tDn = dn.temperature(); //temperature of the downstream SCV (where the cold water is potentially intruding into a steam zone)
- Scalar tUp = up.temperature(); //temperature of the upstream SCV
-
- Scalar sgDn = dn.saturation(gPhaseIdx); //gas phase saturation of the downstream SCV
- Scalar sgUp = up.saturation(gPhaseIdx); //gas phase saturation of the upstream SCV
-
- bool upIsNotSteam = false;
- bool downIsSteam = false;
- bool spuriousFlow = false;
-
- if(sgUp <= 1e-5)
- {upIsNotSteam = true;}
-
- if(sgDn > 1e-5)
- {downIsSteam = true;}
-
- if(upIsNotSteam && downIsSteam && tDn > tUp && phaseIdx == wPhaseIdx)
- {
- spuriousFlow = true;
- spuriousFlowDetected_ = true;
- }
-
- if(spuriousFlow)
- {
- Scalar deltaT = tDn - tUp;
-
- if((deltaT) > 15 )
- {factor = 0.0 ;}
-
- else
- { factor = 1-(deltaT/15); }
-
- }
- return factor;
- }
-
- Implementation *asImp_()
- {
- return static_cast<Implementation *> (this);
- }
-
- const Implementation *asImp_() const
- {
- return static_cast<const Implementation *> (this);
- }
-
-private:
- mutable bool spuriousFlowDetected_;// true if spurious cold-water flow into a steam zone is happening
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/porousmediumflow/2p1c/implicit/model.hh b/dumux/porousmediumflow/2p1c/implicit/model.hh
index 0ec78078c25bc6ad6d20375d57715ea7fcff2146..b66da7eaf00d865092e572b9a06b007799c94c79 100644
--- a/dumux/porousmediumflow/2p1c/implicit/model.hh
+++ b/dumux/porousmediumflow/2p1c/implicit/model.hh
@@ -25,8 +25,24 @@
#ifndef DUMUX_2P1C_MODEL_HH
#define DUMUX_2P1C_MODEL_HH
-#include <dumux/porousmediumflow/implicit/velocityoutput.hh>
-#include "properties.hh"
+#include <dumux/common/properties.hh>
+
+#include <dumux/material/spatialparams/implicit.hh>
+#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
+#include <dumux/material/fluidstates/immiscible.hh>
+#include <dumux/material/fluidstates/compositional.hh>
+
+#include <dumux/porousmediumflow/properties.hh>
+#include <dumux/porousmediumflow/immiscible/localresidual.hh>
+#include <dumux/porousmediumflow/compositional/switchableprimaryvariables.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/properties.hh>
+
+#include <dumux/porousmediumflow/2p/implicit/vtkoutputfields.hh>
+
+#include "indices.hh"
+#include "volumevariables.hh"
+#include "primaryvariableswitch.hh"
namespace Dumux
{
@@ -66,559 +82,114 @@ namespace Dumux
* If only one phase is present, \f$p_g\f$ and \f$T\f$ are the primary variables.
* In the presence of two phases, \f$p_g\f$ and \f$S_w\f$ become the new primary variables.
*/
-template<class TypeTag>
-class TwoPOneCModel: 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),
- numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
-
- pressureIdx = Indices::pressureIdx,
- switch1Idx = Indices::switch1Idx,
-
- wPhaseIdx = Indices::wPhaseIdx,
- gPhaseIdx = Indices::gPhaseIdx,
-
- twoPhases = Indices::twoPhases,
- wPhaseOnly = Indices::wPhaseOnly,
- gPhaseOnly = Indices::gPhaseOnly,
- };
-
- 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::FieldMatrix<Scalar,dim,dim> FieldMatrix;
-
- 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);
-
- staticDat_.resize(this->numDofs());
-
- setSwitched_(false);
-
- if (isBox)
- {
- for(const auto& vertex: vertices(this->gridView_()))
- {
- int globalIdx = this->dofMapper().index(vertex);
- const GlobalPosition &globalPos = vertex.geometry().corner(0);
-
- // initialize phase presence
- staticDat_[globalIdx].phasePresence
- = this->problem_().initialPhasePresence(vertex, globalIdx,
- globalPos);
- staticDat_[globalIdx].wasSwitched = false;
-
- staticDat_[globalIdx].oldPhasePresence
- = staticDat_[globalIdx].phasePresence;
- }
- }
- else
- {
- for (const auto& element : elements(this->gridView_()))
- {
- int globalIdx = this->dofMapper().index(element);
- const GlobalPosition &globalPos = element.geometry().center();
-
- // initialize phase presence
- staticDat_[globalIdx].phasePresence
- = this->problem_().initialPhasePresence(*this->gridView_().template begin<dim> (),
- globalIdx, globalPos);
- staticDat_[globalIdx].wasSwitched = false;
-
- staticDat_[globalIdx].oldPhasePresence
- = staticDat_[globalIdx].phasePresence;
- }
- }
- }
-
- /*!
- * \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 degree of freedom.
- *
- * \param globalIdx The global index of the degree of freedom
- * \param oldSol Evaluate function with solution of current or previous time step
- */
- int phasePresence(int globalIdx, bool oldSol) const
- {
- return
- oldSol
- ? staticDat_[globalIdx].oldPhasePresence
- : staticDat_[globalIdx].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<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];
-
- ScalarField *mobility[numPhases];
- ScalarField *viscosity[numPhases];
-
-
- for (int phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
- saturation[phaseIdx] = writer.allocateManagedBuffer(numDofs);
- pressure[phaseIdx] = writer.allocateManagedBuffer(numDofs);
- density[phaseIdx] = writer.allocateManagedBuffer(numDofs);
-
- mobility[phaseIdx] = writer.allocateManagedBuffer(numDofs);
- viscosity[phaseIdx] = writer.allocateManagedBuffer(numDofs);
- }
-
- ScalarField *phasePresence = writer.allocateManagedBuffer (numDofs);
- ScalarField *poro = writer.allocateManagedBuffer(numDofs);
- ScalarField *permXX = writer.allocateManagedBuffer(numDofs);
- ScalarField *permYY = writer.allocateManagedBuffer(numDofs);
- ScalarField *permZZ = writer.allocateManagedBuffer(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)
- {
- (*velocityW)[i] = Scalar(0);
- (*velocityG)[i] = Scalar(0);
- }
- }
-
- unsigned numElements = this->gridView_().size(0);
- ScalarField *rank = writer.allocateManagedBuffer (numElements);
-
- for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
- {
- int eIdx = this->elementMapper().index(element);
- (*rank)[eIdx] = this->gridView_().comm().rank();
-
- FVElementGeometry fvGeometry;
- fvGeometry.update(this->gridView_(), element);
-
-
- ElementVolumeVariables elemVolVars;
- elemVolVars.update(this->problem_(),
- element,
- fvGeometry,
- false /* oldSol? */);
-
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
- {
- int globalIdx = this->dofMapper().subIndex(element, 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);
-
- (*mobility[phaseIdx])[globalIdx] = elemVolVars[scvIdx].mobility(phaseIdx);
- (*viscosity[phaseIdx])[globalIdx] = elemVolVars[scvIdx].fluidState().viscosity(phaseIdx);
- }
-
- (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
- (*phasePresence)[globalIdx] = staticDat_[globalIdx].phasePresence;
-
- FieldMatrix K = this->problem_().spatialParams().intrinsicPermeability(element,
- fvGeometry,
- scvIdx);
- (*permXX)[globalIdx] = K[0][0];
- if(dimWorld > 1)
- (*permYY)[globalIdx] = K[1][1];
- if(dimWorld > 2)
- (*permZZ)[globalIdx] = K[2][2];
- }
-
- // velocity output
- velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
- velocityOutput.calculateVelocity(*velocityG, elemVolVars, fvGeometry, element, gPhaseIdx);
- }
-
- writer.attachDofData(*saturation[wPhaseIdx], "sw", isBox);
- writer.attachDofData(*saturation[gPhaseIdx], "sg", isBox);
- writer.attachDofData(*pressure[wPhaseIdx], "pw", isBox);
- writer.attachDofData(*pressure[gPhaseIdx], "pg", isBox);
- writer.attachDofData(*density[wPhaseIdx], "rhow", isBox);
- writer.attachDofData(*density[gPhaseIdx], "rhog", isBox);
- writer.attachDofData(*mobility[wPhaseIdx], "MobW", isBox);
- writer.attachDofData(*mobility[gPhaseIdx], "MobG", isBox);
-
- writer.attachDofData(*viscosity[wPhaseIdx], "ViscosW", isBox);
- writer.attachDofData(*viscosity[gPhaseIdx], "ViscosG", isBox);
-
- writer.attachDofData(*poro, "porosity", isBox);
- writer.attachDofData(*permXX, "permeabilityXX", isBox);
- if(dimWorld > 1)
- writer.attachDofData(*permYY, "permeabilityYY", isBox);
- if(dimWorld > 2)
- writer.attachDofData(*permZZ, "permeabilityZZ", isBox);
- writer.attachDofData(*phasePresence, "phase presence", isBox);
-
- if (velocityOutput.enableOutput()) // check if velocity output is demanded
- {
- writer.attachDofData(*velocityW, "velocityW", isBox, dim);
- writer.attachDofData(*velocityG, "velocityG", 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, either a vertex or an element
- */
- template<class Entity>
- void serializeEntity(std::ostream &outStream, const Entity &entity)
- {
- // write primary variables
- ParentType::serializeEntity(outStream, entity);
-
- int dofIdxGlobal = this->dofMapper().index(entity);
-
- if (!outStream.good())
- DUNE_THROW(Dune::IOError, "Could not serialize entity " << dofIdxGlobal);
-
- outStream << staticDat_[dofIdxGlobal].phasePresence << " ";
- }
-
- /*!
- * \brief Reads the current solution from a restart file.
- *
- * \param inStream The input stream of one vertex from the restart file
- * \param entity The entity, either a vertex or an element
- */
- template<class Entity>
- void deserializeEntity(std::istream &inStream, const Entity &entity)
- {
- // read primary variables
- ParentType::deserializeEntity(inStream, entity);
-
- // read phase presence
- int dofIdxGlobal = this->dofMapper().index(entity);
-
- if (!inStream.good())
- DUNE_THROW(Dune::IOError,
- "Could not deserialize entity " << dofIdxGlobal);
-
- inStream >> staticDat_[dofIdxGlobal].phasePresence;
- staticDat_[dofIdxGlobal].oldPhasePresence
- = staticDat_[dofIdxGlobal].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;
- int succeeded;
- try {
- for (unsigned i = 0; i < staticDat_.size(); ++i)
- staticDat_[i].visited = false;
-
- FVElementGeometry fvGeometry;
- static VolumeVariables volVars;
- for (const auto& element : elements(this->gridView_()))
- {
- fvGeometry.update(this->gridView_(), element);
- for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
- {
- int globalIdx = this->dofMapper().subIndex(element, scvIdx, dofCodim);
-
- if (staticDat_[globalIdx].visited)
- continue;
-
- staticDat_[globalIdx].visited = true;
- volVars.update(curGlobalSol[globalIdx],
- this->problem_(),
- element,
- fvGeometry,
- scvIdx,
- false);
- const GlobalPosition &globalPos = fvGeometry.subContVol[scvIdx].global;
-
- if (primaryVarSwitch_(curGlobalSol,
- volVars,
- globalIdx,
- globalPos))
- {
- this->jacobianAssembler().markDofRed(globalIdx);
- wasSwitched = true;
- }
- }
- }
- succeeded = 1;
- }
-
- catch (Dumux::NumericalProblem &e)
- {
- std::cout << "\n"
- << "Rank " << this->problem_().gridView().comm().rank()
- << " caught an exception while updating the static data." << e.what()
- << "\n";
- succeeded = 0;
- }
-
- //make sure that all processes succeeded. If not throw a NumericalProblem to decrease the time step size.
- if (this->gridView_().comm().size() > 1)
- succeeded = this->gridView_().comm().min(succeeded);
-
- if (!succeeded) {
- updateFailed();
-
- DUNE_THROW(NumericalProblem,
- "A process did not succeed in updating the static data.");
-
- return;
- }
-
- // 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;
- };
-
- /*!
- * \brief Reset the current phase presence of all vertices to the old one.
- *
- * This is done after an update failed.
- */
- void resetPhasePresence_()
- {
- for (unsigned int i = 0; i < this->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_()
- {
- for (unsigned int i = 0; i < this->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;
-
- globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(gPhaseIdx);
-
- // check if a primary var switch is necessary
- if (phasePresence == twoPhases)
- {
- Scalar Smin = 0;
- if (staticDat_[globalIdx].wasSwitched)
- Smin = -0.01;
-
- if (volVars.saturation(gPhaseIdx) <= Smin)
- {
- wouldSwitch = true;
- // gas phase disappears
- std::cout << "Gas phase disappears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", sg: "
- << volVars.saturation(gPhaseIdx) << std::endl;
- newPhasePresence = wPhaseOnly;
-
- globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
- }
- else if (volVars.saturation(wPhaseIdx) <= Smin)
- {
- wouldSwitch = true;
- // water phase disappears
- std::cout << "Wetting phase disappears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", sw: "
- << volVars.saturation(wPhaseIdx) << std::endl;
- newPhasePresence = gPhaseOnly;
-
- globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
- }
-
- }
- else if (phasePresence == wPhaseOnly)
- {
- Scalar temp = volVars.fluidState().temperature();
- Scalar tempVap = volVars.vaporTemperature();
-
- // if the the temperature would be larger than
- // the vapor temperature at the given pressure, gas phase appears
- if (temp >= tempVap)
- {
- wouldSwitch = true;
- // gas phase appears
- std::cout << "gas phase appears at vertex " << globalIdx
- << ", coordinates: " << globalPos << std::endl;
-
- newPhasePresence = twoPhases;
- globalSol[globalIdx][switch1Idx] = 0.9999; //wetting phase saturation
- }
- }
-
-
- else if (phasePresence == gPhaseOnly)
- {
-
- Scalar temp = volVars.fluidState().temperature();
- Scalar tempVap = volVars.vaporTemperature();
-
- if (temp < tempVap)
- {
- wouldSwitch = true;
- // wetting phase appears
- std::cout << "wetting phase appears at vertex " << globalIdx
- << ", coordinates: " << globalPos << std::endl;
-
-
- newPhasePresence = twoPhases;
- globalSol[globalIdx][switch1Idx] = 0.0001; //arbitrary small value
- }
- }
- staticDat_[globalIdx].phasePresence = newPhasePresence;
- staticDat_[globalIdx].wasSwitched = wouldSwitch;
- return phasePresence != newPhasePresence;
- }
-
- // parameters given in constructor
- std::vector<StaticVars> staticDat_;
- bool switchFlag_;
-};
-
-}
-
-#include "propertydefaults.hh"
-
-#endif
+namespace Properties
+{
+ //////////////////////////////////////////////////////////////////
+ // Type tags
+ //////////////////////////////////////////////////////////////////
+ NEW_TYPE_TAG(TwoPOneCNI, INHERITS_FROM(PorousMediumFlow, NonIsothermal));
+
+ //////////////////////////////////////////////////////////////////
+ // Property tags
+ //////////////////////////////////////////////////////////////////
+
+ NEW_PROP_TAG(UseBlockingOfSpuriousFlow); //!< Determines whether Blocking ofspurious flow is used
+
+ /*!
+ * \brief Set the property for the number of components.
+ *
+ * We just forward the number from the fluid system and use an static
+ * assert to make sure it is 1.
+ */
+ SET_PROP(TwoPOneCNI, NumComponents)
+ {
+ private:
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ public:
+ static constexpr auto value = FluidSystem::numComponents;
+
+ static_assert(value == 1,
+ "Only fluid systems with 1 component are supported by the 2p1cni model!");
+ };
+
+ /*!
+ * \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(TwoPOneCNI, NumPhases)
+ {
+ private:
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ public:
+ static constexpr auto value = FluidSystem::numPhases;
+ static_assert(value == 2,
+ "Only fluid systems with 2 phases are supported by the 2p1cni model!");
+ };
+
+SET_TYPE_PROP(TwoPOneCNI, LocalResidual, ImmiscibleLocalResidual<TypeTag>); //! The local residual function
+
+SET_BOOL_PROP(TwoPOneCNI, EnableAdvection, true); //! The one-phase model considers advection
+SET_BOOL_PROP(TwoPOneCNI, EnableMolecularDiffusion, false); //! The one-phase model has no molecular diffusion
+
+ /*!
+ * \brief The fluid state which is used by the volume variables to
+ * store the thermodynamic state. This should be chosen
+ * appropriately for the model ((non-)isothermal, equilibrium, ...).
+ * This can be done in the problem.
+ */
+ SET_PROP(TwoPOneCNI, FluidState)
+ {
+ private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ public:
+ using type = CompositionalFluidState<Scalar, FluidSystem>;
+ // using type = ImmiscibleFluidState<Scalar, FluidSystem>;
+ };
+
+ //! Determines whether Blocking ofspurious flow is used
+ SET_BOOL_PROP(TwoPOneCNI, UseBlockingOfSpuriousFlow, false);
+
+ SET_TYPE_PROP(TwoPOneCNI, VolumeVariables, TwoPOneCVolumeVariables<TypeTag>);
+
+ //! The primary variable switch for the 2p1c model
+ SET_TYPE_PROP(TwoPOneCNI, PrimaryVariableSwitch, TwoPOneCPrimaryVariableSwitch<TypeTag>);
+
+ //! The primary variables vector for the 2p1c model
+ SET_TYPE_PROP(TwoPOneCNI, PrimaryVariables, SwitchablePrimaryVariables<TypeTag, int>);
+
+ //! Somerton is used as default model to compute the effective thermal heat conductivity
+ SET_PROP(NonIsothermal, ThermalConductivityModel)
+ {
+ private:
+ typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+ typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ public:
+ typedef ThermalConductivitySomerton<Scalar> type;
+ };
+
+ //////////////////////////////////////////////////////////////////
+ // Property values for isothermal model required for the general non-isothermal model
+ //////////////////////////////////////////////////////////////////
+ //set isothermal VolumeVariables
+ SET_TYPE_PROP(TwoPOneCNI, IsothermalVolumeVariables, TwoPOneCVolumeVariables<TypeTag>);
+
+ //set isothermal LocalResidual
+ SET_TYPE_PROP(TwoPOneCNI, IsothermalLocalResidual, ImmiscibleLocalResidual<TypeTag>);
+
+ //set isothermal Indices
+ SET_TYPE_PROP(TwoPOneCNI, IsothermalIndices, TwoPOneCIndices<TypeTag, 0>);
+
+//! the isothermal vtk output fields
+ SET_TYPE_PROP(TwoPOneCNI, IsothermalVtkOutputFields, TwoPVtkOutputFields<TypeTag>);
+
+ //set isothermal NumEq
+ SET_INT_PROP(TwoPOneCNI, IsothermalNumEq, 1);
+
+} // end namespace Properties
+} // end namespace Dumux
+
+#endif // DUMUX_2P1C_MODEL_HH
diff --git a/dumux/porousmediumflow/2p1c/implicit/newtoncontroller.hh b/dumux/porousmediumflow/2p1c/implicit/newtoncontroller.hh
deleted file mode 100644
index bc5c7e32f493186f9dbdb23cdc9e0125f8132d20..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2p1c/implicit/newtoncontroller.hh
+++ /dev/null
@@ -1,86 +0,0 @@
-// -*- 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 2p1c 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.
- */
-#ifndef DUMUX_2P1CNI_NEWTON_CONTROLLER_HH
-#define DUMUX_2P1CNI_NEWTON_CONTROLLER_HH
-
-#include "properties.hh"
-
-#include <dumux/porousmediumflow/2p2c/implicit/newtoncontroller.hh>
-#include <dumux/nonlinear/newtoncontroller.hh>
-
-namespace Dumux {
-/*!
- * \ingroup TwoPOneCModel
- * \brief A 2p1cni 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 TwoPOneCNINewtonController : public TwoPTwoCNewtonController<TypeTag>
-{
- typedef TwoPTwoCNewtonController<TypeTag> ParentType;
- typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
- typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
-
-public:
- TwoPOneCNINewtonController(const Problem &problem)
- : ParentType(problem)
- {}
-
- void newtonBeginStep()
- {
- // call the method of the base class
- ParentType::newtonBeginStep();
-
- //reset the indicator of spurious cold water fluxes into the steam zone
- this->model_().localJacobian().localResidual().resetSpuriousFlowDetected();
- }
-
-
- /*!
- * \brief Called after each Newton update
- *
- * \param uCurrentIter The current global solution vector
- * \param uLastIter The previous global solution vector
- */
- void newtonEndStep(SolutionVector &uCurrentIter,
- const SolutionVector &uLastIter)
- {
- // call the method of the base class
- ParentType::newtonEndStep(uCurrentIter, uLastIter);
-
- //Output message after each Newton Step if spurious fluxes have been blocked
- if (this->model_().localJacobian().localResidual().spuriousFlowDetected())
- std::cout <<"Spurious fluxes blocked!" <<std::endl;
- }
-
-};
-}
-
-#endif
diff --git a/dumux/porousmediumflow/2p1c/implicit/primaryvariableswitch.hh b/dumux/porousmediumflow/2p1c/implicit/primaryvariableswitch.hh
new file mode 100644
index 0000000000000000000000000000000000000000..1cb87b50f522c25b497199a3501876f37760aeda
--- /dev/null
+++ b/dumux/porousmediumflow/2p1c/implicit/primaryvariableswitch.hh
@@ -0,0 +1,156 @@
+// -*- 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 The primary variable switch for the 2p2c model
+ */
+#ifndef DUMUX_2P1C_PRIMARY_VARIABLE_SWITCH_HH
+#define DUMUX_2P1C_PRIMARY_VARIABLE_SWITCH_HH
+
+#include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup TwoPTwoCModel
+ * \brief The primary variable switch controlling the phase presence state variable
+ */
+template<class TypeTag>
+class TwoPOneCPrimaryVariableSwitch : public PrimaryVariableSwitch<TypeTag>
+{
+ friend typename Dumux::PrimaryVariableSwitch<TypeTag>;
+ using ParentType = Dumux::PrimaryVariableSwitch<TypeTag>;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using IndexType = typename GridView::IndexSet::IndexType;
+ using GlobalPosition = Dune::FieldVector<Scalar, GridView::dimensionworld>;
+
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+ enum {
+ switchIdx = Indices::switch1Idx,
+
+ pressureIdx = Indices::pressureIdx,
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+
+ twoPhases = Indices::twoPhases,
+ wPhaseOnly = Indices::wPhaseOnly,
+ nPhaseOnly = Indices::nPhaseOnly
+ };
+
+public:
+ using ParentType::ParentType;
+
+protected:
+
+ // perform variable switch at a degree of freedom location
+ bool update_(PrimaryVariables& priVars,
+ const VolumeVariables& volVars,
+ IndexType dofIdxGlobal,
+ const GlobalPosition& globalPos)
+ {
+ // evaluate primary variable switch
+ bool wouldSwitch = false;
+ int phasePresence = priVars.state();
+ int newPhasePresence = phasePresence;
+
+ // check if a primary var switch is necessary
+ if (phasePresence == twoPhases)
+ {
+ Scalar Smin = 0;
+ if (this->wasSwitched_[dofIdxGlobal])
+ Smin = -0.01;
+
+ if (volVars.saturation(nPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // gas phase disappears
+ std::cout << "Gas phase disappears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", sg: "
+ << volVars.saturation(nPhaseIdx) << std::endl;
+ newPhasePresence = wPhaseOnly;
+
+ priVars[switchIdx] = volVars.fluidState().temperature();
+ }
+ else if (volVars.saturation(wPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // water phase disappears
+ std::cout << "Wetting phase disappears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", sw: "
+ << volVars.saturation(wPhaseIdx) << std::endl;
+ newPhasePresence = nPhaseOnly;
+
+ priVars[switchIdx] = volVars.fluidState().temperature();
+ }
+
+ }
+ else if (phasePresence == wPhaseOnly)
+ {
+ const Scalar temp = volVars.fluidState().temperature();
+ const Scalar tempVap = volVars.vaporTemperature();
+
+ // if the the temperature would be larger than
+ // the vapor temperature at the given pressure, gas phase appears
+ if (temp >= tempVap)
+ {
+ wouldSwitch = true;
+ // gas phase appears
+ std::cout << "gas phase appears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << std::endl;
+
+ newPhasePresence = twoPhases;
+ priVars[switchIdx] = 0.9999; //wetting phase saturation
+ }
+ }
+
+
+ else if (phasePresence == nPhaseOnly)
+ {
+
+ const Scalar temp = volVars.fluidState().temperature();
+ const Scalar tempVap = volVars.vaporTemperature();
+
+ if (temp < tempVap)
+ {
+ wouldSwitch = true;
+ // wetting phase appears
+ std::cout << "wetting phase appears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << std::endl;
+
+
+ newPhasePresence = twoPhases;
+ priVars[switchIdx] = 0.0001; //arbitrary small value
+ }
+ }
+ priVars.setState(newPhasePresence);
+ this->wasSwitched_[dofIdxGlobal] = wouldSwitch;
+ return phasePresence != newPhasePresence;
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/porousmediumflow/2p1c/implicit/propertydefaults.hh b/dumux/porousmediumflow/2p1c/implicit/propertydefaults.hh
deleted file mode 100644
index 19b589d5ae175bd5d915d437eb15ac8a064782ca..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/2p1c/implicit/propertydefaults.hh
+++ /dev/null
@@ -1,171 +0,0 @@
-// -*- 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 TwoPOneCModel
- */
-/*!
- * \file
- *
- * \brief Defines default values for most properties required by the
- * 2p1cni fully implicit model.
- *
- *Important note: The 2p1c model requires the use of the non-isothermal extension found in dumux/implicit/nonisothermal
- */
-#ifndef DUMUX_2P1C_PROPERTY_DEFAULTS_HH
-#define DUMUX_2P1C_PROPERTY_DEFAULTS_HH
-
-#include "model.hh"
-#include "indices.hh"
-#include "volumevariables.hh"
-#include "properties.hh"
-#include "newtoncontroller.hh"
-#include "localresidual.hh"
-
-#include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
-#include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
-#include <dumux/material/fluidsystems/2pliquidvapor.hh>
-#include <dumux/material/components/constant.hh>
-#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
-#include <dumux/material/spatialparams/implicit.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 and use an static
- * assert to make sure it is 1.
- */
-SET_PROP(TwoPOneCNI, NumComponents)
-{
- private:
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-
- public:
- static const int value = FluidSystem::numComponents;
-
- static_assert(value == 1,
- "Only fluid systems with 1 component are supported by the 2p1cni model!");
-};
-
-/*!
- * \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(TwoPOneCNI, NumPhases)
-{
- private:
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-
- public:
- static const int value = FluidSystem::numPhases;
- static_assert(value == 2,
- "Only fluid systems with 2 phases are supported by the 2p1cni model!");
-};
-
-//! Use the 2p2c specific newton controller for dealing with phase switches
-SET_TYPE_PROP(TwoPOneCNI, NewtonController, TwoPOneCNINewtonController<TypeTag>);
-
-//! the FluxVariables property
-// SET_TYPE_PROP(TwoPOneCNI, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
-
-//! the upwind factor for the mobility.
-SET_SCALAR_PROP(TwoPOneCNI, ImplicitMassUpwindWeight, 1.0);
-
-//! set default mobility upwind weight to 1.0, i.e. fully upwind
-SET_SCALAR_PROP(TwoPOneCNI, ImplicitMobilityUpwindWeight, 1.0);
-
-//! The fluid system to use by default
-SET_PROP(TwoPOneCNI, FluidSystem)
-{ private:
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-public:
- typedef FluidSystems::TwoPLiquidVaporFluidsystem<Scalar, Components::Constant<1, Scalar> > type;
-};
-
-//! Determines whether Blocking ofspurious flow is used
-SET_BOOL_PROP(TwoPOneCNI, UseBlockingOfSpuriousFlow, false);
-
-//! The spatial parameters to be employed.
-//! Use ImplicitSpatialParams by default.
-SET_TYPE_PROP(TwoPOneCNI, SpatialParams, ImplicitSpatialParams<TypeTag>);
-
-// disable velocity output by default
-
-// enable gravity by default
-SET_BOOL_PROP(TwoPOneCNI, ProblemEnableGravity, true);
-
-//! default value for the forchheimer coefficient
-// Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90.
-// Actually the Forchheimer coefficient is also a function of the dimensions of the
-// porous medium. Taking it as a constant is only a first approximation
-// (Nield, Bejan, Convection in porous media, 2006, p. 10)
-SET_SCALAR_PROP(BoxModel, SpatialParamsForchCoeff, 0.55);
-
-
-//! Somerton is used as default model to compute the effective thermal heat conductivity
-SET_PROP(NonIsothermal, ThermalConductivityModel)
-{
-private:
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-public:
- typedef ThermalConductivitySomerton<Scalar> type;
-};
-
-//////////////////////////////////////////////////////////////////
-// Property values for isothermal model required for the general non-isothermal model
-//////////////////////////////////////////////////////////////////
-
-// set isothermal Model
-SET_TYPE_PROP(TwoPOneCNI, IsothermalModel, TwoPOneCModel<TypeTag>);
-
-// set isothermal FluxVariables
-SET_TYPE_PROP(TwoPOneCNI, IsothermalFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
-
-//set isothermal VolumeVariables
-SET_TYPE_PROP(TwoPOneCNI, IsothermalVolumeVariables, TwoPOneCVolumeVariables<TypeTag>);
-
-//set isothermal LocalResidual
-SET_TYPE_PROP(TwoPOneCNI, IsothermalLocalResidual, TwoPOneCLocalResidual<TypeTag>);
-
-//set isothermal Indices
-SET_PROP(TwoPOneCNI, IsothermalIndices)
-{
-
-public:
- typedef TwoPOneCIndices<TypeTag, 0> type;
-};
-
-//set isothermal NumEq
-SET_INT_PROP(TwoPOneCNI, IsothermalNumEq, 1);
-
-}
-
-}
-
-#endif
diff --git a/dumux/porousmediumflow/2p1c/implicit/volumevariables.hh b/dumux/porousmediumflow/2p1c/implicit/volumevariables.hh
index cf406f1ff303196525a8db14fbfe9f77f95760e2..019b54ff27e924f1d96a8f87bd11b25beb6f5aa6 100644
--- a/dumux/porousmediumflow/2p1c/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/2p1c/implicit/volumevariables.hh
@@ -27,9 +27,9 @@
#ifndef DUMUX_2P1C_VOLUME_VARIABLES_HH
#define DUMUX_2P1C_VOLUME_VARIABLES_HH
-#include <dumux/implicit/volumevariables.hh>
+#include <dumux/discretization/volumevariables.hh>
#include <dumux/material/fluidstates/compositional.hh>
-#include "properties.hh"
+#include "indices.hh"
namespace Dumux
{
@@ -43,19 +43,21 @@ namespace Dumux
template <class TypeTag>
class TwoPOneCVolumeVariables : 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, MaterialLaw)::Params MaterialLawParams;
-
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+ using ParentType = ImplicitVolumeVariables<TypeTag>;
+ using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+ using PermeabilityType = typename SpatialParams::PermeabilityType;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params;
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
enum {
dim = GridView::dimension,
@@ -63,7 +65,7 @@ class TwoPOneCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
numComponents = GET_PROP_VALUE(TypeTag, NumComponents),
wPhaseIdx = Indices::wPhaseIdx,
- gPhaseIdx = Indices::gPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
switch1Idx = Indices::switch1Idx,
pressureIdx = Indices::pressureIdx
@@ -73,42 +75,71 @@ class TwoPOneCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
enum {
twoPhases = Indices::twoPhases,
wPhaseOnly = Indices::wPhaseOnly,
- gPhaseOnly = Indices::gPhaseOnly,
+ nPhaseOnly = Indices::nPhaseOnly,
};
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
-
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
- enum { dofCodim = isBox ? dim : 0 };
+ using Element = typename GridView::template Codim<0>::Entity;
public:
//! The type of the object returned by the fluidState() method
- typedef Dumux::CompositionalFluidState<Scalar, FluidSystem> FluidState;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
/*!
* \copydoc ImplicitVolumeVariables::update
*/
- void update(const PrimaryVariables &priVars,
+ void update(const ElementSolutionVector &elemSol,
const Problem &problem,
const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx,
- bool isOldSol)
+ const SubControlVolume& scv)
+ {
+ ParentType::update(elemSol, problem, element, scv);
+
+ completeFluidState(elemSol, problem, element, scv, fluidState_);
+
+ /////////////
+ // calculate the remaining quantities
+ /////////////
+ const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
+
+ // 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 wetting phase saturation
+ // as parameter!
+ kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));
+ relativePermeability_[phaseIdx] = kr;
+ Valgrind::CheckDefined(relativePermeability_[phaseIdx]);
+ }
+
+ // porosity & permeability
+ porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
+ permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
+ }
+
+ /*!
+ * \copydoc ImplicitModel::completeFluidState
+ * \param isOldSol Specifies whether this is the previous solution or the current one
+ */
+ static void completeFluidState(const ElementSolutionVector& elemSol,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume& scv,
+ FluidState& fluidState)
{
- ParentType::update(priVars,
- problem,
- element,
- fvGeometry,
- scvIdx,
- isOldSol);
- // capillary pressure parameters
- const MaterialLawParams &materialParams =
- problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
+ // // capillary pressure parameters
+ const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
- int globalIdx = problem.model().dofMapper().subIndex(element, scvIdx, dofCodim);
- int phasePresence = problem.model().phasePresence(globalIdx, isOldSol);
+ const auto& priVars = ParentType::extractDofPriVars(elemSol, scv);
+ const auto phasePresence = priVars.state();
// get saturations
Scalar sw(0.0);
@@ -123,7 +154,7 @@ public:
sw = 1.0;
sg = 0.0;
}
- else if (phasePresence == gPhaseOnly)
+ else if (phasePresence == nPhaseOnly)
{
sw = 0.0;
sg = 1.0;
@@ -131,8 +162,8 @@ public:
else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
Valgrind::CheckDefined(sg);
- fluidState_.setSaturation(wPhaseIdx, sw);
- fluidState_.setSaturation(gPhaseIdx, sg);
+ fluidState.setSaturation(wPhaseIdx, sw);
+ fluidState.setSaturation(nPhaseIdx, sg);
// get gas phase pressure
const Scalar pg = priVars[pressureIdx];
@@ -144,64 +175,45 @@ public:
const Scalar pw = pg - pc;
//set pressures
- fluidState_.setPressure(wPhaseIdx, pw);
- fluidState_.setPressure(gPhaseIdx, pg);
+ fluidState.setPressure(wPhaseIdx, pw);
+ fluidState.setPressure(nPhaseIdx, pg);
// get temperature
Scalar temperature;
- if (phasePresence == wPhaseOnly || phasePresence == gPhaseOnly)
+ if (phasePresence == wPhaseOnly || phasePresence == nPhaseOnly)
temperature = priVars[switch1Idx];
else if (phasePresence == twoPhases)
- temperature = FluidSystem::vaporTemperature(fluidState_, wPhaseIdx);
+ temperature = FluidSystem::vaporTemperature(fluidState, wPhaseIdx);
else
DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
Valgrind::CheckDefined(temperature);
- fluidState_.setTemperature(temperature);
+ fluidState.setTemperature(temperature);
// set the densities
- const Scalar rhoW = FluidSystem::density(fluidState_, wPhaseIdx);
- const Scalar rhoG = FluidSystem::density(fluidState_, gPhaseIdx);
+ const Scalar rhoW = FluidSystem::density(fluidState, wPhaseIdx);
+ const Scalar rhoG = FluidSystem::density(fluidState, nPhaseIdx);
- fluidState_.setDensity(wPhaseIdx, rhoW);
- fluidState_.setDensity(gPhaseIdx, rhoG);
+ fluidState.setDensity(wPhaseIdx, rhoW);
+ fluidState.setDensity(nPhaseIdx, rhoG);
//get the viscosity and mobility
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
// Mobilities
const Scalar mu =
- FluidSystem::viscosity(fluidState_,
+ FluidSystem::viscosity(fluidState,
phaseIdx);
- fluidState_.setViscosity(phaseIdx,mu);
-
- Scalar kr;
- if (phaseIdx == wPhaseIdx)
- kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx));
- else // ATTENTION: krn requires the wetting phase saturation
- // as parameter!
- kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));
-
- mobility_[phaseIdx] = kr / mu;
- Valgrind::CheckDefined(mobility_[phaseIdx]);
+ fluidState.setViscosity(phaseIdx,mu);
}
- // porosity
- porosity_ = problem.spatialParams().porosity(element,
- fvGeometry,
- scvIdx);
- Valgrind::CheckDefined(porosity_);
-
// the enthalpies (internal energies are directly calculated in the fluidstate
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
- Scalar h = FluidSystem::enthalpy(fluidState_, phaseIdx);
- fluidState_.setEnthalpy(phaseIdx, h);
+ const Scalar h = FluidSystem::enthalpy(fluidState, phaseIdx);
+ fluidState.setEnthalpy(phaseIdx, h);
}
-
- // energy related quantities not contained in the fluid state
- asImp_().updateEnergy_(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
}
/*!
@@ -243,7 +255,7 @@ public:
*
* \param phaseIdx The phase index
*/
- Scalar pressure(const int phaseIdx) const
+ Scalar pressure(const int phaseIdx = 0) const
{ return fluidState_.pressure(phaseIdx); }
/*!
@@ -253,8 +265,8 @@ public:
* temperatures of the rock matrix and of all fluid phases are
* identical.
*/
- Scalar temperature() const
- { return fluidState_.temperature(/*phaseIdx=*/0); }
+ Scalar temperature(const int phaseIdx = 0) const
+ { return fluidState_.temperature(phaseIdx); }
/*!
* \brief Returns the effective mobility of a given phase within
@@ -264,14 +276,15 @@ public:
*/
Scalar mobility(const int phaseIdx) const
{
- return mobility_[phaseIdx];
+ return relativePermeability_[phaseIdx]/fluidState_.viscosity(phaseIdx);
}
/*!
- * \brief Returns the effective capillary pressure within the control volume.
+ * \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_.capillaryPressure(); }
+ { return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx); }
/*!
* \brief Returns the average porosity within the control volume.
@@ -279,6 +292,12 @@ public:
Scalar porosity() const
{ return porosity_; }
+ /*!
+ * \brief Returns the average permeability within the control volume in \f$[m^2]\f$.
+ */
+ const PermeabilityType& permeability() const
+ { return permeability_; }
+
/*!
* \brief Returns the vapor temperature (T_{vap}(p_g) of the fluid within the control volume.
*/
@@ -289,9 +308,9 @@ protected:
FluidState fluidState_;
private:
-
- Scalar porosity_;
- Scalar mobility_[numPhases];
+ Scalar porosity_; //!< Effective porosity within the control volume
+ PermeabilityType permeability_; //!> Effective permeability within the control volume
+ std::array<Scalar, numPhases> relativePermeability_;
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
diff --git a/dumux/porousmediumflow/2p1c/implicit/vtkoutputfields.hh b/dumux/porousmediumflow/2p1c/implicit/vtkoutputfields.hh
new file mode 100644
index 0000000000000000000000000000000000000000..12296f473e9b81ab0b6ca393499cb4672db46f6a
--- /dev/null
+++ b/dumux/porousmediumflow/2p1c/implicit/vtkoutputfields.hh
@@ -0,0 +1,47 @@
+// -*- 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 Adds vtk output fields specific to the onep model
+ */
+#ifndef DUMUX_ONEP_VTK_OUTPUT_FIELDS_HH
+#define DUMUX_ONEP_VTK_OUTPUT_FIELDS_HH
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup OneP, InputOutput
+ * \brief Adds vtk output fields specific to the onep model
+ */
+template<class TypeTag>
+class OnePVtkOutputFields
+{
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+public:
+ template <class VtkOutputModule>
+ static void init(VtkOutputModule& vtk)
+ {
+ vtk.addVolumeVariable([](const auto& volVars){ return volVars.pressure(); }, "pressure");
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/2p1c/implicit/CMakeLists.txt b/test/porousmediumflow/2p1c/implicit/CMakeLists.txt
index 436aaaac55c891152be6a21e14741cf0a592031a..4bfebfdff9b167714328cf9f013c201314d61cb5 100644
--- a/test/porousmediumflow/2p1c/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/2p1c/implicit/CMakeLists.txt
@@ -1,16 +1,20 @@
#add links to input files
add_input_file_links()
-add_dumux_test(test_boxsteaminjection test_boxsteaminjection test_boxsteaminjection.cc
- python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/steaminjectionbox-reference.vtu
- ${CMAKE_CURRENT_BINARY_DIR}/test_boxsteaminjection-00008.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_boxsteaminjection")
+dune_add_test(NAME test_2p1cni_box
+ SOURCES test_2p1c_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=TwoPOneCNIBoxProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/1p2cniccmpfaconvection-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/1p2cni_convectiontest_mpfa-00010.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p1cni_box test_boxsteaminjection.input")
-add_dumux_test(test_ccsteaminjection test_ccsteaminjection test_ccsteaminjection.cc
- python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/steaminjectioncc-reference.vtu
- ${CMAKE_CURRENT_BINARY_DIR}/test_ccsteaminjection-00009.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_ccsteaminjection")
+dune_add_test(NAME test_2p1cni_tpfa
+ SOURCES test_2p1c_fv.cc
+ COMPILE_DEFINITIONS TYPETAG=TwoPOneCNICCTpfaProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/1p2cniccmpfaconvection-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/1p2cni_convectiontest_mpfa-00010.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_2p1cni_tpfa test_ccsteaminjection.input")
diff --git a/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh b/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh
index 97a3347bd0eba5a200f2ec79e9a2983240d388ad..9985d0608d33ce945f216b1fe85543f0e7f4c51b 100644
--- a/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh
+++ b/test/porousmediumflow/2p1c/implicit/steaminjectionproblem.hh
@@ -25,13 +25,18 @@
#ifndef DUMUX_STEAM_INJECTIONPROBLEM_HH
#define DUMUX_STEAM_INJECTIONPROBLEM_HH
+#include <dumux/discretization/cellcentered/tpfa/properties.hh>
+#include <dumux/discretization/box/properties.hh>
#include <dumux/porousmediumflow/2p1c/implicit/model.hh>
-#include <dumux/porousmediumflow/implicit/problem.hh>
+#include <dumux/porousmediumflow/problem.hh>
+
+#include <dumux/material/fluidsystems/2pliquidvapor.hh>
-#include "steaminjectionspatialparams.hh"
#include <dumux/material/components/tabulatedcomponent.hh>
#include <dumux/material/components/h2o.hh>
+#include "steaminjectionspatialparams.hh"
+
namespace Dumux
{
template <class TypeTag>
@@ -40,32 +45,28 @@ class InjectionProblem;
namespace Properties
{
NEW_TYPE_TAG(InjectionProblem, INHERITS_FROM(TwoPOneCNI, InjectionProblemSpatialParams));
-NEW_TYPE_TAG(InjectionBoxProblem, INHERITS_FROM(BoxModel, InjectionProblem));
-NEW_TYPE_TAG(InjectionCCProblem, INHERITS_FROM(CCModel, InjectionProblem));
+NEW_TYPE_TAG(TwoPOneCNIBoxProblem, INHERITS_FROM(BoxModel, InjectionProblem));
+NEW_TYPE_TAG(TwoPOneCNICCTpfaProblem, INHERITS_FROM(CCTpfaModel, InjectionProblem));
SET_TYPE_PROP(InjectionProblem, Grid, Dune::YaspGrid<2>);
// Set the problem property
-SET_PROP(InjectionProblem, Problem)
-{
- typedef Dumux::InjectionProblem<TypeTag> type;
-};
+SET_TYPE_PROP(InjectionProblem, Problem, InjectionProblem<TypeTag>);
-// Set the fluid system
+
+// Set fluid configuration
SET_PROP(InjectionProblem, FluidSystem)
{
private:
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar> > H2OType ;
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using H2OType = Dumux::TabulatedComponent<Scalar, Dumux::H2O<Scalar> >;
public:
- typedef Dumux::FluidSystems::TwoPLiquidVaporFluidsystem<Scalar, H2OType > type;
+ using type = Dumux::FluidSystems::TwoPLiquidVaporFluidsystem<Scalar, H2OType >;
};
-// Enable gravity
-SET_BOOL_PROP(InjectionProblem, ProblemEnableGravity, true);
//Define whether spurious cold-water flow into the steam is blocked
-SET_BOOL_PROP(InjectionProblem, UseBlockingOfSpuriousFlow, true);
+SET_BOOL_PROP(InjectionProblem, UseBlockingOfSpuriousFlow, false);
}
//TODO: Names
@@ -78,16 +79,29 @@ SET_BOOL_PROP(InjectionProblem, UseBlockingOfSpuriousFlow, true);
*
* */
template <class TypeTag>
-class InjectionProblem : public ImplicitPorousMediaProblem<TypeTag>
+class InjectionProblem : public PorousMediumFlowProblem<TypeTag>
{
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GridView::Grid Grid;
-
- typedef ImplicitPorousMediaProblem<TypeTag> ParentType;
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
// copy some indices for convenience
- typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
pressureIdx = Indices::pressureIdx,
switch1Idx = Indices::switch1Idx,
@@ -96,32 +110,20 @@ class InjectionProblem : public ImplicitPorousMediaProblem<TypeTag>
// phase and component indices
wPhaseIdx = Indices::wPhaseIdx,
- gPhaseIdx = Indices::gPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
// Phase State
wPhaseOnly = Indices::wPhaseOnly,
- gPhaseOnly = Indices::gPhaseOnly,
+ nPhaseOnly = Indices::nPhaseOnly,
twoPhases = Indices::twoPhases,
-
- // Grid and world dimension
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld
};
- typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
- typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
- typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
-
- typedef typename GridView::template Codim<0>::Entity Element;
- typedef typename GridView::template Codim<dim>::Entity Vertex;
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
- typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+ static constexpr bool isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box;
+ using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
- typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { dofCodim = isBox ? dim : 0 };
public:
@@ -132,15 +134,9 @@ public:
* \param gridView The grid view
*/
- InjectionProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView)
+ InjectionProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry)
{
-
- name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag,
- std::string,
- Problem,
- Name);
-
FluidSystem::init();
}
@@ -149,20 +145,15 @@ public:
*/
// \{
- /*!
- * \brief The problem name.
- *
- * This is used as a prefix for files generated by the simulation.
- */
- const std::string name() const
- { return name_; }
-
- void sourceAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- values = 0.0;
- }
+ //! \copydoc Dumux::FVProblem::source()
+ Sources source(const Element &element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume &scv) const
+ {
+ return Sources(0.0);
+ }
/*!
* \name Boundary conditions
@@ -171,58 +162,71 @@ public:
/*!
* \brief Specifies which kind of boundary condition should be
- * used for which equation on a given boundary segment.
+ * used for which equation on a given boundary segment
*
- * \param bcTypes The boundary types for the conservation equations
- * \param globalPos The position for which the bc type should be evaluated
+ * \param globalPos The global position
*/
- void boundaryTypesAtPos(BoundaryTypes &bcTypes,
- const GlobalPosition &globalPos) const
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
- if(globalPos[1] > this->bBoxMax()[1] - eps_ || globalPos[0] > this->bBoxMax()[0] - eps_)
+ BoundaryTypes bcTypes;
+
+ if(globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - eps_ || globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_)
bcTypes.setAllDirichlet();
else
bcTypes.setAllNeumann();
+
+ return bcTypes;
}
/*!
- * \brief Evaluate the boundary conditions for a dirichlet
- * boundary segment.
+ * \brief Evaluates the boundary conditions for a Dirichlet
+ * boundary segment
*
- * \param values The dirichlet values for the primary variables
- * \param globalPos The position for which the bc type should be evaluated
- *
- * For this method, the \a values parameter stores primary variables.
+ * \param globalPos The global position
*/
- void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
- initial_(values, globalPos);
+ return initialAtPos(globalPos);
}
/*!
* \brief Evaluate the boundary conditions for a neumann
* boundary segment.
*
- * \param values The dirichlet values for the primary variables
- * \param globalPos The position for which the bc type should be evaluated
+ * This is the method for the case where the Neumann condition is
+ * potentially solution dependent and requires some quantities that
+ * are specific to the fully-implicit method.
+ *
+ * \param values The neumann values for the conservation equations in units of
+ * \f$ [ \textnormal{unit of conserved quantity} / (m^2 \cdot s )] \f$
+ * \param element The finite element
+ * \param fvGeometry The finite-volume geometry
+ * \param elemVolVars All volume variables for the element
+ * \param scvf The sub control volume face
*
- * For this method, the \a values parameter stores the mass flux
+ * For this method, the \a values parameter stores the flux
* in normal direction of each phase. Negative values mean influx.
+ * E.g. for the mass balance that would the mass flux in \f$ [ kg / (m^2 \cdot s)] \f$.
*/
- void neumannAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
+ ResidualVector neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
{
- values = 0.0;
+ ResidualVector values(0.0);
- if (globalPos[0] < eps_)
+ const auto& ipGlobal = scvf.ipGlobal();
+
+ if (ipGlobal[0] < eps_)
{
- if(globalPos[1] > 2.0 - eps_ && globalPos[1] < 3.0 + eps_)
+ if(ipGlobal[1] > 2.0 - eps_ && ipGlobal[1] < 3.0 + eps_)
{
- Scalar massRate = 1e-1;
- values[Indices::conti0EqIdx] = -massRate;
- values[Indices::energyEqIdx] = -massRate * 2690e3;
+ const Scalar massRate = 1e-1;
+ values[Indices::conti0EqIdx] = -massRate;
+ values[Indices::energyEqIdx] = -massRate * 2690e3;
}
}
+ return values;
}
// \}
@@ -233,48 +237,26 @@ public:
// \{
/*!
- * \brief Evaluate the initial value for a control volume.
- *
- * \param values The initial values for the primary variables
- * \param globalPos The position for which the initial condition should be evaluated
+ * \brief Evaluates the initial values for a control volume
*
- * For this method, the \a values parameter stores primary
- * variables.
- */
- void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
- {
- initial_(values, globalPos);
- }
-
- /*!
- * \brief Return the initial phase state inside a control volume.
- *
- * \param vert The vertex
- * \param globalIdx The index of the global vertex
* \param globalPos The global position
*/
- int initialPhasePresence(const Vertex &vert,
- int &globalIdx,
- const GlobalPosition &globalPos) const
+ PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
- return wPhaseOnly;
- }
+ PrimaryVariables values(0.0);
+
+ const Scalar densityW = 1000.0;
+ values[pressureIdx] = 101300.0 + (this->fvGridGeometry().bBoxMax()[1] - globalPos[1])*densityW*9.81; // hydrostatic pressure
+ values[switch1Idx] = 283.13;
+ values.setState(wPhaseOnly);
+ return values;
+ }
private:
- // internal method for the initial condition (reused for the
- // dirichlet conditions!)
- void initial_(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- Scalar densityW = 1000.0;
- values[pressureIdx] = 101300.0 + (this->bBoxMax()[1] - globalPos[1])*densityW*9.81; // hydrostatic pressure
- values[switch1Idx] = 283.13;
- }
static constexpr Scalar eps_ = 1e-6;
- std::string name_;
};
} //end namespace
diff --git a/test/porousmediumflow/2p1c/implicit/steaminjectionspatialparams.hh b/test/porousmediumflow/2p1c/implicit/steaminjectionspatialparams.hh
index 43f02fc690c822b56ef00a58b59f5808f7a0a756..48a00ca88b29986d2095d11a6d8cc7ddcc25478c 100644
--- a/test/porousmediumflow/2p1c/implicit/steaminjectionspatialparams.hh
+++ b/test/porousmediumflow/2p1c/implicit/steaminjectionspatialparams.hh
@@ -64,34 +64,36 @@ SET_PROP(InjectionProblemSpatialParams, MaterialLaw)
template<class TypeTag>
class InjectionProblemSpatialParams : public ImplicitSpatialParams<TypeTag>
{
- typedef ImplicitSpatialParams<TypeTag> ParentType;
-
- typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
- typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
- typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
- typedef typename Grid::ctype CoordScalar;
- enum {
- dim=GridView::dimension,
- dimWorld=GridView::dimensionworld
- };
-
- typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
- typedef typename GridView::template Codim<0>::Entity Element;
-
- typedef Dune::FieldMatrix<Scalar,dim,dim> FieldMatrix;
-
-
+ using ParentType = ImplicitSpatialParams<TypeTag>;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+
+ using MaterialLawParams = typename MaterialLaw::Params;
+ using CoordScalar = typename GridView::ctype;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+
+ static constexpr int dim = GridView::dimension;
+ static constexpr int dimWorld = GridView::dimensionworld;
+ static constexpr int wPhaseIdx = FluidSystem::wPhaseIdx;
+
+ using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
+ using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
public:
- typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
- typedef typename MaterialLaw::Params MaterialLawParams;
+ using PermeabilityType = DimWorldMatrix;
/*!
* \brief The constructor
*
* \param gridView The grid view
*/
- InjectionProblemSpatialParams(const GridView &gridView)
- : ParentType(gridView)
+ InjectionProblemSpatialParams(const Problem& problem)
+ : ParentType(problem)
{
// set Van Genuchten Parameters
materialParams_.setSwr(0.1);
@@ -101,18 +103,13 @@ public:
}
/*!
- * \brief Apply the intrinsic permeability tensor to a pressure
- * potential gradient.
+ * \brief Returns the hydraulic conductivity \f$[m^2]\f$
*
- * \param element The current finite element
- * \param fvElemGeom The current finite volume geometry of the element
- * \param scvIdx The index of the sub-control volume
+ * \param globalPos The global position
*/
- const FieldMatrix intrinsicPermeability(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
+ DimWorldMatrix permeabilityAtPos(const GlobalPosition& globalPos) const
{
- FieldMatrix permMatrix(0.0);
+ DimWorldMatrix permMatrix(0.0);
// intrinsic permeability
permMatrix[0][0] = 1e-9;
@@ -124,78 +121,62 @@ public:
/*!
* \brief Define the porosity \f$[-]\f$ of the spatial parameters
*
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub-control volume where
- * the porosity needs to be defined
+ * \param globalPos The global position
*/
- double porosity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
{
return 0.4;
}
-
/*!
- * \brief return the parameter object for the Brooks-Corey material law which depends on the position
+ * \brief Returns the parameter object for the capillary-pressure/
+ * saturation material law
*
- * \param element The current finite element
- * \param fvGeometry The current finite volume geometry of the element
- * \param scvIdx The index of the sub-control volume
+ * \param globalPos The global position
*/
- const MaterialLawParams& materialLawParams(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ const MaterialLawParams& materialLawParamsAtPos(const GlobalPosition& globalPos) const
{
return materialParams_;
}
-
/*!
* \brief Returns the heat capacity \f$[J / (kg K)]\f$ of the rock matrix.
*
* This is only required for non-isothermal models.
*
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub-control volume
+ * \param element The element
+ * \param scv The sub control volume
+ * \param elemSol The element solution vector
*/
- Scalar solidHeatCapacity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
- {
- return 850.0; // specific heat capacity of granite [J / (kg K)]
- }
+ Scalar solidHeatCapacity(const Element &element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ { return 850.0; /*specific heat capacity of granite [J / (kg K)]*/ }
/*!
* \brief Returns the mass density \f$[kg / m^3]\f$ of the rock matrix.
*
* This is only required for non-isothermal models.
*
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub-control volume
+ * \param element The element
+ * \param scv The sub control volume
+ * \param elemSol The element solution vector
*/
Scalar solidDensity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
- {
- return 2650.0; // density of granite [kg/m^3]
- }
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ { return 2650; /*density of granite [kg/m^3]*/ }
/*!
- * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
+ * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
*
- * \param element The finite element
- * \param fvGeometry The finite volume geometry
- * \param scvIdx The local index of the sub-control volume
+ * \param element The element
+ * \param scv The sub control volume
+ * \param elemSol The element solution vector
*/
Scalar solidThermalConductivity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
- {
- return 2.8;
- }
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ { return 2.8; }
private:
MaterialLawParams materialParams_;
diff --git a/test/porousmediumflow/2p1c/implicit/test_2p1c_fv.cc b/test/porousmediumflow/2p1c/implicit/test_2p1c_fv.cc
new file mode 100644
index 0000000000000000000000000000000000000000..94aafde415fa6cbe7771fdff0bf8dd5ad5607499
--- /dev/null
+++ b/test/porousmediumflow/2p1c/implicit/test_2p1c_fv.cc
@@ -0,0 +1,218 @@
+// -*- 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 test for the 1pnc model
+ */
+ #include <config.h>
+
+ #include "steaminjectionproblem.hh"
+
+ #include <ctime>
+ #include <iostream>
+
+ #include <dune/common/parallel/mpihelper.hh>
+ #include <dune/common/timer.hh>
+ #include <dune/grid/io/file/dgfparser/dgfexception.hh>
+ #include <dune/grid/io/file/vtk.hh>
+ #include <dune/istl/io.hh>
+
+ #include <dumux/discretization/methods.hh>
+
+ #include <dumux/common/properties.hh>
+ #include <dumux/common/parameters.hh>
+ #include <dumux/common/valgrind.hh>
+ #include <dumux/common/dumuxmessage.hh>
+ #include <dumux/common/defaultusagemessage.hh>
+ #include <dumux/common/parameterparser.hh>
+
+ #include <dumux/porousmediumflow/compositional/privarswitchnewtoncontroller.hh>
+ #include <dumux/linear/seqsolverbackend.hh>
+ #include <dumux/nonlinear/newtonmethod.hh>
+
+ #include <dumux/assembly/fvassembler.hh>
+
+ #include <dumux/io/vtkoutputmodule.hh>
+
+ int main(int argc, char** argv) try
+ {
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = TTAG(TYPETAG);
+
+ ////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // print dumux start message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/true);
+
+ // initialize parameter tree
+ Parameters::init(argc, argv);
+
+ //////////////////////////////////////////////////////////////////////
+ // try to create a grid (from the given grid file or the input file)
+ /////////////////////////////////////////////////////////////////////
+
+ using GridCreator = typename GET_PROP_TYPE(TypeTag, GridCreator);
+ GridCreator::makeGrid();
+ GridCreator::loadBalance();
+
+ ////////////////////////////////////////////////////////////
+ // run instationary non-linear problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = GridCreator::grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ auto problem = std::make_shared<Problem>(fvGridGeometry);
+
+ // the solution vector
+ using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+ SolutionVector x(fvGridGeometry->numDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+
+ // the grid variables
+ using GridVariables = typename GET_PROP_TYPE(TypeTag, GridVariables);
+ auto gridVariables = std::make_shared<GridVariables>(problem, fvGridGeometry);
+ gridVariables->init(x, xOld);
+
+ // get some time loop parameters
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+ auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+ auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+
+ // intialize the vtk output module
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0.0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ // the assembler with time loop for instationary problem
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, fvGridGeometry, gridVariables, timeLoop);
+
+ // the linear solver
+ using LinearSolver = ILU0BiCGSTABBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the non-linear solver
+ using NewtonController = Dumux::PriVarSwitchNewtonController<TypeTag>;
+ auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+ NewtonMethod<NewtonController, Assembler, LinearSolver> nonLinearSolver(newtonController, assembler, linearSolver);
+
+ // time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ // try solving the non-linear system
+ for (int i = 0; i < maxDivisions; ++i)
+ {
+ // linearize & solve
+ auto converged = nonLinearSolver.solve(x);
+
+ if (converged)
+ break;
+
+ if (!converged && i == maxDivisions-1)
+ DUNE_THROW(Dune::MathError,
+ "Newton solver didn't converge after "
+ << maxDivisions
+ << " time-step divisions. dt="
+ << timeLoop->timeStepSize()
+ << ".\nThe solutions of the current and the previous time steps "
+ << "have been saved to restart files.");
+ }
+
+ // make the new solution the old solution
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt as suggested by newton controller
+ timeLoop->setTimeStepSize(newtonController->suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/false);
+
+ return 0;
+
+ }
+ catch (Dumux::ParameterException &e)
+ {
+ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+ return 1;
+ }
+ catch (Dune::DGFException & e)
+ {
+ std::cerr << "DGF exception thrown (" << e <<
+ "). Most likely, the DGF file name is wrong "
+ "or the DGF file is corrupted, "
+ "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+ << " ---> Abort!" << std::endl;
+ return 2;
+ }
+ catch (Dune::Exception &e)
+ {
+ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+ return 3;
+ }
+ catch (...)
+ {
+ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+ return 4;
+ }
diff --git a/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.cc b/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.cc
deleted file mode 100644
index bd83890846bb39f76a3e9233e83fa289156e7bfe..0000000000000000000000000000000000000000
--- a/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.cc
+++ /dev/null
@@ -1,62 +0,0 @@
-// -*- 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 Test for the 2p1cni box model
- */
-#include "config.h"
-#include "steaminjectionproblem.hh"
-#include <dumux/common/start.hh>
-
-/*!
- * \brief Provides an interface for customizing error messages associated with
- * reading in parameters.
- *
- * \param progName The name of the program, that was tried to be started.
- * \param errorMsg The error message that was issued by the start function.
- * Comprises the thing that went wrong and a general help message.
- */
-void usage(const char *progName, const std::string &errorMsg)
-{
- if (errorMsg.size() > 0) {
- std::string errorMessageOut = "\nUsage: ";
- errorMessageOut += progName;
- errorMessageOut += " [options]\n";
- errorMessageOut += errorMsg;
- errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
- "\t-TimeManager.TEnd End of the simulation [s] \n"
- "\t-TimeManager.DtInitial Initial timestep size [s] \n"
- "\t-Grid.UpperRight coordinates of the upper right corner of the grid [m] \n"
- "\t-Grid.Cells Number of cells in respective coordinate directions\n"
- "\t-Problem.Name String for naming of the output files \n"
- "\n";
-
- std::cout << errorMessageOut << std::endl;
- }
-}
-
-////////////////////////
-// the main function
-////////////////////////
-int main(int argc, char** argv)
-{
- typedef TTAG(InjectionBoxProblem) TypeTag;
- return Dumux::start<TypeTag>(argc, argv, usage);
-}
diff --git a/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.input b/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.input
index 4696345154e120977b7f4c2d8377e28ec4cdfc07..59aa81a23af877876c9ffba67dd3c141602d507b 100644
--- a/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.input
+++ b/test/porousmediumflow/2p1c/implicit/test_boxsteaminjection.input
@@ -1,4 +1,4 @@
-[TimeManager]
+[TimeLoop]
DtInitial = 50 # [s]
TEnd = 300 # [s]
diff --git a/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.cc b/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.cc
deleted file mode 100644
index 1273bd90ede4491753202a2f9cfef6b72e669adc..0000000000000000000000000000000000000000
--- a/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.cc
+++ /dev/null
@@ -1,62 +0,0 @@
-// -*- 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 Test for the 2p1cni cell-centered model
- */
-#include "config.h"
-#include "steaminjectionproblem.hh"
-#include <dumux/common/start.hh>
-
-/*!
- * \brief Provides an interface for customizing error messages associated with
- * reading in parameters.
- *
- * \param progName The name of the program, that was tried to be started.
- * \param errorMsg The error message that was issued by the start function.
- * Comprises the thing that went wrong and a general help message.
- */
-void usage(const char *progName, const std::string &errorMsg)
-{
- if (errorMsg.size() > 0) {
- std::string errorMessageOut = "\nUsage: ";
- errorMessageOut += progName;
- errorMessageOut += " [options]\n";
- errorMessageOut += errorMsg;
- errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
- "\t-TimeManager.TEnd End of the simulation [s] \n"
- "\t-TimeManager.DtInitial Initial timestep size [s] \n"
- "\t-Grid.UpperRight coordinates of the upper right corner of the grid [m] \n"
- "\t-Grid.Cells Number of cells in respective coordinate directions\n"
- "\t-Problem.Name String for naming of the output files \n"
- "\n";
-
- std::cout << errorMessageOut << std::endl;
- }
-}
-
-////////////////////////
-// the main function
-////////////////////////
-int main(int argc, char** argv)
-{
- typedef TTAG(InjectionCCProblem) TypeTag;
- return Dumux::start<TypeTag>(argc, argv, usage);
-}
diff --git a/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.input b/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.input
index c5b8d8556b1fd4b9c149baa9212b1684d3db858d..fabf21bed2ff303018d5c1d26fc55ec3f9f00e87 100644
--- a/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.input
+++ b/test/porousmediumflow/2p1c/implicit/test_ccsteaminjection.input
@@ -1,4 +1,4 @@
-[TimeManager]
+[TimeLoop]
DtInitial = 50 # [s]
TEnd = 300 # [s]