SimulatorFullyImplicit_impl.hpp

/*
  Copyright 2013, 2015, 2020 SINTEF Digital, Mathematics and Cybernetics.
  Copyright 2015 Andreas Lauser
  Copyright 2017 IRIS

  This file is part of the Open Porous Media project (OPM).

  OPM is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  OPM is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef OPM_SIMULATOR_FULLY_IMPLICIT_IMPL_HEADER_INCLUDED
#define OPM_SIMULATOR_FULLY_IMPLICIT_IMPL_HEADER_INCLUDED

// Improve IDE experience
#ifndef OPM_SIMULATOR_FULLY_IMPLICIT_HEADER_INCLUDED
#include <config.h>
#include <opm/simulators/flow/SimulatorFullyImplicit.hpp>
#endif

#include <opm/input/eclipse/Units/UnitSystem.hpp>

#include <opm/models/tpsa/tpsanewtonmethodparams.hpp>

#include <opm/simulators/linalg/TPSALinearSolverParameters.hpp>

#include <fmt/format.h>

#include <filesystem>
#include <sstream>

namespace Opm {

template<class TypeTag>
SimulatorFullyImplicit<TypeTag>::
SimulatorFullyImplicit(Simulator& simulator)
    : simulator_(simulator)
    , serializer_(*this,
                  FlowGenericVanguard::comm(),
                  simulator_.vanguard().eclState().getIOConfig(),
                  Parameters::Get<Parameters::SaveStep>(),
                  Parameters::Get<Parameters::LoadStep>(),
                  Parameters::Get<Parameters::SaveFile>(),
                  Parameters::Get<Parameters::LoadFile>())
{
    // Only rank 0 does print to std::cout, and only if specifically requested.
    this->terminalOutput_ = false;
    if (this->grid().comm().rank() == 0) {
        this->terminalOutput_ = Parameters::Get<Parameters::EnableTerminalOutput>();

        auto getPhaseName = ConvergenceOutputThread::ComponentToPhaseName {
            [compNames = typename Model::ComponentName{}](const int compIdx)
            { return std::string_view { compNames.name(compIdx) }; }
        };

        if (!simulator_.vanguard().eclState().getIOConfig().initOnly()) {
            this->convergence_output_.
                startThread(this->simulator_.vanguard().eclState(),
                            Parameters::Get<Parameters::OutputExtraConvergenceInfo>(),
                            R"(OutputExtraConvergenceInfo (--output-extra-convergence-info))",
                            getPhaseName);
        }
    }
}

template<class TypeTag>
SimulatorFullyImplicit<TypeTag>::
~SimulatorFullyImplicit()
{
    // Safe to call on all ranks, not just the I/O rank.
    convergence_output_.endThread();
}

template<class TypeTag>
void
SimulatorFullyImplicit<TypeTag>::
registerParameters()
{
    ModelParameters::registerParameters();
    SolverParameters::registerParameters();
    TimeStepper::registerParameters();
    detail::registerSimulatorParameters();

    TpsaNewtonMethodParams<Scalar>::registerParameters();
    TpsaLinearSolverParameters::registerParameters();
}

#ifdef RESERVOIR_COUPLING_ENABLED
template<class TypeTag>
SimulatorReport
SimulatorFullyImplicit<TypeTag>::
run(SimulatorTimer& timer, int argc, char** argv)
{
    init(timer, argc, argv);
#else
template<class TypeTag>
SimulatorReport
SimulatorFullyImplicit<TypeTag>::
run(SimulatorTimer& timer)
{
    init(timer);
#endif
    // Make cache up to date. No need for updating it in elementCtx.
    // NB! Need to be at the correct step in case of restart
    simulator_.setEpisodeIndex(timer.currentStepNum());
    simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
    // Main simulation loop.
    while (!timer.done()) {
        simulator_.problem().writeReports(timer);
        bool continue_looping = runStep(timer);
        if (!continue_looping) break;
    }
    simulator_.problem().writeReports(timer);

#ifdef RESERVOIR_COUPLING_ENABLED
    // Clean up MPI intercommunicators before MPI_Finalize()
    // Master sends terminate=1 signal; slave receives it and both call MPI_Comm_disconnect()
    if (this->reservoirCouplingMaster_) {
        this->reservoirCouplingMaster_->sendTerminateAndDisconnect();
    }
    else if (this->reservoirCouplingSlave_ && !this->reservoirCouplingSlave_->terminated()) {
        // TODO: Implement GECON item 8: stop master process when a slave finishes
        // Only call if not already terminated via maybeReceiveTerminateSignalFromMaster()
        // (which happens when master finishes before slave reaches end of its loop)
        this->reservoirCouplingSlave_->receiveTerminateAndDisconnect();
    }
#endif

    return finalize();
}

#ifdef RESERVOIR_COUPLING_ENABLED
template<class TypeTag>
bool
SimulatorFullyImplicit<TypeTag>::
checkRunningAsReservoirCouplingMaster()
{
    for (std::size_t report_step = 0; report_step < this->schedule().size(); ++report_step) {
        auto rescoup = this->schedule()[report_step].rescoup();
        auto slave_count = rescoup.slaveCount();
        auto master_group_count = rescoup.masterGroupCount();
        // Master mode is enabled when SLAVES keyword is present.
        // - Prediction mode: SLAVES + GRUPMAST (master allocates rates)
        // - History mode: SLAVES only (master synchronizes time-stepping)
        if (slave_count > 0) {
            return true;
        }
        else if (master_group_count > 0) {
            // GRUPMAST without SLAVES is invalid
            throw ReservoirCouplingError(
                "Inconsistent reservoir coupling master schedule: "
                "Master group count is greater than 0 but slave count is 0"
            );
        }
    }
    return false;
}
#endif

#ifdef RESERVOIR_COUPLING_ENABLED
template<class TypeTag>
void
SimulatorFullyImplicit<TypeTag>::
init(const SimulatorTimer& timer, int argc, char** argv)
{
    auto slave_mode = Parameters::Get<Parameters::Slave>();
    if (slave_mode) {
        this->reservoirCouplingSlave_ =
            std::make_unique<ReservoirCouplingSlave<Scalar>>(
                FlowGenericVanguard::comm(),
                this->schedule(), timer
            );
        this->reservoirCouplingSlave_->sendAndReceiveInitialData();
        this->simulator_.setReservoirCouplingSlave(this->reservoirCouplingSlave_.get());
        wellModel_().setReservoirCouplingSlave(this->reservoirCouplingSlave_.get());
    }
    else {
        auto master_mode = checkRunningAsReservoirCouplingMaster();
        if (master_mode) {
            this->reservoirCouplingMaster_ =
                std::make_unique<ReservoirCouplingMaster<Scalar>>(
                    FlowGenericVanguard::comm(),
                    this->schedule(),
                    argc, argv
                );
            this->simulator_.setReservoirCouplingMaster(this->reservoirCouplingMaster_.get());
            wellModel_().setReservoirCouplingMaster(this->reservoirCouplingMaster_.get());
        }
    }
#else
template<class TypeTag>
void
SimulatorFullyImplicit<TypeTag>::
init(const SimulatorTimer& timer)
{
#endif
    simulator_.setEpisodeIndex(-1);

    // Create timers and file for writing timing info.
    solverTimer_ = std::make_unique<time::StopWatch>();
    totalTimer_ = std::make_unique<time::StopWatch>();
    totalTimer_->start();

    // adaptive time stepping
    bool enableAdaptive = Parameters::Get<Parameters::EnableAdaptiveTimeStepping>();
    bool enableTUNING = Parameters::Get<Parameters::EnableTuning>();
    if (enableAdaptive) {
        const UnitSystem& unitSystem = this->simulator_.vanguard().eclState().getUnits();
        const auto& sched_state = schedule()[timer.currentStepNum()];
        auto max_next_tstep = sched_state.max_next_tstep(enableTUNING);
        if (enableTUNING) {
            adaptiveTimeStepping_ = std::make_unique<TimeStepper>(max_next_tstep,
                                                                  sched_state.tuning(),
                                                                  unitSystem, report_, terminalOutput_);
        }
        else {
            adaptiveTimeStepping_ = std::make_unique<TimeStepper>(unitSystem, report_, max_next_tstep, terminalOutput_);
        }
        if (isRestart()) {
            // For restarts the simulator may have gotten some information
            // about the next timestep size from the OPMEXTRA field
            adaptiveTimeStepping_->setSuggestedNextStep(simulator_.timeStepSize());
        }
    }
}

template<class TypeTag>
void
SimulatorFullyImplicit<TypeTag>::
updateTUNING(const Tuning& tuning)
{
    modelParam_.tolerance_cnv_ = tuning.TRGCNV;
    modelParam_.tolerance_cnv_relaxed_ = tuning.XXXCNV;
    modelParam_.tolerance_mb_ = tuning.TRGMBE;
    modelParam_.tolerance_mb_relaxed_ = tuning.XXXMBE;
    modelParam_.newton_max_iter_ = tuning.NEWTMX;
    modelParam_.newton_min_iter_ = tuning.NEWTMN;
    if (terminalOutput_) {
        detail::logTuning(tuning);
    }
}

template<class TypeTag>
void
SimulatorFullyImplicit<TypeTag>::
updateTUNINGDP(const TuningDp& tuning_dp)
{
    // NOTE: If TUNINGDP item is _not_ set it should be 0.0
    modelParam_.tolerance_max_dp_ = tuning_dp.TRGDDP;
    modelParam_.tolerance_max_ds_ = tuning_dp.TRGDDS;
    modelParam_.tolerance_max_drs_ = tuning_dp.TRGDDRS;
    modelParam_.tolerance_max_drv_ = tuning_dp.TRGDDRV;

    // Terminal warnings
    if (terminalOutput_) {
        // Warnings unsupported items
        if (tuning_dp.TRGLCV_has_value) {
            OpmLog::warning("TUNINGDP item 1 (TRGLCV) is not supported.");
        }
        if (tuning_dp.XXXLCV_has_value) {
            OpmLog::warning("TUNINGDP item 2 (XXXLCV) is not supported.");
        }
    }
}

template<class TypeTag>
bool
SimulatorFullyImplicit<TypeTag>::
runStep(SimulatorTimer& timer)
{
    if (schedule().exitStatus().has_value()) {
        if (terminalOutput_) {
            OpmLog::info("Stopping simulation since EXIT was triggered by an action keyword.");
        }
        report_.success.exit_status = schedule().exitStatus().value();
        return false;
    }

    if (serializer_.shouldLoad()) {
        serializer_.loadTimerInfo(timer);
    }

    // Report timestep.
    if (terminalOutput_) {
        std::ostringstream ss;
        timer.report(ss);
        OpmLog::debug(ss.str());
        details::outputReportStep(timer);
    }

    // write the inital state at the report stage
    if (timer.initialStep()) {
        Dune::Timer perfTimer;
        perfTimer.start();

        simulator_.setEpisodeIndex(-1);
        simulator_.setEpisodeLength(0.0);
        simulator_.setTimeStepSize(0.0);
        wellModel_().beginReportStep(timer.currentStepNum());
        simulator_.problem().writeOutput(true);

        report_.success.output_write_time += perfTimer.stop();
    }

    // Run a multiple steps of the solver depending on the time step control.
    solverTimer_->start();

    if (!solver_) {
        solver_ = createSolver(wellModel_());
    }

    simulator_.startNextEpisode(
        simulator_.startTime()
           + schedule().seconds(timer.currentStepNum()),
        timer.currentStepLength());
    simulator_.setEpisodeIndex(timer.currentStepNum());

    if (serializer_.shouldLoad()) {
        wellModel_().prepareDeserialize(serializer_.loadStep() - 1);
        serializer_.loadState();
        simulator_.model().invalidateAndUpdateIntensiveQuantities(/*timeIdx=*/0);
    }

    this->solver_->model().beginReportStep();

    const bool enableTUNING = Parameters::Get<Parameters::EnableTuning>();

    // If sub stepping is enabled allow the solver to sub cycle
    // in case the report steps are too large for the solver to converge
    //
    // \Note: The report steps are met in any case
    // \Note: The sub stepping will require a copy of the state variables
    if (adaptiveTimeStepping_) {
        auto tuningUpdater = [enableTUNING, this,
                              reportStep = timer.currentStepNum()](const double curr_time,
                                                                   double substep_length,
                                                                   const int sub_step_number)
        {
            auto& schedule = this->simulator_.vanguard().schedule();
            auto& events = this->schedule()[reportStep].events();

            bool result = false;
            if (events.hasEvent(ScheduleEvents::TUNING_CHANGE)) {
                // Unset the event to not trigger it again on the next sub step
                schedule.clear_event(ScheduleEvents::TUNING_CHANGE, reportStep);
                const auto& sched_state = schedule[reportStep];
                const auto& max_next_tstep = sched_state.max_next_tstep(enableTUNING);
                const auto& tuning = sched_state.tuning();

                if (enableTUNING) {
                    adaptiveTimeStepping_->updateTUNING(max_next_tstep, tuning);
                    // \Note: Assumes TUNING is only used with adaptive time-stepping
                    // \Note: Need to update both solver (model) and simulator since solver is re-created each report step.
                    solver_->model().updateTUNING(tuning);
                    this->updateTUNING(tuning);
                    substep_length = this->adaptiveTimeStepping_->suggestedNextStep();
                } else {
                    substep_length = max_next_tstep;
                    this->adaptiveTimeStepping_->updateNEXTSTEP(max_next_tstep);
                }
                result = max_next_tstep > 0;
            }

            if (events.hasEvent(ScheduleEvents::TUNINGDP_CHANGE)) {
                // Unset the event to not trigger it again on the next sub step
                schedule.clear_event(ScheduleEvents::TUNINGDP_CHANGE, reportStep);

                // Update TUNINGDP parameters
                // NOTE: Need to update both solver (model) and simulator since solver is re-created each report
                // step.
                const auto& sched_state = schedule[reportStep];
                const auto& tuning_dp = sched_state.tuning_dp();
                solver_->model().updateTUNINGDP(tuning_dp);
                this->updateTUNINGDP(tuning_dp);
            }

            const auto& wcycle = schedule[reportStep].wcycle.get();
            if (wcycle.empty()) {
                return result;
            }

            const auto& wmatcher = schedule.wellMatcher(reportStep);
            double wcycle_time_step =
                wcycle.nextTimeStep(curr_time,
                                    substep_length,
                                    wmatcher,
                                    this->wellModel_().wellOpenTimes(),
                                    this->wellModel_().wellCloseTimes(),
                                    [sub_step_number,
                                     &wg_events = this->wellModel_().reportStepStartEvents()]
                                    (const std::string& name)
                                    {
                                        if (sub_step_number != 0) {
                                            return false;
                                        }
                                        return wg_events.hasEvent(name, ScheduleEvents::REQUEST_OPEN_WELL);
                                    });

            wcycle_time_step = this->grid().comm().min(wcycle_time_step);
            if (substep_length != wcycle_time_step) {
                this->adaptiveTimeStepping_->updateNEXTSTEP(wcycle_time_step);
                return true;
            }

            return result;
        };

        tuningUpdater(timer.simulationTimeElapsed(),
                      this->adaptiveTimeStepping_->suggestedNextStep(), 0);

#ifdef RESERVOIR_COUPLING_ENABLED
        if (this->reservoirCouplingMaster_) {
            this->reservoirCouplingMaster_->maybeSpawnSlaveProcesses(timer.currentStepNum());
            this->reservoirCouplingMaster_->maybeActivate(timer.currentStepNum());
        }
        else if (this->reservoirCouplingSlave_) {
            this->reservoirCouplingSlave_->maybeActivate(timer.currentStepNum());
        }
#endif
        const auto& events = schedule()[timer.currentStepNum()].events();
        bool event = events.hasEvent(ScheduleEvents::NEW_WELL) ||
            events.hasEvent(ScheduleEvents::INJECTION_TYPE_CHANGED) ||
            events.hasEvent(ScheduleEvents::WELL_SWITCHED_INJECTOR_PRODUCER) ||
            events.hasEvent(ScheduleEvents::PRODUCTION_UPDATE) ||
            events.hasEvent(ScheduleEvents::INJECTION_UPDATE) ||
            events.hasEvent(ScheduleEvents::WELL_STATUS_CHANGE);
        auto stepReport = adaptiveTimeStepping_->step(timer, *solver_, event, tuningUpdater);
        report_ += stepReport;
#ifdef RESERVOIR_COUPLING_ENABLED
        // If the master ended its schedule first (e.g. an END keyword truncates the master
        // SCHEDULE), the slave received the terminate signal and disconnected its
        // intercommunicator inside step() above. The coupled run is over for this slave:
        // finish the current report step cleanly and stop the run loop instead of advancing
        // to another report step, which would issue an MPI_Recv on the now-disconnected
        // communicator and abort the job.
        if (this->reservoirCouplingSlave_ && this->reservoirCouplingSlave_->terminated()) {
            this->handleSlaveTerminated_();
            return false;   // breaks the while(!timer.done()) loop in run()
        }
#endif
    } else {
        // solve for complete report step
        auto stepReport = solver_->step(timer, nullptr);
        report_ += stepReport;
        // Pass simulation report to eclwriter for summary output
        simulator_.problem().setSubStepReport(stepReport);
        simulator_.problem().setSimulationReport(report_);
        simulator_.problem().endTimeStep();
        if (terminalOutput_) {
            std::ostringstream ss;
            stepReport.reportStep(ss);
            OpmLog::info(ss.str());
        }
    }

    // write simulation state at the report stage
    Dune::Timer perfTimer;
    perfTimer.start();
    const double nextstep = adaptiveTimeStepping_ ? adaptiveTimeStepping_->suggestedNextStep() : -1.0;
    simulator_.problem().setNextTimeStepSize(nextstep);
    simulator_.problem().writeOutput(true);
    report_.success.output_write_time += perfTimer.stop();

    solver_->model().endReportStep();

    // take time that was used to solve system for this reportStep
    solverTimer_->stop();

    // update timing.
    report_.success.solver_time += solverTimer_->secsSinceStart();

    if (this->grid().comm().rank() == 0) {
        // Grab the step convergence reports that are new since last we
        // were here.
        const auto& reps = this->solver_->model().stepReports();
        convergence_output_.write(reps);
    }

    // Increment timer, remember well state.
    ++timer;

    if (terminalOutput_) {
        std::string msg =
            "Time step took " + std::to_string(solverTimer_->secsSinceStart()) + " seconds; "
            "total solver time " + std::to_string(report_.success.solver_time) + " seconds.";
        OpmLog::debug(msg);
    }

    serializer_.save(timer);

    return true;
}

