removing the current argument in updateWellStateWithTarget

and some other cleaning up.

opm/autodiff/BlackoilWellModel.hpp

@@ -217,7 +217,6 @@ namespace Opm {
                                          std::vector<double>& voidage_conversion_coeffs) const;
 
             // Calculating well potentials for each well
-            // TODO: getBhp() will be refactored to reduce the duplication of the code calculating the bhp from THP.
             void computeWellPotentials(std::vector<double>& well_potentials);
 
             const std::vector<double>& wellPerfEfficiencyFactors() const;

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -704,7 +704,7 @@ namespace Opm {
             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(control, well_state_);
+            well_container_[w]->updateWellStateWithTarget(well_state_);
 
             // The wells are not considered to be newly added
             // for next time step
@@ -950,12 +950,7 @@ namespace Opm {
             wellCollection().updateWellTargets(well_state_.wellRates());
 
             for (int w = 0; w < numWells(); ++w) {
-                // TODO: check whether we need current argument in updateWellStateWithTarget
+                well_container_[w]->updateWellStateWithTarget(well_state_);
-                // maybe there is some circumstances that the current is different from the one
-                // in the WellState.
-                // while probalby, the current argument can be removed
-                const int current = well_state_.currentControls()[w];
-                well_container_[w]->updateWellStateWithTarget(current, well_state_);
             }
         }
     }

opm/autodiff/MultisegmentWell.hpp

@@ -121,8 +121,7 @@ namespace Opm
 
         /// updating the well state based the control mode specified with current
         // TODO: later will check wheter we need current
-        virtual void updateWellStateWithTarget(const int current,
+        virtual void updateWellStateWithTarget(WellState& well_state) const;
-                                               WellState& well_state) const;
 
         /// check whether the well equations get converged for this well
         virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const;

opm/autodiff/MultisegmentWell_impl.hpp

