Merge pull request #3236 from akva2/wellinterface_generic_split

changed: split out non-typetag dependent code from WellInterface

CMakeLists_files.cmake

@@ -59,6 +59,7 @@ list (APPEND MAIN_SOURCE_FILES
   opm/simulators/wells/VFPProdProperties.cpp
   opm/simulators/wells/VFPInjProperties.cpp
   opm/simulators/wells/WellGroupHelpers.cpp
+  opm/simulators/wells/WellInterfaceGeneric.cpp
   opm/simulators/wells/WellProdIndexCalculator.cpp
   opm/simulators/wells/WellState.cpp
   opm/simulators/wells/WellStateFullyImplicitBlackoil.cpp

opm/simulators/wells/WellInterface.hpp

@@ -60,6 +60,8 @@ namespace Opm {
 #include <opm/material/densead/Math.hpp>
 #include <opm/material/densead/Evaluation.hpp>
 
+#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
+
 #include <string>
 #include <memory>
 #include <vector>
@@ -70,7 +72,7 @@ namespace Opm
 
 
     template<typename TypeTag>
-    class WellInterface
+    class WellInterface : public WellInterfaceGeneric
     {
     public:
 
@@ -146,27 +148,8 @@ namespace Opm
                       const int first_perf_index,
                       const std::vector<PerforationData>& perf_data);
 
-        /// Virutal destructor
+        /// Virtual destructor
-        virtual ~WellInterface() {}
+        virtual ~WellInterface() = default;
-
-        /// Well name.
-        const std::string& name() const;
-
-        /// True if the well is an injector.
-        bool isInjector() const;
-
-        /// True if the well is a producer.
-        bool isProducer() const;
-
-        /// Index of well in the wells struct and wellState
-        int indexOfWell() const;
-
-        /// Well cells.
-        const std::vector<int>& cells() const {return well_cells_; }
-
-        void setVFPProperties(const VFPProperties* vfp_properties_arg);
-
-        void setGuideRate(const GuideRate* guide_rate_arg);
 
         virtual void init(const PhaseUsage* phase_usage_arg,
                           const std::vector<double>& depth_arg,
@@ -201,10 +184,6 @@ namespace Opm
                                  WellTestState& wellTestState,
                                  DeferredLogger& deferred_logger) const;
 
-        void setWellEfficiencyFactor(const double efficiency_factor);
-
-        void setRepRadiusPerfLength(const std::vector<int>& cartesian_to_compressed);
-
         /// using the solution x to recover the solution xw for wells and applying
         /// xw to update Well State
         virtual void recoverWellSolutionAndUpdateWellState(const BVector& x,
@@ -273,10 +252,6 @@ namespace Opm
             return out;
         }
 
-        void closeCompletions(WellTestState& wellTestState);
-
-        const Well& wellEcl() const;
-
         // TODO: theoretically, it should be a const function
         // Simulator is not const is because that assembleWellEq is non-const Simulator
         void wellTesting(const Simulator& simulator,
@@ -285,8 +260,6 @@ namespace Opm
                          /* const */ WellState& well_state, const GroupState& group_state, WellTestState& welltest_state,
                          DeferredLogger& deferred_logger);
 
-        void updatePerforatedCell(std::vector<bool>& is_cell_perforated);
-
         void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger);
 
         // check whether the well is operable under the current reservoir condition
@@ -296,15 +269,6 @@ namespace Opm
                                    DeferredLogger& deferred_logger);
 
 
-        // whether the well is operable
-        bool isOperable() const;
-
-        /// Returns true if the well has one or more THP limits/constraints.
-        bool wellHasTHPConstraints(const SummaryState& summaryState) const;
-
-        /// Returns true if the well is currently in prediction mode (i.e. not history mode).
-        bool underPredictionMode() const;
-
         // update perforation water throughput based on solved water rate
         virtual void updateWaterThroughput(const double dt, WellState& well_state) const = 0;
 
@@ -320,29 +284,11 @@ namespace Opm
                                   WellState& well_state,
                                   DeferredLogger& deferred_logger) const;
 
-        void stopWell() {
-            this->wellStatus_ = Well::Status::STOP;
-        }
-
-        void openWell() {
-            this->wellStatus_ = Well::Status::OPEN;
-        }
-
-        bool wellIsStopped() const {
-            return this->wellStatus_ == Well::Status::STOP;
-        }
-
-        void setWsolvent(const double wsolvent);
-
-        void setDynamicThpLimit(const double thp_limit);
-
         void solveWellEquation(const Simulator& ebosSimulator,
                                WellState& well_state,
                                const GroupState& group_state,
                                DeferredLogger& deferred_logger);
 
-        const PhaseUsage& phaseUsage() const;
-
         virtual bool useInnerIterations() const = 0;
 
     protected:
@@ -350,109 +296,16 @@ namespace Opm
         // to indicate a invalid completion
         static const int INVALIDCOMPLETION = INT_MAX;
 
-        Well well_ecl_;
-
-        const ParallelWellInfo& parallel_well_info_;
-
-        const int current_step_;
-
         // simulation parameters
         const ModelParameters& param_;
 
