Merge pull request #1657 from GitPaean/welltesting_physical_reason

first version of the well testing for physical reason

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -324,8 +324,6 @@ namespace Opm {
 
         for (const auto& testWell : wellsForTesting) {
             const std::string& well_name = testWell.first;
-            const std::string msg = std::string("well ") + well_name + std::string(" is tested");
-            OpmLog::info(msg);
 
             // this is the well we will test
             WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx);

opm/autodiff/MultisegmentWell.hpp

@@ -117,8 +117,8 @@ namespace Opm
                                     WellState& well_state) override;
 
         /// updating the well state based the current control mode
-        virtual void updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
-                                               WellState& well_state) /* const */ override;
+        virtual void updateWellStateWithTarget(const Simulator& ebos_simulator,
+                                               WellState& well_state) const override;
 
         /// check whether the well equations get converged for this well
         virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const override;
@@ -349,6 +349,11 @@ namespace Opm
         void assembleWellEqWithoutIteration(const Simulator& ebosSimulator,
                                             const double dt,
                                             WellState& well_state);
+
+        virtual void wellTestingPhysical(Simulator& simulator, const std::vector<double>& B_avg,
+                                         const double simulation_time, const int report_step,
+                                         const bool terminal_output,
+                                         WellState& well_state, WellTestState& welltest_state) override;
     };
 
 }

opm/autodiff/MultisegmentWell_impl.hpp

@@ -252,8 +252,8 @@ namespace Opm
     template <typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
-                              WellState& well_state) /* const */
+    updateWellStateWithTarget(const Simulator& ebos_simulator,
+                              WellState& well_state) const
     {
         // Updating well state bas on well control
         // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
@@ -1648,7 +1648,7 @@ namespace Opm
     checkWellOperability(const Simulator& /* ebos_simulator */,
                          const WellState& /* well_state */)
     {
-        const std::string msg = "Support of well operatability checking for mutlisegment wells is not implemented "
+        const std::string msg = "Support of well operability checking for multisegment wells is not implemented "
                                 "yet, checkWellOperability() for " + name() + " will do nothing";
         OpmLog::warning("NO_OPERATABILITY_CHECKING_MS_WELLS", msg);
     }
@@ -1888,4 +1888,20 @@ namespace Opm
         }
     }
 
+
+
+
+
+    template<typename TypeTag>
+    void
+    MultisegmentWell<TypeTag>::
+    wellTestingPhysical(Simulator& simulator, const std::vector<double>& B_avg,
+                        const double simulation_time, const int report_step, const bool terminal_output,
+                        WellState& well_state, WellTestState& welltest_state)
+    {
+        const std::string msg = "Support of well testing for physical limits for multisegment wells is not "
+                                "implemented yet, wellTestingPhysical() for " + name() + " will do nothing";
+        OpmLog::warning("NO_WELLTESTPHYSICAL_CHECKING_MS_WELLS", msg);
+    }
+
 }

opm/autodiff/StandardWell.hpp

@@ -139,8 +139,8 @@ namespace Opm
                                     const double dt,
                                     WellState& well_state) override;
 
-        virtual void updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
-                                               WellState& well_state) /* const */ override;
+        virtual void updateWellStateWithTarget(const Simulator& ebos_simulator,
+                                               WellState& well_state) const override;
 
         /// check whether the well equations get converged for this well
         virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const override;
@@ -348,11 +348,16 @@ namespace Opm
         // updating the inflow based on the current reservoir condition
         void updateIPR(const Simulator& ebos_simulator) const;
 
-        // check whether the well is operable under the current reservoir condition
+        // update the operability status of the well is operable under the current reservoir condition
         // mostly related to BHP limit and THP limit
         virtual void checkWellOperability(const Simulator& ebos_simulator,
                                           const WellState& well_state) override;
 
+        // check whether the well is operable under the current reservoir condition
+        // mostly related to BHP limit and THP limit
+        void updateWellOperability(const Simulator& ebos_simulator,
+                                   const WellState& well_state);
+
         // check whether the well is operable under BHP limit with current reservoir condition
         void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator);
 
@@ -398,6 +403,10 @@ namespace Opm
         // calculate a relaxation factor to avoid overshoot of total rates
         static double relaxationFactorRate(const std::vector<double>& primary_variables,
                                            const BVectorWell& dwells);
