Merge pull request #181 from atgeirr/num-wells-can-change

Allow number of wells to increase during simulation

opm/autodiff/SimulatorFullyImplicitBlackoil_impl.hpp

@@ -253,7 +253,7 @@ namespace Opm
     SimulatorReport SimulatorFullyImplicitBlackoil<T>::Impl::run(SimulatorTimer& timer,
                                                                  BlackoilState& state)
     {
-        WellStateFullyImplicitBlackoil well_state;
+        WellStateFullyImplicitBlackoil prev_well_state;
 
         // Create timers and file for writing timing info.
         Opm::time::StopWatch solver_timer;
@@ -281,13 +281,15 @@ namespace Opm
                                        Opm::UgGridHelpers::beginFaceCentroids(grid_),
                                        props_.permeability());
 
-            const Wells *wells = wells_manager.c_wells();
+            const Wells* wells = wells_manager.c_wells();
+            WellStateFullyImplicitBlackoil well_state;
+            well_state.init(wells, state);
 
             if (timer.currentStepNum() == 0) {
-                well_state.init(wells, state);
                 output_writer_.writeInit(timer);
             } else {
-                // TODO: add a function to update the well_state here.
+                // Transfer previous well state tu current.
+                well_state.partialCopy(prev_well_state, *wells, prev_well_state.numWells());
             }
 
             if (output_ && (timer.currentStepNum() % output_interval_ == 0)) {
@@ -334,6 +336,7 @@ namespace Opm
             output_writer_.writeTimeStep(timer, state, well_state.basicWellState());
 
             ++timer;
+            prev_well_state = well_state;
         }
 
         total_timer.stop();

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -24,6 +24,7 @@
 #include <opm/core/wells.h>
 #include <opm/core/well_controls.h>
 #include <opm/core/simulator/WellState.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
 #include <vector>
 #include <cassert>
 
@@ -102,6 +103,53 @@ namespace Opm
             return basic_well_state_;
         }
 
+        /// The number of wells present.
+        int numWells() const
+        {
+            return bhp().size();
+        }
+
+        /// The number of phases present.
+        int numPhases() const
+        {
+            return wellRates().size() / numWells();
+        }
+
+        /// 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,
+        /// that the number of perforations associated with the wells
+        /// is unchanging, and that both objects contain at least
+        /// num_wells_to_copy wells.
+        void partialCopy(const WellStateFullyImplicitBlackoil& source,
+                         const Wells& wells,
+                         const int num_wells_to_copy)
+        {
+            if (numPhases() != source.numPhases()) {
+                OPM_THROW(std::logic_error, "partialCopy(): source and destination have different number of phases.");
+            }
+            if (num_wells_to_copy > numWells() || num_wells_to_copy > source.numWells()) {
+                OPM_THROW(std::logic_error, "partialCopy(): trying to copy too many wells.");
+            }
+            // bhp
+            std::copy(source.bhp().begin(),
+                      source.bhp().begin() + num_wells_to_copy,
+                      bhp().begin());
+            // wellRates
+            std::copy(source.wellRates().begin(),
+                      source.wellRates().begin() + numPhases()*num_wells_to_copy,
+                      wellRates().begin());
+            // perfPhaseRates
+            const int num_perfs_to_copy = wells.well_connpos[num_wells_to_copy];
+            std::copy(source.perfPhaseRates().begin(),
+                      source.perfPhaseRates().begin() + numPhases()*num_perfs_to_copy,
+                      perfPhaseRates().begin());
+            // currentControls
+            std::copy(source.currentControls().begin(),
+                      source.currentControls().begin() + num_wells_to_copy,
+                      currentControls().begin());
+        }
+
     private:
         WellState basic_well_state_;
         std::vector<double> perfphaserates_;