diff --git a/dumux/common/properties.hh b/dumux/common/properties.hh
index def2fb54f87c219d8170bc2d5fe24a0b489c3980..95dbf79f5aef913114207e9686bae63fce6f134d 100644
--- a/dumux/common/properties.hh
+++ b/dumux/common/properties.hh
@@ -174,6 +174,11 @@ NEW_PROP_TAG(SetMoleFractionsForWettingPhase); //!< Set the mole fraction in
//////////////////////////////////////////////////////////////
NEW_PROP_TAG(EnableWaterDiffusionInAir); //!< Property for turning Richards into extended Richards
+//////////////////////////////////////////////////////////////
+// Additional properties used by the 3pwateroil model:
+//////////////////////////////////////////////////////////////
+NEW_PROP_TAG(OnlyGasPhaseCanDisappear); //!< reduces the phasestates to threePhases and wnPhaseOnly
+
/////////////////////////////////////////////////////////////
// Properties used by the staggered-grid discretization method
/////////////////////////////////////////////////////////////
diff --git a/dumux/porousmediumflow/3pwateroil/fluxvariables.hh b/dumux/porousmediumflow/3pwateroil/fluxvariables.hh
deleted file mode 100644
index 96a6bd3df740cc2c8007d0feaa2a5685f61fc0ca..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/3pwateroil/fluxvariables.hh
+++ /dev/null
@@ -1,302 +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 This file contains the data which is required to calculate
- * all fluxes of components over a face of a finite volume for
- * the three-phase, two-component model.
- */
-#ifndef DUMUX_3P2CNI_FLUX_VARIABLES_HH
-#define DUMUX_3P2CNI_FLUX_VARIABLES_HH
-
-#include <dumux/common/math.hh>
-#include <dumux/common/spline.hh>
-
-#include "properties.hh"
-
-namespace Dumux
-{
-
-/*!
- * \ingroup ThreePWaterOilModel
- * \ingroup ImplicitFluxVariables
- * \brief This template class contains the data which is required to
- * calculate all fluxes of components over a face of a finite
- * volume for the three-phase, two-component model.
- *
- * This means pressure and concentration gradients, phase densities at
- * the integration point, etc.
- */
-template <class TypeTag>
-class ThreePWaterOilFluxVariables : public GET_PROP_TYPE(TypeTag, BaseFluxVariables)
-{
- friend typename GET_PROP_TYPE(TypeTag, BaseFluxVariables); // be friends with base class
- using BaseFluxVariables = typename GET_PROP_TYPE(TypeTag, BaseFluxVariables);
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using EffectiveDiffusivityModel = typename GET_PROP_TYPE(TypeTag, EffectiveDiffusivityModel);
-
- using Element = typename GridView::template Codim<0>::Entity;
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-
- enum {
- dim = GridView::dimension,
- dimWorld = GridView::dimensionworld,
- numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
- numComponents = GET_PROP_VALUE(TypeTag, NumComponents)
- };
-
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
- using SCVFace = typename FVElementGeometry::SubControlVolumeFace;
-
- using DimVector = Dune::FieldVector<Scalar, dim>;
- using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
-
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
- enum {
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx,
- gPhaseIdx = Indices::gPhaseIdx,
-
- wCompIdx = Indices::wCompIdx,
- nCompIdx = Indices::nCompIdx,
- };
-
-public:
- /*!
- * \brief Compute / update the flux variables
- *
- * \param problem The problem
- * \param element The finite element
- * \param fvGeometry The finite-volume geometry
- * \param fIdx The local index of the SCV (sub-control-volume) face
- * \param elemVolVars The volume variables of the current element
- * \param onBoundary A boolean variable to specify whether the flux variables
- * are calculated for interior SCV faces or boundary faces, default=false
- */
- void update(const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const int fIdx,
- const ElementVolumeVariables &elemVolVars,
- const bool onBoundary = false)
- {
- BaseFluxVariables::update(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary);
- calculatePorousDiffCoeff_(problem, element, elemVolVars);
- }
-
-private:
- void calculateGradients_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
- // initialize to zero
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- density_[phaseIdx] = Scalar(0);
- molarDensity_[phaseIdx] = Scalar(0);
- massFractionCompWGrad_[phaseIdx] = Scalar(0);
- massFractionCompNGrad_[phaseIdx] = Scalar(0);
- moleFractionCompWGrad_[phaseIdx] = Scalar(0);
- moleFractionCompNGrad_[phaseIdx] = Scalar(0);
- }
-
- BaseFluxVariables::calculateGradients_(problem, element, elemVolVars);
-
- // calculate gradients
- DimVector tmp(0.0);
- for (int idx = 0; idx < this->face().numFap; idx++) // loop over adjacent vertices
- {
- // FE gradient at vertex idx
- const DimVector &feGrad = this->face().grad[idx];
-
- // index for the element volume variables
- int volVarsIdx = this->face().fapIndices[idx];
-
- // the concentration gradient of the components
- // component in the phases
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(wPhaseIdx, wCompIdx);
- massFractionCompWGrad_[wPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(nPhaseIdx, wCompIdx);
- massFractionCompWGrad_[nPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(gPhaseIdx, wCompIdx);
- massFractionCompWGrad_[gPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(wPhaseIdx, nCompIdx);
- massFractionCompNGrad_[wPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(nPhaseIdx, nCompIdx);
- massFractionCompNGrad_[nPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().massFraction(gPhaseIdx, nCompIdx);
- massFractionCompNGrad_[gPhaseIdx] += tmp;
-
- // the molar concentration gradients of the components
- // in the phases
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(wPhaseIdx, wCompIdx);
- moleFractionCompWGrad_[wPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, wCompIdx);
- moleFractionCompWGrad_[nPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(gPhaseIdx, wCompIdx);
- moleFractionCompWGrad_[gPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(wPhaseIdx, nCompIdx);
- moleFractionCompNGrad_[wPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(nPhaseIdx, nCompIdx);
- moleFractionCompNGrad_[nPhaseIdx] += tmp;
-
- tmp = feGrad;
- tmp *= elemVolVars[volVarsIdx].fluidState().moleFraction(gPhaseIdx, nCompIdx);
- moleFractionCompNGrad_[gPhaseIdx] += tmp;
-
- }
- // calculate temperature gradient using finite element
- // gradients
- DimVector temperatureGrad(0);
- for (int idx = 0; idx < this->face().numFap; idx++)
- {
- tmp = this->face().grad[idx];
-
- // index for the element volume variables
- int volVarsIdx = this->face().fapIndices[idx];
-
- tmp *= elemVolVars[volVarsIdx].temperature();
- temperatureGrad += tmp;
- }
- }
-
- void calculatePorousDiffCoeff_(const Problem &problem,
- const Element &element,
- const ElementVolumeVariables &elemVolVars)
- {
-
- const VolumeVariables &volVarsI = elemVolVars[this->face().i];
- const VolumeVariables &volVarsJ = elemVolVars[this->face().j];
-
- // the effective diffusion coefficients at vertex i and j
- Scalar diffCoeffI;
- Scalar diffCoeffJ;
-
-
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
- {
- // make sure to calculate only diffusion coefficients
- // for phases which exist in both finite volumes
- /* \todo take care: This should be discussed once again
- * as long as a meaningful value can be found for the required mole fraction
- * diffusion should work even without this one here */
- if (volVarsI.saturation(phaseIdx) <= 0 ||
- volVarsJ.saturation(phaseIdx) <= 0)
- {
- porousDiffCoeff_[phaseIdx] = 0.0;
- }
- else
- {
-
- diffCoeffI = EffectiveDiffusivityModel::effectiveDiffusivity(volVarsI.porosity(),
- volVarsI.saturation(phaseIdx),
- volVarsI.diffCoeff(phaseIdx));
-
- diffCoeffJ = EffectiveDiffusivityModel::effectiveDiffusivity(volVarsJ.porosity(),
- volVarsJ.saturation(phaseIdx),
- volVarsJ.diffCoeff(phaseIdx));
-
- porousDiffCoeff_[phaseIdx] = harmonicMean(diffCoeffI, diffCoeffJ);
- }
- }
- }
-
-public:
- /*!
- * \brief The binary diffusion coefficient for each fluid phase.
- *
- * \param phaseIdx The phase index
- */
- DUNE_DEPRECATED_MSG("porousDiffCoeff() is deprecated. Use porousDiffCoeff(phaseIdx) instead.")
- Dune::FieldVector<Scalar, numPhases> porousDiffCoeff() const
- { return porousDiffCoeff_; };
-
- /*!
- * \brief The binary diffusion coefficient for each fluid phase.
- *
- * \param phaseIdx The phase index
- */
- Scalar porousDiffCoeff(int phaseIdx) const
- { return porousDiffCoeff_[phaseIdx];}
-
- /*!
- * \brief Return density \f$\mathrm{[kg/m^3]}\f$ of a phase.
- */
- Scalar density(int phaseIdx) const
- { return density_[phaseIdx]; }
-
- /*!
- * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ of a phase.
- */
- Scalar molarDensity(int phaseIdx) const
- { return molarDensity_[phaseIdx]; }
-
- const DimVector &massFractionCompWGrad(int phaseIdx) const
- {return massFractionCompWGrad_[phaseIdx];}
-
- const DimVector &massFractionCompNGrad(int phaseIdx) const
- { return massFractionCompNGrad_[phaseIdx]; };
-
- const DimVector &moleFractionCompWGrad(int phaseIdx) const
- { return moleFractionCompWGrad_[phaseIdx]; };
-
- const DimVector &moleFractionCompNGrad(int phaseIdx) const
- { return moleFractionCompNGrad_[phaseIdx]; };
-
-protected:
- // gradients
- DimVector massFractionCompWGrad_[numPhases];
- DimVector massFractionCompNGrad_[numPhases];
- DimVector moleFractionCompWGrad_[numPhases];
- DimVector moleFractionCompNGrad_[numPhases];
-
- // density of each face at the integration point
- Scalar density_[numPhases], molarDensity_[numPhases];
-
- // the diffusivity matrix for the porous medium
- Dune::FieldVector<Scalar, numPhases> porousDiffCoeff_;
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/3pwateroil/indices.hh b/dumux/porousmediumflow/3pwateroil/indices.hh
index 3ac1b47ae70dc74cf6cc166d2580e5b3d3d97159..52541a3205ba01e39c131a2ab268d5f5bb7dea36 100644
--- a/dumux/porousmediumflow/3pwateroil/indices.hh
+++ b/dumux/porousmediumflow/3pwateroil/indices.hh
@@ -18,19 +18,18 @@
*****************************************************************************/
/*!
* \file
+ * \ingroup ThreePWaterOilModel
* \brief Defines the indices required for the 3p2cni model.
*/
#ifndef DUMUX_3P2CNI_INDICES_HH
#define DUMUX_3P2CNI_INDICES_HH
-#include "properties.hh"
-
+#include <dumux/common/properties.hh>
namespace Dumux
{
/*!
* \ingroup ThreePWaterOilModel
- * \ingroup ImplicitIndices
* \brief The indices for the isothermal 3p2cni model.
*
* \tparam formulation The formulation, only pgSwSn
diff --git a/dumux/porousmediumflow/3pwateroil/localresidual.hh b/dumux/porousmediumflow/3pwateroil/localresidual.hh
index 8aeea475c83e558acdfa18599aa31ffdad3b7ed9..4bb312fd7890e2431b25ff0b0980ad98afafdd59 100644
--- a/dumux/porousmediumflow/3pwateroil/localresidual.hh
+++ b/dumux/porousmediumflow/3pwateroil/localresidual.hh
@@ -18,38 +18,47 @@
*****************************************************************************/
/*!
* \file
- *
+ * \ingroup ThreePWaterOilModel
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the three-phase three-component fully implicit model.
*/
#ifndef DUMUX_3P2CNI_LOCAL_RESIDUAL_HH
#define DUMUX_3P2CNI_LOCAL_RESIDUAL_HH
-#include "properties.hh"
+#include <dumux/common/properties.hh>
#include <dumux/porousmediumflow/3p3c/localresidual.hh>
namespace Dumux
{
/*!
* \ingroup ThreePWaterOilModel
- * \ingroup ImplicitLocalResidual
- * \brief Element-wise calculation of the Jacobian matrix for problems
- * using the three-phase three-component fully implicit model.
+ * \brief Element-wise calculation of the local residual for problems
+ * using the ThreePWaterOil fully implicit model.
*
- * This class is used to fill the gaps in BoxLocalResidual for the 3P2C flow.
+ * This class is used to fill the gaps in the CompositionalLocalResidual for the 3PWaterOil flow.
*/
template<class TypeTag>
class ThreePWaterOilLocalResidual: public ThreePThreeCLocalResidual<TypeTag>
{
protected:
+ using ParentType = ThreePThreeCLocalResidual<TypeTag>;
using Implementation = typename GET_PROP_TYPE(TypeTag, LocalResidual);
-
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using ResidualVector = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
+ using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-
+ using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using ElementBoundaryTypes = typename GET_PROP_TYPE(TypeTag, ElementBoundaryTypes);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
enum {
@@ -67,127 +76,128 @@ protected:
nCompIdx = Indices::nCompIdx,
};
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
-
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Element = typename GridView::template Codim<0>::Entity;
-
//! property that defines whether mole or mass fractions are used
static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
public:
+ using ParentType::ParentType;
/*!
- * \brief Evaluate the storage term of the current solution in a
- * single phase.
+ * \brief Evaluate the amount of all conservation quantities
+ * (e.g. phase mass) within a sub-control volume.
*
- * \param element The element
- * \param phaseIdx The index of the fluid phase
+ * 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 evalPhaseStorage(const Element &element, const int phaseIdx)
+ ResidualVector computeStorage(const Problem& problem,
+ const SubControlVolume& scv,
+ const VolumeVariables& volVars) const
{
- FVElementGeometry fvGeometry;
- fvGeometry.update(this->gridView_(), element);
- ElementBoundaryTypes bcTypes;
- bcTypes.update(this->problem_(), element, fvGeometry);
- ElementVolumeVariables elemVolVars;
- elemVolVars.update(this->problem_(), element, fvGeometry, false);
-
- this->storageTerm_.resize(fvGeometry.numScv);
- this->storageTerm_ = 0;
-
- this->elemPtr_ = &element;
- this->fvElemGeomPtr_ = &fvGeometry;
- this->bcTypesPtr_ = &bcTypes;
- this->prevVolVarsPtr_ = 0;
- this->curVolVarsPtr_ = &elemVolVars;
- evalPhaseStorage_(phaseIdx);
+ ResidualVector storage(0.0);
+
+ const auto massOrMoleDensity = [](const auto& volVars, const int phaseIdx)
+ { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
+
+ const auto massOrMoleFraction= [](const auto& volVars, const int phaseIdx, const int compIdx)
+ { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
+
+ // compute storage term of all components within all phases
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ int eqIdx = (compIdx == wCompIdx) ? conti0EqIdx : conti1EqIdx;
+ storage[eqIdx] += volVars.porosity()
+ * volVars.saturation(phaseIdx)
+ * massOrMoleDensity(volVars, phaseIdx)
+ * massOrMoleFraction(volVars, phaseIdx, compIdx);
+ }
+
+ //! The energy storage in the fluid phase with index phaseIdx
+ EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, phaseIdx);
+ }
+
+ //! The energy storage in the solid matrix
+ EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
+
+ return storage;
}
- /*!
- * \brief Adds the diffusive mass flux of all components over
- * a face of a subcontrol volume.
+ /*!
+ * \brief Evaluates the total flux of all conservation quantities
+ * over a face of a sub-control volume.
*
- * \param flux The diffusive flux over the sub-control-volume face for each component
- * \param fluxVars The flux variables at the current SCV
+ * \param flux The flux over the SCV (sub-control-volume) face for each component
+ * \param fIdx 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 computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
+ ResidualVector computeFlux(const Problem& problem,
+ const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf,
+ const ElementFluxVariablesCache& elemFluxVarsCache) const
{
- // TODO: reference!? Dune::FieldMatrix<Scalar, numPhases, numComponents> averagedPorousDiffCoeffMatrix = fluxVars.porousDiffCoeff();
- // add diffusive flux of gas component in liquid phase
- Scalar tmp;
- tmp = - fluxVars.porousDiffCoeff(wPhaseIdx) * fluxVars.molarDensity(wPhaseIdx);
- tmp *= (fluxVars.moleFractionCompNGrad(wPhaseIdx) * fluxVars.face().normal);
- Scalar jNW = tmp;
+ FluxVariables fluxVars;
+ fluxVars.init(problem, element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
+
+ // get upwind weights into local scope
+ ResidualVector flux(0.0);
+
+ const auto massOrMoleDensity = [](const auto& volVars, const int phaseIdx)
+ { return useMoles ? volVars.molarDensity(phaseIdx) : volVars.density(phaseIdx); };
+
+ const auto massOrMoleFraction= [](const auto& volVars, const int phaseIdx, const int compIdx)
+ { return useMoles ? volVars.moleFraction(phaseIdx, compIdx) : volVars.massFraction(phaseIdx, compIdx); };
- Scalar jWW = -jNW;
+ // advective fluxes
+ for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
+ {
+ for (int compIdx = 0; compIdx < numComponents; ++compIdx)
+ {
+ const auto upwindTerm = [&massOrMoleDensity, &massOrMoleFraction, phaseIdx, compIdx] (const auto& volVars)
+ { return massOrMoleDensity(volVars, phaseIdx)*massOrMoleFraction(volVars, phaseIdx, compIdx)*volVars.mobility(phaseIdx); };
+
+ // get equation index
+ auto eqIdx = conti0EqIdx + compIdx;
+ flux[eqIdx] += fluxVars.advectiveFlux(phaseIdx, upwindTerm);
+ }
+
+ //! Add advective phase energy fluxes. For isothermal model the contribution is zero.
+ EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, phaseIdx);
+ }
- tmp = - fluxVars.porousDiffCoeff(gPhaseIdx) * fluxVars.molarDensity(gPhaseIdx);
- tmp *= (fluxVars.moleFractionCompWGrad(gPhaseIdx) * fluxVars.face().normal);
- Scalar jWG = tmp;
+ //! Add diffusive energy fluxes. For isothermal model the contribution is zero.
+ EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
- Scalar jNG = -jWG;
+ const auto diffusionFluxesWPhase = fluxVars.molecularDiffusionFlux(wPhaseIdx);
- tmp = - fluxVars.porousDiffCoeff(nPhaseIdx) * fluxVars.molarDensity(nPhaseIdx);
- tmp *= (fluxVars.moleFractionCompWGrad(nPhaseIdx) * fluxVars.face().normal);
- Scalar jWN = tmp;
+ // diffusive fluxes
+ Scalar jNW = diffusionFluxesWPhase[nCompIdx];
+ Scalar jWW = -(jNW);
+ const auto diffusionFluxesGPhase = fluxVars.molecularDiffusionFlux(gPhaseIdx);
+
+ Scalar jWG = diffusionFluxesGPhase[wCompIdx];
+ Scalar jNG = -(jWG);
+
+ const auto diffusionFluxesNPhase = fluxVars.molecularDiffusionFlux(nPhaseIdx);
+ // At the moment we do not consider diffusion in the NAPL phase
+ Scalar jWN = diffusionFluxesNPhase[wCompIdx];
Scalar jNN = -jWN;
flux[conti0EqIdx] += jWW+jWG+jWN;
flux[conti1EqIdx] += jNW+jNG+jNN;
- }
-protected:
- void evalPhaseStorage_(const int phaseIdx)
- {
- if(!useMoles) //mass-fraction formulation
- {
- // evaluate the storage terms of a single phase
- for (int i=0; i < this->fvGeometry_().numScv; i++) {
- PrimaryVariables &storage = this->storageTerm_[i];
- const ElementVolumeVariables &elemVolVars = this->curVolVars_();
- const VolumeVariables &volVars = elemVolVars[i];
-
- // compute storage term of all components within all phases
- storage = 0;
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- int eqIdx = (compIdx == wCompIdx) ? conti0EqIdx : conti1EqIdx;
- storage[eqIdx] += volVars.density(phaseIdx)
- * volVars.saturation(phaseIdx)
- * volVars.fluidState().massFraction(phaseIdx, compIdx);
- }
+ return flux;
+ }
- storage *= volVars.porosity();
- storage *= this->fvGeometry_().subContVol[i].volume;
- }
- }
- else //mole-fraction formulation
- {
- // evaluate the storage terms of a single phase
- for (int i=0; i < this->fvGeometry_().numScv; i++) {
- PrimaryVariables &storage = this->storageTerm_[i];
- const ElementVolumeVariables &elemVolVars = this->curVolVars_();
- const VolumeVariables &volVars = elemVolVars[i];
-
- // compute storage term of all components within all phases
- storage = 0;
- for (int compIdx = 0; compIdx < numComponents; ++compIdx)
- {
- int eqIdx = (compIdx == wCompIdx) ? conti0EqIdx : conti1EqIdx;
- storage[eqIdx] += volVars.molarDensity(phaseIdx)
- * volVars.saturation(phaseIdx)
- * volVars.fluidState().moleFraction(phaseIdx, compIdx);
- }
- storage *= volVars.porosity();
- storage *= this->fvGeometry_().subContVol[i].volume;
- }
- }
- }
+protected:
Implementation *asImp_()
{
return static_cast<Implementation *> (this);
diff --git a/dumux/porousmediumflow/3pwateroil/model.hh b/dumux/porousmediumflow/3pwateroil/model.hh
index bbde1e52f266d3772bfda8bbac406dd9dfe93ad0..0506938e8876fecf6420394cbc39591a04281179 100644
--- a/dumux/porousmediumflow/3pwateroil/model.hh
+++ b/dumux/porousmediumflow/3pwateroil/model.hh
@@ -18,23 +18,11 @@
*****************************************************************************/
/*!
* \file
- *
+ * \ingroup ThreePWaterOilModel
* \brief Adaption of the fully implicit scheme to the three-phase three-component flow model.
*
* The model is designed for simulating three fluid phases with water, gas, and
* a liquid contaminant (NAPL - non-aqueous phase liquid)
- */
-#ifndef DUMUX_3PWATEROIL_MODEL_HH
-#define DUMUX_3PWATEROIL_MODEL_HH
-
-#include "properties.hh"
-
-namespace Dumux
-{
-/*!
- * \ingroup ThreePWaterOilModel
- * \brief Adaption of the fully implicit scheme to the three-phase two-component flow model.
- *
* This model implements three-phase two-component flow of three fluid phases
* \f$\alpha \in \{ water, gas, NAPL \}\f$ each composed of up to two components
* \f$\kappa \in \{ water, contaminant \}\f$. The standard multiphase Darcy
@@ -79,996 +67,161 @@ namespace Dumux
* <li> ... to be completed ... </li>
* </ul>
*/
-template<class TypeTag>
-class ThreePWaterOilModel: public GET_PROP_TYPE(TypeTag, BaseModel)
-{
- using ParentType = typename GET_PROP_TYPE(TypeTag, BaseModel);
-
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
- using Indices = typename GET_PROP_TYPE(TypeTag, 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,
- switch2Idx = Indices::switch2Idx,
-
-
- wPhaseIdx = Indices::wPhaseIdx,
- nPhaseIdx = Indices::nPhaseIdx,
- gPhaseIdx = Indices::gPhaseIdx,
-
- wCompIdx = Indices::wCompIdx,
- nCompIdx = Indices::nCompIdx,
-
- threePhases = Indices::threePhases,
- wPhaseOnly = Indices::wPhaseOnly,
- gnPhaseOnly = Indices::gnPhaseOnly,
- wnPhaseOnly = Indices::wnPhaseOnly,
- gPhaseOnly = Indices::gPhaseOnly,
- wgPhaseOnly = Indices::wgPhaseOnly
-
- };
-
-
- using Vertex = typename GridView::template Codim<dim>::Entity;
- using Element = typename GridView::template Codim<0>::Entity;
-
- using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
-
- static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
-
-
- 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 Compute the total storage inside one phase of all
- * conservation quantities.
- *
- * \param storage Contains the storage of each component for one phase
- * \param phaseIdx The phase index
- */
- void globalPhaseStorage(PrimaryVariables &storage, const int phaseIdx)
- {
- storage = 0;
-
- for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
- {
- this->localResidual().evalPhaseStorage(element, phaseIdx);
-
- for (unsigned int i = 0; i < this->localResidual().storageTerm().size(); ++i)
- storage += this->localResidual().storageTerm()[i];
- }
- if (this->gridView_().comm().size() > 1)
- storage = this->gridView_().comm().sum(storage);
- }
-
-
- /*!
- * \brief Called by the update() method if applying the newton
- * method was unsuccessful.
- */
- void updateFailed()
- {
- ParentType::updateFailed();
-
- setSwitched_(false);
- resetPhasePresence_();
- };
-
- /*!
- * \brief Called by the problem if a time integration was
- * successful, post processing of the solution is done and the
- * result has been written to disk.
- *
- * This should prepare the model for the next time integration.
- */
- void advanceTimeLevel()
- {
- ParentType::advanceTimeLevel();
-
- // update the phase state
- updateOldPhasePresence_();
- setSwitched_(false);
- }
-
- /*!
- * \brief Return true if the primary variables were switched for
- * at least one vertex after the last timestep.
- */
- bool switched() const
- {
- return switchFlag_;
- }
-
- /*!
- * \brief Returns the phase presence of the current or the old solution of a 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)
- {
- using ScalarField = Dune::BlockVector< Dune::FieldVector<double, 1> >;
- using VectorField = Dune::BlockVector<Dune::FieldVector<double, dim> >;
-
- // 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 *moleFraction[numPhases][numComponents];
- for (int i = 0; i < numPhases; ++i)
- for (int j = 0; j < numComponents; ++j)
- moleFraction[i][j] = writer.allocateManagedBuffer (numDofs);
- ScalarField *temperature = writer.allocateManagedBuffer (numDofs);
- ScalarField *poro = writer.allocateManagedBuffer(numDofs);
- ScalarField *perm = writer.allocateManagedBuffer(numDofs);
- VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
- VectorField *velocityW = writer.template allocateManagedBuffer<double, dim>(numDofs);
- VectorField *velocityG = writer.template allocateManagedBuffer<double, dim>(numDofs);
- ImplicitVelocityOutput<TypeTag> velocityOutput(this->problem_());
-
- if (velocityOutput.enableOutput()) // check if velocity output is demanded
- {
- // initialize velocity fields
- for (unsigned int i = 0; i < numDofs; ++i)
- {
- (*velocityN)[i] = Scalar(0);
- (*velocityW)[i] = Scalar(0);
- (*velocityG)[i] = Scalar(0);
- }
- }
-
- unsigned numElements = this->gridView_().size(0);
- ScalarField *rank = writer.allocateManagedBuffer (numElements);
-
- for (const auto& element : elements(this->gridView_(), Dune::Partitions::interior))
- {
- int idx = this->problem_().elementMapper().index(element);
- (*rank)[idx] = 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);
- }
-
- for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
- for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
- (*moleFraction[phaseIdx][compIdx])[globalIdx] =
- elemVolVars[scvIdx].fluidState().moleFraction(phaseIdx,
- compIdx);
-
- Valgrind::CheckDefined((*moleFraction[phaseIdx][compIdx])[globalIdx]);
- }
- }
-
- (*poro)[globalIdx] = elemVolVars[scvIdx].porosity();
- (*perm)[globalIdx] = elemVolVars[scvIdx].permeability();
- (*temperature)[globalIdx] = elemVolVars[scvIdx].temperature();
- (*phasePresence)[globalIdx] = staticDat_[globalIdx].phasePresence;
- }
-
- // velocity output
- velocityOutput.calculateVelocity(*velocityW, elemVolVars, fvGeometry, element, wPhaseIdx);
- velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, nPhaseIdx);
- velocityOutput.calculateVelocity(*velocityN, elemVolVars, fvGeometry, element, gPhaseIdx);
-
- }
-
- writer.attachDofData(*saturation[wPhaseIdx], "sw", isBox);
- writer.attachDofData(*saturation[nPhaseIdx], "sn", isBox);
- writer.attachDofData(*saturation[gPhaseIdx], "sg", isBox);
- writer.attachDofData(*pressure[wPhaseIdx], "pw", isBox);
- writer.attachDofData(*pressure[nPhaseIdx], "pn", isBox);
- writer.attachDofData(*pressure[gPhaseIdx], "pg", isBox);
- writer.attachDofData(*density[wPhaseIdx], "rhow", isBox);
- writer.attachDofData(*density[nPhaseIdx], "rhon", isBox);
- writer.attachDofData(*density[gPhaseIdx], "rhog", isBox);
-
- writer.attachDofData(*mobility[wPhaseIdx], "MobW", isBox);
- writer.attachDofData(*mobility[nPhaseIdx], "MobN", isBox);
- writer.attachDofData(*mobility[gPhaseIdx], "MobG", isBox);
-
- writer.attachDofData(*viscosity[wPhaseIdx], "ViscosW", isBox);
- writer.attachDofData(*viscosity[nPhaseIdx], "ViscosN", isBox);
- writer.attachDofData(*viscosity[gPhaseIdx], "ViscosG", isBox);
- for (int i = 0; i < numPhases; ++i)
- {
- for (int j = 0; j < numComponents; ++j)
- {
- std::ostringstream oss;
- oss << "x^"
- << FluidSystem::phaseName(i)
- << "_"
- << FluidSystem::componentName(j);
- writer.attachDofData(*moleFraction[i][j], oss.str().c_str(), isBox);
- }
- }
- writer.attachDofData(*poro, "porosity", isBox);
- writer.attachDofData(*perm, "permeability", isBox);
- writer.attachDofData(*temperature, "temperature", isBox);
- writer.attachDofData(*phasePresence, "phase presence", isBox);
-
- if (velocityOutput.enableOutput()) // check if velocity output is demanded
- {
- writer.attachDofData(*velocityW, "velocityW", isBox, dim);
- writer.attachDofData(*velocityN, "velocityN", isBox, dim);
- writer.attachDofData(*velocityG, "velocityG", isBox, dim);
- }
-
- writer.attachCellData(*rank, "process rank");
- }
-
- /*!
- * \brief Write the current solution to a restart file.
- *
- * \param outStream The output stream of one entity 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 globalIdx = this->dofMapper().index(entity);
- if (!outStream.good())
- DUNE_THROW(Dune::IOError, "Could not serialize entity " << globalIdx);
-
- outStream << staticDat_[globalIdx].phasePresence << " ";
- }
-
- /*!
- * \brief Reads the current solution from a restart file.
- *
- * \param inStream The input stream of one entity 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 globalIdx = this->dofMapper().index(entity);
- if (!inStream.good())
- DUNE_THROW(Dune::IOError,
- "Could not deserialize entity " << globalIdx);
-
- inStream >> staticDat_[globalIdx].phasePresence;
- staticDat_[globalIdx].oldPhasePresence
- = staticDat_[globalIdx].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;
-
- bool useSimpleModel = GET_PROP_VALUE(TypeTag, UseSimpleModel);
-
- 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(useSimpleModel)
- {
- if (primaryVarSwitchSimple_(curGlobalSol,
- volVars,
- globalIdx,
- globalPos))
- {
- this->jacobianAssembler().markDofRed(globalIdx);
- wasSwitched = true;
- }
- }
- else
- {
- if (primaryVarSwitch_(curGlobalSol,
- volVars,
- globalIdx,
- globalPos))
- {
- this->jacobianAssembler().markDofRed(globalIdx);
- wasSwitched = true;
- }
- }
- }
- }
-
- // 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;
-
- // check if a primary var switch is necessary
- if (phasePresence == threePhases)
- {
- 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 = wnPhaseOnly;
-
- globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
- globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
- }
- else if (volVars.saturation(wPhaseIdx) <= Smin)
- {
- wouldSwitch = true;
- // water phase disappears
- std::cout << "Water phase disappears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", sw: "
- << volVars.saturation(wPhaseIdx) << std::endl;
- newPhasePresence = gnPhaseOnly;
-
- globalSol[globalIdx][switch1Idx] = volVars.fluidState().saturation(nPhaseIdx);
- globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
- }
- else if (volVars.saturation(nPhaseIdx) <= Smin)
- {
- wouldSwitch = true;
- // NAPL phase disappears
- std::cout << "NAPL phase disappears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", sn: "
- << volVars.saturation(nPhaseIdx) << std::endl;
- newPhasePresence = wgPhaseOnly;
-
- globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
- }
- }
- else if (phasePresence == wPhaseOnly)
- {
- bool gasFlag = 0;
- bool nonwettingFlag = 0;
- // calculate fractions of the partial pressures in the
- // hypothetical gas phase
- Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
- Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
-
- Scalar xgMax = 1.0;
- if (xwg + xng > xgMax)
- wouldSwitch = true;
- if (staticDat_[globalIdx].wasSwitched)
- xgMax *= 1.02;
-
- // if the sum of the mole fractions would be larger than
- // 100%, gas phase appears
- if (xwg + xng > xgMax)
- {
- // gas phase appears
- std::cout << "gas phase appears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", xwg + xng: "
- << xwg + xng << std::endl;
- gasFlag = 1;
- }
-
- // calculate fractions in the hypothetical NAPL phase
- Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
- /* take care:
- for xnn in case wPhaseOnly we compute xnn=henry_mesitylene*x1w,
- where a hypothetical gas pressure is assumed for the Henry
- x0n is set to NULL (all NAPL phase is dirty)
- x2n is set to NULL (all NAPL phase is dirty)
- */
-
- Scalar xnMax = 1.0;
- if (xnn > xnMax)
- wouldSwitch = true;
- if (staticDat_[globalIdx].wasSwitched)
- xnMax *= 1.02;
-
- // if the sum of the hypothetical mole fractions would be larger than
- // 100%, NAPL phase appears
- if (xnn > xnMax)
- {
- // NAPL phase appears
- std::cout << "NAPL phase appears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", xnn: "
- << xnn << std::endl;
- nonwettingFlag = 1;
- }
-
- if ((gasFlag == 1) && (nonwettingFlag == 0))
- {
- newPhasePresence = wgPhaseOnly;
- globalSol[globalIdx][switch1Idx] = 0.9999;
- globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(gPhaseIdx);
- }
- else if ((gasFlag == 1) && (nonwettingFlag == 1))
- {
- newPhasePresence = threePhases;
- globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(gPhaseIdx);
- globalSol[globalIdx][switch1Idx] = 0.999;
- }
- else if ((gasFlag == 0) && (nonwettingFlag == 1))
- {
- newPhasePresence = wnPhaseOnly;
- globalSol[globalIdx][switch2Idx] = 0.0001;
- }
- }
- else if (phasePresence == gnPhaseOnly)
- {
- bool nonwettingFlag = 0;
- bool wettingFlag = 0;
-
- Scalar Smin = 0.0;
- if (staticDat_[globalIdx].wasSwitched)
- Smin = -0.01;
-
- if (volVars.saturation(nPhaseIdx) <= Smin)
- {
- wouldSwitch = true;
- // NAPL phase disappears
- std::cout << "NAPL phase disappears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", sn: "
- << volVars.saturation(nPhaseIdx) << std::endl;
- nonwettingFlag = 1;
- }
+#ifndef DUMUX_3PWATEROIL_MODEL_HH
+#define DUMUX_3PWATEROIL_MODEL_HH
- // calculate fractions of the hypothetical water phase
- Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
- /*
- take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
- If this is larger than 1, then water appears
- */
- Scalar xwMax = 1.0;
- if (xww > xwMax)
- wouldSwitch = true;
- if (staticDat_[globalIdx].wasSwitched)
- xwMax *= 1.02;
+#include <dumux/common/properties.hh>
- // if the sum of the mole fractions would be larger than
- // 100%, water phase appears
- if (xww > xwMax)
- {
- // water phase appears
- std::cout << "water phase appears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
- << xww << std::endl;
- wettingFlag = 1;
- }
+#include <dumux/material/spatialparams/fv.hh>
+#include <dumux/material/fluidmatrixinteractions/3p/thermalconductivitysomerton3p.hh>
+#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
- if ((wettingFlag == 1) && (nonwettingFlag == 0))
- {
- newPhasePresence = threePhases;
- globalSol[globalIdx][switch1Idx] = 0.0001;
- globalSol[globalIdx][switch2Idx] = volVars.saturation(nPhaseIdx);
- }
- else if ((wettingFlag == 1) && (nonwettingFlag == 1))
- {
- newPhasePresence = wgPhaseOnly;
- globalSol[globalIdx][switch1Idx] = 0.0001;
- globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
- }
- else if ((wettingFlag == 0) && (nonwettingFlag == 1))
- {
- newPhasePresence = gPhaseOnly;
- globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
- globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
- }
- }
- else if (phasePresence == wnPhaseOnly)
- {
- bool nonwettingFlag = 0;
- bool gasFlag = 0;
+#include <dumux/porousmediumflow/properties.hh>
+#include <dumux/porousmediumflow/nonisothermal/model.hh>
+#include <dumux/porousmediumflow/compositional/switchableprimaryvariables.hh>
- Scalar Smin = 0.0;
- if (staticDat_[globalIdx].wasSwitched)
- Smin = -0.01;
+#include "indices.hh"
+#include "model.hh"
+#include "volumevariables.hh"
+#include "localresidual.hh"
+#include "primaryvariableswitch.hh"
+#include "vtkoutputfields.hh"
- if (volVars.saturation(nPhaseIdx) <= Smin)
- {
- wouldSwitch = true;
- // NAPL phase disappears
- std::cout << "NAPL phase disappears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", sn: "
- << volVars.saturation(nPhaseIdx) << std::endl;
- nonwettingFlag = 1;
- }
+namespace Dumux
+{
+/*!
+ * \ingroup ThreePWaterOilModel
+ * \brief Adaption of the fully implicit scheme to the three-phase two-component flow model.
+ *
+ */
+namespace Properties {
- // calculate fractions of the partial pressures in the
- // hypothetical gas phase
- Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
- Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+NEW_TYPE_TAG(ThreePWaterOilNI, INHERITS_FROM(PorousMediumFlow, NonIsothermal));
- Scalar xgMax = 1.0;
- if (xwg + xng > xgMax)
- wouldSwitch = true;
- if (staticDat_[globalIdx].wasSwitched)
- xgMax *= 1.02;
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
- // if the sum of the mole fractions would be larger than
- // 100%, gas phase appears
- if (xwg + xng > xgMax)
- {
- // gas phase appears
- std::cout << "gas phase appears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", xwg + xng: "
- << xwg + xng << std::endl;
- gasFlag = 1;
- }
+/*!
+ * \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 2.
+ */
+SET_PROP(ThreePWaterOilNI, NumComponents)
+{
+ private:
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- if ((gasFlag == 1) && (nonwettingFlag == 0))
- {
- newPhasePresence = threePhases;
- globalSol[globalIdx][pressureIdx] = volVars.pressure(gPhaseIdx);
- globalSol[globalIdx][switch1Idx] = volVars.saturation(wPhaseIdx);
- }
- else if ((gasFlag == 1) && (nonwettingFlag == 1))
- {
- newPhasePresence = wgPhaseOnly;
- globalSol[globalIdx][pressureIdx] = volVars.pressure(gPhaseIdx);
- globalSol[globalIdx][switch1Idx] = volVars.temperature();
- globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
- }
- else if ((gasFlag == 0) && (nonwettingFlag == 1))
- {
- newPhasePresence = wPhaseOnly;
- globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
- }
- }
- else if (phasePresence == gPhaseOnly)
- {
- bool nonwettingFlag = 0;
- bool wettingFlag = 0;
+ public:
+ static const int value = FluidSystem::numComponents;
- // calculate fractions in the hypothetical NAPL phase
- Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
- /*
- take care:, xnn, if no NAPL phase is there, take xnn=xng*pg/pcsat
- if this is larger than 1, then NAPL appears
- */
+ static_assert(value == 2,
+ "Only fluid systems with 2 components are supported by the 3p2cni model!");
+};
- Scalar xnMax = 1.0;
- if (xnn > xnMax)
- wouldSwitch = true;
- if (staticDat_[globalIdx].wasSwitched)
- xnMax *= 1.02;
+/*!
+ * \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(ThreePWaterOilNI, NumPhases)
+{
+ private:
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- // if the sum of the hypothetical mole fraction would be larger than
- // 100%, NAPL phase appears
- if (xnn > xnMax)
- {
- // NAPL phase appears
- std::cout << "NAPL phase appears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", xnn: "
- << xnn << std::endl;
- nonwettingFlag = 1;
- }
- // calculate fractions of the hypothetical water phase
- Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
- /*
- take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
- If this is larger than 1, then water appears
- */
- Scalar xwMax = 1.0;
- if (xww > xwMax)
- wouldSwitch = true;
- if (staticDat_[globalIdx].wasSwitched)
- xwMax *= 1.02;
+ public:
+ static const int value = FluidSystem::numPhases;
+ static_assert(value == 3,
+ "Only fluid systems with 3 phases are supported by the 3p2cni model!");
+};
- // if the sum of the mole fractions would be larger than
- // 100%, gas phase appears
- if (xww > xwMax)
- {
- // water phase appears
- std::cout << "water phase appears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
- << xww << std::endl;
- wettingFlag = 1;
- }
- if ((wettingFlag == 1) && (nonwettingFlag == 0))
- {
- newPhasePresence = wgPhaseOnly;
- globalSol[globalIdx][switch1Idx] = 0.0001;
- globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
- }
- else if ((wettingFlag == 1) && (nonwettingFlag == 1))
- {
- newPhasePresence = threePhases;
- globalSol[globalIdx][switch1Idx] = 0.0001;
- globalSol[globalIdx][switch2Idx] = 0.0001;
- }
- else if ((wettingFlag == 0) && (nonwettingFlag == 1))
- {
- newPhasePresence = gnPhaseOnly;
- globalSol[globalIdx][switch2Idx]
- = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
- }
- }
- else if (phasePresence == wgPhaseOnly)
- {
- bool nonwettingFlag = 0;
- bool gasFlag = 0;
- bool wettingFlag = 0;
+/*!
+ * \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(ThreePWaterOilNI, FluidState){
+ private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ public:
+ using type = CompositionalFluidState<Scalar, FluidSystem>;
+};
- // get the fractions in the hypothetical NAPL phase
- Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
+//! The local residual function of the conservation equations
+SET_TYPE_PROP(ThreePWaterOilNI, LocalResidual, ThreePWaterOilLocalResidual<TypeTag>);
- // take care: if the NAPL phase is not present, take
- // xnn=xng*pg/pcsat if this is larger than 1, then NAPL
- // appears
- Scalar xnMax = 1.0;
- if (xnn > xnMax)
- wouldSwitch = true;
- if (staticDat_[globalIdx].wasSwitched)
- xnMax *= 1.02;
+//! Set as default that no component mass balance is replaced by the total mass balance
+SET_INT_PROP(ThreePWaterOilNI, ReplaceCompEqIdx, GET_PROP_VALUE(TypeTag, NumComponents));
- // if the sum of the hypothetical mole fraction would be larger than
- // 100%, NAPL phase appears
- if (xnn > xnMax)
- {
- // NAPL phase appears
- std::cout << "NAPL phase appears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", xnn: "
- << xnn << std::endl;
- nonwettingFlag = 1;
- }
+//! Enable advection
+SET_BOOL_PROP(ThreePWaterOilNI, EnableAdvection, true);
- Scalar Smin = -1.e-6;
- if (staticDat_[globalIdx].wasSwitched)
- Smin = -0.01;
+//! disable molecular diffusion for the 3p model
+SET_BOOL_PROP(ThreePWaterOilNI, EnableMolecularDiffusion, true);
- 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;
- gasFlag = 1;
- }
+//! Isothermal model by default
+SET_BOOL_PROP(ThreePWaterOilNI, EnableEnergyBalance, true);
- Smin = 0.0;
- if (staticDat_[globalIdx].wasSwitched)
- Smin = -0.01;
+//! The primary variable switch for the 3p3c model
+SET_TYPE_PROP(ThreePWaterOilNI, PrimaryVariableSwitch, ThreePWaterOilPrimaryVariableSwitch<TypeTag>);
- if (volVars.saturation(wPhaseIdx) <= Smin)
- {
- wouldSwitch = true;
- // gas phase disappears
- std::cout << "Water phase disappears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", sw: "
- << volVars.saturation(wPhaseIdx) << std::endl;
- wettingFlag = 1;
- }
+//! The primary variables vector for the 3p3c model
+SET_TYPE_PROP(ThreePWaterOilNI, PrimaryVariables, SwitchablePrimaryVariables<TypeTag, int>);
- if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 1))
- {
- newPhasePresence = gnPhaseOnly;
- globalSol[globalIdx][switch1Idx] = 0.0001;
- globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
-;
- }
- else if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 0))
- {
- newPhasePresence = threePhases;
- globalSol[globalIdx][switch2Idx] = 0.0001;
- }
- else if ((gasFlag == 1) && (nonwettingFlag == 0) && (wettingFlag == 0))
- {
- newPhasePresence = wPhaseOnly;
- globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
- globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
- globalSol[globalIdx][switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
- }
- else if ((gasFlag == 0) && (nonwettingFlag == 0) && (wettingFlag == 1))
- {
- newPhasePresence = gPhaseOnly;
- globalSol[globalIdx][switch1Idx]
- = volVars.fluidState().temperature();
- globalSol[globalIdx][switch2Idx]
- = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
- }
- }
+//! Determines whether a constraint solver should be used explicitly
+SET_BOOL_PROP(ThreePWaterOilNI, OnlyGasPhaseCanDisappear, true);
- staticDat_[globalIdx].phasePresence = newPhasePresence;
- staticDat_[globalIdx].wasSwitched = wouldSwitch;
- return phasePresence != newPhasePresence;
- }
+//! the VolumeVariables property
+SET_TYPE_PROP(ThreePWaterOilNI, VolumeVariables, ThreePWaterOilVolumeVariables<TypeTag>);
- // perform variable switch at a vertex; Returns true if a
- // variable switch was performed.
- // Note that this one is the simplified version with only two phase states
- bool primaryVarSwitchSimple_(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;
+//! The spatial parameters to be employed.
+//! Use FVSpatialParams by default.
+SET_TYPE_PROP(ThreePWaterOilNI, SpatialParams, FVSpatialParams<TypeTag>);
- // check if a primary var switch is necessary
- if (phasePresence == threePhases)
- {
- Scalar Smin = 0;
- if (staticDat_[globalIdx].wasSwitched)
- Smin = -0.01;
+//! Use the model after Millington (1961) for the effective diffusivity
+SET_TYPE_PROP(ThreePWaterOilNI, EffectiveDiffusivityModel,
+ DiffusivityMillingtonQuirk<typename GET_PROP_TYPE(TypeTag, Scalar)>);
- 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 = wnPhaseOnly;
+// disable velocity output by default
- globalSol[globalIdx][pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
- globalSol[globalIdx][switch1Idx] = volVars.fluidState().temperature();
- }
- }
- else if (phasePresence == wnPhaseOnly)
- {
- bool gasFlag = 0;
+SET_BOOL_PROP(ThreePWaterOilNI, UseMoles, true); //!< Define that mole fractions are used in the balance equations per default
- // calculate fractions of the partial pressures in the
- // hypothetical gas phase
- Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
- Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+//! Somerton is used as default model to compute the effective thermal heat conductivity
+SET_PROP(ThreePWaterOilNI, ThermalConductivityModel)
+{
+private:
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+public:
+ using type = ThermalConductivitySomerton<Scalar, Indices>;
+};
- Scalar xgMax = 1.0;
- if (xwg + xng > xgMax)
- wouldSwitch = true;
- if (staticDat_[globalIdx].wasSwitched)
- xgMax *= 1.02;
- // if the sum of the mole fractions would be larger than
- // 100%, gas phase appears
- if (xwg + xng > xgMax)
- {
- // gas phase appears
- std::cout << "gas phase appears at vertex " << globalIdx
- << ", coordinates: " << globalPos << ", xwg + xng: "
- << xwg + xng << std::endl;
- gasFlag = 1;
- }
+//////////////////////////////////////////////////////////////////
+// Property values for isothermal model required for the general non-isothermal model
+//////////////////////////////////////////////////////////////////
- if (gasFlag == 1)
- {
- newPhasePresence = threePhases;
- globalSol[globalIdx][pressureIdx] = volVars.pressure(gPhaseIdx);
- globalSol[globalIdx][switch1Idx] = volVars.saturation(wPhaseIdx);
- }
- }
+//set isothermal output fields
+SET_TYPE_PROP(ThreePWaterOilNI, IsothermalVtkOutputFields, ThreePWaterOilVtkOutputFields<TypeTag>);
- staticDat_[globalIdx].phasePresence = newPhasePresence;
- staticDat_[globalIdx].wasSwitched = wouldSwitch;
- return phasePresence != newPhasePresence;
- }
+//set isothermal Indices
+SET_PROP(ThreePWaterOilNI, IsothermalIndices)
+{
- // parameters given in constructor
- std::vector<StaticVars> staticDat_;
- bool switchFlag_;
+public:
+ using type = ThreePWaterOilIndices<TypeTag, /*PVOffset=*/0>;
};
+//set isothermal NumEq
+SET_INT_PROP(ThreePWaterOilNI, IsothermalNumEq, 2);
+
}
-#include "propertydefaults.hh"
+ }
#endif
diff --git a/dumux/porousmediumflow/3pwateroil/newtoncontroller.hh b/dumux/porousmediumflow/3pwateroil/newtoncontroller.hh
deleted file mode 100644
index 419cf1d45b404d4cbd84c53ae4ac8691fcb444db..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/3pwateroil/newtoncontroller.hh
+++ /dev/null
@@ -1,85 +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 3p2cni 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_3P2CNI_NEWTON_CONTROLLER_HH
-#define DUMUX_3P2CNI_NEWTON_CONTROLLER_HH
-
-#include <dumux/nonlinear/newtoncontroller.hh>
-
-#include "properties.hh"
-
-namespace Dumux {
-/*!
- * \ingroup ThreePWaterOilModel
- * \brief A 3p2cni 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 ThreePWaterOilNewtonController : public NewtonController<TypeTag>
-{
- using ParentType = NewtonController<TypeTag>;
- using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
-
-public:
- ThreePWaterOilNewtonController(const Problem &problem)
- : ParentType(problem)
- {};
-
-
- /*!
- * \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
- this->method().model().updateStaticData(uCurrentIter, uLastIter);
- ParentType::newtonEndStep(uCurrentIter, uLastIter);
- }
-
-
- /*!
- * \brief Returns true if the current solution can be considered to
- * be accurate enough
- */
- bool newtonConverged()
- {
- if (this->method().model().switched())
- return false;
-
- return ParentType::newtonConverged();
- };
-};
-}
-
-#endif
diff --git a/dumux/porousmediumflow/3pwateroil/primaryvariableswitch.hh b/dumux/porousmediumflow/3pwateroil/primaryvariableswitch.hh
new file mode 100644
index 0000000000000000000000000000000000000000..35c945abefb71c7a6d05ca6f7a334d1679dd5d87
--- /dev/null
+++ b/dumux/porousmediumflow/3pwateroil/primaryvariableswitch.hh
@@ -0,0 +1,549 @@
+// -*- 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
+ * \ingroup ThreePWaterOilModel
+ * \brief The primary variable switch for the 3p3c model
+ */
+#ifndef DUMUX_3P3C_PRIMARY_VARIABLE_SWITCH_HH
+#define DUMUX_3P3C_PRIMARY_VARIABLE_SWITCH_HH
+
+#include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh>
+
+namespace Dumux
+{
+/*!
+ * \ingroup ThreePWaterOilModel
+ * \brief The primary variable switch controlling the phase presence state variable
+ */
+template<class TypeTag>
+class ThreePWaterOilPrimaryVariableSwitch : public Dumux::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 {
+ switch1Idx = Indices::switch1Idx,
+ switch2Idx = Indices::switch2Idx,
+ pressureIdx = Indices::pressureIdx,
+
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ gPhaseIdx = Indices::gPhaseIdx,
+
+ wCompIdx = Indices::wCompIdx,
+ nCompIdx = Indices::nCompIdx,
+ gCompIdx = Indices::gCompIdx,
+
+ threePhases = Indices::threePhases,
+ wPhaseOnly = Indices::wPhaseOnly,
+ gnPhaseOnly = Indices::gnPhaseOnly,
+ wnPhaseOnly = Indices::wnPhaseOnly,
+ gPhaseOnly = Indices::gPhaseOnly,
+ wgPhaseOnly = Indices::wgPhaseOnly
+ };
+
+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;
+ auto phasePresence = priVars.state();
+ int newPhasePresence = phasePresence;
+
+ bool onlyGasPhaseCanDisappear = GET_PROP_VALUE(TypeTag, OnlyGasPhaseCanDisappear);
+
+ if(onlyGasPhaseCanDisappear)
+ {
+ // check if a primary var switch is necessary
+ if (phasePresence == threePhases)
+ {
+ Scalar Smin = 0;
+ if (this->wasSwitched_[dofIdxGlobal])
+ Smin = -0.01;
+
+ if (volVars.saturation(gPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // gas phase disappears
+ std::cout << "Gas phase disappears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", sg: "
+ << volVars.saturation(gPhaseIdx) << std::endl;
+ newPhasePresence = wnPhaseOnly;
+
+ priVars[pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
+ priVars[switch1Idx] = volVars.fluidState().temperature();
+ }
+ }
+ else if (phasePresence == wnPhaseOnly)
+ {
+ bool gasFlag = 0;
+
+ // calculate fractions of the partial pressures in the
+ // hypothetical gas phase
+ Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+
+ Scalar xgMax = 1.0;
+ if (xwg + xng > xgMax)
+ wouldSwitch = true;
+ if (this->wasSwitched_[dofIdxGlobal])
+ xgMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, gas phase appears
+ if (xwg + xng > xgMax)
+ {
+ // gas phase appears
+ std::cout << "gas phase appears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", xwg + xng: "
+ << xwg + xng << std::endl;
+ gasFlag = 1;
+ }
+
+ if (gasFlag == 1)
+ {
+ newPhasePresence = threePhases;
+ priVars[pressureIdx] = volVars.pressure(gPhaseIdx);
+ priVars[switch1Idx] = volVars.saturation(wPhaseIdx);
+ }
+ }
+ }
+
+ else
+ {
+ // check if a primary var switch is necessary
+ if (phasePresence == threePhases)
+ {
+ Scalar Smin = 0;
+ if (this->wasSwitched_[dofIdxGlobal])
+ Smin = -0.01;
+
+ if (volVars.saturation(gPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // gas phase disappears
+ std::cout << "Gas phase disappears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", sg: "
+ << volVars.saturation(gPhaseIdx) << std::endl;
+ newPhasePresence = wnPhaseOnly;
+
+ priVars[pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
+ priVars[switch1Idx] = volVars.fluidState().temperature();
+ }
+ else if (volVars.saturation(wPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // water phase disappears
+ std::cout << "Water phase disappears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", sw: "
+ << volVars.saturation(wPhaseIdx) << std::endl;
+ newPhasePresence = gnPhaseOnly;
+
+ priVars[switch1Idx] = volVars.fluidState().saturation(nPhaseIdx);
+ priVars[switch2Idx] = volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
+ }
+ else if (volVars.saturation(nPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // NAPL phase disappears
+ std::cout << "NAPL phase disappears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", sn: "
+ << volVars.saturation(nPhaseIdx) << std::endl;
+ newPhasePresence = wgPhaseOnly;
+
+ priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ }
+ else if (phasePresence == wPhaseOnly)
+ {
+ bool gasFlag = 0;
+ bool nonwettingFlag = 0;
+ // calculate fractions of the partial pressures in the
+ // hypothetical gas phase
+ Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+
+ Scalar xgMax = 1.0;
+ if (xwg + xng > xgMax)
+ wouldSwitch = true;
+ if (this->wasSwitched_[dofIdxGlobal])
+ xgMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, gas phase appears
+ if (xwg + xng > xgMax)
+ {
+ // gas phase appears
+ std::cout << "gas phase appears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", xwg + xng: "
+ << xwg + xng << std::endl;
+ gasFlag = 1;
+ }
+
+ // calculate fractions in the hypothetical NAPL phase
+ Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
+ /* take care:
+ for xnn in case wPhaseOnly we compute xnn=henry_mesitylene*x1w,
+ where a hypothetical gas pressure is assumed for the Henry
+ x0n is set to NULL (all NAPL phase is dirty)
+ x2n is set to NULL (all NAPL phase is dirty)
+ */
+
+ Scalar xnMax = 1.0;
+ if (xnn > xnMax)
+ wouldSwitch = true;
+ if (this->wasSwitched_[dofIdxGlobal])
+ xnMax *= 1.02;
+
+ // if the sum of the hypothetical mole fractions would be larger than
+ // 100%, NAPL phase appears
+ if (xnn > xnMax)
+ {
+ // NAPL phase appears
+ std::cout << "NAPL phase appears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", xnn: "
+ << xnn << std::endl;
+ nonwettingFlag = 1;
+ }
+
+ if ((gasFlag == 1) && (nonwettingFlag == 0))
+ {
+ newPhasePresence = wgPhaseOnly;
+ priVars[switch1Idx] = 0.9999;
+ priVars[pressureIdx] = volVars.fluidState().pressure(gPhaseIdx);
+ }
+ else if ((gasFlag == 1) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = threePhases;
+ priVars[pressureIdx] = volVars.fluidState().pressure(gPhaseIdx);
+ priVars[switch1Idx] = 0.999;
+ }
+ else if ((gasFlag == 0) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = wnPhaseOnly;
+ priVars[switch2Idx] = 0.0001;
+ }
+ }
+ else if (phasePresence == gnPhaseOnly)
+ {
+ bool nonwettingFlag = 0;
+ bool wettingFlag = 0;
+
+ Scalar Smin = 0.0;
+ if (this->wasSwitched_[dofIdxGlobal])
+ Smin = -0.01;
+
+ if (volVars.saturation(nPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // NAPL phase disappears
+ std::cout << "NAPL phase disappears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", sn: "
+ << volVars.saturation(nPhaseIdx) << std::endl;
+ nonwettingFlag = 1;
+ }
+
+
+ // calculate fractions of the hypothetical water phase
+ Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
+ /*
+ take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
+ If this is larger than 1, then water appears
+ */
+ Scalar xwMax = 1.0;
+ if (xww > xwMax)
+ wouldSwitch = true;
+ if (this->wasSwitched_[dofIdxGlobal])
+ xwMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, water phase appears
+ if (xww > xwMax)
+ {
+ // water phase appears
+ std::cout << "water phase appears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
+ << xww << std::endl;
+ wettingFlag = 1;
+ }
+
+ if ((wettingFlag == 1) && (nonwettingFlag == 0))
+ {
+ newPhasePresence = threePhases;
+ priVars[switch1Idx] = 0.0001;
+ priVars[switch2Idx] = volVars.saturation(nPhaseIdx);
+ }
+ else if ((wettingFlag == 1) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = wgPhaseOnly;
+ priVars[switch1Idx] = 0.0001;
+ priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ else if ((wettingFlag == 0) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = gPhaseOnly;
+ priVars[switch1Idx] = volVars.fluidState().temperature();
+ priVars[switch2Idx] = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+ }
+ }
+ else if (phasePresence == wnPhaseOnly)
+ {
+ bool nonwettingFlag = 0;
+ bool gasFlag = 0;
+
+ Scalar Smin = 0.0;
+ if (this->wasSwitched_[dofIdxGlobal])
+ Smin = -0.01;
+
+ if (volVars.saturation(nPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // NAPL phase disappears
+ std::cout << "NAPL phase disappears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", sn: "
+ << volVars.saturation(nPhaseIdx) << std::endl;
+ nonwettingFlag = 1;
+ }
+
+ // calculate fractions of the partial pressures in the
+ // hypothetical gas phase
+ Scalar xwg = volVars.fluidState().moleFraction(gPhaseIdx, wCompIdx);
+ Scalar xng = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+
+ Scalar xgMax = 1.0;
+ if (xwg + xng > xgMax)
+ wouldSwitch = true;
+ if (this->wasSwitched_[dofIdxGlobal])
+ xgMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, gas phase appears
+ if (xwg + xng > xgMax)
+ {
+ // gas phase appears
+ std::cout << "gas phase appears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", xwg + xng: "
+ << xwg + xng << std::endl;
+ gasFlag = 1;
+ }
+
+ if ((gasFlag == 1) && (nonwettingFlag == 0))
+ {
+ newPhasePresence = threePhases;
+ priVars[pressureIdx] = volVars.pressure(gPhaseIdx);
+ priVars[switch1Idx] = volVars.saturation(wPhaseIdx);
+ }
+ else if ((gasFlag == 1) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = wgPhaseOnly;
+ priVars[pressureIdx] = volVars.pressure(gPhaseIdx);
+ priVars[switch1Idx] = volVars.temperature();
+ priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ else if ((gasFlag == 0) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = wPhaseOnly;
+ priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ }
+ else if (phasePresence == gPhaseOnly)
+ {
+ bool nonwettingFlag = 0;
+ bool wettingFlag = 0;
+
+ // calculate fractions in the hypothetical NAPL phase
+ Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
+ /*
+ take care:, xnn, if no NAPL phase is there, take xnn=xng*pg/pcsat
+ if this is larger than 1, then NAPL appears
+ */
+
+ Scalar xnMax = 1.0;
+ if (xnn > xnMax)
+ wouldSwitch = true;
+ if (this->wasSwitched_[dofIdxGlobal])
+ xnMax *= 1.02;
+
+ // if the sum of the hypothetical mole fraction would be larger than
+ // 100%, NAPL phase appears
+ if (xnn > xnMax)
+ {
+ // NAPL phase appears
+ std::cout << "NAPL phase appears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", xnn: "
+ << xnn << std::endl;
+ nonwettingFlag = 1;
+ }
+ // calculate fractions of the hypothetical water phase
+ Scalar xww = volVars.fluidState().moleFraction(wPhaseIdx, wCompIdx);
+ /*
+ take care:, xww, if no water is present, then take xww=xwg*pg/pwsat .
+ If this is larger than 1, then water appears
+ */
+ Scalar xwMax = 1.0;
+ if (xww > xwMax)
+ wouldSwitch = true;
+ if (this->wasSwitched_[dofIdxGlobal])
+ xwMax *= 1.02;
+
+ // if the sum of the mole fractions would be larger than
+ // 100%, gas phase appears
+ if (xww > xwMax)
+ {
+ // water phase appears
+ std::cout << "water phase appears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", xww=xwg*pg/pwsat : "
+ << xww << std::endl;
+ wettingFlag = 1;
+ }
+ if ((wettingFlag == 1) && (nonwettingFlag == 0))
+ {
+ newPhasePresence = wgPhaseOnly;
+ priVars[switch1Idx] = 0.0001;
+ priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ else if ((wettingFlag == 1) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = threePhases;
+ priVars[switch1Idx] = 0.0001;
+ priVars[switch2Idx] = 0.0001;
+ }
+ else if ((wettingFlag == 0) && (nonwettingFlag == 1))
+ {
+ newPhasePresence = gnPhaseOnly;
+ priVars[switch2Idx]
+ = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ }
+ else if (phasePresence == wgPhaseOnly)
+ {
+ bool nonwettingFlag = 0;
+ bool gasFlag = 0;
+ bool wettingFlag = 0;
+
+ // get the fractions in the hypothetical NAPL phase
+ Scalar xnn = volVars.fluidState().moleFraction(nPhaseIdx, nCompIdx);
+
+ // take care: if the NAPL phase is not present, take
+ // xnn=xng*pg/pcsat if this is larger than 1, then NAPL
+ // appears
+ Scalar xnMax = 1.0;
+ if (xnn > xnMax)
+ wouldSwitch = true;
+ if (this->wasSwitched_[dofIdxGlobal])
+ xnMax *= 1.02;
+
+ // if the sum of the hypothetical mole fraction would be larger than
+ // 100%, NAPL phase appears
+ if (xnn > xnMax)
+ {
+ // NAPL phase appears
+ std::cout << "NAPL phase appears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", xnn: "
+ << xnn << std::endl;
+ nonwettingFlag = 1;
+ }
+
+ Scalar Smin = -1.e-6;
+ if (this->wasSwitched_[dofIdxGlobal])
+ Smin = -0.01;
+
+ if (volVars.saturation(gPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // gas phase disappears
+ std::cout << "Gas phase disappears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", sg: "
+ << volVars.saturation(gPhaseIdx) << std::endl;
+ gasFlag = 1;
+ }
+
+ Smin = 0.0;
+ if (this->wasSwitched_[dofIdxGlobal])
+ Smin = -0.01;
+
+ if (volVars.saturation(wPhaseIdx) <= Smin)
+ {
+ wouldSwitch = true;
+ // gas phase disappears
+ std::cout << "Water phase disappears at vertex " << dofIdxGlobal
+ << ", coordinates: " << globalPos << ", sw: "
+ << volVars.saturation(wPhaseIdx) << std::endl;
+ wettingFlag = 1;
+ }
+
+ if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 1))
+ {
+ newPhasePresence = gnPhaseOnly;
+ priVars[switch1Idx] = 0.0001;
+ priVars[switch2Idx] = volVars.fluidState().moleFraction(nPhaseIdx, wCompIdx);
+ ;
+ }
+ else if ((gasFlag == 0) && (nonwettingFlag == 1) && (wettingFlag == 0))
+ {
+ newPhasePresence = threePhases;
+ priVars[switch2Idx] = 0.0001;
+ }
+ else if ((gasFlag == 1) && (nonwettingFlag == 0) && (wettingFlag == 0))
+ {
+ newPhasePresence = wPhaseOnly;
+ priVars[pressureIdx] = volVars.fluidState().pressure(wPhaseIdx);
+ priVars[switch1Idx] = volVars.fluidState().temperature();
+ priVars[switch2Idx] = volVars.fluidState().moleFraction(wPhaseIdx, nCompIdx);
+ }
+ else if ((gasFlag == 0) && (nonwettingFlag == 0) && (wettingFlag == 1))
+ {
+ newPhasePresence = gPhaseOnly;
+ priVars[switch1Idx]
+ = volVars.fluidState().temperature();
+ priVars[switch2Idx]
+ = volVars.fluidState().moleFraction(gPhaseIdx, nCompIdx);
+ }
+ }
+ }
+
+
+ priVars.setState(newPhasePresence);
+ this->wasSwitched_[dofIdxGlobal] = wouldSwitch;
+ return phasePresence != newPhasePresence;
+ }
+
+};
+
+} // end namespace dumux
+
+#endif
diff --git a/dumux/porousmediumflow/3pwateroil/properties.hh b/dumux/porousmediumflow/3pwateroil/properties.hh
deleted file mode 100644
index 2ba4bc7e387fe273caae02915fc2701e5396d8ed..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/3pwateroil/properties.hh
+++ /dev/null
@@ -1,75 +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 ThreePWaterOilModel
- */
-/*!
- * \file
- *
- * \brief Defines the properties required for the three-phase two-component
- * fully implicit model. This model requires the use of the non-isothermal
- * extension found in dumux/implicit/nonisothermal
- */
-#ifndef DUMUX_3P2CNI_PROPERTIES_HH
-#define DUMUX_3P2CNI_PROPERTIES_HH
-
-#include <dumux/implicit/box/properties.hh>
-#include <dumux/implicit/cellcentered/properties.hh>
-#include <dumux/porousmediumflow/nonisothermal/model.hh>
-
-namespace Dumux
-{
-
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
-NEW_TYPE_TAG(ThreePWaterOilNI, INHERITS_FROM(NonIsothermal));
-NEW_TYPE_TAG(BoxThreePWaterOilNI, INHERITS_FROM(BoxModel, ThreePWaterOilNI));
-NEW_TYPE_TAG(CCThreePWaterOilNI, INHERITS_FROM(CCModel, ThreePWaterOilNI));
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-
-NEW_PROP_TAG(NumPhases); //!< Number of fluid phases in the system
-NEW_PROP_TAG(NumComponents); //!< Number of fluid components in the system
-NEW_PROP_TAG(Indices); //!< Enumerations for the model
-NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters
-NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
-
-NEW_PROP_TAG(MaterialLaw); //!< The material law which ought to be used (extracted from the spatial parameters)
-
-NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
-
-NEW_PROP_TAG(UseMoles); //!< Defines whether mole (true) or mass (false) fractions are used
-
-NEW_PROP_TAG(UseMassOutput); //!< Defines whether mole or mass are used for phaseStorage output
-NEW_PROP_TAG(EffectiveDiffusivityModel); //!< The employed model for the computation of the effective diffusivity
-
-NEW_PROP_TAG(ImplicitMassUpwindWeight); //!< The value of the upwind parameter for the mobility
-NEW_PROP_TAG(ImplicitMobilityUpwindWeight); //!< Weight for the upwind mobility in the velocity calculation
-NEW_PROP_TAG(UseSimpleModel); //!< Determines whether a simple model with only two phase states (wng) and (wn) should be used
-NEW_PROP_TAG(BaseFluxVariables); //!< The base flux variables
-NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
-}
-}
-
-#endif
diff --git a/dumux/porousmediumflow/3pwateroil/propertydefaults.hh b/dumux/porousmediumflow/3pwateroil/propertydefaults.hh
deleted file mode 100644
index 3de1dacb7653b8b5e4d6e584050f5d7525acac6a..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/3pwateroil/propertydefaults.hh
+++ /dev/null
@@ -1,172 +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 ThreePWaterOilModel
- */
-/*!
- * \file
- *
- * \brief Defines default values for most properties required by the
- * 3p2cni fully implicit model.
- */
-#ifndef DUMUX_3P2CNI_PROPERTY_DEFAULTS_HH
-#define DUMUX_3P2CNI_PROPERTY_DEFAULTS_HH
-
-#include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
-#include <dumux/material/spatialparams/fv.hh>
-#include <dumux/material/fluidmatrixinteractions/3p/thermalconductivitysomerton3p.hh>
-#include <dumux/porousmediumflow/nonisothermal/model.hh>
-#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
-
-#include "indices.hh"
-#include "model.hh"
-#include "fluxvariables.hh"
-#include "volumevariables.hh"
-#include "properties.hh"
-#include "newtoncontroller.hh"
-#include "localresidual.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 2.
- */
-SET_PROP(ThreePWaterOilNI, NumComponents)
-{
- private:
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-
- public:
- static const int value = FluidSystem::numComponents;
-
- static_assert(value == 2,
- "Only fluid systems with 2 components are supported by the 3p2cni 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(ThreePWaterOilNI, NumPhases)
-{
- private:
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-
- public:
- static const int value = FluidSystem::numPhases;
- static_assert(value == 3,
- "Only fluid systems with 3 phases are supported by the 3p2cni model!");
-};
-
-//! Use the 3p2cni specific newton controller for the 3p2cni model
-SET_TYPE_PROP(ThreePWaterOilNI, NewtonController, ThreePWaterOilNewtonController<TypeTag>);
-
-//! define the base flux variables to realize Darcy flow
-SET_TYPE_PROP(ThreePWaterOilNI, BaseFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
-
-//! the upwind factor for the mobility.
-SET_SCALAR_PROP(ThreePWaterOilNI, ImplicitMassUpwindWeight, 1.0);
-
-//! set default mobility upwind weight to 1.0, i.e. fully upwind
-SET_SCALAR_PROP(ThreePWaterOilNI, ImplicitMobilityUpwindWeight, 1.0);
-
-//! Determines whether a constraint solver should be used explicitly
-SET_BOOL_PROP(ThreePWaterOilNI, UseSimpleModel, true);
-
-//! The spatial parameters to be employed.
-//! Use FVSpatialParams by default.
-SET_TYPE_PROP(ThreePWaterOilNI, SpatialParams, FVSpatialParams<TypeTag>);
-
-//! Use the model after Millington (1961) for the effective diffusivity
-SET_TYPE_PROP(ThreePWaterOilNI, EffectiveDiffusivityModel,
- DiffusivityMillingtonQuirk<typename GET_PROP_TYPE(TypeTag, Scalar)>);
-
-// disable velocity output by default
-
-// enable gravity by default
-SET_BOOL_PROP(ThreePWaterOilNI, ProblemEnableGravity, true);
-
-SET_BOOL_PROP(ThreePWaterOilNI, UseMoles, true); //!< Define that mole fractions are used in the balance equations per default
-
-
-
-//! 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(ThreePWaterOilNI, SpatialParamsForchCoeff, 0.55);
-
-
-//! Somerton is used as default model to compute the effective thermal heat conductivity
-SET_PROP(ThreePWaterOilNI, ThermalConductivityModel)
-{
-private:
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-public:
- using type = ThermalConductivitySomerton<Scalar, Indices>;
-};
-
-//! temperature is already written by the isothermal model
-SET_BOOL_PROP(ThreePWaterOilNI, NiOutputLevel, 0);
-
-//////////////////////////////////////////////////////////////////
-// Property values for isothermal model required for the general non-isothermal model
-//////////////////////////////////////////////////////////////////
-
-// set isothermal Model
-SET_TYPE_PROP(ThreePWaterOilNI, IsothermalModel, ThreePWaterOilModel<TypeTag>);
-
-// set isothermal FluxVariables
-SET_TYPE_PROP(ThreePWaterOilNI, IsothermalFluxVariables, ThreePWaterOilFluxVariables<TypeTag>);
-
-//set isothermal VolumeVariables
-SET_TYPE_PROP(ThreePWaterOilNI, IsothermalVolumeVariables, ThreePWaterOilVolumeVariables<TypeTag>);
-
-//set isothermal LocalResidual
-SET_TYPE_PROP(ThreePWaterOilNI, IsothermalLocalResidual, ThreePWaterOilLocalResidual<TypeTag>);
-
-//set isothermal Indices
-SET_PROP(ThreePWaterOilNI, IsothermalIndices)
-{
-
-public:
- using type = ThreePWaterOilIndices<TypeTag, /*PVOffset=*/0>;
-};
-
-//set isothermal NumEq
-SET_INT_PROP(ThreePWaterOilNI, IsothermalNumEq, 2);
-
-}
-
-}
-
-#endif
diff --git a/dumux/porousmediumflow/3pwateroil/volumevariables.hh b/dumux/porousmediumflow/3pwateroil/volumevariables.hh
index d87d3eb8a95b1761c6b7cd6f739f05305099c729..b48f932adac6f620bc183cd8e034424e689e4c9d 100644
--- a/dumux/porousmediumflow/3pwateroil/volumevariables.hh
+++ b/dumux/porousmediumflow/3pwateroil/volumevariables.hh
@@ -18,57 +18,54 @@
*****************************************************************************/
/*!
* \file
- *
+ * \ingroup ThreePWaterOilModel
* \brief Contains the quantities which are constant within a
* finite volume in the three-phase, two-component model.
*/
#ifndef DUMUX_3P2CNI_VOLUME_VARIABLES_HH
#define DUMUX_3P2CNI_VOLUME_VARIABLES_HH
-#include <dumux/implicit/model.hh>
#include <dumux/common/math.hh>
#include <dune/common/parallel/collectivecommunication.hh>
#include <vector>
#include <iostream>
+#include <dumux/porousmediumflow/volumevariables.hh>
+
#include <dumux/material/constants.hh>
#include <dumux/material/fluidstates/compositional.hh>
#include <dumux/material/constraintsolvers/computefromreferencephase.hh>
#include <dumux/material/constraintsolvers/misciblemultiphasecomposition.hh>
-
#include <dune/common/deprecated.hh>
-#include "properties.hh"
-
namespace Dumux
{
/*!
* \ingroup ThreePWaterOilModel
* \brief Contains the quantities which are are constant within a
- * finite volume in the two-phase, two-component model.
+ * finite volume in the three-phase, two-component model.
*/
template <class TypeTag>
-class ThreePWaterOilVolumeVariables : public ImplicitVolumeVariables<TypeTag>
+class ThreePWaterOilVolumeVariables : public PorousMediumFlowVolumeVariables<TypeTag>
{
- using ParentType = ImplicitVolumeVariables<TypeTag>;
+ using ParentType = PorousMediumFlowVolumeVariables<TypeTag>;
using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
- using MaterialLawParams = typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
// constraint solvers
using MiscibleMultiPhaseComposition = Dumux::MiscibleMultiPhaseComposition<Scalar, FluidSystem>;
using ComputeFromReferencePhase = Dumux::ComputeFromReferencePhase<Scalar, FluidSystem>;
-
- using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum {
dim = GridView::dimension,
@@ -77,6 +74,7 @@ class ThreePWaterOilVolumeVariables : public ImplicitVolumeVariables<TypeTag>
wCompIdx = Indices::wCompIdx,
nCompIdx = Indices::nCompIdx,
+ gCompIdx = Indices::gCompIdx,
wPhaseIdx = Indices::wPhaseIdx,
gPhaseIdx = Indices::gPhaseIdx,
@@ -101,42 +99,33 @@ class ThreePWaterOilVolumeVariables : public ImplicitVolumeVariables<TypeTag>
static const Scalar R; // universial gas constant
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+ enum { isBox = GET_PROP_VALUE(TypeTag, DiscretizationMethod) == DiscretizationMethods::Box };
enum { dofCodim = isBox ? dim : 0 };
public:
//! The type of the object returned by the fluidState() method
- using FluidState = Dumux::CompositionalFluidState<Scalar, FluidSystem>;
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+ using ParentType::enthalpy;
/*!
* \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(priVars,
- problem,
- element,
- fvGeometry,
- scvIdx,
- isOldSol);
+ ParentType::update(elemSol, problem, element, scv);
+ const auto& priVars = ParentType::extractDofPriVars(elemSol, scv);
+ const auto phasePresence = priVars.state();
- bool useSimpleModel = GET_PROP_VALUE(TypeTag, UseSimpleModel);
+ bool onlyGasPhaseCanDisappear = GET_PROP_VALUE(TypeTag, OnlyGasPhaseCanDisappear);
// capillary pressure parameters
- const MaterialLawParams &materialParams =
- problem.spatialParams().materialLawParams(element, fvGeometry, scvIdx);
+ 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);
-
- if(!useSimpleModel)
+ if(!onlyGasPhaseCanDisappear)
{
/* first the saturations */
if (phasePresence == threePhases)
@@ -650,8 +639,7 @@ public:
fluidState_.setDensity(gPhaseIdx, rhoG);
fluidState_.setDensity(nPhaseIdx, rhoN);
}
- else
- assert(false); // unhandled phase state
+ else DUNE_THROW(Dune::InvalidStateException, "phasePresence: " << phasePresence << " is invalid.");
}
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
@@ -678,37 +666,27 @@ public:
* for the porous media happens at another place!
*/
- // diffusivity coefficients
diffusionCoefficient_[gPhaseIdx] = FluidSystem::diffusionCoefficient(fluidState_, gPhaseIdx);
-
diffusionCoefficient_[wPhaseIdx] = FluidSystem::diffusionCoefficient(fluidState_, wPhaseIdx);
-
/* no diffusion in NAPL phase considered at the moment, dummy values */
diffusionCoefficient_[nPhaseIdx] = 1.e-10;
- Valgrind::CheckDefined(diffusionCoefficient_);
// porosity
- porosity_ = problem.spatialParams().porosity(element,
- fvGeometry,
- scvIdx);
+ porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
Valgrind::CheckDefined(porosity_);
// permeability
- permeability_ = problem.spatialParams().intrinsicPermeability(element,
- fvGeometry,
- scvIdx);
+ permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
Valgrind::CheckDefined(permeability_);
+// fluidState_.setTemperature(temp_);
// 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);
- }
-
- // energy related quantities not contained in the fluid state
- asImp_().updateEnergy_(priVars, problem, element, fvGeometry, scvIdx, isOldSol);
+ }
}
/*!
@@ -790,9 +768,10 @@ public:
* \param phaseIdx The phase index
*/
Scalar mobility(const int phaseIdx) const
- {
- return mobility_[phaseIdx];
- }
+ { return mobility_[phaseIdx]; }
+
+ Scalar viscosity(const int phaseIdx) const
+ { return fluidState_.viscosity(phaseIdx); }
/*!
* \brief Returns the effective capillary pressure within the control volume.
@@ -815,17 +794,19 @@ public:
/*!
* \brief Returns the diffusivity coefficient matrix
*/
- DUNE_DEPRECATED_MSG("diffusionCoefficient() is deprecated. Use diffCoeff(int phaseIdx) instead.")
+ DUNE_DEPRECATED_MSG("diffusionCoefficient() is deprecated. Use diffusionCoefficient(int phaseIdx) instead.")
Dune::FieldVector<Scalar, numPhases> diffusionCoefficient() const
{
return diffusionCoefficient_;
}
/*!
- * \brief Returns the diffusivity coefficient
+ * \brief Returns the diffusion coeffiecient
*/
- Scalar diffCoeff(int phaseIdx) const
- { return diffusionCoefficient_[phaseIdx]; }
+ Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
+ {
+ return diffusionCoefficient_[phaseIdx];
+ }
/*!
* \brief Returns the adsorption information
@@ -853,17 +834,6 @@ public:
protected:
- /*!
- * \brief Called by update() to compute the energy related quantities
- */
- void updateEnergy_(const PrimaryVariables &priVars,
- const Problem &problem,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx,
- bool isOldSol)
- { }
-
Scalar sw_, sg_, sn_, pg_, pw_, pn_, temp_;
Scalar moleFrac_[numPhases][numComponents];
@@ -879,6 +849,8 @@ protected:
FluidState fluidState_;
private:
+ std::array<std::array<Scalar, numComponents-1>, numPhases> diffCoefficient_;
+
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
diff --git a/dumux/porousmediumflow/3pwateroil/vtkoutputfields.hh b/dumux/porousmediumflow/3pwateroil/vtkoutputfields.hh
new file mode 100644
index 0000000000000000000000000000000000000000..079e1dce9292a8ae047dae1f3c6540a8f63304fd
--- /dev/null
+++ b/dumux/porousmediumflow/3pwateroil/vtkoutputfields.hh
@@ -0,0 +1,82 @@
+// -*- 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
+ * \ingroup ThreePWaterOilModel
+ * \brief Adds vtk output fields specific to the twop model
+ */
+#ifndef DUMUX_THREEPWATEROIL_VTK_OUTPUT_FIELDS_HH
+#define DUMUX_THREEPWATEROIL_VTK_OUTPUT_FIELDS_HH
+
+#include <dumux/common/properties.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup ThreePWaterOilModel
+ * \brief Adds vtk output fields specific to the three-phase three-component model
+ */
+template<class TypeTag>
+class ThreePWaterOilVtkOutputFields
+{
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+ using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+
+ enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
+ numComponents = GET_PROP_VALUE(TypeTag, NumComponents)
+ };
+public:
+ template <class VtkOutputModule>
+ static void init(VtkOutputModule& vtk)
+ {
+ // register standardized vtk output fields
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.saturation(Indices::wPhaseIdx); }, "sw");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.saturation(Indices::nPhaseIdx); },"sn");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.saturation(Indices::gPhaseIdx); },"sg");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.pressure(Indices::wPhaseIdx); },"pw");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.pressure(Indices::nPhaseIdx); },"pn");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.pressure(Indices::gPhaseIdx); },"pg");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.density(Indices::wPhaseIdx); },"rhow");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.density(Indices::nPhaseIdx); },"rhon");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.density(Indices::gPhaseIdx); },"rhog");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return
+ v.mobility(Indices::wPhaseIdx); },"MobW");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.mobility(Indices::nPhaseIdx); },"MobN");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.mobility(Indices::gPhaseIdx); },"MobG");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return
+ v.viscosity(Indices::wPhaseIdx); },"ViscosW");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.viscosity(Indices::nPhaseIdx); },"ViscosN");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.viscosity(Indices::gPhaseIdx); },"ViscosG");
+
+ for (int i = 0; i < numPhases; ++i)
+ for (int j = 0; j < numComponents; ++j)
+ vtk.addVolumeVariable([i,j](const VolumeVariables& v){ return v.moleFraction(i,j); },
+ "x^" + FluidSystem::phaseName(i) + "_" + FluidSystem::componentName(j));
+
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.porosity(); },"porosity");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.priVars().state(); },"phase presence");
+ vtk.addVolumeVariable( [](const VolumeVariables& v){ return v.permeability(); },"permeability");
+ }
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh b/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh
index 7ff97306dbb0acf66da383c7b1833da95256b589..ec746506441ef4beec445b9af567f51bf78d64db 100644
--- a/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh
+++ b/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdproblem.hh
@@ -24,8 +24,9 @@
#ifndef DUMUX_SAGDPROBLEM_HH
#define DUMUX_SAGDPROBLEM_HH
-#include <dumux/porousmediumflow/implicit/problem.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/discretization/box/properties.hh>
#include <dumux/porousmediumflow/3pwateroil/model.hh>
#include <dumux/material/fluidsystems/h2oheavyoilfluidsystem.hh>
#include "3pwateroilsagdspatialparams.hh"
@@ -56,19 +57,9 @@ SET_TYPE_PROP(SagdProblem,
FluidSystem,
Dumux::FluidSystems::H2OHeavyOil<typename GET_PROP_TYPE(TypeTag, Scalar)>);
+SET_BOOL_PROP(SagdProblem, OnlyGasPhaseCanDisappear, true);
-// Enable gravity
-SET_BOOL_PROP(SagdProblem, ProblemEnableGravity, true);
-
-// Use forward differences instead of central differences
-SET_INT_PROP(SagdProblem, ImplicitNumericDifferenceMethod, +1);
-
-// Write newton convergence
-SET_BOOL_PROP(SagdProblem, NewtonWriteConvergence, false);
-
-SET_BOOL_PROP(SagdProblem, UseSimpleModel, true);
-
-SET_BOOL_PROP(SagdProblem, UseMoles, false);
+SET_BOOL_PROP(SagdProblem, UseMoles, true);
}
@@ -81,15 +72,13 @@ SET_BOOL_PROP(SagdProblem, UseMoles, false);
*
* */
template <class TypeTag >
-class SagdProblem : public ImplicitPorousMediaProblem<TypeTag>
+class SagdProblem : public PorousMediumFlowProblem<TypeTag>
{
- using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
- using Grid = typename GridView::Grid;
-
- using ParentType = ImplicitPorousMediaProblem<TypeTag>;
-
- // copy some indices for convenience
+ using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
enum {
pressureIdx = Indices::pressureIdx,
@@ -97,6 +86,7 @@ class SagdProblem : public ImplicitPorousMediaProblem<TypeTag>
switch2Idx = Indices::switch2Idx,
energyEqIdx = Indices::energyEqIdx,
+ temperatureIdx = Indices::temperatureIdx,
// phase and component indices
wPhaseIdx = Indices::wPhaseIdx,
@@ -111,30 +101,21 @@ class SagdProblem : public ImplicitPorousMediaProblem<TypeTag>
wgPhaseOnly = Indices::wgPhaseOnly,
threePhases = Indices::threePhases,
- //contiWEqIdx = Indices::contiWEqIdx,
- //contiNEqIdx = Indices::contiNEqIdx,
// Grid and world dimension
dim = GridView::dimension,
dimWorld = GridView::dimensionworld
};
-
using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using NeumannFluxes = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using Sources = typename GET_PROP_TYPE(TypeTag, NumEqVector);
+ using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
using BoundaryTypes = typename GET_PROP_TYPE(TypeTag, BoundaryTypes);
- using TimeManager = typename GET_PROP_TYPE(TypeTag, TimeManager);
-
using Element = typename GridView::template Codim<0>::Entity;
- using Vertex = typename GridView::template Codim<dim>::Entity;
- using Intersection = typename GridView::Intersection;
-
using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
- using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-
+ using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
using GlobalPosition = Dune::FieldVector<Scalar, dimWorld>;
- enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
-
public:
/*!
@@ -144,8 +125,8 @@ public:
* \param gridView The grid view
*/
- SagdProblem(TimeManager &timeManager, const GridView &gridView)
- : ParentType(timeManager, gridView), pOut_(4e6)
+ SagdProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry), pOut_(4e6)
{
maxDepth_ = 400.0; // [m]
@@ -153,72 +134,9 @@ public:
totalMassProducedOil_ =0;
totalMassProducedWater_ =0;
- this->timeManager().startNextEpisode(86400);
-
- name_ = GET_RUNTIME_PARAM(TypeTag, std::string, Problem.Name);
+ name_ = getParam<std::string>("Problem.Name");
}
- bool shouldWriteRestartFile() const
- {
- return 0;
- }
-
-
- /*!
- * \brief Called directly after the time integration.
- */
- void postTimeStep()
- {
- double time = this->timeManager().time();
- double dt = this->timeManager().timeStepSize();
-
- // Calculate storage terms
- PrimaryVariables storage;
- this->model().globalStorage(storage);
-
- // Write mass balance information for rank 0
- if (this->gridView().comm().rank() == 0)
- {
- std::cout<<"Storage: " << storage << "Time: " << time+dt << std::endl;
- massBalance.open ("massBalance.txt", std::ios::out | std::ios::app );
- massBalance << " Storage " << storage
- << " Time " << time+dt
- << std::endl;
- massBalance.close();
-
- }
-
- // Calculate storage terms
- PrimaryVariables storageW, storageN;
- //Dune::FieldVector<Scalar, 2> flux(0.0);
- this->model().globalPhaseStorage(storageW, wPhaseIdx);
- this->model().globalPhaseStorage(storageN, nPhaseIdx);
-
- //mass of Oil
- const Scalar newMassProducedOil_ = massProducedOil_;
- std::cout<<" newMassProducedOil_ : "<< newMassProducedOil_ << " Time: " << time+dt << std::endl;
-
- totalMassProducedOil_ += newMassProducedOil_;
- std::cout<<" totalMassProducedOil_ : "<< totalMassProducedOil_ << " Time: " << time+dt << std::endl;
- //mass of Water
- const Scalar newMassProducedWater_ = massProducedWater_;
- //std::cout<<" newMassProducedWater_ : "<< newMassProducedWater_ << " Time: " << time+dt << std::endl;
-
- totalMassProducedWater_ += newMassProducedWater_;
- //std::cout<<" totalMassProducedWater_ : "<< totalMassProducedWater_ << " Time: " << time+dt << std::endl;
-
-
- const int timeStepIndex = this->timeManager().timeStepIndex();
-
- if (timeStepIndex == 0 ||
- timeStepIndex % 100 == 0 || //after every 1000000 secs
- this->timeManager().episodeWillBeFinished() ||
- this->timeManager().willBeFinished())
- {
- std::cout<<" totalMassProducedOil_ : "<< totalMassProducedOil_ << " Time: " << time+dt << std::endl;
- std::cout<<" totalMassProducedWater_ : "<< totalMassProducedWater_ << " Time: " << time+dt << std::endl;
- }
- }
void episodeEnd()
{
@@ -241,13 +159,6 @@ public:
const std::string name() const
{ return name_; }
-
- void sourceAtPos(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
- {
- values = 0.0;
- }
-
// \}
/*!
@@ -262,9 +173,9 @@ public:
* \param bcTypes The boundary types for the conservation equations
* \param globalPos The position for which the bc type should be evaluated
*/
- void boundaryTypesAtPos(BoundaryTypes &bcTypes,
- const GlobalPosition &globalPos) const
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
+ BoundaryTypes bcTypes;
// on bottom
if (globalPos[1] < eps_)
bcTypes.setAllNeumann();
@@ -280,6 +191,7 @@ public:
// on Left
else
bcTypes.setAllNeumann();
+ return bcTypes;
}
/*!
@@ -291,9 +203,9 @@ public:
*
* For this method, the \a values parameter stores primary variables.
*/
- void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+ PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
- initial_(values, globalPos);
+ return initial_(globalPos);
}
/*!
@@ -311,22 +223,14 @@ public:
* For this method, the \a values parameter stores the mass flux
* in normal direction of each phase. Negative values mean influx.
*/
- void solDependentNeumann(PrimaryVariables &values,
- const Element &element,
- const FVElementGeometry &fvGeometry,
- const Intersection &is,
- const int scvIdx,
- const int boundaryFaceIdx,
- const ElementVolumeVariables &elemVolVars) const
+ NeumannFluxes neumann(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolumeFace& scvf) const
{
- values = 0;
-
- GlobalPosition globalPos;
- if (isBox)
- globalPos = element.geometry().corner(scvIdx);
- else
- globalPos = is.geometry().center();
+ NeumannFluxes values(0.0);
+ const auto& globalPos = scvf.ipGlobal();
// negative values for injection at injection well
if (globalPos[1] > 8.5 - eps_ && globalPos[1] < 9.5 + eps_)
{
@@ -337,17 +241,17 @@ public:
else if (globalPos[1] > 2.5 - eps_ && globalPos[1] < 3.5 + eps_) // production well
{
- const Scalar elemPressW = elemVolVars[scvIdx].pressure(wPhaseIdx); //Pressures
- const Scalar elemPressN = elemVolVars[scvIdx].pressure(nPhaseIdx);
+ const Scalar elemPressW = elemVolVars[scvf.insideScvIdx()].pressure(wPhaseIdx); //Pressures
+ const Scalar elemPressN = elemVolVars[scvf.insideScvIdx()].pressure(nPhaseIdx);
- const Scalar densityW = elemVolVars[scvIdx].fluidState().density(wPhaseIdx); //Densities
- const Scalar densityN = elemVolVars[scvIdx].fluidState().density(nPhaseIdx);
+ const Scalar densityW = elemVolVars[scvf.insideScvIdx()].fluidState().density(wPhaseIdx); //Densities
+ const Scalar densityN = elemVolVars[scvf.insideScvIdx()].fluidState().density(nPhaseIdx);
- const Scalar elemMobW = elemVolVars[scvIdx].mobility(wPhaseIdx); //Mobilities
- const Scalar elemMobN = elemVolVars[scvIdx].mobility(nPhaseIdx);
+ const Scalar elemMobW = elemVolVars[scvf.insideScvIdx()].mobility(wPhaseIdx); //Mobilities
+ const Scalar elemMobN = elemVolVars[scvf.insideScvIdx()].mobility(nPhaseIdx);
- const Scalar enthW = elemVolVars[scvIdx].enthalpy(wPhaseIdx); //Enthalpies
- const Scalar enthN = elemVolVars[scvIdx].enthalpy(nPhaseIdx);
+ const Scalar enthW = elemVolVars[scvf.insideScvIdx()].enthalpy(wPhaseIdx); //Enthalpies
+ const Scalar enthN = elemVolVars[scvf.insideScvIdx()].enthalpy(nPhaseIdx);
const Scalar wellRadius = 0.50 * 0.3048; // 0.50 ft as specified by SPE9
@@ -368,7 +272,7 @@ public:
// qE = qW*0.018*enthW + qN*enthN*0.350;
//with cooling: see Diplomarbeit Stefan Roll, Sept. 2015
- Scalar wT = elemVolVars[scvIdx].temperature(); // well temperature
+ Scalar wT = elemVolVars[scvf.insideScvIdx()].temperature(); // well temperature
if ( wT > 495. )
{
qE = qW*0.018*enthW + qN*enthN*0.350 + (wT-495.)*5000.; // ~3x injected enthalpy
@@ -384,6 +288,7 @@ public:
massProducedOil_ = qN;
massProducedWater_ = qW;
}
+ return values;
}
// \}
@@ -402,36 +307,24 @@ public:
* 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 wnPhaseOnly;
+ return initial_(globalPos);
}
private:
// internal method for the initial condition (reused for the
// dirichlet conditions!)
- void initial_(PrimaryVariables &values,
- const GlobalPosition &globalPos) const
+ PrimaryVariables initial_(const GlobalPosition &globalPos) const
{
+ PrimaryVariables values(0.0);
+ values.setState(Indices::wnPhaseOnly);
Scalar densityW = 1000.0;
values[pressureIdx] = 101300.0 + (maxDepth_ - globalPos[1])*densityW*9.81;
values[switch1Idx] = 295.13; // temperature
values[switch2Idx] = 0.3; //NAPL saturation
+ return values;
}
Scalar maxDepth_;
diff --git a/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdspatialparams.hh b/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdspatialparams.hh
index d253b4c9229e2030d2c92b14e9a4d014c985e974..1a217f300d5895d59952dc5d553796602bee8740 100644
--- a/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdspatialparams.hh
+++ b/test/porousmediumflow/3pwateroil/implicit/3pwateroilsagdspatialparams.hh
@@ -74,11 +74,15 @@ template<class TypeTag>
class SagdSpatialParams : public FVSpatialParams<TypeTag>
{
using ParentType = FVSpatialParams<TypeTag>;
-
using Grid = typename GET_PROP_TYPE(TypeTag, Grid);
- using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
- using CoordScalar = typename Grid::ctype;
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using MaterialLawParams = typename MaterialLaw::Params;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
enum {
dim=GridView::dimension,
dimWorld=GridView::dimensionworld
@@ -90,30 +94,21 @@ class SagdSpatialParams : public FVSpatialParams<TypeTag>
nPhaseIdx = Indices::nPhaseIdx
};
- using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
- using DimVector = Dune::FieldVector<CoordScalar, dimWorld>;
-
-
- using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
-
- using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
- using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-
- using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+ using CoordScalar = typename Grid::ctype;
+ using GlobalPosition = Dune::FieldVector<CoordScalar,dimWorld>;
using Element = typename GridView::template Codim<0>::Entity;
-
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
public:
- using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
- using MaterialLawParams = typename MaterialLaw::Params;
+ using PermeabilityType = Scalar;
/*!
* \brief The constructor
*
* \param gridView The grid view
*/
- SagdSpatialParams(const GridView &gridView)
- : ParentType(gridView), eps_(1e-6)
+ SagdSpatialParams(const Problem& problem)
+ : ParentType(problem), eps_(1e-6)
{
layerBottom_ = 35.0;
@@ -152,60 +147,69 @@ public:
fineMaterialParams_.setKrRegardsSnr(false);
}
- ~SagdSpatialParams()
- {}
+ /*!
+ * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$
+ * \note It is possibly solution dependent.
+ *
+ * \param element The current element
+ * \param scv The sub-control volume inside the element.
+ * \param elemSol The solution at the dofs connected to the element.
+ * \return permeability
+ */
+ PermeabilityType permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ { return permeabilityAtPos(scv.dofPosition());}
/*!
- * \brief Apply the intrinsic permeability tensor to a pressure
- * potential gradient.
+ * \brief Returns the intrinsic permeability tensor \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 Scalar intrinsicPermeability(const Element &element,
- const FVElementGeometry &fvElemGeom,
- int scvIdx) const
+ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
{
- const GlobalPosition &pos = fvElemGeom.subContVol[scvIdx].global;
- if (isFineMaterial_(pos))
+ if (isFineMaterial_(globalPos))
return fineK_;
return coarseK_;
}
/*!
- * \brief Define the porosity \f$[-]\f$ of the spatial parameters
+ * \brief Returns the porosity \f$[-]\f$
*
- * \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
*/
- Scalar porosity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
{
- const GlobalPosition &pos = fvGeometry.subContVol[scvIdx].global;
- if (isFineMaterial_(pos))
+ if (isFineMaterial_(globalPos))
return finePorosity_;
else
return coarsePorosity_;
}
+ /*!
+ * \brief Function for defining the parameters needed by constitutive relationships (kr-sw, pc-sw, etc.).
+ *
+ * \param element The current element
+ * \param scv The sub-control volume inside the element.
+ * \param elemSol The solution at the dofs connected to the element.
+ * \return the material parameters object
+ */
+ const MaterialLawParams& materialLawParams(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolutionVector& elemSol) const
+ {
+ return materialLawParamsAtPos(scv.dofPosition());
+ }
/*!
- * \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
{
- const GlobalPosition &pos = fvGeometry.subContVol[scvIdx].global;
- if (isFineMaterial_(pos))
+ if (isFineMaterial_(globalPos))
return fineMaterialParams_;
else
return coarseMaterialParams_;
@@ -220,12 +224,9 @@ public:
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume
*/
- Scalar solidHeatCapacity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Scalar solidHeatCapacityAtPos(const GlobalPosition& globalPos) const
{
- const GlobalPosition &pos = fvGeometry.subContVol[scvIdx].global;
- if (isFineMaterial_(pos))
+ if (isFineMaterial_(globalPos))
return fineHeatCap_ ;
else
return coarseHeatCap_;
@@ -240,9 +241,7 @@ public:
* \param fvGeometry The finite volume geometry
* \param scvIdx The local index of the sub-control volume
*/
- Scalar solidDensity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Scalar solidDensityAtPos(const GlobalPosition& globalPos) const
{
return 2650; // density of sand [kg/m^3]
}
@@ -255,9 +254,7 @@ public:
* \param scvIdx The local index of the sub-control volume where
* the heat capacity needs to be defined
*/
- Scalar solidThermalConductivity(const Element &element,
- const FVElementGeometry &fvGeometry,
- const int scvIdx) const
+ Scalar solidThermalConductivityAtPos(const GlobalPosition& globalPos) const
{
return lambdaSolid_;
}
diff --git a/test/porousmediumflow/3pwateroil/implicit/CMakeLists.txt b/test/porousmediumflow/3pwateroil/implicit/CMakeLists.txt
index b0428e98ea7e02a51fb186faaf2262a2e8ceaf2e..f902e222a3cbd0a145eab5d0d81da4add741f612 100644
--- a/test/porousmediumflow/3pwateroil/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/3pwateroil/implicit/CMakeLists.txt
@@ -1,13 +1,13 @@
add_input_file_links()
-dune_add_test(
- NAME test_box3pwateroil
- SOURCES test_box3pwateroil.cc
- COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
- CMD_ARGS --script fuzzy
- --files ${CMAKE_SOURCE_DIR}/test/references/3pwateroilbox-reference.vtu
+dune_add_test(NAME test_box3pwateroil
+ SOURCES test_box3pwateroil.cc
+ COMPILE_DEFINITIONS TYPETAG=ThreePWaterOilSagdBoxProblem
+ COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ CMD_ARGS --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/3pwateroilbox-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/sagd-00011.vtu
- --command "${CMAKE_CURRENT_BINARY_DIR}/test_box3pwateroil"
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_box3pwateroil"
)
#install sources
diff --git a/test/porousmediumflow/3pwateroil/implicit/test_box3pwateroil.cc b/test/porousmediumflow/3pwateroil/implicit/test_box3pwateroil.cc
index 9d42852d6a2d3d20d6e3291231d7212f3879f435..ef46ef7d475312fb4aa98379c44a422ea3e5b7fd 100644
--- a/test/porousmediumflow/3pwateroil/implicit/test_box3pwateroil.cc
+++ b/test/porousmediumflow/3pwateroil/implicit/test_box3pwateroil.cc
@@ -18,19 +18,41 @@
*****************************************************************************/
/*!
* \file
- *
- * \brief test for the 3p3cni box model
+ * \ingroup OnePTests
+ * \brief Test for the three-phase three-component box model
*/
-#include "config.h"
-#include "3pwateroilsagdproblem.hh"
-#include <dumux/common/start.hh>
+#include <config.h>
+#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 <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+#include <dumux/common/defaultusagemessage.hh>
+#include <dumux/common/timeloop.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/nonlinear/newtonmethod.hh>
+#include <dumux/porousmediumflow/compositional/privarswitchnewtoncontroller.hh>
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/methods.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+
+#include "3pwateroilsagdproblem.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.
+ * \ingroup OnePTests
* \param errorMsg The error message that was issued by the start function.
* Comprises the thing that went wrong and a general help message.
*/
@@ -42,21 +64,188 @@ void usage(const char *progName, const std::string &errorMsg)
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.File Name of the file containing the grid \n"
- "\t definition in DGF format\n";
+ "\t-TimeManager.TEnd End of the simulation [s] \n"
+ "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+ "\t-Grid.File Name of the file containing the grid \n"
+ "\t definition in DGF format\n";
std::cout << errorMessageOut
<< "\n";
}
}
-int main(int argc, char** argv)
+int main(int argc, char** argv) try
{
- using ProblemTypeTag = TTAG(ThreePWaterOilSagdBoxProblem);
- return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+ 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);
+
+ // parse command line arguments and input file
+ Parameters::init(argc, argv, usage);
+
+ // 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);
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ const auto maxDivisions = getParam<int>("TimeLoop.MaxTimeStepDivisions");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+ // check if we are about to restart a previously interrupted simulation
+ Scalar restartTime = 0;
+ if (Parameters::getTree().hasKey("Restart") || Parameters::getTree().hasKey("TimeLoop.Restart"))
+ restartTime = getParam<Scalar>("TimeLoop.Restart");
+
+ // intialize the vtk output module
+ using VtkOutputFields = typename GET_PROP_TYPE(TypeTag, VtkOutputFields);
+ VtkOutputModule<TypeTag> vtkWriter(*problem, *fvGridGeometry, *gridVariables, x, problem->name());
+ VtkOutputFields::init(vtkWriter); //! Add model specific output fields
+ vtkWriter.write(0.0);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+ timeLoop->setPeriodicCheckPoint(getParam<Scalar>("TimeLoop.EpisodeLength", std::numeric_limits<Scalar>::max()));
+
+ // 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 = Dumux::AMGBackend<TypeTag>;
+ auto linearSolver = std::make_shared<LinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+ // the non-linear solver
+ using NewtonController = PriVarSwitchNewtonController<TypeTag>;
+ using NewtonMethod = Dumux::NewtonMethod<NewtonController, Assembler, LinearSolver>;
+ auto newtonController = std::make_shared<NewtonController>(leafGridView.comm(), timeLoop);
+ NewtonMethod nonLinearSolver(newtonController, assembler, linearSolver);
+
+ // time loop
+ timeLoop->start();
+ while (!timeLoop->finished())
+ {
+ // 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
+ // if episode length was specificied output only at the end of episodes
+ if (!haveParam("TimeLoop.EpisodeLength") || timeLoop->isCheckPoint() || timeLoop->finished() || timeLoop->timeStepIndex() == 1)
+ 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()));
+ }
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ // print dumux end message
+ if (mpiHelper.rank() == 0)
+ {
+ Parameters::print();
+ 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/3pwateroil/implicit/test_box3pwateroil.input b/test/porousmediumflow/3pwateroil/implicit/test_box3pwateroil.input
index c0656001abb8d2619a893f5604149a952d14b1ac..d1f6789df9cd3161409a4af7bafd7152793e3097 100644
--- a/test/porousmediumflow/3pwateroil/implicit/test_box3pwateroil.input
+++ b/test/porousmediumflow/3pwateroil/implicit/test_box3pwateroil.input
@@ -1,4 +1,4 @@
-[TimeManager]
+[TimeLoop]
DtInitial = 1800 # [s]
TEnd = 60000 # [s]
@@ -13,6 +13,5 @@ Name = sagd # name passed to the output routines
MaxTimeStepDivisions = 100
[Newton]
-RelTolerance = 1e-02
MaxSteps = 8
WriteConvergence = 0