OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-01-12 01:11:55 -06:00
Code Issues Packages Projects Releases Wiki Activity

changed: ewoms/models/richards -> opm/models/richards

Browse Source
This commit is contained in:
Arne Morten Kvarving 2019-09-16 12:17:18 +02:00
parent f558f5d98b
commit 88a5e1db06
11 changed files with 1515 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/problems/richardslensproblem.hh
View File

@ -28,7 +28,7 @@
#ifndef EWOMS_RICHARDS_LENS_PROBLEM_HH
#define EWOMS_RICHARDS_LENS_PROBLEM_HH
#include <ewoms/models/richards/richardsmodel.hh>
#include <opm/models/richards/richardsmodel.hh>
#include <opm/material/components/SimpleH2O.hpp>
#include <opm/material/fluidsystems/LiquidPhase.hpp>

165
opm/models/richards/richardsboundaryratevector.hh Normal file
View File

@ -0,0 +1,165 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::RichardsBoundaryRateVector
*/
#ifndef EWOMS_RICHARDS_BOUNDARY_RATE_VECTOR_HH
#define EWOMS_RICHARDS_BOUNDARY_RATE_VECTOR_HH
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/constraintsolvers/NcpFlash.hpp>
#include "richardsintensivequantities.hh"
namespace Opm {
/*!
* \ingroup RichardsModel
*
* \brief Implements a boundary vector for the fully implicit Richards model.
*/
template <class TypeTag>
class RichardsBoundaryRateVector : public GET_PROP_TYPE(TypeTag, RateVector)
{
typedef typename GET_PROP_TYPE(TypeTag, RateVector) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
enum { contiEqIdx = Indices::contiEqIdx };
enum { liquidPhaseIdx = GET_PROP_VALUE(TypeTag, LiquidPhaseIndex) };
typedef Opm::MathToolbox<Evaluation> Toolbox;
public:
RichardsBoundaryRateVector() : ParentType()
{}
/*!
* \copydoc
* ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(Scalar)
*/
RichardsBoundaryRateVector(const Evaluation& value)
: ParentType(value)
{}
/*!
* \copydoc ImmiscibleBoundaryRateVector::ImmiscibleBoundaryRateVector(const
* ImmiscibleBoundaryRateVector& )
*/
RichardsBoundaryRateVector(const RichardsBoundaryRateVector& value) = default;
RichardsBoundaryRateVector& operator=(const RichardsBoundaryRateVector& value) = default;
/*!
* \copydoc ImmiscibleBoundaryRateVector::setFreeFlow
*/
template <class Context, class FluidState>
void setFreeFlow(const Context& context, unsigned bfIdx, unsigned timeIdx, const FluidState& fluidState)
{
ExtensiveQuantities extQuants;
extQuants.updateBoundary(context, bfIdx, timeIdx, fluidState);
const auto& insideIntQuants = context.intensiveQuantities(bfIdx, timeIdx);
unsigned focusDofIdx = context.focusDofIndex();
unsigned interiorDofIdx = context.interiorScvIndex(bfIdx, timeIdx);
////////
// advective fluxes of all components in all phases
////////
(*this) = Evaluation(0.0);
unsigned phaseIdx = liquidPhaseIdx;
Evaluation density;
if (fluidState.pressure(phaseIdx) > insideIntQuants.fluidState().pressure(phaseIdx)) {
if (focusDofIdx == interiorDofIdx)
density = fluidState.density(phaseIdx);
else
density = Opm::getValue(fluidState.density(phaseIdx));
}
else if (focusDofIdx == interiorDofIdx)
density = insideIntQuants.fluidState().density(phaseIdx);
else
density = Opm::getValue(insideIntQuants.fluidState().density(phaseIdx));
// add advective flux of current component in current
// phase
(*this)[contiEqIdx] += extQuants.volumeFlux(phaseIdx) * density;
#ifndef NDEBUG
for (unsigned i = 0; i < numEq; ++i) {
Opm::Valgrind::CheckDefined((*this)[i]);
}
Opm::Valgrind::CheckDefined(*this);
#endif
}
/*!
* \copydoc ImmiscibleBoundaryRateVector::setInFlow
*/
template <class Context, class FluidState>
void setInFlow(const Context& context,
unsigned bfIdx,
unsigned timeIdx,
const FluidState& fluidState)
{
this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
// we only allow fluxes in the direction opposite to the outer unit normal
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
Evaluation& val = this->operator[](eqIdx);
val = Toolbox::min(0.0, val);
}
}
/*!
* \copydoc ImmiscibleBoundaryRateVector::setOutFlow
*/
template <class Context, class FluidState>
void setOutFlow(const Context& context,
unsigned bfIdx,
unsigned timeIdx,
const FluidState& fluidState)
{
this->setFreeFlow(context, bfIdx, timeIdx, fluidState);
// we only allow fluxes in the same direction as the outer unit normal
for (unsigned eqIdx = 0; eqIdx < numEq; ++eqIdx) {
Evaluation& val = this->operator[](eqIdx);
val = Toolbox::max(0.0, val);
}
}
/*!
* \copydoc ImmiscibleBoundaryRateVector::setNoFlow
*/
void setNoFlow()
{ (*this) = Evaluation(0.0); }
};
} // namespace Opm
#endif

