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