richardsintensivequantities.hh
Go to the documentation of this file.
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 /*
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 2 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 
19  Consult the COPYING file in the top-level source directory of this
20  module for the precise wording of the license and the list of
21  copyright holders.
22 */
28 #ifndef EWOMS_RICHARDS_INTENSIVE_QUANTITIES_HH
29 #define EWOMS_RICHARDS_INTENSIVE_QUANTITIES_HH
30 
31 #include "richardsproperties.hh"
32 
33 #include <opm/material/fluidstates/ImmiscibleFluidState.hpp>
34 
35 #include <dune/common/fvector.hh>
36 #include <dune/common/fmatrix.hh>
37 
38 namespace Ewoms {
39 
46 template <class TypeTag>
48  : public GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities)
49  , public GET_PROP_TYPE(TypeTag, FluxModule)::FluxIntensiveQuantities
50 {
51  typedef typename GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) ParentType;
52  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
53  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
54  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
55  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
56  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
57  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
58  typedef typename GET_PROP_TYPE(TypeTag, FluxModule) FluxModule;
59 
60  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
61  enum { pressureWIdx = Indices::pressureWIdx };
62  enum { numPhases = FluidSystem::numPhases };
63  enum { liquidPhaseIdx = GET_PROP_VALUE(TypeTag, LiquidPhaseIndex) };
64  enum { gasPhaseIdx = GET_PROP_VALUE(TypeTag, GasPhaseIndex) };
65  enum { dimWorld = GridView::dimensionworld };
66 
67  typedef typename FluxModule::FluxIntensiveQuantities FluxIntensiveQuantities;
68  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
69  typedef Dune::FieldVector<Scalar, numPhases> ScalarPhaseVector;
70  typedef Dune::FieldVector<Evaluation, numPhases> PhaseVector;
71  typedef Opm::MathToolbox<Evaluation> Toolbox;
72 
73 public:
75  typedef Opm::ImmiscibleFluidState<Evaluation, FluidSystem> FluidState;
76 
78  {}
79 
81 
82  RichardsIntensiveQuantities& operator=(const RichardsIntensiveQuantities& other) = default;
83 
87  void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
88  {
89  ParentType::update(elemCtx, dofIdx, timeIdx);
90 
91  const auto& T = elemCtx.problem().temperature(elemCtx, dofIdx, timeIdx);
92  fluidState_.setTemperature(T);
93 
94  // material law parameters
95  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
96  const auto& problem = elemCtx.problem();
97  const typename MaterialLaw::Params& materialParams =
98  problem.materialLawParams(elemCtx, dofIdx, timeIdx);
99  const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
100 
102  // calculate the pressures
104 
105  // first, we have to find the minimum capillary pressure (i.e. Sw = 0)
106  fluidState_.setSaturation(liquidPhaseIdx, 1.0);
107  fluidState_.setSaturation(gasPhaseIdx, 0.0);
108  ScalarPhaseVector pC;
109  MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
110 
111  // non-wetting pressure can be larger than the
112  // reference pressure if the medium is fully
113  // saturated by the wetting phase
114  const Evaluation& pW = priVars.makeEvaluation(pressureWIdx, timeIdx);
115  Evaluation pN =
116  Toolbox::max(elemCtx.problem().referencePressure(elemCtx, dofIdx, /*timeIdx=*/0),
117  pW + (pC[gasPhaseIdx] - pC[liquidPhaseIdx]));
118 
120  // calculate the saturations
122  fluidState_.setPressure(liquidPhaseIdx, pW);
123  fluidState_.setPressure(gasPhaseIdx, pN);
124 
125  PhaseVector sat;
126  MaterialLaw::saturations(sat, materialParams, fluidState_);
127  fluidState_.setSaturation(liquidPhaseIdx, sat[liquidPhaseIdx]);
128  fluidState_.setSaturation(gasPhaseIdx, sat[gasPhaseIdx]);
129 
130  typename FluidSystem::template ParameterCache<Evaluation> paramCache;
131  paramCache.updateAll(fluidState_);
132 
133  // compute and set the wetting phase viscosity
134  const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, liquidPhaseIdx);
135  fluidState_.setViscosity(liquidPhaseIdx, mu);
136  fluidState_.setViscosity(gasPhaseIdx, 1e-20);
137 
138  // compute and set the wetting phase density
139  const Evaluation& rho = FluidSystem::density(fluidState_, paramCache, liquidPhaseIdx);
140  fluidState_.setDensity(liquidPhaseIdx, rho);
141  fluidState_.setDensity(gasPhaseIdx, 1e-20);
142 
143  // relperms
144  MaterialLaw::relativePermeabilities(relativePermeability_, materialParams, fluidState_);
145 
146  // mobilities
147  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
148  mobility_[phaseIdx] = relativePermeability_[phaseIdx]/fluidState_.viscosity(phaseIdx);
149 
150  // porosity
151  porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
152 
153  // intrinsic permeability
154  intrinsicPerm_ = problem.intrinsicPermeability(elemCtx, dofIdx, timeIdx);
155 
156  // update the quantities specific for the velocity model
157  FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
158  }
159 
163  const FluidState& fluidState() const
164  { return fluidState_; }
165 
169  const Evaluation& porosity() const
170  { return porosity_; }
171 
175  const DimMatrix& intrinsicPermeability() const
176  { return intrinsicPerm_; }
177 
181  const Evaluation& relativePermeability(unsigned phaseIdx) const
182  { return relativePermeability_[phaseIdx]; }
183 
187  const Evaluation& mobility(unsigned phaseIdx) const
188  { return mobility_[phaseIdx]; }
189 
190 private:
191  FluidState fluidState_;
192  DimMatrix intrinsicPerm_;
193  Evaluation relativePermeability_[numPhases];
194  Evaluation mobility_[numPhases];
195  Evaluation porosity_;
196 };
197 
198 } // namespace Ewoms
199 
200 #endif
Definition: baseauxiliarymodule.hh:37
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
const Evaluation & relativePermeability(unsigned phaseIdx) const
Returns the relative permeability of a given phase within the control volume.
Definition: richardsintensivequantities.hh:181
const DimMatrix & intrinsicPermeability() const
Returns the intrinsic permeability tensor a degree of freedom.
Definition: richardsintensivequantities.hh:175
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
const Evaluation & porosity() const
Returns the average porosity within the control volume.
Definition: richardsintensivequantities.hh:169
const Evaluation & mobility(unsigned phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: richardsintensivequantities.hh:187
Intensive quantities required by the Richards model.
Definition: richardsintensivequantities.hh:47
Contains the property declarations for the Richards model.
Opm::ImmiscibleFluidState< Evaluation, FluidSystem > FluidState
The type returned by the fluidState() method.
Definition: richardsintensivequantities.hh:75
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: richardsintensivequantities.hh:163
void update(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Definition: richardsintensivequantities.hh:87