96 struct MPI_Comm_Deleter;
116 using ModelParameters =
typename Model::ModelParameters;
117 using SolverParameters =
typename Solver::SolverParameters;
123 : simulator_(simulator)
126 simulator_.vanguard().eclState().getIOConfig(),
127 Parameters::Get<Parameters::SaveStep>(),
128 Parameters::Get<Parameters::LoadStep>(),
129 Parameters::Get<Parameters::SaveFile>(),
130 Parameters::Get<Parameters::LoadFile>())
134 this->terminalOutput_ =
false;
135 if (this->grid().comm().rank() == 0) {
139 [compNames =
typename Model::ComponentName{}](
const int compIdx)
140 {
return std::string_view { compNames.name(compIdx) }; }
143 if (!simulator_.vanguard().eclState().getIOConfig().initOnly()) {
144 this->convergence_output_.
145 startThread(this->simulator_.vanguard().eclState(),
147 R
"(OutputExtraConvergenceInfo (--output-extra-convergence-info))",
159 static void registerParameters()
161 ModelParameters::registerParameters();
162 SolverParameters::registerParameters();
163 TimeStepper::registerParameters();
164 detail::registerSimulatorParameters();
176#ifdef RESERVOIR_COUPLING_ENABLED
177 SimulatorReport
run(SimulatorTimer& timer,
int argc,
char** argv)
179 init(timer, argc, argv);
188 simulator_.model().invalidateAndUpdateIntensiveQuantities(0);
190 while (!timer.
done()) {
191 simulator_.problem().writeReports(timer);
192 bool continue_looping = runStep(timer);
193 if (!continue_looping)
break;
195 simulator_.problem().writeReports(timer);
197#ifdef RESERVOIR_COUPLING_ENABLED
200 if (this->reservoirCouplingMaster_) {
201 this->reservoirCouplingMaster_->sendTerminateAndDisconnect();
203 else if (this->reservoirCouplingSlave_ && !this->reservoirCouplingSlave_->terminated()) {
207 this->reservoirCouplingSlave_->receiveTerminateAndDisconnect();
214#ifdef RESERVOIR_COUPLING_ENABLED
219 bool checkRunningAsReservoirCouplingMaster()
221 for (std::size_t report_step = 0; report_step < this->schedule().size(); ++report_step) {
222 auto rescoup = this->schedule()[report_step].rescoup();
223 auto slave_count = rescoup.slaveCount();
224 auto master_group_count = rescoup.masterGroupCount();
228 if (slave_count > 0) {
231 else if (master_group_count > 0) {
233 throw ReservoirCouplingError(
234 "Inconsistent reservoir coupling master schedule: "
235 "Master group count is greater than 0 but slave count is 0"
243#ifdef RESERVOIR_COUPLING_ENABLED
245 void init(
const SimulatorTimer& timer,
int argc,
char** argv)
247 auto slave_mode = Parameters::Get<Parameters::Slave>();
249 this->reservoirCouplingSlave_ =
250 std::make_unique<ReservoirCouplingSlave<Scalar>>(
252 this->schedule(), timer
254 this->reservoirCouplingSlave_->sendAndReceiveInitialData();
255 this->simulator_.setReservoirCouplingSlave(this->reservoirCouplingSlave_.get());
256 wellModel_().setReservoirCouplingSlave(this->reservoirCouplingSlave_.get());
259 auto master_mode = checkRunningAsReservoirCouplingMaster();
261 this->reservoirCouplingMaster_ =
262 std::make_unique<ReservoirCouplingMaster<Scalar>>(
267 this->simulator_.setReservoirCouplingMaster(this->reservoirCouplingMaster_.get());
268 wellModel_().setReservoirCouplingMaster(this->reservoirCouplingMaster_.get());
272 void init(
const SimulatorTimer& timer)
275 simulator_.setEpisodeIndex(-1);
278 solverTimer_ = std::make_unique<time::StopWatch>();
279 totalTimer_ = std::make_unique<time::StopWatch>();
280 totalTimer_->start();
283 bool enableAdaptive = Parameters::Get<Parameters::EnableAdaptiveTimeStepping>();
284 bool enableTUNING = Parameters::Get<Parameters::EnableTuning>();
285 if (enableAdaptive) {
286 const UnitSystem& unitSystem = this->simulator_.vanguard().eclState().getUnits();
287 const auto& sched_state = schedule()[timer.currentStepNum()];
288 auto max_next_tstep = sched_state.max_next_tstep(enableTUNING);
290 adaptiveTimeStepping_ = std::make_unique<TimeStepper>(max_next_tstep,
291 sched_state.tuning(),
292 unitSystem, report_, terminalOutput_);
295 adaptiveTimeStepping_ = std::make_unique<TimeStepper>(unitSystem, report_, max_next_tstep, terminalOutput_);
300 adaptiveTimeStepping_->setSuggestedNextStep(simulator_.timeStepSize());
305 void updateTUNING(
const Tuning& tuning)
307 modelParam_.tolerance_cnv_ = tuning.TRGCNV;
308 modelParam_.tolerance_cnv_relaxed_ = tuning.XXXCNV;
309 modelParam_.tolerance_mb_ = tuning.TRGMBE;
310 modelParam_.tolerance_mb_relaxed_ = tuning.XXXMBE;
311 modelParam_.newton_max_iter_ = tuning.NEWTMX;
312 modelParam_.newton_min_iter_ = tuning.NEWTMN;
313 if (terminalOutput_) {
314 const auto msg = fmt::format(fmt::runtime(
"Tuning values: "
315 "MB: {:.2e}, CNV: {:.2e}, NEWTMN: {}, NEWTMX: {}"),
316 tuning.TRGMBE, tuning.TRGCNV, tuning.NEWTMN, tuning.NEWTMX);
318 if (tuning.TRGTTE_has_value) {
319 OpmLog::warning(
"Tuning item 2-1 (TRGTTE) is not supported.");
321 if (tuning.TRGLCV_has_value) {
322 OpmLog::warning(
"Tuning item 2-4 (TRGLCV) is not supported.");
324 if (tuning.XXXTTE_has_value) {
325 OpmLog::warning(
"Tuning item 2-5 (XXXTTE) is not supported.");
327 if (tuning.XXXLCV_has_value) {
328 OpmLog::warning(
"Tuning item 2-8 (XXXLCV) is not supported.");
330 if (tuning.XXXWFL_has_value) {
331 OpmLog::warning(
"Tuning item 2-9 (XXXWFL) is not supported.");
333 if (tuning.TRGFIP_has_value) {
334 OpmLog::warning(
"Tuning item 2-10 (TRGFIP) is not supported.");
336 if (tuning.TRGSFT_has_value) {
337 OpmLog::warning(
"Tuning item 2-11 (TRGSFT) is not supported.");
339 if (tuning.THIONX_has_value) {
340 OpmLog::warning(
"Tuning item 2-12 (THIONX) is not supported.");
342 if (tuning.TRWGHT_has_value) {
343 OpmLog::warning(
"Tuning item 2-13 (TRWGHT) is not supported.");
345 if (tuning.LITMAX_has_value) {
346 OpmLog::warning(
"Tuning item 3-3 (LITMAX) is not supported.");
348 if (tuning.LITMIN_has_value) {
349 OpmLog::warning(
"Tuning item 3-4 (LITMIN) is not supported.");
351 if (tuning.MXWSIT_has_value) {
352 OpmLog::warning(
"Tuning item 3-5 (MXWSIT) is not supported.");
354 if (tuning.MXWPIT_has_value) {
355 OpmLog::warning(
"Tuning item 3-6 (MXWPIT) is not supported.");
357 if (tuning.DDPLIM_has_value) {
358 OpmLog::warning(
"Tuning item 3-7 (DDPLIM) is not supported.");
360 if (tuning.DDSLIM_has_value) {
361 OpmLog::warning(
"Tuning item 3-8 (DDSLIM) is not supported.");
363 if (tuning.TRGDPR_has_value) {
364 OpmLog::warning(
"Tuning item 3-9 (TRGDPR) is not supported.");
366 if (tuning.XXXDPR_has_value) {
367 OpmLog::warning(
"Tuning item 3-10 (XXXDPR) is not supported.");
369 if (tuning.MNWRFP_has_value) {
370 OpmLog::warning(
"Tuning item 3-11 (MNWRFP) is not supported.");
375 void updateTUNINGDP(
const TuningDp& tuning_dp)
378 modelParam_.tolerance_max_dp_ = tuning_dp.TRGDDP;
379 modelParam_.tolerance_max_ds_ = tuning_dp.TRGDDS;
380 modelParam_.tolerance_max_drs_ = tuning_dp.TRGDDRS;
381 modelParam_.tolerance_max_drv_ = tuning_dp.TRGDDRV;
384 if (terminalOutput_) {
386 if (tuning_dp.TRGLCV_has_value) {
387 OpmLog::warning(
"TUNINGDP item 1 (TRGLCV) is not supported.");
389 if (tuning_dp.XXXLCV_has_value) {
390 OpmLog::warning(
"TUNINGDP item 2 (XXXLCV) is not supported.");
395 bool runStep(SimulatorTimer& timer)
397 if (schedule().exitStatus().has_value()) {
398 if (terminalOutput_) {
399 OpmLog::info(
"Stopping simulation since EXIT was triggered by an action keyword.");
401 report_.success.exit_status = schedule().exitStatus().value();
405 if (serializer_.shouldLoad()) {
406 serializer_.loadTimerInfo(timer);
410 if (terminalOutput_) {
411 std::ostringstream ss;
413 OpmLog::debug(ss.str());
414 details::outputReportStep(timer);
418 if (timer.initialStep()) {
419 Dune::Timer perfTimer;
422 simulator_.setEpisodeIndex(-1);
423 simulator_.setEpisodeLength(0.0);
424 simulator_.setTimeStepSize(0.0);
425 wellModel_().beginReportStep(timer.currentStepNum());
426 simulator_.problem().writeOutput(
true);
428 report_.success.output_write_time += perfTimer.stop();
432 solverTimer_->start();
435 solver_ = createSolver(wellModel_());
438 simulator_.startNextEpisode(
439 simulator_.startTime()
440 + schedule().seconds(timer.currentStepNum()),
441 timer.currentStepLength());
442 simulator_.setEpisodeIndex(timer.currentStepNum());
444 if (serializer_.shouldLoad()) {
445 wellModel_().prepareDeserialize(serializer_.loadStep() - 1);
446 serializer_.loadState();
447 simulator_.model().invalidateAndUpdateIntensiveQuantities(0);
450 this->solver_->model().beginReportStep();
452 const bool enableTUNING = Parameters::Get<Parameters::EnableTuning>();
459 if (adaptiveTimeStepping_) {
460 auto tuningUpdater = [enableTUNING,
this,
461 reportStep = timer.currentStepNum()](
const double curr_time,
462 double dt,
const int timeStep)
464 auto& schedule = this->simulator_.vanguard().schedule();
465 auto& events = this->schedule()[reportStep].events();
468 if (events.hasEvent(ScheduleEvents::TUNING_CHANGE)) {
470 schedule.clear_event(ScheduleEvents::TUNING_CHANGE, reportStep);
471 const auto& sched_state = schedule[reportStep];
472 const auto& max_next_tstep = sched_state.max_next_tstep(enableTUNING);
473 const auto& tuning = sched_state.tuning();
476 adaptiveTimeStepping_->updateTUNING(max_next_tstep, tuning);
479 solver_->model().updateTUNING(tuning);
480 this->updateTUNING(tuning);
481 dt = this->adaptiveTimeStepping_->suggestedNextStep();
484 this->adaptiveTimeStepping_->updateNEXTSTEP(max_next_tstep);
486 result = max_next_tstep > 0;
489 if (events.hasEvent(ScheduleEvents::TUNINGDP_CHANGE)) {
491 schedule.clear_event(ScheduleEvents::TUNINGDP_CHANGE, reportStep);
496 const auto& sched_state = schedule[reportStep];
497 const auto& tuning_dp = sched_state.tuning_dp();
498 solver_->model().updateTUNINGDP(tuning_dp);
499 this->updateTUNINGDP(tuning_dp);
502 const auto& wcycle = schedule[reportStep].wcycle.get();
503 if (wcycle.empty()) {
507 const auto& wmatcher = schedule.wellMatcher(reportStep);
508 double wcycle_time_step =
509 wcycle.nextTimeStep(curr_time,
512 this->wellModel_().wellOpenTimes(),
513 this->wellModel_().wellCloseTimes(),
515 &wg_events = this->wellModel_().reportStepStartEvents()]
516 (
const std::string& name)
521 return wg_events.hasEvent(name, ScheduleEvents::REQUEST_OPEN_WELL);
524 wcycle_time_step = this->grid().comm().min(wcycle_time_step);
525 if (dt != wcycle_time_step) {
526 this->adaptiveTimeStepping_->updateNEXTSTEP(wcycle_time_step);
533 tuningUpdater(timer.simulationTimeElapsed(),
534 this->adaptiveTimeStepping_->suggestedNextStep(), 0);
536#ifdef RESERVOIR_COUPLING_ENABLED
537 if (this->reservoirCouplingMaster_) {
538 this->reservoirCouplingMaster_->maybeSpawnSlaveProcesses(timer.currentStepNum());
539 this->reservoirCouplingMaster_->maybeActivate(timer.currentStepNum());
541 else if (this->reservoirCouplingSlave_) {
542 this->reservoirCouplingSlave_->maybeActivate(timer.currentStepNum());
545 const auto& events = schedule()[timer.currentStepNum()].events();
546 bool event = events.hasEvent(ScheduleEvents::NEW_WELL) ||
547 events.hasEvent(ScheduleEvents::INJECTION_TYPE_CHANGED) ||
548 events.hasEvent(ScheduleEvents::WELL_SWITCHED_INJECTOR_PRODUCER) ||
549 events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE) ||
550 events.hasEvent(ScheduleEvents::INJECTION_UPDATE) ||
551 events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
552 auto stepReport = adaptiveTimeStepping_->step(timer, *solver_, event, tuningUpdater);
553 report_ += stepReport;
556 auto stepReport = solver_->step(timer,
nullptr);
557 report_ += stepReport;
559 simulator_.problem().setSubStepReport(stepReport);
560 simulator_.problem().setSimulationReport(report_);
561 simulator_.problem().endTimeStep();
562 if (terminalOutput_) {
563 std::ostringstream ss;
564 stepReport.reportStep(ss);
565 OpmLog::info(ss.str());
570 Dune::Timer perfTimer;
572 const double nextstep = adaptiveTimeStepping_ ? adaptiveTimeStepping_->suggestedNextStep() : -1.0;
573 simulator_.problem().setNextTimeStepSize(nextstep);
574 simulator_.problem().writeOutput(
true);
575 report_.success.output_write_time += perfTimer.stop();
577 solver_->model().endReportStep();
580 solverTimer_->stop();
583 report_.success.solver_time += solverTimer_->secsSinceStart();
585 if (this->grid().comm().rank() == 0) {
588 const auto& reps = this->solver_->model().stepReports();
589 convergence_output_.write(reps);
595 if (terminalOutput_) {
597 "Time step took " + std::to_string(solverTimer_->secsSinceStart()) +
" seconds; "
598 "total solver time " + std::to_string(report_.success.solver_time) +
" seconds.";
602 serializer_.save(timer);
607 SimulatorReport finalize()
611 Dune::Timer finalOutputTimer;
612 finalOutputTimer.start();
614 simulator_.problem().finalizeOutput();
615 report_.success.output_write_time += finalOutputTimer.stop();
620 report_.success.total_time = totalTimer_->secsSinceStart();
621 report_.success.converged =
true;
626 const Grid& grid()
const
627 {
return simulator_.vanguard().grid(); }
629 template<
class Serializer>
630 void serializeOp(Serializer& serializer)
632 serializer(simulator_);
634 serializer(adaptiveTimeStepping_);
637 const Model& model()
const
638 {
return solver_->model(); }
643 [[maybe_unused]]
const std::string& groupName)
override
646 serializer.read(*
this, groupName,
"simulator_data");
652 [[maybe_unused]]
const std::string& groupName)
const override
655 serializer.write(*
this, groupName,
"simulator_data");
662 std::ostringstream str;
667 simulator_.vanguard().caseName(),
674 return simulator_.vanguard().globalCell();
677 std::unique_ptr<Solver> createSolver(WellModel& wellModel)
679 auto model = std::make_unique<Model>(simulator_,
684 if (this->modelParam_.write_partitions_) {
685 const auto& iocfg = this->eclState().cfg().io();
687 const auto odir = iocfg.getOutputDir()
688 / std::filesystem::path {
"partition" }
689 / iocfg.getBaseName();
691 if (this->grid().comm().rank() == 0) {
692 create_directories(odir);
695 this->grid().comm().barrier();
697 model->writePartitions(odir);
699 this->modelParam_.write_partitions_ =
false;
702 return std::make_unique<Solver>(solverParam_, std::move(model));
705 const EclipseState& eclState()
const
706 {
return simulator_.vanguard().eclState(); }
709 const Schedule& schedule()
const
710 {
return simulator_.vanguard().schedule(); }
712 bool isRestart()
const
714 const auto& initconfig = eclState().getInitConfig();
715 return initconfig.restartRequested();
718 WellModel& wellModel_()
719 {
return simulator_.problem().wellModel(); }
721 const WellModel& wellModel_()
const
722 {
return simulator_.problem().wellModel(); }
725 Simulator& simulator_;
727 ModelParameters modelParam_;
728 SolverParameters solverParam_;
730 std::unique_ptr<Solver> solver_;
734 bool terminalOutput_;
736 SimulatorReport report_;
737 std::unique_ptr<time::StopWatch> solverTimer_;
738 std::unique_ptr<time::StopWatch> totalTimer_;
739 std::unique_ptr<TimeStepper> adaptiveTimeStepping_;
741 SimulatorConvergenceOutput convergence_output_{};
743#ifdef RESERVOIR_COUPLING_ENABLED
744 bool slaveMode_{
false};
745 std::unique_ptr<ReservoirCouplingMaster<Scalar>> reservoirCouplingMaster_{
nullptr};
746 std::unique_ptr<ReservoirCouplingSlave<Scalar>> reservoirCouplingSlave_{
nullptr};
749 SimulatorSerializer serializer_;