#ifdef RESERVOIR_COUPLING_ENABLED
template<class TypeTag>
void
SimulatorFullyImplicit<TypeTag>::
handleSlaveTerminated_()
{
    if (terminalOutput_) {
        OpmLog::info("Reservoir coupling: master simulation has ended; "
                     "stopping slave simulation gracefully.");
    }
    // The slave stops as soon as it is terminated, so no further substeps run for this report
    // step. Per-substep SUMMARY output for any substeps already completed in this report step
    // was written by the adaptive substep loop; the report-step-level restart write is
    // intentionally skipped because a terminated report step is partial - the last completed
    // report step is the clean restart point. finalize() flushes the remaining output.
    // Balance the beginReportStep() issued earlier in runStep().
    this->solver_->model().endReportStep();
    // Account for this report step's solver time and leave the timer stopped: runStep()'s own
    // solverTimer_->stop()/accumulation below is bypassed by the early return on termination.
    this->solverTimer_->stop();
    this->report_.success.solver_time += this->solverTimer_->secsSinceStart();
}
#endif

template<class TypeTag>
SimulatorReport
SimulatorFullyImplicit<TypeTag>::
finalize()
{
    // make sure all output is written to disk before run is finished
    {
        Dune::Timer finalOutputTimer;
        finalOutputTimer.start();

        simulator_.problem().finalizeOutput();
        report_.success.output_write_time += finalOutputTimer.stop();
    }

    // Stop timer and create timing report
    totalTimer_->stop();
    report_.success.total_time = totalTimer_->secsSinceStart();
    report_.success.converged = true;

    return report_;
}

