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 
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/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 Ewoms {
54 template <class TypeTag>
56 
57 namespace Properties {
58 
59 NEW_TYPE_TAG(ReservoirBaseProblem);
60 
61 // Maximum depth of the reservoir
62 NEW_PROP_TAG(MaxDepth);
63 // The temperature inside the reservoir
64 NEW_PROP_TAG(Temperature);
65 // The width of producer/injector wells as a fraction of the width of the spatial domain
66 NEW_PROP_TAG(WellWidth);
67 
68 // Set the grid type
69 SET_TYPE_PROP(ReservoirBaseProblem, Grid, Dune::YaspGrid<2>);
70 
71 // Set the problem property
72 SET_TYPE_PROP(ReservoirBaseProblem, Problem, Ewoms::ReservoirProblem<TypeTag>);
73 
74 // Set the material Law
75 SET_PROP(ReservoirBaseProblem, MaterialLaw)
76 {
77 private:
78  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
79  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
80 
81  typedef Opm::
82  ThreePhaseMaterialTraits<Scalar,
83  /*wettingPhaseIdx=*/FluidSystem::waterPhaseIdx,
84  /*nonWettingPhaseIdx=*/FluidSystem::oilPhaseIdx,
85  /*gasPhaseIdx=*/FluidSystem::gasPhaseIdx> Traits;
86 
87 public:
88  typedef Opm::LinearMaterial<Traits> type;
89 };
90 
91 // Write the Newton convergence behavior to disk?
92 SET_BOOL_PROP(ReservoirBaseProblem, NewtonWriteConvergence, false);
93 
94 // Enable gravity
95 SET_BOOL_PROP(ReservoirBaseProblem, EnableGravity, true);
96 
97 // Enable constraint DOFs?
98 SET_BOOL_PROP(ReservoirBaseProblem, EnableConstraints, true);
99 
100 // set the defaults for some problem specific properties
101 SET_SCALAR_PROP(ReservoirBaseProblem, MaxDepth, 2500);
102 SET_SCALAR_PROP(ReservoirBaseProblem, Temperature, 293.15);
103 
108 SET_SCALAR_PROP(ReservoirBaseProblem, EndTime, 1000.0*24*60*60);
109 
110 // The default for the initial time step size of the simulation [s]
111 SET_SCALAR_PROP(ReservoirBaseProblem, InitialTimeStepSize, 100e3);
112 
113 // The width of producer/injector wells as a fraction of the width of the spatial domain
114 SET_SCALAR_PROP(ReservoirBaseProblem, WellWidth, 0.01);
115 
124 SET_PROP(ReservoirBaseProblem, FluidSystem)
125 {
126 private:
127  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
128 
129 public:
130  typedef Opm::FluidSystems::BlackOil<Scalar> type;
131 };
132 
133 // The default DGF file to load
134 SET_STRING_PROP(ReservoirBaseProblem, GridFile, "data/reservoir.dgf");
135 
136 // increase the tolerance for this problem to get larger time steps
137 SET_SCALAR_PROP(ReservoirBaseProblem, NewtonRawTolerance, 1e-6);
138 } // namespace Properties
139 
156 template <class TypeTag>
157 class ReservoirProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
158 {
159  typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
160 
161  typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
162  typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
163  typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
164  typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
165 
166  // Grid and world dimension
167  enum { dim = GridView::dimension };
168  enum { dimWorld = GridView::dimensionworld };
169 
170  // copy some indices for convenience
171  enum { numPhases = FluidSystem::numPhases };
172  enum { numComponents = FluidSystem::numComponents };
173  enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
174  enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
175  enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
176  enum { gasCompIdx = FluidSystem::gasCompIdx };
177  enum { oilCompIdx = FluidSystem::oilCompIdx };
178  enum { waterCompIdx = FluidSystem::waterCompIdx };
179 
180  typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
181  typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
182  typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
183  typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
184  typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
185  typedef typename GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
186  typedef typename GET_PROP_TYPE(TypeTag, Constraints) Constraints;
187  typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
188  typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
189  typedef typename GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
190 
191  typedef typename GridView::ctype CoordScalar;
192  typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
193  typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
194  typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
195 
196  typedef Opm::CompositionalFluidState<Scalar,
197  FluidSystem,
198  /*enableEnthalpy=*/true> InitialFluidState;
199 
200 public:
205  : ParentType(simulator)
206  { }
207 
211  void finishInit()
212  {
213  ParentType::finishInit();
214 
215  temperature_ = EWOMS_GET_PARAM(TypeTag, Scalar, Temperature);
216  maxDepth_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxDepth);
217  wellWidth_ = EWOMS_GET_PARAM(TypeTag, Scalar, WellWidth);
218 
219  std::vector<std::pair<Scalar, Scalar> > Bo = {
220  { 101353, 1.062 },
221  { 1.82504e+06, 1.15 },
222  { 3.54873e+06, 1.207 },
223  { 6.99611e+06, 1.295 },
224  { 1.38909e+07, 1.435 },
225  { 1.73382e+07, 1.5 },
226  { 2.07856e+07, 1.565 },
227  { 2.76804e+07, 1.695 },
228  { 3.45751e+07, 1.827 }
229  };
230  std::vector<std::pair<Scalar, Scalar> > muo = {
231  { 101353, 0.00104 },
232  { 1.82504e+06, 0.000975 },
233  { 3.54873e+06, 0.00091 },
234  { 6.99611e+06, 0.00083 },
235  { 1.38909e+07, 0.000695 },
236  { 1.73382e+07, 0.000641 },
237  { 2.07856e+07, 0.000594 },
238  { 2.76804e+07, 0.00051 },
239  { 3.45751e+07, 0.000449 }
240  };
241  std::vector<std::pair<Scalar, Scalar> > Rs = {
242  { 101353, 0.178108 },
243  { 1.82504e+06, 16.1187 },
244  { 3.54873e+06, 32.0594 },
245  { 6.99611e+06, 66.0779 },
246  { 1.38909e+07, 113.276 },
247  { 1.73382e+07, 138.033 },
248  { 2.07856e+07, 165.64 },
249  { 2.76804e+07, 226.197 },
250  { 3.45751e+07, 288.178 }
251  };
252  std::vector<std::pair<Scalar, Scalar> > Bg = {
253  { 101353, 0.93576 },
254  { 1.82504e+06, 0.0678972 },
255  { 3.54873e+06, 0.0352259 },
256  { 6.99611e+06, 0.0179498 },
257  { 1.38909e+07, 0.00906194 },
258  { 1.73382e+07, 0.00726527 },
259  { 2.07856e+07, 0.00606375 },
260  { 2.76804e+07, 0.00455343 },
261  { 3.45751e+07, 0.00364386 },
262  { 6.21542e+07, 0.00216723 }
263  };
264  std::vector<std::pair<Scalar, Scalar> > mug = {
265  { 101353, 8e-06 },
266  { 1.82504e+06, 9.6e-06 },
267  { 3.54873e+06, 1.12e-05 },
268  { 6.99611e+06, 1.4e-05 },
269  { 1.38909e+07, 1.89e-05 },
270  { 1.73382e+07, 2.08e-05 },
271  { 2.07856e+07, 2.28e-05 },
272  { 2.76804e+07, 2.68e-05 },
273  { 3.45751e+07, 3.09e-05 },
274  { 6.21542e+07, 4.7e-05 }
275  };
276 
277  Scalar rhoRefO = 786.0; // [kg]
278  Scalar rhoRefG = 0.97; // [kg]
279  Scalar rhoRefW = 1037.0; // [kg]
280  FluidSystem::initBegin(/*numPvtRegions=*/1);
281  FluidSystem::setEnableDissolvedGas(true);
282  FluidSystem::setEnableVaporizedOil(false);
283  FluidSystem::setReferenceDensities(rhoRefO, rhoRefW, rhoRefG, /*regionIdx=*/0);
284 
285  Opm::GasPvtMultiplexer<Scalar> *gasPvt = new Opm::GasPvtMultiplexer<Scalar>;
286  gasPvt->setApproach(Opm::GasPvtMultiplexer<Scalar>::DryGasPvt);
287  auto& dryGasPvt = gasPvt->template getRealPvt<Opm::GasPvtMultiplexer<Scalar>::DryGasPvt>();
288  dryGasPvt.setNumRegions(/*numPvtRegion=*/1);
289  dryGasPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW);
290  dryGasPvt.setGasFormationVolumeFactor(/*regionIdx=*/0, Bg);
291  dryGasPvt.setGasViscosity(/*regionIdx=*/0, mug);
292 
293  Opm::OilPvtMultiplexer<Scalar> *oilPvt = new Opm::OilPvtMultiplexer<Scalar>;
294  oilPvt->setApproach(Opm::OilPvtMultiplexer<Scalar>::LiveOilPvt);
295  auto& liveOilPvt = oilPvt->template getRealPvt<Opm::OilPvtMultiplexer<Scalar>::LiveOilPvt>();
296  liveOilPvt.setNumRegions(/*numPvtRegion=*/1);
297  liveOilPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW);
298  liveOilPvt.setSaturatedOilGasDissolutionFactor(/*regionIdx=*/0, Rs);
299  liveOilPvt.setSaturatedOilFormationVolumeFactor(/*regionIdx=*/0, Bo);
300  liveOilPvt.setSaturatedOilViscosity(/*regionIdx=*/0, muo);
301 
302  Opm::WaterPvtMultiplexer<Scalar> *waterPvt = new Opm::WaterPvtMultiplexer<Scalar>;
303  waterPvt->setApproach(Opm::WaterPvtMultiplexer<Scalar>::ConstantCompressibilityWaterPvt);
304  auto& ccWaterPvt = waterPvt->template getRealPvt<Opm::WaterPvtMultiplexer<Scalar>::ConstantCompressibilityWaterPvt>();
305  ccWaterPvt.setNumRegions(/*numPvtRegions=*/1);
306  ccWaterPvt.setReferenceDensities(/*regionIdx=*/0, rhoRefO, rhoRefG, rhoRefW);
307  ccWaterPvt.setViscosity(/*regionIdx=*/0, 9.6e-4);
308  ccWaterPvt.setCompressibility(/*regionIdx=*/0, 1.450377e-10);
309 
310  gasPvt->initEnd();
311  oilPvt->initEnd();
312  waterPvt->initEnd();
313 
314  typedef std::shared_ptr<Opm::GasPvtMultiplexer<Scalar> > GasPvtSharedPtr;
315  FluidSystem::setGasPvt(GasPvtSharedPtr(gasPvt));
316 
317  typedef std::shared_ptr<Opm::OilPvtMultiplexer<Scalar> > OilPvtSharedPtr;
318  FluidSystem::setOilPvt(OilPvtSharedPtr(oilPvt));
319 
320  typedef std::shared_ptr<Opm::WaterPvtMultiplexer<Scalar> > WaterPvtSharedPtr;
321  FluidSystem::setWaterPvt(WaterPvtSharedPtr(waterPvt));
322 
323  FluidSystem::initEnd();
324 
325  pReservoir_ = 330e5;
326  layerBottom_ = 22.0;
327 
328  // intrinsic permeabilities
329  fineK_ = this->toDimMatrix_(1e-12);
330  coarseK_ = this->toDimMatrix_(1e-11);
331 
332  // porosities
333  finePorosity_ = 0.2;
334  coarsePorosity_ = 0.3;
335 
336  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
337  fineMaterialParams_.setPcMinSat(phaseIdx, 0.0);
338  fineMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
339 
340  coarseMaterialParams_.setPcMinSat(phaseIdx, 0.0);
341  coarseMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
342  }
343 
344  // wrap up the initialization of the material law's parameters
345  fineMaterialParams_.finalize();
346  coarseMaterialParams_.finalize();
347 
348  materialParams_.resize(this->model().numGridDof());
349  ElementContext elemCtx(this->simulator());
350  auto eIt = this->simulator().gridView().template begin<0>();
351  const auto& eEndIt = this->simulator().gridView().template end<0>();
352  for (; eIt != eEndIt; ++eIt) {
353  elemCtx.updateStencil(*eIt);
354  size_t nDof = elemCtx.numPrimaryDof(/*timeIdx=*/0);
355  for (unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
356  unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
357  const GlobalPosition& pos = elemCtx.pos(dofIdx, /*timeIdx=*/0);
358 
359  if (isFineMaterial_(pos))
360  materialParams_[globalDofIdx] = &fineMaterialParams_;
361  else
362  materialParams_[globalDofIdx] = &coarseMaterialParams_;
363  }
364  }
365 
366  initFluidState_();
367 
368  // start the first ("settle down") episode for 100 days
369  this->simulator().startNextEpisode(100.0*24*60*60);
370  }
371 
375  static void registerParameters()
376  {
377  ParentType::registerParameters();
378 
379  EWOMS_REGISTER_PARAM(TypeTag, Scalar, Temperature,
380  "The temperature [K] in the reservoir");
381  EWOMS_REGISTER_PARAM(TypeTag, Scalar, MaxDepth,
382  "The maximum depth [m] of the reservoir");
383  EWOMS_REGISTER_PARAM(TypeTag, Scalar, WellWidth,
384  "The width of producer/injector wells as a fraction of the width"
385  " of the spatial domain");
386  }
387 
391  std::string name() const
392  { return std::string("reservoir_") + Model::name() + "_" + Model::discretizationName(); }
393 
397  void endEpisode()
398  {
399  // in the second episode, the actual work is done (the first is "settle down"
400  // episode). we need to use a pretty short initial time step here as the change
401  // in conditions is quite abrupt.
402  this->simulator().startNextEpisode(1e100);
403  this->simulator().setTimeStepSize(5.0);
404  }
405 
409  void endTimeStep()
410  {
411 #ifndef NDEBUG
412  // checkConservativeness() does not include the effect of constraints, so we
413  // disable it for this problem...
414  //this->model().checkConservativeness();
415 
416  // Calculate storage terms
417  EqVector storage;
418  this->model().globalStorage(storage);
419 
420  // Write mass balance information for rank 0
421  if (this->gridView().comm().rank() == 0) {
422  std::cout << "Storage: " << storage << std::endl << std::flush;
423  }
424 #endif // NDEBUG
425  }
426 
433  template <class Context>
434  const DimMatrix& intrinsicPermeability(const Context& context, unsigned spaceIdx,
435  unsigned timeIdx) const
436  {
437  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
438  if (isFineMaterial_(pos))
439  return fineK_;
440  return coarseK_;
441  }
442 
446  template <class Context>
447  Scalar porosity(const Context& context, unsigned spaceIdx, unsigned timeIdx) const
448  {
449  const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
450  if (isFineMaterial_(pos))
451  return finePorosity_;
452  return coarsePorosity_;
453  }
454 
458  template <class Context>
459  const MaterialLawParams& materialLawParams(const Context& context,
460  unsigned spaceIdx, unsigned timeIdx) const
461  {
462  unsigned globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
463  return *materialParams_[globalIdx];
464  }
465 
466  const MaterialLawParams& materialLawParams(unsigned globalIdx) const
467  { return *materialParams_[globalIdx]; }
468 
472 
474 
483  template <class Context>
484  Scalar temperature(const Context& context OPM_UNUSED,
485  unsigned spaceIdx OPM_UNUSED,
486  unsigned timeIdx OPM_UNUSED) const
487  { return temperature_; }
488 
489  // \}
490 
494 
502  template <class Context>
503  void boundary(BoundaryRateVector& values,
504  const Context& context OPM_UNUSED,
505  unsigned spaceIdx OPM_UNUSED,
506  unsigned timeIdx OPM_UNUSED) const
507  {
508  // no flow on top and bottom
509  values.setNoFlow();
510  }
511 
513 
517 
525  template <class Context>
526  void initial(PrimaryVariables& values,
527  const Context& context OPM_UNUSED,
528  unsigned spaceIdx OPM_UNUSED,
529  unsigned timeIdx OPM_UNUSED) const
530  {
531  values.assignNaive(initialFluidState_);
532 
533 #ifndef NDEBUG
534  for (unsigned pvIdx = 0; pvIdx < values.size(); ++ pvIdx)
535  assert(std::isfinite(values[pvIdx]));
536 #endif
537  }
538 
547  template <class Context>
548  void constraints(Constraints& constraints, const Context& context,
549  unsigned spaceIdx, unsigned timeIdx) const
550  {
551  if (this->simulator().episodeIndex() == 1)
552  return; // no constraints during the "settle down" episode
553 
554  const auto& pos = context.pos(spaceIdx, timeIdx);
555  if (isInjector_(pos)) {
556  constraints.setActive(true);
557  constraints.assignNaive(injectorFluidState_);
558  }
559  else if (isProducer_(pos)) {
560  constraints.setActive(true);
561  constraints.assignNaive(producerFluidState_);
562  }
563  }
564 
570  template <class Context>
571  void source(RateVector& rate,
572  const Context& context OPM_UNUSED,
573  unsigned spaceIdx OPM_UNUSED,
574  unsigned timeIdx OPM_UNUSED) const
575  { rate = Scalar(0.0); }
576 
578 
579 private:
580  void initFluidState_()
581  {
582  auto& fs = initialFluidState_;
583 
585  // set temperatures
587  fs.setTemperature(temperature_);
588 
590  // set saturations
592  fs.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
593  fs.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
594  fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
595 
597  // set pressures
599  Scalar pw = pReservoir_;
600 
601  PhaseVector pC;
602  const auto& matParams = fineMaterialParams_;
603  MaterialLaw::capillaryPressures(pC, matParams, fs);
604 
605  fs.setPressure(oilPhaseIdx, pw + (pC[oilPhaseIdx] - pC[waterPhaseIdx]));
606  fs.setPressure(waterPhaseIdx, pw + (pC[waterPhaseIdx] - pC[waterPhaseIdx]));
607  fs.setPressure(gasPhaseIdx, pw + (pC[gasPhaseIdx] - pC[waterPhaseIdx]));
608 
609  // reset all mole fractions to 0
610  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
611  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
612  fs.setMoleFraction(phaseIdx, compIdx, 0.0);
613 
615  // set composition of the gas and water phases
617  fs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
618  fs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
619 
621  // set composition of the oil phase
623  Scalar RsSat =
624  FluidSystem::saturatedDissolutionFactor(fs, oilPhaseIdx, /*pvtRegionIdx=*/0);
625  Scalar XoGSat = FluidSystem::convertRsToXoG(RsSat, /*pvtRegionIdx=*/0);
626  Scalar xoGSat = FluidSystem::convertXoGToxoG(XoGSat, /*pvtRegionIdx=*/0);
627  Scalar xoG = 0.95*xoGSat;
628  Scalar xoO = 1.0 - xoG;
629 
630  // finally set the oil-phase composition
631  fs.setMoleFraction(oilPhaseIdx, gasCompIdx, xoG);
632  fs.setMoleFraction(oilPhaseIdx, oilCompIdx, xoO);
633 
634  typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> CFRP;
635  typename FluidSystem::template ParameterCache<Scalar> paramCache;
636  CFRP::solve(fs,
637  paramCache,
638  /*refPhaseIdx=*/oilPhaseIdx,
639  /*setViscosities=*/false,
640  /*setEnthalpies=*/false);
641 
642  // set up the fluid state used for the injectors
643  auto& injFs = injectorFluidState_;
644  injFs = initialFluidState_;
645 
646  Scalar pInj = pReservoir_ * 1.5;
647  injFs.setPressure(waterPhaseIdx, pInj);
648  injFs.setPressure(oilPhaseIdx, pInj);
649  injFs.setPressure(gasPhaseIdx, pInj);
650  injFs.setSaturation(waterPhaseIdx, 1.0);
651  injFs.setSaturation(oilPhaseIdx, 0.0);
652  injFs.setSaturation(gasPhaseIdx, 0.0);
653 
654  // set the composition of the phases to immiscible
655  for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
656  for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
657  injFs.setMoleFraction(phaseIdx, compIdx, 0.0);
658 
659  injFs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
660  injFs.setMoleFraction(oilPhaseIdx, oilCompIdx, 1.0);
661  injFs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
662 
663  CFRP::solve(injFs,
664  paramCache,
665  /*refPhaseIdx=*/waterPhaseIdx,
666  /*setViscosities=*/false,
667  /*setEnthalpies=*/false);
668 
669  // set up the fluid state used for the producer
670  auto& prodFs = producerFluidState_;
671  prodFs = initialFluidState_;
672 
673  Scalar pProd = pReservoir_ / 1.5;
674  prodFs.setPressure(waterPhaseIdx, pProd);
675  prodFs.setPressure(oilPhaseIdx, pProd);
676  prodFs.setPressure(gasPhaseIdx, pProd);
677  prodFs.setSaturation(waterPhaseIdx, 0.0);
678  prodFs.setSaturation(oilPhaseIdx, 1.0);
679  prodFs.setSaturation(gasPhaseIdx, 0.0);
680 
681  CFRP::solve(prodFs,
682  paramCache,
683  /*refPhaseIdx=*/oilPhaseIdx,
684  /*setViscosities=*/false,
685  /*setEnthalpies=*/false);
686  }
687 
688  bool isProducer_(const GlobalPosition& pos) const
689  {
690  Scalar x = pos[0] - this->boundingBoxMin()[0];
691  Scalar y = pos[dim - 1] - this->boundingBoxMin()[dim - 1];
692  Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
693  Scalar height = this->boundingBoxMax()[dim - 1] - this->boundingBoxMin()[dim - 1];
694 
695  // only the upper half of the center section of the spatial domain is assumed to
696  // be the producer
697  if (y <= height/2.0)
698  return false;
699 
700  return width/2.0 - width*1e-5 < x && x < width/2.0 + width*(wellWidth_ + 1e-5);
701  }
702 
703  bool isInjector_(const GlobalPosition& pos) const
704  {
705  Scalar x = pos[0] - this->boundingBoxMin()[0];
706  Scalar y = pos[dim - 1] - this->boundingBoxMin()[dim - 1];
707  Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
708  Scalar height = this->boundingBoxMax()[dim - 1] - this->boundingBoxMin()[dim - 1];
709 
710  // only the lower half of the leftmost and rightmost part of the spatial domain
711  // are assumed to be the water injectors
712  if (y > height/2.0)
713  return false;
714 
715  return x < width*wellWidth_ - width*1e-5 || x > width*(1.0 - wellWidth_) + width*1e-5;
716  }
717 
718  bool isFineMaterial_(const GlobalPosition& pos) const
719  { return pos[dim - 1] > layerBottom_; }
720 
721  DimMatrix fineK_;
722  DimMatrix coarseK_;
723  Scalar layerBottom_;
724  Scalar pReservoir_;
725 
726  Scalar finePorosity_;
727  Scalar coarsePorosity_;
728 
729  MaterialLawParams fineMaterialParams_;
730  MaterialLawParams coarseMaterialParams_;
731  std::vector<const MaterialLawParams*> materialParams_;
732 
733  InitialFluidState initialFluidState_;
734  InitialFluidState injectorFluidState_;
735  InitialFluidState producerFluidState_;
736 
737  Scalar temperature_;
738  Scalar maxDepth_;
739  Scalar wellWidth_;
740 };
741 } // namespace Ewoms
742 
743 #endif
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:459
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
void endTimeStep()
Called by the simulator after each time integration.
Definition: reservoirproblem.hh:409
Definition: baseauxiliarymodule.hh:37
Some simple test problem for the black-oil VCVF discretization inspired by an oil reservoir...
Definition: reservoirproblem.hh:55
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
Declares the properties required by the black oil model.
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: reservoirproblem.hh:484
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: reservoirproblem.hh:571
void boundary(BoundaryRateVector &values, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the boundary conditions for a boundary segment.
Definition: reservoirproblem.hh:503
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the constraints for a control volume.
Definition: reservoirproblem.hh:548
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:447
static void registerParameters()
Definition: reservoirproblem.hh:375
std::string name() const
The problem name.
Definition: reservoirproblem.hh:391
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:434
ReservoirProblem(Simulator &simulator)
Definition: reservoirproblem.hh:204
void initial(PrimaryVariables &values, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the initial value for a control volume.
Definition: reservoirproblem.hh:526
Manages the initializing and running of time dependent problems.
Definition: simulator.hh:75
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
void endEpisode()
Called when the end of an simulation episode is reached.
Definition: reservoirproblem.hh:397
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
#define SET_STRING_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant string value.
Definition: propertysystem.hh:416
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: reservoirproblem.hh:211
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394