Merge pull request #3705 from totto82/fixSignPot

guard against negative and trivial potentials

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -2374,4 +2374,19 @@ runWellPIScaling(const int timeStepIdx,
     this->last_run_wellpi_ = timeStepIdx;
 }
 
+
+bool
+BlackoilWellModelGeneric::
+guideRateUpdateIsNeeded() const {
+    const auto need_update =
+    std::any_of(this->well_container_generic_.begin(),
+                this->well_container_generic_.end(),
+    [](const WellInterfaceGeneric* well)
+    {
+        return well->changedToOpenThisStep();
+    });
+    return this->comm_.max(static_cast<int>(need_update));
+}
+
+
 }

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -396,6 +396,8 @@ protected:
                               const SummaryConfig& summaryConfig,
                               DeferredLogger& deferred_logger);
 
+    bool guideRateUpdateIsNeeded() const;
+
     // create the well container
     virtual void createWellContainer(const int time_step) = 0;
     virtual void initWellContainer() = 0;

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -840,6 +840,18 @@ namespace Opm {
         }
         OPM_END_PARALLEL_TRY_CATCH_LOG(local_deferredLogger, "assemble() failed: ",
                                        terminal_output_, grid().comm());
+
+        //update guide rates
+        if (guideRateUpdateIsNeeded()) {
+            const int reportStepIdx = ebosSimulator_.episodeIndex();
+            const double simulationTime = ebosSimulator_.time();
+            const auto& comm = ebosSimulator_.vanguard().grid().comm();
+            const auto& summaryState = ebosSimulator_.vanguard().summaryState();
+            std::vector<double> pot(numPhases(), 0.0);
+            const Group& fieldGroup = schedule().getGroup("FIELD", reportStepIdx);
+            WellGroupHelpers::updateGuideRates(fieldGroup, schedule(), summaryState, this->phase_usage_, reportStepIdx, simulationTime,
+                                        this->wellState(), this->groupState(), comm, &this->guideRate_, pot, local_deferredLogger);
+        }
         last_report_.converged = true;
         last_report_.assemble_time_well += perfTimer.stop();
     }
@@ -1343,7 +1355,8 @@ namespace Opm {
         // and updated only if sucessfull. i.e. the potentials are zero for exceptions
         auto& ws = this->wellState().well(well->indexOfWell());
         for (int p = 0; p < np; ++p) {
-            ws.well_potentials[p] = std::abs(potentials[p]);
+            // make sure the potentials are positive
+            ws.well_potentials[p] = std::max(0.0, potentials[p]);
         }
     }
 
@@ -1424,8 +1437,13 @@ namespace Opm {
                 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);
+            if (param_.solve_welleq_initially_ && well->isOperableAndSolvable()) {
+                try {
+                    well->solveWellEquation(ebosSimulator_, this->wellState(), this->groupState(), deferred_logger);
+                } catch (const std::exception& e) {
+                    const std::string msg = "Compute initial well solution for " + well->name() + " initially failed. Continue with the privious rates";
+                    deferred_logger.warning("WELL_INITIAL_SOLVE_FAILED", msg);
+                }
             }
         }
         updatePrimaryVariables(deferred_logger);

opm/simulators/wells/GasLiftSingleWell_impl.hpp

@@ -112,6 +112,10 @@ computeWellRates_(
             -potentials[this->water_pos_]);
         displayDebugMessage_(msg);
     }
+
+    for (auto& potential : potentials) {
+        potential = std::min(0.0, potential);
+    }
 }
 
 template<typename TypeTag>

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -253,6 +253,7 @@ namespace Opm
         if (this->wellIsStopped()) {
             return;
         }
+        this->operability_status_.has_negative_potentials = false;
 
         // If the well is pressure controlled the potential equals the rate.
         bool thp_controlled_well = false;
@@ -278,15 +279,16 @@ namespace Opm
         if (thp_controlled_well || bhp_controlled_well) {
 
             double total_rate = 0.0;
+            const double sign = this->isInjector() ? 1.0:-1.0;
             for (int phase = 0; phase < np; ++phase){
-                total_rate += ws.surface_rates[phase];
+                total_rate += sign * ws.surface_rates[phase];
             }
             // for pressure controlled wells the well rates are the potentials
             // if the rates are trivial we are most probably looking at the newly
             // opened well and we therefore make the affort of computing the potentials anyway.
-            if (std::abs(total_rate) > 0) {
+            if (total_rate > 0) {
                 for (int phase = 0; phase < np; ++phase){
-                    well_potentials[phase] = ws.surface_rates[phase];
+                    well_potentials[phase] = sign * ws.surface_rates[phase];
                 }
                 return;
             }
@@ -302,6 +304,19 @@ namespace Opm
         }
         deferred_logger.debug("Cost in iterations of finding well potential for well "
                               + this->name() + ": " + std::to_string(debug_cost_counter_));
+
+        const double sign = this->isInjector() ? 1.0:-1.0;
+        double total_potential = 0.0;
+        for (int phase = 0; phase < np; ++phase){
+            well_potentials[phase] *= sign;
+            total_potential += well_potentials[phase];
+        }
+        if (total_potential < 0.0) {
+            // wells with negative potentials are not operable
+            this->operability_status_.has_negative_potentials = true;
+            const std::string msg = std::string("well ") + this->name() + std::string(": has non negative potentials is not operable");
+            deferred_logger.info(msg);
+        }
     }
 
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -1868,6 +1868,7 @@ namespace Opm
             return;
         }
 
