Use well to name mapping when copying the old well state to the new

With the introduction of shut wells the same ordering in the well states
can no longer be assumed. The well names is instead used to map the old
well state to the new. Also the partial_copying is moved into the
initialization.

Tested on SPE1, SPE3, SPE9 and Norne. (Do not change the SPEs and is
nessesary for the Norne)

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -323,11 +323,7 @@ namespace Opm
                                        props_.permeability());
             const Wells* wells = wells_manager.c_wells();
             WellStateFullyImplicitBlackoil well_state;
-            well_state.init(wells, state);
-            if (timer.currentStepNum() != 0) {
-                // Transfer previous well state to current.
-                well_state.partialCopy(prev_well_state, *wells, prev_well_state.numWells());
-            }
+            well_state.init(wells, state, prev_well_state);
 
             // Output state at start of time step.
             if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -27,6 +27,10 @@
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <vector>
 #include <cassert>
+#include <string>
+#include <utility>
+#include <map>
+#include <algorithm>
 
 namespace Opm
 {
@@ -36,12 +40,18 @@ namespace Opm
     class WellStateFullyImplicitBlackoil
     {
     public:
+        typedef std::pair<int,int> mapentry_t;
+        typedef std::map< std::string, mapentry_t > WellMapType;
+
         /// 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
-        template <class State>
-        void init(const Wells* wells, const State& state)
+        template <class State, class PrevState>
+        void init(const Wells* wells, const State& state, const PrevState& prevState)
         {
+            // clear old name mapping
+            wellMap_.clear();
+
             if (wells == 0) {
                 return;
             }
@@ -58,6 +68,12 @@ namespace Opm
             for (int w = 0; w < nw; ++w) {
                 assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
                 const WellControls* ctrl = wells->ctrls[w];
+                std::string name( wells->name[ w ] );
+                assert( name.size() > 0 );
+                mapentry_t& wellMapEntry = wellMap_[ name ];
+                assert( wellMapEntry.size() == 0 );
+                wellMapEntry.first  = w ;
+                wellMapEntry.second = wells->well_connpos[w ] ;
                 if (well_controls_well_is_shut(ctrl)) {
                     // Shut well: perfphaserates_ are all zero.
                 } else {
@@ -79,6 +95,43 @@ namespace Opm
             for (int w = 0; w < nw; ++w) {
                 current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
             }
+
+            // intialize wells that have been there before
+            // order may change so the mapping is based on the well name
+            if( prevState.wellMap().size() > 0 ) 
+            {
+                typedef typename WellMapType :: iterator iterator;
+                const iterator end = prevState.wellMap().end();
+                for (int w = 0; w < nw; ++w) {
+                    std::string name( wells->name[ w ] );
+                    iterator it = prevState.wellMap().find( name );
+                    if( it != end )
+                    {
+                        const int oldIndex = (*it).second.first;
+                        const int newIndex = w;
+
+                        // bhp
+                        bhp()[ newIndex ] = prevState.bhp()[ oldIndex ];
+
+                        // wellrates
+                        for( int i=0, idx=newIndex*np, oldidx=oldIndex*np; i<np; ++i, ++idx, ++oldidx )
+                        {
+                            wellRates()[ idx ] = prevState.wellRates()[ oldidx ];
+                        }
+
+                        // perfPhaseRates
+                        int oldPerf = (*it).second.second * np;
+                        for (int perf = wells->well_connpos[ newIndex ]*np;
+                             perf < wells->well_connpos[ newIndex + 1]*np; ++perf, ++oldPerf )
+                        {
+                            perfPhaseRates()[ perf ] = prevState.perfPhaseRates()[ oldPerf ];
+                        }
+
+                        // currentControls
+                        currentControls()[ newIndex ] = prevState.currentControls()[ oldIndex ];
+                    }
+                }
+            }
         }
 
         /// One bhp pressure per well.
@@ -115,6 +168,8 @@ namespace Opm
             return wellRates().size() / numWells();
         }
 
+        WellMapType& wellMap() const { return wellMap_; }
+
         /// Copy data for the first num_wells_to_copy from source,
         /// overwriting any data in this object associated with those
         /// wells. Assumes that the number of phases are the same,
@@ -125,6 +180,8 @@ namespace Opm
                          const Wells& wells,
                          const int num_wells_to_copy)
         {
+            std::abort();
+
             if (numPhases() != source.numPhases()) {
                 OPM_THROW(std::logic_error, "partialCopy(): source and destination have different number of phases.");
             }
@@ -154,6 +211,7 @@ namespace Opm
         WellState basic_well_state_;
         std::vector<double> perfphaserates_;
         std::vector<int> current_controls_;
+        mutable WellMapType wellMap_;
     };
 
 } // namespace Opm