Use std::variant<> to hold ICD alternatives

opm/parser/eclipse/EclipseState/Schedule/MSW/Segment.hpp

@@ -21,13 +21,14 @@
 #define SEGMENT_HPP_HEADER_INCLUDED
 
 #include <optional>
-#include <memory>
+#include <variant>
 #include <vector>
 
-namespace Opm {
+#include <opm/parser/eclipse/EclipseState/Schedule/MSW/Valve.hpp>
-    class SICD;
+#include <opm/parser/eclipse/EclipseState/Schedule/MSW/SICD.hpp>
-    class Valve;
 
+
+namespace Opm {
     namespace RestartIO {
         struct RstSegment;
     }
@@ -35,6 +36,29 @@ namespace Opm {
 
 namespace Opm {
 
+    /*
+      The current serialization of std::variant<> requires that all the types in
+      the variant have serializeOp() method, that is why this RegularSegment
+      type is introduced, ideally the icd variant should just have
+      std::monostate to represent the regular non ICD segment.
+    */
+    struct RegularSegment : std::monostate {
+
+        template<class Serializer>
+        void serializeOp(Serializer& serializer) {
+        }
+
+        static RegularSegment serializeObject() {
+            return RegularSegment();
+        }
+
+
+        bool operator==(const RegularSegment& ) {
+            return true;
+        }
+    };
+
+
     class Segment {
     public:
 
@@ -51,7 +75,7 @@ namespace Opm {
         Segment(const Segment& src, double new_depth, double new_length);
         Segment(const Segment& src, double new_volume);
         Segment(int segment_number_in, int branch_in, int outlet_segment_in, double length_in, double depth_in,
-                double internal_diameter_in, double roughness_in, double cross_area_in, double volume_in, bool data_ready_in, SegmentType segment_type_in);
+                double internal_diameter_in, double roughness_in, double cross_area_in, double volume_in, bool data_ready_in);
         Segment(const RestartIO::RstSegment& rst_segment);
 
         static Segment serializeObject();
@@ -80,7 +104,7 @@ namespace Opm {
         bool operator==( const Segment& ) const;
         bool operator!=( const Segment& ) const;
 
-        const std::shared_ptr<SICD>& spiralICD() const;
+        const SICD& spiralICD() const;
         const Valve* valve() const;
 
         void updatePerfLength(double perf_length);
@@ -89,6 +113,21 @@ namespace Opm {
         void updateValve(const Valve& valve);
         void addInletSegment(const int segment_number);
 
+        bool isRegular() const
+        {
+            return std::holds_alternative<RegularSegment>(this->m_icd);
+        }
+
+        inline bool isSpiralICD() const
+        {
+            return std::holds_alternative<SICD>(this->m_icd);
+        }
+
+        inline bool isValve() const
+        {
+            return std::holds_alternative<Valve>(this->m_icd);
+        }
+
         template<class Serializer>
         void serializeOp(Serializer& serializer)
         {
@@ -104,9 +143,7 @@ namespace Opm {
             serializer(m_volume);
             serializer(m_data_ready);
             serializer(m_perf_length);
-            serializer(m_segment_type);
+            serializer(m_icd);
-            serializer(m_spiral_icd);
-            serializer(m_valve);
         }
 
     private:
@@ -157,17 +194,7 @@ namespace Opm {
         bool m_data_ready;
 
         std::optional<double> m_perf_length;
-        // indicate the type of the segment
+        std::variant<RegularSegment, SICD, Valve> m_icd;
-        // regular, spiral ICD, or Valve.
-        SegmentType m_segment_type;
-
-        // information related to SpiralICD. It is nullptr for segments are not
-        // spiral ICD type
-        std::shared_ptr<SICD> m_spiral_icd;
-
-        // information related to sub-critical valve. It is nullptr for segments are not
-        // of type of Valve
-        std::shared_ptr<Valve> m_valve;
 
         // We are not handling the length of segment projected onto the X-axis and Y-axis.
         // They are not used in the simulations and we are not supporting the plotting.
@@ -175,20 +202,6 @@ namespace Opm {
         // while they are not supported by the keyword at the moment.
     };
 
-    inline bool isRegular(const Segment& segment)
-    {
-        return segment.segmentType() == Segment::SegmentType::REGULAR;
-    }
-
-    inline bool isSpiralICD(const Segment& segment)
-    {
-        return segment.segmentType() == Segment::SegmentType::SICD;
-    }
-
-    inline bool isValve(const Segment& segment)
-    {
-        return segment.segmentType() == Segment::SegmentType::VALVE;
-    }
 }
 
 #endif

src/opm/output/eclipse/AggregateMSWData.cpp

@@ -424,8 +424,8 @@ namespace {
                 VectorItems::ISeg::index;
 
             const auto& sicd = segment.spiralICD();
-            iSeg[baseIndex + Ix::ICDScalingMode] = sicd->methodFlowScaling();
+            iSeg[baseIndex + Ix::ICDScalingMode] = sicd.methodFlowScaling();
-            iSeg[baseIndex + Ix::ICDOpenShutFlag] = sicd->ecl_status();
+            iSeg[baseIndex + Ix::ICDOpenShutFlag] = sicd.ecl_status();
         }
 
         template <class ISegArray>
@@ -433,7 +433,7 @@ namespace {
                                               const std::size_t   baseIndex,
                                               ISegArray&          iSeg)
         {
-            if (isSpiralICD(segment)) {
+            if (segment.isSpiralICD()) {
                 assignSpiralICDCharacteristics(segment, baseIndex, iSeg);
             }
         }
@@ -474,7 +474,7 @@ namespace {
                     iSeg[iS + 8] = seg_reorder[ind];
 
                     iSeg[iS + Ix::SegmentType] = segment.ecl_type_id();
-                    if (! isRegular(segment)) {
+                    if (! segment.isRegular()) {
                         assignSegmentTypeCharacteristics(segment, iS, iSeg);
                     }
                 }
@@ -577,27 +577,27 @@ namespace {
             const auto& sicd = segment.spiralICD();
 
             rSeg[baseIndex + Ix::DeviceBaseStrength] =
-                usys.from_si(M::icd_strength, sicd->strength());
+                usys.from_si(M::icd_strength, sicd.strength());
 
             rSeg[baseIndex + Ix::CalibrFluidDensity] =
-                usys.from_si(M::density, sicd->densityCalibration());
+                usys.from_si(M::density, sicd.densityCalibration());
 
             rSeg[baseIndex + Ix::CalibrFluidViscosity] =
-                usys.from_si(M::viscosity, sicd->viscosityCalibration());
+                usys.from_si(M::viscosity, sicd.viscosityCalibration());
 
-            rSeg[baseIndex + Ix::CriticalWaterFraction] = sicd->criticalValue();
+            rSeg[baseIndex + Ix::CriticalWaterFraction] = sicd.criticalValue();
 
             rSeg[baseIndex + Ix::TransitionRegWidth] =
-                sicd->widthTransitionRegion();
+                sicd.widthTransitionRegion();
 
             rSeg[baseIndex + Ix::MaxEmulsionRatio] =
-                sicd->maxViscosityRatio();
+                sicd.maxViscosityRatio();
 
             rSeg[baseIndex + Ix::MaxValidFlowRate] =
-                usys.from_si(M::geometric_volume_rate, sicd->maxAbsoluteRate());
+                usys.from_si(M::geometric_volume_rate, sicd.maxAbsoluteRate());
 
             rSeg[baseIndex + Ix::ICDLength] =
-                usys.from_si(M::length, sicd->length());
+                usys.from_si(M::length, sicd.length());
         }
 
         template <class RSegArray>
@@ -606,11 +606,11 @@ namespace {
                                               const int                baseIndex,
                                               RSegArray&               rSeg)
         {
-            if (isSpiralICD(segment)) {
+            if (segment.isSpiralICD()) {
                 assignSpiralICDCharacteristics(segment, usys, baseIndex, rSeg);
             }
 
-            if (isValve(segment)) {
+            if (segment.isValve()) {
                 assignValveCharacteristics(segment, usys, baseIndex, rSeg);
             }
         }
@@ -760,7 +760,7 @@ namespace {
                     rSeg[iS + Ix::item110] = 1.0;
                     rSeg[iS + Ix::item111] = 1.0;
 
-                    if (! isRegular(segment)) {
+                    if (! segment.isRegular()) {
                         assignSegmentTypeCharacteristics(segment, units, iS, rSeg);
                     }
                 }

src/opm/parser/eclipse/EclipseState/Schedule/MSW/Segment.cpp

@@ -71,10 +71,9 @@ namespace {
         m_roughness(if_invalid_value(rst_segment.roughness)),
         m_cross_area(if_invalid_value(rst_segment.area)),
         m_volume(rst_segment.volume),
-        m_data_ready(true),
+        m_data_ready(true)
-        m_segment_type(rst_segment.segment_type)
     {
-        if (this->m_segment_type == SegmentType::SICD) {
+        if (rst_segment.segment_type == SegmentType::SICD) {
             double scalingFactor = -1;  // The scaling factor will be and updated from the simulator.
 
             SICD icd(rst_segment.base_strength,
@@ -92,7 +91,7 @@ namespace {
             this->updateSpiralICD(icd);
         }
 
-        if (this->m_segment_type == SegmentType::VALVE) {
+        if (rst_segment.segment_type == SegmentType::VALVE) {
             /*
               These three variables are currently not stored in the restart
               file; here we initialize with the default values, but if they have
@@ -121,7 +120,7 @@ namespace {
 
     Segment::Segment(int segment_number_in, int branch_in, int outlet_segment_in, double length_in, double depth_in,
                      double internal_diameter_in, double roughness_in, double cross_area_in,
-                     double volume_in, bool data_ready_in, SegmentType segment_type_in)
+                     double volume_in, bool data_ready_in)
     : m_segment_number(segment_number_in),
       m_branch(branch_in),
       m_outlet_segment(outlet_segment_in),
@@ -131,8 +130,7 @@ namespace {
       m_roughness(roughness_in),
       m_cross_area(cross_area_in),
       m_volume(volume_in),
-      m_data_ready(data_ready_in),
+      m_data_ready(data_ready_in)
-      m_segment_type(segment_type_in)
     {
     }
 
@@ -170,10 +168,7 @@ namespace {
       result.m_cross_area = 10.0;
       result.m_volume = 11.0;
       result.m_data_ready = true;
-      result.m_segment_type = SegmentType::SICD;
+      result.m_icd = SICD::serializeObject();
-      result.m_spiral_icd = std::make_shared<SICD>(SICD::serializeObject());
-      result.m_valve = std::make_shared<Valve>(Valve::serializeObject());
-
       return result;
   }
 
@@ -228,7 +223,16 @@ namespace {
     }
 
     Segment::SegmentType Segment::segmentType() const {
-        return m_segment_type;
+        if (this->isRegular())
+            return SegmentType::REGULAR;
+
+        if (this->isSpiralICD())
+            return SegmentType::SICD;
+
+        if (this->isValve())
+            return SegmentType::VALVE;
+
+        throw std::logic_error("This just should not happen ");
     }
 
     const std::vector<int>& Segment::inletSegments() const {
@@ -254,6 +258,7 @@ namespace {
             && this->m_cross_area        == rhs.m_cross_area
             && this->m_volume            == rhs.m_volume
             && this->m_perf_length       == rhs.m_perf_length
+            && this->m_icd               == rhs.m_icd
             && this->m_data_ready        == rhs.m_data_ready;
     }
 
@@ -262,46 +267,48 @@ namespace {
     }
 
     void Segment::updateSpiralICD(const SICD& spiral_icd) {
-        m_segment_type = SegmentType::SICD;
+        this->m_icd = spiral_icd;
-        m_spiral_icd = std::make_shared<SICD>(spiral_icd);
     }
 
-    const std::shared_ptr<SICD>& Segment::spiralICD() const {
+    const SICD& Segment::spiralICD() const {
-        return m_spiral_icd;
+        return std::get<SICD>(this->m_icd);
     }
 
     void Segment::updateValve(const Valve& valve) {
         if (valve.pipeAdditionalLength() < 0)
             throw std::logic_error("Bug in handling of pipe length for valves");
 
+    void Segment::updateValve(const Valve& input_valve, const double segment_length) {
         // we need to update some values for the vale
-        auto valve_ptr = std::make_shared<Valve>(valve);
+        auto valve = input_valve;
 
-        if (valve_ptr->pipeDiameter() < 0.) {
+        if (valve.pipeAdditionalLength() < 0.) { // defaulted for this
-            valve_ptr->setPipeDiameter(m_internal_diameter);
+            valve.setPipeAdditionalLength(segment_length);
+        }
+
+        if (valve.pipeDiameter() < 0.) {
+            valve.setPipeDiameter(m_internal_diameter);
         } else {
-            this->m_internal_diameter = valve_ptr->pipeDiameter();
+            this->m_internal_diameter = valve.pipeDiameter();
         }
 
-        if (valve_ptr->pipeRoughness() < 0.) {
+        if (valve.pipeRoughness() < 0.) {
-            valve_ptr->setPipeRoughness(m_roughness);
+            valve.setPipeRoughness(m_roughness);
         } else {
-            this->m_roughness = valve_ptr->pipeRoughness();
+            this->m_roughness = valve.pipeRoughness();
         }
 
-        if (valve_ptr->pipeCrossArea() < 0.) {
+        if (valve.pipeCrossArea() < 0.) {
-            valve_ptr->setPipeCrossArea(m_cross_area);
+            valve.setPipeCrossArea(m_cross_area);
         } else {
-            this->m_cross_area = valve_ptr->pipeCrossArea();
+            this->m_cross_area = valve.pipeCrossArea();
         }
 
-        if (valve_ptr->conMaxCrossArea() < 0.) {
+        if (valve.conMaxCrossArea() < 0.) {
-            valve_ptr->setConMaxCrossArea(valve_ptr->pipeCrossArea());
+            valve.setConMaxCrossArea(valve.pipeCrossArea());
         }
 
-        this->m_valve = valve_ptr;
+        this->m_icd= valve;
-
-        m_segment_type = SegmentType::VALVE;
     }
 
 
@@ -324,11 +331,12 @@ namespace {
 
 
     const Valve* Segment::valve() const {
-        return m_valve.get();
+        return &std::get<Valve>(this->m_icd);
+        //return m_valve.get();
     }
 
     int Segment::ecl_type_id() const {
-        switch (this->m_segment_type) {
+        switch (this->segmentType()) {
         case SegmentType::REGULAR:
             return -1;
         case SegmentType::SICD:

src/opm/parser/eclipse/EclipseState/Schedule/MSW/WellSegments.cpp

@@ -142,12 +142,12 @@ namespace Opm {
         if (length_depth_type == LengthDepth::INC) {
             m_segments.emplace_back( 1, 1, 0, 0., 0.,
                                      invalid_value, invalid_value, invalid_value,
-                                     volume_top, false , Segment::SegmentType::REGULAR);
+                                     volume_top, false);
 
         } else if (length_depth_type == LengthDepth::ABS) {
             m_segments.emplace_back( 1, 1, 0, length_top, depth_top,
                                      invalid_value, invalid_value, invalid_value,
-                                     volume_top, true , Segment::SegmentType::REGULAR);
+                                     volume_top, true);
         }
 
         // read all the information out from the DECK first then process to get all the required information
@@ -207,13 +207,13 @@ namespace Opm {
 
                 if (length_depth_type == LengthDepth::INC) {
                     m_segments.emplace_back( i, branch, outlet_segment, segment_length, depth_change,
-                                             diameter, roughness, area, volume, false , Segment::SegmentType::REGULAR);
+                                             diameter, roughness, area, volume, false);
                 } else if (i == segment2) {
                     m_segments.emplace_back( i, branch, outlet_segment, segment_length, depth_change,
-                                             diameter, roughness, area, volume, true , Segment::SegmentType::REGULAR);
+                                             diameter, roughness, area, volume, true);
                 } else {
                     m_segments.emplace_back( i, branch, outlet_segment, invalid_value, invalid_value,
-                                             diameter, roughness, area, volume, false , Segment::SegmentType::REGULAR);
+                                             diameter, roughness, area, volume, false);
                 }
             }
         }

tests/parser/MultisegmentWellTests.cpp

@@ -126,30 +126,30 @@ BOOST_AUTO_TEST_CASE(MultisegmentWellTest) {
     const auto& sicd_vector = it->second;
     BOOST_CHECK_EQUAL(1U, sicd_vector.size());
     const int segment_number = sicd_vector[0].first;
-    const Opm::SICD& sicd = sicd_vector[0].second;
+    const Opm::SICD& sicd0 = sicd_vector[0].second;
 
     BOOST_CHECK_EQUAL(8, segment_number);
 
     Opm::Segment segment = segment_set.getFromSegmentNumber(segment_number);
-    segment.updateSpiralICD(sicd);
+    segment.updateSpiralICD(sicd0);
 
     BOOST_CHECK(Opm::Segment::SegmentType::SICD==segment.segmentType());
 
-    const std::shared_ptr<Opm::SICD> sicd_ptr = segment.spiralICD();
+    auto sicd = segment.spiralICD();
-    BOOST_CHECK_GT(sicd_ptr->maxAbsoluteRate(), 1.e99);
+    BOOST_CHECK_GT(sicd.maxAbsoluteRate(), 1.e99);
-    BOOST_CHECK(sicd_ptr->status()==Opm::ICDStatus::SHUT);
+    BOOST_CHECK(sicd.status()==Opm::ICDStatus::SHUT);
     // 0.002 bars*day*day/Volume^2
-    BOOST_CHECK_EQUAL(sicd_ptr->strength(), 0.002*1.e5*86400.*86400.);
+    BOOST_CHECK_EQUAL(sicd.strength(), 0.002*1.e5*86400.*86400.);
-    BOOST_CHECK_EQUAL(sicd_ptr->length(), -0.7);
+    BOOST_CHECK_EQUAL(sicd.length(), -0.7);
-    BOOST_CHECK_EQUAL(sicd_ptr->densityCalibration(), 1000.25);
+    BOOST_CHECK_EQUAL(sicd.densityCalibration(), 1000.25);
     // 1.45 cp
-    BOOST_CHECK_EQUAL(sicd_ptr->viscosityCalibration(), 1.45 * 0.001);
+    BOOST_CHECK_EQUAL(sicd.viscosityCalibration(), 1.45 * 0.001);
-    BOOST_CHECK_EQUAL(sicd_ptr->criticalValue(), 0.6);
+    BOOST_CHECK_EQUAL(sicd.criticalValue(), 0.6);
-    BOOST_CHECK_EQUAL(sicd_ptr->widthTransitionRegion(), 0.05);
+    BOOST_CHECK_EQUAL(sicd.widthTransitionRegion(), 0.05);
-    BOOST_CHECK_EQUAL(sicd_ptr->maxViscosityRatio(), 5.0);
+    BOOST_CHECK_EQUAL(sicd.maxViscosityRatio(), 5.0);
-    BOOST_CHECK_EQUAL(sicd_ptr->methodFlowScaling(), -1);
+    BOOST_CHECK_EQUAL(sicd.methodFlowScaling(), -1);
     // the scaling factor has not been updated properly, so it will throw
-    BOOST_CHECK_THROW(sicd_ptr->scalingFactor(), std::runtime_error);
+    BOOST_CHECK_THROW(sicd.scalingFactor(), std::runtime_error);
 
     const int outlet_segment_number = segment.outletSegment();
     const double outlet_segment_length = segment_set.segmentLength(outlet_segment_number);
@@ -157,11 +157,11 @@ BOOST_AUTO_TEST_CASE(MultisegmentWellTest) {
     const Opm::Connection& connection = new_connection_set.getFromIJK(15, 0, 1);
     const auto& perf_range = connection.perf_range();
     const auto connection_length = perf_range->second - perf_range->first;
-    sicd_ptr->updateScalingFactor(outlet_segment_length, connection_length);
+    sicd.updateScalingFactor(outlet_segment_length, connection_length);
 
     // updated, so it should not throw
-    BOOST_CHECK_NO_THROW(sicd_ptr->scalingFactor());
+    BOOST_CHECK_NO_THROW(sicd.scalingFactor());
-    BOOST_CHECK_EQUAL(0.7, sicd_ptr->scalingFactor());
+    BOOST_CHECK_EQUAL(0.7, sicd.scalingFactor());
 
     BOOST_CHECK_EQUAL(7U, new_connection_set.size());