28 #ifndef EWOMS_RESERVOIR_PROBLEM_HH
29 #define EWOMS_RESERVOIR_PROBLEM_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/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>
57 namespace Properties {
75 SET_PROP(ReservoirBaseProblem, MaterialLaw)
79 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
82 ThreePhaseMaterialTraits<Scalar,
83 FluidSystem::waterPhaseIdx,
84 FluidSystem::oilPhaseIdx,
85 FluidSystem::gasPhaseIdx> Traits;
88 typedef Opm::LinearMaterial<Traits> type;
92 SET_BOOL_PROP(ReservoirBaseProblem, NewtonWriteConvergence, false);
124 SET_PROP(ReservoirBaseProblem, FluidSystem)
130 typedef Opm::FluidSystems::BlackOil<Scalar> type;
156 template <class TypeTag>
159 typedef typename GET_PROP_TYPE(TypeTag, BaseProblem) ParentType;
163 typedef typename
GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
164 typedef typename
GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
167 enum { dim = GridView::dimension };
168 enum { dimWorld = GridView::dimensionworld };
171 enum { numPhases = FluidSystem::numPhases };
172 enum { numComponents = FluidSystem::numComponents };
173 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
174 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
175 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
176 enum { gasCompIdx = FluidSystem::gasCompIdx };
177 enum { oilCompIdx = FluidSystem::oilCompIdx };
178 enum { waterCompIdx = FluidSystem::waterCompIdx };
181 typedef typename
GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
182 typedef typename
GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
184 typedef typename
GET_PROP_TYPE(TypeTag, RateVector) RateVector;
185 typedef typename
GET_PROP_TYPE(TypeTag, BoundaryRateVector) BoundaryRateVector;
186 typedef typename
GET_PROP_TYPE(TypeTag, Constraints) Constraints;
187 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
188 typedef typename
GET_PROP_TYPE(TypeTag, Simulator) Simulator;
189 typedef typename
GET_PROP_TYPE(TypeTag, MaterialLawParams) MaterialLawParams;
191 typedef typename GridView::ctype CoordScalar;
192 typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
193 typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> DimMatrix;
194 typedef Dune::FieldVector<Scalar, numPhases> PhaseVector;
196 typedef Opm::CompositionalFluidState<Scalar,
198 true> InitialFluidState;
205 : ParentType(simulator)
213 ParentType::finishInit();
219 std::vector<std::pair<Scalar, Scalar> > Bo = {
221 { 1.82504e+06, 1.15 },
222 { 3.54873e+06, 1.207 },
223 { 6.99611e+06, 1.295 },
224 { 1.38909e+07, 1.435 },
225 { 1.73382e+07, 1.5 },
226 { 2.07856e+07, 1.565 },
227 { 2.76804e+07, 1.695 },
228 { 3.45751e+07, 1.827 }
230 std::vector<std::pair<Scalar, Scalar> > muo = {
232 { 1.82504e+06, 0.000975 },
233 { 3.54873e+06, 0.00091 },
234 { 6.99611e+06, 0.00083 },
235 { 1.38909e+07, 0.000695 },
236 { 1.73382e+07, 0.000641 },
237 { 2.07856e+07, 0.000594 },
238 { 2.76804e+07, 0.00051 },
239 { 3.45751e+07, 0.000449 }
241 std::vector<std::pair<Scalar, Scalar> > Rs = {
242 { 101353, 0.178108 },
243 { 1.82504e+06, 16.1187 },
244 { 3.54873e+06, 32.0594 },
245 { 6.99611e+06, 66.0779 },
246 { 1.38909e+07, 113.276 },
247 { 1.73382e+07, 138.033 },
248 { 2.07856e+07, 165.64 },
249 { 2.76804e+07, 226.197 },
250 { 3.45751e+07, 288.178 }
252 std::vector<std::pair<Scalar, Scalar> > Bg = {
254 { 1.82504e+06, 0.0678972 },
255 { 3.54873e+06, 0.0352259 },
256 { 6.99611e+06, 0.0179498 },
257 { 1.38909e+07, 0.00906194 },
258 { 1.73382e+07, 0.00726527 },
259 { 2.07856e+07, 0.00606375 },
260 { 2.76804e+07, 0.00455343 },
261 { 3.45751e+07, 0.00364386 },
262 { 6.21542e+07, 0.00216723 }
264 std::vector<std::pair<Scalar, Scalar> > mug = {
266 { 1.82504e+06, 9.6e-06 },
267 { 3.54873e+06, 1.12e-05 },
268 { 6.99611e+06, 1.4e-05 },
269 { 1.38909e+07, 1.89e-05 },
270 { 1.73382e+07, 2.08e-05 },
271 { 2.07856e+07, 2.28e-05 },
272 { 2.76804e+07, 2.68e-05 },
273 { 3.45751e+07, 3.09e-05 },
274 { 6.21542e+07, 4.7e-05 }
277 Scalar rhoRefO = 786.0;
278 Scalar rhoRefG = 0.97;
279 Scalar rhoRefW = 1037.0;
280 FluidSystem::initBegin(1);
281 FluidSystem::setEnableDissolvedGas(
true);
282 FluidSystem::setEnableVaporizedOil(
false);
283 FluidSystem::setReferenceDensities(rhoRefO, rhoRefW, rhoRefG, 0);
285 Opm::GasPvtMultiplexer<Scalar> *gasPvt =
new Opm::GasPvtMultiplexer<Scalar>;
286 gasPvt->setApproach(Opm::GasPvtMultiplexer<Scalar>::DryGasPvt);
287 auto& dryGasPvt = gasPvt->template getRealPvt<Opm::GasPvtMultiplexer<Scalar>::DryGasPvt>();
288 dryGasPvt.setNumRegions(1);
289 dryGasPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
290 dryGasPvt.setGasFormationVolumeFactor(0, Bg);
291 dryGasPvt.setGasViscosity(0, mug);
293 Opm::OilPvtMultiplexer<Scalar> *oilPvt =
new Opm::OilPvtMultiplexer<Scalar>;
294 oilPvt->setApproach(Opm::OilPvtMultiplexer<Scalar>::LiveOilPvt);
295 auto& liveOilPvt = oilPvt->template getRealPvt<Opm::OilPvtMultiplexer<Scalar>::LiveOilPvt>();
296 liveOilPvt.setNumRegions(1);
297 liveOilPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
298 liveOilPvt.setSaturatedOilGasDissolutionFactor(0, Rs);
299 liveOilPvt.setSaturatedOilFormationVolumeFactor(0, Bo);
300 liveOilPvt.setSaturatedOilViscosity(0, muo);
302 Opm::WaterPvtMultiplexer<Scalar> *waterPvt =
new Opm::WaterPvtMultiplexer<Scalar>;
303 waterPvt->setApproach(Opm::WaterPvtMultiplexer<Scalar>::ConstantCompressibilityWaterPvt);
304 auto& ccWaterPvt = waterPvt->template getRealPvt<Opm::WaterPvtMultiplexer<Scalar>::ConstantCompressibilityWaterPvt>();
305 ccWaterPvt.setNumRegions(1);
306 ccWaterPvt.setReferenceDensities(0, rhoRefO, rhoRefG, rhoRefW);
307 ccWaterPvt.setViscosity(0, 9.6e-4);
308 ccWaterPvt.setCompressibility(0, 1.450377e-10);
314 typedef std::shared_ptr<Opm::GasPvtMultiplexer<Scalar> > GasPvtSharedPtr;
315 FluidSystem::setGasPvt(GasPvtSharedPtr(gasPvt));
317 typedef std::shared_ptr<Opm::OilPvtMultiplexer<Scalar> > OilPvtSharedPtr;
318 FluidSystem::setOilPvt(OilPvtSharedPtr(oilPvt));
320 typedef std::shared_ptr<Opm::WaterPvtMultiplexer<Scalar> > WaterPvtSharedPtr;
321 FluidSystem::setWaterPvt(WaterPvtSharedPtr(waterPvt));
323 FluidSystem::initEnd();
329 fineK_ = this->toDimMatrix_(1e-12);
330 coarseK_ = this->toDimMatrix_(1e-11);
334 coarsePorosity_ = 0.3;
336 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
337 fineMaterialParams_.setPcMinSat(phaseIdx, 0.0);
338 fineMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
340 coarseMaterialParams_.setPcMinSat(phaseIdx, 0.0);
341 coarseMaterialParams_.setPcMaxSat(phaseIdx, 0.0);
345 fineMaterialParams_.finalize();
346 coarseMaterialParams_.finalize();
348 materialParams_.resize(this->model().numGridDof());
349 ElementContext elemCtx(this->simulator());
350 auto eIt = this->simulator().gridView().template begin<0>();
351 const auto& eEndIt = this->simulator().gridView().template end<0>();
352 for (; eIt != eEndIt; ++eIt) {
353 elemCtx.updateStencil(*eIt);
354 size_t nDof = elemCtx.numPrimaryDof(0);
355 for (
unsigned dofIdx = 0; dofIdx < nDof; ++ dofIdx) {
356 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, 0);
357 const GlobalPosition& pos = elemCtx.pos(dofIdx, 0);
359 if (isFineMaterial_(pos))
360 materialParams_[globalDofIdx] = &fineMaterialParams_;
362 materialParams_[globalDofIdx] = &coarseMaterialParams_;
369 this->simulator().startNextEpisode(100.0*24*60*60);
377 ParentType::registerParameters();
380 "The temperature [K] in the reservoir");
382 "The maximum depth [m] of the reservoir");
384 "The width of producer/injector wells as a fraction of the width"
385 " of the spatial domain");
392 {
return std::string(
"reservoir_") + Model::name() +
"_" + Model::discretizationName(); }
402 this->simulator().startNextEpisode(1e100);
403 this->simulator().setTimeStepSize(5.0);
418 this->model().globalStorage(storage);
421 if (this->gridView().comm().rank() == 0) {
422 std::cout <<
"Storage: " << storage << std::endl << std::flush;
433 template <
class Context>
435 unsigned timeIdx)
const
437 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
438 if (isFineMaterial_(pos))
446 template <
class Context>
447 Scalar
porosity(
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
449 const GlobalPosition& pos = context.pos(spaceIdx, timeIdx);
450 if (isFineMaterial_(pos))
451 return finePorosity_;
452 return coarsePorosity_;
458 template <
class Context>
460 unsigned spaceIdx,
unsigned timeIdx)
const
462 unsigned globalIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
463 return *materialParams_[globalIdx];
466 const MaterialLawParams& materialLawParams(
unsigned globalIdx)
const
467 {
return *materialParams_[globalIdx]; }
483 template <
class Context>
485 unsigned spaceIdx OPM_UNUSED,
486 unsigned timeIdx OPM_UNUSED)
const
487 {
return temperature_; }
502 template <
class Context>
504 const Context& context OPM_UNUSED,
505 unsigned spaceIdx OPM_UNUSED,
506 unsigned timeIdx OPM_UNUSED)
const
525 template <
class Context>
527 const Context& context OPM_UNUSED,
528 unsigned spaceIdx OPM_UNUSED,
529 unsigned timeIdx OPM_UNUSED)
const
531 values.assignNaive(initialFluidState_);
534 for (
unsigned pvIdx = 0; pvIdx < values.size(); ++ pvIdx)
535 assert(std::isfinite(values[pvIdx]));
547 template <
class Context>
548 void constraints(Constraints& constraints,
const Context& context,
549 unsigned spaceIdx,
unsigned timeIdx)
const
551 if (this->simulator().episodeIndex() == 1)
554 const auto& pos = context.pos(spaceIdx, timeIdx);
555 if (isInjector_(pos)) {
556 constraints.setActive(
true);
557 constraints.assignNaive(injectorFluidState_);
559 else if (isProducer_(pos)) {
560 constraints.setActive(
true);
561 constraints.assignNaive(producerFluidState_);
570 template <
class Context>
572 const Context& context OPM_UNUSED,
573 unsigned spaceIdx OPM_UNUSED,
574 unsigned timeIdx OPM_UNUSED)
const
575 { rate = Scalar(0.0); }
580 void initFluidState_()
582 auto& fs = initialFluidState_;
587 fs.setTemperature(temperature_);
592 fs.setSaturation(FluidSystem::oilPhaseIdx, 1.0);
593 fs.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
594 fs.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
599 Scalar pw = pReservoir_;
602 const auto& matParams = fineMaterialParams_;
603 MaterialLaw::capillaryPressures(pC, matParams, fs);
605 fs.setPressure(oilPhaseIdx, pw + (pC[oilPhaseIdx] - pC[waterPhaseIdx]));
606 fs.setPressure(waterPhaseIdx, pw + (pC[waterPhaseIdx] - pC[waterPhaseIdx]));
607 fs.setPressure(gasPhaseIdx, pw + (pC[gasPhaseIdx] - pC[waterPhaseIdx]));
610 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
611 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
612 fs.setMoleFraction(phaseIdx, compIdx, 0.0);
617 fs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
618 fs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
624 FluidSystem::saturatedDissolutionFactor(fs, oilPhaseIdx, 0);
625 Scalar XoGSat = FluidSystem::convertRsToXoG(RsSat, 0);
626 Scalar xoGSat = FluidSystem::convertXoGToxoG(XoGSat, 0);
627 Scalar xoG = 0.95*xoGSat;
628 Scalar xoO = 1.0 - xoG;
631 fs.setMoleFraction(oilPhaseIdx, gasCompIdx, xoG);
632 fs.setMoleFraction(oilPhaseIdx, oilCompIdx, xoO);
634 typedef Opm::ComputeFromReferencePhase<Scalar, FluidSystem> CFRP;
635 typename FluidSystem::template ParameterCache<Scalar> paramCache;
643 auto& injFs = injectorFluidState_;
644 injFs = initialFluidState_;
646 Scalar pInj = pReservoir_ * 1.5;
647 injFs.setPressure(waterPhaseIdx, pInj);
648 injFs.setPressure(oilPhaseIdx, pInj);
649 injFs.setPressure(gasPhaseIdx, pInj);
650 injFs.setSaturation(waterPhaseIdx, 1.0);
651 injFs.setSaturation(oilPhaseIdx, 0.0);
652 injFs.setSaturation(gasPhaseIdx, 0.0);
655 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
656 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
657 injFs.setMoleFraction(phaseIdx, compIdx, 0.0);
659 injFs.setMoleFraction(gasPhaseIdx, gasCompIdx, 1.0);
660 injFs.setMoleFraction(oilPhaseIdx, oilCompIdx, 1.0);
661 injFs.setMoleFraction(waterPhaseIdx, waterCompIdx, 1.0);
670 auto& prodFs = producerFluidState_;
671 prodFs = initialFluidState_;
673 Scalar pProd = pReservoir_ / 1.5;
674 prodFs.setPressure(waterPhaseIdx, pProd);
675 prodFs.setPressure(oilPhaseIdx, pProd);
676 prodFs.setPressure(gasPhaseIdx, pProd);
677 prodFs.setSaturation(waterPhaseIdx, 0.0);
678 prodFs.setSaturation(oilPhaseIdx, 1.0);
679 prodFs.setSaturation(gasPhaseIdx, 0.0);
688 bool isProducer_(
const GlobalPosition& pos)
const
690 Scalar x = pos[0] - this->boundingBoxMin()[0];
691 Scalar y = pos[dim - 1] - this->boundingBoxMin()[dim - 1];
692 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
693 Scalar height = this->boundingBoxMax()[dim - 1] - this->boundingBoxMin()[dim - 1];
700 return width/2.0 - width*1e-5 < x && x < width/2.0 + width*(wellWidth_ + 1e-5);
703 bool isInjector_(
const GlobalPosition& pos)
const
705 Scalar x = pos[0] - this->boundingBoxMin()[0];
706 Scalar y = pos[dim - 1] - this->boundingBoxMin()[dim - 1];
707 Scalar width = this->boundingBoxMax()[0] - this->boundingBoxMin()[0];
708 Scalar height = this->boundingBoxMax()[dim - 1] - this->boundingBoxMin()[dim - 1];
715 return x < width*wellWidth_ - width*1e-5 || x > width*(1.0 - wellWidth_) + width*1e-5;
718 bool isFineMaterial_(
const GlobalPosition& pos)
const
719 {
return pos[dim - 1] > layerBottom_; }
726 Scalar finePorosity_;
727 Scalar coarsePorosity_;
729 MaterialLawParams fineMaterialParams_;
730 MaterialLawParams coarseMaterialParams_;
731 std::vector<const MaterialLawParams*> materialParams_;
733 InitialFluidState initialFluidState_;
734 InitialFluidState injectorFluidState_;
735 InitialFluidState producerFluidState_;
#define SET_BOOL_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant boolean value.
Definition: propertysystem.hh:361
void endTimeStep()
Called by the simulator after each time integration.
Definition: reservoirproblem.hh:409
void source(RateVector &rate, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the source term for all phases within a given sub-control-volume.
Definition: reservoirproblem.hh:571
Some simple test problem for the black-oil VCVF discretization inspired by an oil reservoir...
Definition: reservoirproblem.hh:55
void boundary(BoundaryRateVector &values, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the boundary conditions for a boundary segment.
Definition: reservoirproblem.hh:503
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:447
#define NEW_TYPE_TAG(...)
Define a new type tag.
Definition: propertysystem.hh:169
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
#define EWOMS_REGISTER_PARAM(TypeTag, ParamType, ParamName, Description)
Register a run-time parameter.
Definition: parametersystem.hh:68
Declares the properties required by the black oil model.
Scalar temperature(const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Definition: reservoirproblem.hh:484
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
#define SET_PROP(EffTypeTagName, PropTagName)
Set a property for a specific type tag.
Definition: propertysystem.hh:297
static void registerParameters()
Definition: reservoirproblem.hh:375
const MaterialLawParams & materialLawParams(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:459
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
void endEpisode()
Called when the end of an simulation episode is reached.
Definition: reservoirproblem.hh:397
#define SET_TYPE_PROP(EffTypeTagName, PropTagName,...)
Set a property which defines a type.
Definition: propertysystem.hh:377
#define SET_STRING_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant string value.
Definition: propertysystem.hh:416
void constraints(Constraints &constraints, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the constraints for a control volume.
Definition: reservoirproblem.hh:548
void finishInit()
Called by the Ewoms::Simulator in order to initialize the problem.
Definition: reservoirproblem.hh:211
std::string name() const
The problem name.
Definition: reservoirproblem.hh:391
#define SET_SCALAR_PROP(EffTypeTagName, PropTagName,...)
Set a property to a simple constant scalar value.
Definition: propertysystem.hh:394
void initial(PrimaryVariables &values, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED) const
Evaluate the initial value for a control volume.
Definition: reservoirproblem.hh:526
const DimMatrix & intrinsicPermeability(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: reservoirproblem.hh:434