-        // number of the perforations for this well
-        int number_of_perforations_;
-
-        // well index for each perforation
-        std::vector<double> well_index_;
-
-        // depth for each perforation
-        std::vector<double> perf_depth_;
-
-        // reference depth for the BHP
-        double ref_depth_;
-
-        double well_efficiency_factor_;
-
-        // cell index for each well perforation
-        std::vector<int> well_cells_;
-
-        // saturation table nubmer for each well perforation
-        std::vector<int> saturation_table_number_;
-
-        // representative radius of the perforations, used in shear calculation
-        std::vector<double> perf_rep_radius_;
-
-        // length of the perforations, use in shear calculation
-        std::vector<double> perf_length_;
-
-        // well bore diameter
-        std::vector<double> bore_diameters_;
-
-        /*
-         *  completions_ contains the mapping from completion id to connection indices
-         *  {
-         *      2 : [ConnectionIndex, ConnectionIndex],
-         *      1 : [ConnectionIndex, ConnectionIndex, ConnectionIndex],
-         *      5 : [ConnectionIndex],
-         *      7 : [ConnectionIndex]
-         *      ...
-         *   }
-         *   The integer IDs correspond to the COMPLETION id given by the COMPLUMP keyword.
-         *   When there is no COMPLUMP keyword used, a default completion number will be assigned
-         *   based on the order of the declaration of the connections.
-         *   Since the connections not OPEN is not included in the Wells, so they will not be considered
-         *   in this mapping relation.
-         */
-        std::map<int, std::vector<int>> completions_;
-
-        const PhaseUsage* phase_usage_;
-
-        bool getAllowCrossFlow() const;
-
-        const VFPProperties* vfp_properties_;
-
-        const GuideRate* guide_rate_;
-
-        double gravity_;
-
         // For the conversion between the surface volume rate and resrevoir voidage rate
         const RateConverterType& rateConverter_;
 
-        // The pvt region of the well. We assume
-        // We assume a well to not penetrate more than one pvt region.
-        const int pvtRegionIdx_;
-
-        const int num_components_;
-
-        // number of phases
-        int number_of_phases_;
-
-        // the index of well in Wells struct
-        int index_of_well_;
-
-        // record the index of the first perforation
-        // of states of individual well.
-        int first_perf_;
-
-        const std::vector<PerforationData>* perf_data_;
-
         std::vector<RateVector> connectionRates_;
 
-        Well::Status wellStatus_;
-
-        double wsolvent_;
-
-        std::optional<double> dynamic_thp_limit_;
-
         std::vector< Scalar > B_avg_;
 
-        // the vectors used to describe the inflow performance relationship (IPR)
-        // Q = IPR_A - BHP * IPR_B
-        // TODO: it minght need to go to WellInterface, let us implement it in StandardWell first
-        // it is only updated and used for producers for now
-        mutable std::vector<double> ipr_a_;
-        mutable std::vector<double> ipr_b_;
-
         bool changed_to_stopped_this_step_ = false;
 
         int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
@@ -461,8 +314,6 @@ namespace Opm
 
         int ebosCompIdxToFlowCompIdx( const unsigned compIdx ) const;
 
-        double wsolvent() const;
-
         double wpolymer() const;
 
         double wfoam() const;
@@ -473,20 +324,14 @@ namespace Opm
                                  const std::vector<double>& well_rates,
                                  DeferredLogger& deferred_logger) const;
 
-        double getTHPConstraint(const SummaryState& summaryState) const;
-
         template <class ValueType>
         ValueType calculateBhpFromThp(const WellState& well_state, const std::vector<ValueType>& rates, const Well& well, const SummaryState& summaryState, DeferredLogger& deferred_logger) const;
 
-        double getALQ(const WellState& well_state) const;
-
         virtual double getRefDensity() const = 0;
 
         // Component fractions for each phase for the well
         const std::vector<double>& compFrac() const;
 
-        double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
-
         struct RatioLimitCheckReport;
 
         void checkMaxWaterCutLimit(const WellEconProductionLimits& econ_production_limits,
@@ -522,13 +367,6 @@ namespace Opm
 
         std::vector<double> initialWellRateFractions(const Simulator& ebosSimulator, const std::vector<double>& potentials) const;
 
-        // whether a well is specified with a non-zero and valid VFP table number
-        bool isVFPActive(DeferredLogger& deferred_logger) const;
-
-        struct OperabilityStatus;
-
-        OperabilityStatus operability_status_;
-
         // check whether the well is operable under BHP limit with current reservoir condition
         virtual void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, DeferredLogger& deferred_logger) =0;
 
@@ -578,17 +416,9 @@ namespace Opm
                                          WellTestState& well_test_state,
                                          DeferredLogger& deferred_logger) const;
 
-        void updateWellTestStatePhysical(const WellState& well_state,
-                                         const double simulation_time,
-                                         const bool write_message_to_opmlog,
-                                         WellTestState& well_test_state,
-                                         DeferredLogger& deferred_logger) const;
-
         void solveWellForTesting(const Simulator& ebosSimulator, WellState& well_state, const GroupState& group_state,
                                  DeferredLogger& deferred_logger);
 
-        void initCompletions();
-
         bool checkConstraints(WellState& well_state,
                               const GroupState& group_state,
                               const Schedule& schedule,
@@ -669,53 +499,6 @@ namespace Opm
                                    DeferredLogger& deferred_logger);
     };
 