+
+        virtual void wellTestingPhysical(Simulator& simulator, const std::vector<double>& B_avg,
+                                         const double simulation_time, const int report_step, const bool terminal_output,
+                                         WellState& well_state, WellTestState& welltest_state);
     };
 
 }

opm/autodiff/StandardWell_impl.hpp

@@ -1088,8 +1088,8 @@ namespace Opm
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
-    updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
-                              WellState& well_state) /* const */
+    updateWellStateWithTarget(const Simulator& ebos_simulator,
+                              WellState& well_state) const
     {
         // number of phases
         const int np = number_of_phases_;
@@ -1112,8 +1112,6 @@ namespace Opm
             break;
 
         case THP: {
-            assert(this->isOperable() );
-
             // when a well can not work under THP target, it switches to BHP control
             if (this->operability_status_.isOperableUnderTHPLimit() ) {
                 updateWellStateWithTHPTargetIPR(ebos_simulator, well_state);
@@ -1307,6 +1305,27 @@ namespace Opm
 
         const bool old_well_operable = this->operability_status_.isOperable();
 
+        updateWellOperability(ebos_simulator, well_state);
+
+        const bool well_operable = this->operability_status_.isOperable();
+
+        if (!well_operable && old_well_operable) {
+            OpmLog::info(" well " + name() + " gets SHUT during iteration ");
+        } else if (well_operable && !old_well_operable) {
+            OpmLog::info(" well " + name() + " gets REVIVED during iteration ");
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    updateWellOperability(const Simulator& ebos_simulator,
+                          const WellState& well_state)
+    {
         this->operability_status_.reset();
 
         updateIPR(ebos_simulator);
@@ -1316,18 +1335,7 @@ namespace Opm
 
         // checking whether the well can operate under the THP constraints.
         if (this->wellHasTHPConstraints()) {
-            this->operability_status_.has_thp_constaint = true;
             checkOperabilityUnderTHPLimitProducer(ebos_simulator);
-            this->operability_status_.can_produce_inject_with_current_bhp =
-                canProduceInjectWithCurrentBhp(ebos_simulator, well_state);
-        }
-
-        const bool well_operable = this->operability_status_.isOperable();
-
-        if (!well_operable && old_well_operable) {
-            OpmLog::info(" well " + name() + " gets SHUT during iteration ");
-        } else if (well_operable && !old_well_operable) {
-            OpmLog::info(" well " + name() + " gets REVIVED during iteration ");
         }
     }
 
@@ -2706,4 +2714,67 @@ namespace Opm
 
         return relaxation_factor;
     }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    wellTestingPhysical(Simulator& ebos_simulator, const std::vector<double>& B_avg,
+                        const double simulation_time, const int report_step, const bool terminal_output,
+                        WellState& well_state, WellTestState& welltest_state)
+    {
+        OpmLog::debug(" well " + name() + " is being tested for physical limits");
+
+        // some most difficult things are the explicit quantities, since there is no information
+        // in the WellState to do a decent initialization
+
+        // TODO: Let us assume that the simulator is updated
+
+        // Let us try to do a normal simualtion running, to keep checking the operability status
+        // If the well is not operable during any of the time. It means it does not pass the physical
+        // limit test.
+
+        // create a copy of the well_state to use. If the operability checking is sucessful, we use this one
+        // to replace the original one
+        WellState well_state_copy = well_state;
+
+        // TODO: well state for this well is kind of all zero status
+        // we should be able to provide a better initialization
+        calculateExplicitQuantities(ebos_simulator, well_state_copy);
+
+        updateWellOperability(ebos_simulator, well_state_copy);
+
+        if ( !this->isOperable() ) {
+            const std::string msg = " well " + name() + " is not operable during well testing for physical reason";
+            OpmLog::debug(msg);
+            return;
+        }
+
+        updateWellStateWithTarget(ebos_simulator, well_state_copy);
+
+        calculateExplicitQuantities(ebos_simulator, well_state_copy);
+        updatePrimaryVariables(well_state_copy);
+        initPrimaryVariablesEvaluation();
+
+        const bool converged = this->solveWellEqUntilConverged(ebos_simulator, B_avg, well_state_copy);
+
+        if (!converged) {
+            const std::string msg = " well " + name() + " did not get converged during well testing for physical reason";
+            OpmLog::debug(msg);
+            return;
+        }
+
+        if (this->isOperable() ) {
+            welltest_state.openWell(name() );
+            const std::string msg = " well " + name() + " is re-opened through well testing for physical reason";
+            OpmLog::info(msg);
+            well_state = well_state_copy;
+        } else {
+            const std::string msg = " well " + name() + " is not operable during well testing for physical reason";
+            OpmLog::debug(msg);
+        }
+    }
 }

opm/autodiff/WellInterface.hpp

@@ -172,8 +172,8 @@ namespace Opm
                                            const WellState& well_state,
                                            std::vector<double>& well_potentials) = 0;
 
-        virtual void updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
-                                               WellState& well_state) /* const */ = 0;
+        virtual void updateWellStateWithTarget(const Simulator& ebos_simulator,
+                                               WellState& well_state) const = 0;
 
         void updateWellControl(/* const */ Simulator& ebos_simulator,
                                WellState& well_state,
@@ -218,8 +218,8 @@ namespace Opm
         // Simulator is not const is because that assembleWellEq is non-const Simulator
         void wellTesting(Simulator& simulator, const std::vector<double>& B_avg,
                          const double simulation_time, const int report_step,  const bool terminal_output,
-                         const WellTestConfig::Reason testing_reason, const WellState& well_state,
-                         WellTestState& welltest_state);
+                         const WellTestConfig::Reason testing_reason,
+                         /* const */ WellState& well_state, WellTestState& welltest_state);
 
         void updatePerforatedCell(std::vector<bool>& is_cell_perforated);
 
