Merge pull request #4094 from totto82/checkControl

Check control

opm/simulators/flow/BlackoilModelEbos.hpp

@@ -982,7 +982,7 @@ namespace Opm {
             // Get convergence reports for reservoir and wells.
             std::vector<Scalar> B_avg(numEq, 0.0);
             auto report = getReservoirConvergence(timer.currentStepLength(), iteration, B_avg, residual_norms);
-            report += wellModel().getWellConvergence(B_avg);
+            report += wellModel().getWellConvergence(B_avg, /*checkWellGroupControls*/report.converged());
 
             return report;
         }

opm/simulators/timestepping/ConvergenceReport.hpp

@@ -86,7 +86,7 @@ namespace Opm
             : status_{AllGood}
             , res_failures_{}
             , well_failures_{}
-            , groupConverged_(true)
+            , wellGroupTargetsViolated_(false)
         {
         }
 
@@ -95,7 +95,7 @@ namespace Opm
             status_ = AllGood;
             res_failures_.clear();
             well_failures_.clear();
-            groupConverged_ = true;
+            wellGroupTargetsViolated_ = false;
         }
 
         void setReservoirFailed(const ReservoirFailure& rf)
@@ -110,9 +110,9 @@ namespace Opm
             well_failures_.push_back(wf);
         }
 
-        void setGroupConverged(const bool groupConverged)
+        void setWellGroupTargetsViolated(const bool wellGroupTargetsViolated)
         {
-            groupConverged_ = groupConverged;
+            wellGroupTargetsViolated_ = wellGroupTargetsViolated;
         }
 
         ConvergenceReport& operator+=(const ConvergenceReport& other)
@@ -122,6 +122,7 @@ namespace Opm
             well_failures_.insert(well_failures_.end(), other.well_failures_.begin(), other.well_failures_.end());
             assert(reservoirFailed() != res_failures_.empty());
             assert(wellFailed() != well_failures_.empty());
+            wellGroupTargetsViolated_ = (wellGroupTargetsViolated_ || other.wellGroupTargetsViolated_);
             return *this;
         }
 
@@ -129,7 +130,7 @@ namespace Opm
 
         bool converged() const
         {
-            return status_ == AllGood && groupConverged_;
+            return (status_ == AllGood) && !wellGroupTargetsViolated_;
         }
 
         bool reservoirFailed() const
@@ -174,7 +175,7 @@ namespace Opm
         Status status_;
         std::vector<ReservoirFailure> res_failures_;
         std::vector<WellFailure> well_failures_;
-        bool groupConverged_;
+        bool wellGroupTargetsViolated_;
     };
 
 } // namespace Opm

opm/simulators/timestepping/SimulatorReport.cpp

@@ -50,6 +50,7 @@ namespace Opm
           min_linear_iterations ( std::numeric_limits<unsigned int>::max() ),
           max_linear_iterations ( 0 ),
           converged(false),
+          well_group_control_changed(false),
           exit_status(EXIT_SUCCESS),
           global_time(0),
           timestep_length(0.0)

opm/simulators/timestepping/SimulatorReport.hpp

@@ -50,6 +50,7 @@ namespace Opm
 
 
         bool converged;
+        bool well_group_control_changed;
         int exit_status;
 
         double global_time;

opm/simulators/wells/BlackoilWellModel.hpp

