tracking if wells are newly added in WellStateFullyImplicitBlackoil

For a newly added well, the initialized rate cans cause big problem when
applied to VFP interploation when THP control is involved.

opm/autodiff/StandardWellsDense.hpp

@@ -249,7 +249,7 @@ enum WellVariablePositions {
             void
             computeWellPotentials(const Simulator& ebosSimulator,
                                   const WellState& well_state,
-                                  std::vector<double>& well_potentials)  const;
+                                  std::vector<double>& well_potentials) const;
 
             // TODO: some preparation work, mostly related to group control and RESV,
             // at the beginning of each time step (Not report step)

opm/autodiff/StandardWellsDense_impl.hpp

@@ -1729,7 +1729,7 @@ namespace Opm {
     StandardWellsDense<FluidSystem, BlackoilIndices, ElementContext, MaterialLaw>::
     computeWellPotentials(const Simulator& ebosSimulator,
                           const WellState& well_state,
-                          std::vector<double>& well_potentials)  const
+                          std::vector<double>& well_potentials) const
     {
 
         // number of wells and phases
@@ -1755,10 +1755,21 @@ namespace Opm {
                 computeWellRatesWithBhp(ebosSimulator, bhp, w, potentials);
 
             } else { // the well has a THP related constraint
-                for (int p = 0; p < np; ++p) {
+                // checking whether a well is newly added, it only happens at the beginning of the report step
-                    // TODO: this is dangerous for new added well
+                if ( !well_state.isNewWell(w) ) {
-                    // since we are not handling the initialization correctly for now
+                    for (int p = 0; p < np; ++p) {
-                    potentials[p] = well_state.wellRates()[w * np + p];
+                        // TODO: this is dangerous for new added well
+                        // since we are not handling the initialization correctly for now
+                        potentials[p] = well_state.wellRates()[w * np + p];
+                    }
+                } else {
+                    // We need to generate a reasonable rates to start the iteration process
+                    computeWellRatesWithBhp(ebosSimulator, bhp, w, potentials);
+                    for (double& value : potentials) {
+                        // make the value a little safer in case the BHP limits are default ones
+                        // TODO: a better way should be a better rescaling based on the investigation of the VFP table.
+                        value *= 0.001;
+                    }
                 }
 
                 potentials = computeWellPotentialWithTHP(ebosSimulator, w, bhp, potentials);
@@ -1875,6 +1886,12 @@ namespace Opm {
             const int control = well_controls_get_current(wc);
             well_state.currentControls()[w] = control;
             updateWellStateWithTarget(wc, control, w, well_state);
+
+            // The wells are not considered to be newly added
+            // for next time step
+            if (well_state.isNewWell(w) ) {
+                well_state.setNewWell(w, false);
+            }
         }
     }
 

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -106,6 +106,8 @@ namespace Opm
             well_potentials_.clear();
             well_potentials_.resize(nperf * np, 0.0);
 
+            is_new_well_.resize(nw, true);
+
             // intialize wells that have been there before
             // order may change so the mapping is based on the well name
             if( ! prevState.wellMap().empty() )
@@ -117,6 +119,9 @@ 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;
 
@@ -273,10 +278,26 @@ namespace Opm
             return res;
         }
 
+
+        bool isNewWell(const int w) const {
+            return is_new_well_[w];
+        }
+
+
+        void setNewWell(const int w, const bool is_new_well) {
+            is_new_well_[w] = is_new_well;
+        }
+
     private:
         std::vector<double> perfphaserates_;
         std::vector<int> current_controls_;
         std::vector<double> well_potentials_;
+
+        // marking whether the well is just added
+        // for newly added well, the current initialized rates from WellState
+        // will have very wrong compsitions for productions wells, will mostly cause
+        // problem with VFP interpolation
+        std::vector<bool> is_new_well_;
     };
 
 } // namespace Opm