30#ifndef OPM_FLOW_PROBLEM_COMP_HPP
31#define OPM_FLOW_PROBLEM_COMP_HPP
38#include <opm/material/fluidstates/CompositionalFluidState.hpp>
40#include <opm/material/thermal/EclThermalLawManager.hpp>
42#include <opm/input/eclipse/EclipseState/Compositional/CompositionalConfig.hpp>
58template <
class TypeTag>
88 using InitialFluidState = CompositionalFluidState<Scalar, FluidSystem>;
105 Parameters::SetDefault<Parameters::NewtonTolerance<Scalar>>(1e-7);
110 auto& simulator = this->simulator();
111 const auto& eclState = simulator.vanguard().eclState();
112 return eclState.compositionalConfig().eosType(0);
120 , thresholdPressures_(simulator)
122 eclWriter_ = std::make_unique<EclWriterType>(simulator);
123 enableEclOutput_ = Parameters::Get<Parameters::EnableEclOutput>();
134 FlowProblemType::finishInit();
136 auto& simulator = this->simulator();
138 auto finishTransmissibilities = [updated =
false,
this]()
mutable {
143 [&vg = this->simulator().vanguard()](
const unsigned int it) {
return vg.gridIdxToEquilGridIdx(it); });
148 finishTransmissibilities();
150 if (enableEclOutput_) {
151 eclWriter_->setTransmissibilities(&simulator.problem().eclTransmissibilities());
152 std::function<
unsigned int(
unsigned int)> equilGridToGrid = [&simulator](
unsigned int i) {
153 return simulator.vanguard().gridEquilIdxToGridIdx(i);
155 eclWriter_->extractOutputTransAndNNC(equilGridToGrid);
158 const auto& eclState = simulator.vanguard().eclState();
159 const auto& schedule = simulator.vanguard().schedule();
162 simulator.setStartTime(schedule.getStartTime());
163 simulator.setEndTime(schedule.simTime(schedule.size() - 1));
169 simulator.setEpisodeIndex(-1);
170 simulator.setEpisodeLength(0.0);
175 this->gravity_ = 0.0;
176 if (Parameters::Get<Parameters::EnableGravity>())
177 this->gravity_[dim - 1] = 9.80665;
178 if (!eclState.getInitConfig().hasGravity())
179 this->gravity_[dim - 1] = 0.0;
184 const auto& tuning = schedule[0].tuning();
191 if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)
192 && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
196 this->
readRockParameters_(simulator.vanguard().cellCenterDepths(), [&simulator](
const unsigned idx) {
197 std::array<int, dim> coords;
198 simulator.vanguard().cartesianCoordinate(idx, coords);
199 std::transform(coords.begin(), coords.end(), coords.begin(),
200 [](const auto c) { return c + 1; });
207 if (enableEclOutput_) {
208 eclWriter_->writeInit();
211 const auto& initconfig = eclState.getInitConfig();
212 if (initconfig.restartRequested())
219 if constexpr (getPropValue<TypeTag, Properties::EnablePolymer>()) {
220 const auto& vanguard = this->simulator().vanguard();
221 const auto& gridView = vanguard.gridView();
222 int numElements = gridView.size(0);
223 this->
polymer_.maxAdsorption.resize(numElements, 0.0);
237 if (this->enableVtkOutput_() && eclState.getIOConfig().initOnly()) {
238 simulator.setTimeStepSize(0.0);
245 if (!initconfig.restartRequested()) {
246 simulator.startNextEpisode(schedule.seconds(1));
247 simulator.setEpisodeIndex(0);
248 simulator.setTimeStepIndex(0);
257 FlowProblemType::endTimeStep();
259 const bool isSubStep = !this->simulator().episodeWillBeOver();
262 this->eclWriter_->mutableOutputModule().invalidateLocalData();
265 const auto& grid = this->simulator().vanguard().gridView().grid();
267 using GridType = std::remove_cv_t<std::remove_reference_t<
decltype(grid)>>;
268 constexpr bool isCpGrid = std::is_same_v<GridType, Dune::CpGrid>;
269 if (!isCpGrid || (grid.maxLevel() == 0)) {
270 this->eclWriter_->evalSummaryState(isSubStep);
276 if (enableEclOutput_){
277 eclWriter_->writeReports(timer);
287 FlowProblemType::writeOutput(verbose);
289 const bool isSubStep = !this->simulator().episodeWillBeOver();
291 data::Solution localCellData = {};
292 if (enableEclOutput_) {
293 if (Parameters::Get<Parameters::EnableWriteAllSolutions>() || !isSubStep) {
294 eclWriter_->writeOutput(std::move(localCellData), isSubStep);
304 template <
class Context>
306 const Context& context,
310 OPM_TIMEBLOCK_LOCAL(eclProblemBoundary);
311 if (!context.intersection(spaceIdx).boundary())
316 if (this->nonTrivialBoundaryConditions()) {
317 throw std::logic_error(
"boundary condition is not supported by compostional modeling yet");
327 template <
class Context>
330 const unsigned globalDofIdx = context.globalSpaceIndex(spaceIdx, timeIdx);
331 const auto& initial_fs = initialFluidStates_[globalDofIdx];
332 Opm::CompositionalFluidState<Scalar, FluidSystem> fs;
333 for (
unsigned p = 0; p < numPhases; ++p) {
335 fs.setPressure(p, initial_fs.pressure(p));
338 fs.setSaturation(p, initial_fs.saturation(p));
341 fs.setTemperature(initial_fs.temperature(p));
345 if (!zmf_initialization_) {
346 for (
unsigned p = 0; p < numPhases; ++p) {
347 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
348 fs.setMoleFraction(p, compIdx, initial_fs.moleFraction(p, compIdx));
353 const auto& eos_type = getEosType();
354 typename FluidSystem::template ParameterCache<Scalar> paramCache(eos_type);
355 paramCache.updatePhase(fs, FluidSystem::oilPhaseIdx);
356 paramCache.updatePhase(fs, FluidSystem::gasPhaseIdx);
357 fs.setDensity(FluidSystem::oilPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::oilPhaseIdx));
358 fs.setDensity(FluidSystem::gasPhaseIdx, FluidSystem::density(fs, paramCache, FluidSystem::gasPhaseIdx));
361 Dune::FieldVector<Scalar, numComponents> z(0.0);
363 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
364 if (Indices::waterEnabled && phaseIdx ==
static_cast<unsigned int>(waterPhaseIdx)){
367 const auto saturation = fs.saturation(phaseIdx);
368 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
369 Scalar tmp = fs.molarity(phaseIdx, compIdx) * saturation;
370 tmp = max(tmp, 1e-8);
376 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
377 fs.setMoleFraction(compIdx, z[compIdx]);
381 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
382 fs.setMoleFraction(compIdx, initial_fs.moleFraction(compIdx));
387 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
388 const auto& Ktmp = fs.wilsonK_(compIdx);
389 fs.setKvalue(compIdx, Ktmp);
392 const Scalar& Ltmp = -1.0;
395 values.assignNaive(fs);
404 {
return initialFluidStates_[globalDofIdx]; }
407 {
return initialFluidStates_; }
410 {
return initialFluidStates_; }
414 assert( !thresholdPressures_.enableThresholdPressure() &&
415 " Threshold Pressures are not supported by compostional simulation ");
416 return thresholdPressures_;
420 {
return *eclWriter_; }
423 {
return *eclWriter_; }
426 template<
class Serializer>
430 serializer(*eclWriter_);
441 throw std::logic_error(
"Equilibration is not supported by compositional modeling yet");
446 throw std::logic_error(
"Restarting is not supported by compositional modeling yet");
451 readExplicitInitialConditionCompositional_();
456 const auto& simulator = this->simulator();
457 const auto& vanguard = simulator.vanguard();
458 const auto& eclState = vanguard.eclState();
459 const auto& fp = eclState.fieldProps();
460 const bool has_pressure = fp.has_double(
"PRESSURE");
462 throw std::runtime_error(
"The ECL input file requires the presence of the PRESSURE "
463 "keyword if the model is initialized explicitly");
465 const bool has_xmf = fp.has_double(
"XMF");
466 const bool has_ymf = fp.has_double(
"YMF");
467 const bool has_zmf = fp.has_double(
"ZMF");
468 if ( !has_zmf && !(has_xmf && has_ymf) ) {
469 throw std::runtime_error(
"The ECL input file requires the presence of ZMF or XMF and YMF "
470 "keyword if the model is initialized explicitly");
473 if (has_zmf && (has_xmf || has_ymf)) {
474 throw std::runtime_error(
"The ECL input file can not handle explicit initialization "
475 "with both ZMF and XMF or YMF");
478 if (has_xmf != has_ymf) {
479 throw std::runtime_error(
"The ECL input file needs XMF and YMF combined to do the explicit "
480 "initializtion when using XMF or YMF");
483 const bool has_temp = fp.has_double(
"TEMPI");
486 assert(fp.has_double(
"SGAS"));
488 std::size_t numDof = this->model().numGridDof();
490 initialFluidStates_.resize(numDof);
492 std::vector<double> waterSaturationData;
493 std::vector<double> gasSaturationData;
494 std::vector<double> soilData;
495 std::vector<double> pressureData;
496 std::vector<double> tempiData;
498 const bool water_active = FluidSystem::phaseIsActive(waterPhaseIdx);
499 const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx);
500 const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx);
502 if (water_active && Indices::numPhases > 2)
503 waterSaturationData = fp.get_double(
"SWAT");
505 waterSaturationData.resize(numDof);
507 pressureData = fp.get_double(
"PRESSURE");
510 tempiData = fp.get_double(
"TEMPI");
516 gasSaturationData = fp.get_double(
"SGAS");
518 gasSaturationData.resize(numDof);
520 for (std::size_t dofIdx = 0; dofIdx < numDof; ++dofIdx) {
521 auto& dofFluidState = initialFluidStates_[dofIdx];
524 Scalar temperatureLoc = tempiData[dofIdx];
525 assert(std::isfinite(temperatureLoc) && temperatureLoc > 0);
526 dofFluidState.setTemperature(temperatureLoc);
529 dofFluidState.setSaturation(FluidSystem::gasPhaseIdx,
530 gasSaturationData[dofIdx]);
533 dofFluidState.setSaturation(FluidSystem::oilPhaseIdx,
535 - waterSaturationData[dofIdx]
536 - gasSaturationData[dofIdx]);
539 dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
540 waterSaturationData[dofIdx]);
546 const Scalar pressure = pressureData[dofIdx];
549 const std::array<Scalar, numPhases> pc = {0};
550 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
551 if (!FluidSystem::phaseIsActive(phaseIdx))
554 if (Indices::oilEnabled)
555 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[oilPhaseIdx]));
556 else if (Indices::gasEnabled)
557 dofFluidState.setPressure(phaseIdx, pressure + (pc[phaseIdx] - pc[gasPhaseIdx]));
558 else if (Indices::waterEnabled)
560 dofFluidState.setPressure(phaseIdx, pressure);
563 if (has_xmf && has_ymf) {
564 const auto& xmfData = fp.get_double(
"XMF");
565 const auto& ymfData = fp.get_double(
"YMF");
566 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
567 const std::size_t data_idx = compIdx * numDof + dofIdx;
568 const Scalar xmf = xmfData[data_idx];
569 const Scalar ymf = ymfData[data_idx];
571 dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf);
572 dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf);
577 zmf_initialization_ =
true;
578 const auto& zmfData = fp.get_double(
"ZMF");
579 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
580 const std::size_t data_idx = compIdx * numDof + dofIdx;
581 const Scalar zmf = zmfData[data_idx];
582 dofFluidState.setMoleFraction(compIdx, zmf);
585 const auto ymf = (dofFluidState.saturation(FluidSystem::gasPhaseIdx) > 0.) ? zmf :
Scalar{0};
586 dofFluidState.setMoleFraction(FluidSystem::gasPhaseIdx, compIdx, ymf);
589 const auto xmf = (dofFluidState.saturation(FluidSystem::oilPhaseIdx) > 0.) ? zmf :
Scalar{0};
590 dofFluidState.setMoleFraction(FluidSystem::oilPhaseIdx, compIdx, xmf);
599 void handleSolventBC(
const BCProp::BCFace& , RateVector& )
const override
601 throw std::logic_error(
"solvent is disabled for compositional modeling and you're trying to add solvent to BC");
604 void handlePolymerBC(
const BCProp::BCFace& , RateVector& )
const override
606 throw std::logic_error(
"polymer is disabled for compositional modeling and you're trying to add polymer to BC");
609 void handleMicrBC(
const BCProp::BCFace& , RateVector& )
const override
611 throw std::logic_error(
"MICP is disabled for compositional modeling and you're trying to add microbes to BC");
614 void handleOxygBC(
const BCProp::BCFace& , RateVector& )
const override
616 throw std::logic_error(
"MICP is disabled for compositional modeling and you're trying to add oxygen to BC");
619 void handleUreaBC(
const BCProp::BCFace& , RateVector& )
const override
621 throw std::logic_error(
"MICP is disabled for compositional modeling and you're trying to add urea to BC");
624 FlowThresholdPressure<TypeTag> thresholdPressures_;
626 std::vector<InitialFluidState> initialFluidStates_;
628 bool zmf_initialization_ {
false};
630 bool enableEclOutput_{
false};
631 std::unique_ptr<EclWriterType> eclWriter_;
Collects necessary output values and pass it to opm-common's ECL output.
Definition: EclWriter.hpp:114
static void registerParameters()
Definition: EclWriter.hpp:137
PolymerSolutionContainer< Scalar > polymer_
Definition: FlowGenericProblem.hpp:334
Scalar initialTimeStepSize_
Definition: FlowGenericProblem.hpp:345
bool enableTuning_
Definition: FlowGenericProblem.hpp:344
void readRockParameters_(const std::vector< Scalar > &cellCenterDepths, std::function< std::array< int, 3 >(const unsigned)> ijkIndex)
Definition: FlowGenericProblem_impl.hpp:133
std::vector< Scalar > maxOilSaturation_
Definition: FlowGenericProblem.hpp:335
void initFluidSystem_()
Definition: FlowGenericProblem_impl.hpp:472
Scalar maxTimeStepAfterWellEvent_
Definition: FlowGenericProblem.hpp:346
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition: FlowProblemComp.hpp:60
void writeOutput(bool verbose) override
Write the requested quantities of the current solution into the output files.
Definition: FlowProblemComp.hpp:285
const std::vector< InitialFluidState > & initialFluidStates() const
Definition: FlowProblemComp.hpp:409
void finishInit()
Called by the Opm::Simulator in order to initialize the problem.
Definition: FlowProblemComp.hpp:129
Opm::CompositionalConfig::EOSType getEosType() const
Definition: FlowProblemComp.hpp:108
FlowProblemComp(Simulator &simulator)
Definition: FlowProblemComp.hpp:118
void writeReports(const SimulatorTimer &timer)
Definition: FlowProblemComp.hpp:275
void readExplicitInitialCondition_() override
Definition: FlowProblemComp.hpp:449
void readExplicitInitialConditionCompositional_()
Definition: FlowProblemComp.hpp:454
void endTimeStep() override
Called by the simulator after each time integration.
Definition: FlowProblemComp.hpp:255
const EclWriterType & eclWriter() const
Definition: FlowProblemComp.hpp:419
std::vector< InitialFluidState > & initialFluidStates()
Definition: FlowProblemComp.hpp:406
const FlowThresholdPressure< TypeTag > & thresholdPressure() const
Definition: FlowProblemComp.hpp:412
void readEclRestartSolution_()
Definition: FlowProblemComp.hpp:444
const InitialFluidState & initialFluidState(unsigned globalDofIdx) const
Definition: FlowProblemComp.hpp:403
void boundary(BoundaryRateVector &values, const Context &context, unsigned spaceIdx, unsigned) const
Evaluate the boundary conditions for a boundary segment.
Definition: FlowProblemComp.hpp:305
void initial(PrimaryVariables &values, const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Evaluate the initial value for a control volume.
Definition: FlowProblemComp.hpp:328
void readEquilInitialCondition_() override
Definition: FlowProblemComp.hpp:439
void serializeOp(Serializer &serializer)
Definition: FlowProblemComp.hpp:427
void updateExplicitQuantities_(int, int, bool) override
Definition: FlowProblemComp.hpp:434
void addToSourceDense(RateVector &, unsigned, unsigned) const override
Definition: FlowProblemComp.hpp:398
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition: FlowProblemComp.hpp:98
EclWriterType & eclWriter()
Definition: FlowProblemComp.hpp:422
This problem simulates an input file given in the data format used by the commercial ECLiPSE simulato...
Definition: FlowProblem.hpp:94
GetPropType< TypeTag, Properties::Evaluation > Evaluation
Definition: FlowProblem.hpp:156
@ dim
Definition: FlowProblem.hpp:110
virtual void writeOutput(bool verbose)
Write the requested quantities of the current solution into the output files.
Definition: FlowProblem.hpp:473
unsigned pvtRegionIndex(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Returns the index of the relevant region for thermodynmic properties.
Definition: FlowProblem.hpp:818
Scalar porosity(const Context &context, unsigned spaceIdx, unsigned timeIdx) const
Definition: FlowProblem.hpp:652
GetPropType< TypeTag, Properties::Vanguard > Vanguard
Definition: FlowProblem.hpp:107
@ oilPhaseIdx
Definition: FlowProblem.hpp:135
@ numPhases
Definition: FlowProblem.hpp:115
GetPropType< TypeTag, Properties::Scalar > Scalar
Definition: FlowProblem.hpp:101
GetPropType< TypeTag, Properties::RateVector > RateVector
Definition: FlowProblem.hpp:145
GetPropType< TypeTag, Properties::Indices > Indices
Definition: FlowProblem.hpp:157
@ gasPhaseIdx
Definition: FlowProblem.hpp:134
GetPropType< TypeTag, Properties::Simulator > Simulator
Definition: FlowProblem.hpp:146
@ dimWorld
Definition: FlowProblem.hpp:111
@ waterPhaseIdx
Definition: FlowProblem.hpp:136
void readThermalParameters_()
Definition: FlowProblem.hpp:1367
GetPropType< TypeTag, Properties::GridView > GridView
Definition: FlowProblem.hpp:102
static void registerParameters()
Registers all available parameters for the problem and the model.
Definition: FlowProblem.hpp:180
void updatePffDofData_()
Definition: FlowProblem.hpp:1480
GetPropType< TypeTag, Properties::PrimaryVariables > PrimaryVariables
Definition: FlowProblem.hpp:144
Vanguard::TransmissibilityType transmissibilities_
Definition: FlowProblem.hpp:1667
virtual void readInitialCondition_()
Definition: FlowProblem.hpp:1408
@ numComponents
Definition: FlowProblem.hpp:116
GetPropType< TypeTag, Properties::FluidSystem > FluidSystem
Definition: FlowProblem.hpp:104
GetPropType< TypeTag, Properties::MaterialLaw > MaterialLaw
Definition: FlowProblem.hpp:154
void readMaterialParameters_()
Definition: FlowProblem.hpp:1333
This class calculates the threshold pressure for grid faces according to the Eclipse Reference Manual...
Definition: FlowThresholdPressure.hpp:59
Definition: SimulatorTimer.hpp:39
Definition: blackoilboundaryratevector.hh:39
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