Comparing fixed timestep only+failed update

compareECLFiles.cmake

@@ -60,7 +60,7 @@ endfunction()
 #   - This test class compares output from a simulation to reference files.
 function(add_test_compareECLFiles)
   set(oneValueArgs CASENAME FILENAME SIMULATOR ABS_TOL REL_TOL DIR DIR_PREFIX PREFIX RESTART_STEP RESTART_SCHED)
-  set(multiValueArgs TEST_ARGS)
+  set(multiValueArgs TEST_ARGS TEST_ARGS_REPLAY)
   cmake_parse_arguments(PARAM "$" "${oneValueArgs}" "${multiValueArgs}" ${ARGN} )
   if(NOT PARAM_DIR)
     set(PARAM_DIR ${PARAM_CASENAME})
@@ -84,6 +84,9 @@ function(add_test_compareECLFiles)
   if(PARAM_RESTART_SCHED STREQUAL "false" OR PARAM_RESTART_SCHED STREQUAL "true")
     list(APPEND DRIVER_ARGS -h ${PARAM_RESTART_SCHED})
   endif()
+  foreach(arg IN LISTS PARAM_TEST_ARGS_REPLAY)
+    list(APPEND DRIVER_ARGS -y ${arg})
+  endforeach()
   opm_add_test(${PARAM_PREFIX}_${PARAM_SIMULATOR}+${PARAM_FILENAME} NO_COMPILE
                EXE_NAME ${PARAM_SIMULATOR}
                DRIVER_ARGS ${DRIVER_ARGS}

flowexperimental/comp/wells/CompWellModel.hpp

@@ -96,6 +96,11 @@ class CompWellModel : WellConnectionAuxiliaryModule<TypeTag, CompWellModel<TypeT
     void beginTimeStep();
     void beginIteration();
 
+     // FlowProblem calls these lifecycle hooks for well models.
+     // Compositional wells currently keep no extra timestep snapshot state.
+     void updateFailed() {}
+     void advanceTimeLevel() {}
+
     void init();
     void endIteration() const {}
     void endTimeStep() {}

opm/models/blackoil/blackoilnewtonmethod.hpp

@@ -90,6 +90,12 @@ class BlackOilNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonM
         wasSwitched_.resize(this->model().numTotalDof(), false);
     }
 
+    void resetPrimaryVariableSwitches()
+    {
+        numPriVarsSwitched_ = 0;
+        std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
+    }
+
     /*!
      * \brief Register all run-time parameters for the blackoil newton method.
      */

opm/models/utils/simulator.hh

@@ -395,7 +395,9 @@ public:
      * \param value The new value for the time step size \f$\mathrm{[s]}\f$
      */
     void setTimeStepSize(Scalar value)
-    { timeStepSize_ = value; }
+    { 
+        timeStepSize_ = float(value); 
+    }
 
     /*!
      * \brief Set the current time step index to a given value.

opm/simulators/flow/BlackoilModel_impl.hpp

@@ -122,10 +122,11 @@ prepareStep(const SimulatorTimerInterface& timer)
                                 "- the previous step succeeded on some ranks but failed on others.");
     }
     if (lastStepFailed) {
-        simulator_.model().updateFailed();
+        simulator_.problem().updateFailed();
+        simulator_.model().newtonMethod().eraseMatrix();
     }
     else {
-        simulator_.model().advanceTimeLevel();
+        simulator_.problem().advanceTimeLevel();
     }
 
     // Set the timestep size and episode index for the model explicitly.
@@ -142,6 +143,7 @@ prepareStep(const SimulatorTimerInterface& timer)
     unsigned numDof = simulator_.model().numGridDof();
     wasSwitched_.resize(numDof);
     std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
+    simulator_.model().newtonMethod().resetPrimaryVariableSwitches();
 
     if (param_.update_equations_scaling_) {
         OpmLog::error("Equation scaling not supported");

opm/simulators/flow/FlowProblem.hpp

@@ -367,6 +367,8 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
         const int episodeIdx = this->episodeIndex();
         const int timeStepSize = this->simulator().timeStepSize();
 
+        this->captureBeginTimeStepState_();
+
         this->beginTimeStep_(enableExperiments,
                              episodeIdx,
                              this->simulator().timeStepIndex(),
@@ -392,6 +394,19 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
 
     }
 
+    void updateFailed()
+    {
+        this->restoreBeginTimeStepState_();
+        wellModel_.updateFailed();
+        this->model().updateFailed();
+    }
+
+    void advanceTimeLevel()
+    {
+        this->model().advanceTimeLevel();
+        wellModel_.advanceTimeLevel();
+    }
+
     /*!
      * \brief Called by the simulator before each Newton-Raphson iteration.
      */
