pvsintensivequantities.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_PVS_INTENSIVE_QUANTITIES_HH
29 #define EWOMS_PVS_INTENSIVE_QUANTITIES_HH
30 
31 #include "pvsproperties.hh"
32 
35 
36 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
37 #include <opm/material/constraintsolvers/MiscibleMultiPhaseComposition.hpp>
38 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
39 #include <opm/common/Valgrind.hpp>
40 
41 #include <dune/common/fvector.hh>
42 #include <dune/common/fmatrix.hh>
43 
44 #include <iostream>
45 
46 namespace Ewoms {
55 template <class TypeTag>
57  : public GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities)
58  , public DiffusionIntensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableDiffusion) >
59  , public EnergyIntensiveQuantities<TypeTag, GET_PROP_VALUE(TypeTag, EnableEnergy) >
60  , public GET_PROP_TYPE(TypeTag, FluxModule)::FluxIntensiveQuantities
61 {
62  typedef typename GET_PROP_TYPE(TypeTag, DiscIntensiveQuantities) ParentType;
63 
64  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
65  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
66  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
67  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
68  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
69  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
70  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
71  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
72  typedef typename GET_PROP_TYPE(TypeTag, FluxModule) FluxModule;
73 
74  enum { switch0Idx = Indices::switch0Idx };
75  enum { pressure0Idx = Indices::pressure0Idx };
76  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
77  enum { numComponents = GET_PROP_VALUE(TypeTag, NumComponents) };
78  enum { enableDiffusion = GET_PROP_VALUE(TypeTag, EnableDiffusion) };
79  enum { enableEnergy = GET_PROP_VALUE(TypeTag, EnableEnergy) };
80  enum { dimWorld = GridView::dimensionworld };
81 
82  typedef Opm::MathToolbox<Evaluation> Toolbox;
83  typedef Opm::MiscibleMultiPhaseComposition<Scalar, FluidSystem, Evaluation> MiscibleMultiPhaseComposition;
84  typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem, Evaluation> ComputeFromReferencePhase;
85 
86  typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
87  typedef Dune::FieldVector<Evaluation, numPhases> EvalPhaseVector;
88  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
89 
90  typedef typename FluxModule::FluxIntensiveQuantities FluxIntensiveQuantities;
93 
94 public:
96  typedef Opm::CompositionalFluidState<Evaluation, FluidSystem> FluidState;
97 
99  { }
100 
101  PvsIntensiveQuantities(const PvsIntensiveQuantities& other) = default;
102 
103  PvsIntensiveQuantities& operator=(const PvsIntensiveQuantities& other) = default;
104 
108  void update(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
109  {
110  ParentType::update(elemCtx, dofIdx, timeIdx);
111  EnergyIntensiveQuantities::updateTemperatures_(fluidState_, elemCtx, dofIdx, timeIdx);
112 
113  const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
114  const auto& problem = elemCtx.problem();
115 
117  // set the saturations
119  Evaluation sumSat = 0.0;
120  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
121  fluidState_.setSaturation(phaseIdx, priVars.explicitSaturationValue(phaseIdx, timeIdx));
122  Opm::Valgrind::CheckDefined(fluidState_.saturation(phaseIdx));
123  sumSat += fluidState_.saturation(phaseIdx);
124  }
125  Opm::Valgrind::CheckDefined(priVars.implicitSaturationIdx());
126  Opm::Valgrind::CheckDefined(sumSat);
127  fluidState_.setSaturation(priVars.implicitSaturationIdx(), 1.0 - sumSat);
128 
130  // set the pressures of the fluid phases
132 
133  // calculate capillary pressure
134  const MaterialLawParams& materialParams =
135  problem.materialLawParams(elemCtx, dofIdx, timeIdx);
136  EvalPhaseVector pC;
137  MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
138 
139  // set the absolute phase pressures in the fluid state
140  const Evaluation& p0 = priVars.makeEvaluation(pressure0Idx, timeIdx);
141  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
142  fluidState_.setPressure(phaseIdx, p0 + (pC[phaseIdx] - pC[0]));
143 
145  // calculate the phase compositions
147 
148  typename FluidSystem::template ParameterCache<Evaluation> paramCache;
149  unsigned lowestPresentPhaseIdx = priVars.lowestPresentPhaseIdx();
150  unsigned numNonPresentPhases = 0;
151  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
152  if (!priVars.phaseIsPresent(phaseIdx))
153  ++numNonPresentPhases;
154  }
155 
156  // now comes the tricky part: calculate phase compositions
157  if (numNonPresentPhases == numPhases - 1) {
158  // only one phase is present, i.e. the primary variables
159  // contain the complete composition of the phase
160  Evaluation sumx = 0.0;
161  for (unsigned compIdx = 1; compIdx < numComponents; ++compIdx) {
162  const Evaluation& x = priVars.makeEvaluation(switch0Idx + compIdx - 1, timeIdx);
163  fluidState_.setMoleFraction(lowestPresentPhaseIdx, compIdx, x);
164  sumx += x;
165  }
166 
167  // set the mole fraction of the first component
168  fluidState_.setMoleFraction(lowestPresentPhaseIdx, 0, 1 - sumx);
169 
170  // calculate the composition of the remaining phases (as
171  // well as the densities of all phases). this is the job
172  // of the "ComputeFromReferencePhase" constraint solver
173  ComputeFromReferencePhase::solve(fluidState_, paramCache,
174  lowestPresentPhaseIdx,
175  /*setViscosity=*/true,
176  /*setEnthalpy=*/false);
177  }
178  else {
179  // create the auxiliary constraints
180  unsigned numAuxConstraints = numComponents + numNonPresentPhases - numPhases;
181  Opm::MMPCAuxConstraint<Evaluation> auxConstraints[numComponents];
182 
183  unsigned auxIdx = 0;
184  unsigned switchIdx = 0;
185  for (; switchIdx < numPhases - 1; ++switchIdx) {
186  unsigned compIdx = switchIdx + 1;
187  unsigned switchPhaseIdx = switchIdx;
188  if (switchIdx >= lowestPresentPhaseIdx)
189  switchPhaseIdx += 1;
190 
191  if (!priVars.phaseIsPresent(switchPhaseIdx)) {
192  auxConstraints[auxIdx].set(lowestPresentPhaseIdx, compIdx,
193  priVars.makeEvaluation(switch0Idx + switchIdx, timeIdx));
194  ++auxIdx;
195  }
196  }
197 
198  for (; auxIdx < numAuxConstraints; ++auxIdx, ++switchIdx) {
199  unsigned compIdx = numPhases - numNonPresentPhases + auxIdx;
200  auxConstraints[auxIdx].set(lowestPresentPhaseIdx, compIdx,
201  priVars.makeEvaluation(switch0Idx + switchIdx, timeIdx));
202  }
203 
204  // both phases are present, i.e. phase compositions are a result of the the
205  // gas <-> liquid equilibrium. This is the job of the
206  // "MiscibleMultiPhaseComposition" constraint solver
207  MiscibleMultiPhaseComposition::solve(fluidState_, paramCache,
208  priVars.phasePresence(),
209  auxConstraints,
210  numAuxConstraints,
211  /*setViscosity=*/true,
212  /*setEnthalpy=*/false);
213  }
214 
215 #ifndef NDEBUG
216  // make valgrind happy and set the enthalpies to NaN
217  if (!enableEnergy) {
218  Scalar myNan = std::numeric_limits<Scalar>::quiet_NaN();
219  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
220  fluidState_.setEnthalpy(phaseIdx, myNan);
221  }
222 #endif
223 
225  // calculate the remaining quantities
227 
228  // calculate relative permeabilities
229  MaterialLaw::relativePermeabilities(relativePermeability_,
230  materialParams, fluidState_);
231  Opm::Valgrind::CheckDefined(relativePermeability_);
232 
233  // mobilities
234  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
235  mobility_[phaseIdx] =
236  relativePermeability_[phaseIdx] / fluidState().viscosity(phaseIdx);
237 
238  // porosity
239  porosity_ = problem.porosity(elemCtx, dofIdx, timeIdx);
240  Opm::Valgrind::CheckDefined(porosity_);
241 
242  // intrinsic permeability
243  intrinsicPerm_ = problem.intrinsicPermeability(elemCtx, dofIdx, timeIdx);
244 
245  // update the quantities specific for the velocity model
246  FluxIntensiveQuantities::update_(elemCtx, dofIdx, timeIdx);
247 
248  // energy related quantities
249  EnergyIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
250 
251  // update the diffusion specific quantities of the intensive quantities
252  DiffusionIntensiveQuantities::update_(fluidState_, paramCache, elemCtx, dofIdx, timeIdx);
253 
254  fluidState_.checkDefined();
255  }
256 
260  const FluidState& fluidState() const
261  { return fluidState_; }
262 
266  const DimMatrix& intrinsicPermeability() const
267  { return intrinsicPerm_; }
268 
272  const Evaluation& relativePermeability(unsigned phaseIdx) const
273  { return relativePermeability_[phaseIdx]; }
274 
278  const Evaluation& mobility(unsigned phaseIdx) const
279  { return mobility_[phaseIdx]; }
280 
284  const Evaluation& porosity() const
285  { return porosity_; }
286 
287 private:
288  FluidState fluidState_;
289  Evaluation porosity_;
290  DimMatrix intrinsicPerm_;
291  Evaluation relativePermeability_[numPhases];
292  Evaluation mobility_[numPhases];
293 };
294 
295 } // namespace Ewoms
296 
297 #endif
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes...
Definition: diffusionmodule.hh:146
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
void update(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Definition: pvsintensivequantities.hh:108
Contains the classes required to consider energy as a conservation quantity in a multi-phase module...
Opm::CompositionalFluidState< Evaluation, FluidSystem > FluidState
The type of the object returned by the fluidState() method.
Definition: pvsintensivequantities.hh:96
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
Provides the volumetric quantities required for the energy equation.
Definition: energymodule.hh:539
const Evaluation & relativePermeability(unsigned phaseIdx) const
Returns the relative permeability of a given phase within the control volume.
Definition: pvsintensivequantities.hh:272
const FluidState & fluidState() const
Returns the phase state for the control-volume.
Definition: pvsintensivequantities.hh:260
const Evaluation & mobility(unsigned phaseIdx) const
Returns the effective mobility of a given phase within the control volume.
Definition: pvsintensivequantities.hh:278
const Evaluation & porosity() const
Returns the average porosity within the control volume.
Definition: pvsintensivequantities.hh:284
const DimMatrix & intrinsicPermeability() const
Returns the intrinsic permeability tensor a degree of freedom.
Definition: pvsintensivequantities.hh:266
Contains the quantities which are are constant within a finite volume in the compositional multi-phas...
Definition: pvsintensivequantities.hh:56
Declares the properties required for the compositional multi-phase primary variable switching model...
Classes required for molecular diffusion.