All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Groups Pages
blackoilintensivequantities.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_BLACK_OIL_INTENSIVE_QUANTITIES_HH
29 #define EWOMS_BLACK_OIL_INTENSIVE_QUANTITIES_HH
30 
31 #include "blackoilproperties.hh"
32 #include "blackoilfluidstate.hh"
35 
36 
37 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
38 #include <opm/common/Valgrind.hpp>
39 
40 #include <dune/common/fmatrix.hh>
41 
42 #include <cstring>
43 #include <utility>
44 
45 namespace Ewoms {
53 template <class TypeTag>
55  : public GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities)
56  , public GET_PROP_TYPE(TypeTag, FluxModule)::FluxIntensiveQuantities
57  , public BlackOilSolventIntensiveQuantities<TypeTag>
58  , public BlackOilPolymerIntensiveQuantities<TypeTag>
59 {
60  typedef typename GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) ParentType;
61  typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) Implementation;
62 
63  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
64  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
65  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
66  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
67  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
68  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
69  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
70  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
71  typedef typename GET_PROP_TYPE(TypeTag, FluxModule) FluxModule;
72 
73  enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
74  enum { enableSolvent = GET_PROP_VALUE(TypeTag, EnableSolvent) };
75  enum { enablePolymer = GET_PROP_VALUE(TypeTag, EnablePolymer) };
76  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
77  enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
78  enum { waterCompIdx = FluidSystem::waterCompIdx };
79  enum { oilCompIdx = FluidSystem::oilCompIdx };
80  enum { gasCompIdx = FluidSystem::gasCompIdx };
81  enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
82  enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
83  enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
84  enum { dimWorld = GridView::dimensionworld };
85  enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
86 
87  static const bool compositionSwitchEnabled = Indices::gasEnabled;
88  static const bool waterEnabled = Indices::waterEnabled;
89 
90 
91  typedef Opm::MathToolbox<Evaluation> Toolbox;
92  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
93  typedef typename FluxModule::FluxIntensiveQuantities FluxIntensiveQuantities;
95 
96 public:
98  {
99  fluidState_.setRs(0.0);
100  fluidState_.setRv(0.0);
101  }
102 
104 
105  BlackOilIntensiveQuantities& operator=(const BlackOilIntensiveQuantities& other) = default;
106 
110  void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
111  {
112  ParentType::update(elemCtx, dofIdx, timeIdx);
113 
114  //fluidState_.setTemperature(elemCtx.problem().temperature(elemCtx, dofIdx, timeIdx));
115 
116  const auto& problem = elemCtx.problem();
117  const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
118 
119  unsigned globalSpaceIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
120  unsigned pvtRegionIdx = priVars.pvtRegionIndex();
121  fluidState_.setPvtRegionIndex(pvtRegionIdx);
122 
123  // extract the water and the gas saturations for convenience
124  Evaluation Sw = 0.0;
125  if (waterEnabled)
126  Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
127 
128  Evaluation Sg = 0.0;
129  if (compositionSwitchEnabled)
130  {
131  if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg)
132  // -> threephase case
133  Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
134  else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
135  // -> gas-water case
136  Sg = 1.0 - Sw;
137 
138  // deal with solvent
139  if (enableSolvent)
140  Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
141  }
142  else
143  {
144  assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs);
145  // -> oil-water case
146  Sg = 0.0;
147  }
148  }
149 
150  Opm::Valgrind::CheckDefined(Sg);
151  Opm::Valgrind::CheckDefined(Sw);
152 
153  Evaluation So = 1.0 - Sw - Sg;
154 
155  // deal with solvent
156  if (enableSolvent)
157  So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
158 
159  fluidState_.setSaturation(waterPhaseIdx, Sw);
160  fluidState_.setSaturation(gasPhaseIdx, Sg);
161  fluidState_.setSaturation(oilPhaseIdx, So);
162 
163  asImp_().solventPreSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
164 
165  // now we compute all phase pressures
166  Evaluation pC[numPhases];
167  const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
168  MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
169 
170  //oil is the reference phase for pressure
171  if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
172  const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
173  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
174  fluidState_.setPressure(phaseIdx, pg + (pC[phaseIdx] - pC[gasPhaseIdx]));
175  }
176 
177  else {
178  const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
179  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
180  fluidState_.setPressure(phaseIdx, po + (pC[phaseIdx] - pC[oilPhaseIdx]));
181  }
182 
183  // calculate relative permeabilities. note that we store the result into the
184  // mobility_ class attribute. the division by the phase viscosity happens later.
185  MaterialLaw::relativePermeabilities(mobility_, materialParams, fluidState_);
186  Opm::Valgrind::CheckDefined(mobility_);
187 
188  // update the Saturation functions for the blackoil solvent module.
189  asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
190 
191  Scalar SoMax = elemCtx.model().maxOilSaturation(globalSpaceIdx);
192 
193  // take the meaning of the switiching primary variable into account for the gas
194  // and oil phase compositions
195  if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
196  // in the threephase case, gas and oil phases are potentially present, i.e.,
197  // we use the compositions of the gas-saturated oil and oil-saturated gas.
198  if (FluidSystem::enableDissolvedGas()) {
199  const Evaluation& RsSat =
200  FluidSystem::saturatedDissolutionFactor(fluidState_,
201  oilPhaseIdx,
202  pvtRegionIdx,
203  SoMax);
204  fluidState_.setRs(RsSat);
205  }
206  else
207  fluidState_.setRs(0.0);
208 
209  if (FluidSystem::enableVaporizedOil()) {
210  const Evaluation& RvSat =
211  FluidSystem::saturatedDissolutionFactor(fluidState_,
212  gasPhaseIdx,
213  pvtRegionIdx,
214  SoMax);
215  fluidState_.setRv(RvSat);
216  }
217  else
218  fluidState_.setRv(0.0);
219  }
220  else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
221  // if the switching variable is the mole fraction of the gas component in the
222  // oil phase, we can directly set the composition of the oil phase
223  const auto& Rs = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
224  fluidState_.setRs(Rs);
225 
226  if (FluidSystem::enableVaporizedOil()) {
227  // the gas phase is not present, but we need to compute its "composition"
228  // for the gravity correction anyway
229  const auto& RvSat =
230  FluidSystem::saturatedDissolutionFactor(fluidState_,
231  gasPhaseIdx,
232  pvtRegionIdx,
233  SoMax);
234 
235  fluidState_.setRv(RvSat);
236  }
237  else
238  fluidState_.setRv(0.0);
239  }
240  else {
241  assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv);
242 
243  const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
244  fluidState_.setRv(Rv);
245 
246  if (FluidSystem::enableDissolvedGas()) {
247  // the oil phase is not present, but we need to compute its "composition" for
248  // the gravity correction anyway
249  const auto& RsSat =
250  FluidSystem::saturatedDissolutionFactor(fluidState_,
251  oilPhaseIdx,
252  pvtRegionIdx,
253  SoMax);
254 
255  fluidState_.setRs(RsSat);
256  }
257  else
258  fluidState_.setRs(0.0);
259  }
260 
261  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
262  typename FluidSystem::template ParameterCache<Evaluation> paramCache;
263  paramCache.setRegionIndex(pvtRegionIdx);
264  paramCache.setMaxOilSat(SoMax);
265  paramCache.updateAll(fluidState_);
266 
267  // compute the phase densities and transform the phase permeabilities into mobilities
268  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
269  if (!FluidSystem::phaseIsActive(phaseIdx))
270  continue;
271 
272  const auto& b = FluidSystem::inverseFormationVolumeFactor(fluidState_, phaseIdx, pvtRegionIdx);
273  fluidState_.setInvB(phaseIdx, b);
274 
275  const auto& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
276  mobility_[phaseIdx] /= mu;
277  }
278  Opm::Valgrind::CheckDefined(mobility_);
279 
280  // calculate the phase densities
281  Evaluation rho;
282  if (FluidSystem::phaseIsActive(waterPhaseIdx)) {
283  rho = fluidState_.invB(waterPhaseIdx);
284  rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
285  fluidState_.setDensity(waterPhaseIdx, rho);
286  }
287 
288  if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
289  rho = fluidState_.invB(gasPhaseIdx);
290  rho *= FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
291  if (FluidSystem::enableVaporizedOil()) {
292  rho +=
293  fluidState_.invB(gasPhaseIdx) *
294  fluidState_.Rv() *
295  FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
296  }
297  fluidState_.setDensity(gasPhaseIdx, rho);
298  }
299 
300  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
301  rho = fluidState_.invB(oilPhaseIdx);
302  rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
303  if (FluidSystem::enableDissolvedGas()) {
304  rho +=
305  fluidState_.invB(oilPhaseIdx) *
306  fluidState_.Rs() *
307  FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
308  }
309  fluidState_.setDensity(oilPhaseIdx, rho);
310  }
311 
312  // retrieve the porosity from the problem
313  porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
314 
315  // the porosity must be modified by the compressibility of the
316  // rock...
317  Scalar rockCompressibility = problem.rockCompressibility(elemCtx, dofIdx, timeIdx);
318  if (rockCompressibility > 0.0) {
319  Scalar rockRefPressure = problem.rockReferencePressure(elemCtx, dofIdx, timeIdx);
320  Evaluation x = rockCompressibility*(fluidState_.pressure(oilPhaseIdx) - rockRefPressure);
321  porosity_ *= 1.0 + x + 0.5*x*x;
322  }
323 
324  asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
325  asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
326 
327 
328  // update the quantities which are required by the chosen
329  // velocity model
330  FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
331 
332 #ifndef NDEBUG
333  // some safety checks in debug mode
334  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++ phaseIdx) {
335  if (!FluidSystem::phaseIsActive(phaseIdx))
336  continue;
337 
338  assert(std::isfinite(Toolbox::value(fluidState_.density(phaseIdx))));
339  assert(std::isfinite(Toolbox::value(fluidState_.saturation(phaseIdx))));
340  assert(std::isfinite(Toolbox::value(fluidState_.temperature(phaseIdx))));
341  assert(std::isfinite(Toolbox::value(fluidState_.pressure(phaseIdx))));
342  assert(std::isfinite(Toolbox::value(fluidState_.invB(phaseIdx))));
343  }
344  assert(std::isfinite(Toolbox::value(fluidState_.Rs())));
345  assert(std::isfinite(Toolbox::value(fluidState_.Rv())));
346 #endif
347  }
348 
352  const FluidState& fluidState() const
353  { return fluidState_; }
354 
358  const Evaluation& mobility(unsigned phaseIdx) const
359  { return mobility_[phaseIdx]; }
360 
364  const Evaluation& porosity() const
365  { return porosity_; }
366 
380  auto pvtRegionIndex() const
381  -> decltype(std::declval<FluidState>().pvtRegionIndex())
382  { return fluidState_.pvtRegionIndex(); }
383 
387  Evaluation relativePermeability(unsigned phaseIdx) const
388  {
389  // warning: slow
390  return fluidState_.viscosity(phaseIdx)*mobility(phaseIdx);
391  }
392 
393 private:
396 
397 
398  Implementation& asImp_()
399  { return *static_cast<Implementation*>(this); }
400 
401  FluidState fluidState_;
402  Evaluation porosity_;
403  Evaluation mobility_[numPhases];
404 };
405 
406 } // namespace Ewoms
407 
408 #endif
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition: blackoilpolymermodules.hh:839
const Evaluation & mobility(unsigned phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: blackoilintensivequantities.hh:358
Provides the volumetric quantities required for the equations needed by the solvents extension of the...
Definition: blackoilsolventmodules.hh:879
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
Contains the classes required to extend the black-oil model by solvents.
Implements a &quot;taylor-made&quot; fluid state class for the black-oil model.
#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: blackoilintensivequantities.hh:364
Declares the properties required by the black oil model.
void update(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Definition: blackoilintensivequantities.hh:110
auto pvtRegionIndex() const -> decltype(std::declval< FluidState >().pvtRegionIndex())
Returns the index of the PVT region used to calculate the thermodynamic quantities.
Definition: blackoilintensivequantities.hh:380
Evaluation relativePermeability(unsigned phaseIdx) const
Returns the relative permeability of a given phase within the control volume.
Definition: blackoilintensivequantities.hh:387
Implements a &quot;taylor-made&quot; fluid state class for the black-oil model.
Definition: blackoilfluidstate.hh:49
Contains the classes required to extend the black-oil model by polymer.
Contains the quantities which are are constant within a finite volume in the black-oil model...
Definition: blackoilintensivequantities.hh:54
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: blackoilintensivequantities.hh:352