52
opm/models/richards/richardsextensivequantities.hh Normal file
View File

@ -0,0 +1,52 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::RichardsExtensiveQuantities
*/
#ifndef EWOMS_RICHARDS_EXTENSIVE_QUANTITIES_HH
#define EWOMS_RICHARDS_EXTENSIVE_QUANTITIES_HH
#include "richardsproperties.hh"
#include <opm/models/common/multiphasebaseextensivequantities.hh>
namespace Opm {
/*!
* \ingroup RichardsModel
* \ingroup ExtensiveQuantities
*
* \brief Calculates and stores the data which is required to
* calculate the flux of fluid over a face of a finite volume.
*/
template <class TypeTag>
class RichardsExtensiveQuantities
: public MultiPhaseBaseExtensiveQuantities<TypeTag>
{
};
} // namespace Opm
#endif

52
opm/models/richards/richardsindices.hh Normal file
View File

@ -0,0 +1,52 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::RichardsIndices
*/
#ifndef EWOMS_RICHARDS_INDICES_HH
#define EWOMS_RICHARDS_INDICES_HH
namespace Opm {
/*!
* \ingroup RichardsModel
* \brief Indices for the primary variables/conservation equations of the
* Richards model.
*/
struct RichardsIndices
{
//! Primary variable index for the wetting phase pressure
static const int pressureWIdx = 0;
//! Equation index for the mass conservation of the wetting phase
static const int contiEqIdx = 0;
//! The number of equations
static const int numEq = 1;
};
} // namespace Opm
#endif

199
opm/models/richards/richardsintensivequantities.hh Normal file
View File

@ -0,0 +1,199 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::RichardsIntensiveQuantities
*/
#ifndef EWOMS_RICHARDS_INTENSIVE_QUANTITIES_HH
#define EWOMS_RICHARDS_INTENSIVE_QUANTITIES_HH
#include "richardsproperties.hh"
#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
#include <dune/common/fvector.hh>
#include <dune/common/fmatrix.hh>
namespace Opm {
/*!
* \ingroup RichardsModel
* \ingroup IntensiveQuantities
*
* \brief Intensive quantities required by the Richards model.
*/
template <class TypeTag>
class RichardsIntensiveQuantities
: public GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities)
, public GET_PROP_TYPE(TypeTag, FluxModule)::FluxIntensiveQuantities
{
typedef typename GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
typedef typename GET_PROP_TYPE(TypeTag, FluxModule) FluxModule;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum { pressureWIdx = Indices::pressureWIdx };
enum { numPhases = FluidSystem::numPhases };
enum { liquidPhaseIdx = GET_PROP_VALUE(TypeTag, LiquidPhaseIndex) };
enum { gasPhaseIdx = GET_PROP_VALUE(TypeTag, GasPhaseIndex) };
enum { dimWorld = GridView::dimensionworld };
typedef typename FluxModule::FluxIntensiveQuantities FluxIntensiveQuantities;
typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
typedef Dune::FieldVector<Scalar, numPhases> ScalarPhaseVector;
typedef Dune::FieldVector<Evaluation, numPhases> PhaseVector;
typedef Opm::MathToolbox<Evaluation> Toolbox;
public:
//! The type returned by the fluidState() method
typedef Opm::ImmiscibleFluidState<Evaluation, FluidSystem> FluidState;
RichardsIntensiveQuantities()
{}
RichardsIntensiveQuantities(const RichardsIntensiveQuantities& other) = default;
RichardsIntensiveQuantities& operator=(const RichardsIntensiveQuantities& other) = default;
/*!
* \copydoc IntensiveQuantities::update
*/
void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
{
ParentType::update(elemCtx, dofIdx, timeIdx);
const auto& T = elemCtx.problem().temperature(elemCtx, dofIdx, timeIdx);
fluidState_.setTemperature(T);
// material law parameters
const auto& problem = elemCtx.problem();
const typename MaterialLaw::Params& materialParams =
problem.materialLawParams(elemCtx, dofIdx, timeIdx);
const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
/////////
// calculate the pressures
/////////
// first, we have to find the minimum capillary pressure (i.e. Sw = 0)
fluidState_.setSaturation(liquidPhaseIdx, 1.0);
fluidState_.setSaturation(gasPhaseIdx, 0.0);
ScalarPhaseVector pC;
MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
// non-wetting pressure can be larger than the
// reference pressure if the medium is fully
// saturated by the wetting phase
const Evaluation& pW = priVars.makeEvaluation(pressureWIdx, timeIdx);
Evaluation pN =
Toolbox::max(elemCtx.problem().referencePressure(elemCtx, dofIdx, /*timeIdx=*/0),
pW + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
/////////
// calculate the saturations
/////////
fluidState_.setPressure(liquidPhaseIdx, pW);
fluidState_.setPressure(gasPhaseIdx, pN);
PhaseVector sat;
MaterialLaw::saturations(sat, materialParams, fluidState_);
fluidState_.setSaturation(liquidPhaseIdx, sat[liquidPhaseIdx]);
fluidState_.setSaturation(gasPhaseIdx, sat[gasPhaseIdx]);
typename FluidSystem::template ParameterCache<Evaluation> paramCache;
paramCache.updateAll(fluidState_);
// compute and set the wetting phase viscosity
const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, liquidPhaseIdx);
fluidState_.setViscosity(liquidPhaseIdx, mu);
fluidState_.setViscosity(gasPhaseIdx, 1e-20);
// compute and set the wetting phase density
const Evaluation& rho = FluidSystem::density(fluidState_, paramCache, liquidPhaseIdx);
fluidState_.setDensity(liquidPhaseIdx, rho);
fluidState_.setDensity(gasPhaseIdx, 1e-20);
// relperms
MaterialLaw::relativePermeabilities(relativePermeability_, materialParams, fluidState_);
// mobilities
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
mobility_[phaseIdx] = relativePermeability_[phaseIdx]/fluidState_.viscosity(phaseIdx);
// porosity
porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
// intrinsic permeability
intrinsicPerm_ = problem.intrinsicPermeability(elemCtx, dofIdx, timeIdx);
// update the quantities specific for the velocity model
FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
}
/*!
* \copydoc ImmiscibleIntensiveQuantities::fluidState
*/
const FluidState& fluidState() const
{ return fluidState_; }
/*!
* \copydoc ImmiscibleIntensiveQuantities::porosity
*/
const Evaluation& porosity() const
{ return porosity_; }
/*!
* \copydoc ImmiscibleIntensiveQuantities::intrinsicPermeability
*/
const DimMatrix& intrinsicPermeability() const
{ return intrinsicPerm_; }
/*!
* \copydoc ImmiscibleIntensiveQuantities::relativePermeability
*/
const Evaluation& relativePermeability(unsigned phaseIdx) const
{ return relativePermeability_[phaseIdx]; }
/*!
* \copydoc ImmiscibleIntensiveQuantities::mobility
*/
const Evaluation& mobility(unsigned phaseIdx) const
{ return mobility_[phaseIdx]; }
private:
FluidState fluidState_;
DimMatrix intrinsicPerm_;
Evaluation relativePermeability_[numPhases];
Evaluation mobility_[numPhases];
Evaluation porosity_;
};
} // namespace Opm
#endif

