From 1c0296960887b251d69a53ff7b3befeb16e2685d Mon Sep 17 00:00:00 2001
From: Andreas Lauser <and@poware.org>
Date: Fri, 13 Nov 2015 13:24:39 +0100
Subject: [PATCH] do not iterate over all wells on DOFs which are not
 penetrated by any wells

for large problems with many wells the performance impact of this is
probably even measureable.

Also, this patch makes it possible to access the well manager outside
of the problem. In the normal case, this should be rarely needed, but
it can come in handy for debugging purposes.
---
 applications/ebos/eclproblem.hh     |  8 ++++++++
 applications/ebos/eclwellmanager.hh | 25 +++++++++++++++++++++----
 2 files changed, 29 insertions(+), 4 deletions(-)

diff --git a/applications/ebos/eclproblem.hh b/applications/ebos/eclproblem.hh
index bfccc7472..1074a3bfd 100644
--- a/applications/ebos/eclproblem.hh
+++ b/applications/ebos/eclproblem.hh
@@ -683,6 +683,14 @@ public:
             rate[eqIdx] /= this->model().dofTotalVolume(globalDofIdx);
     }
 
+    /*!
+     * \brief Returns a reference to the ECL well manager used by the problem.
+     *
+     * This can be used for inspecting wells outside of the problem.
+     */
+    const EclWellManager<TypeTag>& wellManager() const
+    { return wellManager_; }
+
 private:
     static bool enableEclOutput_()
     { return EWOMS_GET_PARAM(TypeTag, bool, EnableEclOutput); }
diff --git a/applications/ebos/eclwellmanager.hh b/applications/ebos/eclwellmanager.hh
index 59d805217..27046bd64 100644
--- a/applications/ebos/eclwellmanager.hh
+++ b/applications/ebos/eclwellmanager.hh
@@ -174,7 +174,7 @@ public:
         computeWellCompletionsMap_(episodeIdx, wellCompMap);
 
         if (wasRestarted || wellTopologyChanged_(eclState, episodeIdx)) {
-            updateWellTopology_(episodeIdx, wellCompMap);
+            updateWellTopology_(episodeIdx, wellCompMap, gridDofIsPenetrated_);
         }
 
         // set those parameters of the wells which do not change the topology of the
@@ -367,6 +367,12 @@ public:
         return wellNameToIndex_.find( wellName ) != wellNameToIndex_.end();
     }
 
+    /*!
+     * \brief Returns true iff a given degree of freedom is currently penetrated by any well.
+     */
+    bool gridDofIsPenetrated(int globalDofIdx) const
+    { return gridDofIsPenetrated_[globalDofIdx]; }
+
     /*!
      * \brief Given a well name, return the corresponding index.
      *
@@ -548,8 +554,10 @@ public:
                                  int dofIdx,
                                  int timeIdx) const
     {
-        for (int eqIdx = 0; eqIdx < numEq; ++ eqIdx)
-            q[eqIdx] = 0.0;
+        q = 0.0;
+
+        if (!gridDofIsPenetrated(context.globalSpaceIndex(dofIdx, timeIdx)))
+            return;
 
         RateVector wellRate;
 
@@ -690,7 +698,9 @@ protected:
         return false;
     }
 
-    void updateWellTopology_(int reportStepIdx, const WellCompletionsMap& wellCompletions)
+    void updateWellTopology_(int reportStepIdx,
+                             const WellCompletionsMap& wellCompletions,
+                             std::vector<bool>& gridDofIsPenetrated) const
     {
         auto& model = simulator_.model();
         const auto& gridManager = simulator_.gridManager();
@@ -704,6 +714,10 @@ protected:
 
         // tell the active wells which DOFs they contain
         const auto gridView = simulator_.gridManager().gridView();
+
+        gridDofIsPenetrated.resize(model.numGridDof());
+        std::fill(gridDofIsPenetrated.begin(), gridDofIsPenetrated.end(), false);
+
         ElementContext elemCtx(simulator_);
         auto elemIt = gridView.template begin</*codim=*/0>();
         const auto elemEndIt = gridView.template end</*codim=*/0>();
@@ -723,6 +737,8 @@ protected:
                     // it...
                     continue;
 
+                gridDofIsPenetrated[globalDofIdx] = true;
+
                 auto eclWell = wellCompletions.at(cartesianDofIdx).second;
                 eclWell->addDof(elemCtx, dofIdx);
 
@@ -868,6 +884,7 @@ protected:
     Simulator &simulator_;
 
     std::vector<std::shared_ptr<Well> > wells_;
+    std::vector<bool> gridDofIsPenetrated_;
     std::map<std::string, int> wellNameToIndex_;
     std::map<std::string, std::array<Scalar, numPhases> > wellTotalInjectedVolume_;
     std::map<std::string, std::array<Scalar, numPhases> > wellTotalProducedVolume_;