Implement multi-segment well handling in COMPTRAJ

Author: verveerpj

CMakeLists_files.cmake

@@ -735,6 +735,7 @@ if(ENABLE_ECL_OUTPUT)
           tests/SPE1CASE1B.DATA
           tests/SPE1CASE1_WELTRAJ.DATA
           tests/SPE1CASE1_WELTRAJ_2.DATA
+          tests/SPE1CASE1_WELTRAJ_MSW.DATA
           tests/props_spe1case1b.inc
           tests/SPE9_CP_PACKED.DATA
           tests/SOFR_TEST.DATA

opm/input/eclipse/Schedule/HandlerContext.cpp

@@ -83,6 +83,13 @@ void HandlerContext::compsegs_handled(const std::string& well_name)
     }
 }
 
+void HandlerContext::comptraj_handled(const std::string& well_name)
+{
+    if (comptraj_wells) {
+        comptraj_wells->insert(well_name);
+    }
+}
+
 ScheduleState& HandlerContext::state()
 {
     return schedule_.snapshots[currentStep];

opm/input/eclipse/Schedule/HandlerContext.hpp

@@ -68,6 +68,7 @@ class HandlerContext
     /// \param wpimult_global_factor_ Global well production index multipliers
     /// \param welsegs_wells_ All wells with a WELSEGS entry for checks.
     /// \param compsegs_wells_ All wells with a COMPSEGS entry for checks.
+    /// \param comptraj_wells_ All wells with a COMPTRAJ entry for checks.
     HandlerContext(Schedule& schedule,
                    const ScheduleBlock& block_,
                    const DeckKeyword& keyword_,
@@ -81,7 +82,8 @@ class HandlerContext
                    const std::unordered_map<std::string, double>* target_wellpi_,
                    std::unordered_map<std::string, double>& wpimult_global_factor_,
                    WelSegsSet* welsegs_wells_,
-                   std::set<std::string>* compsegs_wells_)
+                   std::set<std::string>* compsegs_wells_,
+                   std::set<std::string>* comptraj_wells_)
         : block(block_)
         , keyword(keyword_)
         , currentStep(currentStep_)
@@ -94,6 +96,7 @@ class HandlerContext
         , target_wellpi(target_wellpi_)
         , welsegs_wells(welsegs_wells_)
         , compsegs_wells(compsegs_wells_)
+        , comptraj_wells(comptraj_wells_)
         , sim_update(sim_update_)
         , schedule_(schedule)
     {}
@@ -122,6 +125,9 @@ class HandlerContext
     /// \brief Mark that the well occured in a COMPSEGS keyword.
     void compsegs_handled(const std::string& well_name);
 
+    /// \brief Mark that the well occured in a COMPTRAJ keyword.
+    void comptraj_handled(const std::string& well_name);
+
     //! \brief Set exit code.
     void setExitCode(int code);
 
@@ -208,6 +214,7 @@ class HandlerContext
     const std::unordered_map<std::string, double>* target_wellpi{nullptr};
     WelSegsSet* welsegs_wells{nullptr};
     std::set<std::string>* compsegs_wells{nullptr};
+    std::set<std::string>* comptraj_wells{nullptr};
     SimulatorUpdate* sim_update{nullptr};
     Schedule& schedule_;
 };

opm/input/eclipse/Schedule/MSW/Compsegs.cpp

