Merge pull request #4149 from totto82/improve_conv_thp

Improve convergence for newly opened wells with thp control

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -287,7 +287,7 @@ namespace Opm
                 bhp_controlled_well = true;
             }
         }
-        if (thp_controlled_well || bhp_controlled_well) {
+        if (!this->changed_to_open_this_step_ && (thp_controlled_well || bhp_controlled_well)) {
 
             double total_rate = 0.0;
             const double sign = this->isInjector() ? 1.0:-1.0;

opm/simulators/wells/SingleWellState.cpp

@@ -156,9 +156,6 @@ double SingleWellState::sum_solvent_rates() const {
 void SingleWellState::update_producer_targets(const Well& ecl_well, const SummaryState& st) {
     const double bhp_safety_factor = 0.99;
     const auto& prod_controls = ecl_well.productionControls(st);
-    if (prod_controls.hasControl(Well::ProducerCMode::THP))
-        this->thp = prod_controls.thp_limit;
-
 
     auto cmode_is_bhp = (prod_controls.cmode == Well::ProducerCMode::BHP);
     auto bhp_limit = prod_controls.bhp_limit;
@@ -172,11 +169,6 @@ void SingleWellState::update_producer_targets(const Well& ecl_well, const Summar
         return;
     }
 
-    if (prod_controls.cmode == Well::ProducerCMode::GRUP) {
-        this->bhp = this->perf_data.pressure_first_connection * bhp_safety_factor;
-        return;
-    }
-
     switch (prod_controls.cmode) {
     case Well::ProducerCMode::ORAT:
         assert(this->pu.phase_used[BlackoilPhases::Liquid]);
@@ -190,11 +182,29 @@ void SingleWellState::update_producer_targets(const Well& ecl_well, const Summar
         assert(this->pu.phase_used[BlackoilPhases::Vapour]);
         this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate;
         break;
+    case Well::ProducerCMode::GRUP:
+    case Well::ProducerCMode::THP:
+    case Well::ProducerCMode::BHP:
+        if (this->pu.phase_used[BlackoilPhases::Liquid]) {
+            this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = -1000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
+        }
+        if (this->pu.phase_used[BlackoilPhases::Aqua]) {
+            this->surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = -1000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
+        }
+        if (this->pu.phase_used[BlackoilPhases::Vapour]){
+            this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = -100000.0 * Opm::unit::cubic(Opm::unit::meter) / Opm::unit::day;
+        }
+        break;
+
     default:
         // Keep zero init.
         break;
     }
 
+    if (prod_controls.cmode == Well::ProducerCMode::THP) {
+        this->thp = prod_controls.thp_limit;
+    }
+
     if (cmode_is_bhp)
         this->bhp = bhp_limit;
     else

opm/simulators/wells/StandardWell_impl.hpp

@@ -1994,7 +1994,7 @@ namespace Opm
                 bhp_controlled_well = true;
             }
         }
-        if (thp_controlled_well || bhp_controlled_well) {
+        if (!this->changed_to_open_this_step_ && (thp_controlled_well || bhp_controlled_well)) {
 
             double total_rate = 0.0;
             const double sign = this->isInjector() ? 1.0:-1.0;

opm/simulators/wells/WellInterface_impl.hpp

@@ -459,7 +459,39 @@ namespace Opm
         // keep a copy of the original well state
         const WellState well_state0 = well_state;
         const double dt = ebosSimulator.timeStepSize();
-        const bool converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
+        bool converged = iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
+
+        // Newly opened wells with THP control sometimes struggles to
+        // converge due to bad initial guess. Or due to the simple fact
+        // that the well needs to change to another control.
+        // We therefore try to solve the well with BHP control to get
+        // an better initial guess.
+        // If the well is supposed to operate under THP control
+        // "updateWellControl" will switch it back to THP later.
+        if (!converged) {
+            auto& ws = well_state.well(this->indexOfWell());
+            bool thp_control = false;
+            if (this->well_ecl_.isInjector()) {
+                thp_control = ws.injection_cmode == Well::InjectorCMode::THP;
+                if (thp_control) {
+                    ws.injection_cmode = Well::InjectorCMode::BHP;
+                    this->well_control_log_.push_back(Well::InjectorCMode2String(Well::InjectorCMode::THP));
+                }
+            } else {
+                thp_control = ws.production_cmode == Well::ProducerCMode::THP;
+                if (thp_control) {
+                    ws.production_cmode = Well::ProducerCMode::BHP;
+                    this->well_control_log_.push_back(Well::ProducerCMode2String(Well::ProducerCMode::THP));
+                }
+            }
+            if (thp_control) {
+                const std::string msg = std::string("The newly opened well ") + this->name()
+                    + std::string(" with THP control did not converge during inner iterations, we try again with bhp control");
+                deferred_logger.debug(msg);
+                converged = this->iterateWellEquations(ebosSimulator, dt, well_state, group_state, deferred_logger);
+            }
+        }
+
         if (!converged) {
             const int max_iter = param_.max_welleq_iter_;
             deferred_logger.debug("Compute initial well solution for well " + this->name() + ". Failed to converge in "