111
opm/models/richards/richardslocalresidual.hh Normal file
View File

@ -0,0 +1,111 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::RichardsLocalResidual
*/
#ifndef EWOMS_RICHARDS_LOCAL_RESIDUAL_HH
#define EWOMS_RICHARDS_LOCAL_RESIDUAL_HH
#include "richardsintensivequantities.hh"
#include "richardsextensivequantities.hh"
namespace Opm {
/*!
* \ingroup RichardsModel
* \brief Element-wise calculation of the residual for the Richards model.
*/
template <class TypeTag>
class RichardsLocalResidual : public GET_PROP_TYPE(TypeTag, DiscLocalResidual)
{
typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum { contiEqIdx = Indices::contiEqIdx };
enum { liquidPhaseIdx = GET_PROP_VALUE(TypeTag, LiquidPhaseIndex) };
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
typedef Opm::MathToolbox<Evaluation> Toolbox;
public:
/*!
* \copydoc ImmiscibleLocalResidual::computeStorage
*/
template <class LhsEval>
void computeStorage(Dune::FieldVector<LhsEval, numEq>& storage,
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx) const
{
const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
// partial time derivative of the wetting phase mass
storage[contiEqIdx] =
Toolbox::template decay<LhsEval>(intQuants.fluidState().density(liquidPhaseIdx))
*Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(liquidPhaseIdx))
*Toolbox::template decay<LhsEval>(intQuants.porosity());
}
/*!
* \copydoc ImmiscibleLocalResidual::computeFlux
*/
void computeFlux(RateVector& flux,
const ElementContext& elemCtx,
unsigned scvfIdx,
unsigned timeIdx) const
{
const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
unsigned focusDofIdx = elemCtx.focusDofIndex();
unsigned upIdx = static_cast<unsigned>(extQuants.upstreamIndex(liquidPhaseIdx));
const IntensiveQuantities& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
// compute advective mass flux of the liquid phase. This is slightly hacky
// because it is specific to the element-centered finite volume method.
const Evaluation& rho = up.fluidState().density(liquidPhaseIdx);
if (focusDofIdx == upIdx)
flux[contiEqIdx] = extQuants.volumeFlux(liquidPhaseIdx)*rho;
else
flux[contiEqIdx] = extQuants.volumeFlux(liquidPhaseIdx)*Toolbox::value(rho);
}
/*!
* \copydoc ImmiscibleLocalResidual::computeSource
*/
void computeSource(RateVector& source,
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx) const
{ elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx); }
};
} // namespace Opm
#endif

373
opm/models/richards/richardsmodel.hh Normal file
View File