-
-
-
-
-    // definition of the struct OperabilityStatus
-    template<typename TypeTag>
-    struct
-    WellInterface<TypeTag>::
-    OperabilityStatus {
-        bool isOperable() const {
-            if (!operable_under_only_bhp_limit) {
-                return false;
-            } else {
-                return ( (isOperableUnderBHPLimit() || isOperableUnderTHPLimit()) );
-            }
-        }
-
-        bool isOperableUnderBHPLimit() const {
-            return operable_under_only_bhp_limit && obey_thp_limit_under_bhp_limit;
-        }
-
-        bool isOperableUnderTHPLimit() const {
-            return can_obtain_bhp_with_thp_limit && obey_bhp_limit_with_thp_limit;
-        }
-
-        void reset() {
-            operable_under_only_bhp_limit = true;
-            obey_thp_limit_under_bhp_limit = true;
-            can_obtain_bhp_with_thp_limit = true;
-            obey_bhp_limit_with_thp_limit = true;
-        }
-
-        // whether the well can be operated under bhp limit
-        // without considering other limits.
-        // if it is false, then the well is not operable for sure.
-        bool operable_under_only_bhp_limit = true;
-        // if the well can be operated under bhp limit, will it obey(not violate)
-        // the thp limit when operated under bhp limit
-        bool obey_thp_limit_under_bhp_limit = true;
-        // whether the well operate under the thp limit only
-        bool can_obtain_bhp_with_thp_limit = true;
-        // whether the well obey bhp limit when operated under thp limit
-        bool obey_bhp_limit_with_thp_limit = true;
-
-    };
-
-
     template<typename TypeTag>
     struct
     WellInterface<TypeTag>::

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -0,0 +1,386 @@
+/*
+  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2017 Statoil ASA.
+  Copyright 2018 IRIS
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <config.h>
+#include <opm/simulators/wells/WellInterfaceGeneric.hpp>
+
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/WellTestState.hpp>
+
+#include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
+#include <opm/simulators/wells/PerforationData.hpp>
+#include <opm/simulators/wells/ParallelWellInfo.hpp>
+#include <opm/simulators/wells/VFPProperties.hpp>
+#include <opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp>
+
+#include <cassert>
+#include <cmath>
+#include <cstddef>
+#include <stdexcept>
+
+namespace Opm
+{
+
+WellInterfaceGeneric::WellInterfaceGeneric(const Well& well,
+                                           const ParallelWellInfo& pw_info,
+                                           const int time_step,
+                                           const int pvtRegionIdx,
+                                           const int num_components,
+                                           const int num_phases,
+                                           const int index_of_well,
+                                           const int first_perf_index,
+                                           const std::vector<PerforationData>& perf_data)
+      : well_ecl_(well)
+      , parallel_well_info_(pw_info)
+      , current_step_(time_step)
+      , pvtRegionIdx_(pvtRegionIdx)
+      , num_components_(num_components)
+      , number_of_phases_(num_phases)
+      , index_of_well_(index_of_well)
+      , first_perf_(first_perf_index)
+      , perf_data_(&perf_data)
+      , ipr_a_(number_of_phases_)
+      , ipr_b_(number_of_phases_)
+{
+    assert(well.name()==pw_info.name());
+    assert(std::is_sorted(perf_data.begin(), perf_data.end(),
+                          [](const auto& perf1, const auto& perf2){
+        return perf1.ecl_index < perf2.ecl_index;
+    }));
+    if (time_step < 0) {
+        OPM_THROW(std::invalid_argument, "Negtive time step is used to construct WellInterface");
+    }
+
+    ref_depth_ = well.getRefDepth();
+
+    // We do not want to count SHUT perforations here, so
+    // it would be wrong to use wells.getConnections().size().
+    number_of_perforations_ = perf_data.size();
+
+    // perforations related
+    {
+        well_cells_.resize(number_of_perforations_);
+        well_index_.resize(number_of_perforations_);
+        saturation_table_number_.resize(number_of_perforations_);
+        int perf = 0;
+        for (const auto& pd : perf_data) {
+            well_cells_[perf] = pd.cell_index;
+            well_index_[perf] = pd.connection_transmissibility_factor;
+            saturation_table_number_[perf] = pd.satnum_id;
+            ++perf;
+        }
+    }
+
+    // initialization of the completions mapping
+    initCompletions();
+
+    well_efficiency_factor_ = 1.0;
+
+    this->wellStatus_ = Well::Status::OPEN;
+    if (well.getStatus() == Well::Status::STOP) {
+        this->wellStatus_ = Well::Status::STOP;
+    }
+
+    wsolvent_ = 0.0;
+}
+
+const std::string& WellInterfaceGeneric::name() const
+{
+    return well_ecl_.name();
+}
+
+bool WellInterfaceGeneric::isInjector() const
+{
+    return well_ecl_.isInjector();
+}
+
+bool WellInterfaceGeneric::isProducer() const
+{
+    return well_ecl_.isProducer();
+}
+
+int WellInterfaceGeneric::indexOfWell() const
+{
+    return index_of_well_;
+}
+
+bool WellInterfaceGeneric::getAllowCrossFlow() const
+{
+    return well_ecl_.getAllowCrossFlow();
+}
+
+const Well& WellInterfaceGeneric::wellEcl() const
+{
+    return well_ecl_;
+}
+
+const PhaseUsage& WellInterfaceGeneric::phaseUsage() const
+{
+    assert(phase_usage_ != nullptr);
+
+    return *phase_usage_;
+}
+
+double WellInterfaceGeneric::wsolvent() const
+{
+    return wsolvent_;
+}
+
+bool WellInterfaceGeneric::wellHasTHPConstraints(const SummaryState& summaryState) const
+{
+    if (dynamic_thp_limit_) {
+        return true;
+    }
+
+    if (well_ecl_.isInjector()) {
+        const auto controls = well_ecl_.injectionControls(summaryState);
+        if (controls.hasControl(Well::InjectorCMode::THP))
+            return true;
+    }
+
+    if (well_ecl_.isProducer( )) {
+        const auto controls = well_ecl_.productionControls(summaryState);
+        if (controls.hasControl(Well::ProducerCMode::THP))
+            return true;
+    }
+
+    return false;
+
+}
+
+double WellInterfaceGeneric::mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const
+{
+    if (well_ecl_.isInjector()) {
+        const auto& controls = well_ecl_.injectionControls(summaryState);
+        return controls.bhp_limit;
+    }
+
+    if (well_ecl_.isProducer( )) {
+        const auto& controls = well_ecl_.productionControls(summaryState);
+        return controls.bhp_limit;
+    }
+
+    return 0.0;
+}
+
+double WellInterfaceGeneric::getTHPConstraint(const SummaryState& summaryState) const
+{
+    if (dynamic_thp_limit_) {
+        return *dynamic_thp_limit_;
+    }
+    if (well_ecl_.isInjector()) {
+        const auto& controls = well_ecl_.injectionControls(summaryState);
+        return controls.thp_limit;
+    }
+
+    if (well_ecl_.isProducer( )) {
+        const auto& controls = well_ecl_.productionControls(summaryState);
+        return controls.thp_limit;
+    }
+
+    return 0.0;
+}
+
+bool WellInterfaceGeneric::underPredictionMode() const
+{
+    return well_ecl_.predictionMode();
+}
+
+void WellInterfaceGeneric::initCompletions()
+{
+    assert(completions_.empty() );
+
+    const WellConnections& connections = well_ecl_.getConnections();
+    const std::size_t num_conns = connections.size();
+
+    int num_active_connections = 0;
+    auto my_next_perf = perf_data_->begin();
+    for (std::size_t c = 0; c < num_conns; ++c) {
+        if (my_next_perf == perf_data_->end())
+        {
+            break;
+        }
+        if (my_next_perf->ecl_index > c)
+        {
+            continue;
+        }
+        assert(my_next_perf->ecl_index == c);
+        if (connections[c].state() == Connection::State::OPEN) {
+            completions_[connections[c].complnum()].push_back(num_active_connections++);
+        }
+        ++my_next_perf;
+    }
+    assert(my_next_perf == perf_data_->end());
+}
+
+void WellInterfaceGeneric::closeCompletions(WellTestState& wellTestState)
+{
+    const auto& connections = well_ecl_.getConnections();
+    int perfIdx = 0;
+    for (const auto& connection : connections) {
+        if (connection.state() == Connection::State::OPEN) {
+            if (wellTestState.hasCompletion(name(), connection.complnum())) {
+                well_index_[perfIdx] = 0.0;
+            }
+            perfIdx++;
+        }
+    }
+}
+
+void WellInterfaceGeneric::setVFPProperties(const VFPProperties* vfp_properties_arg)
+{
+    vfp_properties_ = vfp_properties_arg;
+}
+
+void WellInterfaceGeneric::setGuideRate(const GuideRate* guide_rate_arg)
+{
+    guide_rate_ = guide_rate_arg;
+}
+
+void WellInterfaceGeneric::setWellEfficiencyFactor(const double efficiency_factor)
+{
+    well_efficiency_factor_ = efficiency_factor;
+}
+
+void WellInterfaceGeneric::setRepRadiusPerfLength(const std::vector<int>& cartesian_to_compressed)
+{
+    const int nperf = number_of_perforations_;
+
+    perf_rep_radius_.clear();
+    perf_length_.clear();
+    bore_diameters_.clear();
+
+    perf_rep_radius_.reserve(nperf);
+    perf_length_.reserve(nperf);
+    bore_diameters_.reserve(nperf);
+
+    // COMPDAT handling
+    const auto& connectionSet = well_ecl_.getConnections();
+    CheckDistributedWellConnections checker(well_ecl_, parallel_well_info_);
+    for (size_t c=0; c<connectionSet.size(); c++) {
+        const auto& connection = connectionSet.get(c);
+        const int cell =
+                cartesian_to_compressed[connection.global_index()];
+        if (connection.state() != Connection::State::OPEN || cell >= 0)
+        {
+            checker.connectionFound(c);
+        }
+        if (connection.state() == Connection::State::OPEN) {
+
+            if (cell >= 0) {
+                double radius = connection.rw();
+                double re = connection.re(); // area equivalent radius of the grid block
+                double perf_length = connection.connectionLength(); // the length of the well perforation
+                const double repR = std::sqrt(re * radius);
+                perf_rep_radius_.push_back(repR);
+                perf_length_.push_back(perf_length);
+                bore_diameters_.push_back(2. * radius);
+            }
+        }
+    }
+    checker.checkAllConnectionsFound();
+}
+
+void WellInterfaceGeneric::setWsolvent(const double wsolvent)
+{
+    wsolvent_ = wsolvent;
+}
+
+void WellInterfaceGeneric::setDynamicThpLimit(const double thp_limit)
+{
+    dynamic_thp_limit_ = thp_limit;
+}
+
+void WellInterfaceGeneric::updatePerforatedCell(std::vector<bool>& is_cell_perforated)
+{
+
+    for (int perf_idx = 0; perf_idx<number_of_perforations_; ++perf_idx) {
+        is_cell_perforated[well_cells_[perf_idx]] = true;
+    }
+}
+
+bool WellInterfaceGeneric::isVFPActive(DeferredLogger& deferred_logger) 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 (isProducer()) { // producer
+        const int table_id = well_ecl_.vfp_table_number();
+        if (table_id <= 0) {
+            return false;
+        } else {
+            if (vfp_properties_->getProd()->hasTable(table_id)) {
+                return true;
+            } else {
+                OPM_DEFLOG_THROW(std::runtime_error, "VFPPROD table " << std::to_string(table_id) << " is specfied,"
+                              << " for well " << name() << ", while we could not access it during simulation", deferred_logger);
+            }
+        }
+
+    } else { // injector
+        const int table_id = well_ecl_.vfp_table_number();
+        if (table_id <= 0) {
+            return false;
+        } else {
+            if (vfp_properties_->getInj()->hasTable(table_id)) {
+                return true;
+            } else {
+                OPM_DEFLOG_THROW(std::runtime_error, "VFPINJ table " << std::to_string(table_id) << " is specfied,"
+                              << " for well " << name() << ", while we could not access it during simulation", deferred_logger);
+            }
+        }
+    }
+}
+
+void WellInterfaceGeneric::updateWellTestStatePhysical(const WellState& /* well_state */,
+                                                       const double simulation_time,
+                                                       const bool write_message_to_opmlog,
+                                                       WellTestState& well_test_state,
+                                                       DeferredLogger& deferred_logger) const
+{
+    if (!isOperable()) {
+        if (well_test_state.hasWellClosed(name(), WellTestConfig::Reason::ECONOMIC) ||
+            well_test_state.hasWellClosed(name(), WellTestConfig::Reason::PHYSICAL) ) {
+            // Already closed, do nothing.
+        } else {
+            well_test_state.closeWell(name(), WellTestConfig::Reason::PHYSICAL, simulation_time);
+            if (write_message_to_opmlog) {
+                const std::string action = well_ecl_.getAutomaticShutIn() ? "shut" : "stopped";
+                const std::string msg = "Well " + name()
+                    + " will be " + action + " as it can not operate under current reservoir conditions.";
+                deferred_logger.info(msg);
+            }
+        }
+    }
+}
+
+bool WellInterfaceGeneric::isOperable() const
+{
+    return operability_status_.isOperable();
+}
+
+double WellInterfaceGeneric::getALQ(const WellState& well_state) const
+{
+    return well_state.getALQ(name());
+}
+
+
+} // namespace Opm