template<class TypeTag>
template<class Serializer>
void
SimulatorFullyImplicit<TypeTag>::
serializeOp(Serializer& serializer)
{
    serializer(simulator_);
    serializer(report_);
    serializer(adaptiveTimeStepping_);
}

template<class TypeTag>
void
SimulatorFullyImplicit<TypeTag>::
loadState([[maybe_unused]] HDF5Serializer& serializer,
          [[maybe_unused]] const std::string& groupName)
{
#if HAVE_HDF5
    serializer.read(*this, groupName, "simulator_data");
#endif
}

template<class TypeTag>
void
SimulatorFullyImplicit<TypeTag>::
saveState([[maybe_unused]] HDF5Serializer& serializer,
          [[maybe_unused]] const std::string& groupName) const
{
#if HAVE_HDF5
    serializer.write(*this, groupName, "simulator_data");
#endif
}

template<class TypeTag>
std::array<std::string,5>
SimulatorFullyImplicit<TypeTag>::
getHeader() const
{
    std::ostringstream str;
    Parameters::printValues(str);
    return {"OPM Flow",
            moduleVersion(),
            compileTimestamp(),
            simulator_.vanguard().caseName(),
            str.str()};
}

template<class TypeTag>
std::unique_ptr<typename SimulatorFullyImplicit<TypeTag>::Solver>
SimulatorFullyImplicit<TypeTag>::
createSolver(WellModel& wellModel)
{
    auto model = std::make_unique<Model>(simulator_,
                                         modelParam_,
                                         wellModel,
                                         terminalOutput_);

    if (this->modelParam_.write_partitions_) {
        const auto& iocfg = this->eclState().cfg().io();

        const auto odir = iocfg.getOutputDir()
            / std::filesystem::path { "partition" }
            / iocfg.getBaseName();

        if (this->grid().comm().rank() == 0) {
            create_directories(odir);
        }

        this->grid().comm().barrier();

        model->writePartitions(odir);

        this->modelParam_.write_partitions_ = false;
    }

    return std::make_unique<Solver>(solverParam_, std::move(model));
}

} // namespace Opm

#endif // OPM_SIMULATOR_FULLY_IMPLICIT_IMPL_HEADER_INCLUDED