@@ -1271,6 +1286,26 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     { return *static_cast<const Implementation *>(this); }
 
 protected:
+    virtual void captureBeginTimeStepState_()
+    {
+        prev_timestep_first_step_ = first_step_;
+        prev_timestep_max_polymer_adsorption_ = this->polymer_.maxAdsorption;
+        prev_timestep_max_oil_saturation_ = this->maxOilSaturation_;
+        prev_timestep_max_water_saturation_ = this->maxWaterSaturation_;
+        prev_timestep_min_ref_pressure_ = this->minRefPressure_;
+        prev_timestep_rock_comp_trans_mult_val_ = this->rockCompTransMultVal_;
+    }
+
+    virtual void restoreBeginTimeStepState_()
+    {
+        first_step_ = prev_timestep_first_step_;
+        this->polymer_.maxAdsorption = prev_timestep_max_polymer_adsorption_;
+        this->maxOilSaturation_ = prev_timestep_max_oil_saturation_;
+        this->maxWaterSaturation_ = prev_timestep_max_water_saturation_;
+        this->minRefPressure_ = prev_timestep_min_ref_pressure_;
+        this->rockCompTransMultVal_ = prev_timestep_rock_comp_trans_mult_val_;
+    }
+
     template<class UpdateFunc>
     void updateProperty_(const std::string& failureMsg,
                          UpdateFunc func)
@@ -1857,6 +1892,13 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     BCData<int> bcindex_;
     bool nonTrivialBoundaryConditions_ = false;
     bool first_step_ = true;
+    bool prev_timestep_first_step_ = true;
+
+    std::vector<Scalar> prev_timestep_max_polymer_adsorption_;
+    std::vector<Scalar> prev_timestep_max_oil_saturation_;
+    std::vector<Scalar> prev_timestep_max_water_saturation_;
+    std::vector<Scalar> prev_timestep_min_ref_pressure_;
+    std::vector<Scalar> prev_timestep_rock_comp_trans_mult_val_;
 
     /// Whether or not the current episode will end at the end of the
     /// current time step.

opm/simulators/flow/FlowProblemBlackoil.hpp

@@ -188,6 +188,7 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
         : FlowProblemType(simulator)
         , thresholdPressures_(simulator)
         , mixControls_(simulator.vanguard().schedule())
+        , prev_timestep_mixControls_(simulator.vanguard().schedule())
         , actionHandler_(simulator.vanguard().eclState(),
                          simulator.vanguard().schedule(),
                          simulator.vanguard().actionState(),
@@ -1173,6 +1174,18 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
     }
 
 protected:
