blackoilpolymermodules.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_POLYMER_MODULE_HH
29 #define EWOMS_BLACK_OIL_POLYMER_MODULE_HH
30 
31 #include "blackoilproperties.hh"
34 
35 #include <opm/material/common/Tabulated1DFunction.hpp>
36 
37 #if HAVE_OPM_PARSER
38 #include <opm/parser/eclipse/Deck/Deck.hpp>
39 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
40 #include <opm/parser/eclipse/EclipseState/Tables/PlyadsTable.hpp>
41 #include <opm/parser/eclipse/EclipseState/Tables/PlymaxTable.hpp>
42 #include <opm/parser/eclipse/EclipseState/Tables/PlyrockTable.hpp>
43 #include <opm/parser/eclipse/EclipseState/Tables/PlyshlogTable.hpp>
44 #include <opm/parser/eclipse/EclipseState/Tables/PlyviscTable.hpp>
45 #include <opm/parser/eclipse/EclipseState/Tables/PlyshlogTable.hpp>
46 #endif
47 
48 #include <opm/common/Valgrind.hpp>
49 #include <opm/common/Unused.hpp>
50 #include <opm/common/ErrorMacros.hpp>
51 #include <opm/common/Exceptions.hpp>
52 
53 #include <dune/common/fvector.hh>
54 
55 #include <string>
56 
57 namespace Ewoms {
63 template <class TypeTag, bool enablePolymerV = GET_PROP_VALUE(TypeTag, EnablePolymer)>
65 {
66  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
67  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
68  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
69  typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
70  typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
71  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
72  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
73  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
74  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
75  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
76  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
77  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
78 
79  typedef Opm::MathToolbox<Evaluation> Toolbox;
80 
81  typedef typename Opm::Tabulated1DFunction<Scalar> TabulatedFunction;
82 
83  static constexpr unsigned polymerConcentrationIdx = Indices::polymerConcentrationIdx;
84  static constexpr unsigned contiPolymerEqIdx = Indices::contiPolymerEqIdx;
85  static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
86 
87 
88  static constexpr unsigned enablePolymer = enablePolymerV;
89  static constexpr unsigned numEq = GET_PROP_VALUE(TypeTag, NumEq);
90  static constexpr unsigned numPhases = FluidSystem::numPhases;
91 
92 public:
93  enum AdsorptionBehaviour { Desorption = 1, NoDesorption = 2 };
94 
95 #if HAVE_OPM_PARSER
96 
99  static void initFromDeck(const Opm::Deck& deck, const Opm::EclipseState& eclState)
100  {
101  // some sanity checks: if polymers are enabled, the POLYMER keyword must be
102  // present, if polymers are disabled the keyword must not be present.
103  if (enablePolymer && !deck.hasKeyword("POLYMER")) {
104  OPM_THROW(std::runtime_error,
105  "Non-trivial polymer treatment requested at compile time, but "
106  "the deck does not contain the POLYMER keyword");
107  }
108  else if (!enablePolymer && deck.hasKeyword("POLYMER")) {
109  OPM_THROW(std::runtime_error,
110  "Polymer treatment disabled at compile time, but the deck "
111  "contains the POLYMER keyword");
112  }
113 
114  if (!deck.hasKeyword("POLYMER"))
115  return; // polymer treatment is supposed to be disabled
116 
117  const auto& tableManager = eclState.getTableManager();
118 
119  unsigned numSatRegions = tableManager.getTabdims().getNumSatTables();
120  setNumSatRegions(numSatRegions);
121 
122  // initialize the objects which deal with the PLYROCK keyword
123  const auto& plyrockTables = tableManager.getPlyrockTables();
124  if (!plyrockTables.empty()) {
125  assert(numSatRegions == plyrockTables.size());
126  for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) {
127  const auto& plyrockTable = plyrockTables.template getTable<Opm::PlyrockTable>(satRegionIdx);
128  setPlyrock(satRegionIdx,
129  plyrockTable.getDeadPoreVolumeColumn()[satRegionIdx],
130  plyrockTable.getResidualResistanceFactorColumn()[satRegionIdx],
131  plyrockTable.getRockDensityFactorColumn()[satRegionIdx],
132  static_cast<AdsorptionBehaviour>(plyrockTable.getAdsorbtionIndexColumn()[satRegionIdx]),
133  plyrockTable.getMaxAdsorbtionColumn()[satRegionIdx]);
134  }
135  } else {
136  OPM_THROW(std::runtime_error, "PLYROCK must be specified in POLYMER runs\n");
137  }
138 
139  // initialize the objects which deal with the PLYADS keyword
140  const auto& plyadsTables = tableManager.getPlyadsTables();
141  if (!plyadsTables.empty()) {
142  assert(numSatRegions == plyadsTables.size());
143  for (unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) {
144  const auto& plyadsTable = plyadsTables.template getTable<Opm::PlyadsTable>(satRegionIdx);
145  // Copy data
146  const auto& c = plyadsTable.getPolymerConcentrationColumn();
147  const auto& ads = plyadsTable.getAdsorbedPolymerColumn();
148  plyadsAdsorbedPolymer_[satRegionIdx].setXYContainers(c, ads);
149  }
150  } else {
151  OPM_THROW(std::runtime_error, "PLYADS must be specified in POLYMER runs\n");
152  }
153 
154 
155  unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables();
156  setNumPvtRegions(numPvtRegions);
157 
158  // initialize the objects which deal with the PLYVISC keyword
159  const auto& plyviscTables = tableManager.getPlyviscTables();
160  if (!plyviscTables.empty()) {
161 
162  assert(numPvtRegions == plyviscTables.size());
163  for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
164  const auto& plyadsTable = plyviscTables.template getTable<Opm::PlyviscTable>(pvtRegionIdx);
165  // Copy data
166  const auto& c = plyadsTable.getPolymerConcentrationColumn();
167  const auto& visc = plyadsTable.getViscosityMultiplierColumn();
168  plyviscViscosityMultiplierTable_[pvtRegionIdx].setXYContainers(c, visc);
169  }
170  } else {
171  OPM_THROW(std::runtime_error, "PLYVISC must be specified in POLYMER runs\n");
172  }
173 
174  // initialize the objects which deal with the PLYMAX keyword
175  const auto& plymaxTables = tableManager.getPlymaxTables();
176  unsigned numMixRegions = plymaxTables.size();
177  setNumMixRegions(numMixRegions);
178  if (!plymaxTables.empty()) {
179  for (unsigned mixRegionIdx = 0; mixRegionIdx < numMixRegions; ++ mixRegionIdx) {
180  const auto& plymaxTable = plymaxTables.template getTable<Opm::PlymaxTable>(mixRegionIdx);
181  setPlymax(mixRegionIdx, plymaxTable.getPolymerConcentrationColumn()[mixRegionIdx]);
182  }
183  } else {
184  OPM_THROW(std::runtime_error, "PLYMAX must be specified in POLYMER runs\n");
185  }
186 
187  if (deck.hasKeyword("PLMIXPAR")) {
188  // initialize the objects which deal with the PLMIXPAR keyword
189  for (unsigned mixRegionIdx = 0; mixRegionIdx < numMixRegions; ++ mixRegionIdx) {
190  const auto& plmixparRecord = deck.getKeyword("PLMIXPAR").getRecord(mixRegionIdx);
191  setPlmixpar(mixRegionIdx, plmixparRecord.getItem("TODD_LONGSTAFF").getSIDouble(0));
192  }
193  } else {
194  OPM_THROW(std::runtime_error, "PLMIXPAR must be specified in POLYMER runs\n");
195  }
196 
197  hasPlyshlog_ = deck.hasKeyword("PLYSHLOG");
198  hasShrate_ = deck.hasKeyword("SHRATE");
199 
200  if (hasPlyshlog_) {
201  const auto& plyshlogTables = tableManager.getPlyshlogTables();
202  assert(numPvtRegions == plyshlogTables.size());
203  plyshlogShearEffectRefMultiplier_.resize(numPvtRegions);
204  plyshlogShearEffectRefLogVelocity_.resize(numPvtRegions);
205  for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
206  const auto& plyshlogTable = plyshlogTables.template getTable<Opm::PlyshlogTable>(pvtRegionIdx);
207 
208  Scalar plyshlogRefPolymerConcentration = plyshlogTable.getRefPolymerConcentration();
209  std::vector<Scalar> waterVelocity = plyshlogTable.getWaterVelocityColumn().vectorCopy();
210  std::vector<Scalar> shearMultiplier = plyshlogTable.getShearMultiplierColumn().vectorCopy();
211 
212  // do the unit version here for the waterVelocity
213  Opm::UnitSystem unitSystem = deck.getActiveUnitSystem();
214  double siFactor = hasShrate_? unitSystem.parse("1/Time").getSIScaling() : unitSystem.parse("Length/Time").getSIScaling();
215  for (size_t i = 0; i < waterVelocity.size(); ++i ) {
216  waterVelocity[i] *= siFactor;
217  // for plyshlog the input must be stored as logarithms
218  // the interpolation is then done the log-space.
219  waterVelocity[i] = std::log(waterVelocity[i]);
220  }
221 
222  Scalar refViscMult = plyviscViscosityMultiplierTable_[pvtRegionIdx].eval(plyshlogRefPolymerConcentration, /*extrapolate=*/true);
223  // convert the table using referece conditions
224  for (size_t i = 0; i < waterVelocity.size(); ++i ) {
225  shearMultiplier[i] *= refViscMult;
226  shearMultiplier[i] -= 1;
227  shearMultiplier[i] /= (refViscMult - 1);
228  shearMultiplier[i] = shearMultiplier[i];
229  }
230  plyshlogShearEffectRefMultiplier_[pvtRegionIdx].resize(waterVelocity.size());
231  plyshlogShearEffectRefLogVelocity_[pvtRegionIdx].resize(waterVelocity.size());
232 
233  for (size_t i = 0; i < waterVelocity.size(); ++i ) {
234  plyshlogShearEffectRefMultiplier_[pvtRegionIdx][i] = shearMultiplier[i];
235  plyshlogShearEffectRefLogVelocity_[pvtRegionIdx][i] = waterVelocity[i];
236  }
237  }
238  }
239 
240  if (hasShrate_) {
241  if(!hasPlyshlog_) {
242  OPM_THROW(std::runtime_error, "PLYSHLOG must be specified if SHRATE is used in POLYMER runs\n");
243  }
244  const auto& shrateKeyword = deck.getKeyword("SHRATE");
245  const std::vector<double>& shrateFromDeck = shrateKeyword.getSIDoubleData();
246  shrate_.resize(numPvtRegions);
247  for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
248  if (shrateFromDeck.empty()) {
249  shrate_[pvtRegionIdx] = 4.8; //default;
250  } else if (shrateFromDeck.size() == numPvtRegions) {
251  shrate_[pvtRegionIdx] = shrateKeyword.getSIDoubleData()[pvtRegionIdx];
252  } else {
253  OPM_THROW(std::runtime_error, "SHRATE must either have 0 or number of NUMPVT entries\n");
254  }
255  }
256  }
257 
258  }
259 #endif
260 
266  static void setNumSatRegions(unsigned numRegions)
267  {
268  plyrockDeadPoreVolume_.resize(numRegions);
269  plyrockResidualResistanceFactor_.resize(numRegions);
270  plyrockRockDensityFactor_.resize(numRegions);
271  plyrockAdsorbtionIndex_.resize(numRegions);
272  plyrockMaxAdsorbtion_.resize(numRegions);
273  plyadsAdsorbedPolymer_.resize(numRegions);
274  }
275 
281  static void setPlyrock(unsigned satRegionIdx,
282  const Scalar& plyrockDeadPoreVolume,
283  const Scalar& plyrockResidualResistanceFactor,
284  const Scalar& plyrockRockDensityFactor,
285  const Scalar& plyrockAdsorbtionIndex,
286  const Scalar& plyrockMaxAdsorbtion)
287  {
288  plyrockDeadPoreVolume_[satRegionIdx] = plyrockDeadPoreVolume;
289  plyrockResidualResistanceFactor_[satRegionIdx] = plyrockResidualResistanceFactor;
290  plyrockRockDensityFactor_[satRegionIdx] = plyrockRockDensityFactor;
291  plyrockAdsorbtionIndex_[satRegionIdx] = plyrockAdsorbtionIndex;
292  plyrockMaxAdsorbtion_[satRegionIdx] = plyrockMaxAdsorbtion;
293  }
294 
300  static void setAdsrock(unsigned satRegionIdx,
301  const TabulatedFunction& plyadsAdsorbedPolymer)
302  {
303  plyadsAdsorbedPolymer_[satRegionIdx] = plyadsAdsorbedPolymer;
304  }
305 
311  static void setNumPvtRegions(unsigned numRegions)
312  {
313  plyviscViscosityMultiplierTable_.resize(numRegions);
314  }
315 
321  static void setPlyvisc(unsigned satRegionIdx,
322  const TabulatedFunction& plyviscViscosityMultiplierTable)
323  {
324  plyviscViscosityMultiplierTable_[satRegionIdx] = plyviscViscosityMultiplierTable;
325  }
326 
332  static void setNumMixRegions(unsigned numRegions)
333  {
334  plymaxMaxConcentration_.resize(numRegions);
335  plymixparToddLongstaff_.resize(numRegions);
336  }
337 
343  static void setPlymax(unsigned mixRegionIdx,
344  const Scalar& plymaxMaxConcentration)
345  {
346  plymaxMaxConcentration_[mixRegionIdx] = plymaxMaxConcentration;
347  }
348 
354  static void setPlmixpar(unsigned mixRegionIdx,
355  const Scalar& plymixparToddLongstaff)
356  {
357  plymixparToddLongstaff_[mixRegionIdx] = plymixparToddLongstaff;
358  }
359 
360 
361 
365  static void registerParameters()
366  {
367  if (!enablePolymer)
368  // polymers have been disabled at compile time
369  return;
370 
372  }
373 
377  static void registerOutputModules(Model& model,
378  Simulator& simulator)
379  {
380  if (!enablePolymer)
381  // polymers have been disabled at compile time
382  return;
383 
384  model.addOutputModule(new Ewoms::VtkBlackOilPolymerModule<TypeTag>(simulator));
385  }
386 
387  static bool primaryVarApplies(unsigned pvIdx)
388  {
389  if (!enablePolymer)
390  // polymers have been disabled at compile time
391  return false;
392 
393  return pvIdx == polymerConcentrationIdx;
394  }
395 
396  static std::string primaryVarName(unsigned pvIdx OPM_OPTIM_UNUSED)
397  {
398  assert(primaryVarApplies(pvIdx));
399 
400  return "polymer_waterconcentration";
401  }
402 
403  static Scalar primaryVarWeight(unsigned pvIdx OPM_OPTIM_UNUSED)
404  {
405  assert(primaryVarApplies(pvIdx));
406 
407  // TODO: it may be beneficial to chose this differently.
408  return static_cast<Scalar>(1.0);
409  }
410 
411  static bool eqApplies(unsigned eqIdx)
412  {
413  if (!enablePolymer)
414  return false;
415 
416  return eqIdx == contiPolymerEqIdx;
417  }
418 
419  static std::string eqName(unsigned eqIdx OPM_OPTIM_UNUSED)
420  {
421  assert(eqApplies(eqIdx));
422 
423  return "conti^polymer";
424  }
425 
426  static Scalar eqWeight(unsigned eqIdx OPM_OPTIM_UNUSED)
427  {
428  assert(eqApplies(eqIdx));
429 
430  // TODO: it may be beneficial to chose this differently.
431  return static_cast<Scalar>(1.0);
432  }
433 
434  // must be called after water storage is computed
435  template <class LhsEval>
436  static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
437  const IntensiveQuantities& intQuants)
438  {
439  if (!enablePolymer)
440  return;
441 
442  const auto& fs = intQuants.fluidState();
443 
444  LhsEval surfaceVolumeWater =
445  Toolbox::template decay<LhsEval>(fs.saturation(waterPhaseIdx))
446  * Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx))
447  * Toolbox::template decay<LhsEval>(intQuants.porosity());
448 
449  // avoid singular matrix if no water is present.
450  surfaceVolumeWater = Opm::max(surfaceVolumeWater, 1e-10);
451 
452  // polymer in water phase
453  storage[contiPolymerEqIdx] += surfaceVolumeWater
454  * Toolbox::template decay<LhsEval>(intQuants.polymerConcentration())
455  * (1.0 - Toolbox::template decay<LhsEval>(intQuants.polymerDeadPoreVolume()));
456 
457  // polymer in solid phase
458  storage[contiPolymerEqIdx] +=
459  Toolbox::template decay<LhsEval>(1.0 - intQuants.porosity())
460  * Toolbox::template decay<LhsEval>(intQuants.polymerRockDensity())
461  * Toolbox::template decay<LhsEval>(intQuants.polymerAdsorption());
462 
463  }
464 
465  static void computeFlux(RateVector& flux,
466  const ElementContext& elemCtx,
467  unsigned scvfIdx,
468  unsigned timeIdx)
469 
470  {
471  if (!enablePolymer)
472  return;
473 
474  const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
475 
476  unsigned upIdx = extQuants.upstreamIndex(FluidSystem::waterPhaseIdx);
477  unsigned inIdx = extQuants.interiorIndex();
478  const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
479  const unsigned contiWaterEqIdx = Indices::conti0EqIdx + Indices::canonicalToActiveComponentIndex(FluidSystem::waterCompIdx);
480 
481 
482  if (upIdx == inIdx) {
483  flux[contiPolymerEqIdx] =
484  extQuants.volumeFlux(waterPhaseIdx)
485  *up.fluidState().invB(waterPhaseIdx)
486  *up.polymerViscosityCorrection()
487  /extQuants.polymerShearFactor()
488  *up.polymerConcentration();
489 
490  // modify water
491  flux[contiWaterEqIdx] /=
492  extQuants.waterShearFactor();
493  } else {
494  flux[contiPolymerEqIdx] =
495  extQuants.volumeFlux(waterPhaseIdx)
496  *Opm::decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
497  *Opm::decay<Scalar>(up.polymerViscosityCorrection())
498  /Opm::decay<Scalar>(extQuants.polymerShearFactor())
499  *Opm::decay<Scalar>(up.polymerConcentration());
500 
501  // modify water
502  flux[contiWaterEqIdx] /=
503  Opm::decay<Scalar>(extQuants.waterShearFactor());
504  }
505 
506  }
507 
511  static void assignPrimaryVars(PrimaryVariables& priVars,
512  Scalar polymerConcentration)
513  {
514  if (!enablePolymer)
515  return;
516 
517  priVars[polymerConcentrationIdx] = polymerConcentration;
518  }
519 
523  static void updatePrimaryVars(PrimaryVariables& newPv,
524  const PrimaryVariables& oldPv,
525  const EqVector& delta)
526  {
527  if (!enablePolymer)
528  return;
529 
530  // do a plain unchopped Newton update
531  newPv[polymerConcentrationIdx] = oldPv[polymerConcentrationIdx] - delta[polymerConcentrationIdx];
532  }
533 
537  static Scalar computeUpdateError(const PrimaryVariables& oldPv OPM_UNUSED,
538  const EqVector& delta OPM_UNUSED)
539  {
540  // do not consider consider the cange of polymer primary variables for
541  // convergence
542  // TODO: maybe this should be changed
543  return static_cast<Scalar>(0.0);
544  }
545 
549  static Scalar computeResidualError(const EqVector& resid)
550  {
551  // do not weight the residual of polymers when it comes to convergence
552  return std::abs(Toolbox::scalarValue(resid[contiPolymerEqIdx]));
553  }
554 
555  template <class DofEntity>
556  static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
557  {
558  if (!enablePolymer)
559  return;
560 
561  unsigned dofIdx = model.dofMapper().index(dof);
562  const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
563  outstream << priVars[polymerConcentrationIdx];
564  }
565 
566  template <class DofEntity>
567  static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
568  {
569  if (!enablePolymer)
570  return;
571 
572  unsigned dofIdx = model.dofMapper().index(dof);
573  PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
574  PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
575 
576  instream >> priVars0[polymerConcentrationIdx];
577 
578  // set the primary variables for the beginning of the current time step.
579  priVars1 = priVars0[polymerConcentrationIdx];
580  }
581 
582  static const Scalar plyrockDeadPoreVolume(const ElementContext& elemCtx,
583  unsigned scvIdx,
584  unsigned timeIdx)
585  {
586  unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
587  return plyrockDeadPoreVolume_[satnumRegionIdx];
588  }
589 
590  static const Scalar plyrockResidualResistanceFactor(const ElementContext& elemCtx,
591  unsigned scvIdx,
592  unsigned timeIdx)
593  {
594  unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
595  return plyrockResidualResistanceFactor_[satnumRegionIdx];
596  }
597 
598  static const Scalar plyrockRockDensityFactor(const ElementContext& elemCtx,
599  unsigned scvIdx,
600  unsigned timeIdx)
601  {
602  unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
603  return plyrockRockDensityFactor_[satnumRegionIdx];
604  }
605 
606  static const Scalar plyrockAdsorbtionIndex(const ElementContext& elemCtx,
607  unsigned scvIdx,
608  unsigned timeIdx)
609  {
610  unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
611  return plyrockAdsorbtionIndex_[satnumRegionIdx];
612  }
613 
614  static const Scalar plyrockMaxAdsorbtion(const ElementContext& elemCtx,
615  unsigned scvIdx,
616  unsigned timeIdx)
617  {
618  unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
619  return plyrockMaxAdsorbtion_[satnumRegionIdx];
620  }
621 
622  static const TabulatedFunction& plyadsAdsorbedPolymer(const ElementContext& elemCtx,
623  unsigned scvIdx,
624  unsigned timeIdx)
625  {
626  unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
627  return plyadsAdsorbedPolymer_[satnumRegionIdx];
628  }
629 
630  static const TabulatedFunction& plyviscViscosityMultiplierTable(const ElementContext& elemCtx,
631  unsigned scvIdx,
632  unsigned timeIdx)
633  {
634  unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
635  return plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
636  }
637 
638  static const TabulatedFunction& plyviscViscosityMultiplierTable(unsigned pvtnumRegionIdx)
639  {
640  return plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
641  }
642 
643  static const Scalar plymaxMaxConcentration(const ElementContext& elemCtx,
644  unsigned scvIdx,
645  unsigned timeIdx)
646  {
647  unsigned polymerMixRegionIdx = elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
648  return plymaxMaxConcentration_[polymerMixRegionIdx];
649  }
650 
651  static const Scalar plymixparToddLongstaff(const ElementContext& elemCtx,
652  unsigned scvIdx,
653  unsigned timeIdx)
654  {
655  unsigned polymerMixRegionIdx = elemCtx.problem().plmixnumRegionIndex(elemCtx, scvIdx, timeIdx);
656  return plymixparToddLongstaff_[polymerMixRegionIdx];
657  }
658 
659  static bool hasPlyshlog()
660  {
661  return hasPlyshlog_;
662  }
663 
664  static bool hasShrate()
665  {
666  return hasShrate_;
667  }
668 
669  static const Scalar shrate(unsigned pvtnumRegionIdx)
670  {
671  return shrate_[pvtnumRegionIdx];
672  }
673 
680  template <class Evaluation>
681  static Evaluation computeShearFactor(const Evaluation& polymerConcentration,
682  unsigned pvtnumRegionIdx,
683  const Evaluation& v0) {
684 
685  const auto& viscosityMultiplierTable = plyviscViscosityMultiplierTable_[pvtnumRegionIdx];
686  Scalar viscosityMultiplier = viscosityMultiplierTable.eval( Opm::scalarValue(polymerConcentration), /*extrapolate=*/true);
687 
688  const Scalar eps = 1e-14;
689  // return 1.0 if the polymer has no effect on the water.
690  if ( std::abs( (viscosityMultiplier - 1.0) ) < eps){
691  return 1.0;
692  }
693 
694  const std::vector<Scalar>& shearEffectRefLogVelocity = plyshlogShearEffectRefLogVelocity_[pvtnumRegionIdx];
695  auto v0AbsLog = Opm::log(Opm::abs(v0));
696  // return 1.0 if the velocity /sharte is smaller than the first velocity entry.
697  if (v0AbsLog < shearEffectRefLogVelocity[0])
698  return 1.0;
699 
700  // compute shear factor from input
701  // Z = (1 + (P - 1) * M(v) ) / P
702  // where M(v) is computed from user input
703  // and P = viscosityMultiplier
704  const std::vector<Scalar>& shearEffectRefMultiplier = plyshlogShearEffectRefMultiplier_[pvtnumRegionIdx];
705  size_t numTableEntries = shearEffectRefLogVelocity.size();
706  assert(shearEffectRefMultiplier.size() == numTableEntries);
707 
708  std::vector<Scalar> shearEffectMultiplier(numTableEntries, 1.0);
709  for (size_t i = 0; i < numTableEntries; ++i) {
710  shearEffectMultiplier[i] = (1.0 + (viscosityMultiplier - 1.0)*shearEffectRefMultiplier[i]) / viscosityMultiplier;
711  shearEffectMultiplier[i] = Opm::log(shearEffectMultiplier[i]);
712  }
713  // store the logarithmic velocity and logarithmic multipliers in a table for easy look up and
714  // linear interpolation in the logarithmic space.
715  TabulatedFunction logShearEffectMultiplier = TabulatedFunction(numTableEntries, shearEffectRefLogVelocity, shearEffectMultiplier, /*bool sortInputs =*/ false );
716 
717  // Find sheared velocity (v) that satisfies
718  // F = log(v) + log (Z) - log(v0) = 0;
719 
720  // Set up the function
721  // u = log(v)
722  auto F = [&logShearEffectMultiplier, &v0AbsLog](const Evaluation& u) {
723  return u + logShearEffectMultiplier.eval(u, true) - v0AbsLog;
724  };
725  // and its derivative
726  auto dF = [&logShearEffectMultiplier](const Evaluation& u) {
727  return 1 + logShearEffectMultiplier.evalDerivative(u, true);
728  };
729 
730  // Solve F = 0 using Newton
731  // Use log(v0) as initial value for u
732  auto u = v0AbsLog;
733  bool converged = false;
734  for (int i = 0; i < 20; ++i ) {
735  auto f = F(u);
736  auto df = dF(u);
737  u -= f/df;
738  if (std::abs(Opm::scalarValue(f)) < 1e-12) {
739  converged = true;
740  break;
741  }
742  }
743  if (!converged) {
744  OPM_THROW(std::runtime_error, "Not able to compute shear velocity. \n");
745  }
746 
747  // return the shear factor
748  return Opm::exp( logShearEffectMultiplier.eval(u, /*extrapolate=*/true) );
749 
750  }
751 
752 
753 
754 
755 private:
756  static std::vector<Scalar> plyrockDeadPoreVolume_;
757  static std::vector<Scalar> plyrockResidualResistanceFactor_;
758  static std::vector<Scalar> plyrockRockDensityFactor_;
759  static std::vector<Scalar> plyrockAdsorbtionIndex_;
760  static std::vector<Scalar> plyrockMaxAdsorbtion_;
761  static std::vector<TabulatedFunction> plyadsAdsorbedPolymer_;
762  static std::vector<TabulatedFunction> plyviscViscosityMultiplierTable_;
763  static std::vector<Scalar> plymaxMaxConcentration_;
764  static std::vector<Scalar> plymixparToddLongstaff_;
765  static std::vector<std::vector<Scalar>> plyshlogShearEffectRefMultiplier_;
766  static std::vector<std::vector<Scalar>> plyshlogShearEffectRefLogVelocity_;
767  static std::vector<Scalar> shrate_;
768  static bool hasShrate_;
769  static bool hasPlyshlog_;
770 
771 };
772 
773 
774 
775 template <class TypeTag, bool enablePolymerV>
776 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
777 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyrockDeadPoreVolume_;
778 
779 template <class TypeTag, bool enablePolymerV>
780 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
781 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyrockResidualResistanceFactor_;
782 
783 template <class TypeTag, bool enablePolymerV>
784 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
785 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyrockRockDensityFactor_;
786 
787 template <class TypeTag, bool enablePolymerV>
788 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
789 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyrockAdsorbtionIndex_;
790 
791 template <class TypeTag, bool enablePolymerV>
792 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
793 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyrockMaxAdsorbtion_;
794 
795 template <class TypeTag, bool enablePolymerV>
796 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::TabulatedFunction>
797 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyadsAdsorbedPolymer_;
798 
799 template <class TypeTag, bool enablePolymerV>
800 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::TabulatedFunction>
801 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyviscViscosityMultiplierTable_;
802 
803 template <class TypeTag, bool enablePolymerV>
804 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
805 BlackOilPolymerModule<TypeTag, enablePolymerV>::plymaxMaxConcentration_;
806 
807 template <class TypeTag, bool enablePolymerV>
808 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
809 BlackOilPolymerModule<TypeTag, enablePolymerV>::plymixparToddLongstaff_;
810 
811 template <class TypeTag, bool enablePolymerV>
812 std::vector<std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>>
813 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyshlogShearEffectRefMultiplier_;
814 
815 template <class TypeTag, bool enablePolymerV>
816 std::vector<std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>>
817 BlackOilPolymerModule<TypeTag, enablePolymerV>::plyshlogShearEffectRefLogVelocity_;
818 
819 template <class TypeTag, bool enablePolymerV>
820 std::vector<typename BlackOilPolymerModule<TypeTag, enablePolymerV>::Scalar>
821 BlackOilPolymerModule<TypeTag, enablePolymerV>::shrate_;
822 
823 template <class TypeTag, bool enablePolymerV>
824 bool
825 BlackOilPolymerModule<TypeTag, enablePolymerV>::hasShrate_;
826 
827 template <class TypeTag, bool enablePolymerV>
828 bool
829 BlackOilPolymerModule<TypeTag, enablePolymerV>::hasPlyshlog_;
830 
838 template <class TypeTag, bool enablePolymerV = GET_PROP_VALUE(TypeTag, EnablePolymer)>
840 {
841  typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) Implementation;
842 
843  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
844  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
845  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
846  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
847  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
848  typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
849  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
850 
852 
853  enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
854  static constexpr int polymerConcentrationIdx = Indices::polymerConcentrationIdx;
855  static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
856 
857 
858 public:
859 
865  void polymerPropertiesUpdate_(const ElementContext& elemCtx,
866  unsigned dofIdx,
867  unsigned timeIdx)
868  {
869  const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
870  polymerConcentration_ = priVars.makeEvaluation(polymerConcentrationIdx, timeIdx);
871  const Scalar cmax = PolymerModule::plymaxMaxConcentration(elemCtx, dofIdx, timeIdx);
872 
873  // permeability reduction due to polymer
874  const Scalar& maxAdsorbtion = PolymerModule::plyrockMaxAdsorbtion(elemCtx, dofIdx, timeIdx);
875  const auto& plyadsAdsorbedPolymer = PolymerModule::plyadsAdsorbedPolymer(elemCtx, dofIdx, timeIdx);
876  polymerAdsorption_ = plyadsAdsorbedPolymer.eval(polymerConcentration_, /*extrapolate=*/true);
877  if (PolymerModule::plyrockAdsorbtionIndex(elemCtx, dofIdx, timeIdx) == PolymerModule::NoDesorption ) {
878  const Scalar& maxPolymerAdsorption = elemCtx.problem().maxPolymerAdsorption(elemCtx, dofIdx, timeIdx);
879  polymerAdsorption_ = std::max(Evaluation(maxPolymerAdsorption) , polymerAdsorption_);
880  }
881 
882  // compute resitanceFactor
883  const Scalar& residualResistanceFactor = PolymerModule::plyrockResidualResistanceFactor(elemCtx, dofIdx, timeIdx);
884  const Evaluation resistanceFactor = 1.0 + (residualResistanceFactor - 1.0) * polymerAdsorption_ / maxAdsorbtion;
885 
886  // compute effective viscosities
887  auto& fs = asImp_().fluidState_;
888  const Evaluation& muWater = fs.viscosity(waterPhaseIdx);
889  const auto& viscosityMultiplier = PolymerModule::plyviscViscosityMultiplierTable(elemCtx, dofIdx, timeIdx);
890  Evaluation viscosityMixture = viscosityMultiplier.eval(polymerConcentration_, /*extrapolate=*/true) * muWater;
891 
892  // Do the Todd-Longstaff mixing
893  const Scalar plymixparToddLongstaff = PolymerModule::plymixparToddLongstaff(elemCtx, dofIdx, timeIdx);
894  Evaluation viscosityPolymer = viscosityMultiplier.eval(cmax, /*extrapolate=*/true) * muWater;
895  Evaluation viscosityPolymerEffective = pow(viscosityMixture, plymixparToddLongstaff) * pow(viscosityPolymer, 1.0 - plymixparToddLongstaff);
896  Evaluation viscosityWaterEffective = pow(viscosityMixture, plymixparToddLongstaff) * pow(muWater, 1.0 - plymixparToddLongstaff);
897 
898  Evaluation cbar = polymerConcentration_ / cmax;
899  // waterViscosity / effectiveWaterViscosity
900  waterViscosityCorrection_ = muWater * ( (1.0 - cbar) / viscosityWaterEffective + cbar / viscosityPolymerEffective );
901 
902  // adjust water mobility
903  asImp_().mobility_[waterPhaseIdx] *= waterViscosityCorrection_ / resistanceFactor;
904 
905  // effectiveWaterViscosity / effectivePolymerViscosity
906  polymerViscosityCorrection_ = (muWater / waterViscosityCorrection_) / viscosityPolymerEffective;
907 
908  // update rock properties
909  polymerDeadPoreVolume_ = PolymerModule::plyrockDeadPoreVolume(elemCtx, dofIdx, timeIdx);
910  polymerRockDensity_ = PolymerModule::plyrockRockDensityFactor(elemCtx, dofIdx, timeIdx);
911  }
912 
913  const Evaluation& polymerConcentration() const
914  { return polymerConcentration_; }
915 
916  const Scalar& polymerDeadPoreVolume() const
917  { return polymerDeadPoreVolume_; }
918 
919  const Evaluation& polymerAdsorption() const
920  { return polymerAdsorption_; }
921 
922  const Scalar& polymerRockDensity() const
923  { return polymerRockDensity_; }
924 
925  // effectiveWaterViscosity / effectivePolymerViscosity
926  const Evaluation& polymerViscosityCorrection() const
927  { return polymerViscosityCorrection_; }
928 
929  // waterViscosity / effectiveWaterViscosity
930  const Evaluation& waterViscosityCorrection() const
931  { return waterViscosityCorrection_; }
932 
933 
934 protected:
935  Implementation& asImp_()
936  { return *static_cast<Implementation*>(this); }
937 
938  Evaluation polymerConcentration_;
939  Scalar polymerDeadPoreVolume_;
940  Scalar polymerRockDensity_;
941  Evaluation polymerAdsorption_;
942  Evaluation polymerViscosityCorrection_;
943  Evaluation waterViscosityCorrection_;
944 
945 
946 };
947 
948 template <class TypeTag>
950 {
951  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
952  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
953  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
954 
955 public:
956  void polymerPropertiesUpdate_(const ElementContext& elemCtx OPM_UNUSED,
957  unsigned scvIdx OPM_UNUSED,
958  unsigned timeIdx OPM_UNUSED)
959  { }
960 
961  const Evaluation& polymerConcentration() const
962  { OPM_THROW(std::runtime_error, "polymerConcentration() called but polymers are disabled"); }
963 
964  const Evaluation& polymerDeadPoreVolume() const
965  { OPM_THROW(std::runtime_error, "polymerDeadPoreVolume() called but polymers are disabled"); }
966 
967  const Evaluation& polymerAdsorption() const
968  { OPM_THROW(std::runtime_error, "polymerAdsorption() called but polymers are disabled"); }
969 
970  const Evaluation& polymerRockDensity() const
971  { OPM_THROW(std::runtime_error, "polymerRockDensity() called but polymers are disabled"); }
972 
973  const Evaluation& polymerViscosityCorrection() const
974  { OPM_THROW(std::runtime_error, "polymerViscosityCorrection() called but polymers are disabled"); }
975 
976  const Evaluation& waterViscosityCorrection() const
977  { OPM_THROW(std::runtime_error, "waterViscosityCorrection() called but polymers are disabled"); }
978 
979 
980 };
981 
982 
990 template <class TypeTag, bool enablePolymerV = GET_PROP_VALUE(TypeTag, EnablePolymer)>
992 {
993  typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
994 
995  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
996  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
997  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
998  typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
999  typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
1000  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
1001  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
1002 
1003  typedef Opm::MathToolbox<Evaluation> Toolbox;
1004 
1006 
1007 
1008  static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
1009  static constexpr int dimWorld = GridView::dimensionworld;
1010  static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
1011 
1012  typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
1013  typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
1014 
1015 public:
1022  template <class Dummy = bool> // we need to make this method a template to avoid
1023  // compiler errors if it is not instantiated!
1024  void updateShearMultipliersPerm(const ElementContext& elemCtx OPM_UNUSED,
1025  unsigned scvfIdx OPM_UNUSED,
1026  unsigned timeIdx OPM_UNUSED)
1027  {
1028  OPM_THROW(std::runtime_error,
1029  "The extension of the blackoil model for polymers is not yet "
1030  "implemented for problems specified using permeabilities.");
1031  }
1032 
1039  template <class Dummy = bool> // we need to make this method a template to avoid
1040  // compiler errors if it is not instantiated!
1041  void updateShearMultipliers(const ElementContext& elemCtx,
1042  unsigned scvfIdx,
1043  unsigned timeIdx)
1044  {
1045 
1046  waterShearFactor_ = 1.0;
1047  polymerShearFactor_ = 1.0;
1048 
1049  if (!PolymerModule::hasPlyshlog())
1050  return;
1051 
1052  const ExtensiveQuantities& extQuants = asImp_();
1053  unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
1054  unsigned interiorDofIdx = extQuants.interiorIndex();
1055  unsigned exteriorDofIdx = extQuants.exteriorIndex();
1056  const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
1057  const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
1058  const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
1059 
1060  // compute water velocity from flux
1061  Evaluation poroAvg = intQuantsIn.porosity()*0.5 + intQuantsEx.porosity()*0.5;
1062  unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvfIdx, timeIdx);
1063  const Evaluation& Sw = up.fluidState().saturation(waterPhaseIdx);
1064  unsigned cellIdx = elemCtx.globalSpaceIndex(scvfIdx, timeIdx);
1065  const auto& materialLawManager = elemCtx.problem().materialLawManager();
1066  const auto& scaledDrainageInfo =
1067  materialLawManager->oilWaterScaledEpsInfoDrainage(cellIdx);
1068  const Scalar& Swcr = scaledDrainageInfo.Swcr;
1069 
1070  // guard against zero porosity and no mobile water
1071  Evaluation denom = Opm::max(poroAvg * (Sw - Swcr), 1e-12);
1072  Evaluation waterVolumeVelocity = extQuants.volumeFlux(waterPhaseIdx) / denom;
1073 
1074  // if shrate is specified. Compute shrate based on the water velocity
1075  if (PolymerModule::hasShrate()) {
1076  const Evaluation& relWater = up.relativePermeability(waterPhaseIdx);
1077  Scalar trans = elemCtx.problem().transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
1078  if (trans > 0.0) {
1079  Scalar faceArea = elemCtx.stencil(timeIdx).interiorFace(scvfIdx).area();
1080  auto dist = elemCtx.pos(interiorDofIdx, timeIdx ) - elemCtx.pos(exteriorDofIdx, timeIdx );
1081  // compute permeability from transmissibility.
1082  Scalar absPerm = trans / faceArea * dist.two_norm();
1083  waterVolumeVelocity *= PolymerModule::shrate(pvtnumRegionIdx) * Opm::sqrt( poroAvg * Sw / (relWater * absPerm ) ) ;
1084  assert(Opm::isfinite(waterVolumeVelocity));
1085  }
1086  }
1087 
1088  // compute share factors for water and polymer
1089  waterShearFactor_ = PolymerModule::computeShearFactor(up.polymerConcentration(), pvtnumRegionIdx, waterVolumeVelocity);
1090  polymerShearFactor_ = PolymerModule::computeShearFactor(up.polymerConcentration(), pvtnumRegionIdx, waterVolumeVelocity * up.polymerViscosityCorrection());
1091 
1092  }
1093 
1094  const Evaluation& polymerShearFactor() const
1095  { return polymerShearFactor_; }
1096 
1097  const Evaluation& waterShearFactor() const
1098  { return waterShearFactor_; }
1099 
1100 
1101 private:
1102  Implementation& asImp_()
1103  { return *static_cast<Implementation*>(this); }
1104 
1105  Evaluation polymerShearFactor_;
1106  Evaluation waterShearFactor_;
1107 
1108 };
1109 
1110 template <class TypeTag>
1112 {
1113  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
1114  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
1115 
1116 public:
1117  void updateShearMultipliers(const ElementContext& elemCtx OPM_UNUSED,
1118  unsigned scvfIdx OPM_UNUSED,
1119  unsigned timeIdx OPM_UNUSED)
1120  { }
1121 
1122  void updateShearMultipliersPerm(const ElementContext& elemCtx OPM_UNUSED,
1123  unsigned scvfIdx OPM_UNUSED,
1124  unsigned timeIdx OPM_UNUSED)
1125  { }
1126 
1127  const Evaluation& polymerShearFactor() const
1128  { OPM_THROW(std::runtime_error, "polymerShearFactor() called but polymers are disabled"); }
1129 
1130  const Evaluation& waterShearFactor() const
1131  { OPM_THROW(std::runtime_error, "waterShearFactor() called but polymers are disabled"); }
1132 };
1133 
1134 
1135 } // namespace Ewoms
1136 
1137 #endif
Provides the volumetric quantities required for the equations needed by the polymers extension of the...
Definition: blackoilpolymermodules.hh:839
static Scalar computeResidualError(const EqVector &resid)
Return how much a residual is considered an error.
Definition: blackoilpolymermodules.hh:549
static void setAdsrock(unsigned satRegionIdx, const TabulatedFunction &plyadsAdsorbedPolymer)
Specify the polymer rock properties a single region.
Definition: blackoilpolymermodules.hh:300
Definition: baseauxiliarymodule.hh:37
VTK output module for the black oil model&#39;s polymer related quantities.
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: blackoilpolymermodules.hh:537
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
void updateShearMultipliersPerm(const ElementContext &elemCtx OPM_UNUSED, unsigned scvfIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
Method which calculates the shear factor based on flow velocity.
Definition: blackoilpolymermodules.hh:1024
static void assignPrimaryVars(PrimaryVariables &priVars, Scalar polymerConcentration)
Assign the polymer specific primary variables to a PrimaryVariables object.
Definition: blackoilpolymermodules.hh:511
Provides the polymer specific extensive quantities to the generic black-oil module&#39;s extensive quanti...
Definition: blackoilpolymermodules.hh:991
static void setPlmixpar(unsigned mixRegionIdx, const Scalar &plymixparToddLongstaff)
Specify the maximum polymer concentration a single region.
Definition: blackoilpolymermodules.hh:354
void polymerPropertiesUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the intensive properties needed to handle polymers from the primary variables.
Definition: blackoilpolymermodules.hh:865
static void registerOutputModules(Model &model, Simulator &simulator)
Register all polymer specific VTK and ECL output modules.
Definition: blackoilpolymermodules.hh:377
Declares the properties required by the black oil model.
static void setNumSatRegions(unsigned numRegions)
Specify the number of satuation regions.
Definition: blackoilpolymermodules.hh:266
static Evaluation computeShearFactor(const Evaluation &polymerConcentration, unsigned pvtnumRegionIdx, const Evaluation &v0)
Computes the shear factor.
Definition: blackoilpolymermodules.hh:681
void updateShearMultipliers(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Method which calculates the shear factor based on flow velocity.
Definition: blackoilpolymermodules.hh:1041
static void updatePrimaryVars(PrimaryVariables &newPv, const PrimaryVariables &oldPv, const EqVector &delta)
Do a Newton-Raphson update the primary variables of the polymers.
Definition: blackoilpolymermodules.hh:523
This method contains all callback classes for quantities that are required by some extensive quantiti...
static void setNumPvtRegions(unsigned numRegions)
Specify the number of pvt regions.
Definition: blackoilpolymermodules.hh:311
static void setPlyvisc(unsigned satRegionIdx, const TabulatedFunction &plyviscViscosityMultiplierTable)
Specify the polymer viscosity a single region.
Definition: blackoilpolymermodules.hh:321
static void registerParameters()
Register all run-time parameters for the black-oil polymer module.
Definition: blackoilpolymermodules.hh:365
static void setPlyrock(unsigned satRegionIdx, const Scalar &plyrockDeadPoreVolume, const Scalar &plyrockResidualResistanceFactor, const Scalar &plyrockRockDensityFactor, const Scalar &plyrockAdsorbtionIndex, const Scalar &plyrockMaxAdsorbtion)
Specify the polymer rock properties a single region.
Definition: blackoilpolymermodules.hh:281
static void setNumMixRegions(unsigned numRegions)
Specify the number of mix regions.
Definition: blackoilpolymermodules.hh:332
Contains the high level supplements required to extend the black oil model by polymer.
Definition: blackoilpolymermodules.hh:64
static void registerParameters()
Register all run-time parameters for the multi-phase VTK output module.
Definition: vtkblackoilpolymermodule.hh:101
static void setPlymax(unsigned mixRegionIdx, const Scalar &plymaxMaxConcentration)
Specify the maximum polymer concentration a single region.
Definition: blackoilpolymermodules.hh:343
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
VTK output module for the black oil model&#39;s polymer related quantities.
Definition: vtkblackoilpolymermodule.hh:75