opm/simulators/wells/WellInterfaceGeneric.hpp

@@ -0,0 +1,255 @@
+/*
+  Copyright 2017 SINTEF Digital, Mathematics and Cybernetics.
+  Copyright 2017 Statoil ASA.
+  Copyright 2017 IRIS
+  Copyright 2019 Norce
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#ifndef OPM_WELLINTERFACE_GENERIC_HEADER_INCLUDED
+#define OPM_WELLINTERFACE_GENERIC_HEADER_INCLUDED
+
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
+
+#include <map>
+#include <optional>
+#include <string>
+#include <vector>
+
+namespace Opm
+{
+
+class DeferredLogger;
+class GuideRate;
+class ParallelWellInfo;
+struct PerforationData;
+struct PhaseUsage;
+class SummaryState;
+class VFPProperties;
+class WellTestState;
+class WellStateFullyImplicitBlackoil;
+
+class WellInterfaceGeneric {
+public:
+    using WellState = WellStateFullyImplicitBlackoil;
+
+    WellInterfaceGeneric(const Well& well,
+                         const ParallelWellInfo& parallel_well_info,
+                         const int time_step,
+                         const int pvtRegionIdx,
+                         const int num_components,
+                         const int num_phases,
+                         const int index_of_well,
+                         const int first_perf_index,
+                         const std::vector<PerforationData>& perf_data);
+
+    /// Well name.
+    const std::string& name() const;
+
+    /// True if the well is an injector.
+    bool isInjector() const;
+
+    /// True if the well is a producer.
+    bool isProducer() const;
+
+    /// Well cells.
+    const std::vector<int>& cells() const { return well_cells_; }
+
+    /// Index of well in the wells struct and wellState
+    int indexOfWell() const;
+
+    const Well& wellEcl() const;
+    const PhaseUsage& phaseUsage() const;
+
+    /// Returns true if the well is currently in prediction mode (i.e. not history mode).
+    bool underPredictionMode() const;
+
+    // whether the well is operable
+    bool isOperable() const;
+
+    void initCompletions();
+    void closeCompletions(WellTestState& wellTestState);
+
+    void setVFPProperties(const VFPProperties* vfp_properties_arg);
+    void setGuideRate(const GuideRate* guide_rate_arg);
+    void setWellEfficiencyFactor(const double efficiency_factor);
+    void setRepRadiusPerfLength(const std::vector<int>& cartesian_to_compressed);
+    void setWsolvent(const double wsolvent);
+    void setDynamicThpLimit(const double thp_limit);
+    void updatePerforatedCell(std::vector<bool>& is_cell_perforated);
+
+    /// Returns true if the well has one or more THP limits/constraints.
+    bool wellHasTHPConstraints(const SummaryState& summaryState) const;
+
+    void stopWell() {
+        this->wellStatus_ = Well::Status::STOP;
+    }
+
+    void openWell() {
+        this->wellStatus_ = Well::Status::OPEN;
+    }
+
+    bool wellIsStopped() const {
+        return this->wellStatus_ == Well::Status::STOP;
+    }
+
+protected:
+    // whether a well is specified with a non-zero and valid VFP table number
+    bool isVFPActive(DeferredLogger& deferred_logger) const;
+
+    bool getAllowCrossFlow() const;
+    double wsolvent() const;
+    double mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const;
+    double getTHPConstraint(const SummaryState& summaryState) const;
+    void updateWellTestStatePhysical(const WellState& well_state,
+                                     const double simulation_time,
+                                     const bool write_message_to_opmlog,
+                                     WellTestState& well_test_state,
+                                     DeferredLogger& deferred_logger) const;
+    double getALQ(const WellState& well_state) const;
+
+    // definition of the struct OperabilityStatus
+    struct OperabilityStatus {
+        bool isOperable() const {
+            if (!operable_under_only_bhp_limit) {
+                return false;
+            } else {
+                return ( (isOperableUnderBHPLimit() || isOperableUnderTHPLimit()) );
+            }
+        }
+
+        bool isOperableUnderBHPLimit() const {
+            return operable_under_only_bhp_limit && obey_thp_limit_under_bhp_limit;
+        }
+
+        bool isOperableUnderTHPLimit() const {
+            return can_obtain_bhp_with_thp_limit && obey_bhp_limit_with_thp_limit;
+        }
+
+        void reset() {
+            operable_under_only_bhp_limit = true;
+            obey_thp_limit_under_bhp_limit = true;
+            can_obtain_bhp_with_thp_limit = true;
+            obey_bhp_limit_with_thp_limit = true;
+        }
+
+        // whether the well can be operated under bhp limit
+        // without considering other limits.
+        // if it is false, then the well is not operable for sure.
+        bool operable_under_only_bhp_limit = true;
+        // if the well can be operated under bhp limit, will it obey(not violate)
+        // the thp limit when operated under bhp limit
+        bool obey_thp_limit_under_bhp_limit = true;
+        // whether the well operate under the thp limit only
+        bool can_obtain_bhp_with_thp_limit = true;
+        // whether the well obey bhp limit when operated under thp limit
+        bool obey_bhp_limit_with_thp_limit = true;
+    };
+
+    OperabilityStatus operability_status_;
+
+    Well well_ecl_;
+
+    const ParallelWellInfo& parallel_well_info_;
+    const int current_step_;
+
+    // The pvt region of the well. We assume
+    // We assume a well to not penetrate more than one pvt region.
+    const int pvtRegionIdx_;
+
+    const int num_components_;
+
+    // number of phases
+    int number_of_phases_;
+
+    // the index of well in Wells struct
+    int index_of_well_;
+
+    // record the index of the first perforation
+    // of states of individual well.
+    int first_perf_;
+
+    const std::vector<PerforationData>* perf_data_;
+
+    // the vectors used to describe the inflow performance relationship (IPR)
+    // Q = IPR_A - BHP * IPR_B
+    // TODO: it minght need to go to WellInterface, let us implement it in StandardWell first
+    // it is only updated and used for producers for now
+    mutable std::vector<double> ipr_a_;
+    mutable std::vector<double> ipr_b_;
+
+    // cell index for each well perforation
+    std::vector<int> well_cells_;
+
+    // well index for each perforation
+    std::vector<double> well_index_;
+
+    // number of the perforations for this well
+    int number_of_perforations_;
+
+    // depth for each perforation
+    std::vector<double> perf_depth_;
+
+    // representative radius of the perforations, used in shear calculation
+    std::vector<double> perf_rep_radius_;
+
+    // length of the perforations, use in shear calculation
+    std::vector<double> perf_length_;
+
+    // well bore diameter
+    std::vector<double> bore_diameters_;
+
+    /*
+     *  completions_ contains the mapping from completion id to connection indices
+     *  {
+     *      2 : [ConnectionIndex, ConnectionIndex],
+     *      1 : [ConnectionIndex, ConnectionIndex, ConnectionIndex],
+     *      5 : [ConnectionIndex],
+     *      7 : [ConnectionIndex]
+     *      ...
+     *   }
+     *   The integer IDs correspond to the COMPLETION id given by the COMPLUMP keyword.
+     *   When there is no COMPLUMP keyword used, a default completion number will be assigned
+     *   based on the order of the declaration of the connections.
+     *   Since the connections not OPEN is not included in the Wells, so they will not be considered
+     *   in this mapping relation.
+     */
+    std::map<int, std::vector<int>> completions_;
+
+    // reference depth for the BHP
+    double ref_depth_;
+
+    // saturation table nubmer for each well perforation
+    std::vector<int> saturation_table_number_;
+
+    Well::Status wellStatus_;
+
+    const PhaseUsage* phase_usage_;
+
+    double gravity_;
+    double wsolvent_;
+    std::optional<double> dynamic_thp_limit_;
+
+    double well_efficiency_factor_;
+    const VFPProperties* vfp_properties_;
+    const GuideRate* guide_rate_;
+};
+
+}
+
+#endif // OPM_WELLINTERFACE_HEADER_INCLUDED

