28 #ifndef EWOMS_RESERVOIR_PROBLEM_HH
29 #define EWOMS_RESERVOIR_PROBLEM_HH
31 #include <opm/models/blackoil/blackoilproperties.hh>
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>
43 #include <dune/grid/yaspgrid.hh>
44 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
46 #include <dune/common/version.hh>
47 #include <dune/common/fvector.hh>
48 #include <dune/common/fmatrix.hh>
54 template <
class TypeTag>
55 class ReservoirProblem;
59 namespace Opm::Properties {
69 template<
class TypeTag,
class MyTypeTag>
70 struct MaxDepth {
using type = UndefinedProperty; };
72 template<
class TypeTag,
class MyTypeTag>
73 struct Temperature {
using type = UndefinedProperty; };
75 template<
class TypeTag,
class MyTypeTag>
76 struct WellWidth {
using type = UndefinedProperty; };
79 template<
class TypeTag>
80 struct Grid<TypeTag, TTag::ReservoirBaseProblem> {
using type = Dune::YaspGrid<2>; };
83 template<
class TypeTag>
87 template<
class TypeTag>
88 struct MaterialLaw<TypeTag, TTag::ReservoirBaseProblem>
91 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
92 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
95 ThreePhaseMaterialTraits<Scalar,
96 FluidSystem::waterPhaseIdx,
97 FluidSystem::oilPhaseIdx,
98 FluidSystem::gasPhaseIdx>;
101 using type = Opm::LinearMaterial<Traits>;
105 template<
class TypeTag>
106 struct NewtonWriteConvergence<TypeTag, TTag::ReservoirBaseProblem> {
static constexpr
bool value =
false; };
109 template<
class TypeTag>
110 struct EnableGravity<TypeTag, TTag::ReservoirBaseProblem> {
static constexpr
bool value =
true; };
113 template<
class TypeTag>
114 struct EnableConstraints<TypeTag, TTag::ReservoirBaseProblem> {
static constexpr
bool value =
true; };
117 template<
class TypeTag>
118 struct MaxDepth<TypeTag, TTag::ReservoirBaseProblem>
120 using type = GetPropType<TypeTag, Scalar>;
121 static constexpr type value = 2500;
123 template<
class TypeTag>
126 using type = GetPropType<TypeTag, Scalar>;
127 static constexpr type value = 293.15;
134 template<
class TypeTag>
135 struct EndTime<TypeTag, TTag::ReservoirBaseProblem>
137 using type = GetPropType<TypeTag, Scalar>;
138 static constexpr type value = 1000.0*24*60*60;
142 template<
class TypeTag>
143 struct InitialTimeStepSize<TypeTag, TTag::ReservoirBaseProblem>
145 using type = GetPropType<TypeTag, Scalar>;
146 static constexpr type value = 100e3;
150 template<
class TypeTag>
153 using type = GetPropType<TypeTag, Scalar>;
154 static constexpr type value = 0.01;
165 template<
class TypeTag>
166 struct FluidSystem<TypeTag, TTag::ReservoirBaseProblem>
169 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
172 using type = Opm::BlackOilFluidSystem<Scalar>;
176 template<
class TypeTag>
177 struct GridFile<TypeTag, TTag::ReservoirBaseProblem> {
static constexpr
auto value =
"data/reservoir.dgf"; };
180 template<
class TypeTag>
181 struct NewtonTolerance<TypeTag, TTag::ReservoirBaseProblem>
183 using type = GetPropType<TypeTag, Scalar>;
184 static constexpr type value = 1e-6;
207 template <
class TypeTag>
210 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
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>;
218 enum { dim = GridView::dimension };
219 enum { dimWorld = GridView::dimensionworld };
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 };
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>;
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>;
247 using InitialFluidState = Opm::CompositionalFluidState<Scalar,
256 : ParentType(simulator)
264 ParentType::finishInit();
266 temperature_ = EWOMS_GET_PARAM(TypeTag, Scalar, Temperature);
267 maxDepth_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxDepth);
268 wellWidth_ = EWOMS_GET_PARAM(TypeTag, Scalar, WellWidth);
270 std::vector<std::pair<Scalar, Scalar> > Bo = {
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 }
281 std::vector<std::pair<Scalar, Scalar> > muo = {
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 }
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 }
303 std::vector<std::pair<Scalar, Scalar> > Bg = {
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 }
315 std::vector<std::pair<Scalar, Scalar> > mug = {
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 }
328 Scalar rhoRefO = 786.0;
329 Scalar rhoRefG = 0.97;
330 Scalar rhoRefW = 1037.0;
331 FluidSystem::initBegin(1);
332 FluidSystem::setEnableDissolvedGas(
true);
333 FluidSystem::setEnableVaporizedOil(
false);
334 FluidSystem::setReferenceDensities(rhoRefO, rhoRefW, rhoRefG, 0);
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(1);
340 dryGasPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
341 dryGasPvt.setGasFormationVolumeFactor(0, Bg);
342 dryGasPvt.setGasViscosity(0, mug);
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(1);
348 liveOilPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
349 liveOilPvt.setSaturatedOilGasDissolutionFactor(0, Rs);
350 liveOilPvt.setSaturatedOilFormationVolumeFactor(0, Bo);
351 liveOilPvt.setSaturatedOilViscosity(0, muo);
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(1);
357 ccWaterPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
358 ccWaterPvt.setViscosity(0, 9.6e-4);
359 ccWaterPvt.setCompressibility(0, 1.450377e-10);
365 using GasPvtSharedPtr = std::shared_ptr<Opm::GasPvtMultiplexer<Scalar> >;
366 FluidSystem::setGasPvt(GasPvtSharedPtr(gasPvt));
368 using OilPvtSharedPtr = std::shared_ptr<Opm::OilPvtMultiplexer<Scalar> >;
369 FluidSystem::setOilPvt(OilPvtSharedPtr(oilPvt));
371 using WaterPvtSharedPtr = std::shared_ptr<Opm::WaterPvtMultiplexer<Scalar> >;
372 FluidSystem::setWaterPvt(WaterPvtSharedPtr(waterPvt));
374 FluidSystem::initEnd();
380 fineK_ = this->toDimMatrix_(1e-12);
381 coarseK_ = this->toDimMatrix_(1e-11);
385 coarsePorosity_ = 0.3;
387 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
388 fineMaterialParams_.setPcMinSat(phaseIdx, 0.0);
389 fineMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
391 coarseMaterialParams_.setPcMinSat(phaseIdx, 0.0);
392 coarseMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
396 fineMaterialParams_.finalize();
397 coarseMaterialParams_.finalize();
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(0);
406 for (
unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
407 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
408 const GlobalPosition& pos = elemCtx.pos(dofIdx, 0);
410 if (isFineMaterial_(pos))
411 materialParams_[globalDofIdx] = &fineMaterialParams_;
413 materialParams_[globalDofIdx] = &coarseMaterialParams_;
420 this->simulator().startNextEpisode(100.0*24*60*60);
428 ParentType::registerParameters();
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");
443 {
return std::string(
"reservoir_") + Model::name() +
"_" + Model::discretizationName(); }
453 this->simulator().startNextEpisode(1e100);
454 this->simulator().setTimeStepSize(5.0);
469 this->model().globalStorage(storage);
472 if (this->gridView().comm().rank() == 0) {
473 std::cout <<
"Storage: " << storage << std::endl << std::flush;
484 template <
class Context>
486 unsigned timeIdx)
const
488 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
489 if (isFineMaterial_(pos))
497 template <
class Context>
498 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
500 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
501 if (isFineMaterial_(pos))
502 return finePorosity_;
503 return coarsePorosity_;
509 template <
class Context>
511 unsigned spaceIdx,
unsigned timeIdx)
const
513 unsigned globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
514 return *materialParams_[globalIdx];
518 {
return *materialParams_[globalIdx]; }
534 template <
class Context>
536 unsigned spaceIdx OPM_UNUSED,
537 unsigned timeIdx OPM_UNUSED)
const
538 {
return temperature_; }
553 template <
class Context>
555 const Context& context OPM_UNUSED,
556 unsigned spaceIdx OPM_UNUSED,
557 unsigned timeIdx OPM_UNUSED)
const
576 template <
class Context>
578 const Context& context OPM_UNUSED,
579 unsigned spaceIdx OPM_UNUSED,
580 unsigned timeIdx OPM_UNUSED)
const
582 values.assignNaive(initialFluidState_);
585 for (
unsigned pvIdx = 0; pvIdx < values.size(); ++ pvIdx)
586 assert(std::isfinite(values[pvIdx]));
598 template <
class Context>
600 const Context& context,
602 unsigned timeIdx)
const
604 if (this->simulator().episodeIndex() == 1)
607 const auto& pos = context.pos(spaceIdx, timeIdx);
608 if (isInjector_(pos)) {
609 constraintValues.setActive(
true);
610 constraintValues.assignNaive(injectorFluidState_);
612 else if (isProducer_(pos)) {
613 constraintValues.setActive(
true);
614 constraintValues.assignNaive(producerFluidState_);
623 template <
class Context>
625 const Context& context OPM_UNUSED,
626 unsigned spaceIdx OPM_UNUSED,
627 unsigned timeIdx OPM_UNUSED)
const
628 { rate = Scalar(0.0); }
633 void initFluidState_()
635 auto& fs = initialFluidState_;
640 fs.setTemperature(temperature_);
645 fs.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
646 fs.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
647 fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
652 Scalar pw = pReservoir_;
655 const auto& matParams = fineMaterialParams_;
656 MaterialLaw::capillaryPressures(pC, matParams, fs);
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]));
663 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
664 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
665 fs.setMoleFraction(phaseIdx, compIdx, 0.0);
670 fs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
671 fs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
677 FluidSystem::saturatedDissolutionFactor(fs, oilPhaseIdx, 0);
678 Scalar XoGSat = FluidSystem::convertRsToXoG(RsSat, 0);
679 Scalar xoGSat = FluidSystem::convertXoGToxoG(XoGSat, 0);
680 Scalar xoG = 0.95*xoGSat;
681 Scalar xoO = 1.0 - xoG;
684 fs.setMoleFraction(oilPhaseIdx, gasCompIdx, xoG);
685 fs.setMoleFraction(oilPhaseIdx, oilCompIdx, xoO);
687 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
688 typename FluidSystem::template ParameterCache<Scalar> paramCache;
696 auto& injFs = injectorFluidState_;
697 injFs = initialFluidState_;
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);
708 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
709 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
710 injFs.setMoleFraction(phaseIdx, compIdx, 0.0);
712 injFs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
713 injFs.setMoleFraction(oilPhaseIdx, oilCompIdx, 1.0);
714 injFs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
723 auto& prodFs = producerFluidState_;
724 prodFs = initialFluidState_;
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);
741 bool isProducer_(
const GlobalPosition& pos)
const
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];
753 return width/2.0 - width*1e-5 < x && x < width/2.0 + width*(wellWidth_ + 1e-5);
756 bool isInjector_(
const GlobalPosition& pos)
const
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];
768 return x < width*wellWidth_ - width*1e-5 || x > width*(1.0 - wellWidth_) + width*1e-5;
771 bool isFineMaterial_(
const GlobalPosition& pos)
const
772 {
return pos[dim - 1] > layerBottom_; }
779 Scalar finePorosity_;
780 Scalar coarsePorosity_;
782 MaterialLawParams fineMaterialParams_;
783 MaterialLawParams coarseMaterialParams_;
784 std::vector<const MaterialLawParams*> materialParams_;
786 InitialFluidState initialFluidState_;
787 InitialFluidState injectorFluidState_;
788 InitialFluidState producerFluidState_;
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