Merge pull request #3539 from totto82/wellOperIter

stop updating operability during iterations

opm/simulators/flow/BlackoilModelParametersEbos.hpp

@@ -107,7 +107,10 @@ template<class TypeTag, class MyTypeTag>
 struct EnableWellOperabilityCheck {
     using type = UndefinedProperty;
 };
-
+template<class TypeTag, class MyTypeTag>
+struct EnableWellOperabilityCheckIter {
+    using type = UndefinedProperty;
+};
 // parameters for multisegment wells
 template<class TypeTag, class MyTypeTag>
 struct TolerancePressureMsWells {
@@ -284,6 +287,10 @@ struct EnableWellOperabilityCheck<TypeTag, TTag::FlowModelParameters> {
     static constexpr bool value = true;
 };
 template<class TypeTag>
+struct EnableWellOperabilityCheckIter<TypeTag, TTag::FlowModelParameters> {
+    static constexpr bool value = false;
+};
+template<class TypeTag>
 struct RelaxedWellFlowTol<TypeTag, TTag::FlowModelParameters> {
     using type = GetPropType<TypeTag, Scalar>;
     static constexpr type value = 1;
@@ -398,6 +405,12 @@ namespace Opm
         // Whether to add influences of wells between cells to the matrix and preconditioner matrix
         bool matrix_add_well_contributions_;
 
+        // Whether to check well operability
+        bool check_well_operability_;
+        // Whether to check well operability during iterations
+        bool check_well_operability_iter_;
+
+
         /// Construct from user parameters or defaults.
         BlackoilModelParametersEbos()
         {
@@ -429,6 +442,8 @@ namespace Opm
             update_equations_scaling_ = EWOMS_GET_PARAM(TypeTag, bool, UpdateEquationsScaling);
             use_update_stabilization_ = EWOMS_GET_PARAM(TypeTag, bool, UseUpdateStabilization);
             matrix_add_well_contributions_ = EWOMS_GET_PARAM(TypeTag, bool, MatrixAddWellContributions);
+            check_well_operability_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheck);
+            check_well_operability_iter_ = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter);
 
             deck_file_name_ = EWOMS_GET_PARAM(TypeTag, std::string, EclDeckFileName);
         }
@@ -466,6 +481,7 @@ namespace Opm
             EWOMS_REGISTER_PARAM(TypeTag, bool, UseUpdateStabilization, "Try to detect and correct oscillations or stagnation during the Newton method");
             EWOMS_REGISTER_PARAM(TypeTag, bool, MatrixAddWellContributions, "Explicitly specify the influences of wells between cells in the Jacobian and preconditioner matrices");
             EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheck, "Enable the well operability checking");