@ -0,0 +1,373 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::RichardsModel
*/
#ifndef EWOMS_RICHARDS_MODEL_HH
#define EWOMS_RICHARDS_MODEL_HH
#include <opm/material/densead/Math.hpp>
#include "richardsproperties.hh"
#include "richardsindices.hh"
#include "richardslocalresidual.hh"
#include "richardsextensivequantities.hh"
#include "richardsratevector.hh"
#include "richardsboundaryratevector.hh"
#include "richardsprimaryvariables.hh"
#include "richardsintensivequantities.hh"
#include "richardsnewtonmethod.hh"
#include <opm/models/common/multiphasebasemodel.hh>
#include <opm/material/components/NullComponent.hpp>
#include <opm/material/fluidsystems/LiquidPhase.hpp>
#include <opm/material/fluidsystems/GasPhase.hpp>
#include <opm/material/fluidsystems/TwoPhaseImmiscibleFluidSystem.hpp>
#include <opm/material/common/Unused.hpp>
#include <sstream>
#include <string>
namespace Opm {
template <class TypeTag>
class RichardsModel;
}
BEGIN_PROPERTIES
//! The type tag for problems discretized using the Richards model
NEW_TYPE_TAG(Richards, INHERITS_FROM(MultiPhaseBaseModel));
//! By default, assume that the first phase is the liquid one
SET_INT_PROP(Richards, LiquidPhaseIndex, 0);
//! By default, assume that the non-liquid phase is gaseos
SET_INT_PROP(Richards, GasPhaseIndex, 1 - GET_PROP_VALUE(TypeTag, LiquidPhaseIndex));
/*!
* \brief By default, assume that component which the liquid is made of has
* the same index as the liquid phase.
*
* This is a convention which works for most fluid systems shipped
* with eWoms by default, but it cannot generally correct because the
* liquid can be composed of different components. (e.g., do you
* prefer Ethanol of H2O??)
*/
SET_INT_PROP(Richards, LiquidComponentIndex, GET_PROP_VALUE(TypeTag, LiquidPhaseIndex));
//! By default, assume that the gas component is the other than the liquid one
SET_INT_PROP(Richards, GasComponentIndex, 1 - GET_PROP_VALUE(TypeTag, LiquidComponentIndex));
//! The local residual operator
SET_TYPE_PROP(Richards,
LocalResidual,
Opm::RichardsLocalResidual<TypeTag>);
//! The global model used
SET_TYPE_PROP(Richards, Model, Opm::RichardsModel<TypeTag>);
//! the RateVector property
SET_TYPE_PROP(Richards, RateVector, Opm::RichardsRateVector<TypeTag>);
//! the BoundaryRateVector property
SET_TYPE_PROP(Richards, BoundaryRateVector, Opm::RichardsBoundaryRateVector<TypeTag>);
//! the PrimaryVariables property
SET_TYPE_PROP(Richards, PrimaryVariables, Opm::RichardsPrimaryVariables<TypeTag>);
//! The class for the intensive quantities
SET_TYPE_PROP(Richards, IntensiveQuantities, Opm::RichardsIntensiveQuantities<TypeTag>);
//! The class for the quantities required for the flux calculation
SET_TYPE_PROP(Richards, ExtensiveQuantities, Opm::RichardsExtensiveQuantities<TypeTag>);
//! The class of the Newton method
SET_TYPE_PROP(Richards, NewtonMethod, Opm::RichardsNewtonMethod<TypeTag>);
//! The class with all index definitions for the model
SET_TYPE_PROP(Richards, Indices, Opm::RichardsIndices);
/*!
* \brief The wetting phase used.
*
* By default we use the null-phase, i.e. this has to be defined by
* the problem for the program to work. Please be aware that you
* should be careful to use the Richards model in conjunction with
* liquid non-wetting phases. This is only meaningful if the viscosity
* of the liquid phase is _much_ lower than the viscosity of the
* wetting phase.
*/
SET_PROP(Richards, WettingFluid)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef Opm::LiquidPhase<Scalar, Opm::NullComponent<Scalar> > type;
};
/*!
* \brief The non-wetting phase used.
*
* By default we use the null-phase, i.e. this has to be defined by
* the problem for the program to work. This doed not need to be
* specified by the problem for the Richards model to work because the
* Richards model does not conserve the non-wetting phase.
*/
SET_PROP(Richards, NonWettingFluid)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
public:
typedef Opm::GasPhase<Scalar, Opm::NullComponent<Scalar> > type;
};
/*!
*\brief The fluid system used by the model.
*
* By default this uses the immiscible twophase fluid system. The
* actual fluids used are specified using in the problem definition by
* the WettingFluid and NonWettingFluid properties. Be aware that
* using different fluid systems in conjunction with the Richards
* model only makes very limited sense.
*/
SET_PROP(Richards, FluidSystem)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, WettingFluid) WettingFluid;
typedef typename GET_PROP_TYPE(TypeTag, NonWettingFluid) NonWettingFluid;
public:
typedef Opm::TwoPhaseImmiscibleFluidSystem<Scalar, WettingFluid, NonWettingFluid> type;
};
END_PROPERTIES
namespace Opm {
/*!
* \ingroup RichardsModel
*
* \brief This model implements a variant of the Richards equation for
* quasi-twophase flow.
*
* In the unsaturated zone, Richards' equation is frequently used to
* approximate the water distribution above the groundwater level. It
* can be derived from the two-phase equations, i.e.
* \f[
* \frac{\partial\;\phi S_\alpha \rho_\alpha}{\partial t}
* -
* \mathrm{div} \left\{
* \rho_\alpha \frac{k_{r\alpha}}{\mu_\alpha}\; \mathbf{K}\;
* \mathbf{grad}\left[
* p_\alpha - g\rho_\alpha
* \right]
* \right\}
* =
* q_\alpha,
* \f]
* where \f$\alpha \in \{w, n\}\f$ is the index of the fluid phase,
* \f$\rho_\alpha\f$ is the fluid density, \f$S_\alpha\f$ is the fluid
* saturation, \f$\phi\f$ is the porosity of the soil,
* \f$k_{r\alpha}\f$ is the relative permeability for the fluid,
* \f$\mu_\alpha\f$ is the fluid's dynamic viscosity, \f$\mathbf{K}\f$
* is the intrinsic permeability tensor, \f$p_\alpha\f$ is the fluid
* phase pressure and \f$g\f$ is the potential of the gravity field.
*
* In contrast to the "full" two-phase model, the Richards model
* assumes that the non-wetting fluid is gas and that it thus exhibits
* a much lower viscosity than the (liquid) wetting phase. (This
* assumption is quite realistic in many applications: For example, at
* atmospheric pressure and at room temperature, the viscosity of air
* is only about \f$1\%\f$ of the viscosity of liquid water.) As a
* consequence, the \f$\frac{k_{r\alpha}}{\mu_\alpha}\f$ term
* typically is much larger for the gas phase than for the wetting
* phase. Using this reasoning, the Richards model assumes that
* \f$\frac{k_{rn}}{\mu_n}\f$ is infinitely large compared to the same
* term of the liquid phase. This implies that the pressure of the gas
* phase is equivalent to the static pressure distribution and that
* therefore, mass conservation only needs to be considered for the
* liquid phase.
*
* The model thus choses the absolute pressure of the wetting phase
* \f$p_w\f$ as its only primary variable. The wetting phase
* saturation is calculated using the inverse of the capillary
* pressure, i.e.
* \f[
* S_w = p_c^{-1}(p_n - p_w)
* \f]
* holds, where \f$p_n\f$ is a reference pressure given by the
* problem's \c referencePressure() method. Nota bene, that the last
* step assumes that the capillary pressure-saturation curve can be
* uniquely inverted, i.e. it is not possible to set the capillary
* pressure to zero if the Richards model ought to be used!
*/
template <class TypeTag>
class RichardsModel
: public MultiPhaseBaseModel<TypeTag>
{
typedef MultiPhaseBaseModel<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
static const unsigned numPhases = FluidSystem::numPhases;
static const unsigned numComponents = FluidSystem::numComponents;
static const unsigned liquidPhaseIdx = GET_PROP_VALUE(TypeTag, LiquidPhaseIndex);
static const unsigned gasPhaseIdx = GET_PROP_VALUE(TypeTag, GasPhaseIndex);
static const unsigned liquidCompIdx = GET_PROP_VALUE(TypeTag, LiquidComponentIndex);
static const unsigned gasCompIdx = GET_PROP_VALUE(TypeTag, GasComponentIndex);
// some consistency checks
static_assert(numPhases == 2,
"Exactly two fluids are required for this model");
static_assert(numComponents == 2,
"Exactly two components are required for this model");
static_assert(liquidPhaseIdx != gasPhaseIdx,
"The liquid and the gas phases must be different");
static_assert(liquidCompIdx != gasCompIdx,
"The liquid and the gas components must be different");
public:
RichardsModel(Simulator& simulator)
: ParentType(simulator)
{
// the liquid phase must be liquid, the gas phase must be
// gaseous. Think about it!
assert(FluidSystem::isLiquid(liquidPhaseIdx));
assert(!FluidSystem::isLiquid(gasPhaseIdx));
}
/*!
* \copydoc FvBaseDiscretization::registerParameters
*/
static void registerParameters()
{
ParentType::registerParameters();
}
/*!
* \copydoc FvBaseDiscretization::name
*/
static std::string name()
{ return "richards"; }
/*!
* \copydoc FvBaseDiscretization::primaryVarName
*/
std::string primaryVarName(unsigned pvIdx) const
{
std::ostringstream oss;
if (pvIdx == Indices::pressureWIdx)
oss << "pressure_" << FluidSystem::phaseName(liquidPhaseIdx);
else
assert(0);
return oss.str();
}
/*!
* \copydoc FvBaseDiscretization::eqName
*/
std::string eqName(unsigned eqIdx) const
{
std::ostringstream oss;
if (eqIdx == Indices::contiEqIdx)
oss << "continuity_" << FluidSystem::phaseName(liquidPhaseIdx);
else
assert(0);
return oss.str();
}
/*!
* \copydoc FvBaseDiscretization::primaryVarWeight
*/
Scalar primaryVarWeight(unsigned globalDofIdx OPM_UNUSED, unsigned pvIdx) const
{
if (Indices::pressureWIdx == pvIdx) {
return 10 / referencePressure_;
}
return 1;
}
/*!
* \copydoc FvBaseDiscretization::eqWeight
*/
Scalar eqWeight(unsigned globalDofIdx OPM_UNUSED, unsigned OPM_OPTIM_UNUSED eqIdx) const
{
unsigned OPM_OPTIM_UNUSED compIdx = eqIdx - Indices::contiEqIdx;
assert(0 <= compIdx && compIdx <= FluidSystem::numPhases);
// make all kg equal
return 1.0;
}
/*!
* \copydoc FvBaseDiscretization::updateBegin
*/
void updateBegin()
{
ParentType::updateBegin();
// find the a reference pressure. The first degree of freedom
// might correspond to non-interior entities which would lead
// to an undefined value, so we have to iterate...
for (unsigned dofIdx = 0; dofIdx < this->numGridDof(); ++ dofIdx) {
if (this->isLocalDof(dofIdx)) {
referencePressure_ =
this->solution(/*timeIdx=*/0)[dofIdx][/*pvIdx=*/Indices::pressureWIdx];
break;
}
}
}
/*!
* \copydoc FvBaseDiscretization::phaseIsConsidered
*/
bool phaseIsConsidered(unsigned phaseIdx) const
{ return phaseIdx == liquidPhaseIdx; }
void registerOutputModules_()
{
ParentType::registerOutputModules_();
}
mutable Scalar referencePressure_;
};
} // namespace Opm
#endif

