244 throw std::runtime_error(
"ROCKCOMP is activated." + std::to_string(numRocktabTables)
245 +
" ROCKTAB tables is expected, but " + std::to_string(rocktabTables.size()) +
" is provided");
247 rockCompPoroMult_.resize(numRocktabTables);
248 rockCompTransMult_.resize(numRocktabTables);
249 for (std::size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
250 const auto& rocktabTable = rocktabTables.template getTable<RocktabTable>(regionIdx);
251 const auto& pressureColumn = rocktabTable.getPressureColumn();
252 const auto& poroColumn = rocktabTable.getPoreVolumeMultiplierColumn();
253 const auto& transColumn = rocktabTable.getTransmissibilityMultiplierColumn();
254 rockCompPoroMult_[regionIdx].setXYContainers(pressureColumn, poroColumn);
255 rockCompTransMult_[regionIdx].setXYContainers(pressureColumn, transColumn);
258 const auto& rock2dTables = eclState_.getTableManager().getRock2dTables();
259 const auto& rock2dtrTables = eclState_.getTableManager().getRock2dtrTables();
260 const auto& rockwnodTables = eclState_.getTableManager().getRockwnodTables();
261 maxWaterSaturation_.resize(numElem, 0.0);
263 if (rock2dTables.size() != numRocktabTables)
264 throw std::runtime_error(
"Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables)
265 +
" ROCK2D tables is expected, but " + std::to_string(rock2dTables.size()) +
" is provided");
267 if (rockwnodTables.size() != numRocktabTables)
268 throw std::runtime_error(
"Water compation option is selected in ROCKCOMP." + std::to_string(numRocktabTables)
269 +
" ROCKWNOD tables is expected, but " + std::to_string(rockwnodTables.size()) +
" is provided");
271 rockCompPoroMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical));
272 for (std::size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
273 const RockwnodTable& rockwnodTable = rockwnodTables.template getTable<RockwnodTable>(regionIdx);
274 const auto& rock2dTable = rock2dTables[regionIdx];
276 if (rockwnodTable.getSaturationColumn().size() != rock2dTable.sizeMultValues())
277 throw std::runtime_error(
"Number of entries in ROCKWNOD and ROCK2D needs to match.");
279 for (std::size_t xIdx = 0; xIdx < rock2dTable.size(); ++xIdx) {
280 rockCompPoroMultWc_[regionIdx].appendXPos(rock2dTable.getPressureValue(xIdx));
281 for (std::size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx)
282 rockCompPoroMultWc_[regionIdx].appendSamplePoint(xIdx,
283 rockwnodTable.getSaturationColumn()[yIdx],
284 rock2dTable.getPvmultValue(xIdx, yIdx));
288 if (!rock2dtrTables.empty()) {
289 rockCompTransMultWc_.resize(numRocktabTables, TabulatedTwoDFunction(TabulatedTwoDFunction::InterpolationPolicy::Vertical));
290 for (std::size_t regionIdx = 0; regionIdx < numRocktabTables; ++regionIdx) {
291 const RockwnodTable& rockwnodTable = rockwnodTables.template getTable<RockwnodTable>(regionIdx);
292 const auto& rock2dtrTable = rock2dtrTables[regionIdx];
294 if (rockwnodTable.getSaturationColumn().size() != rock2dtrTable.sizeMultValues())
295 throw std::runtime_error(
"Number of entries in ROCKWNOD and ROCK2DTR needs to match.");
297 for (std::size_t xIdx = 0; xIdx < rock2dtrTable.size(); ++xIdx) {
298 rockCompTransMultWc_[regionIdx].appendXPos(rock2dtrTable.getPressureValue(xIdx));
299 for (std::size_t yIdx = 0; yIdx < rockwnodTable.getSaturationColumn().size(); ++yIdx)
300 rockCompTransMultWc_[regionIdx].appendSamplePoint(xIdx,
301 rockwnodTable.getSaturationColumn()[yIdx],
302 rock2dtrTable.getTransMultValue(xIdx, yIdx));
309template<
class Gr
idView,
class Flu
idSystem>
310typename FlowGenericProblem<GridView,FluidSystem>::Scalar
314 if (this->rockParams_.empty())
317 unsigned tableIdx = 0;
318 if (!this->rockTableIdx_.empty()) {
319 tableIdx = this->rockTableIdx_[globalSpaceIdx];
321 return this->rockParams_[tableIdx].compressibility;
324template<
class Gr
idView,
class Flu
idSystem>
325typename FlowGenericProblem<GridView,FluidSystem>::Scalar
327porosity(
unsigned globalSpaceIdx,
unsigned timeIdx)
const
329 return this->referencePorosity_[timeIdx][globalSpaceIdx];
332template<
class Gr
idView,
class Flu
idSystem>
333typename FlowGenericProblem<GridView,FluidSystem>::Scalar
335rockFraction(
unsigned elementIdx,
unsigned timeIdx)
const
341 auto porosity = this->lookUpData_.fieldPropDouble(eclState_.fieldProps(),
"PORO", elementIdx);
342 return referencePorosity(elementIdx, timeIdx) / porosity * (1 - porosity);
345template<
class Gr
idView,
class Flu
idSystem>
348updateNum(
const std::string& name, std::vector<T>& numbers, std::size_t num_regions)
350 if (!eclState_.fieldProps().has_int(name))
353 std::function<void(T,
int)> valueCheck = [num_regions,name](T fieldPropValue, [[maybe_unused]]
int fieldPropIdx) {
354 if ( fieldPropValue > (
int)num_regions) {
355 throw std::runtime_error(
"Values larger than maximum number of regions "
356 + std::to_string(num_regions) +
" provided in " + name);
358 if ( fieldPropValue <= 0) {
359 throw std::runtime_error(
"zero or negative values provided for region array: " + name);
363 numbers = this->lookUpData_.template assignFieldPropsIntOnLeaf<T>(eclState_.fieldProps(), name,
367template<
class Gr
idView,
class Flu
idSystem>
368void FlowGenericProblem<GridView,FluidSystem>::
371 const auto num_regions = eclState_.getTableManager().getTabdims().getNumPVTTables();
372 updateNum(
"PVTNUM", pvtnum_, num_regions);
375template<
class Gr
idView,
class Flu
idSystem>
376void FlowGenericProblem<GridView,FluidSystem>::
379 const auto num_regions = eclState_.getTableManager().getTabdims().getNumSatTables();
380 updateNum(
"SATNUM", satnum_, num_regions);
383template<
class Gr
idView,
class Flu
idSystem>
384void FlowGenericProblem<GridView,FluidSystem>::
387 const auto num_regions = 1;
388 updateNum(
"MISCNUM", miscnum_, num_regions);
391template<
class Gr
idView,
class Flu
idSystem>
392void FlowGenericProblem<GridView,FluidSystem>::
395 const auto num_regions = 1;
396 updateNum(
"PLMIXNUM", plmixnum_, num_regions);
399template<
class Gr
idView,
class Flu
idSystem>
400bool FlowGenericProblem<GridView,FluidSystem>::
401vapparsActive(
int episodeIdx)
const
403 const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
404 return (oilVaporizationControl.getType() == OilVaporizationProperties::OilVaporization::VAPPARS);
407template<
class Gr
idView,
class Flu
idSystem>
408bool FlowGenericProblem<GridView,FluidSystem>::
409beginEpisode_(
bool enableExperiments,
412 if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
415 std::ostringstream ss;
416 boost::posix_time::time_facet* facet =
new boost::posix_time::time_facet(
"%d-%b-%Y");
417 boost::posix_time::ptime curDateTime =
418 boost::posix_time::from_time_t(schedule_.simTime(episodeIdx));
419 ss.imbue(std::locale(std::locale::classic(), facet));
420 ss <<
"Report step " << episodeIdx + 1
421 <<
"/" << schedule_.size() - 1
422 <<
" at day " << schedule_.seconds(episodeIdx)/(24*3600)
423 <<
"/" << schedule_.seconds(schedule_.size() - 1)/(24*3600)
424 <<
", date = " << curDateTime.date()
426 OpmLog::info(ss.str());
429 const auto& events = schedule_[episodeIdx].events();
432 if (episodeIdx > 0 && enableTuning_ && events.hasEvent(ScheduleEvents::TUNING_CHANGE))
434 const auto& sched_state = schedule_[episodeIdx];
435 const auto& tuning = sched_state.tuning();
436 initialTimeStepSize_ = sched_state.max_next_tstep(enableTuning_);
437 maxTimeStepAfterWellEvent_ = tuning.TMAXWC;
444template<
class Gr
idView,
class Flu
idSystem>
445void FlowGenericProblem<GridView,FluidSystem>::
446beginTimeStep_(
bool enableExperiments,
454 if (enableExperiments && gridView_.comm().rank() == 0 && episodeIdx >= 0) {
455 std::ostringstream ss;
456 boost::posix_time::time_facet* facet =
new boost::posix_time::time_facet(
"%d-%b-%Y");
457 boost::posix_time::ptime date = boost::posix_time::from_time_t(startTime) + boost::posix_time::milliseconds(
static_cast<long long>(time / prefix::milli));
458 ss.imbue(std::locale(std::locale::classic(), facet));
459 ss <<
"\nTime step " << timeStepIndex <<
", stepsize "
460 << unit::convert::to(timeStepSize, unit::day) <<
" days,"
461 <<
" at day " << (double)unit::convert::to(time, unit::day)
462 <<
"/" << (double)unit::convert::to(endTime, unit::day)
463 <<
", date = " << date;
464 OpmLog::info(ss.str());
468template<
class Gr
idView,
class Flu
idSystem>
469void FlowGenericProblem<GridView,FluidSystem>::
472 FluidSystem::initFromState(eclState_, schedule_);
475template<
class Gr
idView,
class Flu
idSystem>
476void FlowGenericProblem<GridView,FluidSystem>::
477readBlackoilExtentionsInitialConditions_(std::size_t numDof,
480 bool enablePolymerMolarWeight,
483 auto getArray = [](
const std::vector<double>& input)
485 if constexpr (std::is_same_v<Scalar,double>) {
488 return std::vector<Scalar>{input.begin(), input.end()};
493 if (eclState_.fieldProps().has_double(
"SSOL")) {
494 solventSaturation_ = getArray(eclState_.fieldProps().get_double(
"SSOL"));
496 solventSaturation_.resize(numDof, 0.0);
499 solventRsw_.resize(numDof, 0.0);
503 if (eclState_.fieldProps().has_double(
"SPOLY")) {
504 polymer_.concentration = getArray(eclState_.fieldProps().get_double(
"SPOLY"));
506 polymer_.concentration.resize(numDof, 0.0);
510 if (enablePolymerMolarWeight) {
511 if (eclState_.fieldProps().has_double(
"SPOLYMW")) {
512 polymer_.moleWeight = getArray(eclState_.fieldProps().get_double(
"SPOLYMW"));
514 polymer_.moleWeight.resize(numDof, 0.0);
519 if (eclState_.fieldProps().has_double(
"SMICR")) {
520 micp_.microbialConcentration = getArray(eclState_.fieldProps().get_double(
"SMICR"));
522 micp_.microbialConcentration.resize(numDof, 0.0);
524 if (eclState_.fieldProps().has_double(
"SOXYG")) {
525 micp_.oxygenConcentration = getArray(eclState_.fieldProps().get_double(
"SOXYG"));
527 micp_.oxygenConcentration.resize(numDof, 0.0);
529 if (eclState_.fieldProps().has_double(
"SUREA")) {
530 micp_.ureaConcentration = getArray(eclState_.fieldProps().get_double(
"SUREA"));
532 micp_.ureaConcentration.resize(numDof, 0.0);
534 if (eclState_.fieldProps().has_double(
"SBIOF")) {
535 micp_.biofilmConcentration = getArray(eclState_.fieldProps().get_double(
"SBIOF"));
537 micp_.biofilmConcentration.resize(numDof, 0.0);
539 if (eclState_.fieldProps().has_double(
"SCALC")) {
540 micp_.calciteConcentration = getArray(eclState_.fieldProps().get_double(
"SCALC"));
542 micp_.calciteConcentration.resize(numDof, 0.0);
547template<
class Gr
idView,
class Flu
idSystem>
548typename FlowGenericProblem<GridView,FluidSystem>::Scalar
552 if (maxWaterSaturation_.empty())
555 return maxWaterSaturation_[globalDofIdx];
558template<
class Gr
idView,
class Flu
idSystem>
559typename FlowGenericProblem<GridView,FluidSystem>::Scalar
563 if (minRefPressure_.empty())
566 return minRefPressure_[globalDofIdx];
569template<
class Gr
idView,
class Flu
idSystem>
570typename FlowGenericProblem<GridView,FluidSystem>::Scalar
574 if (overburdenPressure_.empty())
577 return overburdenPressure_[elementIdx];
580template<
class Gr
idView,
class Flu
idSystem>
581typename FlowGenericProblem<GridView,FluidSystem>::Scalar
585 if (solventSaturation_.empty())
588 return solventSaturation_[elemIdx];
591template<
class Gr
idView,
class Flu
idSystem>
592typename FlowGenericProblem<GridView,FluidSystem>::Scalar
596 if (solventRsw_.empty())
599 return solventRsw_[elemIdx];
604template<
class Gr
idView,
class Flu
idSystem>
605typename FlowGenericProblem<GridView,FluidSystem>::Scalar
609 if (polymer_.concentration.empty()) {
613 return polymer_.concentration[elemIdx];
616template<
class Gr
idView,
class Flu
idSystem>
617typename FlowGenericProblem<GridView,FluidSystem>::Scalar
621 if (polymer_.moleWeight.empty()) {
625 return polymer_.moleWeight[elemIdx];
628template<
class Gr
idView,
class Flu
idSystem>
629typename FlowGenericProblem<GridView,FluidSystem>::Scalar
633 if (micp_.microbialConcentration.empty()) {
640template<
class Gr
idView,
class Flu
idSystem>
641typename FlowGenericProblem<GridView,FluidSystem>::Scalar
645 if (micp_.oxygenConcentration.empty()) {
652template<
class Gr
idView,
class Flu
idSystem>
653typename FlowGenericProblem<GridView,FluidSystem>::Scalar
657 if (micp_.ureaConcentration.empty()) {
664template<
class Gr
idView,
class Flu
idSystem>
665typename FlowGenericProblem<GridView,FluidSystem>::Scalar
669 if (micp_.biofilmConcentration.empty()) {
676template<
class Gr
idView,
class Flu
idSystem>
677typename FlowGenericProblem<GridView,FluidSystem>::Scalar
681 if (micp_.calciteConcentration.empty()) {
688template<
class Gr
idView,
class Flu
idSystem>
695 return pvtnum_[elemIdx];
698template<
class Gr
idView,
class Flu
idSystem>
705 return satnum_[elemIdx];
708template<
class Gr
idView,
class Flu
idSystem>
712 if (miscnum_.empty())
715 return miscnum_[elemIdx];
718template<
class Gr
idView,
class Flu
idSystem>
722 if (plmixnum_.empty())
725 return plmixnum_[elemIdx];
728template<
class Gr
idView,
class Flu
idSystem>
729typename FlowGenericProblem<GridView,FluidSystem>::Scalar
733 if (polymer_.maxAdsorption.empty()) {
737 return polymer_.maxAdsorption[elemIdx];
740template<
class Gr
idView,
class Flu
idSystem>
744 return this->maxWaterSaturation_ == rhs.maxWaterSaturation_ &&
745 this->minRefPressure_ == rhs.minRefPressure_ &&
746 this->overburdenPressure_ == rhs.overburdenPressure_ &&
747 this->solventSaturation_ == rhs.solventSaturation_ &&
748 this->solventRsw_ == rhs.solventRsw_ &&
749 this->polymer_ == rhs.polymer_ &&
750 this->micp_ == rhs.micp_;