AdaptiveTimeStepping.hpp

Go to the documentation of this file.

/*
*/
#ifndef OPM_ADAPTIVE_TIME_STEPPING_HPP
#define OPM_ADAPTIVE_TIME_STEPPING_HPP
 
#include <dune/common/version.hh>
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 8)
#include <dune/istl/istlexception.hh>
#else
#include <dune/istl/ilu.hh>
#endif
 
#include <opm/common/Exceptions.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/OpmLog/OpmLog.hpp>
 
#include <opm/core/props/phaseUsageFromDeck.hpp>
 
#include <opm/grid/utility/StopWatch.hpp>
 
#include <opm/input/eclipse/Units/Units.hpp>
#include <opm/input/eclipse/Units/UnitSystem.hpp>
 
#include <opm/input/eclipse/Schedule/Tuning.hpp>
 
#include <opm/models/utils/basicproperties.hh>
#include <opm/models/utils/parametersystem.hh>
#include <opm/models/utils/propertysystem.hh>
 
#include <opm/simulators/timestepping/AdaptiveSimulatorTimer.hpp>
#include <opm/simulators/timestepping/EclTimeSteppingParams.hpp>
#include <opm/simulators/timestepping/SimulatorReport.hpp>
#include <opm/simulators/timestepping/SimulatorTimer.hpp>
#include <opm/simulators/timestepping/TimeStepControl.hpp>
#include <opm/simulators/timestepping/TimeStepControlInterface.hpp>
 
#include <fmt/format.h>
 
#include <algorithm>
#include <cassert>
#include <cmath>
#include <functional>
#include <memory>
#include <set>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
 
namespace Opm::Properties {
 
namespace TTag {
struct FlowTimeSteppingParameters {
  using InheritsFrom = std::tuple<EclTimeSteppingParameters>;
};
}
 
template<class TypeTag, class MyTypeTag>
struct SolverContinueOnConvergenceFailure {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct SolverMaxRestarts {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct SolverVerbosity {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct TimeStepVerbosity {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct InitialTimeStepInDays {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct FullTimeStepInitially {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct TimeStepControl {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct TimeStepControlTolerance {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct TimeStepControlTargetIterations {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct TimeStepControlTargetNewtonIterations {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct TimeStepControlDecayRate {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct TimeStepControlGrowthRate {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct TimeStepControlDecayDampingFactor {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct TimeStepControlGrowthDampingFactor {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct TimeStepControlFileName {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct MinTimeStepBeforeShuttingProblematicWellsInDays {
    using type = UndefinedProperty;
};
template<class TypeTag, class MyTypeTag>
struct MinTimeStepBasedOnNewtonIterations {
    using type = UndefinedProperty;
};
 
template<class TypeTag>
struct SolverContinueOnConvergenceFailure<TypeTag, TTag::FlowTimeSteppingParameters> {
    static constexpr bool value = false;
};
template<class TypeTag>
struct SolverMaxRestarts<TypeTag, TTag::FlowTimeSteppingParameters> {
    static constexpr int value = 10;
};
template<class TypeTag>
struct SolverVerbosity<TypeTag, TTag::FlowTimeSteppingParameters> {
    static constexpr int value = 1;
};
template<class TypeTag>
struct TimeStepVerbosity<TypeTag, TTag::FlowTimeSteppingParameters> {
    static constexpr int value = 1;
};
template<class TypeTag>
struct InitialTimeStepInDays<TypeTag, TTag::FlowTimeSteppingParameters> {
    using type = GetPropType<TypeTag, Scalar>;
    static constexpr type value = 1.0;
};
template<class TypeTag>
struct FullTimeStepInitially<TypeTag, TTag::FlowTimeSteppingParameters> {
    static constexpr bool value = false;
};
template<class TypeTag>
struct TimeStepControl<TypeTag, TTag::FlowTimeSteppingParameters> {
    static constexpr auto value = "pid+newtoniteration";
};
template<class TypeTag>
struct TimeStepControlTolerance<TypeTag, TTag::FlowTimeSteppingParameters> {
    using type = GetPropType<TypeTag, Scalar>;
    static constexpr type value = 1e-1;
};
template<class TypeTag>
struct TimeStepControlTargetIterations<TypeTag, TTag::FlowTimeSteppingParameters> {
    static constexpr int value = 30;
};
template<class TypeTag>
struct TimeStepControlTargetNewtonIterations<TypeTag, TTag::FlowTimeSteppingParameters> {
    static constexpr int value = 8;
};
template<class TypeTag>
struct TimeStepControlDecayRate<TypeTag, TTag::FlowTimeSteppingParameters> {
    using type = GetPropType<TypeTag, Scalar>;
    static constexpr type value = 0.75;
};
template<class TypeTag>
struct TimeStepControlGrowthRate<TypeTag, TTag::FlowTimeSteppingParameters> {
    using type = GetPropType<TypeTag, Scalar>;
    static constexpr type value = 1.25;
};
template<class TypeTag>
struct TimeStepControlDecayDampingFactor<TypeTag, TTag::FlowTimeSteppingParameters> {
    using type = GetPropType<TypeTag, Scalar>;
    static constexpr type value = 1.0;
};
template<class TypeTag>
struct TimeStepControlGrowthDampingFactor<TypeTag, TTag::FlowTimeSteppingParameters> {
    using type = GetPropType<TypeTag, Scalar>;
    static constexpr type value = 3.2;
};
template<class TypeTag>
struct TimeStepControlFileName<TypeTag, TTag::FlowTimeSteppingParameters> {
    static constexpr auto value = "timesteps";
};
template<class TypeTag>
struct MinTimeStepBeforeShuttingProblematicWellsInDays<TypeTag, TTag::FlowTimeSteppingParameters> {
    using type = GetPropType<TypeTag, Scalar>;
    static constexpr type value = 0.01;
};
 
template<class TypeTag>
struct MinTimeStepBasedOnNewtonIterations<TypeTag, TTag::FlowTimeSteppingParameters> {
    using type = GetPropType<TypeTag, Scalar>;
    static constexpr type value = 0.0;
};
 
} // namespace Opm::Properties
 
namespace Opm {
 
struct StepReport;
 
namespace detail {
 
void logTimer(const AdaptiveSimulatorTimer& substepTimer);
 
std::set<std::string> consistentlyFailingWells(const std::vector<StepReport>& sr);
 
}
 