159
opm/models/richards/richardsnewtonmethod.hh Normal file
View File

@ -0,0 +1,159 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::RichardsNewtonMethod
*/
#ifndef EWOMS_RICHARDS_NEWTON_METHOD_HH
#define EWOMS_RICHARDS_NEWTON_METHOD_HH
#include "richardsproperties.hh"
#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
#include <opm/material/common/Unused.hpp>
#include <dune/common/fvector.hh>
namespace Opm {
/*!
* \ingroup RichardsModel
*
* \brief A Richards model specific Newton method.
*/
template <class TypeTag>
class RichardsNewtonMethod : public GET_PROP_TYPE(TypeTag, DiscNewtonMethod)
{
typedef typename GET_PROP_TYPE(TypeTag, DiscNewtonMethod) ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
typedef typename GET_PROP_TYPE(TypeTag, Linearizer) Linearizer;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum { pressureWIdx = Indices::pressureWIdx };
enum { numPhases = FluidSystem::numPhases };
enum { liquidPhaseIdx = GET_PROP_VALUE(TypeTag, LiquidPhaseIndex) };
enum { gasPhaseIdx = GET_PROP_VALUE(TypeTag, GasPhaseIndex) };
typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
public:
RichardsNewtonMethod(Simulator& simulator) : ParentType(simulator)
{}
protected:
friend NewtonMethod<TypeTag>;
friend ParentType;
/*!
* \copydoc FvBaseNewtonMethod::updatePrimaryVariables_
*/
void updatePrimaryVariables_(unsigned globalDofIdx,
PrimaryVariables& nextValue,
const PrimaryVariables& currentValue,
const EqVector& update,
const EqVector& currentResidual OPM_UNUSED)
{
// normal Newton-Raphson update
nextValue = currentValue;
nextValue -= update;
// do not clamp anything after 4 iterations
if (this->numIterations_ > 4)
return;
const auto& problem = this->simulator_.problem();
// calculate the old wetting phase saturation
const MaterialLawParams& matParams =
problem.materialLawParams(globalDofIdx, /*timeIdx=*/0);
Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
// set the temperature
Scalar T = problem.temperature(globalDofIdx, /*timeIdx=*/0);
fs.setTemperature(T);
/////////
// calculate the phase pressures of the previous iteration
/////////
// first, we have to find the minimum capillary pressure
// (i.e. Sw = 0)
fs.setSaturation(liquidPhaseIdx, 1.0);
fs.setSaturation(gasPhaseIdx, 0.0);
PhaseVector pC;
MaterialLaw::capillaryPressures(pC, matParams, fs);
// non-wetting pressure can be larger than the
// reference pressure if the medium is fully
// saturated by the wetting phase
Scalar pWOld = currentValue[pressureWIdx];
Scalar pNOld =
std::max(problem.referencePressure(globalDofIdx, /*timeIdx=*/0),
pWOld + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
/////////
// find the saturations of the previous iteration
/////////
fs.setPressure(liquidPhaseIdx, pWOld);
fs.setPressure(gasPhaseIdx, pNOld);
PhaseVector satOld;
MaterialLaw::saturations(satOld, matParams, fs);
satOld[liquidPhaseIdx] = std::max<Scalar>(0.0, satOld[liquidPhaseIdx]);
/////////
// find the wetting phase pressures which
// corrospond to a 20% increase and a 20% decrease
// of the wetting saturation
/////////
fs.setSaturation(liquidPhaseIdx, satOld[liquidPhaseIdx] - 0.2);
fs.setSaturation(gasPhaseIdx, 1.0 - (satOld[liquidPhaseIdx] - 0.2));
MaterialLaw::capillaryPressures(pC, matParams, fs);
Scalar pwMin = pNOld - (pC[gasPhaseIdx] - pC[liquidPhaseIdx]);
fs.setSaturation(liquidPhaseIdx, satOld[liquidPhaseIdx] + 0.2);
fs.setSaturation(gasPhaseIdx, 1.0 - (satOld[liquidPhaseIdx] + 0.2));
MaterialLaw::capillaryPressures(pC, matParams, fs);
Scalar pwMax = pNOld - (pC[gasPhaseIdx] - pC[liquidPhaseIdx]);
/////////
// clamp the result to the minimum and the maximum
// pressures we just calculated
/////////
Scalar pW = nextValue[pressureWIdx];
pW = std::max(pwMin, std::min(pW, pwMax));
nextValue[pressureWIdx] = pW;
}
};
} // namespace Opm
#endif

193
opm/models/richards/richardsprimaryvariables.hh Normal file
View File

@ -0,0 +1,193 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::RichardsPrimaryVariables
*/
#ifndef EWOMS_RICHARDS_PRIMARY_VARIABLES_HH
#define EWOMS_RICHARDS_PRIMARY_VARIABLES_HH
#include "richardsproperties.hh"
#include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
#include <opm/material/constraintsolvers/ImmiscibleFlash.hpp>
#include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/Unused.hpp>
#include <dune/common/fvector.hh>
namespace Opm {
/*!
* \ingroup RichardsModel
*
* \brief Represents the primary variables used in the Richards model.
*
* This class is basically a Dune::FieldVector which can retrieve its
* contents from an aribitatry fluid state.
*/
template <class TypeTag>
class RichardsPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
{
typedef FvBasePrimaryVariables<TypeTag> ParentType;
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) EnergyModule;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
// primary variable indices
enum { pressureWIdx = Indices::pressureWIdx };
enum { liquidPhaseIdx = GET_PROP_VALUE(TypeTag, LiquidPhaseIndex) };
enum { gasPhaseIdx = GET_PROP_VALUE(TypeTag, GasPhaseIndex) };
enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
typedef Dune::FieldVector<Scalar, numComponents> ComponentVector;
typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
typedef typename Opm::MathToolbox<Evaluation> Toolbox;
typedef Opm::ImmiscibleFlash<Scalar, FluidSystem> ImmiscibleFlash;
public:
RichardsPrimaryVariables() : ParentType()
{ Opm::Valgrind::SetUndefined(*this); }
/*!
* \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(Scalar)
*/
RichardsPrimaryVariables(Scalar value) : ParentType(value)
{}
/*!
* \copydoc ImmisciblePrimaryVariables::ImmisciblePrimaryVariables(const
* ImmisciblePrimaryVariables& )
*/
RichardsPrimaryVariables(const RichardsPrimaryVariables& value) = default;
RichardsPrimaryVariables& operator=(const RichardsPrimaryVariables& value) = default;
/*!
* \brief Set the primary variables with the wetting phase
* pressure, saturation and temperature.
*
* \param T The temperature [K]
* \param pw The pressure of the wetting phase [Pa]
* \param Sw The saturation of the wetting phase []
* \param matParams The capillary pressure law parameters
*/
void assignImmiscibleFromWetting(Scalar T, Scalar pw, Scalar Sw,
const MaterialLawParams& matParams)
{
Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
fs.setTemperature(T);
fs.setSaturation(liquidPhaseIdx, Sw);
fs.setSaturation(gasPhaseIdx, 1 - Sw);
// set phase pressures
PhaseVector pC;
MaterialLaw::capillaryPressures(pC, matParams, fs);
fs.setPressure(liquidPhaseIdx, pw);
fs.setPressure(gasPhaseIdx, pw + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
assignNaive(fs);
}
/*!
* \brief Set the primary variables with the non-wetting phase
* pressure, saturation and temperature.
*
* \param T The temperature [K]
* \param pn The pressure of the non-wetting phase [Pa]
* \param Sn The saturation of the non-wetting phase []
* \param matParams The capillary pressure law parameters
*/
void assignImmiscibleFromNonWetting(Scalar T, Scalar pn, Scalar Sn,
const MaterialLawParams& matParams)
{
Opm::ImmiscibleFluidState<Scalar, FluidSystem> fs;
fs.setTemperature(T);
fs.setSaturation(liquidPhaseIdx, 1 - Sn);
fs.setSaturation(gasPhaseIdx, Sn);
// set phase pressures
PhaseVector pC;
MaterialLaw::capillaryPressures(pC, matParams, fs);
fs.setPressure(gasPhaseIdx, pn);
fs.setPressure(gasPhaseIdx, pn + (pC[liquidPhaseIdx] - pC[gasPhaseIdx]));
assignNaive(fs);
}
/*!
* \copydoc ImmisciblePrimaryVariables::assignMassConservative
*/
template <class FluidState>
void assignMassConservative(const FluidState& fluidState,
const MaterialLawParams& matParams,
bool isInEquilibrium OPM_UNUSED= false)
{
ComponentVector globalMolarities(0.0);
for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
globalMolarities[compIdx] +=
fluidState.molarity(phaseIdx, compIdx) * fluidState.saturation(phaseIdx);
}
}
Opm::ImmiscibleFluidState<Scalar, FluidSystem> fsFlash;
fsFlash.assign(fluidState);
typename FluidSystem::ParameterCache paramCache;
ImmiscibleFlash::template solve<MaterialLaw>(fsFlash, paramCache,
matParams,
globalMolarities);
assignNaive(fsFlash);
}
/*!
* \copydoc ImmisciblePrimaryVariables::assignNaive
*/
template <class FluidState>
void assignNaive(const FluidState& fluidState)
{
// assign the phase temperatures. this is out-sourced to
// the energy module
EnergyModule::setPriVarTemperatures(*this, fluidState);
(*this)[pressureWIdx] = fluidState.pressure(liquidPhaseIdx);
}
};
} // namespace Opm
#endif

