106 GetPropType<TypeTag, Properties::EquilGrid>,
107 GetPropType<TypeTag, Properties::GridView>,
108 GetPropType<TypeTag, Properties::ElementMapper>,
109 GetPropType<TypeTag, Properties::Scalar>>
119 using Element =
typename GridView::template Codim<0>::Entity;
121 using ElementIterator =
typename GridView::template Codim<0>::Iterator;
124 typedef Dune::MultipleCodimMultipleGeomTypeMapper< GridView > VertexMapper;
133 static void registerParameters()
137 Parameters::Register<Parameters::EnableAsyncEclOutput>
138 (
"Write the ECL-formated results in a non-blocking way "
139 "(i.e., using a separate thread).");
140 Parameters::Register<Parameters::EnableEsmry>
141 (
"Write ESMRY file for fast loading of summary data.");
148 :
BaseType(simulator.vanguard().schedule(),
149 simulator.vanguard().eclState(),
150 simulator.vanguard().summaryConfig(),
151 simulator.vanguard().grid(),
152 ((simulator.vanguard().grid().comm().rank() == 0)
153 ? &simulator.vanguard().equilGrid()
155 simulator.vanguard().gridView(),
156 simulator.vanguard().cartesianIndexMapper(),
157 ((simulator.vanguard().grid().comm().rank() == 0)
158 ? &simulator.vanguard().equilCartesianIndexMapper()
160 Parameters::Get<Parameters::EnableAsyncEclOutput>(),
161 Parameters::Get<Parameters::EnableEsmry>())
162 , simulator_(simulator)
165 if (this->simulator_.vanguard().grid().comm().size() > 1) {
166 auto smryCfg = (this->simulator_.vanguard().grid().comm().rank() == 0)
167 ? this->eclIO_->finalSummaryConfig()
170 eclBroadcast(this->simulator_.vanguard().grid().comm(), smryCfg);
172 this->outputModule_ = std::make_unique<OutputBlackOilModule<TypeTag>>
173 (simulator, smryCfg, this->collectOnIORank_);
178 this->outputModule_ = std::make_unique<OutputBlackOilModule<TypeTag>>
179 (simulator, this->eclIO_->finalSummaryConfig(), this->collectOnIORank_);
182 this->rank_ = this->simulator_.vanguard().grid().comm().rank();
184 this->simulator_.vanguard().eclState().computeFipRegionStatistics();
190 const EquilGrid& globalGrid()
const
192 return simulator_.vanguard().equilGrid();
201 const int reportStepNum = simulator_.episodeIndex() + 1;
221 if (reportStepNum == 0)
224 const Scalar curTime = simulator_.time() + simulator_.timeStepSize();
225 const Scalar totalCpuTime =
226 simulator_.executionTimer().realTimeElapsed() +
227 simulator_.setupTimer().realTimeElapsed() +
228 simulator_.vanguard().setupTime();
230 const auto localWellData = simulator_.problem().wellModel().wellData();
231 const auto localWBP = simulator_.problem().wellModel().wellBlockAveragePressures();
232 const auto localGroupAndNetworkData = simulator_.problem().wellModel()
233 .groupAndNetworkData(reportStepNum);
235 const auto localAquiferData = simulator_.problem().aquiferModel().aquiferData();
236 const auto localWellTestState = simulator_.problem().wellModel().wellTestState();
237 this->prepareLocalCellData(isSubStep, reportStepNum);
239 if (this->outputModule_->needInterfaceFluxes(isSubStep)) {
240 this->captureLocalFluxData();
243 if (this->collectOnIORank_.isParallel()) {
244 OPM_BEGIN_PARALLEL_TRY_CATCH()
246 this->collectOnIORank_.collect({},
247 outputModule_->getBlockData(),
250 localGroupAndNetworkData,
253 this->outputModule_->getInterRegFlows(),
257 if (this->collectOnIORank_.isIORank()) {
258 auto& iregFlows = this->collectOnIORank_.globalInterRegFlows();
260 if (! iregFlows.readIsConsistent()) {
261 throw std::runtime_error {
262 "Inconsistent inter-region flow "
263 "region set names in parallel"
270 OPM_END_PARALLEL_TRY_CATCH(
"Collect to I/O rank: ",
271 this->simulator_.vanguard().grid().comm());
275 std::map<std::string, double> miscSummaryData;
276 std::map<std::string, std::vector<double>> regionData;
280 OPM_TIMEBLOCK(outputFipLogAndFipresvLog);
282 inplace = outputModule_->calc_inplace(miscSummaryData, regionData, simulator_.gridView().comm());
284 if (this->collectOnIORank_.isIORank()){
290 if (totalCpuTime != 0.0) {
291 miscSummaryData[
"TCPU"] = totalCpuTime;
293 if (this->sub_step_report_.total_newton_iterations != 0) {
294 miscSummaryData[
"NEWTON"] = this->sub_step_report_.total_newton_iterations;
296 if (this->sub_step_report_.total_linear_iterations != 0) {
297 miscSummaryData[
"MLINEARS"] = this->sub_step_report_.total_linear_iterations;
299 if (this->sub_step_report_.total_newton_iterations != 0) {
300 miscSummaryData[
"NLINEARS"] =
static_cast<float>(this->sub_step_report_.total_linear_iterations) / this->sub_step_report_.total_newton_iterations;
302 if (this->sub_step_report_.min_linear_iterations != std::numeric_limits<unsigned int>::max()) {
303 miscSummaryData[
"NLINSMIN"] = this->sub_step_report_.min_linear_iterations;
305 if (this->sub_step_report_.max_linear_iterations != 0) {
306 miscSummaryData[
"NLINSMAX"] = this->sub_step_report_.max_linear_iterations;
308 if (this->simulation_report_.success.total_newton_iterations != 0) {
309 miscSummaryData[
"MSUMLINS"] = this->simulation_report_.success.total_linear_iterations;
311 if (this->simulation_report_.success.total_newton_iterations != 0) {
312 miscSummaryData[
"MSUMNEWT"] = this->simulation_report_.success.total_newton_iterations;
316 OPM_TIMEBLOCK(evalSummary);
318 const auto& blockData = this->collectOnIORank_.isParallel()
319 ? this->collectOnIORank_.globalBlockData()
320 : this->outputModule_->getBlockData();
322 const auto& interRegFlows = this->collectOnIORank_.isParallel()
323 ? this->collectOnIORank_.globalInterRegFlows()
324 : this->outputModule_->getInterRegFlows();
326 this->evalSummary(reportStepNum,
330 localGroupAndNetworkData,
336 this->outputModule_->initialInplace(),
338 this->summaryState(),
346 const auto& fip = simulator_.vanguard().eclState().getEclipseConfig().fip();
347 if (!fip.output(FIPConfig::OutputField::FIELD) &&
348 !fip.output(FIPConfig::OutputField::RESV)) {
352 const auto& gridView = simulator_.vanguard().gridView();
353 const int num_interior = detail::
354 countLocalInteriorCellsGridView(gridView);
356 this->outputModule_->
357 allocBuffers(num_interior, 0,
false,
false,
false);
360#pragma omp parallel for
362 for (
int dofIdx = 0; dofIdx < num_interior; ++dofIdx) {
363 const auto& intQuants = *simulator_.model().cachedIntensiveQuantities(dofIdx, 0);
364 const auto totVolume = simulator_.model().dofTotalVolume(dofIdx);
366 this->outputModule_->updateFluidInPlace(dofIdx, intQuants, totVolume);
369 std::map<std::string, double> miscSummaryData;
370 std::map<std::string, std::vector<double>> regionData;
373 OPM_TIMEBLOCK(outputFipLogAndFipresvLog);
375 boost::posix_time::ptime start_time = boost::posix_time::from_time_t(simulator_.vanguard().schedule().getStartTime());
377 inplace = outputModule_->calc_inplace(miscSummaryData, regionData, simulator_.gridView().comm());
379 if (this->collectOnIORank_.isIORank()){
382 outputModule_->outputFipAndResvLog(inplace_, 0, 0.0, start_time,
383 false, simulator_.gridView().comm());
391 if ((rstep > 0) && (this->collectOnIORank_.isIORank())){
393 const auto& rpt = this->schedule_[rstep-1].rpt_config.get();
394 if (rpt.contains(
"WELLS") && rpt.at(
"WELLS") > 0) {
396 outputModule_->outputProdLog(rstep-1);
397 outputModule_->outputInjLog(rstep-1);
398 outputModule_->outputCumLog(rstep-1);
409 void writeOutput(data::Solution&& localCellData,
bool isSubStep)
411 OPM_TIMEBLOCK(writeOutput);
413 const int reportStepNum = simulator_.episodeIndex() + 1;
414 this->prepareLocalCellData(isSubStep, reportStepNum);
415 this->outputModule_->outputErrorLog(simulator_.gridView().comm());
418 auto localWellData = simulator_.problem().wellModel().wellData();
419 auto localGroupAndNetworkData = simulator_.problem().wellModel()
420 .groupAndNetworkData(reportStepNum);
422 auto localAquiferData = simulator_.problem().aquiferModel().aquiferData();
423 auto localWellTestState = simulator_.problem().wellModel().wellTestState();
425 const bool isFlowsn = this->outputModule_->hasFlowsn();
426 auto flowsn = this->outputModule_->getFlowsn();
428 const bool isFloresn = this->outputModule_->hasFloresn();
429 auto floresn = this->outputModule_->getFloresn();
432 if (! isSubStep || Parameters::Get<Parameters::EnableWriteAllSolutions>()) {
434 if (localCellData.empty()) {
435 this->outputModule_->assignToSolution(localCellData);
439 this->outputModule_->addRftDataToWells(localWellData, reportStepNum);
442 if (this->collectOnIORank_.isParallel() ||
443 this->collectOnIORank_.doesNeedReordering())
450 this->collectOnIORank_.collect(localCellData,
451 this->outputModule_->getBlockData(),
454 localGroupAndNetworkData,
460 if (this->collectOnIORank_.isIORank()) {
461 this->outputModule_->assignGlobalFieldsToSolution(this->collectOnIORank_.globalCellData());
464 this->outputModule_->assignGlobalFieldsToSolution(localCellData);
467 if (this->collectOnIORank_.isIORank()) {
468 const Scalar curTime = simulator_.time() + simulator_.timeStepSize();
469 const Scalar nextStepSize = simulator_.problem().nextTimeStepSize();
470 std::optional<int> timeStepIdx;
471 if (Parameters::Get<Parameters::EnableWriteAllSolutions>()) {
472 timeStepIdx = simulator_.timeStepIndex();
474 this->doWriteOutput(reportStepNum, timeStepIdx, isSubStep,
475 std::move(localCellData),
476 std::move(localWellData),
477 std::move(localGroupAndNetworkData),
478 std::move(localAquiferData),
479 std::move(localWellTestState),
482 this->summaryState(),
483 this->simulator_.problem().thresholdPressure().getRestartVector(),
484 curTime, nextStepSize,
485 Parameters::Get<Parameters::EclOutputDoublePrecision>(),
486 isFlowsn, std::move(flowsn),
487 isFloresn, std::move(floresn));
493 const auto enablePCHysteresis = simulator_.problem().materialLawManager()->enablePCHysteresis();
494 const auto enableNonWettingHysteresis = simulator_.problem().materialLawManager()->enableNonWettingHysteresis();
495 const auto enableWettingHysteresis = simulator_.problem().materialLawManager()->enableWettingHysteresis();
496 const auto oilActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx);
497 const auto gasActive = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
498 const auto waterActive = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx);
499 const auto enableSwatinit = simulator_.vanguard().eclState().fieldProps().has_double(
"SWATINIT");
501 std::vector<RestartKey> solutionKeys {
502 {
"PRESSURE", UnitSystem::measure::pressure},
503 {
"SWAT", UnitSystem::measure::identity, waterActive},
504 {
"SGAS", UnitSystem::measure::identity, gasActive},
505 {
"TEMP", UnitSystem::measure::temperature, enableEnergy},
506 {
"SSOLVENT", UnitSystem::measure::identity, enableSolvent},
508 {
"RS", UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGas()},
509 {
"RV", UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedOil()},
510 {
"RVW", UnitSystem::measure::oil_gas_ratio, FluidSystem::enableVaporizedWater()},
511 {
"RSW", UnitSystem::measure::gas_oil_ratio, FluidSystem::enableDissolvedGasInWater()},
513 {
"SGMAX", UnitSystem::measure::identity, enableNonWettingHysteresis && oilActive && gasActive},
514 {
"SHMAX", UnitSystem::measure::identity, enableWettingHysteresis && oilActive && gasActive},
516 {
"SOMAX", UnitSystem::measure::identity,
517 (enableNonWettingHysteresis && oilActive && waterActive)
518 || simulator_.problem().vapparsActive(simulator_.episodeIndex())},
520 {
"SOMIN", UnitSystem::measure::identity, enablePCHysteresis && oilActive && gasActive},
521 {
"SWHY1", UnitSystem::measure::identity, enablePCHysteresis && oilActive && waterActive},
522 {
"SWMAX", UnitSystem::measure::identity, enableWettingHysteresis && oilActive && waterActive},
524 {
"PPCW", UnitSystem::measure::pressure, enableSwatinit},
528 const auto& tracers = simulator_.vanguard().eclState().tracer();
530 for (
const auto& tracer : tracers) {
531 const auto enableSolTracer =
532 ((tracer.phase == Phase::GAS) && FluidSystem::enableDissolvedGas()) ||
533 ((tracer.phase == Phase::OIL) && FluidSystem::enableVaporizedOil());
535 solutionKeys.emplace_back(tracer.fname(), UnitSystem::measure::identity,
true);
536 solutionKeys.emplace_back(tracer.sname(), UnitSystem::measure::identity, enableSolTracer);
540 const auto& inputThpres = eclState().getSimulationConfig().getThresholdPressure();
541 const std::vector<RestartKey> extraKeys {
542 {
"OPMEXTRA", UnitSystem::measure::identity,
false},
543 {
"THRESHPR", UnitSystem::measure::pressure, inputThpres.active()},
546 const auto& gridView = this->simulator_.vanguard().gridView();
547 const auto numElements = gridView.size(0);
553 const auto restartStepIdx = this->simulator_.vanguard()
554 .eclState().getInitConfig().getRestartStep();
556 this->outputModule_->allocBuffers(numElements,
564 const auto restartValues =
565 loadParallelRestart(this->eclIO_.get(),
567 this->summaryState(),
568 solutionKeys, extraKeys, gridView.comm());
570 for (
auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
571 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
572 this->outputModule_->setRestart(restartValues.solution, elemIdx, globalIdx);
575 auto& tracer_model = simulator_.problem().tracerModel();
576 for (
int tracer_index = 0; tracer_index < tracer_model.numTracers(); ++tracer_index) {
579 const auto& free_tracer_name = tracer_model.fname(tracer_index);
580 const auto& free_tracer_solution = restartValues.solution
581 .template data<double>(free_tracer_name);
583 for (
auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
584 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
585 tracer_model.setFreeTracerConcentration
586 (tracer_index, elemIdx, free_tracer_solution[globalIdx]);
591 if ((tracer_model.phase(tracer_index) == Phase::GAS && FluidSystem::enableDissolvedGas()) ||
592 (tracer_model.phase(tracer_index) == Phase::OIL && FluidSystem::enableVaporizedOil()))
594 tracer_model.setEnableSolTracers(tracer_index,
true);
596 const auto& sol_tracer_name = tracer_model.sname(tracer_index);
597 const auto& sol_tracer_solution = restartValues.solution
598 .template data<double>(sol_tracer_name);
600 for (
auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
601 const auto globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
602 tracer_model.setSolTracerConcentration
603 (tracer_index, elemIdx, sol_tracer_solution[globalIdx]);
607 tracer_model.setEnableSolTracers(tracer_index,
false);
609 for (
auto elemIdx = 0*numElements; elemIdx < numElements; ++elemIdx) {
610 tracer_model.setSolTracerConcentration(tracer_index, elemIdx, 0.0);
615 if (inputThpres.active()) {
616 const_cast<Simulator&
>(this->simulator_)
617 .problem().thresholdPressure()
618 .setFromRestart(restartValues.getExtra(
"THRESHPR"));
621 restartTimeStepSize_ = restartValues.getExtra(
"OPMEXTRA")[0];
622 if (restartTimeStepSize_ <= 0) {
623 restartTimeStepSize_ = std::numeric_limits<double>::max();
627 this->simulator_.problem().wellModel()
628 .initFromRestartFile(restartValues);
630 if (!restartValues.aquifer.empty()) {
631 this->simulator_.problem().mutableAquiferModel()
632 .initFromRestart(restartValues.aquifer);
642 auto miscSummaryData = std::map<std::string, double>{};
643 auto regionData = std::map<std::string, std::vector<double>>{};
653 auto inplace = this->outputModule_
654 ->calc_inplace(miscSummaryData, regionData,
655 this->simulator_.gridView().comm());
657 if (this->collectOnIORank_.isIORank()) {
658 this->inplace_ = std::move(inplace);
662 const OutputBlackOilModule<TypeTag>& outputModule()
const
663 {
return *outputModule_; }
665 OutputBlackOilModule<TypeTag>& mutableOutputModule()
const
666 {
return *outputModule_; }
668 Scalar restartTimeStepSize()
const
669 {
return restartTimeStepSize_; }
671 template <
class Serializer>
672 void serializeOp(Serializer& serializer)
674 serializer(*outputModule_);
678 static bool enableEclOutput_()
680 static bool enable = Parameters::Get<Parameters::EnableEclOutput>();
684 const EclipseState& eclState()
const
685 {
return simulator_.vanguard().eclState(); }
687 SummaryState& summaryState()
688 {
return simulator_.vanguard().summaryState(); }
690 Action::State& actionState()
691 {
return simulator_.vanguard().actionState(); }
694 {
return simulator_.vanguard().udqState(); }
696 const Schedule& schedule()
const
697 {
return simulator_.vanguard().schedule(); }
699 void prepareLocalCellData(
const bool isSubStep,
700 const int reportStepNum)
702 OPM_TIMEBLOCK(prepareLocalCellData);
704 if (this->outputModule_->localDataValid()) {
708 const auto& gridView = simulator_.vanguard().gridView();
709 const bool log = this->collectOnIORank_.isIORank();
711 const int num_interior = detail::
712 countLocalInteriorCellsGridView(gridView);
713 this->outputModule_->
714 allocBuffers(num_interior, reportStepNum,
715 isSubStep && !Parameters::Get<Parameters::EnableWriteAllSolutions>(),
718 ElementContext elemCtx(simulator_);
720 OPM_BEGIN_PARALLEL_TRY_CATCH();
723 OPM_TIMEBLOCK(prepareCellBasedData);
725 this->outputModule_->prepareDensityAccumulation();
727 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
728 elemCtx.updatePrimaryStencil(elem);
729 elemCtx.updatePrimaryIntensiveQuantities(0);
731 this->outputModule_->processElement(elemCtx);
734 this->outputModule_->accumulateDensityParallel();
737 if constexpr (enableMech) {
738 if (simulator_.vanguard().eclState().runspec().mech()) {
739 OPM_TIMEBLOCK(prepareMechData);
740 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
741 elemCtx.updatePrimaryStencil(elem);
742 elemCtx.updatePrimaryIntensiveQuantities(0);
743 outputModule_->processElementMech(elemCtx);
748 if (! this->simulator_.model().linearizer().getFlowsInfo().empty()) {
749 OPM_TIMEBLOCK(prepareFlowsData);
750 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
751 elemCtx.updatePrimaryStencil(elem);
752 elemCtx.updatePrimaryIntensiveQuantities(0);
754 this->outputModule_->processElementFlows(elemCtx);
759 OPM_TIMEBLOCK(prepareBlockData);
760 for (
const auto& elem : elements(gridView, Dune::Partitions::interior)) {
761 elemCtx.updatePrimaryStencil(elem);
762 elemCtx.updatePrimaryIntensiveQuantities(0);
764 this->outputModule_->processElementBlockData(elemCtx);
769 OPM_TIMEBLOCK(prepareFluidInPlace);
772#pragma omp parallel for
774 for (
int dofIdx = 0; dofIdx < num_interior; ++dofIdx) {
775 const auto& intQuants = *simulator_.model().cachedIntensiveQuantities(dofIdx, 0);
776 const auto totVolume = simulator_.model().dofTotalVolume(dofIdx);
778 this->outputModule_->updateFluidInPlace(dofIdx, intQuants, totVolume);
782 this->outputModule_->validateLocalData();
784 OPM_END_PARALLEL_TRY_CATCH(
"EclWriter::prepareLocalCellData() failed: ",
785 this->simulator_.vanguard().grid().comm());
788 void captureLocalFluxData()
790 OPM_TIMEBLOCK(captureLocalData);
792 const auto& gridView = this->simulator_.vanguard().gridView();
793 const auto timeIdx = 0u;
795 auto elemCtx = ElementContext { this->simulator_ };
797 const auto elemMapper = ElementMapper { gridView, Dune::mcmgElementLayout() };
798 const auto activeIndex = [&elemMapper](
const Element& e)
800 return elemMapper.index(e);
803 const auto cartesianIndex = [
this](
const int elemIndex)
805 return this->cartMapper_.cartesianIndex(elemIndex);
808 this->outputModule_->initializeFluxData();
810 OPM_BEGIN_PARALLEL_TRY_CATCH();
812 for (
const auto& elem : elements(gridView, Dune::Partitions::interiorBorder)) {
813 elemCtx.updateStencil(elem);
814 elemCtx.updateIntensiveQuantities(timeIdx);
815 elemCtx.updateExtensiveQuantities(timeIdx);
817 this->outputModule_->processFluxes(elemCtx, activeIndex, cartesianIndex);
820 OPM_END_PARALLEL_TRY_CATCH(
"EclWriter::captureLocalFluxData() failed: ",
821 this->simulator_.vanguard().grid().comm())
823 this->outputModule_->finalizeFluxData();
826 Simulator& simulator_;
827 std::unique_ptr<OutputBlackOilModule<TypeTag> > outputModule_;
828 Scalar restartTimeStepSize_;