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>
42#include <dune/grid/yaspgrid.hh>
43#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
45#include <dune/common/version.hh>
46#include <dune/common/fvector.hh>
47#include <dune/common/fmatrix.hh>
53template <
class TypeTag>
54class ReservoirProblem;
58namespace Opm::Properties {
68template<
class TypeTag,
class MyTypeTag>
69struct MaxDepth {
using type = UndefinedProperty; };
71template<
class TypeTag,
class MyTypeTag>
72struct Temperature {
using type = UndefinedProperty; };
74template<
class TypeTag,
class MyTypeTag>
75struct WellWidth {
using type = UndefinedProperty; };
78template<
class TypeTag>
79struct Grid<TypeTag, TTag::ReservoirBaseProblem> {
using type = Dune::YaspGrid<2>; };
82template<
class TypeTag>
86template<
class TypeTag>
87struct MaterialLaw<TypeTag, TTag::ReservoirBaseProblem>
90 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
91 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
94 ThreePhaseMaterialTraits<Scalar,
95 FluidSystem::waterPhaseIdx,
96 FluidSystem::oilPhaseIdx,
97 FluidSystem::gasPhaseIdx>;
100 using type = Opm::LinearMaterial<Traits>;
104template<
class TypeTag>
105struct NewtonWriteConvergence<TypeTag, TTag::ReservoirBaseProblem> {
static constexpr bool value =
false; };
108template<
class TypeTag>
109struct EnableGravity<TypeTag, TTag::ReservoirBaseProblem> {
static constexpr bool value =
true; };
112template<
class TypeTag>
113struct EnableConstraints<TypeTag, TTag::ReservoirBaseProblem> {
static constexpr bool value =
true; };
116template<
class TypeTag>
117struct MaxDepth<TypeTag, TTag::ReservoirBaseProblem>
119 using type = GetPropType<TypeTag, Scalar>;
120 static constexpr type value = 2500;
122template<
class TypeTag>
125 using type = GetPropType<TypeTag, Scalar>;
126 static constexpr type value = 293.15;
133template<
class TypeTag>
134struct EndTime<TypeTag, TTag::ReservoirBaseProblem>
136 using type = GetPropType<TypeTag, Scalar>;
137 static constexpr type value = 1000.0*24*60*60;
141template<
class TypeTag>
142struct InitialTimeStepSize<TypeTag, TTag::ReservoirBaseProblem>
144 using type = GetPropType<TypeTag, Scalar>;
145 static constexpr type value = 100e3;
149template<
class TypeTag>
152 using type = GetPropType<TypeTag, Scalar>;
153 static constexpr type value = 0.01;
164template<
class TypeTag>
165struct FluidSystem<TypeTag, TTag::ReservoirBaseProblem>
168 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
171 using type = Opm::BlackOilFluidSystem<Scalar>;
175template<
class TypeTag>
176struct GridFile<TypeTag, TTag::ReservoirBaseProblem> {
static constexpr auto value =
"data/reservoir.dgf"; };
179template<
class TypeTag>
180struct NewtonTolerance<TypeTag, TTag::ReservoirBaseProblem>
182 using type = GetPropType<TypeTag, Scalar>;
183 static constexpr type value = 1e-6;
206template <
class TypeTag>
209 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
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>;
217 enum { dim = GridView::dimension };
218 enum { dimWorld = GridView::dimensionworld };
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 };
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>;
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>;
246 using InitialFluidState = Opm::CompositionalFluidState<Scalar,
255 : ParentType(simulator)
263 ParentType::finishInit();
265 temperature_ = EWOMS_GET_PARAM(TypeTag, Scalar, Temperature);
266 maxDepth_ = EWOMS_GET_PARAM(TypeTag, Scalar, MaxDepth);
267 wellWidth_ = EWOMS_GET_PARAM(TypeTag, Scalar, WellWidth);
269 std::vector<std::pair<Scalar, Scalar> > Bo = {
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 }
280 std::vector<std::pair<Scalar, Scalar> > muo = {
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 }
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 }
302 std::vector<std::pair<Scalar, Scalar> > Bg = {
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 }
314 std::vector<std::pair<Scalar, Scalar> > mug = {
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 }
327 Scalar rhoRefO = 786.0;
328 Scalar rhoRefG = 0.97;
329 Scalar rhoRefW = 1037.0;
330 FluidSystem::initBegin(1);
331 FluidSystem::setEnableDissolvedGas(
true);
332 FluidSystem::setEnableVaporizedOil(
false);
333 FluidSystem::setReferenceDensities(rhoRefO, rhoRefW, rhoRefG, 0);
335 Opm::GasPvtMultiplexer<Scalar> *gasPvt =
new Opm::GasPvtMultiplexer<Scalar>;
336 gasPvt->setApproach(GasPvtApproach::DryGas);
337 auto& dryGasPvt = gasPvt->template getRealPvt<GasPvtApproach::DryGas>();
338 dryGasPvt.setNumRegions(1);
339 dryGasPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
340 dryGasPvt.setGasFormationVolumeFactor(0, Bg);
341 dryGasPvt.setGasViscosity(0, mug);
343 Opm::OilPvtMultiplexer<Scalar> *oilPvt =
new Opm::OilPvtMultiplexer<Scalar>;
344 oilPvt->setApproach(OilPvtApproach::LiveOil);
345 auto& liveOilPvt = oilPvt->template getRealPvt<OilPvtApproach::LiveOil>();
346 liveOilPvt.setNumRegions(1);
347 liveOilPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
348 liveOilPvt.setSaturatedOilGasDissolutionFactor(0, Rs);
349 liveOilPvt.setSaturatedOilFormationVolumeFactor(0, Bo);
350 liveOilPvt.setSaturatedOilViscosity(0, muo);
352 Opm::WaterPvtMultiplexer<Scalar> *waterPvt =
new Opm::WaterPvtMultiplexer<Scalar>;
353 waterPvt->setApproach(WaterPvtApproach::ConstantCompressibilityWater);
354 auto& ccWaterPvt = waterPvt->template getRealPvt<WaterPvtApproach::ConstantCompressibilityWater>();
355 ccWaterPvt.setNumRegions(1);
356 ccWaterPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
357 ccWaterPvt.setViscosity(0, 9.6e-4);
358 ccWaterPvt.setCompressibility(0, 1.450377e-10);
364 using GasPvtSharedPtr = std::shared_ptr<Opm::GasPvtMultiplexer<Scalar> >;
365 FluidSystem::setGasPvt(GasPvtSharedPtr(gasPvt));
367 using OilPvtSharedPtr = std::shared_ptr<Opm::OilPvtMultiplexer<Scalar> >;
368 FluidSystem::setOilPvt(OilPvtSharedPtr(oilPvt));
370 using WaterPvtSharedPtr = std::shared_ptr<Opm::WaterPvtMultiplexer<Scalar> >;
371 FluidSystem::setWaterPvt(WaterPvtSharedPtr(waterPvt));
373 FluidSystem::initEnd();
379 fineK_ = this->toDimMatrix_(1e-12);
380 coarseK_ = this->toDimMatrix_(1e-11);
384 coarsePorosity_ = 0.3;
386 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
387 fineMaterialParams_.setPcMinSat(phaseIdx, 0.0);
388 fineMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
390 coarseMaterialParams_.setPcMinSat(phaseIdx, 0.0);
391 coarseMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
395 fineMaterialParams_.finalize();
396 coarseMaterialParams_.finalize();
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(0);
405 for (
unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
406 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
407 const GlobalPosition& pos = elemCtx.pos(dofIdx, 0);
409 if (isFineMaterial_(pos))
410 materialParams_[globalDofIdx] = &fineMaterialParams_;
412 materialParams_[globalDofIdx] = &coarseMaterialParams_;
419 this->simulator().startNextEpisode(100.0*24*60*60);
427 ParentType::registerParameters();
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");
442 {
return std::string(
"reservoir_") + Model::name() +
"_" + Model::discretizationName(); }
452 this->simulator().startNextEpisode(1e100);
453 this->simulator().setTimeStepSize(5.0);
468 this->model().globalStorage(storage);
471 if (this->gridView().comm().rank() == 0) {
472 std::cout <<
"Storage: " << storage << std::endl << std::flush;
483 template <
class Context>
485 unsigned timeIdx)
const
487 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
488 if (isFineMaterial_(pos))
496 template <
class Context>
497 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
499 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
500 if (isFineMaterial_(pos))
501 return finePorosity_;
502 return coarsePorosity_;
508 template <
class Context>
510 unsigned spaceIdx,
unsigned timeIdx)
const
512 unsigned globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
513 return *materialParams_[globalIdx];
517 {
return *materialParams_[globalIdx]; }
533 template <
class Context>
537 {
return temperature_; }
552 template <
class Context>
575 template <
class Context>
581 values.assignNaive(initialFluidState_);
584 for (
unsigned pvIdx = 0; pvIdx < values.size(); ++ pvIdx)
585 assert(std::isfinite(values[pvIdx]));
597 template <
class Context>
599 const Context& context,
601 unsigned timeIdx)
const
603 if (this->simulator().episodeIndex() == 1)
606 const auto& pos = context.pos(spaceIdx, timeIdx);
607 if (isInjector_(pos)) {
608 constraintValues.setActive(
true);
609 constraintValues.assignNaive(injectorFluidState_);
611 else if (isProducer_(pos)) {
612 constraintValues.setActive(
true);
613 constraintValues.assignNaive(producerFluidState_);
622 template <
class Context>
627 { rate = Scalar(0.0); }
632 void initFluidState_()
634 auto& fs = initialFluidState_;
639 fs.setTemperature(temperature_);
644 fs.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
645 fs.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
646 fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
651 Scalar pw = pReservoir_;
654 const auto& matParams = fineMaterialParams_;
655 MaterialLaw::capillaryPressures(pC, matParams, fs);
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]));
662 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
663 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
664 fs.setMoleFraction(phaseIdx, compIdx, 0.0);
669 fs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
670 fs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
676 FluidSystem::saturatedDissolutionFactor(fs, oilPhaseIdx, 0);
677 Scalar XoGSat = FluidSystem::convertRsToXoG(RsSat, 0);
678 Scalar xoGSat = FluidSystem::convertXoGToxoG(XoGSat, 0);
679 Scalar xoG = 0.95*xoGSat;
680 Scalar xoO = 1.0 - xoG;
683 fs.setMoleFraction(oilPhaseIdx, gasCompIdx, xoG);
684 fs.setMoleFraction(oilPhaseIdx, oilCompIdx, xoO);
686 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
687 typename FluidSystem::template ParameterCache<Scalar> paramCache;
695 auto& injFs = injectorFluidState_;
696 injFs = initialFluidState_;
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);
707 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
708 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
709 injFs.setMoleFraction(phaseIdx, compIdx, 0.0);
711 injFs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
712 injFs.setMoleFraction(oilPhaseIdx, oilCompIdx, 1.0);
713 injFs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
722 auto& prodFs = producerFluidState_;
723 prodFs = initialFluidState_;
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);
740 bool isProducer_(
const GlobalPosition& pos)
const
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];
752 return width/2.0 - width*1e-5 < x && x < width/2.0 + width*(wellWidth_ + 1e-5);
755 bool isInjector_(
const GlobalPosition& pos)
const
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];
767 return x < width*wellWidth_ - width*1e-5 || x > width*(1.0 - wellWidth_) + width*1e-5;
770 bool isFineMaterial_(
const GlobalPosition& pos)
const
771 {
return pos[dim - 1] > layerBottom_; }
778 Scalar finePorosity_;
779 Scalar coarsePorosity_;
781 MaterialLawParams fineMaterialParams_;
782 MaterialLawParams coarseMaterialParams_;
783 std::vector<const MaterialLawParams*> materialParams_;
785 InitialFluidState initialFluidState_;
786 InitialFluidState injectorFluidState_;
787 InitialFluidState producerFluidState_;
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:88
Definition: reservoirproblem.hh:63
Definition: co2injectionproblem.hh:90
Definition: reservoirproblem.hh:75