64
opm/models/richards/richardsproperties.hh Normal file
View File

@ -0,0 +1,64 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
* \ingroup RichardsModel
*
* \brief Contains the property declarations for the Richards model.
*/
#ifndef EWOMS_RICHARDS_PROPERTIES_HH
#define EWOMS_RICHARDS_PROPERTIES_HH
#include <opm/models/common/multiphasebaseproperties.hh>
// \{
BEGIN_PROPERTIES
//! The fluid system used for the problem
NEW_PROP_TAG(FluidSystem);
//! The fluid used as the wetting phase (by default, we set the fluid
//! system to the immiscible one, which requires this property.)
NEW_PROP_TAG(WettingFluid);
//! The fluid used as the non-wetting phase (by default, we set the
//! fluid system to the immiscible one, which requires this property.)
NEW_PROP_TAG(NonWettingFluid);
//! Index of the fluid which represents the wetting phase
NEW_PROP_TAG(LiquidPhaseIndex);
//! Index of the fluid which represents the non-wetting phase
NEW_PROP_TAG(GasPhaseIndex);
//! Index of the component which constitutes the liquid
NEW_PROP_TAG(LiquidComponentIndex);
//! Index of the component which constitutes the gas
NEW_PROP_TAG(GasComponentIndex);
// \}
END_PROPERTIES
#endif

