Merge pull request #1622 from GitPaean/updating_thp_when_VFPActive

Updating thp when valid VFP table is specified

opm/autodiff/AutoDiffHelpers.hpp

@@ -365,7 +365,7 @@ spdiag(const AutoDiffBlock<double>::V& d)
     public:
         typedef AutoDiffBlock<Scalar> ADB;
 
-        enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero, NotNaN };
+        enum CriterionForLeftElement { GreaterEqualZero, GreaterZero, Zero, NotEqualZero, LessZero, LessEqualZero, NotNaN, NotNaNInf };
 
         Selector(const typename ADB::V& selection_basis,
                  CriterionForLeftElement crit = GreaterEqualZero)
@@ -400,6 +400,9 @@ spdiag(const AutoDiffBlock<double>::V& d)
                 case NotNaN:
                     chooseleft = !isnan(selection_basis[i]);
                     break;
+                case NotNaNInf:
+                    chooseleft = !isnan(selection_basis[i]) && !std::isinf(selection_basis[i]);
+                    break;
                 default:
                     OPM_THROW(std::logic_error, "No such criterion: " << crit);
                 }

opm/autodiff/BlackoilWellModel_impl.hpp

@@ -844,8 +844,6 @@ namespace Opm {
         // we simply return.
         if( !wellsActive() ) return ;
 
-#if HAVE_OPENMP
-#endif // HAVE_OPENMP
         wellhelpers::WellSwitchingLogger logger;
 
         for (const auto& well : well_container_) {

opm/autodiff/StandardWell.hpp

@@ -331,7 +331,7 @@ namespace Opm
         template <class ValueType>
         ValueType calculateBhpFromThp(const std::vector<ValueType>& rates, const int control_index) const;
 
-        double calculateThpFromBhp(const std::vector<double>& rates, const int control_index, const double bhp) const;
+        double calculateThpFromBhp(const std::vector<double>& rates, const double bhp) const;
 
         // get the mobility for specific perforation
         void getMobility(const Simulator& ebosSimulator,

opm/autodiff/StandardWell_impl.hpp

@@ -1084,43 +1084,35 @@ namespace Opm
     StandardWell<TypeTag>::
     updateThp(WellState& well_state) const
     {
-        // for the wells having a THP constaint, we should update their thp value
+        // When there is no vaild VFP table provided, we set the thp to be zero.
-        // If it is under THP control, it will be set to be the target value.
+        if (!this->isVFPActive()) {
-        // TODO: a better standard is probably whether we have the table to calculate the THP value
+            well_state.thp()[index_of_well_] = 0.;
-        // TODO: it is something we need to check the output to decide.
+            return;
-        const WellControls* wc = well_controls_;
+        }
-        // TODO: we should only maintain one current control either from the well_state or from well_controls struct.
-        // Either one can be more favored depending on the final strategy for the initilzation of the well control
-        const int nwc = well_controls_get_num(wc);
-        // Looping over all controls until we find a THP constraint
-        for (int ctrl_index = 0; ctrl_index < nwc; ++ctrl_index) {
-            if (well_controls_iget_type(wc, ctrl_index) == THP) {
-                // the current control
-                const int current = well_state.currentControls()[index_of_well_];
-                // if well under THP control at the moment
-                if (current == ctrl_index) {
-                    const double thp_target = well_controls_iget_target(wc, current);
-                    well_state.thp()[index_of_well_] = thp_target;
-                } else { // otherwise we calculate the thp from the bhp value
-                    const Opm::PhaseUsage& pu = phaseUsage();
-                    std::vector<double> rates(3, 0.0);
 
-                    if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
+        // avaiable VFP table is provided, we should update the thp value
-                        rates[ Water ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Water ] ];
-                    }
-                    if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
-                        rates[ Oil ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Oil ] ];
-                    }
-                    if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
-                        rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ];
-                    }
 
-                    const double bhp = well_state.bhp()[index_of_well_];
+        // if the well is under THP control, we should use its target value
+        if (well_controls_get_current_type(well_controls_) == THP) {
+            well_state.thp()[index_of_well_] = well_controls_get_current_target(well_controls_);
+        } else {
+            // the well is under other control types, we calculate the thp based on bhp and rates
+            std::vector<double> rates(3, 0.0);
 
-                    well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, ctrl_index, bhp);
+            const Opm::PhaseUsage& pu = phaseUsage();
-                }
+            if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) {
-                break;
+                rates[ Water ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Water ] ];
             }
