detecting whether some well control events happens to a well

When there is some events happen to a well, we use the control mode
 from the DECK, and update the WellState based on the new control model.
 Otherwise, we can use the control mode from the previous well state,
 and keep the values from the previous well state as an intial guess.

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -172,7 +172,7 @@ namespace Opm {
             const double p = fs.pressure(FluidSystem::oilPhaseIdx).value();
             cellPressures[cellIdx] = p;
         }
-        well_state_.init(wells(), cellPressures, &previous_well_state_, phase_usage_);
+        well_state_.init(wells(), cellPressures, wells_ecl_, timeStepIdx, &previous_well_state_, phase_usage_);
 
         // handling MS well related
         if (param_.use_multisegment_well_) { // if we use MultisegmentWell model
@@ -364,7 +364,7 @@ namespace Opm {
         if (nw > 0) {
             auto phaseUsage = phaseUsageFromDeck(eclState());
             size_t numCells = Opm::UgGridHelpers::numCells(grid());
-            well_state_.resize(wells, numCells, phaseUsage); //Resize for restart step
+            well_state_.resize(wells, numCells); //Resize for restart step
             wellsToState(restartValues.wells, phaseUsage, well_state_);
             previous_well_state_ = well_state_;
         }
@@ -938,18 +938,19 @@ namespace Opm {
             WellControls* wc = well->wellControls();
             const int control = well_controls_get_current(wc);
             well_state_.currentControls()[w] = control;
-            // TODO: for VFP control, the perf_densities are still zero here, investigate better
-            // way to handle it later.
-            well->updateWellStateWithTarget(well_state_);
 
-            // The wells are not considered to be newly added
+            if (well_state_.effectiveEventsHappen(w) ) {
-            // for next time step
+                well->updateWellStateWithTarget(well_state_);
-            if (well_state_.isNewWell(w) ) {
+            }
-                well_state_.setNewWell(w, false);
+
+            // there is no new well control change input within a report step,
+            // so next time step, the well does not consider to have effective events anymore
+            if (well_state_.effectiveEventsHappen(w) ) {
+                well_state_.setEffeciveEventHappen(w, false);
             }
         }  // end of for (int w = 0; w < nw; ++w)
 
-
+        updatePrimaryVariables();
     }
 
 
@@ -1185,8 +1186,10 @@ namespace Opm {
             // it will not change the control mode, only update the targets
             wellCollection().updateWellTargets(well_state_.wellRates());
 
+            // TODO: we should only do the well is involved in the update group targets
             for (auto& well : well_container_) {
                 well->updateWellStateWithTarget(well_state_);
+                well->updatePrimaryVariables(well_state_);
             }
         }
     }

opm/autodiff/MultisegmentWell_impl.hpp

@@ -373,8 +373,6 @@ namespace Opm
 
             break;
         } // end of switch
-
-        updatePrimaryVariables(well_state);
     }
 
 

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

@@ -460,7 +460,7 @@ namespace Opm
 
         const Wells* wells = wellsmanager.c_wells();
         size_t numCells = Opm::UgGridHelpers::numCells(grid);
-        wellstate.resize(wells, numCells, phaseUsage ); //Resize for restart step
+        wellstate.resize(wells, numCells); //Resize for restart step
         auto restart_values = eclIO_->loadRestart(solution_keys, extra_keys);
 
         solutionToSim( restart_values, phaseUsage, simulatorstate );

opm/autodiff/StandardWell_impl.hpp

@@ -1214,8 +1214,6 @@ namespace Opm
 
             break;
         } // end of switch
-
-        updatePrimaryVariables(well_state);
     }
 
 
