My Project
BlackOilFluidSystem.hpp
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 */
27 #ifndef OPM_BLACK_OIL_FLUID_SYSTEM_HPP
28 #define OPM_BLACK_OIL_FLUID_SYSTEM_HPP
29 
35 
38 
43 
44 #if HAVE_ECL_INPUT
45 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
46 #include <opm/parser/eclipse/EclipseState/Tables/FlatTable.hpp>
47 #include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
48 #endif
49 
50 #include <memory>
51 #include <vector>
52 #include <array>
53 
54 namespace Opm {
55 namespace BlackOil {
56 OPM_GENERATE_HAS_MEMBER(Rs, ) // Creates 'HasMember_Rs<T>'.
57 OPM_GENERATE_HAS_MEMBER(Rv, ) // Creates 'HasMember_Rv<T>'.
58 OPM_GENERATE_HAS_MEMBER(saltConcentration, )
59 
60 template <class FluidSystem, class FluidState, class LhsEval>
61 LhsEval getRs_(typename std::enable_if<!HasMember_Rs<FluidState>::value, const FluidState&>::type fluidState,
62  unsigned regionIdx)
63 {
64  const auto& XoG =
65  decay<LhsEval>(fluidState.massFraction(FluidSystem::oilPhaseIdx, FluidSystem::gasCompIdx));
66  return FluidSystem::convertXoGToRs(XoG, regionIdx);
67 }
68 
69 template <class FluidSystem, class FluidState, class LhsEval>
70 auto getRs_(typename std::enable_if<HasMember_Rs<FluidState>::value, const FluidState&>::type fluidState,
71  unsigned)
72  -> decltype(decay<LhsEval>(fluidState.Rs()))
73 { return decay<LhsEval>(fluidState.Rs()); }
74 
75 template <class FluidSystem, class FluidState, class LhsEval>
76 LhsEval getRv_(typename std::enable_if<!HasMember_Rv<FluidState>::value, const FluidState&>::type fluidState,
77  unsigned regionIdx)
78 {
79  const auto& XgO =
80  decay<LhsEval>(fluidState.massFraction(FluidSystem::gasPhaseIdx, FluidSystem::oilCompIdx));
81  return FluidSystem::convertXgOToRv(XgO, regionIdx);
82 }
83 
84 template <class FluidSystem, class FluidState, class LhsEval>
85 auto getRv_(typename std::enable_if<HasMember_Rv<FluidState>::value, const FluidState&>::type fluidState,
86  unsigned)
87  -> decltype(decay<LhsEval>(fluidState.Rv()))
88 { return decay<LhsEval>(fluidState.Rv()); }
89 
90 template <class FluidSystem, class FluidState, class LhsEval>
91 LhsEval getSaltConcentration_(typename std::enable_if<!HasMember_saltConcentration<FluidState>::value,
92  const FluidState&>::type,
93  unsigned)
94 {return 0.0;}
95 
96 template <class FluidSystem, class FluidState, class LhsEval>
97 auto getSaltConcentration_(typename std::enable_if<HasMember_saltConcentration<FluidState>::value, const FluidState&>::type fluidState,
98  unsigned)
99  -> decltype(decay<LhsEval>(fluidState.saltConcentration()))
100 { return decay<LhsEval>(fluidState.saltConcentration()); }
101 
102 }
103 
110 template <class Scalar, class IndexTraits = BlackOilDefaultIndexTraits>
111 class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<Scalar, IndexTraits> >
112 {
114 
115 public:
119 
121  template <class EvaluationT>
122  struct ParameterCache : public NullParameterCache<EvaluationT>
123  {
124  typedef EvaluationT Evaluation;
125 
126  public:
127  ParameterCache(Scalar maxOilSat = 1.0, unsigned regionIdx=0)
128  {
129  maxOilSat_ = maxOilSat;
130  regionIdx_ = regionIdx;
131  }
132 
140  template <class OtherCache>
141  void assignPersistentData(const OtherCache& other)
142  {
143  regionIdx_ = other.regionIndex();
144  maxOilSat_ = other.maxOilSat();
145  }
146 
154  unsigned regionIndex() const
155  { return regionIdx_; }
156 
164  void setRegionIndex(unsigned val)
165  { regionIdx_ = val; }
166 
167  const Evaluation& maxOilSat() const
168  { return maxOilSat_; }
169 
170  void setMaxOilSat(const Evaluation& val)
171  { maxOilSat_ = val; }
172 
173  private:
174  Evaluation maxOilSat_;
175  unsigned regionIdx_;
176  };
177 
178  /****************************************
179  * Initialization
180  ****************************************/
181 #if HAVE_ECL_INPUT
185  static void initFromState(const EclipseState& eclState, const Schedule& schedule)
186  {
187  size_t numRegions = eclState.runspec().tabdims().getNumPVTTables();
189 
190  numActivePhases_ = 0;
191  std::fill_n(&phaseIsActive_[0], numPhases, false);
192 
193 
194  if (eclState.runspec().phases().active(Phase::OIL)) {
195  phaseIsActive_[oilPhaseIdx] = true;
196  ++ numActivePhases_;
197  }
198 
199  if (eclState.runspec().phases().active(Phase::GAS)) {
200  phaseIsActive_[gasPhaseIdx] = true;
201  ++ numActivePhases_;
202  }
203 
204  if (eclState.runspec().phases().active(Phase::WATER)) {
205  phaseIsActive_[waterPhaseIdx] = true;
206  ++ numActivePhases_;
207  }
208 
209  // set the surface conditions using the STCOND keyword
210  surfaceTemperature = eclState.getTableManager().stCond().temperature;
211  surfacePressure = eclState.getTableManager().stCond().pressure;
212 
213  // The reservoir temperature does not really belong into the table manager. TODO:
214  // change this in opm-parser
215  setReservoirTemperature(eclState.getTableManager().rtemp());
216 
217  // this fluidsystem only supports two or three phases
218  assert(numActivePhases_ >= 1 && numActivePhases_ <= 3);
219 
220  setEnableDissolvedGas(eclState.getSimulationConfig().hasDISGAS());
221  setEnableVaporizedOil(eclState.getSimulationConfig().hasVAPOIL());
222 
223  if (phaseIsActive(gasPhaseIdx)) {
224  gasPvt_ = std::make_shared<GasPvt>();
225  gasPvt_->initFromState(eclState, schedule);
226  }
227 
228  if (phaseIsActive(oilPhaseIdx)) {
229  oilPvt_ = std::make_shared<OilPvt>();
230  oilPvt_->initFromState(eclState, schedule);
231  }
232 
234  waterPvt_ = std::make_shared<WaterPvt>();
235  waterPvt_->initFromState(eclState, schedule);
236  }
237 
238  // set the reference densities of all PVT regions
239  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
240  setReferenceDensities(phaseIsActive(oilPhaseIdx)? oilPvt_->oilReferenceDensity(regionIdx):700.,
241  phaseIsActive(waterPhaseIdx)? waterPvt_->waterReferenceDensity(regionIdx):1000.,
242  phaseIsActive(gasPhaseIdx)? gasPvt_->gasReferenceDensity(regionIdx):2.,
243  regionIdx);
244  }
245 
246  // set default molarMass and mappings
247  initEnd();
248 
249  // use molarMass of CO2 and Brine as default
250  // when we are using the the CO2STORE option
251  // NB the oil component is used internally for
252  // brine
253  if (eclState.runspec().co2Storage()) {
254  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
255  molarMass_[regionIdx][oilCompIdx] = BrineCo2Pvt<Scalar>::Brine::molarMass();
256  molarMass_[regionIdx][gasCompIdx] = BrineCo2Pvt<Scalar>::CO2::molarMass();
257  }
258  }
259 
260  setEnableDiffusion(eclState.getSimulationConfig().isDiffusive());
261  if(enableDiffusion()) {
262  const auto& diffCoeffTables = eclState.getTableManager().getDiffusionCoefficientTable();
263  if(!diffCoeffTables.empty()) {
264  // if diffusion coefficient table is empty we relay on the PVT model to
265  // to give us the coefficients.
266  diffusionCoefficients_.resize(numRegions,{0,0,0,0,0,0,0,0,0});
267  assert(diffCoeffTables.size() == numRegions);
268  for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
269  const auto& diffCoeffTable = diffCoeffTables[regionIdx];
270  molarMass_[regionIdx][oilCompIdx] = diffCoeffTable.oil_mw;
271  molarMass_[regionIdx][gasCompIdx] = diffCoeffTable.gas_mw;
272  setDiffusionCoefficient(diffCoeffTable.gas_in_gas, gasCompIdx, gasPhaseIdx, regionIdx);
273  setDiffusionCoefficient(diffCoeffTable.oil_in_gas, oilCompIdx, gasPhaseIdx, regionIdx);
274  setDiffusionCoefficient(diffCoeffTable.gas_in_oil, gasCompIdx, oilPhaseIdx, regionIdx);
275  setDiffusionCoefficient(diffCoeffTable.oil_in_oil, oilCompIdx, oilPhaseIdx, regionIdx);
276  if(diffCoeffTable.gas_in_oil_cross_phase > 0 || diffCoeffTable.oil_in_oil_cross_phase > 0) {
277  throw std::runtime_error("Cross phase diffusion is set in the deck, but not implemented in Flow. "
278  "Please default DIFFC item 7 and item 8 or set it to zero.");
279  }
280  }
281  }
282  }
283  }
284 #endif // HAVE_ECL_INPUT
285 
294  static void initBegin(size_t numPvtRegions)
295  {
296  isInitialized_ = false;
297 
298  enableDissolvedGas_ = true;
299  enableVaporizedOil_ = false;
300  enableDiffusion_ = false;
301 
302  oilPvt_ = nullptr;
303  gasPvt_ = nullptr;
304  waterPvt_ = nullptr;
305 
306  surfaceTemperature = 273.15 + 15.56; // [K]
307  surfacePressure = 1.01325e5; // [Pa]
309 
310  numActivePhases_ = numPhases;
311  std::fill_n(&phaseIsActive_[0], numPhases, true);
312 
313  resizeArrays_(numPvtRegions);
314  }
315 
322  static void setEnableDissolvedGas(bool yesno)
323  { enableDissolvedGas_ = yesno; }
324 
331  static void setEnableVaporizedOil(bool yesno)
332  { enableVaporizedOil_ = yesno; }
333 
339  static void setEnableDiffusion(bool yesno)
340  { enableDiffusion_ = yesno; }
341 
342 
346  static void setGasPvt(std::shared_ptr<GasPvt> pvtObj)
347  { gasPvt_ = pvtObj; }
348 
352  static void setOilPvt(std::shared_ptr<OilPvt> pvtObj)
353  { oilPvt_ = pvtObj; }
354 
358  static void setWaterPvt(std::shared_ptr<WaterPvt> pvtObj)
359  { waterPvt_ = pvtObj; }
360 
368  static void setReferenceDensities(Scalar rhoOil,
369  Scalar rhoWater,
370  Scalar rhoGas,
371  unsigned regionIdx)
372  {
373  referenceDensity_[regionIdx][oilPhaseIdx] = rhoOil;
374  referenceDensity_[regionIdx][waterPhaseIdx] = rhoWater;
375  referenceDensity_[regionIdx][gasPhaseIdx] = rhoGas;
376  }
377 
378 
382  static void initEnd()
383  {
384  // calculate the final 2D functions which are used for interpolation.
385  size_t numRegions = molarMass_.size();
386  for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
387  // calculate molar masses
388 
389  // water is simple: 18 g/mol
390  molarMass_[regionIdx][waterCompIdx] = 18e-3;
391 
392  if (phaseIsActive(gasPhaseIdx)) {
393  // for gas, we take the density at standard conditions and assume it to be ideal
396  Scalar rho_g = referenceDensity_[/*regionIdx=*/0][gasPhaseIdx];
397  molarMass_[regionIdx][gasCompIdx] = Constants<Scalar>::R*T*rho_g / p;
398  }
399  else
400  // hydrogen gas. we just set this do avoid NaNs later
401  molarMass_[regionIdx][gasCompIdx] = 2e-3;
402 
403  // finally, for oil phase, we take the molar mass from the spe9 paper
404  molarMass_[regionIdx][oilCompIdx] = 175e-3; // kg/mol
405  }
406 
407 
408  int activePhaseIdx = 0;
409  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
410  if(phaseIsActive(phaseIdx)){
411  canonicalToActivePhaseIdx_[phaseIdx] = activePhaseIdx;
412  activeToCanonicalPhaseIdx_[activePhaseIdx] = phaseIdx;
413  activePhaseIdx++;
414  }
415  }
416  isInitialized_ = true;
417  }
418 
419  static bool isInitialized()
420  { return isInitialized_; }
421 
422  /****************************************
423  * Generic phase properties
424  ****************************************/
425 
427  static const unsigned numPhases = 3;
428 
430  static const unsigned waterPhaseIdx = IndexTraits::waterPhaseIdx;
432  static const unsigned oilPhaseIdx = IndexTraits::oilPhaseIdx;
434  static const unsigned gasPhaseIdx = IndexTraits::gasPhaseIdx;
435 
438 
441 
443  static const char* phaseName(unsigned phaseIdx)
444  {
445  switch (phaseIdx) {
446  case waterPhaseIdx:
447  return "water";
448  case oilPhaseIdx:
449  return "oil";
450  case gasPhaseIdx:
451  return "gas";
452 
453  default:
454  throw std::logic_error("Phase index " + std::to_string(phaseIdx) + " is unknown");
455  }
456  }
457 
459  static bool isLiquid(unsigned phaseIdx)
460  {
461  assert(phaseIdx < numPhases);
462  return phaseIdx != gasPhaseIdx;
463  }
464 
465  /****************************************
466  * Generic component related properties
467  ****************************************/
468 
470  static const unsigned numComponents = 3;
471 
473  static const unsigned oilCompIdx = IndexTraits::oilCompIdx;
475  static const unsigned waterCompIdx = IndexTraits::waterCompIdx;
477  static const unsigned gasCompIdx = IndexTraits::gasCompIdx;
478 
479 protected:
480  static unsigned char numActivePhases_;
481  static std::array<bool,numPhases> phaseIsActive_;
482 
483 public:
485  static unsigned numActivePhases()
486  { return numActivePhases_; }
487 
489  static unsigned phaseIsActive(unsigned phaseIdx)
490  {
491  assert(phaseIdx < numPhases);
492  return phaseIsActive_[phaseIdx];
493  }
494 
496  static constexpr unsigned solventComponentIndex(unsigned phaseIdx)
497  {
498  switch (phaseIdx) {
499  case waterPhaseIdx:
500  return waterCompIdx;
501  case oilPhaseIdx:
502  return oilCompIdx;
503  case gasPhaseIdx:
504  return gasCompIdx;
505 
506  default:
507  throw std::logic_error("Phase index " + std::to_string(phaseIdx) + " is unknown");
508  }
509  }
510 
512  static constexpr unsigned soluteComponentIndex(unsigned phaseIdx)
513  {
514  switch (phaseIdx) {
515  case waterPhaseIdx:
516  throw std::logic_error("The water phase does not have any solutes in the black oil model!");
517  case oilPhaseIdx:
518  return gasCompIdx;
519  case gasPhaseIdx:
520  return oilCompIdx;
521 
522  default:
523  throw std::logic_error("Phase index " + std::to_string(phaseIdx) + " is unknown");
524  }
525  }
526 
528  static const char* componentName(unsigned compIdx)
529  {
530  switch (compIdx) {
531  case waterCompIdx:
532  return "Water";
533  case oilCompIdx:
534  return "Oil";
535  case gasCompIdx:
536  return "Gas";
537 
538  default:
539  throw std::logic_error("Component index " + std::to_string(compIdx) + " is unknown");
540  }
541  }
542 
544  static Scalar molarMass(unsigned compIdx, unsigned regionIdx = 0)
545  { return molarMass_[regionIdx][compIdx]; }
546 
548  static bool isIdealMixture(unsigned /*phaseIdx*/)
549  {
550  // fugacity coefficients are only pressure dependent -> we
551  // have an ideal mixture
552  return true;
553  }
554 
556  static bool isCompressible(unsigned /*phaseIdx*/)
557  { return true; /* all phases are compressible */ }
558 
560  static bool isIdealGas(unsigned /*phaseIdx*/)
561  { return false; }
562 
563 
564  /****************************************
565  * Black-oil specific properties
566  ****************************************/
572  static size_t numRegions()
573  { return molarMass_.size(); }
574 
581  static bool enableDissolvedGas()
582  { return enableDissolvedGas_; }
583 
590  static bool enableVaporizedOil()
591  { return enableVaporizedOil_; }
592 
598  static bool enableDiffusion()
599  { return enableDiffusion_; }
600 
606  static Scalar referenceDensity(unsigned phaseIdx, unsigned regionIdx)
607  { return referenceDensity_[regionIdx][phaseIdx]; }
608 
609  /****************************************
610  * thermodynamic quantities (generic version)
611  ****************************************/
613  template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
614  static LhsEval density(const FluidState& fluidState,
615  const ParameterCache<ParamCacheEval>& paramCache,
616  unsigned phaseIdx)
617  { return density<FluidState, LhsEval>(fluidState, phaseIdx, paramCache.regionIndex()); }
618 
620  template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
621  static LhsEval fugacityCoefficient(const FluidState& fluidState,
622  const ParameterCache<ParamCacheEval>& paramCache,
623  unsigned phaseIdx,
624  unsigned compIdx)
625  {
626  return fugacityCoefficient<FluidState, LhsEval>(fluidState,
627  phaseIdx,
628  compIdx,
629  paramCache.regionIndex());
630  }
631 
633  template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
634  static LhsEval viscosity(const FluidState& fluidState,
635  const ParameterCache<ParamCacheEval>& paramCache,
636  unsigned phaseIdx)
637  { return viscosity<FluidState, LhsEval>(fluidState, phaseIdx, paramCache.regionIndex()); }
638 
640  template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
641  static LhsEval enthalpy(const FluidState& fluidState,
642  const ParameterCache<ParamCacheEval>& paramCache,
643  unsigned phaseIdx)
644  { return enthalpy<FluidState, LhsEval>(fluidState, phaseIdx, paramCache.regionIndex()); }
645 
646  /****************************************
647  * thermodynamic quantities (black-oil specific version: Note that the PVT region
648  * index is explicitly passed instead of a parameter cache object)
649  ****************************************/
651  template <class FluidState, class LhsEval = typename FluidState::Scalar>
652  static LhsEval density(const FluidState& fluidState,
653  unsigned phaseIdx,
654  unsigned regionIdx)
655  {
656  assert(phaseIdx <= numPhases);
657  assert(regionIdx <= numRegions());
658 
659  const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
660  const LhsEval& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
661  const LhsEval& saltConcentration = BlackOil::template getSaltConcentration_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
662 
663  switch (phaseIdx) {
664  case oilPhaseIdx: {
665  if (enableDissolvedGas()) {
666  // miscible oil
667  const LhsEval& Rs = BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
668  const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
669 
670  return
671  bo*referenceDensity(oilPhaseIdx, regionIdx)
672  + Rs*bo*referenceDensity(gasPhaseIdx, regionIdx);
673  }
674 
675  // immiscible oil
676  const LhsEval Rs(0.0);
677  const auto& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
678 
679  return referenceDensity(phaseIdx, regionIdx)*bo;
680  }
681 
682  case gasPhaseIdx: {
683  if (enableVaporizedOil()) {
684  // miscible gas
685  const LhsEval& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
686  const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
687 
688  return
689  bg*referenceDensity(gasPhaseIdx, regionIdx)
690  + Rv*bg*referenceDensity(oilPhaseIdx, regionIdx);
691  }
692 
693  // immiscible gas
694  const LhsEval Rv(0.0);
695  const auto& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
696  return bg*referenceDensity(phaseIdx, regionIdx);
697  }
698 
699  case waterPhaseIdx:
700  return
701  referenceDensity(waterPhaseIdx, regionIdx)
702  * waterPvt_->inverseFormationVolumeFactor(regionIdx, T, p, saltConcentration);
703  }
704 
705  throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
706  }
707 
715  template <class FluidState, class LhsEval = typename FluidState::Scalar>
716  static LhsEval saturatedDensity(const FluidState& fluidState,
717  unsigned phaseIdx,
718  unsigned regionIdx)
719  {
720  assert(phaseIdx <= numPhases);
721  assert(regionIdx <= numRegions());
722 
723  const auto& p = fluidState.pressure(phaseIdx);
724  const auto& T = fluidState.temperature(phaseIdx);
725 
726  switch (phaseIdx) {
727  case oilPhaseIdx: {
728  if (enableDissolvedGas()) {
729  // miscible oil
730  const LhsEval& Rs = saturatedDissolutionFactor<FluidState, LhsEval>(fluidState, oilPhaseIdx, regionIdx);
731  const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
732 
733  return
734  bo*referenceDensity(oilPhaseIdx, regionIdx)
735  + Rs*bo*referenceDensity(gasPhaseIdx, regionIdx);
736  }
737 
738  // immiscible oil
739  const LhsEval Rs(0.0);
740  const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
741  return referenceDensity(phaseIdx, regionIdx)*bo;
742  }
743 
744  case gasPhaseIdx: {
745  if (enableVaporizedOil()) {
746  // miscible gas
747  const LhsEval& Rv = saturatedDissolutionFactor<FluidState, LhsEval>(fluidState, gasPhaseIdx, regionIdx);
748  const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
749 
750  return
751  bg*referenceDensity(gasPhaseIdx, regionIdx)
752  + Rv*bg*referenceDensity(oilPhaseIdx, regionIdx);
753  }
754 
755  // immiscible gas
756  const LhsEval Rv(0.0);
757  const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
758 
759  return referenceDensity(phaseIdx, regionIdx)*bg;
760 
761  }
762 
763  case waterPhaseIdx:
764  return
765  referenceDensity(waterPhaseIdx, regionIdx)
766  *inverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, waterPhaseIdx, regionIdx);
767  }
768 
769  throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
770  }
771 
780  template <class FluidState, class LhsEval = typename FluidState::Scalar>
781  static LhsEval inverseFormationVolumeFactor(const FluidState& fluidState,
782  unsigned phaseIdx,
783  unsigned regionIdx)
784  {
785  assert(phaseIdx <= numPhases);
786  assert(regionIdx <= numRegions());
787 
788  const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
789  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
790  const auto& saltConcentration = decay<LhsEval>(fluidState.saltConcentration());
791 
792  switch (phaseIdx) {
793  case oilPhaseIdx: {
794  if (enableDissolvedGas()) {
795  const auto& Rs = BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
796  if (fluidState.saturation(gasPhaseIdx) > 0.0
797  && Rs >= (1.0 - 1e-10)*oilPvt_->saturatedGasDissolutionFactor(regionIdx, scalarValue(T), scalarValue(p)))
798  {
799  return oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
800  } else {
801  return oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
802  }
803  }
804 
805  const LhsEval Rs(0.0);
806  return oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
807  }
808  case gasPhaseIdx: {
809  if (enableVaporizedOil()) {
810  const auto& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
811  if (fluidState.saturation(oilPhaseIdx) > 0.0
812  && Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
813  {
814  return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
815  } else {
816  return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
817  }
818  }
819 
820  const LhsEval Rv(0.0);
821  return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv);
822  }
823  case waterPhaseIdx:
824  return waterPvt_->inverseFormationVolumeFactor(regionIdx, T, p, saltConcentration);
825  default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
826  }
827  }
828 
836  template <class FluidState, class LhsEval = typename FluidState::Scalar>
837  static LhsEval saturatedInverseFormationVolumeFactor(const FluidState& fluidState,
838  unsigned phaseIdx,
839  unsigned regionIdx)
840  {
841  assert(phaseIdx <= numPhases);
842  assert(regionIdx <= numRegions());
843 
844  const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
845  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
846  const auto& saltConcentration = decay<LhsEval>(fluidState.saltConcentration());
847 
848  switch (phaseIdx) {
849  case oilPhaseIdx: return oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
850  case gasPhaseIdx: return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
851  case waterPhaseIdx: return waterPvt_->inverseFormationVolumeFactor(regionIdx, T, p, saltConcentration);
852  default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
853  }
854  }
855 
857  template <class FluidState, class LhsEval = typename FluidState::Scalar>
858  static LhsEval fugacityCoefficient(const FluidState& fluidState,
859  unsigned phaseIdx,
860  unsigned compIdx,
861  unsigned regionIdx)
862  {
863  assert(phaseIdx <= numPhases);
864  assert(compIdx <= numComponents);
865  assert(regionIdx <= numRegions());
866 
867  const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
868  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
869 
870  // for the fugacity coefficient of the oil component in the oil phase, we use
871  // some pseudo-realistic value for the vapor pressure to ease physical
872  // interpretation of the results
873  const LhsEval phi_oO = 20e3/p;
874 
875  // for the gas component in the gas phase, assume it to be an ideal gas
876  const Scalar phi_gG = 1.0;
877 
878  // for the fugacity coefficient of the water component in the water phase, we use
879  // the same approach as for the oil component in the oil phase
880  const LhsEval phi_wW = 30e3/p;
881 
882  switch (phaseIdx) {
883  case gasPhaseIdx: // fugacity coefficients for all components in the gas phase
884  switch (compIdx) {
885  case gasCompIdx:
886  return phi_gG;
887 
888  // for the oil component, we calculate the Rv value for saturated gas and Rs
889  // for saturated oil, and compute the fugacity coefficients at the
890  // equilibrium. for this, we assume that in equilibrium the fugacities of the
891  // oil component is the same in both phases.
892  case oilCompIdx: {
893  if (!enableVaporizedOil())
894  // if there's no vaporized oil, the gas phase is assumed to be
895  // immiscible with the oil component
896  return phi_gG*1e6;
897 
898  const auto& R_vSat = gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p);
899  const auto& X_gOSat = convertRvToXgO(R_vSat, regionIdx);
900  const auto& x_gOSat = convertXgOToxgO(X_gOSat, regionIdx);
901 
902  const auto& R_sSat = oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p);
903  const auto& X_oGSat = convertRsToXoG(R_sSat, regionIdx);
904  const auto& x_oGSat = convertXoGToxoG(X_oGSat, regionIdx);
905  const auto& x_oOSat = 1.0 - x_oGSat;
906 
907  const auto& p_o = decay<LhsEval>(fluidState.pressure(oilPhaseIdx));
908  const auto& p_g = decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
909 
910  return phi_oO*p_o*x_oOSat / (p_g*x_gOSat);
911  }
912 
913  case waterCompIdx:
914  // the water component is assumed to be never miscible with the gas phase
915  return phi_gG*1e6;
916 
917  default:
918  throw std::logic_error("Invalid component index "+std::to_string(compIdx));
919  }
920 
921  case oilPhaseIdx: // fugacity coefficients for all components in the oil phase
922  switch (compIdx) {
923  case oilCompIdx:
924  return phi_oO;
925 
926  // for the oil and water components, we proceed analogous to the gas and
927  // water components in the gas phase
928  case gasCompIdx: {
929  if (!enableDissolvedGas())
930  // if there's no dissolved gas, the oil phase is assumed to be
931  // immiscible with the gas component
932  return phi_oO*1e6;
933 
934  const auto& R_vSat = gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p);
935  const auto& X_gOSat = convertRvToXgO(R_vSat, regionIdx);
936  const auto& x_gOSat = convertXgOToxgO(X_gOSat, regionIdx);
937  const auto& x_gGSat = 1.0 - x_gOSat;
938 
939  const auto& R_sSat = oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p);
940  const auto& X_oGSat = convertRsToXoG(R_sSat, regionIdx);
941  const auto& x_oGSat = convertXoGToxoG(X_oGSat, regionIdx);
942 
943  const auto& p_o = decay<LhsEval>(fluidState.pressure(oilPhaseIdx));
944  const auto& p_g = decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
945 
946  return phi_gG*p_g*x_gGSat / (p_o*x_oGSat);
947  }
948 
949  case waterCompIdx:
950  return phi_oO*1e6;
951 
952  default:
953  throw std::logic_error("Invalid component index "+std::to_string(compIdx));
954  }
955 
956  case waterPhaseIdx: // fugacity coefficients for all components in the water phase
957  // the water phase fugacity coefficients are pretty simple: because the water
958  // phase is assumed to consist entirely from the water component, we just
959  // need to make sure that the fugacity coefficients for the other components
960  // are a few orders of magnitude larger than the one of the water
961  // component. (i.e., the affinity of the gas and oil components to the water
962  // phase is lower by a few orders of magnitude)
963  switch (compIdx) {
964  case waterCompIdx: return phi_wW;
965  case oilCompIdx: return 1.1e6*phi_wW;
966  case gasCompIdx: return 1e6*phi_wW;
967  default:
968  throw std::logic_error("Invalid component index "+std::to_string(compIdx));
969  }
970 
971  default:
972  throw std::logic_error("Invalid phase index "+std::to_string(phaseIdx));
973  }
974 
975  throw std::logic_error("Unhandled phase or component index");
976  }
977 
979  template <class FluidState, class LhsEval = typename FluidState::Scalar>
980  static LhsEval viscosity(const FluidState& fluidState,
981  unsigned phaseIdx,
982  unsigned regionIdx)
983  {
984  assert(phaseIdx <= numPhases);
985  assert(regionIdx <= numRegions());
986 
987  const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
988  const LhsEval& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
989  const LhsEval& saltConcentration = BlackOil::template getSaltConcentration_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
990 
991  switch (phaseIdx) {
992  case oilPhaseIdx: {
993  if (enableDissolvedGas()) {
994  const auto& Rs = BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
995  if (fluidState.saturation(gasPhaseIdx) > 0.0
996  && Rs >= (1.0 - 1e-10)*oilPvt_->saturatedGasDissolutionFactor(regionIdx, scalarValue(T), scalarValue(p)))
997  {
998  return oilPvt_->saturatedViscosity(regionIdx, T, p);
999  } else {
1000  return oilPvt_->viscosity(regionIdx, T, p, Rs);
1001  }
1002  }
1003 
1004  const LhsEval Rs(0.0);
1005  return oilPvt_->viscosity(regionIdx, T, p, Rs);
1006  }
1007 
1008  case gasPhaseIdx: {
1009  if (enableVaporizedOil()) {
1010  const auto& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
1011  if (fluidState.saturation(oilPhaseIdx) > 0.0
1012  && Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
1013  {
1014  return gasPvt_->saturatedViscosity(regionIdx, T, p);
1015  } else {
1016  return gasPvt_->viscosity(regionIdx, T, p, Rv);
1017  }
1018  }
1019 
1020  const LhsEval Rv(0.0);
1021  return gasPvt_->viscosity(regionIdx, T, p, Rv);
1022  }
1023 
1024  case waterPhaseIdx:
1025  // since water is always assumed to be immiscible in the black-oil model,
1026  // there is no "saturated water"
1027  return waterPvt_->viscosity(regionIdx, T, p, saltConcentration);
1028  }
1029 
1030  throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
1031  }
1032 
1034  template <class FluidState, class LhsEval = typename FluidState::Scalar>
1035  static LhsEval enthalpy(const FluidState& fluidState,
1036  unsigned phaseIdx,
1037  unsigned regionIdx)
1038  {
1039  assert(phaseIdx <= numPhases);
1040  assert(regionIdx <= numRegions());
1041 
1042  const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
1043  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
1044 
1045  switch (phaseIdx) {
1046  case oilPhaseIdx:
1047  return
1048  oilPvt_->internalEnergy(regionIdx, T, p, BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
1049  + p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
1050 
1051  case gasPhaseIdx:
1052  return
1053  gasPvt_->internalEnergy(regionIdx, T, p, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
1054  + p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
1055 
1056  case waterPhaseIdx:
1057  return
1058  waterPvt_->internalEnergy(regionIdx, T, p)
1059  + p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
1060 
1061  default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
1062  }
1063 
1064  throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
1065  }
1066 
1073  template <class FluidState, class LhsEval = typename FluidState::Scalar>
1074  static LhsEval saturatedDissolutionFactor(const FluidState& fluidState,
1075  unsigned phaseIdx,
1076  unsigned regionIdx,
1077  const LhsEval& maxOilSaturation)
1078  {
1079  assert(phaseIdx <= numPhases);
1080  assert(regionIdx <= numRegions());
1081 
1082  const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
1083  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
1084  const auto& So = decay<LhsEval>(fluidState.saturation(oilPhaseIdx));
1085 
1086  switch (phaseIdx) {
1087  case oilPhaseIdx: return oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p, So, maxOilSaturation);
1088  case gasPhaseIdx: return gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p, So, maxOilSaturation);
1089  case waterPhaseIdx: return 0.0;
1090  default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
1091  }
1092  }
1093 
1102  template <class FluidState, class LhsEval = typename FluidState::Scalar>
1103  static LhsEval saturatedDissolutionFactor(const FluidState& fluidState,
1104  unsigned phaseIdx,
1105  unsigned regionIdx)
1106  {
1107  assert(phaseIdx <= numPhases);
1108  assert(regionIdx <= numRegions());
1109 
1110  const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
1111  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
1112 
1113  switch (phaseIdx) {
1114  case oilPhaseIdx: return oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p);
1115  case gasPhaseIdx: return gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p);
1116  case waterPhaseIdx: return 0.0;
1117  default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
1118  }
1119  }
1120 
1124  template <class FluidState, class LhsEval = typename FluidState::Scalar>
1125  static LhsEval bubblePointPressure(const FluidState& fluidState,
1126  unsigned regionIdx)
1127  {
1128  return saturationPressure(fluidState, oilPhaseIdx, regionIdx);
1129  }
1130 
1131 
1135  template <class FluidState, class LhsEval = typename FluidState::Scalar>
1136  static LhsEval dewPointPressure(const FluidState& fluidState,
1137  unsigned regionIdx)
1138  {
1139  return saturationPressure(fluidState, gasPhaseIdx, regionIdx);
1140  }
1141 
1152  template <class FluidState, class LhsEval = typename FluidState::Scalar>
1153  static LhsEval saturationPressure(const FluidState& fluidState,
1154  unsigned phaseIdx,
1155  unsigned regionIdx)
1156  {
1157  assert(phaseIdx <= numPhases);
1158  assert(regionIdx <= numRegions());
1159 
1160  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
1161 
1162  switch (phaseIdx) {
1163  case oilPhaseIdx: return oilPvt_->saturationPressure(regionIdx, T, BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx));
1164  case gasPhaseIdx: return gasPvt_->saturationPressure(regionIdx, T, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx));
1165  case waterPhaseIdx: return 0.0;
1166  default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
1167  }
1168  }
1169 
1170  /****************************************
1171  * Auxiliary and convenience methods for the black-oil model
1172  ****************************************/
1177  template <class LhsEval>
1178  static LhsEval convertXoGToRs(const LhsEval& XoG, unsigned regionIdx)
1179  {
1180  Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
1181  Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
1182 
1183  return XoG/(1.0 - XoG)*(rho_oRef/rho_gRef);
1184  }
1185 
1190  template <class LhsEval>
1191  static LhsEval convertXgOToRv(const LhsEval& XgO, unsigned regionIdx)
1192  {
1193  Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
1194  Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
1195 
1196  return XgO/(1.0 - XgO)*(rho_gRef/rho_oRef);
1197  }
1198 
1203  template <class LhsEval>
1204  static LhsEval convertRsToXoG(const LhsEval& Rs, unsigned regionIdx)
1205  {
1206  Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
1207  Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
1208 
1209  const LhsEval& rho_oG = Rs*rho_gRef;
1210  return rho_oG/(rho_oRef + rho_oG);
1211  }
1212 
1217  template <class LhsEval>
1218  static LhsEval convertRvToXgO(const LhsEval& Rv, unsigned regionIdx)
1219  {
1220  Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
1221  Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
1222 
1223  const LhsEval& rho_gO = Rv*rho_oRef;
1224  return rho_gO/(rho_gRef + rho_gO);
1225  }
1226 
1230  template <class LhsEval>
1231  static LhsEval convertXoGToxoG(const LhsEval& XoG, unsigned regionIdx)
1232  {
1233  Scalar MO = molarMass_[regionIdx][oilCompIdx];
1234  Scalar MG = molarMass_[regionIdx][gasCompIdx];
1235 
1236  return XoG*MO / (MG*(1 - XoG) + XoG*MO);
1237  }
1238 
1242  template <class LhsEval>
1243  static LhsEval convertxoGToXoG(const LhsEval& xoG, unsigned regionIdx)
1244  {
1245  Scalar MO = molarMass_[regionIdx][oilCompIdx];
1246  Scalar MG = molarMass_[regionIdx][gasCompIdx];
1247 
1248  return xoG*MG / (xoG*(MG - MO) + MO);
1249  }
1250 
1254  template <class LhsEval>
1255  static LhsEval convertXgOToxgO(const LhsEval& XgO, unsigned regionIdx)
1256  {
1257  Scalar MO = molarMass_[regionIdx][oilCompIdx];
1258  Scalar MG = molarMass_[regionIdx][gasCompIdx];
1259 
1260  return XgO*MG / (MO*(1 - XgO) + XgO*MG);
1261  }
1262 
1266  template <class LhsEval>
1267  static LhsEval convertxgOToXgO(const LhsEval& xgO, unsigned regionIdx)
1268  {
1269  Scalar MO = molarMass_[regionIdx][oilCompIdx];
1270  Scalar MG = molarMass_[regionIdx][gasCompIdx];
1271 
1272  return xgO*MO / (xgO*(MO - MG) + MG);
1273  }
1274 
1282  static const GasPvt& gasPvt()
1283  { return *gasPvt_; }
1284 
1292  static const OilPvt& oilPvt()
1293  { return *oilPvt_; }
1294 
1302  static const WaterPvt& waterPvt()
1303  { return *waterPvt_; }
1304 
1310  static Scalar reservoirTemperature(unsigned = 0)
1311  { return reservoirTemperature_; }
1312 
1318  static void setReservoirTemperature(Scalar value)
1319  { reservoirTemperature_ = value; }
1320 
1321  static short activeToCanonicalPhaseIdx(unsigned activePhaseIdx) {
1322  assert(activePhaseIdx<numActivePhases());
1323  return activeToCanonicalPhaseIdx_[activePhaseIdx];
1324  }
1325 
1326  static short canonicalToActivePhaseIdx(unsigned phaseIdx) {
1327  assert(phaseIdx<numPhases);
1328  assert(phaseIsActive(phaseIdx));
1329  return canonicalToActivePhaseIdx_[phaseIdx];
1330  }
1331 
1333  static Scalar diffusionCoefficient(unsigned compIdx, unsigned phaseIdx, unsigned regionIdx = 0)
1334  { return diffusionCoefficients_[regionIdx][numPhases*compIdx + phaseIdx]; }
1335 
1337  static void setDiffusionCoefficient(Scalar coefficient, unsigned compIdx, unsigned phaseIdx, unsigned regionIdx = 0)
1338  { diffusionCoefficients_[regionIdx][numPhases*compIdx + phaseIdx] = coefficient ; }
1339 
1343  template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
1344  static LhsEval diffusionCoefficient(const FluidState& fluidState,
1345  const ParameterCache<ParamCacheEval>& paramCache,
1346  unsigned phaseIdx,
1347  unsigned compIdx)
1348  {
1349  // diffusion is disabled by the user
1350  if(!enableDiffusion())
1351  return 0.0;
1352 
1353  // diffusion coefficients are set, and we use them
1354  if(!diffusionCoefficients_.empty()) {
1355  return diffusionCoefficient(compIdx, phaseIdx, paramCache.regionIndex());
1356  }
1357 
1358  const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
1359  const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
1360 
1361  switch (phaseIdx) {
1362  case oilPhaseIdx: return oilPvt().diffusionCoefficient(T, p, compIdx);
1363  case gasPhaseIdx: return gasPvt().diffusionCoefficient(T, p, compIdx);
1364  case waterPhaseIdx: return 0.0;
1365  default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
1366  }
1367  }
1368 
1369 private:
1370  static void resizeArrays_(size_t numRegions)
1371  {
1372  molarMass_.resize(numRegions);
1373  referenceDensity_.resize(numRegions);
1374  }
1375 
1376  static Scalar reservoirTemperature_;
1377 
1378  static std::shared_ptr<GasPvt> gasPvt_;
1379  static std::shared_ptr<OilPvt> oilPvt_;
1380  static std::shared_ptr<WaterPvt> waterPvt_;
1381 
1382  static bool enableDissolvedGas_;
1383  static bool enableVaporizedOil_;
1384  static bool enableDiffusion_;
1385 
1386  // HACK for GCC 4.4: the array size has to be specified using the literal value '3'
1387  // here, because GCC 4.4 seems to be unable to determine the number of phases from
1388  // the BlackOil fluid system in the attribute declaration below...
1389  static std::vector<std::array<Scalar, /*numPhases=*/3> > referenceDensity_;
1390  static std::vector<std::array<Scalar, /*numComponents=*/3> > molarMass_;
1391  static std::vector<std::array<Scalar, /*numComponents=*/3 * /*numPhases=*/3> > diffusionCoefficients_;
1392 
1393  static std::array<short, numPhases> activeToCanonicalPhaseIdx_;
1394  static std::array<short, numPhases> canonicalToActivePhaseIdx_;
1395 
1396  static bool isInitialized_;
1397 };
1398 
1399 template <class Scalar, class IndexTraits>
1400 unsigned char BlackOilFluidSystem<Scalar, IndexTraits>::numActivePhases_;
1401 
1402 template <class Scalar, class IndexTraits>
1403 std::array<bool, BlackOilFluidSystem<Scalar, IndexTraits>::numPhases> BlackOilFluidSystem<Scalar, IndexTraits>::phaseIsActive_;
1404 
1405 template <class Scalar, class IndexTraits>
1406 std::array<short, BlackOilFluidSystem<Scalar, IndexTraits>::numPhases> BlackOilFluidSystem<Scalar, IndexTraits>::activeToCanonicalPhaseIdx_;
1407 
1408 template <class Scalar, class IndexTraits>
1409 std::array<short, BlackOilFluidSystem<Scalar, IndexTraits>::numPhases> BlackOilFluidSystem<Scalar, IndexTraits>::canonicalToActivePhaseIdx_;
1410 
1411 template <class Scalar, class IndexTraits>
1412 Scalar
1414 
1415 template <class Scalar, class IndexTraits>
1416 Scalar
1418 
1419 template <class Scalar, class IndexTraits>
1420 Scalar
1421 BlackOilFluidSystem<Scalar, IndexTraits>::reservoirTemperature_;
1422 
1423 template <class Scalar, class IndexTraits>
1424 bool BlackOilFluidSystem<Scalar, IndexTraits>::enableDissolvedGas_;
1425 
1426 template <class Scalar, class IndexTraits>
1427 bool BlackOilFluidSystem<Scalar, IndexTraits>::enableVaporizedOil_;
1428 
1429 template <class Scalar, class IndexTraits>
1430 bool BlackOilFluidSystem<Scalar, IndexTraits>::enableDiffusion_;
1431 
1432 template <class Scalar, class IndexTraits>
1433 std::shared_ptr<OilPvtMultiplexer<Scalar> >
1434 BlackOilFluidSystem<Scalar, IndexTraits>::oilPvt_;
1435 
1436 template <class Scalar, class IndexTraits>
1437 std::shared_ptr<GasPvtMultiplexer<Scalar> >
1438 BlackOilFluidSystem<Scalar, IndexTraits>::gasPvt_;
1439 
1440 template <class Scalar, class IndexTraits>
1441 std::shared_ptr<WaterPvtMultiplexer<Scalar> >
1442 BlackOilFluidSystem<Scalar, IndexTraits>::waterPvt_;
1443 
1444 template <class Scalar, class IndexTraits>
1445 std::vector<std::array<Scalar, 3> >
1446 BlackOilFluidSystem<Scalar, IndexTraits>::referenceDensity_;
1447 
1448 template <class Scalar, class IndexTraits>
1449 std::vector<std::array<Scalar, 3> >
1450 BlackOilFluidSystem<Scalar, IndexTraits>::molarMass_;
1451 
1452 template <class Scalar, class IndexTraits>
1453 std::vector<std::array<Scalar, 9> >
1454 BlackOilFluidSystem<Scalar, IndexTraits>::diffusionCoefficients_;
1455 
1456 template <class Scalar, class IndexTraits>
1457 bool BlackOilFluidSystem<Scalar, IndexTraits>::isInitialized_ = false;
1458 
1459 } // namespace Opm
1460 
1461 #endif
The base class for all fluid systems.
The class which specifies the default phase and component indices for the black-oil fluid system.
This class represents the Pressure-Volume-Temperature relations of the liquid phase for a CO2-Brine s...
A central place for various physical constants occuring in some equations.
Provides the opm-material specific exception classes.
This class represents the Pressure-Volume-Temperature relations of the gas phase in the black-oil mod...
This macro generates a class HasMember_${MEMBER_NAME} which can be used for template specialization.
#define OPM_GENERATE_HAS_MEMBER(MEMBER_NAME,...)
This macro generates a class HasMember_${MEMBER_NAME} which can be used for template specialization.
Definition: HasMemberGeneratorMacros.hpp:49
A traits class which provides basic mathematical functions for arbitrary scalar floating point values...
This class represents the Pressure-Volume-Temperature relations of the oil phase in the black-oil mod...
Some templates to wrap the valgrind client request macros.
This class represents the Pressure-Volume-Temperature relations of the water phase in the black-oil m...
The base class for all fluid systems.
Definition: BaseFluidSystem.hpp:44
Scalar Scalar
The type used for scalar quantities.
Definition: BaseFluidSystem.hpp:49
A fluid system which uses the black-oil model assumptions to calculate termodynamically meaningful qu...
Definition: BlackOilFluidSystem.hpp:112
static unsigned numActivePhases()
Returns the number of active fluid phases (i.e., usually three)
Definition: BlackOilFluidSystem.hpp:485
static Scalar diffusionCoefficient(unsigned compIdx, unsigned phaseIdx, unsigned regionIdx=0)
Calculate the binary molecular diffusion coefficient for a component in a fluid phase [mol^2 * s / (k...
Definition: BlackOilFluidSystem.hpp:1333
static void setReferenceDensities(Scalar rhoOil, Scalar rhoWater, Scalar rhoGas, unsigned regionIdx)
Initialize the values of the reference densities.
Definition: BlackOilFluidSystem.hpp:368
static LhsEval convertXoGToRs(const LhsEval &XoG, unsigned regionIdx)
Convert the mass fraction of the gas component in the oil phase to the corresponding gas dissolution ...
Definition: BlackOilFluidSystem.hpp:1178
static LhsEval convertRsToXoG(const LhsEval &Rs, unsigned regionIdx)
Convert a gas dissolution factor to the the corresponding mass fraction of the gas component in the o...
Definition: BlackOilFluidSystem.hpp:1204
static void initEnd()
Finish initializing the black oil fluid system.
Definition: BlackOilFluidSystem.hpp:382
static constexpr unsigned soluteComponentIndex(unsigned phaseIdx)
returns the index of "secondary" component of a phase (solute)
Definition: BlackOilFluidSystem.hpp:512
static LhsEval saturatedInverseFormationVolumeFactor(const FluidState &fluidState, unsigned phaseIdx, unsigned regionIdx)
Returns the formation volume factor of a "saturated" fluid phase.
Definition: BlackOilFluidSystem.hpp:837
static LhsEval convertXoGToxoG(const LhsEval &XoG, unsigned regionIdx)
Convert a gas mass fraction in the oil phase the corresponding mole fraction.
Definition: BlackOilFluidSystem.hpp:1231
static LhsEval viscosity(const FluidState &fluidState, unsigned phaseIdx, unsigned regionIdx)
Calculate the dynamic viscosity of a fluid phase [Pa*s].
Definition: BlackOilFluidSystem.hpp:980
static const unsigned oilPhaseIdx
Index of the oil phase.
Definition: BlackOilFluidSystem.hpp:432
static Scalar reservoirTemperature(unsigned=0)
Set the temperature of the reservoir.
Definition: BlackOilFluidSystem.hpp:1310
static void setEnableDiffusion(bool yesno)
Specify whether the fluid system should consider diffusion.
Definition: BlackOilFluidSystem.hpp:339
static const unsigned waterPhaseIdx
Index of the water phase.
Definition: BlackOilFluidSystem.hpp:430
static LhsEval bubblePointPressure(const FluidState &fluidState, unsigned regionIdx)
Returns the bubble point pressure $P_b$ using the current Rs.
Definition: BlackOilFluidSystem.hpp:1125
static LhsEval density(const FluidState &fluidState, unsigned phaseIdx, unsigned regionIdx)
Calculate the density [kg/m^3] of a fluid phase.
Definition: BlackOilFluidSystem.hpp:652
static LhsEval saturatedDensity(const FluidState &fluidState, unsigned phaseIdx, unsigned regionIdx)
Compute the density of a saturated fluid phase.
Definition: BlackOilFluidSystem.hpp:716
static LhsEval viscosity(const FluidState &fluidState, const ParameterCache< ParamCacheEval > &paramCache, unsigned phaseIdx)
Calculate the dynamic viscosity of a fluid phase [Pa*s].
Definition: BlackOilFluidSystem.hpp:634
static bool isIdealMixture(unsigned)
Returns true if and only if a fluid phase is assumed to be an ideal mixture.
Definition: BlackOilFluidSystem.hpp:548
static const unsigned numComponents
Number of chemical species in the fluid system.
Definition: BlackOilFluidSystem.hpp:470
static const char * componentName(unsigned compIdx)
Return the human readable name of a component.
Definition: BlackOilFluidSystem.hpp:528
static const GasPvt & gasPvt()
Return a reference to the low-level object which calculates the gas phase quantities.
Definition: BlackOilFluidSystem.hpp:1282
static constexpr unsigned solventComponentIndex(unsigned phaseIdx)
returns the index of "primary" component of a phase (solvent)
Definition: BlackOilFluidSystem.hpp:496
static void setWaterPvt(std::shared_ptr< WaterPvt > pvtObj)
Set the pressure-volume-saturation (PVT) relations for the water phase.
Definition: BlackOilFluidSystem.hpp:358
static void setReservoirTemperature(Scalar value)
Return the temperature of the reservoir.
Definition: BlackOilFluidSystem.hpp:1318
static LhsEval saturatedDissolutionFactor(const FluidState &fluidState, unsigned phaseIdx, unsigned regionIdx, const LhsEval &maxOilSaturation)
Returns the dissolution factor of a saturated fluid phase.
Definition: BlackOilFluidSystem.hpp:1074
static size_t numRegions()
Returns the number of PVT regions which are considered.
Definition: BlackOilFluidSystem.hpp:572
static LhsEval enthalpy(const FluidState &fluidState, unsigned phaseIdx, unsigned regionIdx)
Given a phase's composition, temperature, pressure and density, calculate its specific enthalpy [J/kg...
Definition: BlackOilFluidSystem.hpp:1035
static bool isLiquid(unsigned phaseIdx)
Return whether a phase is liquid.
Definition: BlackOilFluidSystem.hpp:459
static LhsEval enthalpy(const FluidState &fluidState, const ParameterCache< ParamCacheEval > &paramCache, unsigned phaseIdx)
Given a phase's composition, temperature, pressure and density, calculate its specific enthalpy [J/kg...
Definition: BlackOilFluidSystem.hpp:641
static LhsEval fugacityCoefficient(const FluidState &fluidState, unsigned phaseIdx, unsigned compIdx, unsigned regionIdx)
Calculate the fugacity coefficient [Pa] of an individual component in a fluid phase.
Definition: BlackOilFluidSystem.hpp:858
static LhsEval convertxoGToXoG(const LhsEval &xoG, unsigned regionIdx)
Convert a gas mole fraction in the oil phase the corresponding mass fraction.
Definition: BlackOilFluidSystem.hpp:1243
static Scalar molarMass(unsigned compIdx, unsigned regionIdx=0)
Return the molar mass of a component in [kg/mol].
Definition: BlackOilFluidSystem.hpp:544
static LhsEval saturatedDissolutionFactor(const FluidState &fluidState, unsigned phaseIdx, unsigned regionIdx)
Returns the dissolution factor of a saturated fluid phase.
Definition: BlackOilFluidSystem.hpp:1103
static LhsEval dewPointPressure(const FluidState &fluidState, unsigned regionIdx)
Returns the dew point pressure $P_d$ using the current Rv.
Definition: BlackOilFluidSystem.hpp:1136
static LhsEval saturationPressure(const FluidState &fluidState, unsigned phaseIdx, unsigned regionIdx)
Returns the saturation pressure of a given phase [Pa] depending on its composition.
Definition: BlackOilFluidSystem.hpp:1153
static LhsEval convertXgOToxgO(const LhsEval &XgO, unsigned regionIdx)
Convert a oil mass fraction in the gas phase the corresponding mole fraction.
Definition: BlackOilFluidSystem.hpp:1255
static LhsEval diffusionCoefficient(const FluidState &fluidState, const ParameterCache< ParamCacheEval > &paramCache, unsigned phaseIdx, unsigned compIdx)
Calculate the binary molecular diffusion coefficient for a component in a fluid phase [mol^2 * s / (k...
Definition: BlackOilFluidSystem.hpp:1344
static bool enableDissolvedGas()
Returns whether the fluid system should consider that the gas component can dissolve in the oil phase...
Definition: BlackOilFluidSystem.hpp:581
static Scalar referenceDensity(unsigned phaseIdx, unsigned regionIdx)
Returns the density of a fluid phase at surface pressure [kg/m^3].
Definition: BlackOilFluidSystem.hpp:606
static Scalar surfaceTemperature
The temperature at the surface.
Definition: BlackOilFluidSystem.hpp:440
static void setOilPvt(std::shared_ptr< OilPvt > pvtObj)
Set the pressure-volume-saturation (PVT) relations for the oil phase.
Definition: BlackOilFluidSystem.hpp:352
static LhsEval convertxgOToXgO(const LhsEval &xgO, unsigned regionIdx)
Convert a oil mole fraction in the gas phase the corresponding mass fraction.
Definition: BlackOilFluidSystem.hpp:1267
static LhsEval inverseFormationVolumeFactor(const FluidState &fluidState, unsigned phaseIdx, unsigned regionIdx)
Returns the formation volume factor of an "undersaturated" fluid phase.
Definition: BlackOilFluidSystem.hpp:781
static void setEnableDissolvedGas(bool yesno)
Specify whether the fluid system should consider that the gas component can dissolve in the oil phase...
Definition: BlackOilFluidSystem.hpp:322
static const char * phaseName(unsigned phaseIdx)
Return the human readable name of a fluid phase.
Definition: BlackOilFluidSystem.hpp:443
static const unsigned numPhases
Number of fluid phases in the fluid system.
Definition: BlackOilFluidSystem.hpp:427
static LhsEval convertXgOToRv(const LhsEval &XgO, unsigned regionIdx)
Convert the mass fraction of the oil component in the gas phase to the corresponding oil vaporization...
Definition: BlackOilFluidSystem.hpp:1191
static LhsEval fugacityCoefficient(const FluidState &fluidState, const ParameterCache< ParamCacheEval > &paramCache, unsigned phaseIdx, unsigned compIdx)
Calculate the fugacity coefficient [Pa] of an individual component in a fluid phase.
Definition: BlackOilFluidSystem.hpp:621
static const unsigned gasPhaseIdx
Index of the gas phase.
Definition: BlackOilFluidSystem.hpp:434
static LhsEval density(const FluidState &fluidState, const ParameterCache< ParamCacheEval > &paramCache, unsigned phaseIdx)
Calculate the density [kg/m^3] of a fluid phase.
Definition: BlackOilFluidSystem.hpp:614
static LhsEval convertRvToXgO(const LhsEval &Rv, unsigned regionIdx)
Convert an oil vaporization factor to the corresponding mass fraction of the oil component in the gas...
Definition: BlackOilFluidSystem.hpp:1218
static bool enableVaporizedOil()
Returns whether the fluid system should consider that the oil component can dissolve in the gas phase...
Definition: BlackOilFluidSystem.hpp:590
static bool isIdealGas(unsigned)
Returns true if and only if a fluid phase is assumed to be an ideal gas.
Definition: BlackOilFluidSystem.hpp:560
static unsigned phaseIsActive(unsigned phaseIdx)
Returns whether a fluid phase is active.
Definition: BlackOilFluidSystem.hpp:489
static Scalar surfacePressure
The pressure at the surface.
Definition: BlackOilFluidSystem.hpp:437
static const WaterPvt & waterPvt()
Return a reference to the low-level object which calculates the water phase quantities.
Definition: BlackOilFluidSystem.hpp:1302
static void initBegin(size_t numPvtRegions)
Begin the initialization of the black oil fluid system.
Definition: BlackOilFluidSystem.hpp:294
static bool enableDiffusion()
Returns whether the fluid system should consider diffusion.
Definition: BlackOilFluidSystem.hpp:598
static const unsigned gasCompIdx
Index of the gas component.
Definition: BlackOilFluidSystem.hpp:477
static const unsigned oilCompIdx
Index of the oil component.
Definition: BlackOilFluidSystem.hpp:473
static void setEnableVaporizedOil(bool yesno)
Specify whether the fluid system should consider that the oil component can dissolve in the gas phase...
Definition: BlackOilFluidSystem.hpp:331
static const unsigned waterCompIdx
Index of the water component.
Definition: BlackOilFluidSystem.hpp:475
static void setDiffusionCoefficient(Scalar coefficient, unsigned compIdx, unsigned phaseIdx, unsigned regionIdx=0)
Definition: BlackOilFluidSystem.hpp:1337
static const OilPvt & oilPvt()
Return a reference to the low-level object which calculates the oil phase quantities.
Definition: BlackOilFluidSystem.hpp:1292
static bool isCompressible(unsigned)
Returns true if and only if a fluid phase is assumed to be compressible.
Definition: BlackOilFluidSystem.hpp:556
static void setGasPvt(std::shared_ptr< GasPvt > pvtObj)
Set the pressure-volume-saturation (PVT) relations for the gas phase.
Definition: BlackOilFluidSystem.hpp:346
static Scalar molarMass()
The molar mass in of the component.
Definition: Brine.hpp:80
static Scalar molarMass()
The mass in [kg] of one mole of CO2.
Definition: CO2.hpp:66
A central place for various physical constants occuring in some equations.
Definition: Constants.hpp:41
This class represents the Pressure-Volume-Temperature relations of the gas phase in the black-oil mod...
Definition: GasPvtMultiplexer.hpp:86
Evaluation diffusionCoefficient(const Evaluation &temperature, const Evaluation &pressure, unsigned compIdx) const
Calculate the binary molecular diffusion coefficient for a component in a fluid phase [mol^2 * s / (k...
Definition: GasPvtMultiplexer.hpp:276
A parameter cache which does nothing.
Definition: NullParameterCache.hpp:40
This class represents the Pressure-Volume-Temperature relations of the oil phase in the black-oil mod...
Definition: OilPvtMultiplexer.hpp:96
Evaluation diffusionCoefficient(const Evaluation &temperature, const Evaluation &pressure, unsigned compIdx) const
Calculate the binary molecular diffusion coefficient for a component in a fluid phase [mol^2 * s / (k...
Definition: OilPvtMultiplexer.hpp:271
This class represents the Pressure-Volume-Temperature relations of the water phase in the black-oil m...
Definition: WaterPvtMultiplexer.hpp:75
The type of the fluid system's parameter cache.
Definition: BlackOilFluidSystem.hpp:123
void assignPersistentData(const OtherCache &other)
Copy the data which is not dependent on the type of the Scalars from another parameter cache.
Definition: BlackOilFluidSystem.hpp:141
void setRegionIndex(unsigned val)
Set the index of the region which should be used to determine the thermodynamic properties.
Definition: BlackOilFluidSystem.hpp:164
unsigned regionIndex() const
Return the index of the region which should be used to determine the thermodynamic properties.
Definition: BlackOilFluidSystem.hpp:154