+            if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
+                rates[ Oil ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Oil ] ];
+            }
+            if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {
+                rates[ Gas ] = well_state.wellRates()[index_of_well_ * number_of_phases_ + pu.phase_pos[ Gas ] ];
+            }
+
+            const double bhp = well_state.bhp()[index_of_well_];
+
+            well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, bhp);
         }
     }
 
@@ -2096,7 +2088,6 @@ namespace Opm
     double
     StandardWell<TypeTag>::
     calculateThpFromBhp(const std::vector<double>& rates,
-                        const int control_index,
                         const double bhp) const
     {
         assert(int(rates.size()) == 3); // the vfp related only supports three phases now.
@@ -2105,32 +2096,30 @@ namespace Opm
         const double liquid = rates[Oil];
         const double vapour = rates[Gas];
 
-        const int vfp        = well_controls_iget_vfp(well_controls_, control_index);
-        const double& alq    = well_controls_iget_alq(well_controls_, control_index);
-
         // pick the density in the top layer
         const double rho = perf_densities_[0];
 
         double thp = 0.0;
         if (well_type_ == INJECTOR) {
-            const double vfp_ref_depth = vfp_properties_->getInj()->getTable(vfp)->getDatumDepth();
+            const int table_id = well_ecl_->getInjectionProperties(current_step_).VFPTableNumber;
-
+            const double vfp_ref_depth = vfp_properties_->getInj()->getTable(table_id)->getDatumDepth();
             const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
 
-            thp = vfp_properties_->getInj()->thp(vfp, aqua, liquid, vapour, bhp + dp);
+            thp = vfp_properties_->getInj()->thp(table_id, aqua, liquid, vapour, bhp + dp);
-         }
+        }
-         else if (well_type_ == PRODUCER) {
+        else if (well_type_ == PRODUCER) {
-             const double vfp_ref_depth = vfp_properties_->getProd()->getTable(vfp)->getDatumDepth();
+            const int table_id = well_ecl_->getProductionProperties(current_step_).VFPTableNumber;
+            const double alq = well_ecl_->getProductionProperties(current_step_).ALQValue;
+            const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
+            const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
 
-             const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
+            thp = vfp_properties_->getProd()->thp(table_id, aqua, liquid, vapour, bhp + dp, alq);
+        }
+        else {
+            OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
+        }
 
-             thp = vfp_properties_->getProd()->thp(vfp, aqua, liquid, vapour, bhp + dp, alq);
+        return thp;
-         }
-         else {
-             OPM_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well");
-         }
-
-         return thp;
     }
 
 

opm/autodiff/VFPHelpers.hpp

@@ -21,6 +21,7 @@
 #ifndef OPM_AUTODIFF_VFPHELPERS_HPP_
 #define OPM_AUTODIFF_VFPHELPERS_HPP_
 
+#include <opm/common/OpmLog/OpmLog.hpp>
 
 #include <cmath>
 #include <opm/common/ErrorMacros.hpp>
