Merge pull request #1609 from GitPaean/refactoring_well_test

Refactoring well test to prepare for Physical limit related well test.

opm/autodiff/BlackoilWellModel.hpp

@@ -237,6 +237,8 @@ namespace Opm {
             // create the well container
             std::vector<WellInterfacePtr > createWellContainer(const int time_step);
 
+            WellInterfacePtr createWellForWellTest(const std::string& well_name, const int report_step) const;
+
             WellState well_state_;
             WellState previous_well_state_;
 

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -158,7 +158,7 @@ namespace Opm {
             well->setVFPProperties(vfp_properties_.get());
         }
 
-        // Close wells and completions due to economical reasons
+        // Close completions due to economical reasons
         for (auto& well : well_container_) {
             well->closeCompletions(wellTestState_);
         }
@@ -170,102 +170,38 @@ namespace Opm {
     void
     BlackoilWellModel<TypeTag>::wellTesting(const int timeStepIdx, const double simulationTime) {
         const auto& wtest_config = schedule().wtestConfig(timeStepIdx);
+        if (wtest_config.size() == 0) { // there is no WTEST request
+            return;
+        }
+
         const auto& wellsForTesting = wellTestState_.updateWell(wtest_config, simulationTime);
+        if (wellsForTesting.size() == 0) { // there is no well available for WTEST at the moment
+            return;
+        }
 
-        // Do the well testing if enabled
+        // average B factors are required for the convergence checking of well equations
-        if (wtest_config.size() > 0 && wellsForTesting.size() > 0) {
+        std::vector< Scalar > B_avg(numComponents(), Scalar() );
-            // solve the well equation isolated from the reservoir.
+        computeAverageFormationFactor(B_avg);
-            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 (const auto& testWell : wellsForTesting) {
-            for (auto& testWell : wellsForTesting) {
+            const std::string& well_name = testWell.first;
-                const std::string msg = std::string("well ") + testWell.first + std::string(" is tested");
+            const std::string msg = std::string("well ") + well_name + std::string(" is tested");
-                OpmLog::info(msg);
+            OpmLog::info(msg);
 
-                // Finding the location of the well in wells_ecl
+            // this is the well we will test
-                const int nw_wells_ecl = wells_ecl_.size();
+            WellInterfacePtr well = createWellForWellTest(well_name, timeStepIdx);
-                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];
 
-                // Finding the location of the well in wells struct.
+            // some preparation before the well can be used
-                const int nw = numWells();
+            well->init(&phase_usage_, depth_, gravity_, number_of_cells_);
-                int wellidx = -999;
+            const WellNode& well_node = wellCollection().findWellNode(well_name);
-                for (int w = 0; w < nw; ++w) {
+            const double well_efficiency_factor = well_node.getAccumulativeEfficiencyFactor();
-                    if (testWell.first == std::string(wells()->name[w])) {
+            well->setWellEfficiencyFactor(well_efficiency_factor);
-                        wellidx = w;
+            well->setVFPProperties(vfp_properties_.get());
-                        break;
-                    }
-                }                
-                if (wellidx < 0) {
-                    OPM_THROW(std::logic_error, "Could not find the well  " << testWell.first << " in the well struct ");
-                }
 
-                // Use the pvtRegionIdx from the top cell
+            const WellTestConfig::Reason testing_reason = testWell.second;
-                const int well_cell_top = wells()->well_cells[wells()->well_connpos[wellidx]];
-                const int pvtreg = pvt_region_idx_[well_cell_top];
 
-                if ( !well_ecl->isMultiSegment(timeStepIdx) || !param_.use_multisegment_well_) {
+            well->wellTesting(ebosSimulator_, B_avg, simulationTime, timeStepIdx, terminal_output_,
-                    well_container.emplace_back(new StandardWell<TypeTag>(well_ecl, timeStepIdx, wells(),
+                              testing_reason, well_state_, wellTestState_);
-                                                                          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;
-                WellState wellStateCopy = well_state_;
-                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->updatePrimaryVariables(wellStateCopy);
-                well->initPrimaryVariablesEvaluation();
-
-                bool testWell = true;
-		// if a well is closed because all completions are closed, we need to check each completion
-		// individually. We first open all completions, then we close one by one by calling updateWellTestState
-		// untill the number of closed completions do not increase anymore. 
-                while (testWell) {
-                    const size_t numberOfClosedCompletions = wellTestStateForTheWellTest.sizeCompletions();
-                    well->solveWellForTesting(ebosSimulator_, wellStateCopy, B_avg, terminal_output_);
-                    well->updateWellTestState(wellStateCopy, simulationTime, wellTestStateForTheWellTest, /*writeMessageToOPMLog=*/ false);
-                    well->closeCompletions(wellTestStateForTheWellTest);
-
-                    // Stop testing if the well is closed or shut due to all completions shut
-                    // Also check if number of completions has increased. If the number of closed completions do not increased
-                    // we stop the testing.  
-                    if (wellTestStateForTheWellTest.sizeWells() > 0 || numberOfClosedCompletions == wellTestStateForTheWellTest.sizeCompletions())
-                        testWell = false;
-                }
-
-                // 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 (auto& completion : well->wellEcl()->getCompletions(timeStepIdx)) {
-                        if (!wellTestStateForTheWellTest.hasCompletion(well->name(), completion.first))
-                            wellTestState_.dropCompletion(well->name(), completion.first);
-                    }
-                }
-            }
         }
     }
 
@@ -373,6 +309,58 @@ namespace Opm {
 
 
 
+    template<typename TypeTag>
+    typename BlackoilWellModel<TypeTag>::WellInterfacePtr
+    BlackoilWellModel<TypeTag>::
+    createWellForWellTest(const std::string& well_name,
+                          const int report_step) const
+    {
+        // Finding the location of the well in wells_ecl
+        const int nw_wells_ecl = wells_ecl_.size();
+        int index_well_ecl = 0;
+        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];
+
+        // Finding the location of the well in wells struct.
+        const int nw = numWells();
+        int well_index_wells = -999;
+        for (int w = 0; w < nw; ++w) {
+            if (well_name == std::string(wells()->name[w])) {
+                well_index_wells = w;
+                break;
+            }
+        }
+
+        if (well_index_wells < 0) {
+            OPM_THROW(std::logic_error, "Could not find the well  " << well_name << " in the well struct ");
+        }
+
+        // Use the pvtRegionIdx from the top cell
+        const int well_cell_top = wells()->well_cells[wells()->well_connpos[well_index_wells]];
+        const int pvtreg = pvt_region_idx_[well_cell_top];
+
+        if ( !well_ecl->isMultiSegment(report_step) || !param_.use_multisegment_well_) {
+             return WellInterfacePtr(new StandardWell<TypeTag>(well_ecl, report_step, wells(),
+                                                 param_, *rateConverter_, pvtreg, numComponents() ) );
+        } else {
+             return WellInterfacePtr(new MultisegmentWell<TypeTag>(well_ecl, report_step, wells(),
+                                                 param_, *rateConverter_, pvtreg, numComponents() ) );
+        }
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
@@ -598,7 +586,6 @@ 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
@@ -767,7 +754,7 @@ namespace Opm {
     updateWellTestState(const double& simulationTime, WellTestState& wellTestState) const
     {
         for (const auto& well : well_container_) {
-            well->updateWellTestState(well_state_, simulationTime, wellTestState, /*writeMessageToOPMLog=*/ true);
+            well->updateWellTestState(well_state_, simulationTime, /*writeMessageToOPMLog=*/ true, wellTestState);
         }
     }
 

opm/autodiff/WellInterface.hpp

@@ -149,8 +149,9 @@ namespace Opm
 
         void updateWellTestState(const WellState& well_state,
                                  const double& simulationTime,
-                                 WellTestState& wellTestState,
+                                 const bool& writeMessageToOPMLog,
-                                 const bool& writeMessageToOPMLog) const;
+                                 WellTestState& wellTestState) const;
+
 
         void setWellEfficiencyFactor(const double efficiency_factor);
 
@@ -195,12 +196,16 @@ namespace Opm
         virtual void addWellContributions(Mat&) const
         {}
 
-        void solveWellForTesting(Simulator& ebosSimulator, WellState& well_state, const std::vector<double>& B_avg, bool terminal_output);
-
         void closeCompletions(WellTestState& wellTestState);
 
         const Well* wellEcl() const;
 
+        // TODO: theoretically, it should be a const function
+        // Simulator is not const is because that assembleWellEq is non-const Simulator
+        void wellTesting(Simulator& simulator, const std::vector<double>& B_avg,
+                         const double simulation_time, const int report_step,  const bool terminal_output,
+                         const WellTestConfig::Reason testing_reason, const WellState& well_state,
+                         WellTestState& welltest_state);
 
     protected:
 
@@ -317,6 +322,16 @@ namespace Opm
 
         double scalingFactor(const int comp_idx) const;
 
+        void wellTestingEconomic(Simulator& simulator, const std::vector<double>& B_avg,
+                                 const double simulation_time, const int report_step, const bool terminal_output,
+                                 const WellState& well_state, WellTestState& welltest_state);
+
+        void updateWellTestStateEconomic(const WellState& well_state,
+                                         const double simulation_time,
+                                         const bool write_message_to_opmlog,
+                                         WellTestState& well_test_state) const;
+
+        void solveWellForTesting(Simulator& ebosSimulator, WellState& well_state, const std::vector<double>& B_avg, bool terminal_output);
 
     };
 

opm/autodiff/WellInterface_impl.hpp

@@ -641,16 +641,32 @@ namespace Opm
     WellInterface<TypeTag>::
     updateWellTestState(const WellState& well_state,
                         const double& simulationTime,
-                        WellTestState& wellTestState,
+                        const bool& writeMessageToOPMLog,
-                        const bool& writeMessageToOPMLog) const
+                        WellTestState& wellTestState) const
     {
-        // economic limits only apply for production wells.
+        // currently, we only updateWellTestState for producers
         if (wellType() != PRODUCER) {
             return;
         }
 
-        // flag to check if the mim oil/gas rate limit is violated
+        // updating well test state based on Economic limits.
-        bool rate_limit_violated = false;
+        updateWellTestStateEconomic(well_state, simulationTime, writeMessageToOPMLog, wellTestState);
+
+        // TODO: well can be shut/closed due to other reasons
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    updateWellTestStateEconomic(const WellState& well_state,
+                                const double simulation_time,
+                                const bool write_message_to_opmlog,
+                                WellTestState& well_test_state) const
+    {
         const WellEconProductionLimits& econ_production_limits = well_ecl_->getEconProductionLimits(current_step_);
 
         // if no limit is effective here, then continue to the next well
@@ -658,13 +674,14 @@ namespace Opm
             return;
         }
 
-        const std::string well_name = name();
+        // flag to check if the mim oil/gas rate limit is violated
+        bool rate_limit_violated = false;
 
         // for the moment, we only handle rate limits, not handling potential limits
         // the potential limits should not be difficult to add
         const WellEcon::QuantityLimitEnum& quantity_limit = econ_production_limits.quantityLimit();
         if (quantity_limit == WellEcon::POTN) {
-            const std::string msg = std::string("POTN limit for well ") + well_name + std::string(" is not supported for the moment. \n")
+            const std::string msg = std::string("POTN limit for well ") + name() + std::string(" is not supported for the moment. \n")
                                   + std::string("All the limits will be evaluated based on RATE. ");
             OpmLog::warning("NOT_SUPPORTING_POTN", msg);
         }
@@ -675,8 +692,10 @@ namespace Opm
 
         if (rate_limit_violated) {
             if (econ_production_limits.endRun()) {
-                const std::string warning_message = std::string("ending run after well closed due to economic limits is not supported yet \n")
+                const std::string warning_message = std::string("ending run after well closed due to economic limits")
-                                                  + std::string("the program will keep running after ") + well_name + std::string(" is closed");
+                                                  + std::string("is not supported yet \n")
+                                                  + std::string("the program will keep running after ") + name()
+                                                  + std::string(" is closed");
                 OpmLog::warning("NOT_SUPPORTING_ENDRUN", warning_message);
             }
 
@@ -684,13 +703,13 @@ 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);
+            well_test_state.addClosedWell(name(), WellTestConfig::Reason::ECONOMIC, simulation_time);
-            if (writeMessageToOPMLog) {
+            if (write_message_to_opmlog) {
                 if (well_ecl_->getAutomaticShutIn()) {
-                    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 ") + name() + std::string(" will be shut due to rate economic limit");
                     OpmLog::info(msg);
                 } else {
-                    const std::string msg = std::string("well ") + well_name + std::string(" will be stopped due to rate economic limit");
+                    const std::string msg = std::string("well ") + name() + std::string(" will be stopped due to rate economic limit");
                     OpmLog::info(msg);
                 }
             }
@@ -714,15 +733,15 @@ namespace Opm
                 {
                     const int worst_offending_completion = std::get<1>(ratio_check_return);
 
-                    wellTestState.addClosedCompletion(well_name, worst_offending_completion, simulationTime);
+                    well_test_state.addClosedCompletion(name(), worst_offending_completion, simulation_time);
-                    if (writeMessageToOPMLog) {
+                    if (write_message_to_opmlog) {
                         if (worst_offending_completion < 0) {
-                            const std::string msg = std::string("Connection ") + std::to_string(- worst_offending_completion) + std::string(" for well ")
+                            const std::string msg = std::string("Connection ") + std::to_string(- worst_offending_completion)
-                                    + well_name + std::string(" will be closed due to economic limit");
+                                    + std::string(" for well ") + name() + std::string(" will be closed due to economic limit");
                             OpmLog::info(msg);
                         } else {
-                            const std::string msg = std::string("Completion ") + std::to_string(worst_offending_completion) + std::string(" for well ")
+                            const std::string msg = std::string("Completion ") + std::to_string(worst_offending_completion)
-                                    + well_name + std::string(" will be closed due to economic limit");
+                                    + std::string(" for well ") + name() + std::string(" will be closed due to economic limit");
                             OpmLog::info(msg);
                         }
                     }
@@ -730,19 +749,19 @@ namespace Opm
                     bool allCompletionsClosed = true;
                     const auto& connections = well_ecl_->getConnections(current_step_);
                     for (const auto& connection : connections) {
-                        if (!wellTestState.hasCompletion(name(), connection.complnum())) {
+                        if (!well_test_state.hasCompletion(name(), connection.complnum())) {
                             allCompletionsClosed = false;
                         }
                     }
 
                     if (allCompletionsClosed) {
-                        wellTestState.addClosedWell(well_name, WellTestConfig::Reason::ECONOMIC, simulationTime);
+                        well_test_state.addClosedWell(name(), WellTestConfig::Reason::ECONOMIC, simulation_time);
-                        if (writeMessageToOPMLog) {
+                        if (write_message_to_opmlog) {
                             if (well_ecl_->getAutomaticShutIn()) {
-                            	const std::string msg = well_name + std::string(" will be shut due to last completion closed");
+                                const std::string msg = name() + std::string(" will be shut due to last completion closed");
                             	OpmLog::info(msg);
                             } else {
-                            	const std::string msg = well_name + std::string(" will be stopped due to last completion closed");
+                                const std::string msg = name() + std::string(" will be stopped due to last completion closed");
                                 OpmLog::info(msg);
                             }
                         }
@@ -751,14 +770,14 @@ namespace Opm
                 }
                 case WellEcon::WELL:
                 {
-                wellTestState.addClosedWell(well_name, WellTestConfig::Reason::ECONOMIC, 0);
+                well_test_state.addClosedWell(name(), WellTestConfig::Reason::ECONOMIC, 0);
-                if (writeMessageToOPMLog) {
+                if (write_message_to_opmlog) {
                     if (well_ecl_->getAutomaticShutIn()) {
                         // tell the controll that the well is closed
-                        const std::string msg = well_name + std::string(" will be shut due to ratio economic limit");
+                        const std::string msg = name() + std::string(" will be shut due to ratio economic limit");
                         OpmLog::info(msg);
                     } else {
-                        const std::string msg = well_name + std::string(" will be stopped due to ratio economic limit");
+                        const std::string msg = name() + std::string(" will be stopped due to ratio economic limit");
                         OpmLog::info(msg);
                     }
                 }
@@ -768,7 +787,80 @@ namespace Opm
                     break;
                 default:
                 {
-                    OpmLog::warning("NOT_SUPPORTED_WORKOVER_TYPE", "not supporting workover type " + WellEcon::WorkoverEnumToString(workover) );
+                    OpmLog::warning("NOT_SUPPORTED_WORKOVER_TYPE",
+                                    "not supporting workover type " + WellEcon::WorkoverEnumToString(workover) );
+                }
+            }
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    wellTesting(Simulator& simulator, const std::vector<double>& B_avg,
+                const double simulation_time, const int report_step, const bool terminal_output,
+                const WellTestConfig::Reason testing_reason, const WellState& well_state,
+                WellTestState& well_test_state)
+    {
+        if (testing_reason == WellTestConfig::Reason::ECONOMIC) {
+            wellTestingEconomic(simulator, B_avg, simulation_time, report_step,
+                                terminal_output, well_state, well_test_state);
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    wellTestingEconomic(Simulator& simulator, const std::vector<double>& B_avg,
+                        const double simulation_time, const int report_step, const bool terminal_output,
+                        const WellState& well_state, WellTestState& welltest_state)
+    {
+        WellState well_state_copy = well_state;
+
+        updatePrimaryVariables(well_state_copy);
+        initPrimaryVariablesEvaluation();
+
+        // create a well
+        WellTestState welltest_state_temp;
+
+        bool testWell = true;
+        // if a well is closed because all completions are closed, we need to check each completion
+        // individually. We first open all completions, then we close one by one by calling updateWellTestState
+        // untill the number of closed completions do not increase anymore.
+        while (testWell) {
+            const size_t original_number_closed_completions = welltest_state_temp.sizeCompletions();
+            solveWellForTesting(simulator, well_state_copy, B_avg, terminal_output);
+            updateWellTestState(well_state_copy, simulation_time, /*writeMessageToOPMLog=*/ false, welltest_state_temp);
+            closeCompletions(welltest_state_temp);
+
+            // Stop testing if the well is closed or shut due to all completions shut
+            // Also check if number of completions has increased. If the number of closed completions do not increased
+            // we stop the testing.
+            // TODO: it can be tricky here, if the well is shut/closed due to other reasons
+            if ( welltest_state_temp.sizeWells() > 0 ||
+                (original_number_closed_completions == welltest_state_temp.sizeCompletions()) ) {
+                    testWell = false; // this terminates the while loop
+            }
+        }
+
+        // update wellTestState if the well test succeeds
+        if (!welltest_state_temp.hasWell(name(), WellTestConfig::Reason::ECONOMIC)) {
+            welltest_state.openWell(name());
+            const std::string msg = std::string("well ") + name() + std::string(" is re-opened");
+            OpmLog::info(msg);
+
+            // also reopen completions
+            for (auto& completion : well_ecl_->getCompletions(report_step)) {
+                if (!welltest_state_temp.hasCompletion(name(), completion.first)) {
+                    welltest_state.dropCompletion(name(), completion.first);
                 }
             }
         }