opm/simulators/wells/WellInterface_impl.hpp

@@ -40,63 +40,14 @@ namespace Opm
                   const int index_of_well,
                   const int first_perf_index,
                   const std::vector<PerforationData>& perf_data)
-      : well_ecl_(well)
+      : WellInterfaceGeneric(well, pw_info, time_step, pvtRegionIdx,
-      , parallel_well_info_(pw_info)
+                             num_components, num_phases, index_of_well,
-      , current_step_(time_step)
+                             first_perf_index, perf_data)
       , param_(param)
       , rateConverter_(rate_converter)
-      , pvtRegionIdx_(pvtRegionIdx)
-      , num_components_(num_components)
-      , number_of_phases_(num_phases)
-      , index_of_well_(index_of_well)
-      , first_perf_(first_perf_index)
-      , perf_data_(&perf_data)
-      , ipr_a_(number_of_phases_)
-      , ipr_b_(number_of_phases_)
     {
-        assert(well.name()==pw_info.name());
-        assert(std::is_sorted(perf_data.begin(), perf_data.end(),
-                              [](const auto& perf1, const auto& perf2){
-                                  return perf1.ecl_index < perf2.ecl_index;
-                              }));
-        if (time_step < 0) {
-            OPM_THROW(std::invalid_argument, "Negtive time step is used to construct WellInterface");
-        }
-
-        ref_depth_ = well.getRefDepth();
-
-        // We do not want to count SHUT perforations here, so
-        // it would be wrong to use wells.getConnections().size().
-        number_of_perforations_ = perf_data.size();
-
-        // perforations related
-        {
-            well_cells_.resize(number_of_perforations_);
-            well_index_.resize(number_of_perforations_);
-            saturation_table_number_.resize(number_of_perforations_);
-            int perf = 0;
-            for (const auto& pd : perf_data) {
-                well_cells_[perf] = pd.cell_index;
-                well_index_[perf] = pd.connection_transmissibility_factor;
-                saturation_table_number_[perf] = pd.satnum_id;
-                ++perf;
-            }
-        }
-
-        // initialization of the completions mapping
-        initCompletions();
-
-        well_efficiency_factor_ = 1.0;
-
         connectionRates_.resize(number_of_perforations_);
 
-        this->wellStatus_ = Well::Status::OPEN;
-        if (well.getStatus() == Well::Status::STOP) {
-            this->wellStatus_ = Well::Status::STOP;
-        }
-
-        wsolvent_ = 0.0;
-
         if constexpr (has_solvent || has_zFraction) {
             if (well.isInjector()) {
                 auto injectorType = well_ecl_.injectorType();
@@ -107,18 +58,6 @@ namespace Opm
         }
     }
 
-    template<typename TypeTag>
-    void
-    WellInterface<TypeTag>::
-    updatePerforatedCell(std::vector<bool>& is_cell_perforated)
-    {
-
-        for (int perf_idx = 0; perf_idx<number_of_perforations_; ++perf_idx) {
-            is_cell_perforated[well_cells_[perf_idx]] = true;
-        }
-    }
-
-
 
     template<typename TypeTag>
     void
@@ -135,153 +74,6 @@ namespace Opm
     }
 
 
