28#ifndef EWOMS_RESERVOIR_PROBLEM_HH
29#define EWOMS_RESERVOIR_PROBLEM_HH
31#include <dune/common/version.hh>
32#include <dune/common/fvector.hh>
33#include <dune/common/fmatrix.hh>
35#include <dune/grid/yaspgrid.hh>
36#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
38#include <opm/material/fluidmatrixinteractions/LinearMaterial.hpp>
39#include <opm/material/fluidmatrixinteractions/MaterialTraits.hpp>
40#include <opm/material/fluidstates/CompositionalFluidState.hpp>
41#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
42#include <opm/material/constraintsolvers/ComputeFromReferencePhase.hpp>
43#include <opm/material/fluidsystems/blackoilpvt/DryGasPvt.hpp>
44#include <opm/material/fluidsystems/blackoilpvt/LiveOilPvt.hpp>
45#include <opm/material/fluidsystems/blackoilpvt/ConstantCompressibilityWaterPvt.hpp>
62template <
class TypeTag>
63class ReservoirProblem;
77template<
class TypeTag>
78struct Grid<TypeTag, TTag::ReservoirBaseProblem> {
using type = Dune::YaspGrid<2>; };
81template<
class TypeTag>
85template<
class TypeTag>
93 ThreePhaseMaterialTraits<Scalar,
94 FluidSystem::waterPhaseIdx,
95 FluidSystem::oilPhaseIdx,
96 FluidSystem::gasPhaseIdx>;
99 using type = Opm::LinearMaterial<Traits>;
103template<
class TypeTag>
104struct EnableConstraints<TypeTag, TTag::ReservoirBaseProblem> {
static constexpr bool value =
true; };
114template<
class TypeTag>
121 using type = Opm::BlackOilFluidSystem<Scalar>;
129template<
class Scalar>
130struct MaxDepth {
static constexpr Scalar
value = 2500.0; };
133template<
class Scalar>
134struct Temperature {
static constexpr Scalar
value = 293.15; };
137template<
class Scalar>
160template <
class TypeTag>
171 enum { dim = GridView::dimension };
172 enum { dimWorld = GridView::dimensionworld };
175 enum { numPhases = FluidSystem::numPhases };
176 enum { numComponents = FluidSystem::numComponents };
177 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
178 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
179 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
180 enum { gasCompIdx = FluidSystem::gasCompIdx };
181 enum { oilCompIdx = FluidSystem::oilCompIdx };
182 enum { waterCompIdx = FluidSystem::waterCompIdx };
195 using CoordScalar =
typename GridView::ctype;
196 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
197 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
198 using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
200 using InitialFluidState = Opm::CompositionalFluidState<Scalar,
209 : ParentType(simulator)
217 ParentType::finishInit();
219 temperature_ = Parameters::Get<Parameters::Temperature<Scalar>>();
220 maxDepth_ = Parameters::Get<Parameters::MaxDepth<Scalar>>();
221 wellWidth_ = Parameters::Get<Parameters::WellWidth<Scalar>>();
223 std::vector<std::pair<Scalar, Scalar> > Bo = {
225 { 1.82504e+06, 1.15 },
226 { 3.54873e+06, 1.207 },
227 { 6.99611e+06, 1.295 },
228 { 1.38909e+07, 1.435 },
229 { 1.73382e+07, 1.5 },
230 { 2.07856e+07, 1.565 },
231 { 2.76804e+07, 1.695 },
232 { 3.45751e+07, 1.827 }
234 std::vector<std::pair<Scalar, Scalar> > muo = {
236 { 1.82504e+06, 0.000975 },
237 { 3.54873e+06, 0.00091 },
238 { 6.99611e+06, 0.00083 },
239 { 1.38909e+07, 0.000695 },
240 { 1.73382e+07, 0.000641 },
241 { 2.07856e+07, 0.000594 },
242 { 2.76804e+07, 0.00051 },
243 { 3.45751e+07, 0.000449 }
245 std::vector<std::pair<Scalar, Scalar> > Rs = {
246 { 101353, 0.178108 },
247 { 1.82504e+06, 16.1187 },
248 { 3.54873e+06, 32.0594 },
249 { 6.99611e+06, 66.0779 },
250 { 1.38909e+07, 113.276 },
251 { 1.73382e+07, 138.033 },
252 { 2.07856e+07, 165.64 },
253 { 2.76804e+07, 226.197 },
254 { 3.45751e+07, 288.178 }
256 std::vector<std::pair<Scalar, Scalar> > Bg = {
258 { 1.82504e+06, 0.0678972 },
259 { 3.54873e+06, 0.0352259 },
260 { 6.99611e+06, 0.0179498 },
261 { 1.38909e+07, 0.00906194 },
262 { 1.73382e+07, 0.00726527 },
263 { 2.07856e+07, 0.00606375 },
264 { 2.76804e+07, 0.00455343 },
265 { 3.45751e+07, 0.00364386 },
266 { 6.21542e+07, 0.00216723 }
268 std::vector<std::pair<Scalar, Scalar> > mug = {
270 { 1.82504e+06, 9.6e-06 },
271 { 3.54873e+06, 1.12e-05 },
272 { 6.99611e+06, 1.4e-05 },
273 { 1.38909e+07, 1.89e-05 },
274 { 1.73382e+07, 2.08e-05 },
275 { 2.07856e+07, 2.28e-05 },
276 { 2.76804e+07, 2.68e-05 },
277 { 3.45751e+07, 3.09e-05 },
278 { 6.21542e+07, 4.7e-05 }
281 Scalar rhoRefO = 786.0;
282 Scalar rhoRefG = 0.97;
283 Scalar rhoRefW = 1037.0;
284 FluidSystem::initBegin(1);
285 FluidSystem::setEnableDissolvedGas(
true);
286 FluidSystem::setEnableVaporizedOil(
false);
287 FluidSystem::setReferenceDensities(rhoRefO, rhoRefW, rhoRefG, 0);
289 Opm::GasPvtMultiplexer<Scalar> *gasPvt =
new Opm::GasPvtMultiplexer<Scalar>;
290 gasPvt->setApproach(GasPvtApproach::DryGas);
291 auto& dryGasPvt = gasPvt->template getRealPvt<GasPvtApproach::DryGas>();
292 dryGasPvt.setNumRegions(1);
293 dryGasPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
294 dryGasPvt.setGasFormationVolumeFactor(0, Bg);
295 dryGasPvt.setGasViscosity(0, mug);
297 Opm::OilPvtMultiplexer<Scalar> *oilPvt =
new Opm::OilPvtMultiplexer<Scalar>;
298 oilPvt->setApproach(OilPvtApproach::LiveOil);
299 auto& liveOilPvt = oilPvt->template getRealPvt<OilPvtApproach::LiveOil>();
300 liveOilPvt.setNumRegions(1);
301 liveOilPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
302 liveOilPvt.setSaturatedOilGasDissolutionFactor(0, Rs);
303 liveOilPvt.setSaturatedOilFormationVolumeFactor(0, Bo);
304 liveOilPvt.setSaturatedOilViscosity(0, muo);
306 Opm::WaterPvtMultiplexer<Scalar> *waterPvt =
new Opm::WaterPvtMultiplexer<Scalar>;
307 waterPvt->setApproach(WaterPvtApproach::ConstantCompressibilityWater);
308 auto& ccWaterPvt = waterPvt->template getRealPvt<WaterPvtApproach::ConstantCompressibilityWater>();
309 ccWaterPvt.setNumRegions(1);
310 ccWaterPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
311 ccWaterPvt.setViscosity(0, 9.6e-4);
312 ccWaterPvt.setCompressibility(0, 1.450377e-10);
318 using GasPvtSharedPtr = std::shared_ptr<Opm::GasPvtMultiplexer<Scalar> >;
319 FluidSystem::setGasPvt(GasPvtSharedPtr(gasPvt));
321 using OilPvtSharedPtr = std::shared_ptr<Opm::OilPvtMultiplexer<Scalar> >;
322 FluidSystem::setOilPvt(OilPvtSharedPtr(oilPvt));
324 using WaterPvtSharedPtr = std::shared_ptr<Opm::WaterPvtMultiplexer<Scalar> >;
325 FluidSystem::setWaterPvt(WaterPvtSharedPtr(waterPvt));
327 FluidSystem::initEnd();
333 fineK_ = this->toDimMatrix_(1e-12);
334 coarseK_ = this->toDimMatrix_(1e-11);
338 coarsePorosity_ = 0.3;
340 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
341 fineMaterialParams_.setPcMinSat(phaseIdx, 0.0);
342 fineMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
344 coarseMaterialParams_.setPcMinSat(phaseIdx, 0.0);
345 coarseMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
349 fineMaterialParams_.finalize();
350 coarseMaterialParams_.finalize();
352 materialParams_.resize(this->model().numGridDof());
353 ElementContext elemCtx(this->simulator());
354 auto eIt = this->simulator().gridView().template begin<0>();
355 const auto& eEndIt = this->simulator().gridView().template end<0>();
356 for (; eIt != eEndIt; ++eIt) {
357 elemCtx.updateStencil(*eIt);
358 size_t nDof = elemCtx.numPrimaryDof(0);
359 for (
unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
360 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
361 const GlobalPosition& pos = elemCtx.pos(dofIdx, 0);
363 if (isFineMaterial_(pos))
364 materialParams_[globalDofIdx] = &fineMaterialParams_;
366 materialParams_[globalDofIdx] = &coarseMaterialParams_;
373 this->simulator().startNextEpisode(100.0*24*60*60);
381 ParentType::registerParameters();
383 Parameters::Register<Parameters::Temperature<Scalar>>
384 (
"The temperature [K] in the reservoir");
385 Parameters::Register<Parameters::MaxDepth<Scalar>>
386 (
"The maximum depth [m] of the reservoir");
387 Parameters::Register<Parameters::WellWidth<Scalar>>
388 (
"The width of producer/injector wells as a fraction of the width"
389 " of the spatial domain");
391 Parameters::SetDefault<Parameters::GridFile>(
"data/reservoir.dgf");
395 Parameters::SetDefault<Parameters::EndTime<Scalar>>(1000.0*24*60*60);
397 Parameters::SetDefault<Parameters::EnableStorageCache>(
true);
398 Parameters::SetDefault<Parameters::GridFile>(
"data/reservoir.dgf");
399 Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(100e3);
401 Parameters::SetDefault<Parameters::NewtonTolerance<Scalar>>(1e-6);
403 Parameters::SetDefault<Parameters::EnableGravity>(
true);
410 {
return std::string(
"reservoir_") + Model::name() +
"_" + Model::discretizationName(); }
420 this->simulator().startNextEpisode(1e100);
421 this->simulator().setTimeStepSize(5.0);
436 this->model().globalStorage(storage);
439 if (this->gridView().comm().rank() == 0) {
440 std::cout <<
"Storage: " << storage << std::endl << std::flush;
451 template <
class Context>
453 unsigned timeIdx)
const
455 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
456 if (isFineMaterial_(pos))
464 template <
class Context>
465 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
467 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
468 if (isFineMaterial_(pos))
469 return finePorosity_;
470 return coarsePorosity_;
476 template <
class Context>
478 unsigned spaceIdx,
unsigned timeIdx)
const
480 unsigned globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
481 return *materialParams_[globalIdx];
485 {
return *materialParams_[globalIdx]; }
501 template <
class Context>
505 {
return temperature_; }
520 template <
class Context>
543 template <
class Context>
549 values.assignNaive(initialFluidState_);
552 for (
unsigned pvIdx = 0; pvIdx < values.size(); ++ pvIdx)
553 assert(std::isfinite(values[pvIdx]));
565 template <
class Context>
567 const Context& context,
569 unsigned timeIdx)
const
571 if (this->simulator().episodeIndex() == 1)
574 const auto& pos = context.pos(spaceIdx, timeIdx);
575 if (isInjector_(pos)) {
576 constraintValues.setActive(
true);
577 constraintValues.assignNaive(injectorFluidState_);
579 else if (isProducer_(pos)) {
580 constraintValues.setActive(
true);
581 constraintValues.assignNaive(producerFluidState_);
590 template <
class Context>
595 { rate = Scalar(0.0); }
600 void initFluidState_()
602 auto& fs = initialFluidState_;
607 fs.setTemperature(temperature_);
612 fs.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
613 fs.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
614 fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
619 Scalar pw = pReservoir_;
622 const auto& matParams = fineMaterialParams_;
623 MaterialLaw::capillaryPressures(pC, matParams, fs);
625 fs.setPressure(oilPhaseIdx, pw + (pC[oilPhaseIdx] - pC[waterPhaseIdx]));
626 fs.setPressure(waterPhaseIdx, pw + (pC[waterPhaseIdx] - pC[waterPhaseIdx]));
627 fs.setPressure(gasPhaseIdx, pw + (pC[gasPhaseIdx] - pC[waterPhaseIdx]));
630 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
631 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
632 fs.setMoleFraction(phaseIdx, compIdx, 0.0);
637 fs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
638 fs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
644 FluidSystem::saturatedDissolutionFactor(fs, oilPhaseIdx, 0);
645 Scalar XoGSat = FluidSystem::convertRsToXoG(RsSat, 0);
646 Scalar xoGSat = FluidSystem::convertXoGToxoG(XoGSat, 0);
647 Scalar xoG = 0.95*xoGSat;
648 Scalar xoO = 1.0 - xoG;
651 fs.setMoleFraction(oilPhaseIdx, gasCompIdx, xoG);
652 fs.setMoleFraction(oilPhaseIdx, oilCompIdx, xoO);
654 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
655 typename FluidSystem::template ParameterCache<Scalar> paramCache;
663 auto& injFs = injectorFluidState_;
664 injFs = initialFluidState_;
666 Scalar pInj = pReservoir_ * 1.5;
667 injFs.setPressure(waterPhaseIdx, pInj);
668 injFs.setPressure(oilPhaseIdx, pInj);
669 injFs.setPressure(gasPhaseIdx, pInj);
670 injFs.setSaturation(waterPhaseIdx, 1.0);
671 injFs.setSaturation(oilPhaseIdx, 0.0);
672 injFs.setSaturation(gasPhaseIdx, 0.0);
675 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
676 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
677 injFs.setMoleFraction(phaseIdx, compIdx, 0.0);
679 injFs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
680 injFs.setMoleFraction(oilPhaseIdx, oilCompIdx, 1.0);
681 injFs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
690 auto& prodFs = producerFluidState_;
691 prodFs = initialFluidState_;
693 Scalar pProd = pReservoir_ / 1.5;
694 prodFs.setPressure(waterPhaseIdx, pProd);
695 prodFs.setPressure(oilPhaseIdx, pProd);
696 prodFs.setPressure(gasPhaseIdx, pProd);
697 prodFs.setSaturation(waterPhaseIdx, 0.0);
698 prodFs.setSaturation(oilPhaseIdx, 1.0);
699 prodFs.setSaturation(gasPhaseIdx, 0.0);
708 bool isProducer_(
const GlobalPosition& pos)
const
710 Scalar x = pos[0] - this->boundingBoxMin()[0];
711 Scalar y = pos[dim - 1] - this->boundingBoxMin()[dim - 1];
712 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
713 Scalar height = this->boundingBoxMax()[dim - 1] - this->boundingBoxMin()[dim - 1];
720 return width/2.0 - width*1e-5 < x && x < width/2.0 + width*(wellWidth_ + 1e-5);
723 bool isInjector_(
const GlobalPosition& pos)
const
725 Scalar x = pos[0] - this->boundingBoxMin()[0];
726 Scalar y = pos[dim - 1] - this->boundingBoxMin()[dim - 1];
727 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
728 Scalar height = this->boundingBoxMax()[dim - 1] - this->boundingBoxMin()[dim - 1];
735 return x < width*wellWidth_ - width*1e-5 || x > width*(1.0 - wellWidth_) + width*1e-5;
738 bool isFineMaterial_(
const GlobalPosition& pos)
const
739 {
return pos[dim - 1] > layerBottom_; }
746 Scalar finePorosity_;
747 Scalar coarsePorosity_;
749 MaterialLawParams fineMaterialParams_;
750 MaterialLawParams coarseMaterialParams_;
751 std::vector<const MaterialLawParams*> materialParams_;
753 InitialFluidState initialFluidState_;
754 InitialFluidState injectorFluidState_;
755 InitialFluidState producerFluidState_;
Defines a type tags and some fundamental properties all models.
Declares the properties required by the black oil model.
Some simple test problem for the black-oil VCVF discretization inspired by an oil reservoir.
Definition: reservoirproblem.hh:162
static void registerParameters()
Definition: reservoirproblem.hh:379
void source(RateVector &rate, const Context &, unsigned, unsigned) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: reservoirproblem.hh:591
void endTimeStep()
Called by the simulator after each time integration.
Definition: reservoirproblem.hh:427
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:465
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Evaluate the initial value for a control volume.
Definition: reservoirproblem.hh:544
ReservoirProblem(Simulator &simulator)
Definition: reservoirproblem.hh:208
void constraints(Constraints &constraintValues, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the constraints for a control volume.
Definition: reservoirproblem.hh:566
void boundary(BoundaryRateVector &values, const Context &, unsigned, unsigned) const
Evaluate the boundary conditions for a boundary segment.
Definition: reservoirproblem.hh:521
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:477
std::string name() const
The problem name.
Definition: reservoirproblem.hh:409
const MaterialLawParams & materialLawParams(unsigned globalIdx) const
Definition: reservoirproblem.hh:484
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: reservoirproblem.hh:215
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:452
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: reservoirproblem.hh:502
void endEpisode()
Called when the end of an simulation episode is reached.
Definition: reservoirproblem.hh:415
Declare the properties used by the infrastructure code of the finite volume discretizations.
Declare the properties used by the infrastructure code of the finite volume discretizations.
Defines the common parameters for the porous medium multi-phase models.
Definition: blackoilnewtonmethodparams.hpp:31
Definition: blackoilmodel.hh:72
Definition: blackoilboundaryratevector.hh:37
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:235
static constexpr Scalar value
Definition: co2injectionproblem.hh:160
static constexpr Scalar value
Definition: co2injectionproblem.hh:165
Definition: reservoirproblem.hh:138
static constexpr Scalar value
Definition: reservoirproblem.hh:138
Specify whether the some degrees of fredom can be constraint.
Definition: fvbaseproperties.hh:199
Opm::BlackOilFluidSystem< Scalar > type
Definition: reservoirproblem.hh:121
The fluid systems including the information about the phases.
Definition: multiphasebaseproperties.hh:69
Dune::YaspGrid< 2 > type
Definition: reservoirproblem.hh:78
The type of the DUNE grid.
Definition: basicproperties.hh:100
Opm::LinearMaterial< Traits > type
Definition: reservoirproblem.hh:99
The material law which ought to be used (extracted from the spatial parameters)
Definition: multiphasebaseproperties.hh:51
The type of the problem.
Definition: fvbaseproperties.hh:81
Definition: reservoirproblem.hh:72