@@ -37,19 +38,31 @@ namespace detail {
 
 
 /**
- * Returns zero if input value is NaN
+ * Returns zero if input value is NaN of INF
  */
-inline double zeroIfNan(const double& value) {
+inline double zeroIfNanInf(const double& value) {
-    return (std::isnan(value)) ? 0.0 : value;
+    const bool nan_or_inf = std::isnan(value) || std::isinf(value);
+
+    if (nan_or_inf) {
+        OpmLog::warning("NAN_OR_INF_VFP", "NAN or INF value encountered during VFP calculation, the value is set to zero");
+    }
+
+    return nan_or_inf ? 0.0 : value;
 }
 
 
 /**
- * Returns zero if input value is NaN
+ * Returns zero if input value is NaN or INF
  */
 template <class EvalWell>
-inline EvalWell zeroIfNan(const EvalWell& value) {
+inline EvalWell zeroIfNanInf(const EvalWell& value) {
-    return (std::isnan(value.value())) ? 0.0 : value;
+    const bool nan_or_inf = std::isnan(value.value()) || std::isinf(value.value());
+
+    if (nan_or_inf) {
+        OpmLog::warning("NAN_OR_INF_VFP_EVAL", "NAN or INF Evalution encountered during VFP calculation, the Evalution is set to zero");
+    }
+
+    return nan_or_inf ? 0.0 : value;
 }
 
 
@@ -120,14 +133,14 @@ static T getWFR(const T& aqua, const T& liquid, const T& vapour,
         case VFPProdTable::WFR_WOR: {
             //Water-oil ratio = water / oil
             T wor = aqua / liquid;
-            return zeroIfNan(wor);
+            return zeroIfNanInf(wor);
         }
         case VFPProdTable::WFR_WCT:
             //Water cut = water / (water + oil)
-            return zeroIfNan(aqua / (aqua + liquid));
+            return zeroIfNanInf(aqua / (aqua + liquid));
         case VFPProdTable::WFR_WGR:
             //Water-gas ratio = water / gas
-            return zeroIfNan(aqua / vapour);
+            return zeroIfNanInf(aqua / vapour);
         case VFPProdTable::WFR_INVALID: //Intentional fall-through
         default:
             OPM_THROW(std::logic_error, "Invalid WFR_TYPE: '" << type << "'");
@@ -149,13 +162,13 @@ static T getGFR(const T& aqua, const T& liquid, const T& vapour,
     switch(type) {
         case VFPProdTable::GFR_GOR:
             // Gas-oil ratio = gas / oil
-            return zeroIfNan(vapour / liquid);
+            return zeroIfNanInf(vapour / liquid);
         case VFPProdTable::GFR_GLR:
             // Gas-liquid ratio = gas / (oil + water)
-            return zeroIfNan(vapour / (liquid + aqua));
+            return zeroIfNanInf(vapour / (liquid + aqua));
         case VFPProdTable::GFR_OGR:
             // Oil-gas ratio = oil / gas
-            return zeroIfNan(liquid / vapour);
+            return zeroIfNanInf(liquid / vapour);
         case VFPProdTable::GFR_INVALID: //Intentional fall-through
         default:
             OPM_THROW(std::logic_error, "Invalid GFR_TYPE: '" << type << "'");
@@ -537,13 +550,22 @@ template <typename T>
 const T* getTable(const std::map<int, T*> tables, int table_id) {
     auto entry = tables.find(table_id);
     if (entry == tables.end()) {
-        OPM_THROW(std::invalid_argument, "Nonexistent table " << table_id << " referenced.");
+        OPM_THROW(std::invalid_argument, "Nonexistent VFP table " << table_id << " referenced.");
     }
     else {
         return entry->second;
     }
 }
 
+/**
+ * Check whether we have a table with the table number
+ */
+template <typename T>
+bool hasTable(const std::map<int, T*> tables, int table_id) {
+    const auto entry = tables.find(table_id);
+    return (entry != tables.end() );
+}
+
 
 /**
  * Returns the type variable for FLO/GFR/WFR for production tables

opm/autodiff/VFPHelpersLegacy.hpp

@@ -37,15 +37,15 @@ typedef AutoDiffBlock<double> ADB;
 
 
 /**
- * Returns zero for every entry in the ADB which is NaN
+ *  * Returns zero for every entry in the ADB which is NaN or INF
- */
+ *   */
-inline ADB zeroIfNan(const ADB& values) {
+inline ADB zeroIfNanInf(const ADB& values) {
-    Selector<ADB::V::Scalar> not_nan_selector(values.value(), Selector<ADB::V::Scalar>::NotNaN);
+    Selector<ADB::V::Scalar> not_nan_inf_selector(values.value(), Selector<ADB::V::Scalar>::NotNaNInf);
 
     const ADB::V z = ADB::V::Zero(values.size());
     const ADB zero = ADB::constant(z, values.blockPattern());
 
-    ADB retval = not_nan_selector.select(values, zero);
+    ADB retval = not_nan_inf_selector.select(values, zero);
     return retval;
 }
 

opm/autodiff/VFPInjProperties.cpp

@@ -91,9 +91,12 @@ double VFPInjProperties::thp(int table_id,
     return retval;
 }
 
-
 const VFPInjTable* VFPInjProperties::getTable(const int table_id) const {
     return detail::getTable(m_tables, table_id);
 }
 
+bool VFPInjProperties::hasTable(const int table_id) const {
+    return detail::hasTable(m_tables, table_id);
+}
+
 } //Namespace Opm

opm/autodiff/VFPInjProperties.hpp

@@ -109,6 +109,11 @@ public:
      */
     const VFPInjTable* getTable(const int table_id) const;
 
+    /**
+     * Check whether there is table associated with ID
+     */
+    bool hasTable(const int table_id) const;
+
     /**
      * Returns true if no vfp tables are in the current map
      */

opm/autodiff/VFPProdProperties.cpp

@@ -105,5 +105,9 @@ const VFPProdTable* VFPProdProperties::getTable(const int table_id) const {
     return detail::getTable(m_tables, table_id);
 }
 
+bool VFPProdProperties::hasTable(const int table_id) const {
+    return detail::hasTable(m_tables, table_id);
+}
+
 
 }

opm/autodiff/VFPProdProperties.hpp

@@ -158,6 +158,11 @@ public:
      */
     const VFPProdTable* getTable(const int table_id) const;
 
+    /**
+     * Check whether there is table associated with ID
+     */
+    bool hasTable(const int table_id) const;
+
     /**
      * Returns true if no vfp tables are in the current map
      */

opm/autodiff/WellInterface.hpp

@@ -334,6 +334,9 @@ namespace Opm
 
         double scalingFactor(const int comp_idx) const;
 
+        // whether a well is specified with a non-zero and valid VFP table number
+        bool isVFPActive() 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);

opm/autodiff/WellInterface_impl.hpp

@@ -979,6 +979,48 @@ namespace Opm
 
 
 
+
+    template<typename TypeTag>
+    bool
+    WellInterface<TypeTag>::isVFPActive() const
+    {
+        // since the well_controls only handles the VFP number when THP constraint/target is there.
+        // we need to get the table number through the parser, in case THP constraint/target is not there.
+        // When THP control/limit is not active, if available VFP table is provided, we will still need to
+        // update THP value. However, it will only used for output purpose.
+
+        if (well_type_ == PRODUCER) { // producer
+            const int table_id = well_ecl_->getProductionProperties(current_step_).VFPTableNumber;
+            if (table_id <= 0) {
+                return false;
+            } else {
+                if (vfp_properties_->getProd()->hasTable(table_id)) {
+                    return true;
+                } else {
+                    OPM_THROW(std::runtime_error, "VFPPROD table " << std::to_string(table_id) << " is specfied,"
+                              << " for well " << name() << ", while we could not access it during simulation");
+                }
+            }
+
+        } else { // injector
+            const int table_id = well_ecl_->getInjectionProperties(current_step_).VFPTableNumber;
+            if (table_id <= 0) {
+                return false;
+            } else {
+                if (vfp_properties_->getInj()->hasTable(table_id)) {
+                    return true;
+                } else {
+                    OPM_THROW(std::runtime_error, "VFPINJ table " << std::to_string(table_id) << " is specfied,"
+                              << " for well " << name() << ", while we could not access it during simulation");
+                }
+            }
+        }
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::calculateReservoirRates(WellState& well_state) const