adding function updateWellStateWithTarget to StandardWell

without dealing with wsolvent function. It can be just a member variable
since we are handling well one by one individually.

opm/autodiff/StandardWell.hpp

@@ -137,6 +137,10 @@ namespace Opm
                              const BlackoilModelParameters& param,
                              WellState& well_state) const;
 
+        // TODO: later will check wheter we need current
+        void updateWellStateWithTarget(const int current,
+                                       WellState& xw) const;
+
         using WellInterface<TypeTag>::phaseUsage;
         using WellInterface<TypeTag>::active;
         using WellInterface<TypeTag>::numberOfPerforations;
@@ -157,15 +161,19 @@ namespace Opm
         using WellInterface<TypeTag>::numPhases;
         using WellInterface<TypeTag>::has_solvent;
         using WellInterface<TypeTag>::wellIndex;
+        using WellInterface<TypeTag>::wsolvent;
 
     protected:
 
         // TODO: maybe this function can go to some helper file.
         void localInvert(Mat& istlA) const;
 
+        // TODO: decide wether to use member function to refer to private member later
         using WellInterface<TypeTag>::vfp_properties_;
         using WellInterface<TypeTag>::gravity_;
         using WellInterface<TypeTag>::well_efficiency_factor_;
+        using WellInterface<TypeTag>::active_;
+        using WellInterface<TypeTag>::phase_usage_;
 
         // densities of the fluid in each perforation
         std::vector<double> perf_densities_;

opm/autodiff/StandardWell_impl.hpp

