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
52namespace Opm {
53template <class TypeTag>
54class ReservoirProblem;
55
56} // namespace Opm
57
58namespace Opm::Properties {
59
60
61namespace TTag {
62
64
65} // namespace TTag
66
67// Maximum depth of the reservoir
68template<class TypeTag, class MyTypeTag>
69struct MaxDepth { using type = UndefinedProperty; };
70// The temperature inside the reservoir
71template<class TypeTag, class MyTypeTag>
72struct Temperature { using type = UndefinedProperty; };
73// The width of producer/injector wells as a fraction of the width of the spatial domain
74template<class TypeTag, class MyTypeTag>
75struct WellWidth { using type = UndefinedProperty; };
76
77// Set the grid type
78template<class TypeTag>
79struct Grid<TypeTag, TTag::ReservoirBaseProblem> { using type = Dune::YaspGrid<2>; };
80
81// Set the problem property
82template<class TypeTag>
83struct Problem<TypeTag, TTag::ReservoirBaseProblem> { using type = Opm::ReservoirProblem<TypeTag>; };
84
85// Set the material Law
86template<class TypeTag>
87struct MaterialLaw<TypeTag, TTag::ReservoirBaseProblem>
88{
89private:
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
99public:
100 using type = Opm::LinearMaterial<Traits>;
101};
102
103// Write the Newton convergence behavior to disk?
104template<class TypeTag>
105struct NewtonWriteConvergence<TypeTag, TTag::ReservoirBaseProblem> { static constexpr bool value = false; };
106
107// Enable gravity
108template<class TypeTag>
109struct EnableGravity<TypeTag, TTag::ReservoirBaseProblem> { static constexpr bool value = true; };
110
111// Enable constraint DOFs?
112template<class TypeTag>
113struct EnableConstraints<TypeTag, TTag::ReservoirBaseProblem> { static constexpr bool value = true; };
114
115// set the defaults for some problem specific properties
116template<class TypeTag>
117struct MaxDepth<TypeTag, TTag::ReservoirBaseProblem>
118{
119 using type = GetPropType<TypeTag, Scalar>;
120 static constexpr type value = 2500;
121};
122template<class TypeTag>
123struct Temperature<TypeTag, TTag::ReservoirBaseProblem>
124{
125 using type = GetPropType<TypeTag, Scalar>;
126 static constexpr type value = 293.15;
127};
128
133template<class TypeTag>
134struct 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]
141template<class TypeTag>
142struct 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
149template<class TypeTag>
150struct WellWidth<TypeTag, TTag::ReservoirBaseProblem>
151{
152 using type = GetPropType<TypeTag, Scalar>;
153 static constexpr type value = 0.01;
154};
155
164template<class TypeTag>
165struct FluidSystem<TypeTag, TTag::ReservoirBaseProblem>
166{
167private:
168 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
169
170public:
171 using type = Opm::BlackOilFluidSystem<Scalar>;
172};
173
174// The default DGF file to load
175template<class TypeTag>
176struct GridFile<TypeTag, TTag::ReservoirBaseProblem> { static constexpr auto value = "data/reservoir.dgf"; };
177
178// increase the tolerance for this problem to get larger time steps
179template<class TypeTag>
180struct 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
188namespace Opm {
189
206template <class TypeTag>
207class 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
250public:
254 ReservoirProblem(Simulator& simulator)
255 : ParentType(simulator)
256 { }
257
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
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
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
631private:
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
std::string name() const
Definition: reservoirproblem.hh:441
void finishInit()
Definition: reservoirproblem.hh:261
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:484
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