262 +
" ROCK2D tables is expected, but " + std::to_string(rock2dTables.size()) +
" is provided");
264 if (rockwnodTables.size() != numRocktabTables)
265 throw std::runtime_error(
"Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables)
266 +
" ROCKWNOD tables is expected, but " + std::to_string(rockwnodTables.size()) +
" is provided");
268 rockCompPoroMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical));
269 for (std::size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
270 const RockwnodTable& rockwnodTable = rockwnodTables.template getTable<RockwnodTable>(regionIdx);
271 const auto& rock2dTable = rock2dTables[regionIdx];
273 if (rockwnodTable.getSaturationColumn().size() != rock2dTable.sizeMultValues())
274 throw std::runtime_error(
"Number of entries in ROCKWNOD and ROCK2D needs to match.");
276 for (std::size_t xIdx = 0; xIdx < rock2dTable.size(); ++xIdx) {
277 rockCompPoroMultWc_[regionIdx].appendXPos(rock2dTable.getPressureValue(xIdx));
278 for (std::size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx)
279 rockCompPoroMultWc_[regionIdx].appendSamplePoint(xIdx,
280 rockwnodTable.getSaturationColumn()[yIdx],
281 rock2dTable.getPvmultValue(xIdx, yIdx));
285 if (!rock2dtrTables.empty()) {
286 rockCompTransMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical));
287 for (std::size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
288 const RockwnodTable& rockwnodTable = rockwnodTables.template getTable<RockwnodTable>(regionIdx);
289 const auto& rock2dtrTable = rock2dtrTables[regionIdx];
291 if (rockwnodTable.getSaturationColumn().size() != rock2dtrTable.sizeMultValues())
292 throw std::runtime_error(
"Number of entries in ROCKWNOD and ROCK2DTR needs to match.");
294 for (std::size_t xIdx = 0; xIdx < rock2dtrTable.size(); ++xIdx) {
295 rockCompTransMultWc_[regionIdx].appendXPos(rock2dtrTable.getPressureValue(xIdx));
296 for (std::size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx)
297 rockCompTransMultWc_[regionIdx].appendSamplePoint(xIdx,
298 rockwnodTable.getSaturationColumn()[yIdx],
299 rock2dtrTable.getTransMultValue(xIdx, yIdx));
306template<
class Gr
idView,
class Flu
idSystem>
307typename FlowGenericProblem<GridView,FluidSystem>::Scalar
311 if (this->rockParams_.empty())
314 unsigned tableIdx = 0;
315 if (!this->rockTableIdx_.empty()) {
316 tableIdx = this->rockTableIdx_[globalSpaceIdx];
318 return this->rockParams_[tableIdx].compressibility;
321template<
class Gr
idView,
class Flu
idSystem>
322typename FlowGenericProblem<GridView,FluidSystem>::Scalar
326 if (this->rockParams_.empty())
329 unsigned tableIdx = 0;
330 if (!this->rockTableIdx_.empty()) {
331 tableIdx = this->rockTableIdx_[globalSpaceIdx];
333 return this->rockParams_[tableIdx].referencePressure;
336template<
class Gr
idView,
class Flu
idSystem>
337typename FlowGenericProblem<GridView,FluidSystem>::Scalar
339porosity(
unsigned globalSpaceIdx,
unsigned timeIdx)
const
341 return this->referencePorosity_[timeIdx][globalSpaceIdx];
344template<
class Gr
idView,
class Flu
idSystem>
345typename FlowGenericProblem<GridView,FluidSystem>::Scalar
347rockFraction(
unsigned elementIdx,
unsigned timeIdx)
const
353 auto porosity = this->lookUpData_.fieldPropDouble(eclState_.fieldProps(),
"PORO", elementIdx);
354 return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity);
357template<
class Gr
idView,
class Flu
idSystem>
360updateNum(
const std::string& name, std::vector<T>& numbers, std::size_t num_regions)
362 if (!eclState_.fieldProps().has_int(name))
365 std::function<void(T,
int)> valueCheck = [num_regions,name](T fieldPropValue, [[maybe_unused]]
int fieldPropIdx) {
366 if ( fieldPropValue > (
int)num_regions) {
367 throw std::runtime_error(
"Values larger than maximum number of regions "
368 + std::to_string(num_regions) +
" provided in " + name);
370 if ( fieldPropValue <= 0) {
371 throw std::runtime_error(
"zero or negative values provided for region array: " + name);
375 numbers = this->lookUpData_.template assignFieldPropsIntOnLeaf<T>(eclState_.fieldProps(), name,
379template<
class Gr
idView,
class Flu
idSystem>
380void FlowGenericProblem<GridView,FluidSystem>::
383 const auto num_regions = eclState_.getTableManager().getTabdims().getNumPVTTables();
384 updateNum(
"PVTNUM", pvtnum_, num_regions);
387template<
class Gr
idView,
class Flu
idSystem>
388void FlowGenericProblem<GridView,FluidSystem>::
391 const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
392 updateNum(
"SATNUM", satnum_, num_regions);
395template<
class Gr
idView,
class Flu
idSystem>
396void FlowGenericProblem<GridView,FluidSystem>::
399 const auto num_regions = 1;
400 updateNum(
"MISCNUM", miscnum_, num_regions);
403template<
class Gr
idView,
class Flu
idSystem>
404void FlowGenericProblem<GridView,FluidSystem>::
407 const auto num_regions = 1;
408 updateNum(
"PLMIXNUM", plmixnum_, num_regions);
411template<
class Gr
idView,
class Flu
idSystem>
412void FlowGenericProblem<GridView,FluidSystem>::
415 const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
416 updateNum(
"KRNUMX", krnumx_, num_regions);
417 updateNum(
"KRNUMY", krnumy_, num_regions);
418 updateNum(
"KRNUMZ", krnumz_, num_regions);
421template<
class Gr
idView,
class Flu
idSystem>
422void FlowGenericProblem<GridView,FluidSystem>::
425 const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
426 updateNum(
"IMBNUMX", imbnumx_, num_regions);
427 updateNum(
"IMBNUMY", imbnumy_, num_regions);
428 updateNum(
"IMBNUMZ", imbnumz_, num_regions);
431template<
class Gr
idView,
class Flu
idSystem>
432bool FlowGenericProblem<GridView,FluidSystem>::
433vapparsActive(
int episodeIdx)
const
435 const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
436 return (oilVaporizationControl.getType() == OilVaporizationProperties::OilVaporization::VAPPARS);
439template<
class Gr
idView,
class Flu
idSystem>
440bool FlowGenericProblem<GridView,FluidSystem>::
441beginEpisode_(
bool enableExperiments,
444 if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
447 std::ostringstream ss;
448 boost::posix_time::time_facet* facet =
new boost::posix_time::time_facet(
"%d-%b-%Y");
449 boost::posix_time::ptime curDateTime =
450 boost::posix_time::from_time_t(schedule_.simTime(episodeIdx));
451 ss.imbue(std::locale(std::locale::classic(), facet));
452 ss <<
"Report step " << episodeIdx + 1
453 <<
"/" << schedule_.size() - 1
454 <<
" at day " << schedule_.seconds(episodeIdx)/(24*3600)
455 <<
"/" << schedule_.seconds(schedule_.size() - 1)/(24*3600)
456 <<
", date = " << curDateTime.date()
458 OpmLog::info(ss.str());
461 const auto& events = schedule_[episodeIdx].events();
464 if (episodeIdx > 0 && enableTuning_ && events.hasEvent(ScheduleEvents::TUNING_CHANGE))
466 const auto& sched_state = schedule_[episodeIdx];
467 const auto& tuning = sched_state.tuning();
468 initialTimeStepSize_ = sched_state.max_next_tstep(enableTuning_);
469 maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
476template<
class Gr
idView,
class Flu
idSystem>
477void FlowGenericProblem<GridView,FluidSystem>::
478beginTimeStep_(
bool enableExperiments,
486 if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
487 std::ostringstream ss;
488 boost::posix_time::time_facet* facet =
new boost::posix_time::time_facet(
"%d-%b-%Y");
489 boost::posix_time::ptime date = boost::posix_time::from_time_t(startTime) + boost::posix_time::milliseconds(
static_cast<long long>(time / prefix::milli));
490 ss.imbue(std::locale(std::locale::classic(), facet));
491 ss <<
"\nTime step " << timeStepIndex <<
", stepsize "
492 << unit::convert::to(timeStepSize, unit::day) <<
" days,"
493 <<
" at day " << (double)unit::convert::to(time, unit::day)
494 <<
"/" << (double)unit::convert::to(endTime, unit::day)
495 <<
", date = " << date;
496 OpmLog::info(ss.str());
500 this->mixControls_.updateExplicitQuantities(episodeIdx, timeStepSize);
503template<
class Gr
idView,
class Flu
idSystem>
504void FlowGenericProblem<GridView,FluidSystem>::
507 FluidSystem::initFromState(eclState_, schedule_);
510template<
class Gr
idView,
class Flu
idSystem>
511void FlowGenericProblem<GridView,FluidSystem>::
512readBlackoilExtentionsInitialConditions_(std::size_t numDof,
515 bool enablePolymerMolarWeight,
519 if (eclState_.fieldProps().has_double(
"SSOL"))
520 solventSaturation_ = eclState_.fieldProps().get_double(
"SSOL");
522 solventSaturation_.resize(numDof, 0.0);
527 solventRsw_.resize(numDof, 0.0);
531 if (eclState_.fieldProps().has_double(
"SPOLY")) {
532 polymer_.concentration = eclState_.fieldProps().get_double(
"SPOLY");
534 polymer_.concentration.resize(numDof, 0.0);
538 if (enablePolymerMolarWeight) {
539 if (eclState_.fieldProps().has_double(
"SPOLYMW")) {
540 polymer_.moleWeight = eclState_.fieldProps().get_double(
"SPOLYMW");
542 polymer_.moleWeight.resize(numDof, 0.0);
547 if (eclState_.fieldProps().has_double(
"SMICR")) {
548 micp_.microbialConcentration = eclState_.fieldProps().get_double(
"SMICR");
550 micp_.microbialConcentration.resize(numDof, 0.0);
552 if (eclState_.fieldProps().has_double(
"SOXYG")) {
553 micp_.oxygenConcentration = eclState_.fieldProps().get_double(
"SOXYG");
555 micp_.oxygenConcentration.resize(numDof, 0.0);
557 if (eclState_.fieldProps().has_double(
"SUREA")) {
558 micp_.ureaConcentration = eclState_.fieldProps().get_double(
"SUREA");
560 micp_.ureaConcentration.resize(numDof, 0.0);
562 if (eclState_.fieldProps().has_double(
"SBIOF")) {
563 micp_.biofilmConcentration = eclState_.fieldProps().get_double(
"SBIOF");
565 micp_.biofilmConcentration.resize(numDof, 0.0);
567 if (eclState_.fieldProps().has_double(
"SCALC")) {
568 micp_.calciteConcentration = eclState_.fieldProps().get_double(
"SCALC");
570 micp_.calciteConcentration.resize(numDof, 0.0);
575template<
class Gr
idView,
class Flu
idSystem>
576typename FlowGenericProblem<GridView,FluidSystem>::Scalar
580 if (maxWaterSaturation_.empty())
583 return maxWaterSaturation_[globalDofIdx];
586template<
class Gr
idView,
class Flu
idSystem>
587typename FlowGenericProblem<GridView,FluidSystem>::Scalar
591 if (minRefPressure_.empty())
594 return minRefPressure_[globalDofIdx];
597template<
class Gr
idView,
class Flu
idSystem>
598typename FlowGenericProblem<GridView,FluidSystem>::Scalar
602 if (overburdenPressure_.empty())
605 return overburdenPressure_[elementIdx];
608template<
class Gr
idView,
class Flu
idSystem>
609typename FlowGenericProblem<GridView,FluidSystem>::Scalar
613 if (solventSaturation_.empty())
616 return solventSaturation_[elemIdx];
619template<
class Gr
idView,
class Flu
idSystem>
620typename FlowGenericProblem<GridView,FluidSystem>::Scalar
624 if (solventRsw_.empty())
627 return solventRsw_[elemIdx];
630template<
class Gr
idView,
class Flu
idSystem>
631typename FlowGenericProblem<GridView,FluidSystem>::Scalar
633drsdtcon(
unsigned elemIdx,
int episodeIdx)
const
635 return this->mixControls_.
drsdtcon(elemIdx, episodeIdx,
636 this->pvtRegionIndex(elemIdx));
639template<
class Gr
idView,
class Flu
idSystem>
640typename FlowGenericProblem<GridView,FluidSystem>::Scalar
644 if (polymer_.concentration.empty()) {
648 return polymer_.concentration[elemIdx];
651template<
class Gr
idView,
class Flu
idSystem>
652typename FlowGenericProblem<GridView,FluidSystem>::Scalar
656 if (polymer_.moleWeight.empty()) {
660 return polymer_.moleWeight[elemIdx];
663template<
class Gr
idView,
class Flu
idSystem>
664typename FlowGenericProblem<GridView,FluidSystem>::Scalar
668 if (micp_.microbialConcentration.empty()) {
675template<
class Gr
idView,
class Flu
idSystem>
676typename FlowGenericProblem<GridView,FluidSystem>::Scalar
680 if (micp_.oxygenConcentration.empty()) {
687template<
class Gr
idView,
class Flu
idSystem>
688typename FlowGenericProblem<GridView,FluidSystem>::Scalar
692 if (micp_.ureaConcentration.empty()) {
699template<
class Gr
idView,
class Flu
idSystem>
700typename FlowGenericProblem<GridView,FluidSystem>::Scalar
704 if (micp_.biofilmConcentration.empty()) {
711template<
class Gr
idView,
class Flu
idSystem>
712typename FlowGenericProblem<GridView,FluidSystem>::Scalar
716 if (micp_.calciteConcentration.empty()) {
723template<
class Gr
idView,
class Flu
idSystem>
730 return pvtnum_[elemIdx];
733template<
class Gr
idView,
class Flu
idSystem>
740 return satnum_[elemIdx];
743template<
class Gr
idView,
class Flu
idSystem>
747 if (miscnum_.empty())
750 return miscnum_[elemIdx];
753template<
class Gr
idView,
class Flu
idSystem>
757 if (plmixnum_.empty())
760 return plmixnum_[elemIdx];
763template<
class Gr
idView,
class Flu
idSystem>
764typename FlowGenericProblem<GridView,FluidSystem>::Scalar
768 if (polymer_.maxAdsorption.empty()) {
772 return polymer_.maxAdsorption[elemIdx];
775template<
class Gr
idView,
class Flu
idSystem>
779 return this->maxWaterSaturation_ == rhs.maxWaterSaturation_ &&
780 this->minRefPressure_ == rhs.minRefPressure_ &&
781 this->overburdenPressure_ == rhs.overburdenPressure_ &&
782 this->solventSaturation_ == rhs.solventSaturation_ &&
783 this->solventRsw_ == rhs.solventRsw_ &&
784 this->polymer_ == rhs.polymer_ &&
785 this->micp_ == rhs.micp_ &&
786 this->mixControls_ == rhs.mixControls_;