My Project
reservoirproblem.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_RESERVOIR_PROBLEM_HH
29 #define EWOMS_RESERVOIR_PROBLEM_HH
30 
31 #include <opm/models/blackoil/blackoilproperties.hh>
32 
33 #include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
34 #include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
35 #include <opm/material/fluidstates/CompositionalFluidState.hpp>
36 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
37 #include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
38 #include <opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp>
39 #include <opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp>
40 #include <opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp>
41 
42 #include <dune/grid/yaspgrid.hh>
43 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
44 
45 #include <dune/common/version.hh>
46 #include <dune/common/fvector.hh>
47 #include <dune/common/fmatrix.hh>
48 
49 #include <vector>
50 #include <string>
51 
52 namespace Opm {
53 template <class TypeTag>
54 class ReservoirProblem;
55 
56 } // namespace Opm
57 
58 namespace Opm::Properties {
59 
60 
61 namespace TTag {
62 
64 
65 } // namespace TTag
66 
67 // Maximum depth of the reservoir
68 template<class TypeTag, class MyTypeTag>
69 struct MaxDepth { using type = UndefinedProperty; };
70 // The temperature inside the reservoir
71 template<class TypeTag, class MyTypeTag>
72 struct Temperature { using type = UndefinedProperty; };
73 // The width of producer/injector wells as a fraction of the width of the spatial domain
74 template<class TypeTag, class MyTypeTag>
75 struct WellWidth { using type = UndefinedProperty; };
76 
77 // Set the grid type
78 template<class TypeTag>
79 struct Grid<TypeTag, TTag::ReservoirBaseProblem> { using type = Dune::YaspGrid<2>; };
80 
81 // Set the problem property
82 template<class TypeTag>
83 struct Problem<TypeTag, TTag::ReservoirBaseProblem> { using type = Opm::ReservoirProblem<TypeTag>; };
84 
85 // Set the material Law
86 template<class TypeTag>
87 struct MaterialLaw<TypeTag, TTag::ReservoirBaseProblem>
88 {
89 private:
90  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
91  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
92 
93  using Traits = Opm::
94  ThreePhaseMaterialTraits<Scalar,
95  /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
96  /*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
97  /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx>;
98 
99 public:
100  using type = Opm::LinearMaterial<Traits>;
101 };
102 
103 // Write the Newton convergence behavior to disk?
104 template<class TypeTag>
105 struct NewtonWriteConvergence<TypeTag, TTag::ReservoirBaseProblem> { static constexpr bool value = false; };
106 
107 // Enable gravity
108 template<class TypeTag>
109 struct EnableGravity<TypeTag, TTag::ReservoirBaseProblem> { static constexpr bool value = true; };
110 
111 // Enable constraint DOFs?
112 template<class TypeTag>
113 struct EnableConstraints<TypeTag, TTag::ReservoirBaseProblem> { static constexpr bool value = true; };
114 
115 // set the defaults for some problem specific properties
116 template<class TypeTag>
117 struct MaxDepth<TypeTag, TTag::ReservoirBaseProblem>
118 {
119  using type = GetPropType<TypeTag, Scalar>;
120  static constexpr type value = 2500;
121 };
122 template<class TypeTag>
123 struct Temperature<TypeTag, TTag::ReservoirBaseProblem>
124 {
125  using type = GetPropType<TypeTag, Scalar>;
126  static constexpr type value = 293.15;
127 };
128 
133 template<class TypeTag>
134 struct EndTime<TypeTag, TTag::ReservoirBaseProblem>
135 {
136  using type = GetPropType<TypeTag, Scalar>;
137  static constexpr type value = 1000.0*24*60*60;
138 };
139 
140 // The default for the initial time step size of the simulation [s]
141 template<class TypeTag>
142 struct InitialTimeStepSize<TypeTag, TTag::ReservoirBaseProblem>
143 {
144  using type = GetPropType<TypeTag, Scalar>;
145  static constexpr type value = 100e3;
146 };
147 
148 // The width of producer/injector wells as a fraction of the width of the spatial domain
149 template<class TypeTag>
150 struct WellWidth<TypeTag, TTag::ReservoirBaseProblem>
151 {
152  using type = GetPropType<TypeTag, Scalar>;
153  static constexpr type value = 0.01;
154 };
155 
164 template<class TypeTag>
165 struct FluidSystem<TypeTag, TTag::ReservoirBaseProblem>
166 {
167 private:
168  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
169 
170 public:
171  using type = Opm::BlackOilFluidSystem<Scalar>;
172 };
173 
174 // The default DGF file to load
175 template<class TypeTag>
176 struct GridFile<TypeTag, TTag::ReservoirBaseProblem> { static constexpr auto value = "data/reservoir.dgf"; };
177 
178 // increase the tolerance for this problem to get larger time steps
179 template<class TypeTag>
180 struct NewtonTolerance<TypeTag, TTag::ReservoirBaseProblem>
181 {
182  using type = GetPropType<TypeTag, Scalar>;
183  static constexpr type value = 1e-6;
184 };
185 
186 } // namespace Opm::Properties
187 
188 namespace Opm {
189 
206 template <class TypeTag>
207 class ReservoirProblem : public GetPropType<TypeTag, Properties::BaseProblem>
208 {
209  using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
210 
211  using GridView = GetPropType<TypeTag, Properties::GridView>;
212  using Scalar = GetPropType<TypeTag, Properties::Scalar>;
213  using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
214  using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
215 
216  // Grid and world dimension
217  enum { dim = GridView::dimension };
218  enum { dimWorld = GridView::dimensionworld };
219 
220  // copy some indices for convenience
221  enum { numPhases = FluidSystem::numPhases };
222  enum { numComponents = FluidSystem::numComponents };
223  enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
224  enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
225  enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
226  enum { gasCompIdx = FluidSystem::gasCompIdx };
227  enum { oilCompIdx = FluidSystem::oilCompIdx };
228  enum { waterCompIdx = FluidSystem::waterCompIdx };
229 
230  using Model = GetPropType<TypeTag, Properties::Model>;
231  using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
232  using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
233  using EqVector = GetPropType<TypeTag, Properties::EqVector>;
234  using RateVector = GetPropType<TypeTag, Properties::RateVector>;
235  using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
236  using Constraints = GetPropType<TypeTag, Properties::Constraints>;
237  using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
238  using Simulator = GetPropType<TypeTag, Properties::Simulator>;
239  using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
240 
241  using CoordScalar = typename GridView::ctype;
242  using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
243  using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
244  using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
245 
246  using InitialFluidState = Opm::CompositionalFluidState<Scalar,
247  FluidSystem,
248  /*enableEnthalpy=*/true>;
249 
250 public:
254  ReservoirProblem(Simulator& simulator)
255  : ParentType(simulator)
256  { }
257 
261  void finishInit()
262  {
263  ParentType::finishInit();
264 
265  temperature_ = EWOMS_GET_PARAM(TypeTag, Scalar, Temperature);
266  maxDepth_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxDepth);
267  wellWidth_ = EWOMS_GET_PARAM(TypeTag, Scalar, WellWidth);
268 
269  std::vector<std::pair<Scalar, Scalar> > Bo = {
270  { 101353, 1.062 },
271  { 1.82504e+06, 1.15 },
272  { 3.54873e+06, 1.207 },
273  { 6.99611e+06, 1.295 },
274  { 1.38909e+07, 1.435 },
275  { 1.73382e+07, 1.5 },
276  { 2.07856e+07, 1.565 },
277  { 2.76804e+07, 1.695 },
278  { 3.45751e+07, 1.827 }
279  };
280  std::vector<std::pair<Scalar, Scalar> > muo = {
281  { 101353, 0.00104 },
282  { 1.82504e+06, 0.000975 },
283  { 3.54873e+06, 0.00091 },
284  { 6.99611e+06, 0.00083 },
285  { 1.38909e+07, 0.000695 },
286  { 1.73382e+07, 0.000641 },
287  { 2.07856e+07, 0.000594 },
288  { 2.76804e+07, 0.00051 },
289  { 3.45751e+07, 0.000449 }
290  };
291  std::vector<std::pair<Scalar, Scalar> > Rs = {
292  { 101353, 0.178108 },
293  { 1.82504e+06, 16.1187 },
294  { 3.54873e+06, 32.0594 },
295  { 6.99611e+06, 66.0779 },
296  { 1.38909e+07, 113.276 },
297  { 1.73382e+07, 138.033 },
298  { 2.07856e+07, 165.64 },
299  { 2.76804e+07, 226.197 },
300  { 3.45751e+07, 288.178 }
301  };
302  std::vector<std::pair<Scalar, Scalar> > Bg = {
303  { 101353, 0.93576 },
304  { 1.82504e+06, 0.0678972 },
305  { 3.54873e+06, 0.0352259 },
306  { 6.99611e+06, 0.0179498 },
307  { 1.38909e+07, 0.00906194 },
308  { 1.73382e+07, 0.00726527 },
309  { 2.07856e+07, 0.00606375 },
310  { 2.76804e+07, 0.00455343 },
311  { 3.45751e+07, 0.00364386 },
312  { 6.21542e+07, 0.00216723 }
313  };
314  std::vector<std::pair<Scalar, Scalar> > mug = {
315  { 101353, 8e-06 },
316  { 1.82504e+06, 9.6e-06 },
317  { 3.54873e+06, 1.12e-05 },
318  { 6.99611e+06, 1.4e-05 },
319  { 1.38909e+07, 1.89e-05 },
320  { 1.73382e+07, 2.08e-05 },
321  { 2.07856e+07, 2.28e-05 },
322  { 2.76804e+07, 2.68e-05 },
323  { 3.45751e+07, 3.09e-05 },
324  { 6.21542e+07, 4.7e-05 }
325  };
326 
327  Scalar rhoRefO = 786.0; // [kg]
328  Scalar rhoRefG = 0.97; // [kg]
329  Scalar rhoRefW = 1037.0; // [kg]
330  FluidSystem::initBegin(/*numPvtRegions=*/1);
331  FluidSystem::setEnableDissolvedGas(true);
332  FluidSystem::setEnableVaporizedOil(false);
333  FluidSystem::setReferenceDensities(rhoRefO, rhoRefW, rhoRefG, /*regionIdx=*/0);
334 
335  Opm::GasPvtMultiplexer<Scalar> *gasPvt = new Opm::GasPvtMultiplexer<Scalar>;
336  gasPvt->setApproach(GasPvtApproach::DryGasPvt);
337  auto& dryGasPvt = gasPvt->template getRealPvt<GasPvtApproach::DryGasPvt>();
338  dryGasPvt.setNumRegions(/*numPvtRegion=*/1);
339  dryGasPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW);
340  dryGasPvt.setGasFormationVolumeFactor(/*regionIdx=*/0, Bg);
341  dryGasPvt.setGasViscosity(/*regionIdx=*/0, mug);
342 
343  Opm::OilPvtMultiplexer<Scalar> *oilPvt = new Opm::OilPvtMultiplexer<Scalar>;
344  oilPvt->setApproach(OilPvtApproach::LiveOilPvt);
345  auto& liveOilPvt = oilPvt->template getRealPvt<OilPvtApproach::LiveOilPvt>();
346  liveOilPvt.setNumRegions(/*numPvtRegion=*/1);
347  liveOilPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW);
348  liveOilPvt.setSaturatedOilGasDissolutionFactor(/*regionIdx=*/0, Rs);
349  liveOilPvt.setSaturatedOilFormationVolumeFactor(/*regionIdx=*/0, Bo);
350  liveOilPvt.setSaturatedOilViscosity(/*regionIdx=*/0, muo);
351 
352  Opm::WaterPvtMultiplexer<Scalar> *waterPvt = new Opm::WaterPvtMultiplexer<Scalar>;
353  waterPvt->setApproach(WaterPvtApproach::ConstantCompressibilityWaterPvt);
354  auto& ccWaterPvt = waterPvt->template getRealPvt<WaterPvtApproach::ConstantCompressibilityWaterPvt>();
355  ccWaterPvt.setNumRegions(/*numPvtRegions=*/1);
356  ccWaterPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW);
357  ccWaterPvt.setViscosity(/*regionIdx=*/0, 9.6e-4);
358  ccWaterPvt.setCompressibility(/*regionIdx=*/0, 1.450377e-10);
359 
360  gasPvt->initEnd();
361  oilPvt->initEnd();
362  waterPvt->initEnd();
363 
364  using GasPvtSharedPtr = std::shared_ptr<Opm::GasPvtMultiplexer<Scalar> >;
365  FluidSystem::setGasPvt(GasPvtSharedPtr(gasPvt));
366 
367  using OilPvtSharedPtr = std::shared_ptr<Opm::OilPvtMultiplexer<Scalar> >;
368  FluidSystem::setOilPvt(OilPvtSharedPtr(oilPvt));
369 
370  using WaterPvtSharedPtr = std::shared_ptr<Opm::WaterPvtMultiplexer<Scalar> >;
371  FluidSystem::setWaterPvt(WaterPvtSharedPtr(waterPvt));
372 
373  FluidSystem::initEnd();
374 
375  pReservoir_ = 330e5;
376  layerBottom_ = 22.0;
377 
378  // intrinsic permeabilities
379  fineK_ = this->toDimMatrix_(1e-12);
380  coarseK_ = this->toDimMatrix_(1e-11);
381 
382  // porosities
383  finePorosity_ = 0.2;
384  coarsePorosity_ = 0.3;
385 
386  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
387  fineMaterialParams_.setPcMinSat(phaseIdx, 0.0);
388  fineMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
389 
390  coarseMaterialParams_.setPcMinSat(phaseIdx, 0.0);
391  coarseMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
392  }
393 
394  // wrap up the initialization of the material law's parameters
395  fineMaterialParams_.finalize();
396  coarseMaterialParams_.finalize();
397 
398  materialParams_.resize(this->model().numGridDof());
399  ElementContext elemCtx(this->simulator());
400  auto eIt = this->simulator().gridView().template begin<0>();
401  const auto& eEndIt = this->simulator().gridView().template end<0>();
402  for (; eIt != eEndIt; ++eIt) {
403  elemCtx.updateStencil(*eIt);
404  size_t nDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
405  for (unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
406  unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
407  const GlobalPosition& pos = elemCtx.pos(dofIdx, /*timeIdx=*/0);
408 
409  if (isFineMaterial_(pos))
410  materialParams_[globalDofIdx] = &fineMaterialParams_;
411  else
412  materialParams_[globalDofIdx] = &coarseMaterialParams_;
413  }
414  }
415 
416  initFluidState_();
417 
418  // start the first ("settle down") episode for 100 days
419  this->simulator().startNextEpisode(100.0*24*60*60);
420  }
421 
425  static void registerParameters()
426  {
427  ParentType::registerParameters();
428 
429  EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature,
430  "The temperature [K] in the reservoir");
431  EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxDepth,
432  "The maximum depth [m] of the reservoir");
433  EWOMS_REGISTER_PARAM(TypeTag, Scalar, WellWidth,
434  "The width of producer/injector wells as a fraction of the width"
435  " of the spatial domain");
436  }
437 
441  std::string name() const
442  { return std::string("reservoir_") + Model::name() + "_" + Model::discretizationName(); }
443 
447  void endEpisode()
448  {
449  // in the second episode, the actual work is done (the first is "settle down"
450  // episode). we need to use a pretty short initial time step here as the change
451  // in conditions is quite abrupt.
452  this->simulator().startNextEpisode(1e100);
453  this->simulator().setTimeStepSize(5.0);
454  }
455 
459  void endTimeStep()
460  {
461 #ifndef NDEBUG
462  // checkConservativeness() does not include the effect of constraints, so we
463  // disable it for this problem...
464  //this->model().checkConservativeness();
465 
466  // Calculate storage terms
467  EqVector storage;
468  this->model().globalStorage(storage);
469 
470  // Write mass balance information for rank 0
471  if (this->gridView().comm().rank() == 0) {
472  std::cout << "Storage: " << storage << std::endl << std::flush;
473  }
474 #endif // NDEBUG
475  }
476 
483  template <class Context>
484  const DimMatrix& intrinsicPermeability(const Context& context, unsigned spaceIdx,
485  unsigned timeIdx) const
486  {
487  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
488  if (isFineMaterial_(pos))
489  return fineK_;
490  return coarseK_;
491  }
492 
496  template <class Context>
497  Scalar porosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
498  {
499  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
500  if (isFineMaterial_(pos))
501  return finePorosity_;
502  return coarsePorosity_;
503  }
504 
508  template <class Context>
509  const MaterialLawParams& materialLawParams(const Context& context,
510  unsigned spaceIdx, unsigned timeIdx) const
511  {
512  unsigned globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
513  return *materialParams_[globalIdx];
514  }
515 
516  const MaterialLawParams& materialLawParams(unsigned globalIdx) const
517  { return *materialParams_[globalIdx]; }
518 
523 
524 
533  template <class Context>
534  Scalar temperature(const Context& /*context*/,
535  unsigned /*spaceIdx*/,
536  unsigned /*timeIdx*/) const
537  { return temperature_; }
538 
539  // \}
540 
545 
552  template <class Context>
553  void boundary(BoundaryRateVector& values,
554  const Context& /*context*/,
555  unsigned /*spaceIdx*/,
556  unsigned /*timeIdx*/) const
557  {
558  // no flow on top and bottom
559  values.setNoFlow();
560  }
561 
563 
568 
575  template <class Context>
576  void initial(PrimaryVariables& values,
577  const Context& /*context*/,
578  unsigned /*spaceIdx*/,
579  unsigned /*timeIdx*/) const
580  {
581  values.assignNaive(initialFluidState_);
582 
583 #ifndef NDEBUG
584  for (unsigned pvIdx = 0; pvIdx < values.size(); ++ pvIdx)
585  assert(std::isfinite(values[pvIdx]));
586 #endif
587  }
588 
597  template <class Context>
598  void constraints(Constraints& constraintValues,
599  const Context& context,
600  unsigned spaceIdx,
601  unsigned timeIdx) const
602  {
603  if (this->simulator().episodeIndex() == 1)
604  return; // no constraints during the "settle down" episode
605 
606  const auto& pos = context.pos(spaceIdx, timeIdx);
607  if (isInjector_(pos)) {
608  constraintValues.setActive(true);
609  constraintValues.assignNaive(injectorFluidState_);
610  }
611  else if (isProducer_(pos)) {
612  constraintValues.setActive(true);
613  constraintValues.assignNaive(producerFluidState_);
614  }
615  }
616 
622  template <class Context>
623  void source(RateVector& rate,
624  const Context& /*context*/,
625  unsigned /*spaceIdx*/,
626  unsigned /*timeIdx*/) const
627  { rate = Scalar(0.0); }
628 
630 
631 private:
632  void initFluidState_()
633  {
634  auto& fs = initialFluidState_;
635 
637  // set temperatures
639  fs.setTemperature(temperature_);
640 
642  // set saturations
644  fs.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
645  fs.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
646  fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
647 
649  // set pressures
651  Scalar pw = pReservoir_;
652 
653  PhaseVector pC;
654  const auto& matParams = fineMaterialParams_;
655  MaterialLaw::capillaryPressures(pC, matParams, fs);
656 
657  fs.setPressure(oilPhaseIdx, pw + (pC[oilPhaseIdx] - pC[waterPhaseIdx]));
658  fs.setPressure(waterPhaseIdx, pw + (pC[waterPhaseIdx] - pC[waterPhaseIdx]));
659  fs.setPressure(gasPhaseIdx, pw + (pC[gasPhaseIdx] - pC[waterPhaseIdx]));
660 
661  // reset all mole fractions to 0
662  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
663  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
664  fs.setMoleFraction(phaseIdx, compIdx, 0.0);
665 
667  // set composition of the gas and water phases
669  fs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
670  fs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
671 
673  // set composition of the oil phase
675  Scalar RsSat =
676  FluidSystem::saturatedDissolutionFactor(fs, oilPhaseIdx, /*pvtRegionIdx=*/0);
677  Scalar XoGSat = FluidSystem::convertRsToXoG(RsSat, /*pvtRegionIdx=*/0);
678  Scalar xoGSat = FluidSystem::convertXoGToxoG(XoGSat, /*pvtRegionIdx=*/0);
679  Scalar xoG = 0.95*xoGSat;
680  Scalar xoO = 1.0 - xoG;
681 
682  // finally set the oil-phase composition
683  fs.setMoleFraction(oilPhaseIdx, gasCompIdx, xoG);
684  fs.setMoleFraction(oilPhaseIdx, oilCompIdx, xoO);
685 
686  using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
687  typename FluidSystem::template ParameterCache<Scalar> paramCache;
688  CFRP::solve(fs,
689  paramCache,
690  /*refPhaseIdx=*/oilPhaseIdx,
691  /*setViscosities=*/false,
692  /*setEnthalpies=*/false);
693 
694  // set up the fluid state used for the injectors
695  auto& injFs = injectorFluidState_;
696  injFs = initialFluidState_;
697 
698  Scalar pInj = pReservoir_ * 1.5;
699  injFs.setPressure(waterPhaseIdx, pInj);
700  injFs.setPressure(oilPhaseIdx, pInj);
701  injFs.setPressure(gasPhaseIdx, pInj);
702  injFs.setSaturation(waterPhaseIdx, 1.0);
703  injFs.setSaturation(oilPhaseIdx, 0.0);
704  injFs.setSaturation(gasPhaseIdx, 0.0);
705 
706  // set the composition of the phases to immiscible
707  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
708  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
709  injFs.setMoleFraction(phaseIdx, compIdx, 0.0);
710 
711  injFs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
712  injFs.setMoleFraction(oilPhaseIdx, oilCompIdx, 1.0);
713  injFs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
714 
715  CFRP::solve(injFs,
716  paramCache,
717  /*refPhaseIdx=*/waterPhaseIdx,
718  /*setViscosities=*/true,
719  /*setEnthalpies=*/false);
720 
721  // set up the fluid state used for the producer
722  auto& prodFs = producerFluidState_;
723  prodFs = initialFluidState_;
724 
725  Scalar pProd = pReservoir_ / 1.5;
726  prodFs.setPressure(waterPhaseIdx, pProd);
727  prodFs.setPressure(oilPhaseIdx, pProd);
728  prodFs.setPressure(gasPhaseIdx, pProd);
729  prodFs.setSaturation(waterPhaseIdx, 0.0);
730  prodFs.setSaturation(oilPhaseIdx, 1.0);
731  prodFs.setSaturation(gasPhaseIdx, 0.0);
732 
733  CFRP::solve(prodFs,
734  paramCache,
735  /*refPhaseIdx=*/oilPhaseIdx,
736  /*setViscosities=*/true,
737  /*setEnthalpies=*/false);
738  }
739 
740  bool isProducer_(const GlobalPosition& pos) const
741  {
742  Scalar x = pos[0] - this->boundingBoxMin()[0];
743  Scalar y = pos[dim - 1] - this->boundingBoxMin()[dim - 1];
744  Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
745  Scalar height = this->boundingBoxMax()[dim - 1] - this->boundingBoxMin()[dim - 1];
746 
747  // only the upper half of the center section of the spatial domain is assumed to
748  // be the producer
749  if (y <= height/2.0)
750  return false;
751 
752  return width/2.0 - width*1e-5 < x && x < width/2.0 + width*(wellWidth_ + 1e-5);
753  }
754 
755  bool isInjector_(const GlobalPosition& pos) const
756  {
757  Scalar x = pos[0] - this->boundingBoxMin()[0];
758  Scalar y = pos[dim - 1] - this->boundingBoxMin()[dim - 1];
759  Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
760  Scalar height = this->boundingBoxMax()[dim - 1] - this->boundingBoxMin()[dim - 1];
761 
762  // only the lower half of the leftmost and rightmost part of the spatial domain
763  // are assumed to be the water injectors
764  if (y > height/2.0)
765  return false;
766 
767  return x < width*wellWidth_ - width*1e-5 || x > width*(1.0 - wellWidth_) + width*1e-5;
768  }
769 
770  bool isFineMaterial_(const GlobalPosition& pos) const
771  { return pos[dim - 1] > layerBottom_; }
772 
773  DimMatrix fineK_;
774  DimMatrix coarseK_;
775  Scalar layerBottom_;
776  Scalar pReservoir_;
777 
778  Scalar finePorosity_;
779  Scalar coarsePorosity_;
780 
781  MaterialLawParams fineMaterialParams_;
782  MaterialLawParams coarseMaterialParams_;
783  std::vector<const MaterialLawParams*> materialParams_;
784 
785  InitialFluidState initialFluidState_;
786  InitialFluidState injectorFluidState_;
787  InitialFluidState producerFluidState_;
788 
789  Scalar temperature_;
790  Scalar maxDepth_;
791  Scalar wellWidth_;
792 };
793 } // namespace Opm
794 
795 #endif
Some simple test problem for the black-oil VCVF discretization inspired by an oil reservoir.
Definition: reservoirproblem.hh:208
static void registerParameters()
Definition: reservoirproblem.hh:425
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Definition: reservoirproblem.hh:623
void endTimeStep()
Definition: reservoirproblem.hh:459
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:497
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Definition: reservoirproblem.hh:576
ReservoirProblem(Simulator &simulator)
Definition: reservoirproblem.hh:254
void constraints(Constraints &constraintValues, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:598
void boundary(BoundaryRateVector &values, const Context &, unsigned, unsigned) const
Definition: reservoirproblem.hh:553
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:509
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:484
std::string name() const
Definition: reservoirproblem.hh:441
void finishInit()
Definition: reservoirproblem.hh:261
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: reservoirproblem.hh:534
void endEpisode()
Definition: reservoirproblem.hh:447
Definition: co2injectionproblem.hh:91
Definition: reservoirproblem.hh:63
Definition: co2injectionproblem.hh:93
Definition: reservoirproblem.hh:75