+        this->operability_status_.has_negative_potentials = false;
         // If the well is pressure controlled the potential equals the rate.
         bool thp_controlled_well = false;
         bool bhp_controlled_well = false;
@@ -1892,15 +1893,16 @@ namespace Opm
         if (thp_controlled_well || bhp_controlled_well) {
 
             double total_rate = 0.0;
+            const double sign = this->isInjector() ? 1.0:-1.0;
             for (int phase = 0; phase < np; ++phase){
-                total_rate += ws.surface_rates[phase];
+                total_rate += sign * ws.surface_rates[phase];
             }
             // for pressure controlled wells the well rates are the potentials
             // if the rates are trivial we are most probably looking at the newly
             // opened well and we therefore make the affort of computing the potentials anyway.
-            if (std::abs(total_rate) > 0) {
+            if (total_rate > 0) {
                 for (int phase = 0; phase < np; ++phase){
-                    well_potentials[phase] = ws.surface_rates[phase];
+                    well_potentials[phase] = sign * ws.surface_rates[phase];
                 }
                 return;
             }
@@ -1917,6 +1919,19 @@ namespace Opm
             // the well has a THP related constraint
             well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger, well_state);
         }
+
+        const double sign = this->isInjector() ? 1.0:-1.0;
+        double total_potential = 0.0;
+        for (int phase = 0; phase < np; ++phase){
+            well_potentials[phase] *= sign;
+            total_potential += well_potentials[phase];
+        }
+        if (total_potential < 0.0) {
+            // wells with negative potentials are not operable
+            this->operability_status_.has_negative_potentials = true;
+            const std::string msg = std::string("well ") + this->name() + std::string(": has negative potentials and is not operable");
+            deferred_logger.info(msg);
+        }
     }
 
 

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -174,6 +174,9 @@ public:
 
     void reportWellSwitching(const SingleWellState& ws, DeferredLogger& deferred_logger) const;
 
+    bool changedToOpenThisStep() const {
+        return this->changed_to_open_this_step_;
+    }
 protected:
     bool getAllowCrossFlow() const;
     double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
@@ -186,7 +189,7 @@ protected:
     // definition of the struct OperabilityStatus
     struct OperabilityStatus {
         bool isOperableAndSolvable() const {
-            if (!operable_under_only_bhp_limit || !solvable) {
+            if (!operable_under_only_bhp_limit || !solvable || has_negative_potentials) {
                 return false;
             } else {
                 return ( (isOperableUnderBHPLimit() || isOperableUnderTHPLimit()) );
@@ -221,6 +224,8 @@ protected:
         bool obey_bhp_limit_with_thp_limit = true;
         // the well is solveable
         bool solvable = true;
+        // the well have non positive potentials
+        bool has_negative_potentials = false;
     };
 
     OperabilityStatus operability_status_;
@@ -308,6 +313,8 @@ protected:
     const GuideRate* guide_rate_;
 
     std::vector< std::string> well_control_log_;
+
+    bool changed_to_open_this_step_ = false;
 };
 
 }

opm/simulators/wells/WellInterface_impl.hpp

@@ -349,7 +349,7 @@ namespace Opm
             }
             const int np = well_state_copy.numPhases();
             for (int p = 0; p < np; ++p) {
-                ws.well_potentials[p] = std::abs(potentials[p]);
+                ws.well_potentials[p] = std::max(0.0, potentials[p]);
             }
             this->updateWellTestState(well_state_copy.well(this->indexOfWell()), simulation_time, /*writeMessageToOPMLog=*/ false, welltest_state_temp, deferred_logger);
             this->closeCompletions(welltest_state_temp);
@@ -482,7 +482,23 @@ namespace Opm
                     this->operability_status_.solvable = false;
             }
         }
-
+        if (this->operability_status_.has_negative_potentials) {
+            auto well_state_copy = well_state;
+            std::vector<double> potentials;
+            try {
+                computeWellPotentials(ebosSimulator, well_state_copy, potentials, deferred_logger);
+            } catch (const std::exception& e) {
+                const std::string msg = std::string("well ") + this->name() + std::string(": computeWellPotentials() failed during attempt to recompute potentials for well : ") + e.what();
+                deferred_logger.info(msg);
+                this->operability_status_.has_negative_potentials = true;
+            }
+            auto& ws = well_state.well(this->indexOfWell());
+            const int np = well_state.numPhases();
+            for (int p = 0; p < np; ++p) {
+                ws.well_potentials[p] = std::max(0.0, potentials[p]);
+            }
+        }
+        this->changed_to_open_this_step_ = false;
         const bool well_operable = this->operability_status_.isOperableAndSolvable();
         if (!well_operable && old_well_operable) {
             if (this->well_ecl_.getAutomaticShutIn()) {
@@ -498,6 +514,7 @@ namespace Opm
             deferred_logger.info(" well " + this->name() + " gets REVIVED during iteration ");
             this->openWell();
             changed_to_stopped_this_step_ = false;
+            this->changed_to_open_this_step_ = true;
         }
 
         const auto& summary_state = ebosSimulator.vanguard().summaryState();