@@ -310,14 +310,11 @@ namespace {
 
 
     std::pair<WellConnections, WellSegments>
-    processCOMPSEGS(const DeckKeyword& compsegs,
+    processCOMPSEGSfromRecords(const std::vector<Record>& compsegs_vector,
                     const WellConnections& input_connections,
                     const WellSegments& input_segments,
-                    const ScheduleGrid& grid,
-                    const ParseContext& parseContext,
-                    ErrorGuard& errors)
+                    const ScheduleGrid& grid)
     {
-            const auto& compsegs_vector = Compsegs::compsegsFromCOMPSEGSKeyword( compsegs, input_segments, grid, parseContext, errors);
             WellSegments new_segment_set = input_segments;
             WellConnections new_connection_set = input_connections;
 
@@ -351,6 +348,43 @@ namespace {
                                    WellSegments( std::move(new_segment_set)));
     }
 
+    std::pair<WellConnections, WellSegments>
+    processCOMPSEGS(const DeckKeyword& compsegs,
+                    const WellConnections& input_connections,
+                    const WellSegments& input_segments,
+                    const ScheduleGrid& grid,
+                    const ParseContext& parseContext,
+                    ErrorGuard& errors)
+    {
+        const auto& compsegs_vector = Compsegs::compsegsFromCOMPSEGSKeyword(compsegs, input_segments, grid, parseContext, errors);
+        return processCOMPSEGSfromRecords(compsegs_vector, input_connections, input_segments, grid);
+    }
+
+    std::pair<WellConnections, WellSegments>
+    processCOMPSEGS(std::vector<std::tuple<double, double, std::array<int, 3>>>& segments_md_and_ijk,
+                    const WellSegments& segments,
+                    const WellConnections& input_connections,
+                    const WellSegments& input_segments,
+                    const ScheduleGrid& grid)
+    {
+        std::vector<Record> compsegs;
+
+        for (auto [startMD, endMD, ijk] : segments_md_and_ijk) {
+            // Defaulted values:
+            auto direction = Connection::Direction::X;
+            const double center_depth = 0.0;
+            int segment_number = 0;
+            const int branch = 1;
+
+            std::size_t seqIndex = compsegs.size();
+            compsegs.emplace_back(
+                ijk[0], ijk[1], ijk[2], branch, startMD, endMD, direction, center_depth, segment_number, seqIndex);
+        }
+
+        processCOMPSEGS__(compsegs, segments);
+        return processCOMPSEGSfromRecords(compsegs, input_connections, input_segments, grid);
+    }
+
 namespace {
     // Duplicated from Well.cpp
     Connection::Order order_from_int(int int_value) {

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

@@ -58,6 +58,13 @@ namespace Compsegs {
                     const ParseContext& parseContext,
                     ErrorGuard& errors);
 
+    std::pair<WellConnections, WellSegments>
+    processCOMPSEGS(std::vector<std::tuple<double, double, std::array<int, 3>>>& segments_md_and_ijk,
+                    const WellSegments& segments,
+                    const WellConnections& input_connections,
+                    const WellSegments& input_segments,
+                    const ScheduleGrid& grid);
+
 
     std::pair<WellConnections, WellSegments>
     rstUpdate(const RestartIO::RstWell& rst_well,

opm/input/eclipse/Schedule/MSW/WelSegsSet.cpp

@@ -58,6 +58,24 @@ WelSegsSet::difference(const std::set<std::string>& compsegs,
     return difference;
 }
 
+std::vector<WelSegsSet::Entry> 
+WelSegsSet::intersection(const std::set<std::string>& wells1,
+                         const std::set<std::string>& wells2) const
+{
+    std::vector<Entry> intersection;
+    intersection.reserve(entries_.size());
+
+    for (const auto& entry : entries_) {
+        if (std::count(wells1.begin(), wells1.end(), std::get<0>(entry)) &&
+            std::count(wells2.begin(), wells2.end(), std::get<0>(entry))) {
+            intersection.push_back(entry);
+        }
+    }
+
+    return intersection;
+}
+
+
 bool WelSegsSet::PairComp::
 operator()(const Entry& pair, const std::string& str) const
 {

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

@@ -40,6 +40,8 @@ class WelSegsSet
 
     std::vector<Entry> difference(const std::set<std::string>& compsegs,
                                   const std::vector<Well>& wells) const;
+    std::vector<Entry> intersection(const std::set<std::string>& compsegs,
+                                    const std::set<std::string>& comptraj) const;
 
 private:
     struct PairComp

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

@@ -63,6 +63,14 @@ namespace Opm {
     }
 
 
+    WellSegments::WellSegments(const std::string &wname,
+        std::vector<std::pair<double, double>>& lengths_and_depths,
+        double diameter, const UnitSystem& unit_system)
+    {
+        this->addWellSegmentsFromLengthsAndDepths(wname, lengths_and_depths, diameter, unit_system);
+    }
+
+
     WellSegments WellSegments::serializationTestObject()
     {
         WellSegments result;
@@ -185,6 +193,44 @@ namespace Opm {
     }
 
 
+    void WellSegments::addWellSegmentsFromLengthsAndDepths(const std::string &wname,
+                                                           std::vector<std::pair<double, double>>& lengths_and_depths,
+                                                           double diameter, const UnitSystem& unit_system)
+    {
+        // Only LengthDepth::ABS is supported:
+        assert (this->m_length_depth_type == WellSegments::LengthDepth::ABS);
+
+        const int branchID = 1;  // Only main branch for now.
+        
+        const double roughness = 0.0;  // Defaulted: ROUGHNESS in WELSEGS.
+        const double area = M_PI * diameter * diameter / 4.0;
+        const double volume = Segment::invalidValue();
+        
+        // Add segments:
+        int segmentID = 2;
+        for (auto [length, depth]: lengths_and_depths) {
+            this->addSegment(
+                segmentID, branchID, segmentID - 1, depth, length, diameter,
+                roughness, area, volume, true, 0.0, 0.0
+            );           
+            segmentID += 1;
+        }
+
+        // Fix inlets:
+        for (const auto& segment : this->m_segments) {
+            const int outlet_segment = segment.outletSegment();
+            if (outlet_segment <= 0) { // no outlet segment
+                continue;
+            }
+
+            const int outlet_segment_index = segment_number_to_index[outlet_segment];
+            m_segments[outlet_segment_index].addInletSegment(segment.segmentNumber());
+        }
+
+        this->process(wname, unit_system, WellSegments::LengthDepth::ABS, this->depthTopSegment(), this->lengthTopSegment());
+    }
+
+
     void WellSegments::loadWELSEGS(const DeckKeyword& welsegsKeyword, const UnitSystem& unit_system)
     {
         // For the first record, which provides the information for the top
@@ -204,6 +250,8 @@ namespace Opm {
         const auto nodeX_top = record1.getItem("TOP_X").getSIDouble(0);
         const auto nodeY_top = record1.getItem("TOP_Y").getSIDouble(0);
 
+        this->m_length_depth_type = length_depth_type;
+
         // The main branch is 1 instead of 0.  The segment number for top
         // segment is also 1.
         if (length_depth_type == LengthDepth::INC) {

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

@@ -71,6 +71,12 @@ namespace Opm {
         WellSegments(CompPressureDrop compDrop,
                      const std::vector<Segment>& segments);
         void loadWELSEGS( const DeckKeyword& welsegsKeyword, const UnitSystem& unit_system);
+        explicit WellSegments(const std::string &wname,
+                              std::vector<std::pair<double, double>>& lengths_and_depths,
+                              double diameter, const UnitSystem& unit_system);
+        void addWellSegmentsFromLengthsAndDepths(const std::string &wname,
+                                                 std::vector<std::pair<double, double>>& lengths_and_depths,
+                                                 double diameter, const UnitSystem& unit_system);
 
         static WellSegments serializationTestObject();
 
@@ -120,6 +126,12 @@ namespace Opm {
             serializer(m_comp_pressure_drop);
             serializer(m_segments);
             serializer(segment_number_to_index);
+            serializer(m_length_depth_type);
+        }
+        
+        LengthDepth getLengthDepthType() const 
+        {
+            return this->m_length_depth_type;
         }
 
     private:
@@ -151,6 +163,9 @@ namespace Opm {
         // the mapping from the segment number to the
         // storage index in the vector
         std::map<int, int> segment_number_to_index{};
+
+        // The length/depth type, incremental or absolute
+        LengthDepth m_length_depth_type{LengthDepth::ABS};
     };
 }
 

opm/input/eclipse/Schedule/Schedule.cpp

@@ -430,12 +430,13 @@ namespace Opm {
                                  const std::unordered_map<std::string, double>* target_wellpi,
                                  std::unordered_map<std::string, double>& wpimult_global_factor,
                                  WelSegsSet* welsegs_wells,
-                                 std::set<std::string>* compsegs_wells)
+                                 std::set<std::string>* compsegs_wells,
+                                 std::set<std::string>* comptraj_wells)
     {
         HandlerContext handlerContext { *this, block, keyword, grid, currentStep,
                                         matches, action_mode,
                                         parseContext, errors, sim_update, target_wellpi,
-                                        wpimult_global_factor, welsegs_wells, compsegs_wells};
+                                        wpimult_global_factor, welsegs_wells, compsegs_wells, comptraj_wells};
 
         if (!KeywordHandlers::getInstance().handleKeyword(handlerContext)) {
             OpmLog::warning(fmt::format("No handler registered for keyword {} "
@@ -580,30 +581,47 @@ class ScheduleLogger
 
 namespace
 {
-/// \brief Check whether each MS well has COMPSEGS entry andissue error if not.
-/// \param welsegs All wells with a WELSEGS entry together with the location.
-/// \param compsegs All wells with a COMPSEGS entry
-void check_compsegs_consistency(const Opm::WelSegsSet& welsegs,
-                                const std::set<std::string>& compsegs,
-                                const std::vector<::Opm::Well>& wells)
-{
-    const auto difference = welsegs.difference(compsegs, wells);
 
-    if (!difference.empty()) {
+void report_welsegs_error(const std::vector<Opm::WelSegsSet::Entry>& segments, std::string format_str) 
+{
+    if (!segments.empty()) {
         std::string well_str = "well";
-        if (difference.size() > 1) {
+        if (segments.size() > 1) {
             well_str.append("s");
         }
         well_str.append(":");
 
-        for(const auto& [name, location] : difference) {
+        for(const auto& [name, location] : segments) {
             well_str.append(fmt::format("\n   {} in {} at line {}",
                                         name, location.filename, location.lineno));
         }
-        auto msg = fmt::format("Missing COMPSEGS keyword for the following multisegment {}.", well_str);
-        throw Opm::OpmInputError(msg, std::get<1>(difference[0]));
+        auto msg = fmt::format(format_str, well_str);
+        throw Opm::OpmInputError(msg, std::get<1>(segments[0]));
     }
 }
+
+/// \brief Check whether each MS well has COMPSEGS entry andissue error if not.
+/// \param welsegs All wells with a WELSEGS entry together with the location.
+/// \param compsegs All wells with a COMPSEGS entry
+/// \param comptraj All wells with a COMPTRAJ entry
+void check_compsegs_and_comptraj_consistency(const Opm::WelSegsSet& welsegs,
+                                             const std::set<std::string>& compsegs,
+                                             const std::set<std::string>& comptraj,
+                                             const std::vector<::Opm::Well>& wells)
+{
+    std::set<std::string> comp = compsegs;
+    comp.insert(comptraj.begin(), comptraj.end());
+    const auto difference = welsegs.difference(comp, wells);
+    report_welsegs_error(
+        difference,"Missing COMPSEGS keyword for the following multisegment {}."
+    );
+
+    const auto intersection = welsegs.intersection(compsegs, comptraj);
+    report_welsegs_error(
+        intersection, "Overlapping COMPSEGS and COMPTRAJ keywords for the following multisegment {}."
+    );
+}
+
 }// end anonymous namespace
 
 namespace Opm
@@ -666,6 +684,7 @@ void Schedule::iterateScheduleSection(std::size_t load_start, std::size_t load_e
         }
 
         std::set<std::string> compsegs_wells;
+        std::set<std::string> comptraj_wells;
         WelSegsSet welsegs_wells;
 
         const auto matches = Action::Result { false }.matches();
@@ -759,12 +778,13 @@ void Schedule::iterateScheduleSection(std::size_t load_start, std::size_t load_e
                                     target_wellpi,
                                     wpimult_global_factor,
                                     &welsegs_wells,
-                                    &compsegs_wells);
+                                    &compsegs_wells,
+                                    &comptraj_wells);
                 keyword_index++;
             }
 
             this->updateICDScalingFactors();
-            check_compsegs_consistency(welsegs_wells, compsegs_wells, this->getWells(report_step));
+            check_compsegs_and_comptraj_consistency(welsegs_wells, compsegs_wells, comptraj_wells, this->getWells(report_step));
             this->applyGlobalWPIMULT(wpimult_global_factor);
             this->end_report(report_step);