blackoilsolventmodules.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_SOLVENT_MODULE_HH
29 #define EWOMS_BLACK_OIL_SOLVENT_MODULE_HH
30 
31 #include "blackoilproperties.hh"
34 
35 #include <opm/material/fluidsystems/blackoilpvt/SolventPvt.hpp>
36 #include <opm/material/common/Tabulated1DFunction.hpp>
37 
38 #if HAVE_OPM_PARSER
39 #include <opm/parser/eclipse/Deck/Deck.hpp>
40 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
41 #include <opm/parser/eclipse/EclipseState/Tables/SsfnTable.hpp>
42 #include <opm/parser/eclipse/EclipseState/Tables/Sof2Table.hpp>
43 #include <opm/parser/eclipse/EclipseState/Tables/MsfnTable.hpp>
44 #include <opm/parser/eclipse/EclipseState/Tables/PmiscTable.hpp>
45 #include <opm/parser/eclipse/EclipseState/Tables/MiscTable.hpp>
46 #include <opm/parser/eclipse/EclipseState/Tables/SorwmisTable.hpp>
47 #include <opm/parser/eclipse/EclipseState/Tables/SgcwmisTable.hpp>
48 #include <opm/parser/eclipse/EclipseState/Tables/TlpmixpaTable.hpp>
49 #endif
50 
51 #include <opm/common/Valgrind.hpp>
52 #include <opm/common/Unused.hpp>
53 #include <opm/common/ErrorMacros.hpp>
54 #include <opm/common/Exceptions.hpp>
55 
56 #include <dune/common/fvector.hh>
57 
58 #include <string>
59 
60 namespace Ewoms {
66 template <class TypeTag, bool enableSolventV = GET_PROP_VALUE(TypeTag, EnableSolvent)>
67 class BlackOilSolventModule
68 {
69  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
70  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
71  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
72  typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
73  typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
74  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
75  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
76  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
77  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
78  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
79  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
80  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
81 
82  typedef Opm::MathToolbox<Evaluation> Toolbox;
83  typedef Opm::SolventPvt<Scalar> SolventPvt;
84 
85  typedef typename Opm::Tabulated1DFunction<Scalar> TabulatedFunction;
86 
87  static constexpr unsigned solventSaturationIdx = Indices::solventSaturationIdx;
88  static constexpr unsigned contiSolventEqIdx = Indices::contiSolventEqIdx;
89  static constexpr unsigned enableSolvent = enableSolventV;
90  static constexpr unsigned numEq = GET_PROP_VALUE(TypeTag, NumEq);
91  static constexpr unsigned numPhases = FluidSystem::numPhases;
92  static constexpr bool blackoilConserveSurfaceVolume = GET_PROP_VALUE(TypeTag, BlackoilConserveSurfaceVolume);
93 
94 
95 public:
96 #if HAVE_OPM_PARSER
97 
100  static void initFromDeck(const Opm::Deck& deck, const Opm::EclipseState& eclState)
101  {
102  // some sanity checks: if solvents are enabled, the SOLVENT keyword must be
103  // present, if solvents are disabled the keyword must not be present.
104  if (enableSolvent && !deck.hasKeyword("SOLVENT")) {
105  OPM_THROW(std::runtime_error,
106  "Non-trivial solvent treatment requested at compile time, but "
107  "the deck does not contain the SOLVENT keyword");
108  }
109  else if (!enableSolvent && deck.hasKeyword("SOLVENT")) {
110  OPM_THROW(std::runtime_error,
111  "Solvent treatment disabled at compile time, but the deck "
112  "contains the SOLVENT keyword");
113  }
114 
115  if (!deck.hasKeyword("SOLVENT"))
116  return; // solvent treatment is supposed to be disabled
117 
118  solventPvt_.initFromDeck(deck, eclState);
119 
120  const auto& tableManager = eclState.getTableManager();
121  // initialize the objects which deal with the SSFN keyword
122  const auto& ssfnTables = tableManager.getSsfnTables();
123  unsigned numSatRegions = tableManager.getTabdims().getNumSatTables();
124  setNumSatRegions(numSatRegions);
125  for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) {
126  const auto& ssfnTable = ssfnTables.template getTable<Opm::SsfnTable>(satRegionIdx);
127  ssfnKrg_[satRegionIdx].setXYContainers(ssfnTable.getSolventFractionColumn(),
128  ssfnTable.getGasRelPermMultiplierColumn(),
129  /*sortInput=*/true);
130  ssfnKrs_[satRegionIdx].setXYContainers(ssfnTable.getSolventFractionColumn(),
131  ssfnTable.getSolventRelPermMultiplierColumn(),
132  /*sortInput=*/true);
133  }
134 
135  // initialize the objects needed for miscible solvent and oil simulations
136  isMiscible_ = false;
137  if (deck.hasKeyword("MISCIBLE")) {
138  isMiscible_ = true;
139 
140  unsigned numSatRegions = tableManager.getTabdims().getNumSatTables();
141  unsigned numMiscRegions = 1;
142 
143  // misicible hydrocabon relative permeability wrt water
144  const auto& sof2Tables = tableManager.getSof2Tables();
145  if (!sof2Tables.empty()) {
146 
147  // resize the attributes of the object
148  sof2Krn_.resize(numSatRegions);
149  for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) {
150  const auto& sof2Table = sof2Tables.template getTable<Opm::Sof2Table>(satRegionIdx);
151  sof2Krn_[satRegionIdx].setXYContainers(sof2Table.getSoColumn(),
152  sof2Table.getKroColumn(),
153  /*sortInput=*/true);
154  }
155 
156  } else {
157  OPM_THROW(std::runtime_error, "SOF2 must be specified in MISCIBLE (SOLVENT) runs\n");
158  }
159 
160  const auto& miscTables = tableManager.getMiscTables();
161  if (!miscTables.empty()) {
162 
163  assert(numMiscRegions == miscTables.size());
164 
165  // resize the attributes of the object
166  misc_.resize(numMiscRegions);
167  for (unsigned miscRegionIdx = 0; miscRegionIdx < numMiscRegions; ++miscRegionIdx) {
168  const auto& miscTable = miscTables.template getTable<Opm::MiscTable>(miscRegionIdx);
169 
170  // solventFraction = Ss / (Ss + Sg);
171  const auto& solventFraction = miscTable.getSolventFractionColumn();
172  const auto& misc = miscTable.getMiscibilityColumn();
173  misc_[miscRegionIdx].setXYContainers(solventFraction, misc);
174 
175  }
176  } else {
177  OPM_THROW(std::runtime_error, "MISC must be specified in MISCIBLE (SOLVENT) runs\n");
178  }
179 
180  // resize the attributes of the object
181  pmisc_.resize(numMiscRegions);
182  const auto& pmiscTables = tableManager.getPmiscTables();
183  if (!pmiscTables.empty()) {
184 
185  assert(numMiscRegions == pmiscTables.size());
186 
187  for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
188  const auto& pmiscTable = pmiscTables.template getTable<Opm::PmiscTable>(regionIdx);
189 
190  // Copy data
191  const auto& po = pmiscTable.getOilPhasePressureColumn();
192  const auto& pmisc = pmiscTable.getMiscibilityColumn();
193 
194  pmisc_[regionIdx].setXYContainers(po, pmisc);
195 
196  }
197  } else {
198  std::vector<double> x = {0.0,1.0e20};
199  std::vector<double> y = {1.0,1.0};
200  TabulatedFunction constant = TabulatedFunction(2, x, y);
201  for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
202  setPmisc(regionIdx, constant);
203  }
204 
205  }
206 
207  // miscible relative permeability multipleiers
208  msfnKrsg_.resize(numSatRegions);
209  msfnKro_.resize(numSatRegions);
210  const auto& msfnTables = tableManager.getMsfnTables();
211  if (!msfnTables.empty()) {
212 
213  assert(numSatRegions == msfnTables.size());
214 
215 
216  for (unsigned regionIdx = 0; regionIdx < numSatRegions; ++regionIdx) {
217  const Opm::MsfnTable& msfnTable = msfnTables.template getTable<Opm::MsfnTable>(regionIdx);
218 
219  // Copy data
220  // Ssg = Ss + Sg;
221  const auto& Ssg = msfnTable.getGasPhaseFractionColumn();
222  const auto& krsg = msfnTable.getGasSolventRelpermMultiplierColumn();
223  const auto& kro = msfnTable.getOilRelpermMultiplierColumn();
224 
225  msfnKrsg_[regionIdx].setXYContainers(Ssg, krsg);
226  msfnKro_[regionIdx].setXYContainers(Ssg, kro);
227 
228  }
229  } else {
230  std::vector<double> x = {0.0,1.0};
231  std::vector<double> y = {1.0,0.0};
232  TabulatedFunction unit = TabulatedFunction(2, x, x);
233  TabulatedFunction inv_unit = TabulatedFunction(2, x, y);
234 
235  for (unsigned regionIdx = 0; regionIdx < numSatRegions; ++regionIdx) {
236  setMsfn(regionIdx, unit, inv_unit);
237  }
238  }
239  // resize the attributes of the object
240  sorwmis_.resize(numMiscRegions);
241  const auto& sorwmisTables = tableManager.getSorwmisTables();
242  if (!sorwmisTables.empty()) {
243  assert(numMiscRegions == sorwmisTables.size());
244 
245  for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
246  const auto& sorwmisTable = sorwmisTables.template getTable<Opm::SorwmisTable>(regionIdx);
247 
248  // Copy data
249  const auto& sw = sorwmisTable.getWaterSaturationColumn();
250  const auto& sorwmis = sorwmisTable.getMiscibleResidualOilColumn();
251 
252  sorwmis_[regionIdx].setXYContainers(sw, sorwmis);
253  }
254  } else {
255  // default
256  std::vector<double> x = {0.0,1.0};
257  std::vector<double> y = {0.0,0.0};
258  TabulatedFunction zero = TabulatedFunction(2, x, y);
259  for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
260  setSorwmis(regionIdx, zero);
261  }
262  }
263 
264  // resize the attributes of the object
265  sgcwmis_.resize(numMiscRegions);
266  const auto& sgcwmisTables = tableManager.getSgcwmisTables();
267  if (!sgcwmisTables.empty()) {
268 
269  assert(numMiscRegions ==sgcwmisTables.size());
270 
271  for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
272  const auto& sgcwmisTable = sgcwmisTables.template getTable<Opm::SgcwmisTable>(regionIdx);
273 
274  // Copy data
275  const auto& sw = sgcwmisTable.getWaterSaturationColumn();
276  const auto& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn();
277 
278  sgcwmis_[regionIdx].setXYContainers(sw, sgcwmis);
279  }
280  } else {
281  // default
282  std::vector<double> x = {0.0,1.0};
283  std::vector<double> y = {0.0,0.0};
284  TabulatedFunction zero = TabulatedFunction(2, x, y);
285  for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
286  setSgcmis(regionIdx, zero);
287  }
288  }
289 
290 
291  if (deck.hasKeyword("TLMIXPAR")) {
292  // resize the attributes of the object
293  tlMixParamViscosity_.resize(numMiscRegions);
294  tlMixParamDensity_.resize(numMiscRegions);
295 
296  assert(numMiscRegions == deck.getKeyword("TLMIXPAR").size() );
297 
298  for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
299  const auto& tlmixparRecord = deck.getKeyword("TLMIXPAR").getRecord(regionIdx);
300  const auto& mix_params_viscosity = tlmixparRecord.getItem("TL_VISCOSITY_PARAMETER").getSIDoubleData();
301  tlMixParamViscosity_[regionIdx] = mix_params_viscosity[0];
302  const auto& mix_params_density = tlmixparRecord.getItem("TL_DENSITY_PARAMETER").getSIDoubleData();
303  const int numDensityItems = mix_params_density.size();
304  if (numDensityItems == 0) {
305  tlMixParamDensity_[regionIdx] = tlMixParamViscosity_[regionIdx];
306  } else if (numDensityItems == 1) {
307  tlMixParamDensity_[regionIdx] = mix_params_density[0];
308  } else {
309  OPM_THROW(std::runtime_error, "Only one value can be entered for the TL parameter pr MISC region.");
310  }
311  }
312  } else {
313  OPM_THROW(std::runtime_error, "TLMIXPAR must be specified in MISCIBLE (SOLVENT) runs\n");
314  }
315 
316  // resize the attributes of the object
317  tlPMixTable_.resize(numMiscRegions);
318  if (deck.hasKeyword("TLPMIXPA")) {
319  const auto& tlpmixparTables = tableManager.getTlpmixpaTables();
320  if (!tlpmixparTables.empty()) {
321 
322  assert(numMiscRegions == tlpmixparTables.size());
323  for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
324  const auto& tlpmixparTable = tlpmixparTables.template getTable<Opm::TlpmixpaTable>(regionIdx);
325 
326  // Copy data
327  const auto& po = tlpmixparTable.getOilPhasePressureColumn();
328  const auto& tlpmixpa = tlpmixparTable.getMiscibilityColumn();
329 
330  tlPMixTable_[regionIdx].setXYContainers(po, tlpmixpa);
331 
332  }
333  } else {
334  // if empty keyword. Try to use the pmisc table as default.
335  if (pmisc_.size() > 0) {
336  tlPMixTable_ = pmisc_;
337  } else {
338  OPM_THROW(std::invalid_argument, "If the pressure dependent TL values in TLPMIXPA is defaulted (no entries), then the PMISC tables must be specified.");
339  }
340  }
341  } else {
342  // default
343  std::vector<double> x = {0.0,1.0e20};
344  std::vector<double> y = {1.0,1.0};
345  TabulatedFunction ones = TabulatedFunction(2, x, y);
346  for (unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
347  setTlpmixpa(regionIdx, ones);
348  }
349  }
350 
351 
352  }
353 
354 
355 
356  }
357 #endif
358 
364  static void setNumSatRegions(unsigned numRegions)
365  {
366  ssfnKrg_.resize(numRegions);
367  ssfnKrs_.resize(numRegions);
368  }
369 
375  static void setSsfn(unsigned satRegionIdx,
376  const TabulatedFunction& ssfnKrg,
377  const TabulatedFunction& ssfnKrs)
378  {
379  ssfnKrg_[satRegionIdx] = ssfnKrg;
380  ssfnKrs_[satRegionIdx] = ssfnKrs;
381  }
382 
388  static void setSof2(unsigned satRegionIdx,
389  const TabulatedFunction& sof2Krn)
390  {
391  sof2Krn_[satRegionIdx] = sof2Krn;
392  }
393 
399  static void setMisc(unsigned miscRegionIdx,
400  const TabulatedFunction& misc)
401  {
402  misc_[miscRegionIdx] = misc;
403  }
404 
410  static void setPmisc(unsigned miscRegionIdx,
411  const TabulatedFunction& pmisc)
412  {
413  pmisc_[miscRegionIdx] = pmisc;
414  }
415 
421  static void setMsfn(unsigned satRegionIdx,
422  const TabulatedFunction& msfnKrsg,
423  const TabulatedFunction& msfnKro)
424  {
425  msfnKrsg_[satRegionIdx] = msfnKrsg;
426  msfnKro_[satRegionIdx] = msfnKro;
427  }
428 
434  static void setSorwmis(unsigned miscRegionIdx,
435  const TabulatedFunction& sorwmis)
436  {
437  sorwmis_[miscRegionIdx] = sorwmis;
438  }
439 
445  static void setSgcmis(unsigned miscRegionIdx,
446  const TabulatedFunction& sgcwmis)
447  {
448  sgcwmis_[miscRegionIdx] = sgcwmis;
449  }
450 
456  static void setTlmixpar(unsigned miscRegionIdx,
457  const Scalar& tlMixParamViscosity,
458  const Scalar& tlMixParamDensity)
459  {
460  tlMixParamViscosity_[miscRegionIdx] = tlMixParamViscosity;
461  tlMixParamDensity_[miscRegionIdx] = tlMixParamDensity;
462  }
463 
469  static void setTlpmixpa(unsigned miscRegionIdx,
470  const TabulatedFunction& tlPMixTable)
471  {
472  tlPMixTable_[miscRegionIdx] = tlPMixTable;
473  }
474 
478  static void setSolventPvt(const SolventPvt& value)
479  { solventPvt_ = value; }
480 
481 
482  static void setIsMiscible(const bool isMiscible)
483  { isMiscible_ = isMiscible; }
484 
488  static void registerParameters()
489  {
490  if (!enableSolvent)
491  // solvents have disabled at compile time
492  return;
493 
495  }
496 
500  static void registerOutputModules(Model& model,
501  Simulator& simulator)
502  {
503  if (!enableSolvent)
504  // solvents have disabled at compile time
505  return;
506 
507  model.addOutputModule(new Ewoms::VtkBlackOilSolventModule<TypeTag>(simulator));
508  }
509 
510  static bool primaryVarApplies(unsigned pvIdx)
511  {
512  if (!enableSolvent)
513  // solvents have disabled at compile time
514  return false;
515 
516  return pvIdx == solventSaturationIdx;
517  }
518 
519  static std::string primaryVarName(unsigned pvIdx OPM_OPTIM_UNUSED)
520  {
521  assert(primaryVarApplies(pvIdx));
522 
523  return "saturation_solvent";
524  }
525 
526  static Scalar primaryVarWeight(unsigned pvIdx OPM_OPTIM_UNUSED)
527  {
528  assert(primaryVarApplies(pvIdx));
529 
530  // TODO: it may be beneficial to chose this differently.
531  return static_cast<Scalar>(1.0);
532  }
533 
534  static bool eqApplies(unsigned eqIdx)
535  {
536  if (!enableSolvent)
537  return false;
538 
539  return eqIdx == contiSolventEqIdx;
540  }
541 
542  static std::string eqName(unsigned eqIdx OPM_OPTIM_UNUSED)
543  {
544  assert(eqApplies(eqIdx));
545 
546  return "conti^solvent";
547  }
548 
549  static Scalar eqWeight(unsigned eqIdx OPM_OPTIM_UNUSED)
550  {
551  assert(eqApplies(eqIdx));
552 
553  // TODO: it may be beneficial to chose this differently.
554  return static_cast<Scalar>(1.0);
555  }
556 
557  template <class LhsEval>
558  static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
559  const IntensiveQuantities& intQuants)
560  {
561  if (!enableSolvent)
562  return;
563  if (blackoilConserveSurfaceVolume) {
564  storage[contiSolventEqIdx] +=
565  Toolbox::template decay<LhsEval>(intQuants.porosity())
566  * Toolbox::template decay<LhsEval>(intQuants.solventSaturation())
567  * Toolbox::template decay<LhsEval>(intQuants.solventInverseFormationVolumeFactor());
568  } else {
569  storage[contiSolventEqIdx] +=
570  Toolbox::template decay<LhsEval>(intQuants.porosity())
571  * Toolbox::template decay<LhsEval>(intQuants.solventSaturation())
572  * Toolbox::template decay<LhsEval>(intQuants.solventDensity());
573  }
574  }
575 
576  static void computeFlux(RateVector& flux,
577  const ElementContext& elemCtx,
578  unsigned scvfIdx,
579  unsigned timeIdx)
580 
581  {
582  if (!enableSolvent)
583  return;
584 
585  const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
586 
587  unsigned upIdx = extQuants.solventUpstreamIndex();
588  unsigned inIdx = extQuants.interiorIndex();
589  const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
590 
591  if (blackoilConserveSurfaceVolume) {
592  if (upIdx == inIdx)
593  flux[contiSolventEqIdx] =
594  extQuants.solventVolumeFlux()
595  *up.solventInverseFormationVolumeFactor();
596  else
597  flux[contiSolventEqIdx] =
598  extQuants.solventVolumeFlux()
599  *Opm::decay<Scalar>(up.solventInverseFormationVolumeFactor());
600  } else {
601  if (upIdx == inIdx)
602  flux[contiSolventEqIdx] =
603  extQuants.solventVolumeFlux()
604  *up.solventDensity();
605  else
606  flux[contiSolventEqIdx] =
607  extQuants.solventVolumeFlux()
608  *Opm::decay<Scalar>(up.solventDensity());
609  }
610  }
611 
615  static void assignPrimaryVars(PrimaryVariables& priVars,
616  Scalar solventSaturation)
617  {
618  if (!enableSolvent)
619  return;
620 
621  priVars[solventSaturationIdx] = solventSaturation;
622  }
623 
627  static void updatePrimaryVars(PrimaryVariables& newPv,
628  const PrimaryVariables& oldPv,
629  const EqVector& delta)
630  {
631  if (!enableSolvent)
632  return;
633 
634  // do a plain unchopped Newton update
635  newPv[solventSaturationIdx] = oldPv[solventSaturationIdx] - delta[solventSaturationIdx];
636  }
637 
641  static Scalar computeUpdateError(const PrimaryVariables& oldPv OPM_UNUSED,
642  const EqVector& delta OPM_UNUSED)
643  {
644  // do not consider consider the cange of solvent primary variables for
645  // convergence
646  // TODO: maybe this should be changed
647  return static_cast<Scalar>(0.0);
648  }
649 
653  static Scalar computeResidualError(const EqVector& resid)
654  {
655  // do not weight the residual of solvents when it comes to convergence
656  return std::abs(Toolbox::scalarValue(resid[contiSolventEqIdx]));
657  }
658 
659  template <class DofEntity>
660  static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
661  {
662  if (!enableSolvent)
663  return;
664 
665  unsigned dofIdx = model.dofMapper().index(dof);
666 
667  const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
668  outstream << priVars[solventSaturationIdx];
669  }
670 
671  template <class DofEntity>
672  static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
673  {
674  if (!enableSolvent)
675  return;
676 
677  unsigned dofIdx = model.dofMapper().index(dof);
678 
679  PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
680  PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
681 
682  instream >> priVars0[solventSaturationIdx];
683 
684  // set the primary variables for the beginning of the current time step.
685  priVars1 = priVars0[solventSaturationIdx];
686  }
687 
688  static const SolventPvt& solventPvt()
689  { return solventPvt_; }
690 
691  static const TabulatedFunction& ssfnKrg(const ElementContext& elemCtx,
692  unsigned scvIdx,
693  unsigned timeIdx)
694  {
695  unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
696  return ssfnKrg_[satnumRegionIdx];
697  }
698 
699  static const TabulatedFunction& ssfnKrs(const ElementContext& elemCtx,
700  unsigned scvIdx,
701  unsigned timeIdx)
702  {
703  unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
704  return ssfnKrs_[satnumRegionIdx];
705  }
706 
707  static const TabulatedFunction& sof2Krn(const ElementContext& elemCtx,
708  unsigned scvIdx,
709  unsigned timeIdx)
710  {
711  unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
712  return sof2Krn_[satnumRegionIdx];
713  }
714 
715  static const TabulatedFunction& misc(const ElementContext& elemCtx,
716  unsigned scvIdx,
717  unsigned timeIdx)
718  {
719  unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
720  return misc_[miscnumRegionIdx];
721  }
722 
723  static const TabulatedFunction& pmisc(const ElementContext& elemCtx,
724  unsigned scvIdx,
725  unsigned timeIdx)
726  {
727  unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
728  return pmisc_[miscnumRegionIdx];
729  }
730 
731  static const TabulatedFunction& msfnKrsg(const ElementContext& elemCtx,
732  unsigned scvIdx,
733  unsigned timeIdx)
734  {
735  unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
736  return msfnKrsg_[satnumRegionIdx];
737  }
738 
739  static const TabulatedFunction& msfnKro(const ElementContext& elemCtx,
740  unsigned scvIdx,
741  unsigned timeIdx)
742  {
743  unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
744  return msfnKro_[satnumRegionIdx];
745  }
746 
747  static const TabulatedFunction& sorwmis(const ElementContext& elemCtx,
748  unsigned scvIdx,
749  unsigned timeIdx)
750  {
751  unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
752  return sorwmis_[miscnumRegionIdx];
753  }
754 
755  static const TabulatedFunction& sgcwmis(const ElementContext& elemCtx,
756  unsigned scvIdx,
757  unsigned timeIdx)
758  {
759  unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
760  return sgcwmis_[miscnumRegionIdx];
761  }
762 
763  static const TabulatedFunction& tlPMixTable(const ElementContext& elemCtx,
764  unsigned scvIdx,
765  unsigned timeIdx)
766  {
767  unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
768  return tlPMixTable_[miscnumRegionIdx];
769  }
770 
771  static const Scalar& tlMixParamViscosity(const ElementContext& elemCtx,
772  unsigned scvIdx,
773  unsigned timeIdx)
774  {
775  unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
776  return tlMixParamViscosity_[miscnumRegionIdx];
777  }
778 
779  static const Scalar& tlMixParamDensity(const ElementContext& elemCtx,
780  unsigned scvIdx,
781  unsigned timeIdx)
782  {
783  unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
784  return tlMixParamDensity_[miscnumRegionIdx];
785  }
786 
787  static bool isMiscible()
788  {
789  return isMiscible_;
790  }
791 
792 
793 private:
794  static SolventPvt solventPvt_;
795 
796  static std::vector<TabulatedFunction> ssfnKrg_; // the krg(Fs) column of the SSFN table
797  static std::vector<TabulatedFunction> ssfnKrs_; // the krs(Fs) column of the SSFN table
798  static std::vector<TabulatedFunction> sof2Krn_; // the krn(Sn) column of the SOF2 table
799  static std::vector<TabulatedFunction> misc_; // the misc(Ss) column of the MISC table
800  static std::vector<TabulatedFunction> pmisc_; // the pmisc(pg) column of the PMISC table
801  static std::vector<TabulatedFunction> msfnKrsg_; // the krsg(Ssg) column of the MSFN table
802  static std::vector<TabulatedFunction> msfnKro_; // the kro(Ssg) column of the MSFN table
803  static std::vector<TabulatedFunction> sorwmis_; // the sorwmis(Sw) column of the SORWMIS table
804  static std::vector<TabulatedFunction> sgcwmis_; // the sgcwmis(Sw) column of the SGCWMIS table
805 
806  static std::vector<Scalar> tlMixParamViscosity_; // Todd-Longstaff mixing parameter for viscosity
807  static std::vector<Scalar> tlMixParamDensity_; // Todd-Longstaff mixing parameter for density
808  static std::vector<TabulatedFunction> tlPMixTable_; // the tlpmixpa(Po) column of the TLPMIXPA table
809 
810  static bool isMiscible_;
811 };
812 
813 template <class TypeTag, bool enableSolventV>
814 typename BlackOilSolventModule<TypeTag, enableSolventV>::SolventPvt
815 BlackOilSolventModule<TypeTag, enableSolventV>::solventPvt_;
816 
817 template <class TypeTag, bool enableSolventV>
818 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
819 BlackOilSolventModule<TypeTag, enableSolventV>::ssfnKrg_;
820 
821 template <class TypeTag, bool enableSolventV>
822 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
823 BlackOilSolventModule<TypeTag, enableSolventV>::ssfnKrs_;
824 
825 template <class TypeTag, bool enableSolventV>
826 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
827 BlackOilSolventModule<TypeTag, enableSolventV>::sof2Krn_;
828 
829 template <class TypeTag, bool enableSolventV>
830 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
831 BlackOilSolventModule<TypeTag, enableSolventV>::misc_;
832 
833 template <class TypeTag, bool enableSolventV>
834 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
835 BlackOilSolventModule<TypeTag, enableSolventV>::pmisc_;
836 
837 template <class TypeTag, bool enableSolventV>
838 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
839 BlackOilSolventModule<TypeTag, enableSolventV>::msfnKrsg_;
840 
841 template <class TypeTag, bool enableSolventV>
842 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
843 BlackOilSolventModule<TypeTag, enableSolventV>::msfnKro_;
844 
845 template <class TypeTag, bool enableSolventV>
846 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
847 BlackOilSolventModule<TypeTag, enableSolventV>::sorwmis_;
848 
849 template <class TypeTag, bool enableSolventV>
850 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
851 BlackOilSolventModule<TypeTag, enableSolventV>::sgcwmis_;
852 
853 template <class TypeTag, bool enableSolventV>
854 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::Scalar>
855 BlackOilSolventModule<TypeTag, enableSolventV>::tlMixParamViscosity_;
856 
857 template <class TypeTag, bool enableSolventV>
858 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::Scalar>
859 BlackOilSolventModule<TypeTag, enableSolventV>::tlMixParamDensity_;
860 
861 
862 template <class TypeTag, bool enableSolventV>
863 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
864 BlackOilSolventModule<TypeTag, enableSolventV>::tlPMixTable_;
865 
866 template <class TypeTag, bool enableSolventV>
867 bool
868 BlackOilSolventModule<TypeTag, enableSolventV>::isMiscible_;
869 
870 
878 template <class TypeTag, bool enableSolventV = GET_PROP_VALUE(TypeTag, EnableSolvent)>
880 {
881  typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) Implementation;
882 
883  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
884  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
885  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
886  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
887  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
888  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
889  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
890 
892 
893  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
894  static constexpr int solventSaturationIdx = Indices::solventSaturationIdx;
895  static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
896  static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
897  static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
898  static constexpr double cutOff = 1e-16;
899 
900 
901 public:
908  void solventPreSatFuncUpdate_(const ElementContext& elemCtx,
909  unsigned dofIdx,
910  unsigned timeIdx)
911  {
912  const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
913  auto& fs = asImp_().fluidState_;
914  solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx);
915  hydrocarbonSaturation_ = fs.saturation(gasPhaseIdx);
916 
917  // apply a cut-off. Don't waste calculations if no solvent
918  if (std::abs(solventSaturation().value()) < cutOff) {
919  return;
920  }
921 
922  // make the saturation of the gas phase which is used by the saturation functions
923  // the sum of the solvent "saturation" and the saturation the hydrocarbon gas.
924  fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_ + solventSaturation_);
925  }
926 
934  void solventPostSatFuncUpdate_(const ElementContext& elemCtx,
935  unsigned dofIdx,
936  unsigned timeIdx)
937  {
938  // revert the gas "saturation" of the fluid state back to the saturation of the
939  // hydrocarbon gas.
940  auto& fs = asImp_().fluidState_;
941  fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_);
942 
943  solventMobility_ = 0.0;
944 
945  // apply a cut-off. Don't waste calculations if no solvent
946  if (std::abs(solventSaturation().value()) < cutOff) {
947  return;
948  }
949 
950  // Pressure effects on capillary pressure miscibility
951  if(SolventModule::isMiscible()) {
952 
953  const Evaluation& p = fs.pressure(oilPhaseIdx); // or gas pressure?
954  const Evaluation pmisc = SolventModule::pmisc(elemCtx, dofIdx, timeIdx).eval(p, /*extrapolate=*/true);
955  const Evaluation& pgImisc = fs.pressure(gasPhaseIdx);
956 
957  // compute capillary pressure for miscible fluid
958  const auto& problem = elemCtx.problem();
959  const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
960  Evaluation pgMisc = 0.0;
961  Evaluation pC[numPhases];
962  const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
963  MaterialLaw::capillaryPressures(pC, materialParams, fs);
964 
965  //oil is the reference phase for pressure
966  if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
967  pgMisc = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
968  } else {
969  const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
970  pgMisc = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
971  }
972 
973  fs.setPressure(gasPhaseIdx, pmisc * pgMisc + (1.0 - pmisc) * pgImisc);
974  }
975 
976 
977  Evaluation gasSolventSat = hydrocarbonSaturation_ + solventSaturation_;
978 
979  if (std::abs(gasSolventSat.value()) < cutOff) { // avoid division by zero
980  return;
981  }
982 
983  Evaluation Fhydgas = hydrocarbonSaturation_/gasSolventSat;
984  Evaluation Fsolgas = solventSaturation_/gasSolventSat;
985 
986  // account for miscibility of oil and solvent
987  if(SolventModule::isMiscible()) {
988  const auto& misc = SolventModule::misc(elemCtx, dofIdx, timeIdx);
989  const auto& pmisc = SolventModule::pmisc(elemCtx, dofIdx, timeIdx);
990  const Evaluation& p = fs.pressure(oilPhaseIdx); // or gas pressure?
991  const Evaluation miscibility = misc.eval(Fsolgas, /*extrapolate=*/true) * pmisc.eval(p, /*extrapolate=*/true);
992 
993  // TODO adjust endpoints of sn and ssg
994  unsigned cellIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
995  const auto& materialLawManager = elemCtx.problem().materialLawManager();
996  const auto& scaledDrainageInfo =
997  materialLawManager->oilWaterScaledEpsInfoDrainage(cellIdx);
998 
999  const Scalar& sgcr = scaledDrainageInfo.Sgcr;
1000  const Scalar& sogcr = scaledDrainageInfo.Sogcr;
1001  const Evaluation& sw = fs.saturation(waterPhaseIdx);
1002  const auto& sorwmis = SolventModule::sorwmis(elemCtx, dofIdx, timeIdx);
1003  const auto& sgcwmis = SolventModule::sgcwmis(elemCtx, dofIdx, timeIdx);
1004 
1005  Evaluation sor = miscibility * sorwmis.eval(sw, /*extrapolate=*/true) + ( 1.0 - miscibility) * sogcr;
1006  Evaluation sgc = miscibility * sgcwmis.eval(sw, /*extrapolate=*/true) + ( 1.0 - miscibility) * sgcr;
1007 
1008  const Evaluation oilGasSolventSat = gasSolventSat + fs.saturation(oilPhaseIdx);
1009  const Evaluation zero = 0.0;
1010  const Evaluation oilGasSolventEffSat = std::max(oilGasSolventSat - sor - sgc, zero);
1011 
1012  Evaluation F_totalGas = 0.0;
1013  if (std::abs(oilGasSolventEffSat.value()) > cutOff) {
1014  const Evaluation gasSolventEffSat = std::max(gasSolventSat - sgc, zero);
1015  F_totalGas = gasSolventEffSat / oilGasSolventEffSat;
1016  }
1017  const auto& msfnKro = SolventModule::msfnKro(elemCtx, dofIdx, timeIdx);
1018  const auto& msfnKrsg = SolventModule::msfnKrsg(elemCtx, dofIdx, timeIdx);
1019  const auto& sof2Krn = SolventModule::sof2Krn(elemCtx, dofIdx, timeIdx);
1020 
1021  const Evaluation mkrgt = msfnKrsg.eval( F_totalGas, /*extrapolate=*/true ) * sof2Krn.eval( oilGasSolventSat, /*extrapolate=*/true);
1022  const Evaluation mkro = msfnKro.eval( F_totalGas, /*extrapolate=*/true ) * sof2Krn.eval( oilGasSolventSat, /*extrapolate=*/true );
1023 
1024  Evaluation& kro = asImp_().mobility_[oilPhaseIdx];
1025  Evaluation& krg = asImp_().mobility_[gasPhaseIdx];
1026 
1027  // combine immiscible and miscible part of the relperm
1028  krg *= (1.0 - miscibility);
1029  krg += miscibility * mkrgt;
1030  kro *= (1.0 - miscibility);
1031  kro += miscibility * mkro;
1032  }
1033 
1034 
1035 
1036  // compute the mobility of the solvent "phase" and modify the gas phase
1037  const auto& ssfnKrg = SolventModule::ssfnKrg(elemCtx, dofIdx, timeIdx);
1038  const auto& ssfnKrs = SolventModule::ssfnKrs(elemCtx, dofIdx, timeIdx);
1039 
1040  Evaluation& krg = asImp_().mobility_[gasPhaseIdx];
1041  solventMobility_ = krg * ssfnKrs.eval(Fsolgas, /*extrapolate=*/true);
1042  krg *= ssfnKrg.eval(Fhydgas, /*extrapolate=*/true);
1043 
1044  }
1045 
1052  void solventPvtUpdate_(const ElementContext& elemCtx,
1053  unsigned scvIdx,
1054  unsigned timeIdx)
1055  {
1056  const auto& iq = asImp_();
1057  const auto& fs = iq.fluidState();
1058  const auto& solventPvt = SolventModule::solventPvt();
1059 
1060  unsigned pvtRegionIdx = iq.pvtRegionIndex();
1061  solventRefDensity_ = solventPvt.referenceDensity(pvtRegionIdx);
1062  const Evaluation& T = fs.temperature(gasPhaseIdx);
1063  const Evaluation& p = fs.pressure(gasPhaseIdx);
1064  solventInvFormationVolumeFactor_ = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
1065 
1066  solventDensity_ = solventInvFormationVolumeFactor_*solventRefDensity_;
1067  solventViscosity_ = solventPvt.viscosity(pvtRegionIdx, T, p);
1068 
1069  effectiveProperties(elemCtx, scvIdx, timeIdx);
1070 
1071  solventMobility_ /= solventViscosity_;
1072 
1073 
1074  }
1075 
1076  const Evaluation& solventSaturation() const
1077  { return solventSaturation_; }
1078 
1079  const Evaluation& solventDensity() const
1080  { return solventDensity_; }
1081 
1082  const Evaluation& solventViscosity() const
1083  { return solventViscosity_; }
1084 
1085  const Evaluation& solventMobility() const
1086  { return solventMobility_; }
1087 
1088  const Evaluation& solventInverseFormationVolumeFactor() const
1089  { return solventInvFormationVolumeFactor_; }
1090 
1091  // This could be stored pr pvtRegion instead
1092  const Scalar& solventRefDensity() const
1093  { return solventRefDensity_; }
1094 
1095 private:
1096  // Computes the effective properties based on
1097  // Todd-Longstaff mixing model.
1098  void effectiveProperties(const ElementContext& elemCtx,
1099  unsigned scvIdx,
1100  unsigned timeIdx)
1101  {
1102  if (!SolventModule::isMiscible()) {
1103  return;
1104  }
1105 
1106  // Don't waste calculations if no solvent
1107  // Apply a cut-off for small and negative solvent saturations
1108  if ( solventSaturation() < cutOff) { //
1109  return;
1110  }
1111 
1112  auto& fs = asImp_().fluidState_;
1113 
1114  // Compute effective saturations
1115  const auto& sorwmis = SolventModule::sorwmis(elemCtx, scvIdx, timeIdx);
1116  const auto& sgcwmis = SolventModule::sgcwmis(elemCtx, scvIdx, timeIdx);
1117  const Evaluation& sw = fs.saturation(waterPhaseIdx);
1118  const Evaluation oilEffSat = fs.saturation(oilPhaseIdx) - sorwmis.eval(sw, /*extrapolate=*/true);
1119  const Evaluation gasEffSat = fs.saturation(gasPhaseIdx) - sgcwmis.eval(sw, /*extrapolate=*/true);
1120  const Evaluation solventEffSat = solventSaturation() - sgcwmis.eval(sw, /*extrapolate=*/true);
1121  const Evaluation oilGasSolventEffSat = oilEffSat + gasEffSat + solventEffSat;
1122  const Evaluation oilSolventEffSat = oilEffSat + solventEffSat;
1123  const Evaluation solventGasEffSat = solventEffSat + gasEffSat;
1124 
1125  // Compute effective viscosities
1126  const Evaluation& muGas = fs.viscosity(gasPhaseIdx);
1127  const Evaluation& muOil = fs.viscosity(oilPhaseIdx);
1128  const Evaluation& muSolvent = solventViscosity_;
1129 
1130  const Evaluation muOilPow = pow(muOil, 0.25);
1131  const Evaluation muGasPow = pow(muGas, 0.25);
1132  const Evaluation muSolventPow = pow(muSolvent, 0.25);
1133 
1134  Evaluation muMixOilSolvent = muOil;
1135  if ( std::abs(oilSolventEffSat.value()) > cutOff)
1136  muMixOilSolvent *= muSolvent / pow( ( (oilEffSat / oilSolventEffSat) * muSolventPow) + ( (solventEffSat / oilSolventEffSat) * muOilPow) , 4.0);
1137  Evaluation muMixSolventGas = muGas;
1138  if ( std::abs(solventGasEffSat.value()) > cutOff)
1139  muMixSolventGas *= muSolvent / pow( ( (gasEffSat / solventGasEffSat) * muSolventPow) + ( (solventEffSat / solventGasEffSat) * muGasPow) , 4.0);
1140 
1141  Evaluation muMixSolventGasOil = muOil;
1142  if (std::abs(oilGasSolventEffSat.value()) > cutOff)
1143  muMixSolventGasOil *= muSolvent * muGas / pow( ( (oilEffSat / oilGasSolventEffSat) * muSolventPow * muGasPow)
1144  + ( (solventEffSat / oilGasSolventEffSat) * muOilPow * muGasPow) + ( (gasEffSat / oilGasSolventEffSat) * muSolventPow * muOilPow), 4.0);
1145 
1146  // Mixing parameter for viscosity
1147  // The pressureMixingParameter represent the miscibility of the solvent while the mixingParameterViscosity the effect of the porous media.
1148  // The pressureMixingParameter is not implemented in ecl100.
1149  const Evaluation& po = fs.pressure(oilPhaseIdx);
1150  const auto& tlPMixTable = SolventModule::tlPMixTable(elemCtx, scvIdx, timeIdx);
1151  const Evaluation tlMixParamMu = SolventModule::tlMixParamViscosity(elemCtx, scvIdx, timeIdx) * tlPMixTable.eval(po, /*extrapolate=*/true);
1152 
1153  Evaluation muOilEff = pow(muOil,1.0 - tlMixParamMu) * pow(muMixOilSolvent, tlMixParamMu);
1154  Evaluation muGasEff = pow(muGas,1.0 - tlMixParamMu) * pow(muMixSolventGas, tlMixParamMu);
1155  Evaluation muSolventEff = pow(muSolvent,1.0 - tlMixParamMu) * pow(muMixSolventGasOil, tlMixParamMu);
1156 
1157  // Compute effective densities
1158  const Evaluation& rhoGas = fs.density(gasPhaseIdx);
1159  const Evaluation& rhoOil = fs.density(oilPhaseIdx);
1160  const Evaluation& rhoSolvent = solventDensity_;
1161 
1162  // Mixing parameter for density
1163  // The pressureMixingParameter represent the miscibility of the solvent while the mixingParameterDenisty the effect of the porous media.
1164  // The pressureMixingParameter is not implemented in ecl100.
1165  const Evaluation tlMixParamRho = SolventModule::tlMixParamDensity(elemCtx, scvIdx, timeIdx) * tlPMixTable.eval(po, /*extrapolate=*/true);
1166 
1167  // compute effective viscosities for density calculations. These have to
1168  // be recomputed as a different mixing parameter may be used.
1169  const Evaluation muOilEffPow = pow(pow(muOil, 1.0 - tlMixParamRho) * pow(muMixOilSolvent, tlMixParamRho), 0.25);
1170  const Evaluation muGasEffPow = pow(pow(muGas, 1.0 - tlMixParamRho) * pow(muMixSolventGas, tlMixParamRho), 0.25);
1171  const Evaluation muSolventEffPow = pow(pow(muSolvent, 1.0 - tlMixParamRho) * pow(muMixSolventGasOil, tlMixParamRho), 0.25);
1172 
1173  const Evaluation oilGasEffSaturation = oilEffSat + gasEffSat;
1174  Evaluation sof = 0.0;
1175  Evaluation sgf = 0.0;
1176  if ( std::abs(oilGasEffSaturation.value()) > cutOff ) {
1177  sof = oilEffSat / oilGasEffSaturation;
1178  sgf = gasEffSat / oilGasEffSaturation;
1179  }
1180 
1181  const Evaluation muSolventOilGasPow = muSolventPow * ( (sgf * muOilPow) + (sof * muGasPow) );
1182 
1183  Evaluation rhoMixSolventGasOil = 0.0;
1184  if (std::abs(oilGasSolventEffSat.value()) > cutOff )
1185  rhoMixSolventGasOil = (rhoOil * oilEffSat / oilGasSolventEffSat) + (rhoGas * gasEffSat / oilGasSolventEffSat) + (rhoSolvent * solventEffSat / oilGasSolventEffSat);
1186 
1187  Evaluation rhoGasEff = 0.0;
1188  if ( std::abs(muSolvent.value() - muGas.value()) < cutOff ) {
1189  rhoGasEff = ( (1.0 - tlMixParamRho) * rhoGas) + (tlMixParamRho * rhoMixSolventGasOil);
1190  } else {
1191  const Evaluation solventGasEffFraction = (muGasPow * (muSolventPow - muGasEffPow)) / (muGasEffPow * (muSolventPow - muGasPow));
1192  rhoGasEff = (rhoGas * solventGasEffFraction) + (rhoSolvent * (1.0 - solventGasEffFraction));
1193  }
1194 
1195  Evaluation rhoOilEff = 0.0;
1196  if ( std::abs(muGas.value() - muOil.value()) < cutOff ) {
1197  rhoOilEff = ( (1.0 - tlMixParamRho) * rhoOil) + (tlMixParamRho * rhoMixSolventGasOil);
1198  } else {
1199  const Evaluation solventOilEffFraction = (muOilPow * (muOilEffPow - muSolventPow)) / (muOilEffPow * (muOilPow - muSolventPow));
1200  rhoOilEff = (rhoOil * solventOilEffFraction) + (rhoSolvent * (1.0 - solventOilEffFraction));
1201  }
1202 
1203  Evaluation rhoSolventEff = 0.0;
1204  if ( std::abs( ( muSolventOilGasPow.value() - (muOilPow.value() * muGasPow.value()))) < cutOff ) {
1205  rhoSolventEff = ( (1.0 - tlMixParamRho) * rhoSolvent) + (tlMixParamRho * rhoMixSolventGasOil);
1206  } else {
1207  const Evaluation sfraction_se = (muSolventOilGasPow - ( muOilPow * muGasPow * muSolventPow / muSolventEffPow) ) / ( muSolventOilGasPow - (muOilPow * muGasPow));
1208  rhoSolventEff = (rhoSolvent * sfraction_se) + (rhoGas * sgf * (1.0 - sfraction_se)) + (rhoOil * sof * (1.0 - sfraction_se));
1209  }
1210 
1211  unsigned pvtRegionIdx = asImp_().pvtRegionIndex();
1212  // compute invB from densities.
1213  const Evaluation bOilEff = rhoOilEff / (FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) + FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) * fs.Rs());
1214  const Evaluation bGasEff = rhoGasEff / (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) + FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) * fs.Rv());
1215  const Evaluation bSolventEff = rhoSolventEff / solventRefDensity();
1216 
1217  // account for pressure effects
1218  const auto& pmiscTable = SolventModule::pmisc(elemCtx, scvIdx, timeIdx);
1219  const Evaluation pmisc = pmiscTable.eval(po, /*extrapolate=*/true);
1220 
1221  // copy the unmodified invB factors
1222  const Evaluation bo = fs.invB(oilPhaseIdx);
1223  const Evaluation bg = fs.invB(gasPhaseIdx);
1224  const Evaluation bs = solventInverseFormationVolumeFactor();
1225 
1226  // Set the effective invB factors
1227  fs.setInvB(oilPhaseIdx, pmisc * bOilEff + (1.0 - pmisc) * bo);
1228  fs.setInvB(gasPhaseIdx, pmisc * bGasEff + (1.0 - pmisc) * bg);
1229  solventInvFormationVolumeFactor_ = pmisc * bSolventEff + (1.0 - pmisc) * bs;
1230 
1231  // set the densities
1232  fs.setDensity(oilPhaseIdx, fs.invB(oilPhaseIdx) * (FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) + FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) * fs.Rs()) );
1233  fs.setDensity(gasPhaseIdx, fs.invB(gasPhaseIdx) * (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) + FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) * fs.Rv()) );
1234  solventDensity_ = solventInverseFormationVolumeFactor() * solventRefDensity();
1235 
1236  // set the viscosity / mobility
1237  // TODO make it possible to store and modify the viscosity in fs directly
1238 
1239  // keep the mu*b interpolation
1240  Evaluation& mobo = asImp_().mobility_[oilPhaseIdx];
1241  muOilEff = fs.invB(oilPhaseIdx) / (pmisc * bOilEff / muOilEff + (1.0 - pmisc) * bo / muOil);
1242  mobo *= muOil / muOilEff;
1243 
1244  Evaluation& mobg = asImp_().mobility_[gasPhaseIdx];
1245  muGasEff = fs.invB(gasPhaseIdx) / (pmisc * bGasEff / muGasEff + (1.0 - pmisc) * bg / muGas);
1246  mobg *= muGas / muGasEff;
1247 
1248  // Update viscosity of solvent
1249  solventViscosity_ = solventInvFormationVolumeFactor_ / (pmisc * bSolventEff / muSolventEff + (1.0 - pmisc) * bs / muSolvent);
1250  }
1251 
1252 protected:
1253  Implementation& asImp_()
1254  { return *static_cast<Implementation*>(this); }
1255 
1256  Evaluation hydrocarbonSaturation_;
1257  Evaluation solventSaturation_;
1258  Evaluation solventDensity_;
1259  Evaluation solventViscosity_;
1260  Evaluation solventMobility_;
1261  Evaluation solventInvFormationVolumeFactor_;
1262 
1263  Scalar solventRefDensity_;
1264 };
1265 
1266 template <class TypeTag>
1268 {
1269  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
1270  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
1271  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
1272 
1273 
1274 public:
1275  void solventPreSatFuncUpdate_(const ElementContext& elemCtx OPM_UNUSED,
1276  unsigned scvIdx OPM_UNUSED,
1277  unsigned timeIdx OPM_UNUSED)
1278  { }
1279 
1280  void solventPostSatFuncUpdate_(const ElementContext& elemCtx OPM_UNUSED,
1281  unsigned scvIdx OPM_UNUSED,
1282  unsigned timeIdx OPM_UNUSED)
1283  { }
1284 
1285  void solventPvtUpdate_(const ElementContext& elemCtx OPM_UNUSED,
1286  unsigned scvIdx OPM_UNUSED,
1287  unsigned timeIdx OPM_UNUSED)
1288  { }
1289 
1290  const Evaluation& solventSaturation() const
1291  { OPM_THROW(std::runtime_error, "solventSaturation() called but solvents are disabled"); }
1292 
1293  const Evaluation& solventDensity() const
1294  { OPM_THROW(std::runtime_error, "solventDensity() called but solvents are disabled"); }
1295 
1296  const Evaluation& solventViscosity() const
1297  { OPM_THROW(std::runtime_error, "solventViscosity() called but solvents are disabled"); }
1298 
1299  const Evaluation& solventMobility() const
1300  { OPM_THROW(std::runtime_error, "solventMobility() called but solvents are disabled"); }
1301 
1302  const Evaluation& solventInverseFormationVolumeFactor() const
1303  { OPM_THROW(std::runtime_error, "solventInverseFormationVolumeFactor() called but solvents are disabled"); }
1304 
1305  const Scalar& solventRefDensity() const
1306  { OPM_THROW(std::runtime_error, "solventRefDensity() called but solvents are disabled"); }
1307 };
1308 
1316 template <class TypeTag, bool enableSolventV = GET_PROP_VALUE(TypeTag, EnableSolvent)>
1318 {
1319  typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
1320 
1321  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
1322  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
1323  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
1324  typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
1325  typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
1326  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
1327  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
1328 
1329  typedef Opm::MathToolbox<Evaluation> Toolbox;
1330 
1331  static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
1332  static constexpr int dimWorld = GridView::dimensionworld;
1333 
1334  typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
1335  typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
1336 
1337 public:
1342  template <class Dummy = bool> // we need to make this method a template to avoid
1343  // compiler errors if it is not instantiated!
1344  void updateVolumeFluxPerm(const ElementContext& elemCtx,
1345  unsigned scvfIdx,
1346  unsigned timeIdx)
1347  {
1348  const auto& gradCalc = elemCtx.gradientCalculator();
1349  Ewoms::PressureCallback<TypeTag> pressureCallback(elemCtx);
1350 
1351  const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
1352  const auto& faceNormal = scvf.normal();
1353 
1354  unsigned i = scvf.interiorIndex();
1355  unsigned j = scvf.exteriorIndex();
1356 
1357  // calculate the "raw" pressure gradient
1358  DimEvalVector solventPGrad;
1359  pressureCallback.setPhaseIndex(gasPhaseIdx);
1360  gradCalc.calculateGradient(solventPGrad,
1361  elemCtx,
1362  scvfIdx,
1363  pressureCallback);
1364  Opm::Valgrind::CheckDefined(solventPGrad);
1365 
1366  // correct the pressure gradients by the gravitational acceleration
1367  if (EWOMS_GET_PARAM(TypeTag, bool, EnableGravity)) {
1368  // estimate the gravitational acceleration at a given SCV face
1369  // using the arithmetic mean
1370  const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
1371  const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx);
1372 
1373  const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
1374  const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
1375 
1376  const auto& posIn = elemCtx.pos(i, timeIdx);
1377  const auto& posEx = elemCtx.pos(j, timeIdx);
1378  const auto& posFace = scvf.integrationPos();
1379 
1380  // the distance between the centers of the control volumes
1381  DimVector distVecIn(posIn);
1382  DimVector distVecEx(posEx);
1383  DimVector distVecTotal(posEx);
1384 
1385  distVecIn -= posFace;
1386  distVecEx -= posFace;
1387  distVecTotal -= posIn;
1388  Scalar absDistTotalSquared = distVecTotal.two_norm2();
1389 
1390  // calculate the hydrostatic pressure at the integration point of the face
1391  auto rhoIn = intQuantsIn.solventDensity();
1392  auto pStatIn = - rhoIn*(gIn*distVecIn);
1393 
1394  // the quantities on the exterior side of the face do not influence the
1395  // result for the TPFA scheme, so they can be treated as scalar values.
1396  Scalar rhoEx = Toolbox::value(intQuantsEx.solventDensity());
1397  Scalar pStatEx = - rhoEx*(gEx*distVecEx);
1398 
1399  // compute the hydrostatic gradient between the two control volumes (this
1400  // gradient exhibitis the same direction as the vector between the two
1401  // control volume centers and the length (pStaticExterior -
1402  // pStaticInterior)/distanceInteriorToExterior
1403  DimEvalVector f(distVecTotal);
1404  f *= (pStatEx - pStatIn)/absDistTotalSquared;
1405 
1406  // calculate the final potential gradient
1407  for (unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
1408  solventPGrad[dimIdx] += f[dimIdx];
1409 
1410  if (!Opm::isfinite(solventPGrad[dimIdx])) {
1411  OPM_THROW(Opm::NumericalProblem,
1412  "Non-finite potential gradient for solvent 'phase'");
1413  }
1414  }
1415  }
1416 
1417  // determine the upstream and downstream DOFs
1418  Evaluation solventPGradNormal = 0.0;
1419  for (unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
1420  solventPGradNormal += solventPGrad[dimIdx]*faceNormal[dimIdx];
1421 
1422  if (solventPGradNormal > 0) {
1423  solventUpstreamDofIdx_ = j;
1424  solventDownstreamDofIdx_ = i;
1425  }
1426  else {
1427  solventUpstreamDofIdx_ = i;
1428  solventDownstreamDofIdx_ = j;
1429  }
1430 
1431  const auto& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
1432 
1433  // this is also slightly hacky because it assumes that the derivative of the
1434  // flux between two DOFs only depends on the primary variables in the
1435  // upstream direction. For non-TPFA flux approximation schemes, this is not
1436  // true...
1437  if (solventUpstreamDofIdx_ == i)
1438  solventVolumeFlux_ = solventPGradNormal*up.solventMobility();
1439  else
1440  solventVolumeFlux_ = solventPGradNormal*Opm::scalarValue(up.solventMobility());
1441  }
1442 
1447  template <class Dummy = bool> // we need to make this method a template to avoid
1448  // compiler errors if it is not instantiated!
1449  void updateVolumeFluxTrans(const ElementContext& elemCtx,
1450  unsigned scvfIdx,
1451  unsigned timeIdx)
1452  {
1453  const ExtensiveQuantities& extQuants = asImp_();
1454 
1455  unsigned interiorDofIdx = extQuants.interiorIndex();
1456  unsigned exteriorDofIdx = extQuants.exteriorIndex();
1457  assert(interiorDofIdx != exteriorDofIdx);
1458 
1459  const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
1460  const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
1461 
1462  unsigned I = elemCtx.globalSpaceIndex(interiorDofIdx, timeIdx);
1463  unsigned J = elemCtx.globalSpaceIndex(exteriorDofIdx, timeIdx);
1464 
1465  Scalar thpres = elemCtx.problem().thresholdPressure(I, J);
1466  Scalar trans = elemCtx.problem().transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
1467  Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
1468 
1469  Scalar zIn = elemCtx.problem().dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
1470  Scalar zEx = elemCtx.problem().dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
1471  Scalar distZ = zIn - zEx;
1472 
1473  const Evaluation& rhoIn = intQuantsIn.solventDensity();
1474  Scalar rhoEx = Toolbox::value(intQuantsEx.solventDensity());
1475  const Evaluation& rhoAvg = rhoIn*0.5 + rhoEx*0.5;
1476 
1477  const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(gasPhaseIdx);
1478  Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(gasPhaseIdx));
1479  pressureExterior += distZ*g*rhoAvg;
1480 
1481  Evaluation pressureDiffSolvent = pressureExterior - pressureInterior;
1482  if (std::abs(Opm::scalarValue(pressureDiffSolvent)) > thpres) {
1483  if (pressureDiffSolvent < 0.0)
1484  pressureDiffSolvent += thpres;
1485  else
1486  pressureDiffSolvent -= thpres;
1487  }
1488  else
1489  pressureDiffSolvent = 0.0;
1490 
1491  if (pressureDiffSolvent > 0.0) {
1492  solventUpstreamDofIdx_ = exteriorDofIdx;
1493  solventDownstreamDofIdx_ = interiorDofIdx;
1494  }
1495  else if (pressureDiffSolvent < 0.0) {
1496  solventUpstreamDofIdx_ = interiorDofIdx;
1497  solventDownstreamDofIdx_ = exteriorDofIdx;
1498  }
1499  else {
1500  // pressure potential gradient is zero; force consistent upstream and
1501  // downstream indices over the intersection regardless of the side which it
1502  // is looked at.
1503  solventUpstreamDofIdx_ = std::min(interiorDofIdx, exteriorDofIdx);
1504  solventDownstreamDofIdx_ = std::max(interiorDofIdx, exteriorDofIdx);
1505  solventVolumeFlux_ = 0.0;
1506  return;
1507  }
1508 
1509  Scalar faceArea = elemCtx.stencil(timeIdx).interiorFace(scvfIdx).area();
1510  const IntensiveQuantities& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
1511  if (solventUpstreamDofIdx_ == interiorDofIdx)
1512  solventVolumeFlux_ =
1513  up.solventMobility()
1514  *(-trans/faceArea)
1515  *pressureDiffSolvent;
1516  else
1517  solventVolumeFlux_ =
1518  Opm::scalarValue(up.solventMobility())
1519  *(-trans/faceArea)
1520  *pressureDiffSolvent;
1521  }
1522 
1523  unsigned solventUpstreamIndex() const
1524  { return solventUpstreamDofIdx_; }
1525 
1526  unsigned solventDownstreamIndex() const
1527  { return solventDownstreamDofIdx_; }
1528 
1529  const Evaluation& solventVolumeFlux() const
1530  { return solventVolumeFlux_; }
1531 
1532 private:
1533  Implementation& asImp_()
1534  { return *static_cast<Implementation*>(this); }
1535 
1536  Evaluation solventVolumeFlux_;
1537  unsigned solventUpstreamDofIdx_;
1538  unsigned solventDownstreamDofIdx_;
1539 };
1540 
1541 template <class TypeTag>
1543 {
1544  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
1545  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
1546 
1547 public:
1548  void updateVolumeFluxPerm(const ElementContext& elemCtx OPM_UNUSED,
1549  unsigned scvfIdx OPM_UNUSED,
1550  unsigned timeIdx OPM_UNUSED)
1551  { }
1552 
1553  void updateVolumeFluxTrans(const ElementContext& elemCtx OPM_UNUSED,
1554  unsigned scvfIdx OPM_UNUSED,
1555  unsigned timeIdx OPM_UNUSED)
1556  { }
1557 
1558  unsigned solventUpstreamIndex() const
1559  { OPM_THROW(std::runtime_error, "solventUpstreamIndex() called but solvents are disabled"); }
1560 
1561  unsigned solventDownstreamIndex() const
1562  { OPM_THROW(std::runtime_error, "solventDownstreamIndex() called but solvents are disabled"); }
1563 
1564  const Evaluation& solventVolumeFlux() const
1565  { OPM_THROW(std::runtime_error, "solventVolumeFlux() called but solvents are disabled"); }
1566 };
1567 
1568 } // namespace Ewoms
1569 
1570 #endif
static Scalar computeResidualError(const EqVector &resid)
Return how much a residual is considered an error.
Definition: blackoilsolventmodules.hh:653
static void registerParameters()
Register all run-time parameters for the multi-phase VTK output module.
Definition: vtkblackoilsolventmodule.hh:97
void updateVolumeFluxTrans(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Method which calculates the volume flux of the polymer "phase" using the gas pressure potential diffe...
Definition: blackoilsolventmodules.hh:1449
Provides the volumetric quantities required for the equations needed by the solvents extension of the...
Definition: blackoilsolventmodules.hh:879
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 solventPvtUpdate_(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Update the intensive PVT properties needed to handle solvents from the primary variables.
Definition: blackoilsolventmodules.hh:1052
static void registerOutputModules(Model &model, Simulator &simulator)
Register all solvent specific VTK and ECL output modules.
Definition: blackoilsolventmodules.hh:500
static void setSgcmis(unsigned miscRegionIdx, const TabulatedFunction &sgcwmis)
Misicibe critical gas saturation function wrt water saturation of a single region.
Definition: blackoilsolventmodules.hh:445
static void setSof2(unsigned satRegionIdx, const TabulatedFunction &sof2Krn)
Specify misicible hydrocabon relative permeability wrt water of a single region.
Definition: blackoilsolventmodules.hh:388
static void setTlmixpar(unsigned miscRegionIdx, const Scalar &tlMixParamViscosity, const Scalar &tlMixParamDensity)
Todd-Longstaff mixing parameters of a single region.
Definition: blackoilsolventmodules.hh:456
static void setSolventPvt(const SolventPvt &value)
Specify the solvent PVT of a all PVT regions.
Definition: blackoilsolventmodules.hh:478
VTK output module for the black oil model&#39;s solvent related quantities.
Definition: vtkblackoilsolventmodule.hh:71
static void registerParameters()
Register all run-time parameters for the black-oil solvent module.
Definition: blackoilsolventmodules.hh:488
Declares the properties required by the black oil model.
static void setMisc(unsigned miscRegionIdx, const TabulatedFunction &misc)
Misicibility function wrt solvent fraction of a single region.
Definition: blackoilsolventmodules.hh:399
static void setSsfn(unsigned satRegionIdx, const TabulatedFunction &ssfnKrg, const TabulatedFunction &ssfnKrs)
Specify the solvent saturation functions of a single region.
Definition: blackoilsolventmodules.hh:375
static void updatePrimaryVars(PrimaryVariables &newPv, const PrimaryVariables &oldPv, const EqVector &delta)
Do a Newton-Raphson update the primary variables of the solvents.
Definition: blackoilsolventmodules.hh:627
Provides the solvent specific extensive quantities to the generic black-oil module&#39;s extensive quanti...
Definition: blackoilsolventmodules.hh:1317
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which the pressure should be returned.
Definition: quantitycallbacks.hh:108
static void setSorwmis(unsigned miscRegionIdx, const TabulatedFunction &sorwmis)
Misicibe residual oil saturation function wrt water saturation of a single region.
Definition: blackoilsolventmodules.hh:434
This method contains all callback classes for quantities that are required by some extensive quantiti...
Contains the high level supplements required to extend the black oil model by solvents.
Definition: blackoilprimaryvariables.hh:47
Callback class for a phase pressure.
Definition: quantitycallbacks.hh:83
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
void updateVolumeFluxPerm(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Method which calculates the volume flux of the polymer "phase" using the pressure potential gradient ...
Definition: blackoilsolventmodules.hh:1344
static void setPmisc(unsigned miscRegionIdx, const TabulatedFunction &pmisc)
Misicibility function wrt pressure of a single region.
Definition: blackoilsolventmodules.hh:410
static Scalar computeUpdateError(const PrimaryVariables &oldPv OPM_UNUSED, const EqVector &delta OPM_UNUSED)
Return how much a Newton-Raphson update is considered an error.
Definition: blackoilsolventmodules.hh:641
static void setMsfn(unsigned satRegionIdx, const TabulatedFunction &msfnKrsg, const TabulatedFunction &msfnKro)
Specify misicible relative permeability multipliers of a single region.
Definition: blackoilsolventmodules.hh:421
static void setTlpmixpa(unsigned miscRegionIdx, const TabulatedFunction &tlPMixTable)
Todd-Longstaff mixing parameter multiplier wrt pressure of a single region.
Definition: blackoilsolventmodules.hh:469
static void setNumSatRegions(unsigned numRegions)
Specify the number of satuation regions.
Definition: blackoilsolventmodules.hh:364
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
static void assignPrimaryVars(PrimaryVariables &priVars, Scalar solventSaturation)
Assign the solvent specific primary variables to a PrimaryVariables object.
Definition: blackoilsolventmodules.hh:615
void solventPreSatFuncUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Called before the saturation functions are doing their magic.
Definition: blackoilsolventmodules.hh:908
VTK output module for the black oil model&#39;s solvent related quantities.
void solventPostSatFuncUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Called after the saturation functions have been doing their magic.
Definition: blackoilsolventmodules.hh:934