-
-
-
-    template<typename TypeTag>
-    void
-    WellInterface<TypeTag>::
-    initCompletions()
-    {
-        assert(completions_.empty() );
-
-        const WellConnections& connections = well_ecl_.getConnections();
-        const std::size_t num_conns = connections.size();
-
-        int num_active_connections = 0;
-        auto my_next_perf = perf_data_->begin();
-        for (std::size_t c = 0; c < num_conns; ++c) {
-            if (my_next_perf == perf_data_->end())
-            {
-                break;
-            }
-            if (my_next_perf->ecl_index > c)
-            {
-                continue;
-            }
-            assert(my_next_perf->ecl_index == c);
-            if (connections[c].state() == Connection::State::OPEN) {
-                completions_[connections[c].complnum()].push_back(num_active_connections++);
-            }
-            ++my_next_perf;
-        }
-        assert(my_next_perf == perf_data_->end());
-    }
-
-
-
-
-
-
-    template<typename TypeTag>
-    void
-    WellInterface<TypeTag>::
-    setVFPProperties(const VFPProperties* vfp_properties_arg)
-    {
-        vfp_properties_ = vfp_properties_arg;
-    }
-
-    template<typename TypeTag>
-    void
-    WellInterface<TypeTag>::
-    setGuideRate(const GuideRate* guide_rate_arg)
-    {
-        guide_rate_ = guide_rate_arg;
-    }
-
-
-    template<typename TypeTag>
-    const std::string&
-    WellInterface<TypeTag>::
-    name() const
-    {
-        return well_ecl_.name();
-    }
-
-
-
-
-
-    template<typename TypeTag>
-    bool
-    WellInterface<TypeTag>::
-    isInjector() const
-    {
-        return well_ecl_.isInjector();
-    }
-
-
-
-
-    template<typename TypeTag>
-    bool
-    WellInterface<TypeTag>::
-    isProducer() const
-    {
-        return well_ecl_.isProducer();
-    }
-
-
-
-
-    template<typename TypeTag>
-    int
-    WellInterface<TypeTag>::
-    indexOfWell() const
-    {
-        return index_of_well_;
-    }
-
-
-
-
-
-
-    template<typename TypeTag>
-    bool
-    WellInterface<TypeTag>::
-    getAllowCrossFlow() const
-    {
-        return well_ecl_.getAllowCrossFlow();
-    }
-
-
-
-
-    template<typename TypeTag>
-    void
-    WellInterface<TypeTag>::
-    setWellEfficiencyFactor(const double efficiency_factor)
-    {
-        well_efficiency_factor_ = efficiency_factor;
-    }
-
-
-
-    template<typename TypeTag>
-    const Well&
-    WellInterface<TypeTag>::
-    wellEcl() const
-    {
-      return well_ecl_;
-    }
-
-
-
-    template<typename TypeTag>
-    const PhaseUsage&
-    WellInterface<TypeTag>::
-    phaseUsage() const
-    {
-        assert(phase_usage_ != nullptr);
-
-        return *phase_usage_;
-    }
-
-
-
-
-
     template<typename TypeTag>
     int
     WellInterface<TypeTag>::
