adding function updateWellStateWithTHPTargetIPR()

this function can update the well state related based on the inflow
perfomance relationship and THP target.

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -847,7 +847,7 @@ namespace Opm {
         wellhelpers::WellSwitchingLogger logger;
 
         for (const auto& well : well_container_) {
-            well->updateWellControl(well_state_, logger);
+            well->updateWellControl(ebosSimulator_, well_state_, logger);
         }
 
         updateGroupControls();
@@ -938,7 +938,7 @@ namespace Opm {
             well_state_.currentControls()[w] = control;
 
             if (well_state_.effectiveEventsOccurred(w) ) {
-                well->updateWellStateWithTarget(well_state_);
+                well->updateWellStateWithTarget(ebosSimulator_, well_state_);
             }
 
             // there is no new well control change input within a report step,
@@ -1186,7 +1186,7 @@ namespace Opm {
 
             // TODO: we should only do the well is involved in the update group targets
             for (auto& well : well_container_) {
-                well->updateWellStateWithTarget(well_state_);
+                well->updateWellStateWithTarget(ebosSimulator_, well_state_);
                 well->updatePrimaryVariables(well_state_);
             }
         }

opm/autodiff/MultisegmentWell.hpp

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

opm/autodiff/MultisegmentWell_impl.hpp

@@ -254,7 +254,8 @@ namespace Opm
     template <typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    updateWellStateWithTarget(WellState& well_state) const
+    updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
+                              WellState& well_state) /* const */
     {
         // Updating well state bas on well control
         // Target values are used as initial conditions for BHP, THP, and SURFACE_RATE

opm/autodiff/StandardWell.hpp

@@ -140,9 +140,8 @@ namespace Opm
                                     const double dt,
                                     WellState& well_state) override;
 
-        /// updating the well state based the control mode specified with current
+        virtual void updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
-        // TODO: later will check wheter we need current
+                                               WellState& well_state) /* const */ override;
-        virtual void updateWellStateWithTarget(WellState& well_state) const override;
 
         /// check whether the well equations get converged for this well
         virtual ConvergenceReport getWellConvergence(const std::vector<double>& B_avg) const override;
@@ -364,6 +363,18 @@ namespace Opm
         // updating the inflow based on the current reservoir condition
         void updateIPR(const Simulator& ebos_simulator) const;
 
+        // update WellState based on IPR and associated VFP table
+        void updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator,
+                                             WellState& well_state) const;
+
+        void updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator,
+                                                     WellState& well_state) const;
+
+        // calculate the BHP from THP target based on IPR
+        // TODO: we need to check the operablility here first, if not operable, then maybe there is
+        // no point to do this
+        double calculateBHPWithTHPTargetIPR() const;
+
         // relaxation factor considering only one fraction value
         static double relaxationFactorFraction(const double old_value,
                                                const double dx);

opm/autodiff/StandardWell_impl.hpp

@@ -1125,7 +1125,8 @@ namespace Opm
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
-    updateWellStateWithTarget(WellState& well_state) const
+    updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
+                              WellState& well_state) /* const */
     {
         // number of phases
         const int np = number_of_phases_;
@@ -1142,32 +1143,23 @@ namespace Opm
             // TODO: similar to the way below to handle THP
             // we should not something related to thp here when there is thp constraint
             // or when can calculate the THP (table avaiable or requested for output?)
+            // TODO: we should address this in a function updateWellStateWithBHPTarget.
+            // TODO: however, the reason that this one minght not be that critical with
+            // TODO: the effects remaining to be investigated.
+
             break;
-
         case THP: {
-            // TODO: this will be the big task here.
+            // TODO: adding the checking for the operability
-            // p_bhp = BHP(THP, rates(p_bhp))
+            // TODO: should we do updateIPR before this or within the related functions
-            // more sophiscated techniques is required to obtain the bhp and rates here
+            updateIPR(ebos_simulator);
-            well_state.thp()[well_index] = target;
+            updateWellStateWithTHPTargetIPR(ebos_simulator, well_state);
-
-            const Opm::PhaseUsage& pu = phaseUsage();
-            std::vector<double> rates(3, 0.0);
-            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-                rates[ Water ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Water ] ];
-            }
-            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
-                 rates[ Oil ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Oil ] ];
-            }
-            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                rates[ Gas ] = well_state.wellRates()[well_index*np + pu.phase_pos[ Gas ] ];
-            }
-
-            well_state.bhp()[well_index] = calculateBhpFromThp(rates, current);
             break;
         }
 
         case RESERVOIR_RATE: // intentional fall-through
         case SURFACE_RATE:
+            // TODO: something needs to be done with BHP and THP here
+            // TODO: they should go to a separate function
             // checking the number of the phases under control
             int numPhasesWithTargetsUnderThisControl = 0;
             for (int phase = 0; phase < np; ++phase) {
@@ -1326,6 +1318,101 @@ namespace Opm
 
 
 
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    updateWellStateWithTHPTargetIPR(const Simulator& ebos_simulator,
+                                    WellState& well_state) const
+    {
+        if (well_type_ == PRODUCER) {
+            updateWellStateWithTHPTargetIPRProducer(ebos_simulator,
+                                                    well_state);
+        }
+
+        if (well_type_ == INJECTOR) {
+            well_state.thp()[index_of_well_] = this->getTHPConstraint();
+            // TODO: more work needs to be done for the injectors here, while injectors
+            // have been okay with the current strategy relying on well control equation directly.
+        }
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    updateWellStateWithTHPTargetIPRProducer(const Simulator& ebos_simulator,
+                                            WellState& well_state) const
+    {
+
+        well_state.thp()[index_of_well_] = this->getTHPConstraint();
+
+        const double bhp = calculateBHPWithTHPTargetIPR();
+
+        well_state.bhp()[index_of_well_] = bhp;
+
+        // TODO: explicit quantities are always tricky for this type of situation
+        updatePrimaryVariables(well_state);
+        initPrimaryVariablesEvaluation();
+
+        std::vector<double> rates;
+        computeWellRatesWithBhp(ebos_simulator, bhp, rates);
+
+        // TODO: double checke the obtained rates
+        // this is another places we might obtain negative rates
+
+        for (size_t p = 0; p < number_of_phases_; ++p) {
+            well_state.wellRates()[number_of_phases_ * index_of_well_ + p] = rates[p];
+        }
+
+        // TODO: there will be something need to be done for the cases not the defaulted 3 phases,
+        // like 2 phases or solvent, polymer, etc. But we are not addressing them with THP control yet.
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    double
+    StandardWell<TypeTag>::
+    calculateBHPWithTHPTargetIPR() const
+    {
+        const double thp_target = this->getTHPConstraint();
+        const double thp_control_index = this->getTHPControlIndex();
+        const  int thp_table_id = well_controls_iget_vfp(well_controls_, thp_control_index);
+        const double alq = well_controls_iget_alq(well_controls_, thp_control_index);
+
+        // not considering injectors for now
+        const double vfp_ref_depth = vfp_properties_->getProd()->getTable(thp_table_id)->getDatumDepth();
+
+        // the density of the top perforation
+        // TODO: make sure this is properly initialized
+        // TODO: with IPR, it should be possible, even this well is a new well, we do not have
+        // TODO: any information of previous rates. However, we are lacking the pressure though.
+        // TODO: but let us not do more work related now to see what will happen
+        const double rho = perf_densities_[0];
+
+        // TODO: call a/the function for dp
+        const double dp = (vfp_ref_depth - ref_depth_) * rho * gravity_;
+
+        const double  bhp_limit = mostStrictBhpFromBhpLimits();
+
+        const double obtain_bhp = vfp_properties_->getProd()->calculateBhpWithTHPTarget(ipr_a_, ipr_b_,
+                                             bhp_limit, thp_table_id, thp_target, alq, dp);
+
+        // we should have made sure that this well should be operable under THP limit now
+        assert(obtain_bhp > 0.);
+
+        return obtain_bhp;
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::

opm/autodiff/WellInterface.hpp

@@ -172,10 +172,12 @@ namespace Opm
                                            const WellState& well_state,
                                            std::vector<double>& well_potentials) = 0;
 
-        virtual void updateWellStateWithTarget(WellState& well_state) const = 0;
+        virtual void updateWellStateWithTarget(/* const */ Simulator& ebos_simulator,
+                                               WellState& well_state) /* const */ = 0;
 
-        void updateWellControl(WellState& well_state,
+        void updateWellControl(/* const */ Simulator& ebos_simulator,
-                               wellhelpers::WellSwitchingLogger& logger) const;
+                               WellState& well_state,
+                               wellhelpers::WellSwitchingLogger& logger) /* const */;
 
         virtual void updatePrimaryVariables(const WellState& well_state) const = 0;
 

opm/autodiff/WellInterface_impl.hpp

@@ -415,8 +415,9 @@ namespace Opm
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::
-    updateWellControl(WellState& well_state,
+    updateWellControl(/* const */ Simulator& ebos_simulator,
-                      wellhelpers::WellSwitchingLogger& logger) const
+                      WellState& well_state,
+                      wellhelpers::WellSwitchingLogger& logger) /* const */
     {
         const int np = number_of_phases_;
         const int w = index_of_well_;
@@ -467,7 +468,7 @@ namespace Opm
         }
 
         if (updated_control_index != old_control_index) { //  || well_collection_->groupControlActive()) {
-            updateWellStateWithTarget(well_state);
+            updateWellStateWithTarget(ebos_simulator, well_state);
             updatePrimaryVariables(well_state);
         }
     }
@@ -1133,7 +1134,7 @@ namespace Opm
             solveEqAndUpdateWellState(well_state);
 
             wellhelpers::WellSwitchingLogger logger;
-            updateWellControl(well_state, logger);
+            updateWellControl(ebosSimulator, well_state, logger);
             initPrimaryVariablesEvaluation();
         } while (it < max_iter);