@@ -252,11 +252,11 @@ namespace Opm
     template <typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    updateWellStateWithTarget(const int current,
+    updateWellStateWithTarget(WellState& well_state) const
-                              WellState& well_state) const
     {
         // Updating well state bas on well control
         // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
+        const int current = well_state.currentControls()[index_of_well_];
         const double target = well_controls_iget_target(well_controls_, current);
         const double* distr = well_controls_iget_distr(well_controls_, current);
         switch (well_controls_iget_type(well_controls_, current)) {

opm/autodiff/StandardWell.hpp

@@ -132,8 +132,7 @@ namespace Opm
 
         /// updating the well state based the control mode specified with current
         // TODO: later will check wheter we need current
-        virtual void updateWellStateWithTarget(const int current,
+        virtual void updateWellStateWithTarget(WellState& well_state) const;
-                                               WellState& xw) const;
 
         /// check whether the well equations get converged for this well
         virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const;
@@ -252,7 +251,7 @@ namespace Opm
         // calculate the properties for the well connections
         // to calulate the pressure difference between well connections.
         void computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
-                                                         const WellState& xw,
+                                                         const WellState& well_state,
                                                          std::vector<double>& b_perf,
                                                          std::vector<double>& rsmax_perf,
                                                          std::vector<double>& rvmax_perf,
@@ -268,7 +267,7 @@ namespace Opm
 
         void computeConnectionPressureDelta();
 
-        void computeWellConnectionDensitesPressures(const WellState& xw,
+        void computeWellConnectionDensitesPressures(const WellState& well_state,
                                                     const std::vector<double>& b_perf,
                                                     const std::vector<double>& rsmax_perf,
                                                     const std::vector<double>& rvmax_perf,
@@ -278,7 +277,7 @@ namespace Opm
         void computeAccumWell();
 
         void computeWellConnectionPressures(const Simulator& ebosSimulator,
-                                                    const WellState& xw);
+                                                    const WellState& well_state);
 
         // TODO: to check whether all the paramters are required
         void computePerfRate(const IntensiveQuantities& intQuants,

opm/autodiff/StandardWell_impl.hpp

@@ -588,14 +588,14 @@ namespace Opm
             }
 
             if (has_polymer) {
-                EvalWell cq_s_poly = cq_s[Water];
+                // TODO: the application of well efficiency factor has not been tested with an example yet
+                EvalWell cq_s_poly = cq_s[Water] * well_efficiency_factor_;
                 if (well_type_ == INJECTOR) {
                     cq_s_poly *= wpolymer();
                 } else {
                     cq_s_poly *= extendEval(intQuants.polymerConcentration() * intQuants.polymerViscosityCorrection());
                 }
                 if (!only_wells) {
-                    // TODO: we need to consider the efficiency here.
                     for (int pvIdx = 0; pvIdx < numEq; ++pvIdx) {
                         ebosJac[cell_idx][cell_idx][contiPolymerEqIdx][pvIdx] -= cq_s_poly.derivative(pvIdx);
                     }
@@ -642,12 +642,12 @@ namespace Opm
             const int cell_idx = well_cells_[perf];
             const auto& intQuants = *(ebosSimulator.model().cachedIntensiveQuantities(cell_idx, /*timeIdx=*/0));
             const auto& fs = intQuants.fluidState();
-            EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
+            const EvalWell pressure = extendEval(fs.pressure(FluidSystem::oilPhaseIdx));
-            EvalWell bhp = getBhp();
+            const EvalWell bhp = getBhp();
 
             // Pressure drawdown (also used to determine direction of flow)
-            EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
+            const EvalWell well_pressure = bhp + perf_pressure_diffs_[perf];
-            EvalWell drawdown = pressure - well_pressure;
+            const EvalWell drawdown = pressure - well_pressure;
 
             if (drawdown.value() < 0 && well_type_ == INJECTOR)  {
                 return false;
@@ -976,40 +976,40 @@ namespace Opm
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
-    updateWellStateWithTarget(const int current,
+    updateWellStateWithTarget(WellState& well_state) const
-                              WellState& xw) const
     {
         // number of phases
         const int np = number_of_phases_;
         const int well_index = index_of_well_;
         const WellControls* wc = well_controls_;
+        const int current = well_state.currentControls()[well_index];
         // Updating well state and primary variables.
         // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE
         const double target = well_controls_iget_target(wc, current);
         const double* distr = well_controls_iget_distr(wc, current);
         switch (well_controls_iget_type(wc, current)) {
         case BHP:
-            xw.bhp()[well_index] = target;
+            well_state.bhp()[well_index] = target;
             // TODO: similar to the way below to handle THP
             // we should not something related to thp here when there is thp constraint
             break;
 
         case THP: {
-            xw.thp()[well_index] = target;
+            well_state.thp()[well_index] = target;
 
             const Opm::PhaseUsage& pu = phaseUsage();
             std::vector<double> rates(3, 0.0);
             if (active()[ Water ]) {
-                rates[ Water ] = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
+                rates[ Water ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
             }
             if (active()[ Oil ]) {
-                 rates[ Oil ] = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
+                 rates[ Oil ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
             }
             if (active()[ Gas ]) {
-                rates[ Gas ] = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
+                rates[ Gas ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
             }
 
-            xw.bhp()[well_index] = calculateBhpFromThp(rates, current);
+            well_state.bhp()[well_index] = calculateBhpFromThp(rates, current);
             break;
         }
 
@@ -1032,9 +1032,9 @@ namespace Opm
 
                 for (int phase = 0; phase < np; ++phase) {
                     if (distr[phase] > 0.) {
-                        xw.wellRates()[np*well_index + phase] = target / distr[phase];
+                        well_state.wellRates()[np*well_index + phase] = target / distr[phase];
                     } else {
-                        xw.wellRates()[np * well_index + phase] = 0.;
+                        well_state.wellRates()[np * well_index + phase] = 0.;
                     }
                 }
             } else if (well_type_ == PRODUCER) {
@@ -1044,7 +1044,7 @@ namespace Opm
                 double original_rates_under_phase_control = 0.0;
                 for (int phase = 0; phase < np; ++phase) {
                     if (distr[phase] > 0.0) {
-                        original_rates_under_phase_control += xw.wellRates()[np * well_index + phase] * distr[phase];
+                        original_rates_under_phase_control += well_state.wellRates()[np * well_index + phase] * distr[phase];
                     }
                 }
 
@@ -1052,17 +1052,17 @@ namespace Opm
                     double scaling_factor = target / original_rates_under_phase_control;
 
                     for (int phase = 0; phase < np; ++phase) {
-                        xw.wellRates()[np * well_index + phase] *= scaling_factor;
+                        well_state.wellRates()[np * well_index + phase] *= scaling_factor;
                     }
                 } else { // scaling factor is not well defied when original_rates_under_phase_control is zero
                     // separating targets equally between phases under control
                     const double target_rate_divided = target / numPhasesWithTargetsUnderThisControl;
                     for (int phase = 0; phase < np; ++phase) {
                         if (distr[phase] > 0.0) {
-                            xw.wellRates()[np * well_index + phase] = target_rate_divided / distr[phase];
+                            well_state.wellRates()[np * well_index + phase] = target_rate_divided / distr[phase];
                         } else {
                             // this only happens for SURFACE_RATE control
-                            xw.wellRates()[np * well_index + phase] = target_rate_divided;
+                            well_state.wellRates()[np * well_index + phase] = target_rate_divided;
                         }
                     }
                 }
@@ -1073,7 +1073,7 @@ namespace Opm
             break;
         } // end of switch
 
-        updatePrimaryVariables(xw);
+        updatePrimaryVariables(well_state);
     }
 
 
@@ -1084,7 +1084,7 @@ namespace Opm
     void
     StandardWell<TypeTag>::
     computePropertiesForWellConnectionPressures(const Simulator& ebosSimulator,
-                                                const WellState& xw,
+                                                const WellState& well_state,
                                                 std::vector<double>& b_perf,
                                                 std::vector<double>& rsmax_perf,
                                                 std::vector<double>& rvmax_perf,
@@ -1110,8 +1110,8 @@ namespace Opm
 
             // TODO: this is another place to show why WellState need to be a vector of WellState.
             // TODO: to check why should be perf - 1
-            const double p_above = perf == 0 ? xw.bhp()[w] : xw.perfPress()[first_perf_ + perf - 1];
+            const double p_above = perf == 0 ? well_state.bhp()[w] : well_state.perfPress()[first_perf_ + perf - 1];
-            const double p_avg = (xw.perfPress()[first_perf_ + perf] + p_above)/2;
+            const double p_avg = (well_state.perfPress()[first_perf_ + perf] + p_above)/2;
             const double temperature = fs.temperature(FluidSystem::oilPhaseIdx).value();
 
             if (pu.phase_used[Water]) {
@@ -1125,10 +1125,10 @@ namespace Opm
 
                 if (pu.phase_used[Oil]) {
                     const int oilpos_well = pu.phase_pos[Oil] + w * pu.num_phases;
-                    const double oilrate = std::abs(xw.wellRates()[oilpos_well]); //in order to handle negative rates in producers
+                    const double oilrate = std::abs(well_state.wellRates()[oilpos_well]); //in order to handle negative rates in producers
                     rvmax_perf[perf] = FluidSystem::gasPvt().saturatedOilVaporizationFactor(fs.pvtRegionIndex(), temperature, p_avg);
                     if (oilrate > 0) {
-                        const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w);
+                        const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - well_state.solventWellRate(w);
                         double rv = 0.0;
                         if (gasrate > 0) {
                             rv = oilrate / gasrate;
@@ -1152,9 +1152,9 @@ namespace Opm
                 if (pu.phase_used[Gas]) {
                     rsmax_perf[perf] = FluidSystem::oilPvt().saturatedGasDissolutionFactor(fs.pvtRegionIndex(), temperature, p_avg);
                     const int gaspos_well = pu.phase_pos[Gas] + w * pu.num_phases;
-                    const double gasrate = std::abs(xw.wellRates()[gaspos_well]) - xw.solventWellRate(w);
+                    const double gasrate = std::abs(well_state.wellRates()[gaspos_well]) - well_state.solventWellRate(w);
                     if (gasrate > 0) {
-                        const double oilrate = std::abs(xw.wellRates()[oilpos_well]);
+                        const double oilrate = std::abs(well_state.wellRates()[oilpos_well]);
                         double rs = 0.0;
                         if (oilrate > 0) {
                             rs = gasrate / oilrate;
@@ -1329,9 +1329,6 @@ namespace Opm
     StandardWell<TypeTag>::
     getWellConvergence(const std::vector<double>& B_avg) const
     {
-        typedef double Scalar;
-        typedef std::vector< Scalar > Vector;
-
         const int np = number_of_phases_;
 
         // the following implementation assume that the polymer is always after the w-o-g phases
@@ -1343,13 +1340,13 @@ namespace Opm
 
         // TODO: it should be the number of numWellEq
         // using num_components_ here for flow_ebos running 2p case.
-        std::vector<Scalar> res(num_components_);
+        std::vector<double> res(num_components_);
         for (int comp = 0; comp < num_components_; ++comp) {
             // magnitude of the residual matters
             res[comp] = std::abs(resWell_[0][comp]);
         }
 
-        Vector well_flux_residual(num_components_);
+        std::vector<double> well_flux_residual(num_components_);
 
         // Finish computation
         for ( int compIdx = 0; compIdx < num_components_; ++compIdx )
@@ -1396,7 +1393,7 @@ namespace Opm
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
-    computeWellConnectionDensitesPressures(const WellState& xw,
+    computeWellConnectionDensitesPressures(const WellState& well_state,
                                            const std::vector<double>& b_perf,
                                            const std::vector<double>& rsmax_perf,
                                            const std::vector<double>& rvmax_perf,
@@ -1409,10 +1406,10 @@ namespace Opm
 
         for (int perf = 0; perf < nperf; ++perf) {
             for (int phase = 0; phase < np; ++phase) {
-                perfRates[perf * num_components_ + phase] =  xw.perfPhaseRates()[(first_perf_ + perf) * np + phase];
+                perfRates[perf * num_components_ + phase] =  well_state.perfPhaseRates()[(first_perf_ + perf) * np + phase];
             }
             if(has_solvent) {
-                perfRates[perf * num_components_ + contiSolventEqIdx] =  xw.perfRateSolvent()[first_perf_ + perf];
+                perfRates[perf * num_components_ + contiSolventEqIdx] =  well_state.perfRateSolvent()[first_perf_ + perf];
             }
         }
 

opm/autodiff/WellInterface.hpp

@@ -193,16 +193,15 @@ namespace Opm
                                            const WellState& well_state,
                                            std::vector<double>& well_potentials) = 0;
 
-        virtual void updateWellStateWithTarget(const int current,
+        virtual void updateWellStateWithTarget(WellState& well_state) const = 0;
-                                               WellState& xw) const = 0;
 
-        void updateWellControl(WellState& xw,
+        void updateWellControl(WellState& well_state,
                                wellhelpers::WellSwitchingLogger& logger) const;
 
         virtual void updatePrimaryVariables(const WellState& well_state) const = 0;
 
         virtual void calculateExplicitQuantities(const Simulator& ebosSimulator,
-                                                 const WellState& xw) = 0; // should be const?
+                                                 const WellState& well_state) = 0; // should be const?
 
     protected:
 

opm/autodiff/WellInterface_impl.hpp

@@ -443,7 +443,7 @@ namespace Opm
         }
 
         if (updated_control_index != old_control_index) { //  || well_collection_->groupControlActive()) {
-            updateWellStateWithTarget(updated_control_index, well_state);
+            updateWellStateWithTarget(well_state);
         }
     }