Merge pull request #3498 from bska/dynamic-segdims-restart

Dimension Segment Related Restart Arrays According to Dynamic Size

opm/input/eclipse/Schedule/MSW/WellSegments.hpp

@@ -77,6 +77,8 @@ namespace Opm {
 
         std::size_t size() const;
         bool empty() const;
+        int maxSegmentID() const;
+        int maxBranchID() const;
         double depthTopSegment() const;
         double lengthTopSegment() const;
         double volumeTopSegment() const;

opm/input/eclipse/Schedule/Well/Well.hpp

@@ -387,6 +387,8 @@ public:
     const std::vector<const Connection *> getConnections(int completion) const;
     const WellConnections& getConnections() const;
     const WellSegments& getSegments() const;
+    int maxSegmentID() const;
+    int maxBranchID() const;
 
     const WellProductionProperties& getProductionProperties() const;
     const WellInjectionProperties& getInjectionProperties() const;

src/opm/input/eclipse/Schedule/MSW/WellSegments.cpp

@@ -28,11 +28,19 @@
 #include <opm/input/eclipse/Deck/DeckKeyword.hpp>
 #include <opm/input/eclipse/Deck/DeckRecord.hpp>
 
+#include <algorithm>
 #include <cassert>
 #include <cmath>
-#include <map>
+#include <cstddef>
 #include <iterator>
+#include <map>
+#include <numeric>
+#include <set>
+#include <stdexcept>
+#include <string>
 #include <unordered_set>
+#include <utility>
+#include <vector>
 
 #ifdef _WIN32
 #define _USE_MATH_DEFINES
@@ -76,6 +84,26 @@ namespace Opm {
         return this->m_segments.empty();
     }
 
+    int WellSegments::maxSegmentID() const
+    {
+        return std::accumulate(this->m_segments.begin(),
+                               this->m_segments.end(), 0,
+            [](const int maxID, const Segment& seg)
+        {
+            return std::max(maxID, seg.segmentNumber());
+        });
+    }
+
+    int WellSegments::maxBranchID() const
+    {
+        return std::accumulate(this->m_segments.begin(),
+                               this->m_segments.end(), 0,
+            [](const int maxID, const Segment& seg)
+        {
+            return std::max(maxID, seg.branchNumber());
+        });
+    }
+
     const Segment& WellSegments::topSegment() const {
         return this->m_segments[0];
     }

src/opm/input/eclipse/Schedule/Well/Well.cpp

@@ -1122,6 +1122,20 @@ const WellSegments& Well::getSegments() const {
         throw std::logic_error("Asked for segment information in not MSW well: " + this->name());
 }
 
+int Well::maxSegmentID() const
+{
+    return (this->segments == nullptr)
+        ? 0
+        : this->segments->maxSegmentID();
+}
+
+int Well::maxBranchID() const
+{
+    return (this->segments == nullptr)
+        ? 0
+        : this->segments->maxBranchID();
+}
+
 const Well::WellInjectionProperties& Well::getInjectionProperties() const {
     return *this->injection;
 }

src/opm/output/eclipse/CreateInteHead.cpp

@@ -237,6 +237,38 @@ namespace {
         });
     }
 
+    int maxNumSegments(const ::Opm::Schedule& sched,
+                       const std::size_t      report_step,
+                       const std::size_t      lookup_step)
+    {
+        if (report_step == 0) { return 0; }
+
+        const auto& wnames = sched.wellNames(lookup_step);
+
+        return std::accumulate(std::begin(wnames), std::end(wnames), 0,
+            [&sched, lookup_step](const int m, const std::string& wname) -> int
+        {
+            // maxSegmentID() returns 0 for standard (non-MS) wells.
+            return std::max(m, sched.getWell(wname, lookup_step).maxSegmentID());
+        });
+    }
+
+    int maxNumLateralBranches(const ::Opm::Schedule& sched,
+                              const std::size_t      report_step,
+                              const std::size_t      lookup_step)
+    {
+        if (report_step == 0) { return 0; }
+
+        const auto& wnames = sched.wellNames(lookup_step);
+
+        return std::accumulate(std::begin(wnames), std::end(wnames), 0,
+            [&sched, lookup_step](const int m, const std::string& wname) -> int
+        {
+            // maxBranchID() returns 0 for standard (non-MS) wells.
+            return std::max(m, sched.getWell(wname, lookup_step).maxBranchID());
+        });
+    }
+
     Opm::RestartIO::InteHEAD::WellTableDim
     getWellTableDims(const int              nwgmax,
                      const int              ngmax,
@@ -384,13 +416,18 @@ namespace {
     {
         const auto& wsd = rspec.wellSegmentDimensions();
 
+        const auto numMSW = numMultiSegWells(sched, report_step, lookup_step);
+        const auto maxNumSeg = maxNumSegments(sched, report_step, lookup_step);
+        const auto maxNumBr = maxNumLateralBranches(sched, report_step, lookup_step);
+
         return {
-            numMultiSegWells(sched, report_step, lookup_step),
-            wsd.maxSegmentedWells(),
-            wsd.maxSegmentsPerWell(),
-            wsd.maxLateralBranchesPerWell(),
+            numMSW,
+            std::max(numMSW, wsd.maxSegmentedWells()),
+            std::max(maxNumSeg, wsd.maxSegmentsPerWell()),
+            std::max(maxNumBr, wsd.maxLateralBranchesPerWell()),
             22,           // Number of entries per segment in ISEG (2017.2)
-            Opm::RestartIO::InteHEAD::numRsegElem(rspec.phases()) + 8*num_water_tracer, // Number of entries per segment in RSEG
+            Opm::RestartIO::InteHEAD::numRsegElem(rspec.phases())
+               + 8*num_water_tracer, // Number of entries per segment in RSEG
             10            // Number of entries per segment in ILBR (2017.2)
         };
     }