@@ -336,40 +128,6 @@ namespace Opm
 
 
 
-    template<typename TypeTag>
-    double
-    WellInterface<TypeTag>::
-    wsolvent() const
-    {
-        return wsolvent_;
-    }
-
-
-
-    template<typename TypeTag>
-    void
-    WellInterface<TypeTag>::
-    setWsolvent(const double wsolvent)
-    {
-       wsolvent_ = wsolvent;
-    }
-
-
-
-
-
-    template<typename TypeTag>
-    void
-    WellInterface<TypeTag>::
-    setDynamicThpLimit(const double thp_limit)
-    {
-       dynamic_thp_limit_ = thp_limit;
-    }
-
-
-
-
-
     template<typename TypeTag>
     double
     WellInterface<TypeTag>::
@@ -438,78 +196,6 @@ namespace Opm
     }
 
 
-
-
-    template<typename TypeTag>
-    bool
-    WellInterface<TypeTag>::
-    wellHasTHPConstraints(const SummaryState& summaryState) const
-    {
-        if (dynamic_thp_limit_) {
-            return true;
-        }
-
-        if (well_ecl_.isInjector()) {
-            const auto controls = well_ecl_.injectionControls(summaryState);
-            if (controls.hasControl(Well::InjectorCMode::THP))
-                return true;
-        }
-
-        if (well_ecl_.isProducer( )) {
-            const auto controls = well_ecl_.productionControls(summaryState);
-            if (controls.hasControl(Well::ProducerCMode::THP))
-                return true;
-        }
-
-        return false;
-
-    }
-
-    template<typename TypeTag>
-    double
-    WellInterface<TypeTag>::
-    mostStrictBhpFromBhpLimits(const SummaryState& summaryState) const
-    {
-        if (well_ecl_.isInjector()) {
-            const auto& controls = well_ecl_.injectionControls(summaryState);
-            return controls.bhp_limit;
-        }
-
-        if (well_ecl_.isProducer( )) {
-            const auto& controls = well_ecl_.productionControls(summaryState);
-            return controls.bhp_limit;
-        }
-
-        return 0.0;
-    }
-
-    template<typename TypeTag>
-    double
-    WellInterface<TypeTag>::
-    getTHPConstraint(const SummaryState& summaryState) const
-    {
-        if (dynamic_thp_limit_) {
-            return *dynamic_thp_limit_;
-        }
-        if (well_ecl_.isInjector()) {
-            const auto& controls = well_ecl_.injectionControls(summaryState);
-            return controls.thp_limit;
-        }
-
-        if (well_ecl_.isProducer( )) {
-            const auto& controls = well_ecl_.productionControls(summaryState);
-            return controls.thp_limit;
-        }
-
-        return 0.0;
-
-    }
-
-
-
-
-
-
     template<typename TypeTag>
     bool
     WellInterface<TypeTag>::
@@ -572,18 +258,6 @@ namespace Opm
 
 
 
-    template<typename TypeTag>
-    bool
-    WellInterface<TypeTag>::
-    underPredictionMode() const
-    {
-        return well_ecl_.predictionMode();
-    }
-
-
-
-
-
     template<typename TypeTag>
     bool
     WellInterface<TypeTag>::
