Commit c6fb1dc6 authored by Timo Koch's avatar Timo Koch
Browse files

[porousmflow] Introduce compositional and immisible folders

* The immiscible folder contains files common to immiscible
  porousmediumflow. It now contains a general local residual for np flow

* The compositional folder contains files common to the compositional
  porousmediumflow models. It now contains the privarswitch base class.
parent 8188805b
......@@ -30,10 +30,10 @@
#include "properties.hh"
#include "model.hh"
#include "localresidual.hh"
#include "volumevariables.hh"
#include "indices.hh"
#include <dumux/porousmediumflow/immiscible/localresidual.hh>
#include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
#include <dumux/material/fluidsystems/gasphase.hh>
#include <dumux/material/fluidsystems/liquidphase.hh>
......@@ -54,7 +54,7 @@ SET_INT_PROP(OneP, NumEq, 1); //!< set the number of equations to 1
SET_INT_PROP(OneP, NumPhases, 1); //!< The number of phases in the 1p model is 1
//! The local residual function
SET_TYPE_PROP(OneP, LocalResidual, OnePLocalResidual<TypeTag>);
SET_TYPE_PROP(OneP, LocalResidual, ImmiscibleLocalResidual<TypeTag>);
//! the Model property
SET_TYPE_PROP(OneP, Model, OnePModel<TypeTag>);
......@@ -142,7 +142,7 @@ SET_TYPE_PROP(OnePNI, IsothermalModel, OnePModel<TypeTag>);
SET_TYPE_PROP(OnePNI, IsothermalVolumeVariables, OnePVolumeVariables<TypeTag>);
//set isothermal LocalResidual
SET_TYPE_PROP(OnePNI, IsothermalLocalResidual, OnePLocalResidual<TypeTag>);
SET_TYPE_PROP(OnePNI, IsothermalLocalResidual, ImmiscibleLocalResidual<TypeTag>);
//set isothermal Indices
SET_TYPE_PROP(OnePNI, IsothermalIndices, OnePIndices);
......
......@@ -128,6 +128,12 @@ public:
Scalar pressure(int phaseIdx = 0) const
{ return fluidState_.pressure(phaseIdx); }
/*!
* \brief Return the saturation
*/
Scalar saturation(int phaseIdx = 0) const
{ return 1.0; }
/*!
* \brief Return the mass density \f$\mathrm{[kg/m^3]}\f$ of a given phase within the
* control volume.
......
......@@ -49,7 +49,7 @@ struct TwoPFormulation
*
* \tparam TypeTag The problem type tag
*/
template <class TypeTag>
template <class TypeTag, int PVOffset = 0>
struct TwoPCommonIndices
{
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
......@@ -57,8 +57,18 @@ struct TwoPCommonIndices
// Phase indices
static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< Index of the wetting phase
static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< Index of the non-wetting phase
// Primary variable indices
static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
static const int saturationIdx = PVOffset + 1; //!< Index of the saturation of the non-wetting/wetting phase
// indices of the equations
static const int conti0EqIdx = PVOffset + 0; //!< Index of the first continuity equation
static const int contiWEqIdx = PVOffset + 0; //!< Index of the continuity equation of the wetting phase
static const int contiNEqIdx = PVOffset + 1; //!< Index of the continuity equation of the non-wetting phase
};
/*!
* \ingroup TwoPBoxModel
* \ingroup ImplicitIndices
......@@ -73,19 +83,11 @@ template <class TypeTag,
int formulation = TwoPFormulation::pwsn,
int PVOffset = 0>
struct TwoPIndices
: public TwoPCommonIndices<TypeTag>, TwoPFormulation
: public TwoPCommonIndices<TypeTag, PVOffset>, TwoPFormulation
{
// Primary variable indices
static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
static const int saturationIdx = PVOffset + 1; //!< Index of the saturation of the non-wetting/wetting phase
// indices of the primary variables
static const int pwIdx = PVOffset + 0; //!< index of the wetting phase pressure
static const int snIdx = PVOffset + 1; //!< index of the nonwetting phase saturation
// indices of the equations
static const int contiWEqIdx = PVOffset + 0; //!< Index of the continuity equation of the wetting phase
static const int contiNEqIdx = PVOffset + 1; //!< Index of the continuity equation of the non-wetting phase
};
/*!
......@@ -98,19 +100,11 @@ struct TwoPIndices
*/
template <class TypeTag, int PVOffset>
struct TwoPIndices<TypeTag, TwoPFormulation::pnsw, PVOffset>
: public TwoPCommonIndices<TypeTag>, TwoPFormulation
: public TwoPCommonIndices<TypeTag, PVOffset>, TwoPFormulation
{
// Primary variable indices
static const int pressureIdx = PVOffset + 0; //!< Index for wetting/non-wetting phase pressure (depending on formulation) in a solution vector
static const int saturationIdx = PVOffset + 1; //!< Index of the saturation of the non-wetting/wetting phase
// indices of the primary variables
static const int pnIdx = PVOffset + 0; //!< index of the nonwetting phase pressure
static const int swIdx = PVOffset + 1; //!< index of the wetting phase saturation
// indices of the equations
static const int contiNEqIdx = PVOffset + 0; //!< Index of the continuity equation of the non-wetting phase
static const int contiWEqIdx = PVOffset + 1; //!< Index of the continuity equation of the wetting phase
};
// \}
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief Element-wise calculation of the residual for the two-phase fully implicit model.
*/
#ifndef DUMUX_TWOP_LOCAL_RESIDUAL_BASE_HH
#define DUMUX_TWOP_LOCAL_RESIDUAL_BASE_HH
#include "properties.hh"
namespace Dumux
{
/*!
* \ingroup TwoPModel
* \ingroup ImplicitLocalResidual
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the two-phase fully implicit model.
*
* This class is also used for the non-isothermal model, which means
* that it uses static polymorphism.
*/
template<class TypeTag>
class TwoPLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
{
protected:
typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum
{
contiWEqIdx = Indices::contiWEqIdx,
contiNEqIdx = Indices::contiNEqIdx,
wPhaseIdx = Indices::wPhaseIdx,
nPhaseIdx = Indices::nPhaseIdx,
numPhases = GET_PROP_VALUE(TypeTag, NumPhases)
};
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
/*!
* \brief Constructor
*
* Sets the upwind weight.
*/
TwoPLocalResidual()
{
// retrieve the upwind weight for the mass conservation equations. Use the value
// specified via the property system as default, and overwrite
// it by the run-time parameter from the Dune::ParameterTree
massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
};
/*!
* \brief Evaluate the amount of all conservation quantities
* (e.g. phase mass) within a sub-control volume.
*
* \param storage The phase mass 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
*
* The result should be averaged over the volume (e.g. phase mass
* inside a sub-control volume divided by the volume)
*/
PrimaryVariables computeStorage(const SubControlVolume& scv,
const VolumeVariables& volVars) const
{
PrimaryVariables storage(0);
// wetting phase mass
storage[contiWEqIdx] = volVars.density(wPhaseIdx) * volVars.porosity()
* volVars.saturation(wPhaseIdx);
// non-wetting phase mass
storage[contiNEqIdx] = volVars.density(nPhaseIdx) * volVars.porosity()
* volVars.saturation(nPhaseIdx);
return storage;
}
/*!
* \brief Evaluates the total flux of all conservation quantities
* over a face of a sub-control volume.
*
* \param flux The flux over the sub-control-volume face for each component
* \param fIdx The index of the sub-control-volume face
* \param onBoundary Evaluate flux at inner sub-control-volume face or on a boundary face
*/
PrimaryVariables computeFlux(const SubControlVolumeFace& scvFace) const
{
const auto& fluxVars = this->model_().fluxVars(scvFace);
auto massWeight = massUpwindWeight_;
auto mobWeight = mobilityUpwindWeight_;
// loop over all phases
PrimaryVariables flux;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
auto upwindRule = [massWeight, mobWeight, phaseIdx](const VolumeVariables& up, const VolumeVariables& dn)
{ return (up.density(phaseIdx)*massWeight + dn.density(phaseIdx)*(1-massWeight))
*(up.mobility(phaseIdx)*mobWeight + dn.mobility(phaseIdx)*(1-mobWeight)); };
auto eqIdx = (phaseIdx == wPhaseIdx) ? contiWEqIdx : contiNEqIdx;
flux[eqIdx] = fluxVars.advection().flux(phaseIdx, upwindRule);
}
return flux;
}
protected:
Implementation *asImp_()
{
return static_cast<Implementation *> (this);
}
const Implementation *asImp_() const
{
return static_cast<const Implementation *> (this);
}
private:
Scalar massUpwindWeight_;
Scalar mobilityUpwindWeight_;
};
}
#endif
......@@ -32,8 +32,8 @@
#include "model.hh"
#include "indices.hh"
#include "volumevariables.hh"
#include "localresidual.hh"
#include <dumux/porousmediumflow/immiscible/localresidual.hh>
#include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
#include <dumux/material/fluidsystems/gasphase.hh>
#include <dumux/material/fluidsystems/liquidphase.hh>
......@@ -64,7 +64,7 @@ SET_INT_PROP(TwoP,
//! Use the 2p local jacobian operator for the 2p model
SET_TYPE_PROP(TwoP,
LocalResidual,
TwoPLocalResidual<TypeTag>);
ImmiscibleLocalResidual<TypeTag>);
//! the Model property
SET_TYPE_PROP(TwoP, Model, TwoPModel<TypeTag>);
......@@ -176,7 +176,7 @@ SET_TYPE_PROP(TwoPNI, IsothermalModel, TwoPModel<TypeTag>);
SET_TYPE_PROP(TwoPNI, IsothermalVolumeVariables, TwoPVolumeVariables<TypeTag>);
//set isothermal LocalResidual
SET_TYPE_PROP(TwoPNI, IsothermalLocalResidual, TwoPLocalResidual<TypeTag>);
SET_TYPE_PROP(TwoPNI, IsothermalLocalResidual, ImmiscibleLocalResidual<TypeTag>);
//set isothermal Indices
SET_PROP(TwoPNI, IsothermalIndices)
......
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 2 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
*
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the three-phase fully implicit model.
*/
#ifndef DUMUX_3P_LOCAL_RESIDUAL_HH
#define DUMUX_3P_LOCAL_RESIDUAL_HH
#include "properties.hh"
namespace Dumux
{
/*!
* \ingroup ThreePModel
* \ingroup ImplicitLocalResidual
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the three-phase fully implicit model.
*
* This class is used to fill the gaps in BoxLocalResidual for three-phase flow.
*/
template<class TypeTag>
class ThreePLocalResidual: public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
{
protected:
typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
typedef typename GET_PROP_TYPE(TypeTag, BaseLocalResidual) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolume) SubControlVolume;
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum {
numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
conti0EqIdx = Indices::conti0EqIdx //!< Index of the mass conservation equation for the water component
};
typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
public:
ThreePLocalResidual() : ParentType()
{
massWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
}
/*!
* \brief Evaluate the amount of all conservation quantities
* (e.g. phase mass) within a sub-control volume.
*
* The result should be averaged over the volume (e.g. phase mass
* inside a sub-control volume divided by the volume).
*
* \param storage The mass of the component within the sub-control volume
* \param scvIdx The SCV (sub-control volume) index
* \param usePrevSol Evaluate function with solution of current or previous time step
*/
PrimaryVariables computeStorage(const SubControlVolume &scv,
const VolumeVariables &volVars) const
{
PrimaryVariables storage(0.0);
// compute storage term of all components within all phases
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
auto eqIdx = conti0EqIdx + phaseIdx;
storage[eqIdx] +=
volVars.porosity()
* volVars.saturation(phaseIdx)
* volVars.density(phaseIdx);
}
return storage;
}
/*!
* \brief Evaluates the total flux of all conservation quantities
* over a face of a sub-control volume.
*
* \param flux The flux over the SCV (sub-control-volume) face for each component
* \param 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
*/
PrimaryVariables computeFlux(const SubControlVolumeFace& scvFace) const
{
auto& fluxVars = this->model_().fluxVars(scvFace);
// get upwind weights into local scope
auto massWeight = massWeight_;
PrimaryVariables flux(0.0);
// advective fluxes of all phases
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
auto upwindRule = [massWeight, phaseIdx](const VolumeVariables& up, const VolumeVariables& dn)
{
return (massWeight)*up.mobility(phaseIdx)*up.density(phaseIdx)
+(1-massWeight)*dn.mobility(phaseIdx)*dn.density(phaseIdx);
};
// add advective flux of current phase
auto eqIdx = conti0EqIdx + phaseIdx;
flux[eqIdx] += fluxVars.advection().flux(phaseIdx, upwindRule);
}
return flux;
}
protected:
Implementation *asImp_()
{
return static_cast<Implementation *> (this);
}
const Implementation *asImp_() const
{
return static_cast<const Implementation *> (this);
}
private:
Scalar massWeight_;
};
} // end namespace
#endif
......@@ -34,8 +34,8 @@
#include "indices.hh"
#include "volumevariables.hh"
#include "properties.hh"
#include "localresidual.hh"
#include <dumux/porousmediumflow/immiscible/localresidual.hh>
#include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
#include <dumux/porousmediumflow/implicit/darcyfluxvariables.hh>
#include <dumux/material/spatialparams/implicit.hh>
......@@ -78,7 +78,7 @@ SET_INT_PROP(ThreeP, NumEq, 3); //!< set the number of equations to 2
SET_TYPE_PROP(ThreeP, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
//! The local residual function of the conservation equations
SET_TYPE_PROP(ThreeP, LocalResidual, ThreePLocalResidual<TypeTag>);
SET_TYPE_PROP(ThreeP, LocalResidual, ImmiscibleLocalResidual<TypeTag>);
//! Enable advection
SET_BOOL_PROP(ThreeP, EnableAdvection, true);
......@@ -147,7 +147,7 @@ SET_TYPE_PROP(ThreePNI, IsothermalModel, ThreePModel<TypeTag>);
SET_TYPE_PROP(ThreePNI, IsothermalVolumeVariables, ThreePVolumeVariables<TypeTag>);
//set isothermal LocalResidual
SET_TYPE_PROP(ThreePNI, IsothermalLocalResidual, ThreePLocalResidual<TypeTag>);
SET_TYPE_PROP(ThreePNI, IsothermalLocalResidual, ImmiscibleLocalResidual<TypeTag>);
//set isothermal Indices
SET_TYPE_PROP(ThreePNI, IsothermalIndices, ThreePIndices<TypeTag,/*PVOffset=*/0>);
......
......@@ -24,7 +24,7 @@
#ifndef DUMUX_3P3C_PRIMARY_VARIABLE_SWITCH_HH
#define DUMUX_3P3C_PRIMARY_VARIABLE_SWITCH_HH
#include <dumux/porousmediumflow/implicit/primaryvariableswitch.hh>
#include <dumux/porousmediumflow/compositional/primaryvariableswitch.hh>
namespace Dumux
{
......
......@@ -20,12 +20,10 @@
* \file
*
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the one-phase fully implicit model.
* using the n-phase immiscible fully implicit models.
*/
#ifndef DUMUX_1P_LOCAL_RESIDUAL_HH
#define DUMUX_1P_LOCAL_RESIDUAL_HH
#include "properties.hh"
#ifndef DUMUX_IMMISCIBLE_LOCAL_RESIDUAL_HH
#define DUMUX_IMMISCIBLE_LOCAL_RESIDUAL_HH
namespace Dumux
{
......@@ -33,10 +31,10 @@ namespace Dumux
* \ingroup OnePModel
* \ingroup ImplicitLocalResidual
* \brief Element-wise calculation of the Jacobian matrix for problems
* using the one-phase fully implicit model.
* using the n-phase immiscible fully implicit models.
*/
template<class TypeTag>
class OnePLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
class ImmiscibleLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
{
typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
......@@ -48,81 +46,78 @@ class OnePLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
typedef typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace) SubControlVolumeFace;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
//index of the mass balance equation
enum {
conti0EqIdx = Indices::conti0EqIdx //index for the mass balance
};
// first index for the mass balance
enum { conti0EqIdx = Indices::conti0EqIdx };
static const int numPhases = GET_PROP_VALUE(TypeTag, NumPhases);
public:
/*!
* \brief Constructor. Sets the upwind weight.
* \brief Constructor. Gets the upwind weight.
*/
OnePLocalResidual()
ImmiscibleLocalResidual()
{
// retrieve the upwind weight for the mass conservation equations. Use the value
// specified via the property system as default, and overwrite
// it by the run-time parameter from the Dune::ParameterTree
massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
upwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
}
/*!
* \brief Evaluate the rate of change of all conservation
* quantites (e.g. phase mass) within a sub-control
* volume of a finite volume element for the OneP
* model.
*
* This function should not include the source and sink terms.
* \param scv The sub control volume
* \param volVars The current or previous volVars
*
* The volVars can be different to allow computing the implicit euler time derivative here
* volume of a finite volume element for the immiscible models.
* \param scv The sub control volume
* \param volVars The current or previous volVars
* \note This function should not include the source and sink terms.
* \note The volVars can be different to allow computing
* the implicit euler time derivative here
*/
PrimaryVariables computeStorage(const SubControlVolume& scv,
const VolumeVariables& volVars) const
{
PrimaryVariables storage(0);
// partial time derivative of the wetting phase mass
storage[conti0EqIdx] = volVars.density() * volVars.porosity();
// partial time derivative of the phase mass
PrimaryVariables storage;
for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
{
auto eqIdx = conti0EqIdx + phaseIdx;
storage[eqIdx] = volVars.porosity()
* volVars.density(phaseIdx)
* volVars.saturation(phaseIdx);
}
return storage;
}
/*!
* \brief Evaluate the mass flux over a face of a sub-control
* volume.
*
* \param flux The flux over the SCV (sub-control-volume) face
* \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
* \brief Evaluate the mass flux over a face of a sub control volume
* \param scvf The sub control volume face to compute the flux on
*/
PrimaryVariables computeFlux(const SubControlVolumeFace& scvFace) const
PrimaryVariables computeFlux(const SubControlVolumeFace& scvf)
{
const auto& fluxVars = this->model_().fluxVars(scvFace);
auto& fluxVars = this->model_().fluxVars_(scvf);
fluxVars.beginFluxComputation();
auto massWeight = massUpwindWeight_;
auto mobWeight = mobilityUpwindWeight_;
// copy weight to local scope for use in lambda expression
auto w = upwindWeight_;