@@ -1104,4 +1104,226 @@ namespace Opm
         }
     }
 
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    updateWellStateWithTarget(const int current,
+                              WellState& xw) const
+    {
+        // number of phases
+        const int np = numberOfPhases();
+        const int well_index = indexOfWell();
+        const WellControls* wc = wellControls();
+        // 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;
+            // 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;
+
+            double aqua = 0.0;
+            double liquid = 0.0;
+            double vapour = 0.0;
+
+            const Opm::PhaseUsage& pu = phase_usage_;
+
+            if (active_[ Water ]) {
+                aqua = xw.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
+            }
+            if (active_[ Oil ]) {
+                 liquid = xw.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
+            }
+            if (active_[ Gas ]) {
+                vapour = xw.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
+            }
+
+            const int table_id = well_controls_iget_vfp(wc, current);
+            const double& thp    = well_controls_iget_target(wc, current);
+            const double& alq    = well_controls_iget_alq(wc, current);
+
+            //Set *BHP* target by calculating bhp from THP
+
+            // pick the density in the top layer
+            const double rho = perf_densities_[0];
+            const double well_ref_depth = perfDepth()[0];
+
+            // TODO: make the following a function and we call it so many times.
+            if (wellType() == INJECTOR) {
+
+                const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
+
+                const double dp = wellhelpers::computeHydrostaticCorrection(well_ref_depth, vfp_ref_depth, rho, gravity_);
+
+                xw.bhp()[well_index] = vfp_properties_->getInj()->bhp(table_id, aqua, liquid, vapour, thp) - dp;
+            }
+            else if (wellType() == PRODUCER) {
+                const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
+
+                const double dp = wellhelpers::computeHydrostaticCorrection(well_ref_depth, vfp_ref_depth, rho, gravity_);
+
+                xw.bhp()[well_index] = vfp_properties_->getProd()->bhp(table_id, aqua, liquid, vapour, thp, alq) - dp;
+            }
+            else {
+                OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
+            }
+            break;
+        }
+
+        case RESERVOIR_RATE: // intentional fall-through
+        case SURFACE_RATE:
+            // checking the number of the phases under control
+            int numPhasesWithTargetsUnderThisControl = 0;
+            for (int phase = 0; phase < np; ++phase) {
+                if (distr[phase] > 0.0) {
+                    numPhasesWithTargetsUnderThisControl += 1;
+                }
+            }
+
+            assert(numPhasesWithTargetsUnderThisControl > 0);
+
+            if (wellType() == INJECTOR) {
+                // assign target value as initial guess for injectors
+                // only handles single phase control at the moment
+                assert(numPhasesWithTargetsUnderThisControl == 1);
+
+                for (int phase = 0; phase < np; ++phase) {
+                    if (distr[phase] > 0.) {
+                        xw.wellRates()[np*well_index + phase] = target / distr[phase];
+                    } else {
+                        xw.wellRates()[np * well_index + phase] = 0.;
+                    }
+                }
+            } else if (wellType() == PRODUCER) {
+                // update the rates of phases under control based on the target,
+                // and also update rates of phases not under control to keep the rate ratio,
+                // assuming the mobility ratio does not change for the production wells
+                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];
+                    }
+                }
+
+                if (original_rates_under_phase_control != 0.0 ) {
+                    double scaling_factor = target / original_rates_under_phase_control;
+
+                    for (int phase = 0; phase < np; ++phase) {
+                        xw.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];
+                        } else {
+                            // this only happens for SURFACE_RATE control
+                            xw.wellRates()[np * well_index + phase] = target_rate_divided;
+                        }
+                    }
+                }
+            } else {
+                OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
+            }
+
+            break;
+        } // end of switch
+
+
+        std::vector<double> g = {1.0, 1.0, 0.01};
+        if (well_controls_iget_type(wc, current) == RESERVOIR_RATE) {
+            for (int phase = 0; phase < np; ++phase) {
+                g[phase] = distr[phase];
+            }
+        }
+
+        // the number of wells
+        const int nw = xw.bhp().size();
+
+        switch (well_controls_iget_type(wc, current)) {
+        case THP:
+        case BHP: {
+            xw.wellSolutions()[nw*XvarWell + well_index] = 0.0;
+            if (wellType() == INJECTOR) {
+                for (int p = 0; p < np; ++p) {
+                    xw.wellSolutions()[nw*XvarWell + well_index] += xw.wellRates()[np*well_index + p] * compFrac()[p];
+                }
+            } else {
+                for (int p = 0; p < np; ++p) {
+                    xw.wellSolutions()[nw*XvarWell + well_index] += g[p] * xw.wellRates()[np*well_index + p];
+                }
+            }
+            break;
+        }
+        case RESERVOIR_RATE: // Intentional fall-through
+        case SURFACE_RATE:
+            xw.wellSolutions()[nw*XvarWell + well_index] = xw.bhp()[well_index];
+            break;
+        } // end of switch
+
+        double tot_well_rate = 0.0;
+        for (int p = 0; p < np; ++p)  {
+            tot_well_rate += g[p] * xw.wellRates()[np*well_index + p];
+        }
+        if(std::abs(tot_well_rate) > 0) {
+            if (active_[ Water ]) {
+                xw.wellSolutions()[WFrac*nw + well_index] = g[Water] * xw.wellRates()[np*well_index + Water] / tot_well_rate;
+            }
+            if (active_[ Gas ]) {
+                xw.wellSolutions()[GFrac*nw + well_index] = g[Gas] * (1.0 - wsolvent()) * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
+            }
+            if (has_solvent) {
+                xw.wellSolutions()[SFrac*nw + well_index] = g[Gas] * wsolvent() * xw.wellRates()[np*well_index + Gas] / tot_well_rate ;
+            }
+        } else { // tot_well_rate == 0
+            if (wellType() == INJECTOR) {
+                // only single phase injection handled
+                if (active_[Water]) {
+                    if (distr[Water] > 0.0) {
+                        xw.wellSolutions()[WFrac * nw + well_index] = 1.0;
+                    } else {
+                        xw.wellSolutions()[WFrac * nw + well_index] = 0.0;
+                    }
+                }
+
+                if (active_[Gas]) {
+                    if (distr[Gas] > 0.0) {
+                        xw.wellSolutions()[GFrac * nw + well_index] = 1.0 - wsolvent();
+                        if (has_solvent) {
+                            xw.wellSolutions()[SFrac * nw + well_index] = wsolvent();
+                        }
+                    } else {
+                        xw.wellSolutions()[GFrac * nw + well_index] = 0.0;
+                    }
+                }
+
+                // TODO: it is possible to leave injector as a oil well,
+                // when F_w and F_g both equals to zero, not sure under what kind of circumstance
+                // this will happen.
+            } else if (wellType() == PRODUCER) { // producers
+                // TODO: the following are not addressed for the solvent case yet
+                if (active_[Water]) {
+                    xw.wellSolutions()[WFrac * nw + well_index] = 1.0 / np;
+                }
+                if (active_[Gas]) {
+                    xw.wellSolutions()[GFrac * nw + well_index] = 1.0 / np;
+                }
+            } else {
+                OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
+            }
+        }
+    }
+
+
 }

opm/autodiff/WellInterface.hpp

@@ -138,6 +138,8 @@ namespace Opm
         // TODO: for this kind of function, maybe can make a function with parameter perf
         const std::vector<int>& saturationTableNumber() const;
 
+        const double wsolvent() const;
+
     protected:
         // TODO: some variables shared by all the wells should be made static
         // well name

opm/autodiff/WellInterface_impl.hpp

@@ -352,4 +352,19 @@ namespace Opm
         return numComp;
     }
 
+
+
+
+    template<typename TypeTag>
+    const double
+    WellInterface<TypeTag>::
+    wsolvent() const
+    {
+        // TODO: not handling it for the moment
+        // TODO: it needs information from the well_ecl
+        // TODO: will decide on well_ecl role later.
+        // It can be just one member variable and no need to deal with well_ecl at all
+        return 0.0;
+    }
+
 }