using self-copying when calculating well potentials

in order to re-use the existing well model in stead of re-creating new
ones. The latter way really complicates the function createWellContainer
in BlackWellModel.

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1068,13 +1068,7 @@ namespace Opm {
         const int np = numPhases();
         well_potentials.resize(nw * np, 0.0);
 
-
-        const int reportStepIdx = ebosSimulator_.episodeIndex();
-        const double invalid_alq = -1e100;
-        const double invalid_vfp = -2147483647;
         auto well_state_copy = well_state_;
-        const Wells* local_wells = clone_wells(wells());
-        std::vector<WellInterfacePtr> well_container_copy = createWellContainer(reportStepIdx, local_wells, /*allow_closing_opening_wells=*/ false, deferred_logger);
 
         // average B factors are required for the convergence checking of well equations
         // Note: this must be done on all processes, even those with
@@ -1087,51 +1081,7 @@ namespace Opm {
         const bool write_restart_file = ebosSimulator_.vanguard().eclState().getRestartConfig().getWriteRestartFile(reportStepIdx);
         int exception_thrown = 0;
         try {
-            for (const auto& well : well_container_copy) {
+            for (const auto& well : well_container_) {
-                // Only compute the well potential when asked for
-                well->init(&phase_usage_, depth_, gravity_, number_of_cells_);
-
-                WellControls* wc = well->wellControls();
-                well_controls_clear(wc);
-                well_controls_assert_number_of_phases( wc , np);
-                if (well->wellType() == INJECTOR) {
-                    const auto controls = well->wellEcl()->injectionControls(summaryState);
-
-                    if (controls.hasControl(WellInjector::THP)) {
-                        const double thp_limit  = controls.thp_limit;
-                        const int    vfp_number = controls.vfp_table_number;
-                        well_controls_add_new(THP, thp_limit, invalid_alq, vfp_number, NULL, wc);
-                    }
-
-                    // we always have a bhp limit
-                    const double bhp_limit = controls.bhp_limit;
-                    well_controls_add_new(BHP, bhp_limit, invalid_alq, invalid_vfp, NULL, wc);
-                } else {
-                    const auto controls = well->wellEcl()->productionControls(summaryState);
-                    if (controls.hasControl(WellProducer::THP)) {
-                        const double thp_limit  = controls.thp_limit;
-                        const double alq_value  = controls.alq_value;
-                        const int    vfp_number = controls.vfp_table_number;
-                        well_controls_add_new(THP, thp_limit, alq_value, vfp_number, NULL, wc);
-                    }
-
-                    // we always have a bhp limit
-                    const double bhp_limit = controls.bhp_limit;
-                    well_controls_add_new(BHP, bhp_limit, invalid_alq, invalid_vfp, NULL, wc);
-
-                    well->setVFPProperties(vfp_properties_.get());
-
-                }
-
-                if (has_polymer_)
-                {
-                    const Grid& grid = ebosSimulator_.vanguard().grid();
-                    if (PolymerModule::hasPlyshlog() || GET_PROP_VALUE(TypeTag, EnablePolymerMW) ) {
-                        well->computeRepRadiusPerfLength(grid, cartesian_to_compressed_, deferred_logger);
-                    }
-                }
-
-
 
                 const bool needed_for_summary = ((summaryConfig.hasSummaryKey( "WWPI:" + well->name()) ||
                                                   summaryConfig.hasSummaryKey( "WOPI:" + well->name()) ||

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -548,23 +548,29 @@ namespace Opm
                           std::vector<double>& well_potentials,
                           Opm::DeferredLogger& deferred_logger)
     {
+        // creating a copy of the well itself, to avoid messing up the explicit informations
+        // during this copy, the only information not copied properly is the well controls
+        MultisegmentWell<TypeTag> well(*this);
+
+        well.well_controls_ = this->createWellControlsWithBHPAndTHP(deferred_logger);
+
         const int np = number_of_phases_;
         well_potentials.resize(np, 0.0);
 
-        updatePrimaryVariables(well_state, deferred_logger);
+        well.updatePrimaryVariables(well_state, deferred_logger);
 
         // initialize the primary variables in Evaluation, which is used in computePerfRate for computeWellPotentials
         // TODO: for computeWellPotentials, no derivative is required actually
-        initPrimaryVariablesEvaluation();
+        well.initPrimaryVariablesEvaluation();
 
         // get the bhp value based on the bhp constraints
-        const double bhp = Base::mostStrictBhpFromBhpLimits(deferred_logger);
+        const double bhp = well.mostStrictBhpFromBhpLimits(deferred_logger);
 
         // does the well have a THP related constraint?
-        if ( !Base::wellHasTHPConstraints() ) {
+        if ( !well.wellHasTHPConstraints() ) {
             assert(std::abs(bhp) != std::numeric_limits<double>::max());
 
-            computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
+            well.computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
         } else {
 
             const std::string msg = std::string("Well potential calculation is not supported for thp controlled multisegment wells \n")
@@ -573,9 +579,10 @@ namespace Opm
                     + "you will have to change the " + name() + " well to a standard well \n";
 
             deferred_logger.warning("WELL_POTENTIAL_FOR_THP_NOT_IMPLEMENTED_FOR_MULTISEG_WELLS", msg);
-            return;
         }
 
+        // destroy the newly created WellControls
+        well_controls_destroy(well.well_controls_);
     }
 
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -2430,23 +2430,30 @@ namespace Opm
                           std::vector<double>& well_potentials,
                           Opm::DeferredLogger& deferred_logger) // const
     {
-        updatePrimaryVariables(well_state, deferred_logger);
+
-        computeWellConnectionPressures(ebosSimulator, well_state);
+        // creating a copy of the well itself, to avoid messing up the explicit informations
+        // during this copy, the only information not copied properly is the well controls
+        StandardWell<TypeTag> well(*this);
+
+        well.well_controls_ = this->createWellControlsWithBHPAndTHP(deferred_logger);
+
+        well.updatePrimaryVariables(well_state, deferred_logger);
+        well.computeWellConnectionPressures(ebosSimulator, well_state);
 
         // initialize the primary variables in Evaluation, which is used in computePerfRate for computeWellPotentials
         // TODO: for computeWellPotentials, no derivative is required actually
-        initPrimaryVariablesEvaluation();
+        well.initPrimaryVariablesEvaluation();
 
         const int np = number_of_phases_;
         well_potentials.resize(np, 0.0);
 
         // get the bhp value based on the bhp constraints
-        const double bhp = mostStrictBhpFromBhpLimits(deferred_logger);
+        const double bhp = well.mostStrictBhpFromBhpLimits(deferred_logger);
 
         // does the well have a THP related constraint?
-        if ( !wellHasTHPConstraints() ) {
+        if ( !well.wellHasTHPConstraints() ) {
             assert(std::abs(bhp) != std::numeric_limits<double>::max());
-            computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
+            well.computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
         } else {
             // the well has a THP related constraint
             // checking whether a well is newly added, it only happens at the beginning of the report step
@@ -2458,7 +2465,7 @@ namespace Opm
                 }
             } else {
                 // We need to generate a reasonable rates to start the iteration process
-                computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
+                well.computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
                 for (double& value : well_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.
@@ -2467,8 +2474,12 @@ namespace Opm
                 }
             }
 
-            well_potentials = computeWellPotentialWithTHP(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
+            well_potentials = well.computeWellPotentialWithTHP(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
         }
+
+
+        // destroy the newly created WellControls
+        well_controls_destroy(well.well_controls_);
     }
 
 

opm/simulators/wells/WellInterface.hpp

@@ -447,6 +447,8 @@ namespace Opm
 
         void initCompletions();
 
+        WellControls* createWellControlsWithBHPAndTHP(DeferredLogger& deferred_logger) const;
+
         // count the number of times an output log message is created in the productivity
         // index calculations
         int well_productivity_index_logger_counter_;

opm/simulators/wells/WellInterface_impl.hpp

@@ -1442,4 +1442,35 @@ namespace Opm
     }
 
 
+
+
+
+
+    template<typename TypeTag>
+    WellControls*
+    WellInterface<TypeTag>::
+    createWellControlsWithBHPAndTHP(DeferredLogger& deferred_logger) const
+    {
+        WellControls* wc = well_controls_create();
+        well_controls_assert_number_of_phases(wc, number_of_phases_);
+
+        // a well always has a bhp limit
+        const double invalid_alq = -1e100;
+        const double invalid_vfp = -2147483647;
+        const double bhp_limit = this->mostStrictBhpFromBhpLimits(deferred_logger);
+        well_controls_add_new(BHP, bhp_limit, invalid_alq, invalid_vfp, NULL, wc);
+
+        if (this->wellHasTHPConstraints()) {
+            // it might be better to do it through EclipseState?
+            const double thp_limit = this->getTHPConstraint(deferred_logger);
+            const double thp_control_index = this->getControlIndex(THP);
+            const  int vfp_number = well_controls_iget_vfp(well_controls_, thp_control_index);
+            const double alq = well_controls_iget_alq(well_controls_, thp_control_index);
+            well_controls_add_new(THP, thp_limit, alq, vfp_number, NULL, wc);
+        }
+
+        return wc;
+    }
+
+
 }