implement WTEST support

compareECLFiles.cmake

@@ -259,6 +259,12 @@ add_test_compareECLFiles(CASENAME spe5
                          REL_TOL ${coarse_rel_tol}
                          TEST_ARGS max_iter=20)
 
+add_test_compareECLFiles(CASENAME wecon
+                         FILENAME 3D_WECON
+                         SIMULATOR flow
+                         ABS_TOL ${abs_tol}
+                         REL_TOL ${coarse_rel_tol})
+
 # Restart tests
 opm_set_test_driver(${PROJECT_SOURCE_DIR}/tests/run-restart-regressionTest.sh "")
 

opm/autodiff/BlackoilModelEbos.hpp

@@ -206,7 +206,7 @@ namespace Opm {
             wasSwitched_.resize(numDof);
             std::fill(wasSwitched_.begin(), wasSwitched_.end(), false);
 
-            wellModel().beginTimeStep();
+            wellModel().beginTimeStep(timer.reportStepNum(), timer.simulationTimeElapsed());
 
             if (param_.update_equations_scaling_) {
                 std::cout << "equation scaling not suported yet" << std::endl;
@@ -340,7 +340,7 @@ namespace Opm {
         /// \param[in] timer                  simulation timer
         void afterStep(const SimulatorTimerInterface& OPM_UNUSED timer)
         {
-            wellModel().timeStepSucceeded();
+            wellModel().timeStepSucceeded(timer.simulationTimeElapsed());
             aquiferModel().timeStepSucceeded(timer);
             ebosSimulator_.problem().endTimeStep();
 

opm/autodiff/BlackoilWellModel.hpp

@@ -2,7 +2,7 @@
   Copyright 2016 SINTEF ICT, Applied Mathematics.
   Copyright 2016 - 2017 Statoil ASA.
   Copyright 2017 Dr. Blatt - HPC-Simulation-Software & Services
-  Copyright 2016 - 2017 IRIS AS
+  Copyright 2016 - 2018 IRIS AS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -33,6 +33,7 @@
 #include <tuple>
 
 #include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/WellTestState.hpp>
 
 #include <opm/core/wells.h>
 #include <opm/core/wells/DynamicListEconLimited.hpp>
@@ -141,7 +142,8 @@ namespace Opm {
             // compute the well fluxes and assemble them in to the reservoir equations as source terms
             // and in the well equations.
             void assemble(const int iterationIdx,
-                                     const double dt);
+                          const double dt,
+                          bool wtest = false);
 
             // substract Binv(D)rw from r;
             void apply( BVector& r) const;
@@ -175,9 +177,9 @@ namespace Opm {
             void setRestartWellState(const WellState& well_state);
 
             // called at the beginning of a time step
-            void beginTimeStep();
+            void beginTimeStep(const int timeStepIdx,const double simulationTime);
             // called at the end of a time step
-            void timeStepSucceeded();
+            void timeStepSucceeded(const double& simulationTime);
 
             // called at the beginning of a report step
             void beginReportStep(const int time_step);
@@ -236,7 +238,7 @@ namespace Opm {
             using ConvergenceReport = typename WellInterface<TypeTag>::ConvergenceReport;
 
             // create the well container
-            std::vector<WellInterfacePtr > createWellContainer(const int time_step) const;
+            std::vector<WellInterfacePtr > createWellContainer(const int time_step);
 
             WellState well_state_;
             WellState previous_well_state_;
@@ -254,12 +256,13 @@ namespace Opm {
             std::vector<double> depth_;
             bool initial_step_;
 
-            DynamicListEconLimited dynamic_list_econ_limited_;
             std::unique_ptr<RateConverterType> rateConverter_;
             std::unique_ptr<VFPProperties> vfp_properties_;
 
             SimulatorReport last_report_;
 
+            WellTestState wellTestState_;
+
             // used to better efficiency of calcuation
             mutable BVector scaleAddRes_;
 
@@ -344,10 +347,12 @@ namespace Opm {
             bool localWellsActive() const;
 
             /// upate the dynamic lists related to economic limits
-            void updateListEconLimited(DynamicListEconLimited& list_econ_limited) const;
+            void updateListEconLimited(const double& simulationTime, WellTestState& wellTestState) const;
 
             void updatePerforationIntensiveQuantities();
 
+            void wellTesting(const int timeStepIdx,const double simulationTime);
+
         };
 
 

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -45,6 +45,10 @@ namespace Opm {
         wells_ecl_ = schedule().getWells(timeStepIdx);
 
         // Create wells and well state.
+        // Pass empty dynamicListEconLimited class
+        // The closing of wells due to limites is
+        // handled by the wellTestState class
+        DynamicListEconLimited dynamic_list_econ_limited;
         wells_manager_.reset( new WellsManager (eclState,
                                                 schedule(),
                                                 timeStepIdx,
@@ -54,7 +58,7 @@ namespace Opm {
                                                 Opm::UgGridHelpers::dimensions(grid),
                                                 Opm::UgGridHelpers::cell2Faces(grid),
                                                 Opm::UgGridHelpers::beginFaceCentroids(grid),
-                                                dynamic_list_econ_limited_,
+                                                dynamic_list_econ_limited,
                                                 grid.comm().size() > 1,
                                                 defunct_well_names) );
 
@@ -103,11 +107,40 @@ namespace Opm {
             }
         }
 
+        // update the previous well state. This is used to restart failed steps.
+        previous_well_state_ = well_state_;
+
+
+        if (wellCollection().havingVREPGroups() ) {
+            rateConverter_->template defineState<ElementContext>(ebosSimulator_);
+        }
+
         // Compute reservoir volumes for RESV controls.
         rateConverter_.reset(new RateConverterType (phase_usage_,
                                          std::vector<int>(number_of_cells_, 0)));
         computeRESV(timeStepIdx);
 
+        // update VFP properties
+        vfp_properties_.reset (new VFPProperties (
+                                   schedule().getVFPInjTables(timeStepIdx),
+                                   schedule().getVFPProdTables(timeStepIdx)) );
+
+
+
+    }
+
+
+    // called at the beginning of a time step
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    beginTimeStep(const int timeStepIdx, const double simulationTime) {
+        well_state_ = previous_well_state_;
+
+
+        // test wells
+        wellTesting(timeStepIdx, simulationTime);
+
         // create the well container
         well_container_ = createWellContainer(timeStepIdx);
 
@@ -121,6 +154,8 @@ namespace Opm {
         // calculate the efficiency factors for each well
         calculateEfficiencyFactors();
 
+        const Grid& grid = ebosSimulator_.vanguard().grid();
+
         if (has_polymer_)
         {
             if (PolymerModule::hasPlyshlog()) {
@@ -128,32 +163,101 @@ namespace Opm {
             }
         }
 
-        // update VFP properties
-        vfp_properties_.reset (new VFPProperties (
-                                   schedule().getVFPInjTables(timeStepIdx),
-                                   schedule().getVFPProdTables(timeStepIdx)) );
-
         for (auto& well : well_container_) {
             well->setVFPProperties(vfp_properties_.get());
         }
 
-        // update the previous well state. This is used to restart failed steps.
+        // Close wells and connections due to economical reasons
-        previous_well_state_ = well_state_;
+        for (auto& well : well_container_) {
-
+            well->closeWellsAndCompletions(wellTestState_);
+        }
 
     }
 
 
-    // called at the beginning of a time step
     template<typename TypeTag>
     void
-    BlackoilWellModel<TypeTag>::
+    BlackoilWellModel<TypeTag>::wellTesting(const int timeStepIdx, const double simulationTime) {
-    beginTimeStep() {
+        const auto& wtest_config = schedule().wtestConfig(timeStepIdx);
-        well_state_ = previous_well_state_;
+        const auto& wellsForTesting = wellTestState_.updateWell(wtest_config, simulationTime);
 
-        if (wellCollection().havingVREPGroups() ) {
+        // Do the well testing if enabled
-            rateConverter_->template defineState<ElementContext>(ebosSimulator_);
+        if (!initial_step_ && wtest_config.size() > 0 && wellsForTesting.size() > 0) {
+            // solve the well equation isolated from the reservoir.
+            const int numComp = numComponents();
+            std::vector< Scalar > B_avg( numComp, Scalar() );
+            computeAverageFormationFactor(B_avg);
+            std::vector<WellInterfacePtr> well_container;
 
+            well_container.reserve(wellsForTesting.size());
+            for (auto& testWell : wellsForTesting) {
+                const std::string msg = std::string("well ") + testWell.first + std::string(" will be tested");
+                OpmLog::info(msg);
+
+                // finding the location of the well in wells_ecl
+                const int nw_wells_ecl = wells_ecl_.size();
+                int index_well = 0;
+                for (; index_well < nw_wells_ecl; ++index_well) {
+                    if (testWell.first == wells_ecl_[index_well]->name()) {
+                        break;
+                    }
+                }
+
+                // It should be able to find in wells_ecl.
+                if (index_well == nw_wells_ecl) {
+                    OPM_THROW(std::logic_error, "Could not find well " << testWell.first << " in wells_ecl ");
+                }
+                const Well* well_ecl = wells_ecl_[index_well];
+
+
+                // Find the index in the wells() struct
+                const int nw = numWells();
+                int wellidx = -999;
+                for (int w = 0; w < nw; ++w) {
+                    if (testWell.first == std::string(wells()->name[w])) {
+                        wellidx = w;
+                        break;
+                    }
+                }
+
+
+                // Use the pvtRegionIdx from the top cell
+                const int well_cell_top = wells()->well_cells[wells()->well_connpos[wellidx]];
+                const int pvtreg = pvt_region_idx_[well_cell_top];
+
+                //WellInterface<TypeTag> well(well_ecl, timeStepIdx, wells(), param_, *rateConverter_, pvtreg, numComponents() );
+
+                if ( !well_ecl->isMultiSegment(timeStepIdx) || !param_.use_multisegment_well_) {
+                    well_container.emplace_back(new StandardWell<TypeTag>(well_ecl, timeStepIdx, wells(),
+                                                                          param_, *rateConverter_, pvtreg, numComponents() ) );
+                } else {
+                    well_container.emplace_back(new MultisegmentWell<TypeTag>(well_ecl, timeStepIdx, wells(),
+                                                                              param_, *rateConverter_, pvtreg, numComponents() ) );
+                }
+            }
+
+            for (auto& well : well_container) {
+                WellTestState wellTestStateForTheWellTest;
+                well->init(&phase_usage_, depth_, gravity_, number_of_cells_);
+                const std::string well_name = well->name();
+                const WellNode& well_node = wellCollection().findWellNode(well_name);
+                const double well_efficiency_factor = well_node.getAccumulativeEfficiencyFactor();
+                well->setWellEfficiencyFactor(well_efficiency_factor);
+                well->setVFPProperties(vfp_properties_.get());
+                well->solveWellEq(ebosSimulator_, well_state_, /*dt (not relevant for well test) =*/ 1.0, B_avg, terminal_output_);
+                well->updateListEconLimited(well_state_, simulationTime, wellTestStateForTheWellTest);
+                // update wellTestState if the well test succeeds
+                if (!wellTestStateForTheWellTest.hasWell(well->name(), WellTestConfig::Reason::ECONOMIC)) {
+                    wellTestState_.openWell(well->name());
+                    const std::string msg = std::string("well ") + well->name() + std::string(" is re-opened");
+                    OpmLog::info(msg);
+                    // also reopen completions
+                    for (int completionIdx = 0;  completionIdx <well->numberOfCompletions(); ++completionIdx) {
+                        if (!wellTestStateForTheWellTest.hasCompletion(well->name(), completionIdx))
+                            wellTestState_.dropCompletion(well->name(), completionIdx);
+                    }
+                }
+            }
         }
     }
 
@@ -168,10 +272,6 @@ namespace Opm {
     void
     BlackoilWellModel<TypeTag>::
     endReportStep() {
-        // update the list contanining information of closed wells
-        // and connections due to economical limits
-        // Used by the wellManager
-        updateListEconLimited(dynamic_list_econ_limited_);
     }
 
     // called at the end of a report step
@@ -184,20 +284,20 @@ namespace Opm {
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
-    timeStepSucceeded() {
+    timeStepSucceeded(const double& simulationTime) {
         // TODO: when necessary
         rateConverter_->template defineState<ElementContext>(ebosSimulator_);
         for (const auto& well : well_container_) {
             well->calculateReservoirRates(well_state_);
         }
-
+        updateListEconLimited(simulationTime, wellTestState_);
         previous_well_state_ = well_state_;
     }
 
     template<typename TypeTag>
     std::vector<typename BlackoilWellModel<TypeTag>::WellInterfacePtr >
     BlackoilWellModel<TypeTag>::
-    createWellContainer(const int time_step) const
+    createWellContainer(const int time_step)
     {
         std::vector<WellInterfacePtr> well_container;
 
@@ -227,6 +327,16 @@ namespace Opm {
 
                 const Well* well_ecl = wells_ecl_[index_well];
 
+                // well is shut due to economical reasons
+                if (wellTestState_.hasWell(well_name, WellTestConfig::Reason::ECONOMIC) && well_ecl->getAutomaticShutIn() ) {
+                    well_state_.bhp()[w] = 0;
+                    const int np = numPhases();
+                    for (int p = 0; p < np; ++p) {
+                        well_state_.wellRates()[np * w + p] = 0;
+                    }
+                    continue;
+                }
+
                 // Use the pvtRegionIdx from the top cell
                 const int well_cell_top = wells()->well_cells[wells()->well_connpos[w]];
                 const int pvtreg = pvt_region_idx_[well_cell_top];
@@ -251,7 +361,8 @@ namespace Opm {
     void
     BlackoilWellModel<TypeTag>::
     assemble(const int iterationIdx,
-             const double dt)
+             const double dt,
+             bool onlyDoTheWellTest)
     {
 
 
@@ -275,6 +386,9 @@ namespace Opm {
         if (param_.solve_welleq_initially_ && iterationIdx == 0) {
             // solve the well equations as a pre-processing step
             last_report_ = solveWellEq(dt);
+            if (onlyDoTheWellTest)
+                return;
+
             if (initial_step_) {
                 // update the explicit quantities to get the initial fluid distribution in the well correct.
                 calculateExplicitQuantities();
@@ -301,8 +415,8 @@ namespace Opm {
     assembleWellEq(const double dt,
                    bool only_wells)
     {
-        for (int w = 0; w < numWells(); ++w) {
+        for (auto& well : well_container_) {
-            well_container_[w]->assembleWellEq(ebosSimulator_, dt, well_state_, only_wells);
+            well->assembleWellEq(ebosSimulator_, dt, well_state_, only_wells);
         }
     }
 
@@ -392,13 +506,9 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     resetWellControlFromState() const
     {
-        const int        nw   = numWells();
+        for (auto& well : well_container_) {
-
+            WellControls* wc = well->wellControls();
-        assert(nw == int(well_container_.size()) );
+            well_controls_set_current( wc, well_state_.currentControls()[well->indexOfWell()]);
-
-        for (int w = 0; w < nw; ++w) {
-            WellControls* wc = well_container_[w]->wellControls();
-            well_controls_set_current( wc, well_state_.currentControls()[w]);
         }
     }
 
@@ -475,6 +585,7 @@ namespace Opm {
         do {
             assembleWellEq(dt, true);
 
+            //std::cout << "well convergence only wells " << std::endl;
             converged = getWellConvergence(B_avg);
 
             // checking whether the group targets are converged
@@ -632,10 +743,10 @@ namespace Opm {
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
-    updateListEconLimited(DynamicListEconLimited& list_econ_limited) const
+    updateListEconLimited(const double& simulationTime, WellTestState& wellTestState) const
     {
         for (const auto& well : well_container_) {
-            well->updateListEconLimited(well_state_, list_econ_limited);
+            well->updateListEconLimited(well_state_, simulationTime, wellTestState);
         }
     }
 
@@ -651,13 +762,13 @@ namespace Opm {
         const int np = numPhases();
         well_potentials.resize(nw * np, 0.0);
 
-        for (int w = 0; w < nw; ++w) {
+        for (const auto& well : well_container_) {
             std::vector<double> potentials;
-            well_container_[w]->computeWellPotentials(ebosSimulator_, well_state_, potentials);
+            well->computeWellPotentials(ebosSimulator_, well_state_, potentials);
 
             // putting the sucessfully calculated potentials to the well_potentials
             for (int p = 0; p < np; ++p) {
-                well_potentials[w * np + p] = std::abs(potentials[p]);
+                well_potentials[well->indexOfWell() * np + p] = std::abs(potentials[p]);
             }
         } // end of for (int w = 0; w < nw; ++w)
     }
@@ -687,13 +798,14 @@ namespace Opm {
         prepareGroupControl();
 
         // since the controls are all updated, we should update well_state accordingly
-        for (int w = 0; w < numWells(); ++w) {
+        for (const auto& well : well_container_) {
-            WellControls* wc = well_container_[w]->wellControls();
+            const int w = well->indexOfWell();
+            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_container_[w]->updateWellStateWithTarget(well_state_);
+            well->updateWellStateWithTarget(well_state_);
 
             // The wells are not considered to be newly added
             // for next time step
@@ -702,6 +814,7 @@ namespace Opm {
             }
         }  // end of for (int w = 0; w < nw; ++w)
 
+
     }
 
 
@@ -718,9 +831,9 @@ namespace Opm {
     {
         // group control related processing
         if (wellCollection().groupControlActive()) {
-            for (int w = 0; w < numWells(); ++w) {
+            for (const auto& well : well_container_) {
-                WellControls* wc = well_container_[w]->wellControls();
+                WellControls* wc = well->wellControls();
-                WellNode& well_node = wellCollection().findWellNode(well_container_[w]->name());
+                WellNode& well_node = wellCollection().findWellNode(well->name());
 
                 // handling the situation that wells do not have a valid control
                 // it happens the well specified with GRUP and restarting due to non-convergencing
@@ -731,7 +844,7 @@ namespace Opm {
                 if (group_control_index >= 0 && ctrl_index < 0) {
                     // put well under group control
                     well_controls_set_current(wc, group_control_index);
-                    well_state_.currentControls()[w] = group_control_index;
+                    well_state_.currentControls()[well->indexOfWell()] = group_control_index;
                 }
 
                 // Final step, update whehter the well is under group control or individual control
@@ -804,15 +917,13 @@ namespace Opm {
             return;
         }
 
-        const int nw = numWells();
+        for (auto& well : well_container_) {
-
+            const std::string well_name = well->name();
-        for (int w = 0; w < nw; ++w) {
-            const std::string well_name = well_container_[w]->name();
             const WellNode& well_node = wellCollection().findWellNode(well_name);
 
             const double well_efficiency_factor = well_node.getAccumulativeEfficiencyFactor();
 
-            well_container_[w]->setWellEfficiencyFactor(well_efficiency_factor);
+            well->setWellEfficiencyFactor(well_efficiency_factor);
         }
     }
 
@@ -846,9 +957,10 @@ namespace Opm {
         std::vector<double> well_rates(np, 0.0);
         std::vector<double> convert_coeff(np, 1.0);
 
-        for (int w = 0; w < nw; ++w) {
+         for (auto& well : well_container_) {
-            const bool is_producer = well_container_[w]->wellType() == PRODUCER;
+            const bool is_producer = well->wellType() == PRODUCER;
-            const int well_cell_top = well_container_[w]->cells()[0];
+            const int well_cell_top =well->cells()[0];
+            const int w = well->indexOfWell();
             const int pvtRegionIdx = pvt_region_idx_[well_cell_top];
 
             // not sure necessary to change all the value to be positive
@@ -917,13 +1029,13 @@ namespace Opm {
     {
 
         if (wellCollection().groupControlActive()) {
-            for (int w = 0; w < numWells(); ++w) {
+           for (auto& well : well_container_) {
                 // update whether well is under group control
                 // get well node in the well collection
-                WellNode& well_node = wellCollection().findWellNode(well_container_[w]->name());
+                WellNode& well_node = wellCollection().findWellNode(well->name());
 
                 // update whehter the well is under group control or individual control
-                const int current = well_state_.currentControls()[w];
+                const int current = well_state_.currentControls()[well->indexOfWell()];
                 if (well_node.groupControlIndex() >= 0 && current == well_node.groupControlIndex()) {
                     // under group control
                     well_node.setIndividualControl(false);
@@ -938,8 +1050,8 @@ namespace Opm {
             // it will not change the control mode, only update the targets
             wellCollection().updateWellTargets(well_state_.wellRates());
 
-            for (int w = 0; w < numWells(); ++w) {
+            for (auto& well : well_container_) {
-                well_container_[w]->updateWellStateWithTarget(well_state_);
+                well->updateWellStateWithTarget(well_state_);
             }
         }
     }

opm/autodiff/StandardWell.hpp

@@ -188,6 +188,7 @@ namespace Opm
         using Base::wellHasTHPConstraints;
         using Base::mostStrictBhpFromBhpLimits;
         using Base::scalingFactor;
+        using Base::updateWellControl;
 
         // protected member variables from the Base class
         using Base::current_step_;
@@ -338,6 +339,9 @@ namespace Opm
 
         // handle the non reasonable fractions due to numerical overshoot
         void processFractions() const;
+
+        SimulatorReport solveWellEq(Simulator& ebosSimulator, WellState& well_state, const double dt, const std::vector<double>& B_avg, bool terminal_output);
+
     };
 
 }

opm/autodiff/StandardWell_impl.hpp

@@ -1570,6 +1570,7 @@ namespace Opm
             for ( int compIdx = 0; compIdx < num_components_; ++compIdx )
             {
                 report.converged = report.converged && (well_flux_residual[compIdx] < tol_wells) && control_eq_converged;
+                //std::cout << name() << " " << well_flux_residual[compIdx] << std::endl;
             }
         } else { // abnormal values found and no need to check the convergence
             report.converged = false;
@@ -2229,4 +2230,57 @@ namespace Opm
             }
         }
     }
+
+    template<typename TypeTag>
+    SimulatorReport
+    StandardWell<TypeTag>::solveWellEq(Simulator& ebosSimulator, WellState& well_state, const double dt, const std::vector<double>& B_avg, bool terminal_output)
+    {
+        const int max_iter = param_.max_welleq_iter_;
+        int it  = 0;
+        bool converged;
+        WellState well_state0 = well_state;
+        do {
+            assembleWellEq(ebosSimulator, dt, well_state, true);
+
+            ConvergenceReport report;
+            report = getWellConvergence(B_avg);
+            converged = report.converged;
+
+            if (converged) {
+                break;
+            }
+
+            ++it;
+            solveEqAndUpdateWellState(well_state);
+
+            wellhelpers::WellSwitchingLogger logger;
+            updateWellControl(well_state, logger);
+            initPrimaryVariablesEvaluation();
+            } while (it < max_iter);
+
+        if (converged) {
+            if ( terminal_output ) {
+                OpmLog::debug("Well equation solution gets converged with " + std::to_string(it) + " iterations");
+            }
+        } else {
+            if ( terminal_output ) {
+                OpmLog::debug("Well equation solution failed in getting converged with " + std::to_string(it) + " iterations");
+                well_state = well_state0;
+                updatePrimaryVariables(well_state);
+                // also recover the old well controls
+                //WellControls* wc = wellControls();
+                //well_controls_set_current(wc, well_state.currentControls()[indexOfWell()]);
+            }
+
+//#warning need the unconverged solution in the wtest. Either add a flag or always use the unconverged solution?
+
+
+        }
+
+        SimulatorReport report;
+        report.converged = converged;
+        report.total_well_iterations = it;
+        return report;
+
+    }
 }

opm/autodiff/WellInterface.hpp

@@ -28,10 +28,13 @@
 #include <opm/common/Exceptions.hpp>
 
 #include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/WellTestState.hpp>
+
 #include <opm/core/wells.h>
 #include <opm/core/well_controls.h>
 #include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/core/wells/WellsManager.hpp>
+#include <opm/core/simulator/SimulatorReport.hpp>
 
 #include <opm/autodiff/VFPProperties.hpp>
 #include <opm/autodiff/VFPInjProperties.hpp>
@@ -113,6 +116,9 @@ namespace Opm
         /// Well name.
         const std::string& name() const;
 
+        /// Index of well in the wells struct and wellState
+        const int indexOfWell() const;
+
         /// Well cells.
         const std::vector<int>& cells() {return well_cells_; }
 
@@ -172,7 +178,8 @@ namespace Opm
                                     bool only_wells) = 0;
 
         void updateListEconLimited(const WellState& well_state,
-                                   DynamicListEconLimited& list_econ_limited) const;
+                                   const double& simulationTime,
+                                   WellTestState& wellTestState) const;
 
         void setWellEfficiencyFactor(const double efficiency_factor);
 
@@ -216,6 +223,20 @@ namespace Opm
         // Add well contributions to matrix
         virtual void addWellContributions(Mat&) const
         {}
+
+        virtual SimulatorReport solveWellEq(Simulator& ebosSimulator, WellState& well_state, const double dt, const std::vector<double>& B_avg, bool terminal_output)
+        {
+#warning need to add this to multisegment wells
+            SimulatorReport report;
+            return report;
+        }
+
+        void closeWellsAndCompletions(WellTestState& wellTestState);
+
+#warning currently just return number of connections
+        int numberOfCompletions(){ return number_of_perforations_;}
+
+
     protected:
 
         // to indicate a invalid connection
@@ -332,6 +353,8 @@ namespace Opm
 
         double scalingFactor(const int comp_idx) const;
 
+
+
     };
 
 }

opm/autodiff/WellInterface_impl.hpp

@@ -1,6 +1,7 @@
 /*
   Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
   Copyright 2017 Statoil ASA.
+  Copyright 2018 IRIS
 
   This file is part of the Open Porous Media project (OPM).
 
@@ -101,8 +102,8 @@ namespace Opm
                       wells->sat_table_id + perf_index_end,
                       saturation_table_number_.begin() );
         }
-
         well_efficiency_factor_ = 1.0;
+
     }
 
 
@@ -169,6 +170,15 @@ namespace Opm
         return well_controls_;
     }
 
+    template<typename TypeTag>
+    const int
+    WellInterface<TypeTag>::
+    indexOfWell() const
+    {
+        return index_of_well_;
+    }
+
+
 
 
 
@@ -618,7 +628,8 @@ namespace Opm
     void
     WellInterface<TypeTag>::
     updateListEconLimited(const WellState& well_state,
-                          DynamicListEconLimited& list_econ_limited) const
+                          const double& simulationTime,
+                          WellTestState& wellTestState) const
     {
         // economic limits only apply for production wells.
         if (wellType() != PRODUCER) {
@@ -660,12 +671,11 @@ namespace Opm
                 OpmLog::warning("NOT_SUPPORTING_FOLLOWONWELL", "opening following on well after well closed is not supported yet");
             }
 
+            wellTestState.addClosedWell(well_name, WellTestConfig::Reason::ECONOMIC, simulationTime);
             if (well_ecl_->getAutomaticShutIn()) {
-                list_econ_limited.addShutWell(well_name);
+                const std::string msg = std::string("well ") + well_name + std::string(" will be shut due to rate economic limit");
-                const std::string msg = std::string("well ") + well_name + std::string(" will be shut in due to rate economic limit");
                 OpmLog::info(msg);
             } else {
-                list_econ_limited.addStoppedWell(well_name);
                 const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to rate economic limit");
                 OpmLog::info(msg);
             }
@@ -687,33 +697,43 @@ namespace Opm
             switch (workover) {
                 case WellEcon::CON:
                 {
-                    const bool last_connection = std::get<1>(ratio_check_return);
                     const int worst_offending_connection = std::get<2>(ratio_check_return);
 
                     assert((worst_offending_connection >= 0) && (worst_offending_connection < number_of_perforations_));
-
+#warning map to completions
-                    const int cell_worst_offending_connection = well_cells_[worst_offending_connection];
+                    wellTestState.addClosedCompletion(well_name, worst_offending_connection, simulationTime);
-                    list_econ_limited.addClosedConnectionsForWell(well_name, cell_worst_offending_connection);
                     const std::string msg = std::string("Connection ") + std::to_string(worst_offending_connection) + std::string(" for well ")
                                             + well_name + std::string(" will be closed due to economic limit");
                     OpmLog::info(msg);
 
-                    if (last_connection) {
+
-                        // TODO: there is more things to check here
+                    bool allCompletionsClosed = true;
-                        list_econ_limited.addShutWell(well_name);
+                    for (int perf = 0; perf < number_of_perforations_ ; ++perf) {
-                        const std::string msg2 = well_name + std::string(" will be shut due to the last connection closed");
+                        if (!wellTestState.hasCompletion(name(), perf)) {
-                        OpmLog::info(msg2);
+                            allCompletionsClosed = false;
+                        }
+                    }
+
+                    if (allCompletionsClosed) {
+                        wellTestState.addClosedWell(well_name, WellTestConfig::Reason::ECONOMIC, simulationTime);
+                        if (well_ecl_->getAutomaticShutIn()) {
+                            const std::string msg = well_name + std::string(" will be shut due to last compleation closed");
+                            OpmLog::info(msg);
+                        } else {
+                            const std::string msg = well_name + std::string(" will be stopped due to last compleation closed");
+                            OpmLog::info(msg);
+                        }
                     }
                     break;
                 }
                 case WellEcon::WELL:
                 {
+                    wellTestState.addClosedWell(well_name, WellTestConfig::Reason::ECONOMIC, 0);
                     if (well_ecl_->getAutomaticShutIn()) {
-                        list_econ_limited.addShutWell(well_name);
+                        // tell the controll that the well is closed
                         const std::string msg = well_name + std::string(" will be shut due to ratio economic limit");
                         OpmLog::info(msg);
                     } else {
-                        list_econ_limited.addStoppedWell(well_name);
                         const std::string msg = well_name + std::string(" will be stopped due to ratio economic limit");
                         OpmLog::info(msg);
                     }
@@ -866,4 +886,20 @@ namespace Opm
         }    
     }
 
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::closeWellsAndCompletions(WellTestState& wellTestState)
+    {
+        if (wellTestState.hasWell(name(), WellTestConfig::Reason::ECONOMIC)) {
+            assert(!well_ecl_->getAutomaticShutIn());
+            well_controls_stop_well(wellControls());
+        }
+
+        for (int perf = 0; perf < number_of_perforations_ ; ++perf) {
+            if (wellTestState.hasCompletion(name(), perf)) {
+                well_index_[perf] = 0.0;
+            }
+        }
+    }
+
 }

opm/core/wells/DynamicListEconLimited.hpp

@@ -80,6 +80,22 @@ namespace Opm
             }
         }
 
+        void removeShutWell(const std::string& well_name) {
+            auto itr = std::find(m_shut_wells.begin(), m_shut_wells.end(), well_name);
+            if (itr != m_shut_wells.end())
+                m_shut_wells.erase(itr);
+        }
+
+        void removeStoppedWell(const std::string& well_name) {
+            auto itr = std::find(m_stopped_wells.begin(), m_stopped_wells.end(), well_name);
+            if (itr != m_stopped_wells.end())
+                m_stopped_wells.erase(itr);
+        }
+
+        void removeClosedConnectionsForWell(const std::string& well_name) {
+             m_cells_closed_connections.erase(well_name);
+        }
+
     private:
         std::vector <std::string> m_shut_wells;
         std::vector <std::string> m_stopped_wells;