Limit timestep when we have THP constraints present.

Using 16 days as the current limit. Could be made a user parameter.

opm/autodiff/BlackoilWellModel.hpp

@@ -236,6 +236,9 @@ namespace Opm {
             // called at the beginning of a report step
             void beginReportStep(const int time_step);
 
+            /// Return true if any well has a THP constraint.
+            bool hasTHPConstraints() const;
+
         protected:
 
             void extractLegacyPressure_(std::vector<double>& cellPressure) const

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -116,6 +116,23 @@ namespace Opm {
     }
 
 
+    /// Return true if any well has a THP constraint.
+    template<typename TypeTag>
+    bool
+    BlackoilWellModel<TypeTag>::
+    hasTHPConstraints() const
+    {
+        for (const auto& well : well_container_) {
+            if (well->wellHasTHPConstraints()) {
+                return true;
+            }
+        }
+        return false;
+    }
+
+
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::

opm/autodiff/WellInterface.hpp

@@ -228,6 +228,9 @@ namespace Opm
         // whether the well is operable
         bool isOperable() const;
 
+        /// Returns true if the well has one or more THP limits/constraints.
+        bool wellHasTHPConstraints() const;
+
     protected:
 
         // to indicate a invalid completion
@@ -324,8 +327,6 @@ namespace Opm
 
         bool underPredictionMode() const;
 
-        bool wellHasTHPConstraints() const;
-
         double getTHPConstraint() const;
 
         int getTHPControlIndex() const;

opm/simulators/timestepping/AdaptiveTimeSteppingEbos.hpp

@@ -298,6 +298,13 @@ namespace Opm {
                         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);
+                    }
+
                     if (timestepVerbose_) {
                         std::ostringstream ss;
                         substepReport.reportStep(ss);