@@ -1908,13 +1906,14 @@ namespace Opm
         } else {
             // the well has a THP related constraint
             // checking whether a well is newly added, it only happens at the beginning of the report step
-            if ( !well_state.isNewWell(index_of_well_) ) {
+            if ( !well_state.effectiveEventsHappen(index_of_well_) ) {
                 for (int p = 0; p < np; ++p) {
                     // This is dangerous for new added well
                     // since we are not handling the initialization correctly for now
                     well_potentials[p] = well_state.wellRates()[index_of_well_ * np + p];
                 }
-            } else {
+            } else
+            {
                 // We need to generate a reasonable rates to start the iteration process
                 computeWellRatesWithBhp(ebosSimulator, bhp, well_potentials);
                 for (double& value : well_potentials) {

opm/autodiff/WellInterface_impl.hpp

@@ -438,6 +438,7 @@ namespace Opm
 
         if (updated_control_index != old_control_index) { //  || well_collection_->groupControlActive()) {
             updateWellStateWithTarget(well_state);
+            updatePrimaryVariables(well_state);
         }
     }
 

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -62,7 +62,9 @@ namespace Opm
         /// Allocate and initialize if wells is non-null.  Also tries
         /// to give useful initial values to the bhp(), wellRates()
         /// and perfPhaseRates() fields, depending on controls
-        void init(const Wells* wells, const std::vector<double>& cellPressures, const WellStateFullyImplicitBlackoil* prevState, const PhaseUsage& pu)
+        void init(const Wells* wells, const std::vector<double>& cellPressures,
+                  const std::vector<const Well*>& wells_ecl, const int report_step,
+                  const WellStateFullyImplicitBlackoil* prevState, const PhaseUsage& pu)
         {
             // call init on base class
             BaseType :: init(wells, cellPressures);
@@ -84,15 +86,27 @@ namespace Opm
             well_dissolved_gas_rates_.resize(nw, 0.0);
             well_vaporized_oil_rates_.resize(nw, 0.0);
 
-            is_new_well_.resize(nw, true);
+            // checking whether some effective well control happens
-            if (prevState && !prevState->wellMap().empty()) {
+            effective_events_happen_.resize(nw);
-                const auto& end = prevState->wellMap().end();
+            for (int w = 0; w <nw; ++w) {
-                for (int w = 0; w < nw; ++w) {
+                const int nw_wells_ecl = wells_ecl.size();
-                    const auto& it = prevState->wellMap().find( wells->name[w]);
+                int index_well_ecl = 0;
-                    if (it != end) {
+                const std::string well_name(wells->name[w]);
-                        is_new_well_[w] = false;
+                for (; index_well_ecl < nw_wells_ecl; ++index_well_ecl) {
+                    if (well_name == wells_ecl[index_well_ecl]->name()) {
+                        break;
                     }
                 }
+
+                // It should be able to find in wells_ecl.
+                if (index_well_ecl == nw_wells_ecl) {
+                    OPM_THROW(std::logic_error, "Could not find well " << well_name << " in wells_ecl ");
+                }
+
+                const Well* well_ecl = wells_ecl[index_well_ecl];
+                // PRODUCTION_UPDATE, INJECTION_UPDATE, WELL_STATUS_CHANGE
+                // 16 + 32 + 128
+                effective_events_happen_[w] = (well_ecl->hasEvent(176, report_step) );
             }
 
             // Ensure that we start out with zero rates by default.
@@ -117,9 +131,11 @@ namespace Opm
                 }
             }
 
-            // Initialize current_controls_.
+            // set up the current well controls, if there are some effective events happened to the well,
-            // The controls set in the Wells object are treated as defaults,
+            // we use the control specified in the DECK
-            // and also used for initial values.
+            // if not, we will just use the control from the prevState
+            // the second part will be done when copying the data from  prevState
+            // here, we initilaize all the current_controls_ based on the control mode from the DECK
             current_controls_.resize(nw);
             for (int w = 0; w < nw; ++w) {
                 current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
@@ -136,7 +152,7 @@ namespace Opm
                 typedef typename WellMapType :: const_iterator const_iterator;
                 const_iterator end = prevState->wellMap().end();
                 for (int w = 0; w < nw; ++w) {
-                    std::string name( wells->name[ w ] );
+                    const std::string name( wells->name[ w ] );
                     const_iterator it = prevState->wellMap().find( name );
                     if( it != end )
                     {
@@ -149,6 +165,14 @@ namespace Opm
                         // thp
                         thp()[ newIndex ] = prevState->thp()[ oldIndex ];
 
+                        // if there is no effective control event happens to the well, we use the current_controls_ from prevState
+                        if (!effective_events_happen_[w]) {
+                            current_controls_[ newIndex ] = prevState->currentControls()[ oldIndex ];
+                            // also change the one in the WellControls
+                            // TODO: checking if this is appropriate
+                            well_controls_set_current(wells->ctrls[w], current_controls_[ newIndex ]);
+                        }
+
                         // wellrates
                         for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
                         {
@@ -230,10 +254,24 @@ namespace Opm
             }
         }
 
-        void resize(const Wells* wells, size_t numCells, const PhaseUsage& pu)
+        void resize(const Wells* wells, size_t numCells)
         {
             std::vector<double> tmp(numCells, 0.0); // <- UGLY HACK to pass the size
-            init(wells, tmp, nullptr, pu);
+            BaseType :: init(wells, tmp);
+
+            const int np = wells->number_of_phases;
+            const int nw = wells->number_of_wells;
+            const int nperf = wells->well_connpos[nw];
+
+            well_reservoir_rates_.resize(nw * np, 0.0);
+            well_dissolved_gas_rates_.resize(nw, 0.0);
+            well_vaporized_oil_rates_.resize(nw, 0.0);
+            effective_events_happen_.resize(nw);
+            perfphaserates_.resize(nperf * np, 0.0);
+            current_controls_.resize(nw);
+            perfRateSolvent_.resize(nperf, 0.0);
+            top_segment_index_.resize(nw);
+            // TODO: segrates_ and segpress_ are not taken care of here
         }
 
         /// Allocate and initialize if wells is non-null.  Also tries
@@ -296,8 +334,6 @@ namespace Opm
                 current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
             }
 
-            is_new_well_.resize(nw, true);
-
             perfRateSolvent_.clear();
             perfRateSolvent_.resize(nperf, 0.0);
 
@@ -312,9 +348,6 @@ namespace Opm
                     const_iterator it = prevState.wellMap().find( name );
                     if( it != end )
                     {
-                        // this is not a new added well
-                        is_new_well_[w] = false;
-
                         const int oldIndex = (*it).second[ 0 ];
                         const int newIndex = w;
 
@@ -707,13 +740,13 @@ namespace Opm
         }
 
 
-        bool isNewWell(const int w) const {
+        bool effectiveEventsHappen(const int w) const {
-            return is_new_well_[w];
+            return effective_events_happen_[w];
         }
 
 
-        void setNewWell(const int w, const bool is_new_well) {
+        void setEffeciveEventHappen(const int w, const bool effective_events_happen) {
-            is_new_well_[w] = is_new_well;
+            effective_events_happen_[w] = effective_events_happen;
         }
 
 
@@ -802,11 +835,10 @@ namespace Opm
         // should be zero for injection wells
         std::vector<double> well_vaporized_oil_rates_;
 
-        // marking whether the well is just added
+        // some events happens to the well, like this well is a new well
-        // for newly added well, the current initialized rates from WellState
+        // or new well control keywords happens
-        // will have very wrong compositions for production wells, will mostly cause
+        // \Note: for now, only WCON* keywords, and well status change is considered
-        // problem with VFP interpolation
+        std::vector<bool> effective_events_happen_;
-        std::vector<bool> is_new_well_;
 
         // MS well related
         // for StandardWell, the number of segments will be one