+            EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWellOperabilityCheckIter, "Enable the well operability checking during iterations");
         }
     };
 } // namespace Opm

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1136,7 +1136,7 @@ namespace Opm {
         ConvergenceReport local_report;
         const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
         for (const auto& well : well_container_) {
-            if (well->isOperable() ) {
+            if (well->isOperableAndSolvable() ) {
                 local_report += well->getWellConvergence(this->wellState(), B_avg, local_deferredLogger, iterationIdx > param_.strict_outer_iter_wells_ );
             }
         }
@@ -1178,7 +1178,7 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     calculateExplicitQuantities(DeferredLogger& deferred_logger) const
     {
-        // TODO: checking isOperable() ?
+        // TODO: checking isOperableAndSolvable() ?
         for (auto& well : well_container_) {
             well->calculateExplicitQuantities(ebosSimulator_, this->wellState(), deferred_logger);
         }
@@ -1285,7 +1285,7 @@ namespace Opm {
         for (const auto& well : well_container_) {
             const auto& wname = well->name();
             const auto wasClosed = wellTestState.hasWellClosed(wname);
-
+            well->checkWellOperability(ebosSimulator_, this->wellState(), local_deferredLogger);
             well->updateWellTestState(this->wellState().well(wname), simulationTime, /*writeMessageToOPMLog=*/ true, wellTestState, local_deferredLogger);
 
             if (!wasClosed && wellTestState.hasWellClosed(wname)) {
@@ -1394,11 +1394,6 @@ namespace Opm {
     prepareTimeStep(DeferredLogger& deferred_logger)
     {
         for (const auto& well : well_container_) {
-            const bool old_well_operable = well->isOperable();
-            well->checkWellOperability(ebosSimulator_, this->wellState(), deferred_logger);
-
-            if (!well->isOperable() ) continue;
-
             auto& events = this->wellState().well(well->indexOfWell()).events;
             if (events.hasEvent(WellState::event_mask)) {
                 well->updateWellStateWithTarget(ebosSimulator_, this->groupState(), this->wellState(), deferred_logger);
@@ -1406,36 +1401,16 @@ namespace Opm {
                 // so next time step, the well does not consider to have effective events anymore.
                 events.clearEvent(WellState::event_mask);
             }
-
             // solve the well equation initially to improve the initial solution of the well model
             if (param_.solve_welleq_initially_) {
                 well->solveWellEquation(ebosSimulator_, this->wellState(), this->groupState(), deferred_logger);
             }
-
+        }
-            const bool well_operable = well->isOperable();
-            if (!well_operable && old_well_operable) {
-                const Well& well_ecl = getWellEcl(well->name());
-                if (well_ecl.getAutomaticShutIn()) {
-                    deferred_logger.info(" well " + well->name() + " gets SHUT at the beginning of the time step ");
-                } else {
-                    if (!well->wellIsStopped()) {
-                        deferred_logger.info(" well " + well->name() + " gets STOPPED at the beginning of the time step ");
-                        well->stopWell();
-                    }
-                }
-            } else if (well_operable && !old_well_operable) {
-                deferred_logger.info(" well " + well->name() + " gets REVIVED at the beginning of the time step ");
-                well->openWell();
-            }
-
-        }  // end of for (const auto& well : well_container_)
         updatePrimaryVariables(deferred_logger);
     }
 
 
 
-
-
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -439,7 +439,7 @@ updatePrimaryVariables(const WellState& well_state) const
 
     // TODO: to test using rate conversion coefficients to see if it will be better than
     // this default one
-    if (!baseif_.isOperable() && !baseif_.wellIsStopped()) return;
+    if (!baseif_.isOperableAndSolvable() && !baseif_.wellIsStopped()) return;
 
     const Well& well = baseif_.wellEcl();
 
@@ -510,7 +510,7 @@ void
 MultisegmentWellEval<FluidSystem,Indices,Scalar>::
 recoverSolutionWell(const BVector& x, BVectorWell& xw) const
 {
-    if (!baseif_.isOperable() && !baseif_.wellIsStopped()) return;
+    if (!baseif_.isOperableAndSolvable() && !baseif_.wellIsStopped()) return;
 
     BVectorWell resWell = resWell_;
     // resWell = resWell - B * x

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -183,7 +183,7 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     apply(const BVector& x, BVector& Ax) const
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         if ( this->param_.matrix_add_well_contributions_ )
         {
@@ -210,7 +210,7 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     apply(BVector& r) const
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         // invDrw_ = duneD^-1 * resWell_
         const BVectorWell invDrw = mswellhelpers::applyUMFPack(this->duneD_, this->duneDSolver_, this->resWell_);
@@ -227,7 +227,7 @@ namespace Opm
                                           WellState& well_state,
                                           DeferredLogger& deferred_logger) const
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         BVectorWell xw(1);
         this->recoverSolutionWell(x, xw);
@@ -484,7 +484,7 @@ namespace Opm
     MultisegmentWell<TypeTag>::
     solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger)
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         // We assemble the well equations, then we check the convergence,
         // which is why we do not put the assembleWellEq here.
@@ -578,7 +578,7 @@ namespace Opm
                     DeferredLogger& deferred_logger,
                     const double relaxation_factor) const
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         const double dFLimit = this->param_.dwell_fraction_max_;
         const double max_pressure_change = this->param_.max_pressure_change_ms_wells_;
@@ -1336,7 +1336,7 @@ namespace Opm
                              const GroupState& group_state,
                              DeferredLogger& deferred_logger)
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return true;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return true;
 
         const int max_iter_number = this->param_.max_inner_iter_ms_wells_;
         const WellState well_state0 = well_state;
@@ -1460,7 +1460,7 @@ namespace Opm
                                    DeferredLogger& deferred_logger)
     {
 
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         // update the upwinding segments
         this->updateUpwindingSegments();

opm/simulators/wells/StandardWellEval.cpp

@@ -250,7 +250,7 @@ updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_log
     static constexpr int Oil = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Oil;
     static constexpr int Water = WellInterfaceIndices<FluidSystem,Indices,Scalar>::Water;
 
-    if (!baseif_.isOperable() && !baseif_.wellIsStopped()) return;
+    if (!baseif_.isOperableAndSolvable() && !baseif_.wellIsStopped()) return;
 
     const int well_index = baseif_.indexOfWell();
     const int np = baseif_.numPhases();

opm/simulators/wells/StandardWell_impl.hpp

@@ -397,7 +397,7 @@ namespace Opm
     {
         // TODO: only_wells should be put back to save some computation
         // for example, the matrices B C does not need to update if only_wells
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         // clear all entries
         this->duneB_ = 0.0;
@@ -852,7 +852,7 @@ namespace Opm
                     WellState& well_state,
                     DeferredLogger& deferred_logger) const
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         updatePrimaryVariablesNewton(dwells, well_state);
 
@@ -1517,7 +1517,7 @@ namespace Opm
     StandardWell<TypeTag>::
     solveEqAndUpdateWellState(WellState& well_state, DeferredLogger& deferred_logger)
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         // We assemble the well equations, then we check the convergence,
         // which is why we do not put the assembleWellEq here.
@@ -1552,7 +1552,7 @@ namespace Opm
     StandardWell<TypeTag>::
     apply(const BVector& x, BVector& Ax) const
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         if (this->param_.matrix_add_well_contributions_)
         {
@@ -1583,7 +1583,7 @@ namespace Opm
     StandardWell<TypeTag>::
     apply(BVector& r) const
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         assert( this->invDrw_.size() == this->invDuneD_.N() );
 
@@ -1598,7 +1598,7 @@ namespace Opm
     StandardWell<TypeTag>::
     recoverSolutionWell(const BVector& x, BVectorWell& xw) const
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         BVectorWell resWell = this->resWell_;
         // resWell = resWell - B * x
@@ -1618,7 +1618,7 @@ namespace Opm
                                           WellState& well_state,
                                           DeferredLogger& deferred_logger) const
     {
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         BVectorWell xw(1);
         xw[0].resize(this->numWellEq_);
@@ -1903,7 +1903,7 @@ namespace Opm
     updatePrimaryVariables(const WellState& well_state, DeferredLogger& deferred_logger) const
     {
         this->StdWellEval::updatePrimaryVariables(well_state, deferred_logger);
-        if (!this->isOperable() && !this->wellIsStopped()) return;
+        if (!this->isOperableAndSolvable() && !this->wellIsStopped()) return;
 
         // other primary variables related to polymer injection
         if constexpr (Base::has_polymermw) {

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -352,7 +352,7 @@ void WellInterfaceGeneric::updateWellTestStatePhysical(const double simulation_t
                                                        WellTestState& well_test_state,
                                                        DeferredLogger& deferred_logger) const
 {
-    if (!isOperable()) {
+    if (!isOperableAndSolvable()) {
         if (well_test_state.hasWellClosed(name(), WellTestConfig::Reason::ECONOMIC) ||
             well_test_state.hasWellClosed(name(), WellTestConfig::Reason::PHYSICAL) ) {
             // Already closed, do nothing.
@@ -368,9 +368,9 @@ void WellInterfaceGeneric::updateWellTestStatePhysical(const double simulation_t
     }
 }
 
-bool WellInterfaceGeneric::isOperable() const
+bool WellInterfaceGeneric::isOperableAndSolvable() const
 {
-    return operability_status_.isOperable();
+    return operability_status_.isOperableAndSolvable();
 }
 
 double WellInterfaceGeneric::getALQ(const WellState& well_state) const

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -82,7 +82,7 @@ public:
     bool underPredictionMode() const;
 
     // whether the well is operable
-    bool isOperable() const;
+    bool isOperableAndSolvable() const;
 
     void initCompletions();
     void closeCompletions(WellTestState& wellTestState);
@@ -180,7 +180,7 @@ protected:
 
     // definition of the struct OperabilityStatus
     struct OperabilityStatus {
-        bool isOperable() const {
+        bool isOperableAndSolvable() const {
             if (!operable_under_only_bhp_limit || !solvable) {
                 return false;
             } else {
@@ -196,12 +196,11 @@ protected:
             return can_obtain_bhp_with_thp_limit && obey_bhp_limit_with_thp_limit;
         }
 
-        void reset() {
+        void resetOperability() {
             operable_under_only_bhp_limit = true;
             obey_thp_limit_under_bhp_limit = true;
             can_obtain_bhp_with_thp_limit = true;
             obey_bhp_limit_with_thp_limit = true;
-            solvable = true;
         }
 
         // whether the well can be operated under bhp limit

opm/simulators/wells/WellInterface_impl.hpp

@@ -246,7 +246,7 @@ namespace Opm
             }
 
             updateWellOperability(simulator, well_state_copy, deferred_logger);
-            if ( !this->isOperable() ) {
+            if ( !this->isOperableAndSolvable() ) {
                 const auto msg = fmt::format("WTEST: Well {} is not operable (physical)", this->name());
                 deferred_logger.debug(msg);
                 return;
@@ -346,7 +346,7 @@ namespace Opm
                       const GroupState& group_state,
                       DeferredLogger& deferred_logger)
     {
-        if (!this->isOperable())
+        if (!this->isOperableAndSolvable())
             return;
 
         // keep a copy of the original well state
@@ -357,10 +357,6 @@ namespace Opm
             const int max_iter = param_.max_welleq_iter_;
             deferred_logger.debug("Compute initial well solution for well " + this->name() + ". Failed to converge in "
                                   + std::to_string(max_iter) + " iterations");
-            // the well operability system currently works only for producers in prediction mode
-            if (this->shutUnsolvableWells())
-                this->operability_status_.solvable = false;
-
             well_state = well_state0;
         }
     }
@@ -376,21 +372,25 @@ namespace Opm
                    const GroupState& group_state,
                    DeferredLogger& deferred_logger)
     {
-        const bool old_well_operable = this->operability_status_.isOperable();
+        const bool old_well_operable = this->operability_status_.isOperableAndSolvable();
-        checkWellOperability(ebosSimulator, well_state, deferred_logger);
+
+        if (param_.check_well_operability_iter_)
+            checkWellOperability(ebosSimulator, well_state, deferred_logger);
 
         // only use inner well iterations for the first newton iterations.
         const int iteration_idx = ebosSimulator.model().newtonMethod().numIterations();
-        bool converged = true;
+        if (iteration_idx < param_.max_niter_inner_well_iter_) {
-        if (iteration_idx < param_.max_niter_inner_well_iter_)
+            this->operability_status_.solvable = true;
-            converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
+            bool converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
 
-        // unsolvable wells are treated as not operable and will not be solved for in this iteration.
+            // unsolvable wells are treated as not operable and will not be solved for in this iteration.
-        if (!converged) {
+            if (!converged) {
-            if (this->shutUnsolvableWells())
+                if (this->shutUnsolvableWells())
-                this->operability_status_.solvable = false;
+                    this->operability_status_.solvable = false;
+            }
         }
-        const bool well_operable = this->operability_status_.isOperable();
+
+        const bool well_operable = this->operability_status_.isOperableAndSolvable();
         if (!well_operable && old_well_operable) {
             if (this->well_ecl_.getAutomaticShutIn()) {
                 deferred_logger.info(" well " + this->name() + " gets SHUT during iteration ");
@@ -454,8 +454,7 @@ namespace Opm
                          DeferredLogger& deferred_logger)
     {
 
-        const bool checkOperability = EWOMS_GET_PARAM(TypeTag, bool, EnableWellOperabilityCheck);
+        if (!param_.check_well_operability_) {
-        if (!checkOperability) {
             return;
         }
 
@@ -497,7 +496,7 @@ namespace Opm
                           const WellState& well_state,
                           DeferredLogger& deferred_logger)
     {
-        this->operability_status_.reset();
+        this->operability_status_.resetOperability();
 
         auto current_control = well_state.well(this->index_of_well_).production_cmode;
         // Operability checking is not free