From 31ac5378cf499ecd3cfd5e17db0e0947ce5822da Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tosa@norceresearch.no>
Date: Thu, 19 Aug 2021 12:15:35 +0200
Subject: [PATCH] fix case where wells are shut behind our back

---
 ebos/ecltracermodel.hh                        | 14 ++++++++--
 opm/simulators/wells/BlackoilWellModel.hpp    |  3 +-
 .../wells/BlackoilWellModel_impl.hpp          | 28 ++++++++-----------
 3 files changed, 26 insertions(+), 19 deletions(-)

diff --git a/ebos/ecltracermodel.hh b/ebos/ecltracermodel.hh
index 925ac3cdda..89b7bdaf67 100644
--- a/ebos/ecltracermodel.hh
+++ b/ebos/ecltracermodel.hh
@@ -302,6 +302,11 @@ protected:
             if (well.getStatus() == Well::Status::SHUT)
                 continue;
 
+            if (!simulator_.problem().wellModel().hasWell(well.name()))
+                continue;
+
+            const auto& wellPtr = simulator_.problem().wellModel().getWell(well.name());
+
             std::vector<double> wtracer(tr.numTracer());
             for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
                 wtracer[tIdx] = well.getTracerProperties().getConcentration(this->tracerNames_[tr.idx_[tIdx]]);
@@ -318,7 +323,7 @@ protected:
                 cartesianCoordinate[2] = connection.getK();
                 const size_t cartIdx = simulator_.vanguard().cartesianIndex(cartesianCoordinate);
                 const int I = this->cartToGlobal_[cartIdx];
-                Scalar rate = simulator_.problem().wellModel().well(well.name())->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
+                Scalar rate = wellPtr->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
                 if (rate > 0) {
                     for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
                         tr.residual_[tIdx][I][0] -= rate*wtracer[tIdx];
@@ -390,6 +395,11 @@ protected:
             if (well.getStatus() == Well::Status::SHUT)
                 continue;
 
+            if (!simulator_.problem().wellModel().hasWell(well.name()))
+                continue;
+
+            const auto& wellPtr = simulator_.problem().wellModel().getWell(well.name());
+
             if (!well.isProducer()) //Injection rates already reported during assembly
                 continue;
 
@@ -406,7 +416,7 @@ protected:
                 cartesianCoordinate[2] = connection.getK();
                 const size_t cartIdx = simulator_.vanguard().cartesianIndex(cartesianCoordinate);
                 const int I = this->cartToGlobal_[cartIdx];
-                Scalar rate = simulator_.problem().wellModel().well(well.name())->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
+                Scalar rate = wellPtr->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
                 if (rate < 0) {
                     rateWellNeg += rate;
                     for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp
index 42c70c9fde..f1e3158272 100644
--- a/opm/simulators/wells/BlackoilWellModel.hpp
+++ b/opm/simulators/wells/BlackoilWellModel.hpp
@@ -210,7 +210,6 @@ namespace Opm {
 
 
             using WellInterfacePtr = std::shared_ptr<WellInterface<TypeTag> >;
-            WellInterfacePtr well(const std::string& wellName) const;
 
             using BlackoilWellModelGeneric::initFromRestartFile;
             void initFromRestartFile(const RestartValue& restartValues)
@@ -279,6 +278,8 @@ namespace Opm {
             void initPrimaryVariablesEvaluation() const;
             void updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls);
             WellInterfacePtr getWell(const std::string& well_name) const;
+            bool hasWell(const std::string& well_name) const;
+
             void initGliftEclWellMap(GLiftEclWells &ecl_well_map);
 
         protected:
diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp
index 80d58456f5..c7f2c43638 100644
--- a/opm/simulators/wells/BlackoilWellModel_impl.hpp
+++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp
@@ -511,22 +511,6 @@ namespace Opm {
 
 
 
-    template<typename TypeTag>
-    typename BlackoilWellModel<TypeTag>::WellInterfacePtr
-    BlackoilWellModel<TypeTag>::
-    well(const std::string& wellName) const
-    {
-        for (const auto& well : well_container_) {
-            if (well->name() == wellName) {
-                return well;
-            }
-        }
-        OPM_THROW(std::invalid_argument, "The well with name " + wellName + " is not in the well Container");
-        return nullptr;
-    }
-
-
-
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
@@ -1645,6 +1629,18 @@ namespace Opm {
         return *well;
     }
 
+    template<typename TypeTag>
+    bool
+    BlackoilWellModel<TypeTag>::
+    hasWell(const std::string& well_name) const
+    {
+        return std::any_of(well_container_.begin(), well_container_.end(),
+            [&well_name](const WellInterfacePtr& elem) -> bool
+        {
+            return elem->name() == well_name;
+        });
+    }
+