@@ -670,7 +344,7 @@ namespace Opm
              const auto& controls = well.productionControls(summaryState);
              const double vfp_ref_depth = vfp_properties_->getProd()->getTable(controls.vfp_table_number).getDatumDepth();
              const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
-             return vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState), getALQ(well_state)) - dp;
+             return vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, this->getTHPConstraint(summaryState), this->getALQ(well_state)) - dp;
          }
          else {
              OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER for well " + name(), deferred_logger);
@@ -681,13 +355,6 @@ namespace Opm
     }
 
 
-    template<typename TypeTag>
-    double
-    WellInterface<TypeTag>::
-    getALQ(const WellState &well_state) const
-    {
-        return well_state.getALQ(name());
-    }
 
 
 
@@ -992,36 +659,6 @@ namespace Opm
 
 
 
-    template<typename TypeTag>
-    void
-    WellInterface<TypeTag>::
-    updateWellTestStatePhysical(const WellState& /* well_state */,
-                                const double simulation_time,
-                                const bool write_message_to_opmlog,
-                                WellTestState& well_test_state,
-                                DeferredLogger& deferred_logger) const
-    {
-        if (!isOperable()) {
-            if (well_test_state.hasWellClosed(name(), WellTestConfig::Reason::ECONOMIC) ||
-                well_test_state.hasWellClosed(name(), WellTestConfig::Reason::PHYSICAL) ) {
-                // Already closed, do nothing.
-            } else {
-                well_test_state.closeWell(name(), WellTestConfig::Reason::PHYSICAL, simulation_time);
-                if (write_message_to_opmlog) {
-                    const std::string action = well_ecl_.getAutomaticShutIn() ? "shut" : "stopped";
-                    const std::string msg = "Well " + name()
-                        + " will be " + action + " as it can not operate under current reservoir conditions.";
-                    deferred_logger.info(msg);
-                }
-            }
-        }
-
-    }
-
-
-
-
-
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::
@@ -1245,50 +882,6 @@ namespace Opm
 
 
 
-
-
-    template<typename TypeTag>
-    void
-    WellInterface<TypeTag>::
-    setRepRadiusPerfLength(const std::vector<int>& cartesian_to_compressed)
-    {
-        const int nperf = number_of_perforations_;
-
-        perf_rep_radius_.clear();
-        perf_length_.clear();
-        bore_diameters_.clear();
-
-        perf_rep_radius_.reserve(nperf);
-        perf_length_.reserve(nperf);
-        bore_diameters_.reserve(nperf);
-
-        // COMPDAT handling
-        const auto& connectionSet = well_ecl_.getConnections();
-        CheckDistributedWellConnections checker(well_ecl_, parallel_well_info_);
-        for (size_t c=0; c<connectionSet.size(); c++) {
-            const auto& connection = connectionSet.get(c);
-            const int cell =
-                cartesian_to_compressed[connection.global_index()];
-            if (connection.state() != Connection::State::OPEN || cell >= 0)
-            {
-                checker.connectionFound(c);
-            }
-            if (connection.state() == Connection::State::OPEN) {
-
-                if (cell >= 0) {
-                    double radius = connection.rw();
-                    double re = connection.re(); // area equivalent radius of the grid block
-                    double perf_length = connection.connectionLength(); // the length of the well perforation
-                    const double repR = std::sqrt(re * radius);
-                    perf_rep_radius_.push_back(repR);
-                    perf_length_.push_back(perf_length);
-                    bore_diameters_.push_back(2. * radius);
-                }
-            }
-        }
-        checker.checkAllConnectionsFound();
-    }
-
     template<typename TypeTag>
     double
     WellInterface<TypeTag>::scalingFactor(const int phaseIdx) const
@@ -1311,47 +904,6 @@ namespace Opm
 
 
 
-
-    template<typename TypeTag>
-    bool
-    WellInterface<TypeTag>::isVFPActive(DeferredLogger& deferred_logger) 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 (isProducer()) { // producer
-            const int table_id = well_ecl_.vfp_table_number();
-            if (table_id <= 0) {
-                return false;
-            } else {
-                if (vfp_properties_->getProd()->hasTable(table_id)) {
-                    return true;
-                } else {
-                    OPM_DEFLOG_THROW(std::runtime_error, "VFPPROD table " << std::to_string(table_id) << " is specfied,"
-                              << " for well " << name() << ", while we could not access it during simulation", deferred_logger);
-                }
-            }
-
-        } else { // injector
-            const int table_id = well_ecl_.vfp_table_number();
-            if (table_id <= 0) {
-                return false;
-            } else {
-                if (vfp_properties_->getInj()->hasTable(table_id)) {
-                    return true;
-                } else {
-                    OPM_DEFLOG_THROW(std::runtime_error, "VFPINJ table " << std::to_string(table_id) << " is specfied,"
-                              << " for well " << name() << ", while we could not access it during simulation", deferred_logger);
-                }
-            }
-        }
-    }
-
-
-
-
-
     template<typename TypeTag>
     bool
     WellInterface<TypeTag>::
@@ -1393,22 +945,6 @@ namespace Opm
         }
     }
 
-    template<typename TypeTag>
-    void
-    WellInterface<TypeTag>::closeCompletions(WellTestState& wellTestState)
-    {
-        const auto& connections = well_ecl_.getConnections();
-        int perfIdx = 0;
-        for (const auto& connection : connections) {
-            if (connection.state() == Connection::State::OPEN) {
-                if (wellTestState.hasCompletion(name(), connection.complnum())) {
-                    well_index_[perfIdx] = 0.0;
-                }
-                perfIdx++;
-            }
-        }
-    }
-
     template<typename TypeTag>
     void
     WellInterface<TypeTag>::
@@ -2004,15 +1540,6 @@ namespace Opm
         return scaling_factor;
     }
 
-    template<typename TypeTag>
-    bool
-    WellInterface<TypeTag>::
-    isOperable() const {
-        return operability_status_.isOperable();
-    }
-
-
-
 
 
     template <typename TypeTag>