311class CO2PTProblem :
public GetPropType<TypeTag, Properties::BaseProblem>
313 using ParentType = GetPropType<TypeTag, Properties::BaseProblem>;
315 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
316 using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
317 using GridView = GetPropType<TypeTag, Properties::GridView>;
318 using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
319 enum { dim = GridView::dimension };
320 enum { dimWorld = GridView::dimensionworld };
321 using Indices = GetPropType<TypeTag, Properties::Indices>;
322 using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
323 using RateVector = GetPropType<TypeTag, Properties::RateVector>;
324 using BoundaryRateVector = GetPropType<TypeTag, Properties::BoundaryRateVector>;
325 using MaterialLaw = GetPropType<TypeTag, Properties::MaterialLaw>;
326 using Simulator = GetPropType<TypeTag, Properties::Simulator>;
327 using Model = GetPropType<TypeTag, Properties::Model>;
328 using MaterialLawParams = GetPropType<TypeTag, Properties::MaterialLawParams>;
330 using Toolbox = Opm::MathToolbox<Evaluation>;
331 using CoordScalar =
typename GridView::ctype;
333 enum { numPhases = FluidSystem::numPhases };
334 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
335 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
336 enum { conti0EqIdx = Indices::conti0EqIdx };
337 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
338 enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
339 enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
340 enum { enableGravity = getPropValue<TypeTag, Properties::EnableGravity>() };
342 using GlobalPosition = Dune::FieldVector<CoordScalar, dimWorld>;
343 using DimMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
344 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
345 using ComponentVector = Dune::FieldVector<Evaluation, numComponents>;
346 using FlashSolver = GetPropType<TypeTag, Properties::FlashSolver>;
349 using FluidState = Opm::CompositionalFluidState<Evaluation, FluidSystem, enableEnergy>;
354 : ParentType(simulator)
356 const Scalar epi_len = Parameters::get<TypeTag, Properties::EpisodeLength>();
357 simulator.setEpisodeLength(epi_len);
359 using CompParm =
typename FluidSystem::ComponentParam;
360 using CO2 = Opm::SimpleCO2<Scalar>;
361 using C1 = Opm::C1<Scalar>;
362 using C10 = Opm::C10<Scalar>;
363 FluidSystem::addComponent(CompParm {CO2::name(), CO2::molarMass(), CO2::criticalTemperature(),
364 CO2::criticalPressure(), CO2::criticalVolume(), CO2::acentricFactor()});
365 FluidSystem::addComponent(CompParm {C1::name(), C1::molarMass(), C1::criticalTemperature(),
366 C1::criticalPressure(), C1::criticalVolume(), C1::acentricFactor()});
367 FluidSystem::addComponent(CompParm{C10::name(), C10::molarMass(), C10::criticalTemperature(),
368 C10::criticalPressure(), C10::criticalVolume(), C10::acentricFactor()});
372 void initPetrophysics()
374 temperature_ = Parameters::get<TypeTag, Properties::Temperature>();
375 K_ = this->toDimMatrix_(9.869232667160131e-14);
380 template <
class Context>
382 gravity([[maybe_unused]]
const Context& context, [[maybe_unused]]
unsigned spaceIdx, [[maybe_unused]]
unsigned timeIdx)
const
387 const DimVector& gravity()
const
397 ParentType::finishInit();
407 ParentType::registerParameters();
409 Parameters::registerParam<TypeTag, Properties::Temperature>
410 (
"The temperature [K] in the reservoir");
411 Parameters::registerParam<TypeTag, Properties::Initialpressure>
412 (
"The initial pressure [Pa s] in the reservoir");
413 Parameters::registerParam<TypeTag, Properties::SimulationName>
414 (
"The name of the simulation used for the output files");
415 Parameters::registerParam<TypeTag, Properties::EpisodeLength>
416 (
"Time interval [s] for episode length");
424 std::ostringstream oss;
425 oss << Parameters::get<TypeTag, Properties::SimulationName>();
434 Scalar epi_len = Parameters::get<TypeTag, Properties::EpisodeLength>();
435 this->simulator().startNextEpisode(epi_len);
440 bool shouldWriteOutput()
442 return this->simulator().episodeWillBeOver() || (this->simulator().timeStepIndex() == -1);
447 bool shouldWriteRestartFile()
457 Scalar tol = this->model().newtonMethod().tolerance() * 1e5;
458 this->model().checkConservativeness(tol);
461 PrimaryVariables storageO, storageW;
462 this->model().globalPhaseStorage(storageO, oilPhaseIdx);
465 PrimaryVariables storageL, storageG;
466 this->model().globalPhaseStorage(storageL, 0);
467 this->model().globalPhaseStorage(storageG, 1);
479 template <
class Context>
480 void initial(PrimaryVariables& values,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
482 Opm::CompositionalFluidState<Evaluation, FluidSystem> fs;
483 initialFluidState(fs, context, spaceIdx, timeIdx);
484 values.assignNaive(fs);
488 template <
class Context>
489 Scalar temperature([[maybe_unused]]
const Context& context, [[maybe_unused]]
unsigned spaceIdx, [[maybe_unused]]
unsigned timeIdx)
const
496 template <
class Context>
497 const DimMatrix& intrinsicPermeability([[maybe_unused]]
const Context& context,
498 [[maybe_unused]]
unsigned spaceIdx,
499 [[maybe_unused]]
unsigned timeIdx)
const
505 template <
class Context>
506 Scalar porosity([[maybe_unused]]
const Context& context, [[maybe_unused]]
unsigned spaceIdx, [[maybe_unused]]
unsigned timeIdx)
const
508 int spatialIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
510 int prod = Parameters::get<TypeTag, Properties::CellsX>() - 1;
511 if (spatialIdx == inj || spatialIdx == prod) {
521 template <
class Context>
523 [[maybe_unused]]
unsigned spaceIdx,
524 [[maybe_unused]]
unsigned timeIdx)
const
531 template <
class Context>
532 void boundary(BoundaryRateVector& values,
536 { values.setNoFlow(); }
539 template <
class Context>
540 void source(RateVector& source_rate,
541 [[maybe_unused]]
const Context& context,
542 [[maybe_unused]]
unsigned spaceIdx,
543 [[maybe_unused]]
unsigned timeIdx)
const
545 source_rate = Scalar(0.0);
553 template <
class Flu
idState,
class Context>
554 void initialFluidState(FluidState& fs,
const Context& context,
unsigned spaceIdx,
unsigned timeIdx)
const
561 int prod = Parameters::get<TypeTag, Properties::CellsX>() - 1;
562 int spatialIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
563 ComponentVector comp;
564 comp[0] = Evaluation::createVariable(0.5, 1);
565 comp[1] = Evaluation::createVariable(0.3, 2);
566 comp[2] = 1. - comp[0] - comp[1];
567 if (spatialIdx == inj) {
568 comp[0] = Evaluation::createVariable(0.99, 1);
569 comp[1] = Evaluation::createVariable(0.01 - 1e-3, 2);
570 comp[2] = 1. - comp[0] - comp[1];
574 sat[1] = 1.0 - sat[0];
576 Scalar p0 = Parameters::get<TypeTag, Properties::Initialpressure>();
580 if (spatialIdx == inj) {
583 if (spatialIdx == prod) {
586 Evaluation p_init = Evaluation::createVariable(p0, 0);
588 fs.setPressure(FluidSystem::oilPhaseIdx, p_init);
589 fs.setPressure(FluidSystem::gasPhaseIdx, p_init);
591 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
592 fs.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, comp[compIdx]);
593 fs.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, comp[compIdx]);
597 fs.setSaturation(FluidSystem::oilPhaseIdx, sat[0]);
598 fs.setSaturation(FluidSystem::gasPhaseIdx, sat[1]);
600 fs.setTemperature(temperature_);
604 typename FluidSystem::template ParameterCache<Evaluation> paramCache;
605 paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
606 paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
607 fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
608 fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
609 fs.setViscosity(FluidSystem::oilPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::oilPhaseIdx));
610 fs.setViscosity(FluidSystem::gasPhaseIdx, FluidSystem::viscosity(fs, paramCache, FluidSystem::gasPhaseIdx));
614 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
615 const Evaluation Ktmp = fs.wilsonK_(compIdx);
616 fs.setKvalue(compIdx, Ktmp);
619 const Evaluation& Ltmp = -1.0;
626 MaterialLawParams mat_;