@@ -364,6 +364,10 @@ namespace Opm
                                  const double simulation_time, const int report_step, const bool terminal_output,
                                  const WellState& well_state, WellTestState& welltest_state);
 
+        virtual void wellTestingPhysical(Simulator& simulator, const std::vector<double>& B_avg,
+                                 const double simulation_time, const int report_step, const bool terminal_output,
+                                 WellState& well_state, WellTestState& welltest_state) = 0;
+
         void updateWellTestStateEconomic(const WellState& well_state,
                                          const double simulation_time,
                                          const bool write_message_to_opmlog,
@@ -374,7 +378,12 @@ namespace Opm
                                          const bool write_message_to_opmlog,
                                          WellTestState& well_test_state) const;
 
-        void solveWellForTesting(Simulator& ebosSimulator, WellState& well_state, const std::vector<double>& B_avg, bool terminal_output);
+        void  solveWellForTesting(Simulator& ebosSimulator, WellState& well_state,
+                                  const std::vector<double>& B_avg, bool terminal_output);
+
+        bool solveWellEqUntilConverged(Simulator& ebosSimulator,
+                                        const std::vector<double>& B_avg,
+                                        WellState& well_state);
 
         void scaleProductivityIndex(const int perfIdx, double& productivity_index) const;
 
@@ -393,8 +402,7 @@ namespace Opm
             if (!operable_under_only_bhp_limit) {
                 return false;
             } else {
-                return ( (isOperableUnderBHPLimit() || isOperableUnderTHPLimit()) &&
-                        !(has_thp_constaint && !can_produce_inject_with_current_bhp) );
+                return ( (isOperableUnderBHPLimit() || isOperableUnderTHPLimit()) );
             }
         }
 
@@ -411,8 +419,6 @@ namespace Opm
             obey_thp_limit_under_bhp_limit = true;
             can_obtain_bhp_with_thp_limit = true;
             obey_bhp_limit_with_thp_limit = true;
-            can_produce_inject_with_current_bhp = true;
-            has_thp_constaint = false;
         }
 
         // whether the well can be operated under bhp limit
@@ -426,15 +432,6 @@ namespace Opm
         bool can_obtain_bhp_with_thp_limit = true;
         // whether the well obey bhp limit when operated under thp limit
         bool obey_bhp_limit_with_thp_limit = true;
-
-        // TODO: the following criterion is based on the current state of
-        // the well, we consider it is a numerical criterion.
-        // at the moment, we only apply it with well has THP constraint.
-        // whether the well can produce / inject with the current bhp of the well
-        // it might be updated with other criterion with investigation with more cases.
-        bool can_produce_inject_with_current_bhp = true;
-        // whether the well has a THP constraint
-        bool has_thp_constaint = false;
     };
 
 }

opm/autodiff/WellInterface_impl.hpp