+    void captureBeginTimeStepState_() override
+    {
+        FlowProblemType::captureBeginTimeStepState_();
+        prev_timestep_mixControls_ = mixControls_;
+    }
+
+    void restoreBeginTimeStepState_() override
+    {
+        FlowProblemType::restoreBeginTimeStepState_();
+        mixControls_ = prev_timestep_mixControls_;
+    }
+
     void updateExplicitQuantities_(int episodeIdx, int timeStepSize, const bool first_step_after_restart) override
     {
         this->updateExplicitQuantities_(first_step_after_restart);
@@ -1734,6 +1747,7 @@ class FlowProblemBlackoil : public FlowProblem<TypeTag>
     std::unique_ptr<DamarisWriterType> damarisWriter_;
 #endif
     MixingRateControls<FluidSystem> mixControls_;
+    MixingRateControls<FluidSystem> prev_timestep_mixControls_;
 
     ActionHandler<Scalar, IndexTraits> actionHandler_;
 

opm/simulators/linalg/ISTLSolver.hpp

@@ -313,9 +313,18 @@ std::unique_ptr<Matrix> blockJacobiAdjacency(const Grid& grid,
             element_chunks_ = std::make_unique<ElementChunksType>(simulator_.vanguard().gridView(), Dune::Partitions::all, ThreadManager::maxThreads());
         }
 
-        // nothing to clean here
         void eraseMatrix() override
         {
+            matrix_ = nullptr;
+            rhs_ = nullptr;
+            iterations_ = 0;
+            solveCount_ = 0;
+
+            for (auto& solverInfo : flexibleSolver_) {
+                solverInfo.pre_ = nullptr;
+                solverInfo.solver_.reset();
+                solverInfo.op_.reset();
+            }
         }
 
         void setActiveSolver(const int num) override

opm/simulators/linalg/ISTLSolverTPSA.hpp

@@ -274,7 +274,15 @@ class ISTLSolverTPSA : public AbstractISTLSolver<GetPropType<TypeTag, Properties
     * \copydoc AbstractISTLSolver::eraseMatrix
     */
     void eraseMatrix() override
-    { }
+    {
+        matrix_ = nullptr;
+        rhs_ = nullptr;
+        iterations_ = 0;
+        solveCount_ = 0;
+        flexibleSolver_.pre_ = nullptr;
+        flexibleSolver_.solver_.reset();
+        flexibleSolver_.op_.reset();
+    }
 
     /*!
     * \copydoc AbstractISTLSolver::setActiveSolver

opm/simulators/linalg/gpuistl/ISTLSolverGPUISTL.hpp

@@ -141,12 +141,21 @@ class ISTLSolverGPUISTL : public AbstractISTLSolver<GetPropType<TypeTag, Propert
 
     /**
      * \copydoc AbstractISTLSolver::eraseMatrix
-     *
-     * \note This method will not do anything.
      */
     void eraseMatrix() override
     {
-        // Nothing, this is the same as the ISTLSolver
+        m_lastSeenIterations = 0;
+        m_solveCount = 0;
+        m_matrix.reset();
+        m_gpuSolver.reset();
+        m_rhs.reset();
+        m_x.reset();
+        m_pinnedMatrixMemory.reset();
+        m_pinnedRhsMemory.reset();
+        m_pinnedXMemory.reset();
+        m_pinnedWeightsMemory.reset();
+        m_weights.reset();
+        m_diagonalIndices.reset();
     }
 
     /**

opm/simulators/timestepping/AdaptiveTimeStepping_impl.hpp

@@ -522,6 +522,34 @@ maybeModifySuggestedTimeStepAtBeginningOfReportStep_(const double original_time_
     this->adaptive_time_stepping_.maybeModifySuggestedTimeStepAtBeginningOfReportStep_(
         original_time_step, this->is_event_
     );
+
+    if (this->adaptive_time_stepping_.time_step_control_type_ != TimeStepControlType::HardCodedTimeStep) {
+        return;
+    }
+
+    struct ZeroRelativeChange final : RelativeChangeInterface {
+        double relativeChange() const override { return 0.0; }
+    } zero_relative_change;
+
+    const auto* hardcoded_control = static_cast<const HardcodedTimeStepControl*>(
+        this->adaptive_time_stepping_.time_step_control_.get());
+    AdaptiveSimulatorTimer report_step_timer{
+        this->simulator_timer_.startDateTime(),
+        original_time_step,
+        this->simulator_timer_.simulationTimeElapsed(),
+        original_time_step,
+        this->simulator_timer_.reportStepNum(),
+        maxTimeStep_()
+    };
+
+    const double hardcoded_initial_step = hardcoded_control->computeTimeStepSize(
+        original_time_step,
+        0,
+        zero_relative_change,
+        report_step_timer);
+    if (std::isfinite(hardcoded_initial_step) && hardcoded_initial_step > 0.0) {
+        this->adaptive_time_stepping_.setSuggestedNextStep(hardcoded_initial_step);
+    }
 }
 
 // The maybeUpdateTuning_() lambda callback is defined in SimulatorFullyImplicitBlackoil::runStep()

opm/simulators/timestepping/SimulatorReport.cpp

@@ -395,7 +395,7 @@ namespace Opm
         os << "  Time(day)  TStep(day)  Assembly    LSetup    LSolve    LocSol    Update    Output WellIt Lins NewtIt LinIt Conv\n";
         for (std::size_t i = 0; i < this->stepreports.size(); ++i) {
             const SimulatorReportSingle& sr = this->stepreports[i];
-            os.precision(10);
+            os.precision(20);
             os << std::defaultfloat;
             os << std::setw(11) << unit::convert::to(sr.global_time, unit::day) << " ";
             os << std::setw(11) << unit::convert::to(sr.timestep_length, unit::day) << " ";

opm/simulators/timestepping/TimeStepControl.cpp

@@ -41,6 +41,14 @@
 
 namespace Opm
 {
+    namespace {
+    double hardcodedTimeTolerance(const double a, const double b, const double c = 0.0)
+    {
+        const double scale = std::max({1.0, std::abs(a), std::abs(b), std::abs(c)});
+        return 64.0 * std::numeric_limits<double>::epsilon() * scale;
+    }
+    }
+
     ////////////////////////////////////////////////////////
     //
     //  InterationCountTimeStepControl Implementation
@@ -128,12 +136,14 @@ namespace Opm
         std::string::size_type sz;
         std::string line;
         while ( std::getline(infile, line)) {
-            if( line[0] != '-') { // ignore lines starting with '-'
+            if (!line.empty() && line[0] != '-') { // ignore lines starting with '-'
                 const double time = std::stod(line,&sz); // read the first number i.e. the actual substep time
                 subStepTime_.push_back( time * unit::day );
             }
 
         }
+
+        std::sort(subStepTime_.begin(), subStepTime_.end());
     }
 
     HardcodedTimeStepControl HardcodedTimeStepControl::serializationTestObject()
@@ -145,14 +155,25 @@ namespace Opm
     }
 
     double
-    HardcodedTimeStepControl::computeTimeStepSize(const double /*dt */,
+    HardcodedTimeStepControl::computeTimeStepSize(const double dt,
                                                   const int /*iterations */,
                                                   const RelativeChangeInterface& /* relativeChange */,
                                                   const AdaptiveSimulatorTimer& substepTimer) const
     {
-        auto nextTime
-            = std::upper_bound(subStepTime_.begin(), subStepTime_.end(), substepTimer.simulationTimeElapsed());
-        return (*nextTime - substepTimer.simulationTimeElapsed());
+        const double currentTime = substepTimer.simulationTimeElapsed();
+        const double remaining = substepTimer.totalTime() - currentTime;
+        const double tol = hardcodedTimeTolerance(currentTime, substepTimer.totalTime(), dt);
+
+        auto nextTime = std::upper_bound(subStepTime_.begin(), subStepTime_.end(), currentTime + tol);
+        while (nextTime != subStepTime_.end() && (*nextTime - currentTime) <= tol) {
+            ++nextTime;
+        }
+
+        if (nextTime == subStepTime_.end()) {
+            return remaining > tol ? remaining : std::max(dt, tol);
+        }
+
+        return std::max(*nextTime - currentTime, tol);
     }
 
     bool HardcodedTimeStepControl::operator==(const HardcodedTimeStepControl& ctrl) const

opm/simulators/wells/BlackoilWellModel.hpp

@@ -154,6 +154,10 @@ template<class Scalar> class WellContributions;
 
             void beginTimeStep();
 
+            void updateFailed();
+
+            void advanceTimeLevel();
+
             void beginIteration()
             {
                 OPM_TIMEBLOCK(beginIteration);

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -108,9 +108,12 @@ BlackoilWellModelGeneric(Schedule& schedule,
     , terminal_output_(comm_.rank() == 0 &&
                        Parameters::Get<Parameters::EnableTerminalOutput>())
     , guideRate_(schedule)
+    , prev_timestep_guideRate_(schedule)
     , active_wgstate_(pu)
     , last_valid_wgstate_(pu)
+    , prev_timestep_wgstate_(pu)
     , nupcol_wgstate_(pu)
+    , prev_timestep_nupcol_wgstate_(pu)
     , group_state_helper_(this->wellState(),
                           this->groupState(),
                           this->schedule(),
@@ -1964,10 +1967,12 @@ operator==(const BlackoilWellModelGeneric& rhs) const
         && this->last_run_wellpi_ == rhs.last_run_wellpi_
         && this->local_shut_wells_ == rhs.local_shut_wells_
         && this->closed_this_step_ == rhs.closed_this_step_
+        && this->prev_timestep_closed_this_step_ == rhs.prev_timestep_closed_this_step_
         && this->genNetwork_ == rhs.genNetwork_
         && this->prev_inj_multipliers_ == rhs.prev_inj_multipliers_
         && this->active_wgstate_ == rhs.active_wgstate_
         && this->last_valid_wgstate_ == rhs.last_valid_wgstate_
+        && this->prev_timestep_wgstate_ == rhs.prev_timestep_wgstate_
         && this->nupcol_wgstate_ == rhs.nupcol_wgstate_
         && this->switched_prod_groups_ == rhs.switched_prod_groups_
         && this->switched_inj_groups_ == rhs.switched_inj_groups_