255 const auto& rock2dtrTable = rock2dtrTables[regionIdx];
257 if (rockwnodTable.getSaturationColumn().size() != rock2dtrTable.sizeMultValues())
258 throw std::runtime_error(
"Number of entries in ROCKWNOD and ROCK2DTR needs to match.");
260 for (std::size_t xIdx = 0; xIdx < rock2dtrTable.size(); ++xIdx) {
261 rockCompTransMultWc_[regionIdx].appendXPos(rock2dtrTable.getPressureValue(xIdx));
262 for (std::size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx)
263 rockCompTransMultWc_[regionIdx].appendSamplePoint(xIdx,
264 rockwnodTable.getSaturationColumn()[yIdx],
265 rock2dtrTable.getTransMultValue(xIdx, yIdx));
272template<
class Gr
idView,
class Flu
idSystem>
273typename FlowGenericProblem<GridView,FluidSystem>::Scalar
277 if (this->rockParams_.empty())
280 unsigned tableIdx = 0;
281 if (!this->rockTableIdx_.empty()) {
282 tableIdx = this->rockTableIdx_[globalSpaceIdx];
284 return this->rockParams_[tableIdx].compressibility;
287template<
class Gr
idView,
class Flu
idSystem>
288typename FlowGenericProblem<GridView,FluidSystem>::Scalar
292 if (this->rockParams_.empty())
295 unsigned tableIdx = 0;
296 if (!this->rockTableIdx_.empty()) {
297 tableIdx = this->rockTableIdx_[globalSpaceIdx];
299 return this->rockParams_[tableIdx].referencePressure;
302template<
class Gr
idView,
class Flu
idSystem>
303typename FlowGenericProblem<GridView,FluidSystem>::Scalar
305porosity(
unsigned globalSpaceIdx,
unsigned timeIdx)
const
307 return this->referencePorosity_[timeIdx][globalSpaceIdx];
310template<
class Gr
idView,
class Flu
idSystem>
311typename FlowGenericProblem<GridView,FluidSystem>::Scalar
313rockFraction(
unsigned elementIdx,
unsigned timeIdx)
const
319 auto porosity = this->lookUpData_.fieldPropDouble(eclState_.fieldProps(),
"PORO", elementIdx);
320 return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity);
323template<
class Gr
idView,
class Flu
idSystem>
326updateNum(
const std::string& name, std::vector<T>& numbers, std::size_t num_regions)
328 if (!eclState_.fieldProps().has_int(name))
331 std::function<void(T,
int)> valueCheck = [num_regions,name](T fieldPropValue, [[maybe_unused]]
int fieldPropIdx) {
332 if ( fieldPropValue > (
int)num_regions) {
333 throw std::runtime_error(
"Values larger than maximum number of regions "
334 + std::to_string(num_regions) +
" provided in " + name);
336 if ( fieldPropValue <= 0) {
337 throw std::runtime_error(
"zero or negative values provided for region array: " + name);
341 numbers = this->lookUpData_.template assignFieldPropsIntOnLeaf<T>(eclState_.fieldProps(), name,
345template<
class Gr
idView,
class Flu
idSystem>
346void FlowGenericProblem<GridView,FluidSystem>::
349 const auto num_regions = eclState_.getTableManager().getTabdims().getNumPVTTables();
350 updateNum(
"PVTNUM", pvtnum_, num_regions);
353template<
class Gr
idView,
class Flu
idSystem>
354void FlowGenericProblem<GridView,FluidSystem>::
357 const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
358 updateNum(
"SATNUM", satnum_, num_regions);
361template<
class Gr
idView,
class Flu
idSystem>
362void FlowGenericProblem<GridView,FluidSystem>::
365 const auto num_regions = 1;
366 updateNum(
"MISCNUM", miscnum_, num_regions);
369template<
class Gr
idView,
class Flu
idSystem>
370void FlowGenericProblem<GridView,FluidSystem>::
373 const auto num_regions = 1;
374 updateNum(
"PLMIXNUM", plmixnum_, num_regions);
377template<
class Gr
idView,
class Flu
idSystem>
378bool FlowGenericProblem<GridView,FluidSystem>::
379vapparsActive(
int episodeIdx)
const
381 const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
382 return (oilVaporizationControl.getType() == OilVaporizationProperties::OilVaporization::VAPPARS);
385template<
class Gr
idView,
class Flu
idSystem>
386bool FlowGenericProblem<GridView,FluidSystem>::
387beginEpisode_(
bool enableExperiments,
390 if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
393 std::ostringstream ss;
394 boost::posix_time::time_facet* facet =
new boost::posix_time::time_facet(
"%d-%b-%Y");
395 boost::posix_time::ptime curDateTime =
396 boost::posix_time::from_time_t(schedule_.simTime(episodeIdx));
397 ss.imbue(std::locale(std::locale::classic(), facet));
398 ss <<
"Report step " << episodeIdx + 1
399 <<
"/" << schedule_.size() - 1
400 <<
" at day " << schedule_.seconds(episodeIdx)/(24*3600)
401 <<
"/" << schedule_.seconds(schedule_.size() - 1)/(24*3600)
402 <<
", date = " << curDateTime.date()
404 OpmLog::info(ss.str());
407 const auto& events = schedule_[episodeIdx].events();
410 if (episodeIdx > 0 && enableTuning_ && events.hasEvent(ScheduleEvents::TUNING_CHANGE))
412 const auto& sched_state = schedule_[episodeIdx];
413 const auto& tuning = sched_state.tuning();
414 initialTimeStepSize_ = sched_state.max_next_tstep(enableTuning_);
415 maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
422template<
class Gr
idView,
class Flu
idSystem>
423void FlowGenericProblem<GridView,FluidSystem>::
424beginTimeStep_(
bool enableExperiments,
432 if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
433 std::ostringstream ss;
434 boost::posix_time::time_facet* facet =
new boost::posix_time::time_facet(
"%d-%b-%Y");
435 boost::posix_time::ptime date = boost::posix_time::from_time_t(startTime) + boost::posix_time::milliseconds(
static_cast<long long>(time / prefix::milli));
436 ss.imbue(std::locale(std::locale::classic(), facet));
437 ss <<
"\nTime step " << timeStepIndex <<
", stepsize "
438 << unit::convert::to(timeStepSize, unit::day) <<
" days,"
439 <<
" at day " << (double)unit::convert::to(time, unit::day)
440 <<
"/" << (double)unit::convert::to(endTime, unit::day)
441 <<
", date = " << date;
442 OpmLog::info(ss.str());
446template<
class Gr
idView,
class Flu
idSystem>
447void FlowGenericProblem<GridView,FluidSystem>::
450 FluidSystem::initFromState(eclState_, schedule_);
453template<
class Gr
idView,
class Flu
idSystem>
454void FlowGenericProblem<GridView,FluidSystem>::
455readBlackoilExtentionsInitialConditions_(std::size_t numDof,
458 bool enablePolymerMolarWeight,
461 auto getArray = [](
const std::vector<double>& input)
463 if constexpr (std::is_same_v<Scalar,double>) {
466 return std::vector<Scalar>{input.begin(), input.end()};
471 if (eclState_.fieldProps().has_double(
"SSOL")) {
472 solventSaturation_ = getArray(eclState_.fieldProps().get_double(
"SSOL"));
474 solventSaturation_.resize(numDof, 0.0);
477 solventRsw_.resize(numDof, 0.0);
481 if (eclState_.fieldProps().has_double(
"SPOLY")) {
482 polymer_.concentration = getArray(eclState_.fieldProps().get_double(
"SPOLY"));
484 polymer_.concentration.resize(numDof, 0.0);
488 if (enablePolymerMolarWeight) {
489 if (eclState_.fieldProps().has_double(
"SPOLYMW")) {
490 polymer_.moleWeight = getArray(eclState_.fieldProps().get_double(
"SPOLYMW"));
492 polymer_.moleWeight.resize(numDof, 0.0);
497 if (eclState_.fieldProps().has_double(
"SMICR")) {
498 micp_.microbialConcentration = getArray(eclState_.fieldProps().get_double(
"SMICR"));
500 micp_.microbialConcentration.resize(numDof, 0.0);
502 if (eclState_.fieldProps().has_double(
"SOXYG")) {
503 micp_.oxygenConcentration = getArray(eclState_.fieldProps().get_double(
"SOXYG"));
505 micp_.oxygenConcentration.resize(numDof, 0.0);
507 if (eclState_.fieldProps().has_double(
"SUREA")) {
508 micp_.ureaConcentration = getArray(eclState_.fieldProps().get_double(
"SUREA"));
510 micp_.ureaConcentration.resize(numDof, 0.0);
512 if (eclState_.fieldProps().has_double(
"SBIOF")) {
513 micp_.biofilmConcentration = getArray(eclState_.fieldProps().get_double(
"SBIOF"));
515 micp_.biofilmConcentration.resize(numDof, 0.0);
517 if (eclState_.fieldProps().has_double(
"SCALC")) {
518 micp_.calciteConcentration = getArray(eclState_.fieldProps().get_double(
"SCALC"));
520 micp_.calciteConcentration.resize(numDof, 0.0);
525template<
class Gr
idView,
class Flu
idSystem>
526typename FlowGenericProblem<GridView,FluidSystem>::Scalar
530 if (maxWaterSaturation_.empty())
533 return maxWaterSaturation_[globalDofIdx];
536template<
class Gr
idView,
class Flu
idSystem>
537typename FlowGenericProblem<GridView,FluidSystem>::Scalar
541 if (minRefPressure_.empty())
544 return minRefPressure_[globalDofIdx];
547template<
class Gr
idView,
class Flu
idSystem>
548typename FlowGenericProblem<GridView,FluidSystem>::Scalar
552 if (overburdenPressure_.empty())
555 return overburdenPressure_[elementIdx];
558template<
class Gr
idView,
class Flu
idSystem>
559typename FlowGenericProblem<GridView,FluidSystem>::Scalar
563 if (solventSaturation_.empty())
566 return solventSaturation_[elemIdx];
569template<
class Gr
idView,
class Flu
idSystem>
570typename FlowGenericProblem<GridView,FluidSystem>::Scalar
574 if (solventRsw_.empty())
577 return solventRsw_[elemIdx];
582template<
class Gr
idView,
class Flu
idSystem>
583typename FlowGenericProblem<GridView,FluidSystem>::Scalar
587 if (polymer_.concentration.empty()) {
591 return polymer_.concentration[elemIdx];
594template<
class Gr
idView,
class Flu
idSystem>
595typename FlowGenericProblem<GridView,FluidSystem>::Scalar
599 if (polymer_.moleWeight.empty()) {
603 return polymer_.moleWeight[elemIdx];
606template<
class Gr
idView,
class Flu
idSystem>
607typename FlowGenericProblem<GridView,FluidSystem>::Scalar
611 if (micp_.microbialConcentration.empty()) {
618template<
class Gr
idView,
class Flu
idSystem>
619typename FlowGenericProblem<GridView,FluidSystem>::Scalar
623 if (micp_.oxygenConcentration.empty()) {
630template<
class Gr
idView,
class Flu
idSystem>
631typename FlowGenericProblem<GridView,FluidSystem>::Scalar
635 if (micp_.ureaConcentration.empty()) {
642template<
class Gr
idView,
class Flu
idSystem>
643typename FlowGenericProblem<GridView,FluidSystem>::Scalar
647 if (micp_.biofilmConcentration.empty()) {
654template<
class Gr
idView,
class Flu
idSystem>
655typename FlowGenericProblem<GridView,FluidSystem>::Scalar
659 if (micp_.calciteConcentration.empty()) {
666template<
class Gr
idView,
class Flu
idSystem>
673 return pvtnum_[elemIdx];
676template<
class Gr
idView,
class Flu
idSystem>
683 return satnum_[elemIdx];
686template<
class Gr
idView,
class Flu
idSystem>
690 if (miscnum_.empty())
693 return miscnum_[elemIdx];
696template<
class Gr
idView,
class Flu
idSystem>
700 if (plmixnum_.empty())
703 return plmixnum_[elemIdx];
706template<
class Gr
idView,
class Flu
idSystem>
707typename FlowGenericProblem<GridView,FluidSystem>::Scalar
711 if (polymer_.maxAdsorption.empty()) {
715 return polymer_.maxAdsorption[elemIdx];
718template<
class Gr
idView,
class Flu
idSystem>
722 return this->maxWaterSaturation_ == rhs.maxWaterSaturation_ &&
723 this->minRefPressure_ == rhs.minRefPressure_ &&
724 this->overburdenPressure_ == rhs.overburdenPressure_ &&
725 this->solventSaturation_ == rhs.solventSaturation_ &&
726 this->solventRsw_ == rhs.solventRsw_ &&
727 this->polymer_ == rhs.polymer_ &&
728 this->micp_ == rhs.micp_;