    // AdaptiveTimeStepping
    //---------------------
    template<class TypeTag>
    class AdaptiveTimeStepping
    {
        template <class Solver>
        class SolutionTimeErrorSolverWrapper : public RelativeChangeInterface
        {
            const Solver& solver_;
        public:
            SolutionTimeErrorSolverWrapper(const Solver& solver)
              : solver_(solver)
            {}
 
            double relativeChange() const
            { return solver_.model().relativeChange(); }
        };
 
        template<class E>
        void logException_(const E& exception, bool verbose)
        {
            if (verbose) {
                std::string message;
                message = "Caught Exception: ";
                message += exception.what();
                OpmLog::debug(message);
            }
        }
 
    public:
        AdaptiveTimeStepping() = default;
 
        AdaptiveTimeStepping(const UnitSystem& unitSystem,
                             const double max_next_tstep = -1.0,
                             const bool terminalOutput = true)
            : timeStepControl_()
            , restartFactor_(Parameters::get<TypeTag, Properties::SolverRestartFactor>()) // 0.33
            , growthFactor_(Parameters::get<TypeTag, Properties::SolverGrowthFactor>()) // 2.0
            , maxGrowth_(Parameters::get<TypeTag, Properties::SolverMaxGrowth>()) // 3.0
            , maxTimeStep_(Parameters::get<TypeTag, Properties::SolverMaxTimeStepInDays>() * 24 * 60 * 60) // 365.25
            , minTimeStep_(unitSystem.to_si(UnitSystem::measure::time, Parameters::get<TypeTag, Properties::SolverMinTimeStep>())) // 1e-12;
            , ignoreConvergenceFailure_(Parameters::get<TypeTag, Properties::SolverContinueOnConvergenceFailure>()) // false;
            , solverRestartMax_(Parameters::get<TypeTag, Properties::SolverMaxRestarts>()) // 10
            , solverVerbose_(Parameters::get<TypeTag, Properties::SolverVerbosity>() > 0 && terminalOutput) // 2
            , timestepVerbose_(Parameters::get<TypeTag, Properties::TimeStepVerbosity>() > 0 && terminalOutput) // 2
            , suggestedNextTimestep_((max_next_tstep <= 0 ? Parameters::get<TypeTag, Properties::InitialTimeStepInDays>() : max_next_tstep) * 24 * 60 * 60) // 1.0
            , fullTimestepInitially_(Parameters::get<TypeTag, Properties::FullTimeStepInitially>()) // false
            , timestepAfterEvent_(Parameters::get<TypeTag, Properties::TimeStepAfterEventInDays>() * 24 * 60 * 60) // 1e30
            , useNewtonIteration_(false)
            , minTimeStepBeforeShuttingProblematicWells_(Parameters::get<TypeTag, Properties::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
 
        {
            init_(unitSystem);
        }
 
 
 
        AdaptiveTimeStepping(double max_next_tstep,
                             const Tuning& tuning,
                             const UnitSystem& unitSystem,
                             const bool terminalOutput = true)
            : timeStepControl_()
            , restartFactor_(tuning.TSFCNV)
            , growthFactor_(tuning.TFDIFF)
            , maxGrowth_(tuning.TSFMAX)
            , maxTimeStep_(tuning.TSMAXZ) // 365.25
            , minTimeStep_(tuning.TSFMIN) // 0.1;
            , ignoreConvergenceFailure_(true)
            , solverRestartMax_(Parameters::get<TypeTag, Properties::SolverMaxRestarts>()) // 10
            , solverVerbose_(Parameters::get<TypeTag, Properties::SolverVerbosity>() > 0 && terminalOutput) // 2
            , timestepVerbose_(Parameters::get<TypeTag, Properties::TimeStepVerbosity>() > 0 && terminalOutput) // 2
            , suggestedNextTimestep_(max_next_tstep <= 0 ? Parameters::get<TypeTag, Properties::InitialTimeStepInDays>() * 24 * 60 * 60 : max_next_tstep) // 1.0
            , fullTimestepInitially_(Parameters::get<TypeTag, Properties::FullTimeStepInitially>()) // false
            , timestepAfterEvent_(tuning.TMAXWC) // 1e30
            , useNewtonIteration_(false)
            , minTimeStepBeforeShuttingProblematicWells_(Parameters::get<TypeTag, Properties::MinTimeStepBeforeShuttingProblematicWellsInDays>() * unit::day)
        {
            init_(unitSystem);
        }
 
        static void registerParameters()
        {
            registerEclTimeSteppingParameters<TypeTag>();
            // TODO: make sure the help messages are correct (and useful)
            Parameters::registerParam<TypeTag, Properties::SolverContinueOnConvergenceFailure>
                ("Continue instead of stop when minimum solver time step is reached");
            Parameters::registerParam<TypeTag, Properties::SolverMaxRestarts>
                ("The maximum number of breakdowns before a substep is given up and "
                 "the simulator is terminated");
            Parameters::registerParam<TypeTag, Properties::SolverVerbosity>
                ("Specify the \"chattiness\" of the non-linear solver itself");
            Parameters::registerParam<TypeTag, Properties::TimeStepVerbosity>
                ("Specify the \"chattiness\" during the time integration");
            Parameters::registerParam<TypeTag, Properties::InitialTimeStepInDays>
                ("The size of the initial time step in days");
            Parameters::registerParam<TypeTag, Properties::FullTimeStepInitially>
                ("Always attempt to finish a report step using a single substep");
            Parameters::registerParam<TypeTag, Properties::TimeStepControl>
                ("The algorithm used to determine time-step sizes. "
                 "Valid options are: "
                 "'pid' (default), "
                 "'pid+iteration', "
                 "'pid+newtoniteration', "
                 "'iterationcount', "
                "'newtoniterationcount' "
                "and 'hardcoded'");
            Parameters::registerParam<TypeTag, Properties::TimeStepControlTolerance>
                ("The tolerance used by the time step size control algorithm");
            Parameters::registerParam<TypeTag, Properties::TimeStepControlTargetIterations>
                ("The number of linear iterations which the time step control scheme "
                 "should aim for (if applicable)");
            Parameters::registerParam<TypeTag, Properties::TimeStepControlTargetNewtonIterations>
                ("The number of Newton iterations which the time step control scheme "
                 "should aim for (if applicable)");
            Parameters::registerParam<TypeTag, Properties::TimeStepControlDecayRate>
                ("The decay rate of the time step size of the number of "
                 "target iterations is exceeded");
            Parameters::registerParam<TypeTag, Properties::TimeStepControlGrowthRate>
                ("The growth rate of the time step size of the number of "
                 "target iterations is undercut");
            Parameters::registerParam<TypeTag, Properties::TimeStepControlDecayDampingFactor>
                ("The decay rate of the time step decrease when the "
                 "target iterations is exceeded");
            Parameters::registerParam<TypeTag, Properties::TimeStepControlGrowthDampingFactor>
                ("The growth rate of the time step increase when the "
                 "target iterations is undercut");
            Parameters::registerParam<TypeTag, Properties::TimeStepControlFileName>
                ("The name of the file which contains the hardcoded time steps sizes");
            Parameters::registerParam<TypeTag, Properties::MinTimeStepBeforeShuttingProblematicWellsInDays>
                ("The minimum time step size in days for which problematic wells are not shut");
            Parameters::registerParam<TypeTag, Properties::MinTimeStepBasedOnNewtonIterations>
                ("The minimum time step size (in days for field and metric unit and hours for lab unit) "
                 "can be reduced to based on newton iteration counts");
        }
 
        template <class Solver>
        SimulatorReport step(const SimulatorTimer& simulatorTimer,
                             Solver& solver,
                             const bool isEvent,
                             const std::vector<int>* fipnum = nullptr,
                             const std::function<bool()> tuningUpdater = [](){return false;})
        {
            // Maybe update tuning
            tuningUpdater();
            SimulatorReport report;
            const double timestep = simulatorTimer.currentStepLength();
 
            // init last time step as a fraction of the given time step
            if (suggestedNextTimestep_ < 0) {
                suggestedNextTimestep_ = restartFactor_ * timestep;
            }
 
            if (fullTimestepInitially_) {
                suggestedNextTimestep_ = timestep;
            }
 
            // use seperate time step after event
            if (isEvent && timestepAfterEvent_ > 0) {
                suggestedNextTimestep_ = timestepAfterEvent_;
            }
 
            auto& simulator = solver.model().simulator();
            auto& problem = simulator.problem();
 
            // create adaptive step timer with previously used sub step size
            AdaptiveSimulatorTimer substepTimer(simulatorTimer, suggestedNextTimestep_, maxTimeStep_);
 
            // counter for solver restarts
            int restarts = 0;
 
            // sub step time loop
            while (!substepTimer.done()) {
                // Maybe update tuning
                // get current delta t
                auto oldValue = suggestedNextTimestep_;
                if (tuningUpdater()) {
                    // Use provideTimeStepEstimate to make we sure don't simulate longer than the report step is.
                    substepTimer.provideTimeStepEstimate(suggestedNextTimestep_);
                    suggestedNextTimestep_ = oldValue;
                }
                const double dt = substepTimer.currentStepLength();
                if (timestepVerbose_) {
                    detail::logTimer(substepTimer);
                }
 
                SimulatorReportSingle substepReport;
                std::string causeOfFailure;
                try {
                    substepReport = solver.step(substepTimer);
 
                    if (solverVerbose_) {
                        // report number of linear iterations
                        OpmLog::debug("Overall linear iterations used: " + std::to_string(substepReport.total_linear_iterations));
                    }
                }
                catch (const TooManyIterations& e) {
                    substepReport = solver.failureReport();
                    causeOfFailure = "Solver convergence failure - Iteration limit reached";
 
                    logException_(e, solverVerbose_);
                    // since linearIterations is < 0 this will restart the solver
                }
                catch (const LinearSolverProblem& e) {
                    substepReport = solver.failureReport();
                    causeOfFailure = "Linear solver convergence failure";
 
                    logException_(e, solverVerbose_);
                    // since linearIterations is < 0 this will restart the solver
                }
                catch (const NumericalProblem& e) {
                    substepReport = solver.failureReport();
                    causeOfFailure = "Solver convergence failure - Numerical problem encountered";
 
                    logException_(e, solverVerbose_);
                    // since linearIterations is < 0 this will restart the solver
                }
                catch (const std::runtime_error& e) {
                    substepReport = solver.failureReport();
 
                    logException_(e, solverVerbose_);
                    // also catch linear solver not converged
                }
                catch (const Dune::ISTLError& e) {
                    substepReport = solver.failureReport();
 
                    logException_(e, solverVerbose_);
                    // also catch errors in ISTL AMG that occur when time step is too large
                }
                catch (const Dune::MatrixBlockError& e) {
                    substepReport = solver.failureReport();
 
                    logException_(e, solverVerbose_);
                    // this can be thrown by ISTL's ILU0 in block mode, yet is not an ISTLError
                }
 
                //Pass substep to eclwriter for summary output
                simulator.problem().setSubStepReport(substepReport);
 
                report += substepReport;
 
                bool continue_on_uncoverged_solution = ignoreConvergenceFailure_ &&
                                                       !substepReport.converged  &&
                                                       dt <= minTimeStep_;
 
                if (continue_on_uncoverged_solution && solverVerbose_) {
                    const auto msg = fmt::format(
                        "Solver failed to converge but timestep "
                        "{} is smaller or equal to {}\n"
                        "which is the minimum threshold given "
                        "by option --solver-min-time-step\n",
                        dt, minTimeStep_
                    );
                    OpmLog::problem(msg);
                }
 
                if (substepReport.converged || continue_on_uncoverged_solution) {
 
                    // advance by current dt
                    ++substepTimer;
 
                    // create object to compute the time error, simply forwards the call to the model
                    SolutionTimeErrorSolverWrapper<Solver> relativeChange(solver);
 
                    // compute new time step estimate
                    const int iterations = useNewtonIteration_ ? substepReport.total_newton_iterations
                        : substepReport.total_linear_iterations;
                    double dtEstimate = timeStepControl_->computeTimeStepSize(dt, iterations, relativeChange,
                                                                               substepTimer.simulationTimeElapsed());
 
                    assert(dtEstimate > 0);
                    // limit the growth of the timestep size by the growth factor
                    dtEstimate = std::min(dtEstimate, double(maxGrowth_ * dt));
                    assert(dtEstimate > 0);
                    // further restrict time step size growth after convergence problems
                    if (restarts > 0) {
                        dtEstimate = std::min(growthFactor_ * dt, dtEstimate);
                        // solver converged, reset restarts counter
                        restarts = 0;
                    }
 
                    // Further restrict time step size if we are in
                    // prediction mode with THP constraints.
                    if (solver.model().wellModel().hasTHPConstraints()) {
                        const double maxPredictionTHPTimestep = 16.0 * unit::day;
                        dtEstimate = std::min(dtEstimate, maxPredictionTHPTimestep);
                    }
                    assert(dtEstimate > 0);
                    if (timestepVerbose_) {
                        std::ostringstream ss;
                        substepReport.reportStep(ss);
                        OpmLog::info(ss.str());
                    }
 
                    // write data if outputWriter was provided
                    // if the time step is done we do not need
                    // to write it as this will be done by the simulator
                    // anyway.
                    if (!substepTimer.done()) {
                        if (fipnum) {
                            solver.computeFluidInPlace(*fipnum);
                        }
                        time::StopWatch perfTimer;
                        perfTimer.start();
 
                        problem.writeOutput(simulatorTimer);
 
                        report.success.output_write_time += perfTimer.secsSinceStart();
                    }
 
                    // set new time step length
                    substepTimer.provideTimeStepEstimate(dtEstimate);
 
                    report.success.converged = substepTimer.done();
                    substepTimer.setLastStepFailed(false);
 
                }
                else { // in case of no convergence
                    substepTimer.setLastStepFailed(true);
 
                    // If we have restarted (i.e. cut the timestep) too
                    // many times, we have failed and throw an exception.
                    if (restarts >= solverRestartMax_) {
                        const auto msg = fmt::format(
                            "Solver failed to converge after cutting timestep {} times.",
                            restarts
                        );
                        if (solverVerbose_) {
                            OpmLog::error(msg);
                        }
                        // Use throw directly to prevent file and line
                        throw TimeSteppingBreakdown{msg};
                    }
 
                    // The new, chopped timestep.
                    const double newTimeStep = restartFactor_ * dt;
 
 
                    // If we have restarted (i.e. cut the timestep) too
                    // much, we have failed and throw an exception.
                    if (newTimeStep < minTimeStep_) {
                        const auto msg = fmt::format(
                            "Solver failed to converge after cutting timestep to {}\n"
                            "which is the minimum threshold given by option --solver-min-time-step\n",
                            minTimeStep_
                        );
                        if (solverVerbose_) {
                            OpmLog::error(msg);
                        }
                        // Use throw directly to prevent file and line
                        throw TimeSteppingBreakdown{msg};
                    }
 
                    // Define utility function for chopping timestep.
                    auto chopTimestep = [&]() {
                        substepTimer.provideTimeStepEstimate(newTimeStep);
                        if (solverVerbose_) {
                            const auto msg = fmt::format(
                                "{}\nTimestep chopped to {} days\n",
                                causeOfFailure,
                                std::to_string(unit::convert::to(substepTimer.currentStepLength(), unit::day))
                            );
                            OpmLog::problem(msg);
                        }
                        ++restarts;
                    };
 
                    const double minimumChoppedTimestep = minTimeStepBeforeShuttingProblematicWells_;
                    if (newTimeStep > minimumChoppedTimestep) {
                        chopTimestep();
                    } else {
                        // We are below the threshold, and will check if there are any
                        // wells we should close rather than chopping again.
                        std::set<std::string> failing_wells = detail::consistentlyFailingWells(solver.model().stepReports());
                        if (failing_wells.empty()) {
                            // Found no wells to close, chop the timestep as above.
                            chopTimestep();
                        } else {
                            // Close all consistently failing wells.
                            int num_shut_wells = 0;
                            for (const auto& well : failing_wells) {
                                bool was_shut = solver.model().wellModel().forceShutWellByName(well, substepTimer.simulationTimeElapsed());
                                if (was_shut) {
                                    ++num_shut_wells;
                                }
                            }
                            if (num_shut_wells == 0) {
                                // None of the problematic wells were shut.
                                // We must fall back to chopping again.
                                chopTimestep();
                            } else {
                                substepTimer.provideTimeStepEstimate(dt);
                                if (solverVerbose_) {
                                    std::string msg;
                                    msg = "\nProblematic well(s) were shut: ";
                                    for (const auto& well : failing_wells) {
                                        msg += well;
                                        msg += " ";
                                    }
                                    msg += "(retrying timestep)\n";
                                    OpmLog::problem(msg);
                                }
                            }
                        }
                    }
                }
                problem.setNextTimeStepSize(substepTimer.currentStepLength());
            }
 
            // store estimated time step for next reportStep
            suggestedNextTimestep_ = substepTimer.currentStepLength();
            if (timestepVerbose_) {
                std::ostringstream ss;
                substepTimer.report(ss);
                ss << "Suggested next step size = " << unit::convert::to(suggestedNextTimestep_, unit::day) << " (days)" << std::endl;
                OpmLog::debug(ss.str());
            }
 
            if (! std::isfinite(suggestedNextTimestep_)) { // check for NaN
                suggestedNextTimestep_ = timestep;
            }
            return report;
        }
 
        double suggestedNextStep() const
        { return suggestedNextTimestep_; }
 
        void setSuggestedNextStep(const double x)
        { suggestedNextTimestep_ = x; }
 
        void updateTUNING(double max_next_tstep, const Tuning& tuning)
        {
            restartFactor_ = tuning.TSFCNV;
            growthFactor_ = tuning.TFDIFF;
            maxGrowth_ = tuning.TSFMAX;
            maxTimeStep_ = tuning.TSMAXZ;
            updateNEXTSTEP(max_next_tstep);
            timestepAfterEvent_ = tuning.TMAXWC;
        }
 
        void updateNEXTSTEP(double max_next_tstep)
        {
             // \Note Only update next suggested step if TSINIT was explicitly set in TUNING or NEXTSTEP is active. 
            if (max_next_tstep > 0) {
                suggestedNextTimestep_ = max_next_tstep;
            }
        }
 
        template<class Serializer>
        void serializeOp(Serializer& serializer)
        {
            serializer(timeStepControlType_);
            switch (timeStepControlType_) {
            case TimeStepControlType::HardCodedTimeStep:
                allocAndSerialize<HardcodedTimeStepControl>(serializer);
                break;
            case TimeStepControlType::PIDAndIterationCount:
                allocAndSerialize<PIDAndIterationCountTimeStepControl>(serializer);
                break;
            case TimeStepControlType::SimpleIterationCount:
                allocAndSerialize<SimpleIterationCountTimeStepControl>(serializer);
                break;
            case TimeStepControlType::PID:
                allocAndSerialize<PIDTimeStepControl>(serializer);
                break;
            }
            serializer(restartFactor_);
            serializer(growthFactor_);
            serializer(maxGrowth_);
            serializer(maxTimeStep_);
            serializer(minTimeStep_);
            serializer(ignoreConvergenceFailure_);
            serializer(solverRestartMax_);
            serializer(solverVerbose_);
            serializer(timestepVerbose_);
            serializer(suggestedNextTimestep_);
            serializer(fullTimestepInitially_);
            serializer(timestepAfterEvent_);
            serializer(useNewtonIteration_);
            serializer(minTimeStepBeforeShuttingProblematicWells_);
        }
 
        static AdaptiveTimeStepping<TypeTag> serializationTestObjectHardcoded()
        {
            return serializationTestObject_<HardcodedTimeStepControl>();
        }
 
        static AdaptiveTimeStepping<TypeTag> serializationTestObjectPID()
        {
            return serializationTestObject_<PIDTimeStepControl>();
        }
 
        static AdaptiveTimeStepping<TypeTag> serializationTestObjectPIDIt()
        {
            return serializationTestObject_<PIDAndIterationCountTimeStepControl>();
        }
 
        static AdaptiveTimeStepping<TypeTag> serializationTestObjectSimple()
        {
            return serializationTestObject_<SimpleIterationCountTimeStepControl>();
        }
 
        bool operator==(const AdaptiveTimeStepping<TypeTag>& rhs)
        {
            if (timeStepControlType_ != rhs.timeStepControlType_ ||
                (timeStepControl_ && !rhs.timeStepControl_) ||
                (!timeStepControl_ && rhs.timeStepControl_)) {
                return false;
            }
 
            bool result = false;
            switch (timeStepControlType_) {
            case TimeStepControlType::HardCodedTimeStep:
                result = castAndComp<HardcodedTimeStepControl>(rhs);
                break;
            case TimeStepControlType::PIDAndIterationCount:
                result = castAndComp<PIDAndIterationCountTimeStepControl>(rhs);
                break;
            case TimeStepControlType::SimpleIterationCount:
                result = castAndComp<SimpleIterationCountTimeStepControl>(rhs);
                break;
            case TimeStepControlType::PID:
                result = castAndComp<PIDTimeStepControl>(rhs);
                break;
            }
 
            return result &&
                   this->restartFactor_ == rhs.restartFactor_ &&
                   this->growthFactor_ == rhs.growthFactor_ &&
                   this->maxGrowth_ == rhs.maxGrowth_ &&
                   this->maxTimeStep_ == rhs.maxTimeStep_ &&
                   this->minTimeStep_ == rhs.minTimeStep_ &&
                   this->ignoreConvergenceFailure_ == rhs.ignoreConvergenceFailure_ &&
                   this->solverRestartMax_== rhs.solverRestartMax_ &&
                   this->solverVerbose_ == rhs.solverVerbose_ &&
                   this->fullTimestepInitially_ == rhs.fullTimestepInitially_ &&
                   this->timestepAfterEvent_ == rhs.timestepAfterEvent_ &&
                   this->useNewtonIteration_ == rhs.useNewtonIteration_ &&
                   this->minTimeStepBeforeShuttingProblematicWells_ ==
                       rhs.minTimeStepBeforeShuttingProblematicWells_;
        }
 
    private:
        template<class Controller>
        static AdaptiveTimeStepping<TypeTag> serializationTestObject_()
        {
            AdaptiveTimeStepping<TypeTag> result;
 
            result.restartFactor_ = 1.0;
            result.growthFactor_ = 2.0;
            result.maxGrowth_ = 3.0;
            result.maxTimeStep_ = 4.0;
            result.minTimeStep_ = 5.0;
            result.ignoreConvergenceFailure_ = true;
            result.solverRestartMax_ = 6;
            result.solverVerbose_ = true;
            result.timestepVerbose_ = true;
            result.suggestedNextTimestep_ = 7.0;
            result.fullTimestepInitially_ = true;
            result.useNewtonIteration_ = true;
            result.minTimeStepBeforeShuttingProblematicWells_ = 9.0;
            result.timeStepControlType_ = Controller::Type;
            result.timeStepControl_ = std::make_unique<Controller>(Controller::serializationTestObject());
 
            return result;
        }
        template<class T, class Serializer>
        void allocAndSerialize(Serializer& serializer)
        {
            if (!serializer.isSerializing()) {
                timeStepControl_ = std::make_unique<T>();
            }
            serializer(*static_cast<T*>(timeStepControl_.get()));
        }
 
        template<class T>
        bool castAndComp(const AdaptiveTimeStepping<TypeTag>& Rhs) const
        {
            const T* lhs = static_cast<const T*>(timeStepControl_.get());
            const T* rhs = static_cast<const T*>(Rhs.timeStepControl_.get());
            return *lhs == *rhs;
        }
 
    protected:
        void init_(const UnitSystem& unitSystem)
        {
            // valid are "pid" and "pid+iteration"
            std::string control = Parameters::get<TypeTag, Properties::TimeStepControl>(); // "pid"
 
            const double tol =  Parameters::get<TypeTag, Properties::TimeStepControlTolerance>(); // 1e-1
            if (control == "pid") {
                timeStepControl_ = std::make_unique<PIDTimeStepControl>(tol);
                timeStepControlType_ = TimeStepControlType::PID;
            }
            else if (control == "pid+iteration") {
                const int iterations =  Parameters::get<TypeTag, Properties::TimeStepControlTargetIterations>(); // 30
                const double decayDampingFactor = Parameters::get<TypeTag, Properties::TimeStepControlDecayDampingFactor>(); // 1.0
                const double growthDampingFactor = Parameters::get<TypeTag, Properties::TimeStepControlGrowthDampingFactor>(); // 3.2
                timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor, growthDampingFactor, tol);
                timeStepControlType_ = TimeStepControlType::PIDAndIterationCount;
            }
            else if (control == "pid+newtoniteration") {
                const int iterations =  Parameters::get<TypeTag, Properties::TimeStepControlTargetNewtonIterations>(); // 8
                const double decayDampingFactor = Parameters::get<TypeTag, Properties::TimeStepControlDecayDampingFactor>(); // 1.0
                const double growthDampingFactor = Parameters::get<TypeTag, Properties::TimeStepControlGrowthDampingFactor>(); // 3.2
                const double nonDimensionalMinTimeStepIterations = Parameters::get<TypeTag, Properties::MinTimeStepBasedOnNewtonIterations>(); // 0.0 by default
                // the min time step can be reduced by the newton iteration numbers
                double minTimeStepReducedByIterations = unitSystem.to_si(UnitSystem::measure::time, nonDimensionalMinTimeStepIterations);
                timeStepControl_ = std::make_unique<PIDAndIterationCountTimeStepControl>(iterations, decayDampingFactor,
                                                                                         growthDampingFactor, tol, minTimeStepReducedByIterations);
                timeStepControlType_ = TimeStepControlType::PIDAndIterationCount;
                useNewtonIteration_ = true;
            }
            else if (control == "iterationcount") {
                const int iterations =  Parameters::get<TypeTag, Properties::TimeStepControlTargetIterations>(); // 30
                const double decayrate = Parameters::get<TypeTag, Properties::TimeStepControlDecayRate>(); // 0.75
                const double growthrate = Parameters::get<TypeTag, Properties::TimeStepControlGrowthRate>(); // 1.25
                timeStepControl_ = std::make_unique<SimpleIterationCountTimeStepControl>(iterations, decayrate, growthrate);
                timeStepControlType_ = TimeStepControlType::SimpleIterationCount;
            }
            else if (control == "newtoniterationcount") {
                const int iterations =  Parameters::get<TypeTag, Properties::TimeStepControlTargetNewtonIterations>(); // 8
                const double decayrate = Parameters::get<TypeTag, Properties::TimeStepControlDecayRate>(); // 0.75
                const double growthrate = Parameters::get<TypeTag, Properties::TimeStepControlGrowthRate>(); // 1.25
                timeStepControl_ = std::make_unique<SimpleIterationCountTimeStepControl>(iterations, decayrate, growthrate);
                useNewtonIteration_ = true;
                timeStepControlType_ = TimeStepControlType::SimpleIterationCount;
            }
            else if (control == "hardcoded") {
                const std::string filename = Parameters::get<TypeTag, Properties::TimeStepControlFileName>(); // "timesteps"
                timeStepControl_ = std::make_unique<HardcodedTimeStepControl>(filename);
                timeStepControlType_ = TimeStepControlType::HardCodedTimeStep;
            }
            else
                OPM_THROW(std::runtime_error,
                          "Unsupported time step control selected " + control);
 
            // make sure growth factor is something reasonable
            assert(growthFactor_ >= 1.0);
        }
 
        using TimeStepController = std::unique_ptr<TimeStepControlInterface>;
 
        TimeStepControlType timeStepControlType_; 
        TimeStepController timeStepControl_; 
        double restartFactor_;               
        double growthFactor_;                
        double maxGrowth_;                   
        double maxTimeStep_;                 
        double minTimeStep_;                 
        bool ignoreConvergenceFailure_;      
        int solverRestartMax_;               
        bool solverVerbose_;                 
        bool timestepVerbose_;               
        double suggestedNextTimestep_;       
        bool fullTimestepInitially_;         
        double timestepAfterEvent_;          
        bool useNewtonIteration_;            
        double minTimeStepBeforeShuttingProblematicWells_; 
    };
}
 
#endif // OPM_ADAPTIVE_TIME_STEPPING_HPP