146
opm/models/richards/richardsratevector.hh Normal file
View File

@ -0,0 +1,146 @@
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*
This file is part of the Open Porous Media project (OPM).
OPM 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.
OPM 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 OPM. If not, see <http://www.gnu.org/licenses/>.
Consult the COPYING file in the top-level source directory of this
module for the precise wording of the license and the list of
copyright holders.
*/
/*!
* \file
*
* \copydoc Opm::RichardsRateVector
*/
#ifndef EWOMS_RICHARDS_RATE_VECTOR_HH
#define EWOMS_RICHARDS_RATE_VECTOR_HH
#include <dune/common/fvector.hh>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/constraintsolvers/NcpFlash.hpp>
#include "richardsintensivequantities.hh"
namespace Opm {
/*!
* \ingroup RichardsModel
*
* \brief Implements a vector representing mass, molar or volumetric rates.
*
* This class is basically a Dune::FieldVector which can be set using either mass, molar
* or volumetric rates.
*/
template <class TypeTag>
class RichardsRateVector
: public Dune::FieldVector<typename GET_PROP_TYPE(TypeTag, Evaluation),
GET_PROP_VALUE(TypeTag, NumEq)>
{
typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) EnergyModule;
typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
enum { contiEqIdx = Indices::contiEqIdx };
enum { liquidCompIdx = GET_PROP_VALUE(TypeTag, LiquidComponentIndex) };
enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
typedef Dune::FieldVector<Evaluation, numEq> ParentType;
public:
RichardsRateVector() : ParentType()
{ Opm::Valgrind::SetUndefined(*this); }
/*!
* \copydoc ImmiscibleRateVector::ImmiscibleRateVector(Scalar)
*/
RichardsRateVector(const Evaluation& value)
: ParentType(value)
{}
/*!
* \copydoc ImmiscibleRateVector::ImmiscibleRateVector(const
* ImmiscibleRateVector& )
*/
RichardsRateVector(const RichardsRateVector& value)
: ParentType(value)
{}
/*!
* \copydoc ImmiscibleRateVector::setMassRate
*/
void setMassRate(const ParentType& value)
{ ParentType::operator=(value); }
/*!
* \copydoc ImmiscibleRateVector::setMolarRate
*/
void setMolarRate(const ParentType& value)
{
// convert to mass rates
ParentType::operator[](contiEqIdx) =
value[contiEqIdx]*FluidSystem::molarMass(liquidCompIdx);
}
/*!
* \copydoc ImmiscibleRateVector::setEnthalpyRate
*/
template <class RhsEval>
void setEnthalpyRate(const RhsEval& rate)
{ EnergyModule::setEnthalpyRate(*this, rate); }
/*!
* \copydoc ImmiscibleRateVector::setVolumetricRate
*/
template <class FluidState, class RhsEval>
void setVolumetricRate(const FluidState& fluidState, unsigned phaseIdx, const RhsEval& volume)
{
(*this)[contiEqIdx] =
fluidState.density(phaseIdx)
* fluidState.massFraction(phaseIdx, liquidCompIdx)
* volume;
EnergyModule::setEnthalpyRate(*this, fluidState, phaseIdx, volume);
}
/*!
* \brief Assignment operator from a scalar or a function evaluation
*/
template <class RhsEval>
RichardsRateVector& operator=(const RhsEval& value)
{
for (unsigned i=0; i < this->size(); ++i)
(*this)[i] = value;
return *this;
}
/*!
* \brief Assignment operator from another rate vector
*/
RichardsRateVector& operator=(const RichardsRateVector& other)
{
for (unsigned i=0; i < this->size(); ++i)
(*this)[i] = other[i];
return *this;
}
};
} // namespace Opm
#endif