@@ -923,9 +923,15 @@ namespace Opm
     WellInterface<TypeTag>::
     wellTesting(Simulator& simulator, const std::vector<double>& B_avg,
                 const double simulation_time, const int report_step, const bool terminal_output,
-                const WellTestConfig::Reason testing_reason, const WellState& well_state,
+                const WellTestConfig::Reason testing_reason,
+                /* const */ WellState& well_state,
                 WellTestState& well_test_state)
     {
+        if (testing_reason == WellTestConfig::Reason::PHYSICAL) {
+            wellTestingPhysical(simulator, B_avg, simulation_time, report_step,
+                                terminal_output, well_state, well_test_state);
+        }
+
         if (testing_reason == WellTestConfig::Reason::ECONOMIC) {
             wellTestingEconomic(simulator, B_avg, simulation_time, report_step,
                                 terminal_output, well_state, well_test_state);
@@ -943,6 +949,8 @@ namespace Opm
                         const double simulation_time, const int report_step, const bool terminal_output,
                         const WellState& well_state, WellTestState& welltest_state)
     {
+        OpmLog::debug(" well " + name() + " is being tested for economic limits");
+
         WellState well_state_copy = well_state;
 
         updatePrimaryVariables(well_state_copy);
@@ -1139,6 +1147,42 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    bool
+    WellInterface<TypeTag>::
+    solveWellEqUntilConverged(Simulator& ebosSimulator,
+                              const std::vector<double>& B_avg,
+                              WellState& well_state)
+    {
+        const int max_iter = param_.max_welleq_iter_;
+        int it = 0;
+        const double dt = 1.0; //not used for the well tests
+        bool converged;
+        WellState well_state0 = well_state;
+        do {
+            assembleWellEq(ebosSimulator, dt, well_state);
+
+            auto report = getWellConvergence(B_avg);
+            converged = report.converged();
+            if (converged) {
+                break;
+            }
+
+            ++it;
+            solveEqAndUpdateWellState(well_state);
+
+            wellhelpers::WellSwitchingLogger logger;
+            updateWellControl(ebosSimulator, well_state, logger);
+            initPrimaryVariablesEvaluation();
+        } while (it < max_iter);
+
+        return converged;
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::calculateReservoirRates(WellState& well_state) const
@@ -1176,37 +1220,24 @@ namespace Opm
 
     template<typename TypeTag>
     void
-    WellInterface<TypeTag>::solveWellForTesting(Simulator& ebosSimulator, WellState& well_state, const std::vector<double>& B_avg, bool terminal_output)
+    WellInterface<TypeTag>::
+    solveWellForTesting(Simulator& ebosSimulator, WellState& well_state,
+                        const std::vector<double>& B_avg, bool terminal_output)
     {
-        const int max_iter = param_.max_welleq_iter_;
-        int it = 0;
-        const double dt = 1.0; //not used for the well tests
-        bool converged;
-        WellState well_state0 = well_state;
-        do {
-            assembleWellEq(ebosSimulator, dt, well_state);
+        // keep a copy of the original well state
+        const WellState well_state0 = well_state;
 
-            auto report = getWellConvergence(B_avg);
-            converged = report.converged();
-            if (converged) {
-                break;
-            }
-
-            ++it;
-            solveEqAndUpdateWellState(well_state);
-
-            wellhelpers::WellSwitchingLogger logger;
-            updateWellControl(ebosSimulator, well_state, logger);
-            initPrimaryVariablesEvaluation();
-        } while (it < max_iter);
+        const bool converged = solveWellEqUntilConverged(ebosSimulator, B_avg, well_state);
 
         if (converged) {
             if ( terminal_output ) {
-                OpmLog::debug("WellTest: Well equation for well " + name() +  " solution gets converged with " + std::to_string(it) + " iterations");
+                OpmLog::debug("WellTest: Well equation for well " + name() +  " solution gets converged");
             }
         } else {
             if ( terminal_output ) {
-                OpmLog::debug("WellTest: Well equation for well" +name() + " solution failed in getting converged with " + std::to_string(it) + " iterations");
+                const int max_iter = param_.max_welleq_iter_;
+                OpmLog::debug("WellTest: Well equation for well" +name() + " solution failed in getting converged with "
+                              + std::to_string(max_iter) + " iterations");
             }
             well_state = well_state0;
         }