@@ -260,7 +260,7 @@ namespace Opm {
             void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const;
 
             // Check if well equations is converged.
-            ConvergenceReport getWellConvergence(const std::vector<Scalar>& B_avg, const bool checkGroupConvergence = false) const;
+            ConvergenceReport getWellConvergence(const std::vector<Scalar>& B_avg, const bool checkWellGroupControls = false) const;
 
             const SimulatorReportSingle& lastReport() const;
 
@@ -281,7 +281,7 @@ namespace Opm {
             void prepareTimeStep(DeferredLogger& deferred_logger);
             void initPrimaryVariablesEvaluation() const;
             bool shouldBalanceNetwork(const int reportStepIndex, const int iterationIdx) const;
-            std::pair<bool, double> updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls);
+            std::tuple<bool, bool, double> updateWellControls(DeferredLogger& deferred_logger);
 
             void updateAndCommunicate(const int reportStepIdx,
                                       const int iterationIdx,
@@ -373,7 +373,7 @@ namespace Opm {
             // and in the well equations.
             void assemble(const int iterationIdx,
                           const double dt);
-            void assembleImpl(const int iterationIdx,
+            bool assembleImpl(const int iterationIdx,
                               const double dt,
                               const std::size_t recursion_level,
                               DeferredLogger& local_deferredLogger);

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -880,9 +880,11 @@ namespace Opm {
                                                terminal_output_, grid().comm());
         }
 
-        assembleImpl(iterationIdx, dt, 0, local_deferredLogger);
+        const bool well_group_control_changed = assembleImpl(iterationIdx, dt, 0, local_deferredLogger);
 
-        last_report_.converged = true;
+        // if group or well control changes we don't consider the
+        // case converged
+        last_report_.well_group_control_changed = well_group_control_changed;
         last_report_.assemble_time_well += perfTimer.stop();
     }
 
@@ -891,7 +893,7 @@ namespace Opm {
 
 
     template<typename TypeTag>
-    void
+    bool
     BlackoilWellModel<TypeTag>::
     assembleImpl(const int iterationIdx,
                  const double dt,
@@ -899,7 +901,7 @@ namespace Opm {
                  DeferredLogger& local_deferredLogger)
     {
 
-        const auto [network_changed, network_imbalance] = updateWellControls(local_deferredLogger, /* check group controls */ true);
+        auto [well_group_control_changed, network_changed, network_imbalance] = updateWellControls(local_deferredLogger);
 
         bool alq_updated = false;
         OPM_BEGIN_PARALLEL_TRY_CATCH();
@@ -931,10 +933,11 @@ namespace Opm {
                 const auto& balance = schedule()[reportStepIdx].network_balance();
                 // Iterate if not converged, and number of iterations is not yet max (NETBALAN item 3).
                 if (recursion_level < balance.pressure_max_iter() && network_imbalance > balance.pressure_tolerance()) {
-                    assembleImpl(iterationIdx, dt, recursion_level + 1, local_deferredLogger);
+                    well_group_control_changed = assembleImpl(iterationIdx, dt, recursion_level + 1, local_deferredLogger);
                 }
             }
         }
+        return well_group_control_changed;
     }
 
 
@@ -1389,7 +1392,7 @@ namespace Opm {
     template<typename TypeTag>
     ConvergenceReport
     BlackoilWellModel<TypeTag>::
-    getWellConvergence(const std::vector<Scalar>& B_avg, bool checkGroupConvergence) const
+    getWellConvergence(const std::vector<Scalar>& B_avg, bool checkWellGroupControls) const
     {
 
         DeferredLogger local_deferredLogger;
@@ -1406,15 +1409,19 @@ namespace Opm {
                 local_report += report;
             }
         }
-        
+
         const Opm::Parallel::Communication comm = grid().comm();
         DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger, comm);
+        ConvergenceReport report = gatherConvergenceReport(local_report, comm);
+
+        // the well_group_control_changed info is already communicated
+        if (checkWellGroupControls) {
+            report.setWellGroupTargetsViolated(this->lastReport().well_group_control_changed);
+        }
+
         if (terminal_output_) {
             global_deferredLogger.logMessages();
         }
-
-        ConvergenceReport report = gatherConvergenceReport(local_report, comm);
-
         // Log debug messages for NaN or too large residuals.
         if (terminal_output_) {
             for (const auto& f : report.wellFailures()) {
@@ -1425,15 +1432,6 @@ namespace Opm {
                 }
             }
         }
-
-        if (checkGroupConvergence) {
-            const int reportStepIdx = ebosSimulator_.episodeIndex();
-            const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
-            bool violated = checkGroupConstraints(fieldGroup,
-                                                  ebosSimulator_.episodeIndex(),
-                                                  global_deferredLogger);
-            report.setGroupConverged(!violated);
-        }
         return report;
     }
 
@@ -1481,14 +1479,14 @@ namespace Opm {
 
 
     template<typename TypeTag>
-    std::pair<bool, double>
+    std::tuple<bool, bool, double>
     BlackoilWellModel<TypeTag>::
-    updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls)
+    updateWellControls(DeferredLogger& deferred_logger)
     {
         // Even if there are no wells active locally, we cannot
         // return as the DeferredLogger uses global communication.
         // For no well active globally we simply return.
-        if( !wellsActive() ) return { false, 0.0 };
+        if( !wellsActive() ) return { false, false, 0.0 };
 
         const int episodeIdx = ebosSimulator_.episodeIndex();
         const int iterationIdx = ebosSimulator_.model().newtonMethod().numIterations();
@@ -1501,39 +1499,34 @@ namespace Opm {
         const bool network_changed = comm.sum(local_network_changed);
         const double network_imbalance = comm.max(local_network_imbalance);
 
-        std::set<std::string> switched_wells;
+        bool changed_well_group = false;
-
+        // Check group individual constraints.
-        if (checkGroupControls) {
+        const int nupcol = schedule()[episodeIdx].nupcol();
-
+        // don't switch group control when iterationIdx > nupcol
-            // Check group individual constraints.
+        // to avoid oscilations between group controls
-            const int nupcol = schedule()[episodeIdx].nupcol();
+        if (iterationIdx <= nupcol) {
-            // don't switch group control when iterationIdx > nupcol
+            const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
-            // to avoid oscilations between group controls
+            changed_well_group = updateGroupControls(fieldGroup, deferred_logger, episodeIdx, iterationIdx);
-            if (iterationIdx <= nupcol) {
+        }
-                const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
+        // Check wells' group constraints and communicate.
-                updateGroupControls(fieldGroup, deferred_logger, episodeIdx, iterationIdx);
+        bool changed_well_to_group = false;
+        for (const auto& well : well_container_) {
+            const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Group;
+            const bool changed_well = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
+            if (changed_well) {
+                changed_well_to_group = changed_well || changed_well_to_group;
             }
-            // Check wells' group constraints and communicate.
+        }
-            bool changed_well_group = false;
+
-            for (const auto& well : well_container_) {
+        changed_well_to_group = comm.sum(changed_well_to_group);
-                const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Group;
+        if (changed_well_to_group) {
-                const bool changed_well = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
+            updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
-                if (changed_well) {
+            changed_well_group = true;
-                    switched_wells.insert(well->name());
-                    changed_well_group = changed_well || changed_well_group;
-                }
-            }
-            changed_well_group = comm.sum(changed_well_group);
-            if (changed_well_group)
-                updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
         }
 
         // Check individual well constraints and communicate.
         bool changed_well_individual = false;
         for (const auto& well : well_container_) {
-            if (switched_wells.count(well->name())) {
-                continue;
-            }
             const auto mode = WellInterface<TypeTag>::IndividualOrGroup::Individual;
             const bool changed_well = well->updateWellControl(ebosSimulator_, mode, this->wellState(), this->groupState(), deferred_logger);
             if (changed_well) {
@@ -1541,14 +1534,16 @@ namespace Opm {
             }
         }
         changed_well_individual = comm.sum(changed_well_individual);
-        if (changed_well_individual)
+        if (changed_well_individual) {
             updateAndCommunicate(episodeIdx, iterationIdx, deferred_logger);
+            changed_well_group = true;
+        }
 
         // update wsolvent fraction for REIN wells
         const Group& fieldGroup = schedule().getGroup("FIELD", episodeIdx);
         updateWsolvent(fieldGroup, episodeIdx,  this->nupcolWellState());
 
-        return { network_changed, network_imbalance };
+        return { changed_well_group, network_changed, network_imbalance };
     }
 
 

opm/simulators/wells/MultisegmentWellEval.cpp

@@ -1363,6 +1363,12 @@ updateThp(WellState& well_state,
         ws.thp = 0;
         return;
     }
+    // For THP controlled wells, we know the thp value
+    bool thp_controlled = baseif_.isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP:
+                                              ws.production_cmode == Well::ProducerCMode::THP;
+    if (thp_controlled) {
+        return;
+    }
 
     // the well is under other control types, we calculate the thp based on bhp and rates
     std::vector<double> rates(3, 0.0);

opm/simulators/wells/StandardWellEval.cpp

@@ -573,6 +573,13 @@ updateThp(WellState& well_state,
         return;
     }
 
+    // For THP controlled wells, we know the thp value
+    bool thp_controlled = baseif_.isInjector() ? ws.injection_cmode == Well::InjectorCMode::THP:
+                                              ws.production_cmode == Well::ProducerCMode::THP;
+    if (thp_controlled) {
+        return;
+    }
+
     // the well is under other control types, we calculate the thp based on bhp and rates
     std::vector<double> rates(3, 0.0);
 

opm/simulators/wells/WellInterface_impl.hpp

@@ -710,15 +710,15 @@ namespace Opm
     {
         this->operability_status_.resetOperability();
 
-        bool thp_controled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
+        bool thp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::THP:
                                               well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::THP;
-        bool bhp_controled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
+        bool bhp_controlled = this->isInjector() ? well_state.well(this->index_of_well_).injection_cmode == Well::InjectorCMode::BHP:
                                               well_state.well(this->index_of_well_).production_cmode == Well::ProducerCMode::BHP;
 
         // Operability checking is not free
         // Only check wells under BHP and THP control
-        bool check_thp = thp_controled || this->operability_status_.thp_limit_violated_but_not_switched;
+        bool check_thp = thp_controlled || this->operability_status_.thp_limit_violated_but_not_switched;
-        if (check_thp || bhp_controled) {
+        if (check_thp || bhp_controlled) {
             updateIPR(ebos_simulator, deferred_logger);
             checkOperabilityUnderBHPLimit(well_state, ebos_simulator, deferred_logger);
         }
@@ -813,6 +813,7 @@ namespace Opm
                 auto rates = ws.surface_rates;
                 double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
                 ws.bhp = bhp;
+                ws.thp = this->getTHPConstraint(summaryState);
 
                 // if the total rates are negative or zero
                 // we try to provide a better intial well rate
@@ -1037,6 +1038,7 @@ namespace Opm
                 this->adaptRatesForVFP(rates);
                 double bhp = this->calculateBhpFromThp(well_state, rates, well, summaryState, this->getRefDensity(), deferred_logger);
                 ws.bhp = bhp;
+                ws.thp = this->getTHPConstraint(summaryState);
 
                 // if the total rates are negative or zero
                 // we try to provide a better intial well rate

tests/test_glift1.cpp

@@ -146,7 +146,7 @@ BOOST_AUTO_TEST_CASE(G1)
     Opm::DeferredLogger deferred_logger;
     well_model.calculateExplicitQuantities(deferred_logger);
     well_model.prepareTimeStep(deferred_logger);
-    well_model.updateWellControls(deferred_logger, /* check group controls */ true);
+    well_model.updateWellControls(deferred_logger);
     well_model.initPrimaryVariablesEvaluation();
     Opm::WellInterface<TypeTag> *well_ptr = well_model.getWell("B-1H").get();
     StdWell *std_well = dynamic_cast<StdWell *>(well_ptr);