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,
101 using type = Opm::LinearMaterial<Traits>;
105template<
class TypeTag>
106struct EnableConstraints<TypeTag, TTag::ReservoirBaseProblem> {
static constexpr bool value =
true; };
116template<
class TypeTag>
123 using type = Opm::BlackOilFluidSystem<Scalar>;
131template<
class Scalar>
132struct MaxDepth {
static constexpr Scalar
value = 2500.0; };
135template<
class Scalar>
136struct Temperature {
static constexpr Scalar
value = 293.15; };
139template<
class Scalar>
162template <
class TypeTag>
173 enum { dim = GridView::dimension };
174 enum { dimWorld = GridView::dimensionworld };
177 enum { numPhases = FluidSystem::numPhases };
178 enum { numComponents = FluidSystem::numComponents };
179 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
180 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
181 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
182 enum { gasCompIdx = FluidSystem::gasCompIdx };
183 enum { oilCompIdx = FluidSystem::oilCompIdx };
184 enum { waterCompIdx = FluidSystem::waterCompIdx };
197 using CoordScalar =
typename GridView::ctype;
198 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
199 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
200 using PhaseVector = Dune::FieldVector<Scalar, numPhases>;
202 using InitialFluidState = Opm::CompositionalFluidState<Scalar,
211 : ParentType(simulator)
219 ParentType::finishInit();
221 temperature_ = Parameters::Get<Parameters::Temperature<Scalar>>();
222 maxDepth_ = Parameters::Get<Parameters::MaxDepth<Scalar>>();
223 wellWidth_ = Parameters::Get<Parameters::WellWidth<Scalar>>();
225 std::vector<std::pair<Scalar, Scalar> > Bo = {
227 { 1.82504e+06, 1.15 },
228 { 3.54873e+06, 1.207 },
229 { 6.99611e+06, 1.295 },
230 { 1.38909e+07, 1.435 },
231 { 1.73382e+07, 1.5 },
232 { 2.07856e+07, 1.565 },
233 { 2.76804e+07, 1.695 },
234 { 3.45751e+07, 1.827 }
236 std::vector<std::pair<Scalar, Scalar> > muo = {
238 { 1.82504e+06, 0.000975 },
239 { 3.54873e+06, 0.00091 },
240 { 6.99611e+06, 0.00083 },
241 { 1.38909e+07, 0.000695 },
242 { 1.73382e+07, 0.000641 },
243 { 2.07856e+07, 0.000594 },
244 { 2.76804e+07, 0.00051 },
245 { 3.45751e+07, 0.000449 }
247 std::vector<std::pair<Scalar, Scalar> > Rs = {
248 { 101353, 0.178108 },
249 { 1.82504e+06, 16.1187 },
250 { 3.54873e+06, 32.0594 },
251 { 6.99611e+06, 66.0779 },
252 { 1.38909e+07, 113.276 },
253 { 1.73382e+07, 138.033 },
254 { 2.07856e+07, 165.64 },
255 { 2.76804e+07, 226.197 },
256 { 3.45751e+07, 288.178 }
258 std::vector<std::pair<Scalar, Scalar> > Bg = {
260 { 1.82504e+06, 0.0678972 },
261 { 3.54873e+06, 0.0352259 },
262 { 6.99611e+06, 0.0179498 },
263 { 1.38909e+07, 0.00906194 },
264 { 1.73382e+07, 0.00726527 },
265 { 2.07856e+07, 0.00606375 },
266 { 2.76804e+07, 0.00455343 },
267 { 3.45751e+07, 0.00364386 },
268 { 6.21542e+07, 0.00216723 }
270 std::vector<std::pair<Scalar, Scalar> > mug = {
272 { 1.82504e+06, 9.6e-06 },
273 { 3.54873e+06, 1.12e-05 },
274 { 6.99611e+06, 1.4e-05 },
275 { 1.38909e+07, 1.89e-05 },
276 { 1.73382e+07, 2.08e-05 },
277 { 2.07856e+07, 2.28e-05 },
278 { 2.76804e+07, 2.68e-05 },
279 { 3.45751e+07, 3.09e-05 },
280 { 6.21542e+07, 4.7e-05 }
283 Scalar rhoRefO = 786.0;
284 Scalar rhoRefG = 0.97;
285 Scalar rhoRefW = 1037.0;
286 FluidSystem::initBegin(1);
287 FluidSystem::setEnableDissolvedGas(
true);
288 FluidSystem::setEnableVaporizedOil(
false);
289 FluidSystem::setReferenceDensities(rhoRefO, rhoRefW, rhoRefG, 0);
291 Opm::GasPvtMultiplexer<Scalar> *gasPvt =
new Opm::GasPvtMultiplexer<Scalar>;
292 gasPvt->setApproach(GasPvtApproach::DryGas);
293 auto& dryGasPvt = gasPvt->template getRealPvt<GasPvtApproach::DryGas>();
294 dryGasPvt.setNumRegions(1);
295 dryGasPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
296 dryGasPvt.setGasFormationVolumeFactor(0, Bg);
297 dryGasPvt.setGasViscosity(0, Opm::Tabulated1DFunction<Scalar>(mug));
299 Opm::OilPvtMultiplexer<Scalar> *oilPvt =
new Opm::OilPvtMultiplexer<Scalar>;
300 oilPvt->setApproach(OilPvtApproach::LiveOil);
301 auto& liveOilPvt = oilPvt->template getRealPvt<OilPvtApproach::LiveOil>();
302 liveOilPvt.setNumRegions(1);
303 liveOilPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
304 liveOilPvt.setSaturatedOilGasDissolutionFactor(0, Rs);
305 liveOilPvt.setSaturatedOilFormationVolumeFactor(0, Bo);
306 liveOilPvt.setSaturatedOilViscosity(0, muo);
308 Opm::WaterPvtMultiplexer<Scalar> *waterPvt =
new Opm::WaterPvtMultiplexer<Scalar>;
309 waterPvt->setApproach(WaterPvtApproach::ConstantCompressibilityWater);
310 auto& ccWaterPvt = waterPvt->template getRealPvt<WaterPvtApproach::ConstantCompressibilityWater>();
311 ccWaterPvt.setNumRegions(1);
312 ccWaterPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
313 ccWaterPvt.setViscosity(0, 9.6e-4);
314 ccWaterPvt.setCompressibility(0, 1.450377e-10);
320 using GasPvtSharedPtr = std::shared_ptr<Opm::GasPvtMultiplexer<Scalar> >;
321 FluidSystem::setGasPvt(GasPvtSharedPtr(gasPvt));
323 using OilPvtSharedPtr = std::shared_ptr<Opm::OilPvtMultiplexer<Scalar> >;
324 FluidSystem::setOilPvt(OilPvtSharedPtr(oilPvt));
326 using WaterPvtSharedPtr = std::shared_ptr<Opm::WaterPvtMultiplexer<Scalar> >;
327 FluidSystem::setWaterPvt(WaterPvtSharedPtr(waterPvt));
329 FluidSystem::initEnd();
335 fineK_ = this->toDimMatrix_(1e-12);
336 coarseK_ = this->toDimMatrix_(1e-11);
340 coarsePorosity_ = 0.3;
342 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
343 fineMaterialParams_.setPcMinSat(phaseIdx, 0.0);
344 fineMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
346 coarseMaterialParams_.setPcMinSat(phaseIdx, 0.0);
347 coarseMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
351 fineMaterialParams_.finalize();
352 coarseMaterialParams_.finalize();
354 materialParams_.resize(this->model().numGridDof());
355 ElementContext elemCtx(this->simulator());
356 auto eIt = this->simulator().gridView().template begin<0>();
357 const auto& eEndIt = this->simulator().gridView().template end<0>();
358 for (; eIt != eEndIt; ++eIt) {
359 elemCtx.updateStencil(*eIt);
360 size_t nDof = elemCtx.numPrimaryDof(0);
361 for (
unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
362 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
363 const GlobalPosition& pos = elemCtx.pos(dofIdx, 0);
365 if (isFineMaterial_(pos))
366 materialParams_[globalDofIdx] = &fineMaterialParams_;
368 materialParams_[globalDofIdx] = &coarseMaterialParams_;
375 this->simulator().startNextEpisode(100.0*24*60*60);
383 ParentType::registerParameters();
385 Parameters::Register<Parameters::Temperature<Scalar>>
386 (
"The temperature [K] in the reservoir");
387 Parameters::Register<Parameters::MaxDepth<Scalar>>
388 (
"The maximum depth [m] of the reservoir");
389 Parameters::Register<Parameters::WellWidth<Scalar>>
390 (
"The width of producer/injector wells as a fraction of the width"
391 " of the spatial domain");
393 Parameters::SetDefault<Parameters::GridFile>(
"data/reservoir.dgf");
397 Parameters::SetDefault<Parameters::EndTime<Scalar>>(1000.0*24*60*60);
399 Parameters::SetDefault<Parameters::EnableStorageCache>(
true);
400 Parameters::SetDefault<Parameters::GridFile>(
"data/reservoir.dgf");
401 Parameters::SetDefault<Parameters::InitialTimeStepSize<Scalar>>(100e3);
403 Parameters::SetDefault<Parameters::NewtonTolerance<Scalar>>(1e-6);
405 Parameters::SetDefault<Parameters::EnableGravity>(
true);
412 {
return std::string(
"reservoir_") + Model::name() +
"_" + Model::discretizationName(); }
422 this->simulator().startNextEpisode(1e100);
423 this->simulator().setTimeStepSize(5.0);
438 this->model().globalStorage(storage);
441 if (this->gridView().comm().rank() == 0) {
442 std::cout <<
"Storage: " << storage << std::endl << std::flush;
453 template <
class Context>
455 unsigned timeIdx)
const
457 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
458 if (isFineMaterial_(pos))
466 template <
class Context>
467 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
469 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
470 if (isFineMaterial_(pos))
471 return finePorosity_;
472 return coarsePorosity_;
478 template <
class Context>
480 unsigned spaceIdx,
unsigned timeIdx)
const
482 unsigned globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
483 return *materialParams_[globalIdx];
487 {
return *materialParams_[globalIdx]; }
503 template <
class Context>
507 {
return temperature_; }
522 template <
class Context>
545 template <
class Context>
551 values.assignNaive(initialFluidState_);
554 for (
unsigned pvIdx = 0; pvIdx < values.size(); ++ pvIdx)
555 assert(std::isfinite(values[pvIdx]));
567 template <
class Context>
569 const Context& context,
571 unsigned timeIdx)
const
573 if (this->simulator().episodeIndex() == 1)
576 const auto& pos = context.pos(spaceIdx, timeIdx);
577 if (isInjector_(pos)) {
578 constraintValues.setActive(
true);
579 constraintValues.assignNaive(injectorFluidState_);
581 else if (isProducer_(pos)) {
582 constraintValues.setActive(
true);
583 constraintValues.assignNaive(producerFluidState_);
592 template <
class Context>
597 { rate = Scalar(0.0); }
602 void initFluidState_()
604 auto& fs = initialFluidState_;
609 fs.setTemperature(temperature_);
614 fs.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
615 fs.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
616 fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
621 Scalar pw = pReservoir_;
624 const auto& matParams = fineMaterialParams_;
625 MaterialLaw::capillaryPressures(pC, matParams, fs);
627 fs.setPressure(oilPhaseIdx, pw + (pC[oilPhaseIdx] - pC[waterPhaseIdx]));
628 fs.setPressure(waterPhaseIdx, pw + (pC[waterPhaseIdx] - pC[waterPhaseIdx]));
629 fs.setPressure(gasPhaseIdx, pw + (pC[gasPhaseIdx] - pC[waterPhaseIdx]));
632 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
633 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
634 fs.setMoleFraction(phaseIdx, compIdx, 0.0);
639 fs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
640 fs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
646 FluidSystem::saturatedDissolutionFactor(fs, oilPhaseIdx, 0);
647 Scalar XoGSat = FluidSystem::convertRsToXoG(RsSat, 0);
648 Scalar xoGSat = FluidSystem::convertXoGToxoG(XoGSat, 0);
649 Scalar xoG = 0.95*xoGSat;
650 Scalar xoO = 1.0 - xoG;
653 fs.setMoleFraction(oilPhaseIdx, gasCompIdx, xoG);
654 fs.setMoleFraction(oilPhaseIdx, oilCompIdx, xoO);
656 using CFRP = Opm::ComputeFromReferencePhase<Scalar, FluidSystem>;
657 typename FluidSystem::template ParameterCache<Scalar> paramCache;
665 auto& injFs = injectorFluidState_;
666 injFs = initialFluidState_;
668 Scalar pInj = pReservoir_ * 1.5;
669 injFs.setPressure(waterPhaseIdx, pInj);
670 injFs.setPressure(oilPhaseIdx, pInj);
671 injFs.setPressure(gasPhaseIdx, pInj);
672 injFs.setSaturation(waterPhaseIdx, 1.0);
673 injFs.setSaturation(oilPhaseIdx, 0.0);
674 injFs.setSaturation(gasPhaseIdx, 0.0);
677 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
678 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
679 injFs.setMoleFraction(phaseIdx, compIdx, 0.0);
681 injFs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
682 injFs.setMoleFraction(oilPhaseIdx, oilCompIdx, 1.0);
683 injFs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
692 auto& prodFs = producerFluidState_;
693 prodFs = initialFluidState_;
695 Scalar pProd = pReservoir_ / 1.5;
696 prodFs.setPressure(waterPhaseIdx, pProd);
697 prodFs.setPressure(oilPhaseIdx, pProd);
698 prodFs.setPressure(gasPhaseIdx, pProd);
699 prodFs.setSaturation(waterPhaseIdx, 0.0);
700 prodFs.setSaturation(oilPhaseIdx, 1.0);
701 prodFs.setSaturation(gasPhaseIdx, 0.0);
710 bool isProducer_(
const GlobalPosition& pos)
const
712 Scalar x = pos[0] - this->boundingBoxMin()[0];
713 Scalar y = pos[dim - 1] - this->boundingBoxMin()[dim - 1];
714 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
715 Scalar height = this->boundingBoxMax()[dim - 1] - this->boundingBoxMin()[dim - 1];
722 return width/2.0 - width*1e-5 < x && x < width/2.0 + width*(wellWidth_ + 1e-5);
725 bool isInjector_(
const GlobalPosition& pos)
const
727 Scalar x = pos[0] - this->boundingBoxMin()[0];
728 Scalar y = pos[dim - 1] - this->boundingBoxMin()[dim - 1];
729 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
730 Scalar height = this->boundingBoxMax()[dim - 1] - this->boundingBoxMin()[dim - 1];
737 return x < width*wellWidth_ - width*1e-5 || x > width*(1.0 - wellWidth_) + width*1e-5;
740 bool isFineMaterial_(
const GlobalPosition& pos)
const
741 {
return pos[dim - 1] > layerBottom_; }
748 Scalar finePorosity_;
749 Scalar coarsePorosity_;
751 MaterialLawParams fineMaterialParams_;
752 MaterialLawParams coarseMaterialParams_;
753 std::vector<const MaterialLawParams*> materialParams_;
755 InitialFluidState initialFluidState_;
756 InitialFluidState injectorFluidState_;
757 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:164
static void registerParameters()
Definition: reservoirproblem.hh:381
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:593
void endTimeStep()
Called by the simulator after each time integration.
Definition: reservoirproblem.hh:429
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:467
void initial(PrimaryVariables &values, const Context &, unsigned, unsigned) const
Evaluate the initial value for a control volume.
Definition: reservoirproblem.hh:546
ReservoirProblem(Simulator &simulator)
Definition: reservoirproblem.hh:210
void constraints(Constraints &constraintValues, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the constraints for a control volume.
Definition: reservoirproblem.hh:568
void boundary(BoundaryRateVector &values, const Context &, unsigned, unsigned) const
Evaluate the boundary conditions for a boundary segment.
Definition: reservoirproblem.hh:523
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:479
std::string name() const
The problem name.
Definition: reservoirproblem.hh:411
const MaterialLawParams & materialLawParams(unsigned globalIdx) const
Definition: reservoirproblem.hh:486
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: reservoirproblem.hh:217
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:454
Scalar temperature(const Context &, unsigned, unsigned) const
Definition: reservoirproblem.hh:504
void endEpisode()
Called when the end of an simulation episode is reached.
Definition: reservoirproblem.hh:417
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:79
Definition: blackoilbioeffectsmodules.hh:43
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:233
static constexpr Scalar value
Definition: co2injectionproblem.hh:160
static constexpr Scalar value
Definition: co2injectionproblem.hh:165
Definition: reservoirproblem.hh:140
static constexpr Scalar value
Definition: reservoirproblem.hh:140
Specify whether the some degrees of fredom can be constraint.
Definition: fvbaseproperties.hh:199
Opm::BlackOilFluidSystem< Scalar > type
Definition: reservoirproblem.hh:123
The fluid systems including the information about the phases.
Definition: multiphasebaseproperties.hh:79
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:101
The material law which ought to be used (extracted from the spatial parameters)
Definition: multiphasebaseproperties.hh:55
The type of the problem.
Definition: fvbaseproperties.hh:81
Definition: reservoirproblem.hh:72