Merge pull request #2864 from joakim-hove/tracer-zero

Pad XCON, XGRP and RSEG vectors with tracer zero

opm/output/eclipse/InteHEAD.hpp

@@ -31,7 +31,7 @@ namespace Opm {
 class EclipseGrid;
 class EclipseState;
 class UnitSystem;
-
+class Phases;
 }
 
 namespace Opm { namespace RestartIO {
@@ -228,6 +228,7 @@ namespace Opm { namespace RestartIO {
         InteHEAD& whistControlMode(int mode);
         InteHEAD& liftOptParam(int in_enc);
 
+        static int numRsegElem(const Opm::Phases& phase);
         const std::vector<int>& data() const
         {
             return this->data_;

opm/output/eclipse/VectorItems/connection.hpp

@@ -86,6 +86,8 @@ namespace Opm { namespace RestartIO { namespace Helpers { namespace VectorItems
             ResVRate     = 49,  // Reservoir voidage rate
             VoidPrTotal  = 50,  // Total cumulative reservoir voidage volume production
             VoidInjTotal = 51,  // Total cumulative reservoir voidage volume injection
+
+            TracerOffset = 58,  // Tracer data starts after this index
         };
     } // XConn
 }}}} // Opm::RestartIO::Helpers::VectorItems

opm/output/eclipse/VectorItems/group.hpp

@@ -158,6 +158,8 @@ namespace Opm { namespace RestartIO { namespace Helpers { namespace VectorItems
                                    // production (observed/historical rates)
             HistGasInjTotal = 144, // Group's total cumulative gas injection
                                    // (observed/historical rates)
+
+            TracerOffset = 180,    // Tracer data starts here
         };
     } // XGroup
 

src/opm/output/eclipse/AggregateConnectionData.cpp

@@ -295,6 +295,9 @@ namespace {
             xConn[Ix::OilRate_Copy] = xConn[Ix::OilRate];
             xConn[Ix::GasRate_Copy] = xConn[Ix::GasRate];
             xConn[Ix::WaterRate_Copy] = xConn[Ix::WaterRate];
+
+            // Pad the connection array with tracer values
+            std::fill(xConn.begin() + Ix::TracerOffset, xConn.end(), 0);
         }
     } // XConn
 } // Anonymous

src/opm/output/eclipse/AggregateGroupData.cpp

@@ -975,6 +975,8 @@ void dynamicContrib(const std::vector<std::string>&      restart_group_keys,
     xGrp[Ix::VoidPrGuideRate_2] = xGrp[Ix::VoidPrGuideRate];
 
     xGrp[Ix::WatInjGuideRate_2] = xGrp[Ix::WatInjGuideRate];
+
+    std::fill(xGrp.begin() + Ix::TracerOffset, xGrp.end(), 0);
 }
 } // XGrp
 

src/opm/output/eclipse/AggregateMSWData.cpp

@@ -18,7 +18,7 @@
 */
 
 #include <opm/output/eclipse/AggregateMSWData.hpp>
-
+#include <opm/output/eclipse/InteHEAD.hpp>
 #include <opm/output/eclipse/VectorItems/msw.hpp>
 
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
@@ -722,7 +722,19 @@ namespace {
         }
 
         template <class RSegArray>
-        void staticContrib_useMSW(const Opm::Well&           well,
+        void assignTracerData(std::size_t segment_offset,
+                              const Opm::Runspec& runspec,
+                              RSegArray& rSeg)
+        {
+            auto tracer_offset = segment_offset + Opm::RestartIO::InteHEAD::numRsegElem(runspec.phases());
+            auto tracer_end    = tracer_offset + runspec.tracers().water_tracers() * 8;
+
+            std::fill(rSeg.begin() + tracer_offset, rSeg.begin() + tracer_end, 0.0);
+        }
+
+        template <class RSegArray>
+        void staticContrib_useMSW(const Opm::Runspec&         runspec,
+                                  const Opm::Well&            well,
                                   const std::vector<int>&     inteHead,
                                   const Opm::EclipseGrid&     grid,
                                   const Opm::UnitSystem&      units,
@@ -763,70 +775,76 @@ namespace {
                     sSFR = getSegmentSetFlowRates(welSegSet, wRatesIt->second.connections, welConns, units);
                 }
 
-                std::string stringSegNum = std::to_string(segment0.segmentNumber());
+                auto get = [&smry, &wname](const std::string& vector, const std::string& segment_nr)
-                auto get = [&smry, &wname, &stringSegNum](const std::string& vector)
                 {
-                    // 'stringSegNum' is one-based (1 .. #segments inclusive)
+                    const auto key = vector + ':' + wname + ':' + segment_nr;
-                    const auto key = vector + ':' + wname + ':' + stringSegNum;
+                    return smry.get(key, 0.0);
-                    return smry.has(key) ? smry.get(key) : 0.0;
                 };
 
-                auto iS = (segment0.segmentNumber() - 1)*noElmSeg;
+
                 // Treat the top segment individually
-                rSeg[iS + Ix::DistOutlet]      = units.from_si(M::length, welSegSet.lengthTopSegment());
+                {
-                rSeg[iS + Ix::OutletDepthDiff] = units.from_si(M::length, welSegSet.depthTopSegment());
+                    const int segNumber = segment0.segmentNumber();
-                rSeg[iS + Ix::SegVolume]       = volFromLengthUnitConv*welSegSet.volumeTopSegment();
+                    const auto& segment_string = std::to_string(segNumber);
-                rSeg[iS + Ix::DistBHPRef]      = rSeg[iS + Ix::DistOutlet];
+                    auto iS = (segNumber - 1)*noElmSeg;
-                rSeg[iS + Ix::DepthBHPRef]     = rSeg[iS + Ix::OutletDepthDiff];
+
-                //
+                    rSeg[iS + Ix::DistOutlet]      = units.from_si(M::length, welSegSet.lengthTopSegment());
-                // branch according to whether multisegment well calculations are switched on or not
+                    rSeg[iS + Ix::OutletDepthDiff] = units.from_si(M::length, welSegSet.depthTopSegment());
+                    rSeg[iS + Ix::SegVolume]       = volFromLengthUnitConv*welSegSet.volumeTopSegment();
+                    rSeg[iS + Ix::DistBHPRef]      = rSeg[iS + Ix::DistOutlet];
+                    rSeg[iS + Ix::DepthBHPRef]     = rSeg[iS + Ix::OutletDepthDiff];
+                    //
+                    // branch according to whether multisegment well calculations are switched on or not
 
 
-                if (haveWellRes && wRatesIt->second.segments.size() < 2) {
+                    if (haveWellRes && wRatesIt->second.segments.size() < 2) {
-                    // Note: Segment flow rates and pressure from 'smry' have correct
+                        // Note: Segment flow rates and pressure from 'smry' have correct
-                    // output units and sign conventions.
+                        // output units and sign conventions.
-                    temp_o = sSFR.sofr[0];
+                        temp_o = sSFR.sofr[0];
-                    temp_w = sSFR.swfr[0]*0.1;
+                        temp_w = sSFR.swfr[0]*0.1;
-                    temp_g = sSFR.sgfr[0]*gfactor;
+                        temp_g = sSFR.sgfr[0]*gfactor;
-                    //Item 12 Segment pressure - use well flow bhp
+                        //Item 12 Segment pressure - use well flow bhp
-                    rSeg[iS + Ix::Pressure] = smry.get_well_var(wname, "WBHP", 0);
+                        rSeg[iS + Ix::Pressure] = smry.get_well_var(wname, "WBHP", 0);
+                    }
+                    else {
+                        // Note: Segment flow rates and pressure from 'smry' have correct
+                        // output units and sign conventions.
+                        temp_o = get("SOFR", segment_string);
+                        temp_w = get("SWFR", segment_string)*0.1;
+                        temp_g = get("SGFR", segment_string)*gfactor;
+                        //Item 12 Segment pressure
+                        rSeg[iS + Ix::Pressure] = get("SPR", segment_string);
+                    }
+
+                    rSeg[iS + Ix::TotFlowRate] = temp_o + temp_w + temp_g;
+                    rSeg[iS + Ix::WatFlowFract] = (std::abs(temp_w) > 0) ? temp_w / rSeg[8] : 0.;
+                    rSeg[iS + Ix::GasFlowFract] = (std::abs(temp_g) > 0) ? temp_g / rSeg[8] : 0.;
+
+
+                    rSeg[iS + Ix::item31] = rSeg[iS + Ix::WatFlowFract];
+
+                    //  value is 1. based on tests on several data sets
+                    rSeg[iS + Ix::item40] = 1.;
+
+                    rSeg[iS + Ix::flowFractionOilDensityExponent]       = 1.0;
+                    rSeg[iS + Ix::flowFractionWaterDensityExponent]     = 1.0;
+                    rSeg[iS + Ix::flowFractionGasDensityExponent]       = 1.0;
+                    rSeg[iS + Ix::flowFractionOilViscosityExponent]     = 1.0;
+                    rSeg[iS + Ix::flowFractionWaterViscosityExponent]   = 1.0;
+                    rSeg[iS + Ix::flowFractionGasViscosityExponent]     = 1.0;
+
+                    assignTracerData(iS, runspec, rSeg);
                 }
-                else {
-                    // Note: Segment flow rates and pressure from 'smry' have correct
-                    // output units and sign conventions.
-                    temp_o = get("SOFR");
-                    temp_w = get("SWFR")*0.1;
-                    temp_g = get("SGFR")*gfactor;
-                    //Item 12 Segment pressure
-                    rSeg[iS + Ix::Pressure] = get("SPR");
-                }
-
-                rSeg[iS + Ix::TotFlowRate] = temp_o + temp_w + temp_g;
-                rSeg[iS + Ix::WatFlowFract] = (std::abs(temp_w) > 0) ? temp_w / rSeg[8] : 0.;
-                rSeg[iS + Ix::GasFlowFract] = (std::abs(temp_g) > 0) ? temp_g / rSeg[8] : 0.;
-
-
-                rSeg[iS + Ix::item31] = rSeg[iS + Ix::WatFlowFract];
-
-                //  value is 1. based on tests on several data sets
-                rSeg[iS + Ix::item40] = 1.;
-
-                rSeg[iS + Ix::flowFractionOilDensityExponent]       = 1.0;
-                rSeg[iS + Ix::flowFractionWaterDensityExponent]     = 1.0;
-                rSeg[iS + Ix::flowFractionGasDensityExponent]       = 1.0;
-                rSeg[iS + Ix::flowFractionOilViscosityExponent]     = 1.0;
-                rSeg[iS + Ix::flowFractionWaterViscosityExponent]   = 1.0;
-                rSeg[iS + Ix::flowFractionGasViscosityExponent]     = 1.0;
 
                 //Treat subsequent segments
                 for (std::size_t segIndex = 1; segIndex < welSegSet.size(); segIndex++) {
                     const auto& segment = welSegSet[segIndex];
                     const auto& outlet_segment = welSegSet.getFromSegmentNumber( segment.outletSegment() );
                     const int segNumber = segment.segmentNumber();
-                    stringSegNum = std::to_string(segNumber);
+                    const auto& segment_string = std::to_string(segNumber);
 
                     // set the elements of the rSeg array
-                    iS = (segNumber - 1)*noElmSeg;
+                    auto iS = (segNumber - 1)*noElmSeg;
                     rSeg[iS + Ix::DistOutlet]      = units.from_si(M::length, (segment.totalLength() - outlet_segment.totalLength()));
                     rSeg[iS + Ix::OutletDepthDiff] = units.from_si(M::length, (segment.depth() - outlet_segment.depth()));
                     rSeg[iS + Ix::SegDiam]         = units.from_si(M::length, (segment.internalDiameter()));
@@ -850,11 +868,11 @@ namespace {
                     else {
                         // Note: Segment flow rates and pressure from 'smry' have correct
                         // output units and sign conventions.
-                        temp_o = get("SOFR");
+                        temp_o = get("SOFR", segment_string);
-                        temp_w = get("SWFR")*0.1;
+                        temp_w = get("SWFR", segment_string)*0.1;
-                        temp_g = get("SGFR")*gfactor;
+                        temp_g = get("SGFR", segment_string)*gfactor;
                         //Item 12 Segment pressure
-                        rSeg[iS +  Ix::Pressure] = get("SPR");
+                        rSeg[iS +  Ix::Pressure] = get("SPR", segment_string);
                     }
 
                     rSeg[iS + Ix::TotFlowRate] = temp_o + temp_w + temp_g;
@@ -875,6 +893,8 @@ namespace {
                     if (! segment.isRegular()) {
                         assignSegmentTypeCharacteristics(segment, units, iS, rSeg);
                     }
+
+                    assignTracerData(iS, runspec, rSeg);
                 }
             }
             else {
@@ -1012,12 +1032,11 @@ captureDeclaredMSWData(const Schedule&         sched,
     }
     // Extract Contributions to RSeg Array
     {
-        MSWLoop(msw, [&units, &inteHead, &grid, &smry, this, &wr]
+        MSWLoop(msw, [&units, &inteHead, &sched, &grid, &smry, this, &wr]
             (const Well& well, const std::size_t mswID) -> void
         {
             auto rmsw = this->rSeg_[mswID];
-
+            RSeg::staticContrib_useMSW(sched.runspec(), well, inteHead, grid, units, smry, wr, rmsw);
-            RSeg::staticContrib_useMSW(well, inteHead, grid, units, smry, wr, rmsw);
         });
     }
     // Extract Contributions to ILBS Array

src/opm/output/eclipse/CreateInteHead.cpp

@@ -367,24 +367,6 @@ namespace {
         };
     }
 
-    int numRsegElem(const ::Opm::Phases& phase)
-    {
-        const auto nact = phase.active(::Opm::Phase::OIL)
-            + phase.active(::Opm::Phase::GAS)
-            + phase.active(::Opm::Phase::WATER);
-
-        switch (nact) {
-        case 1: return 126;
-        case 2: return 134;
-        case 3: return 146;
-        }
-
-        throw std::invalid_argument {
-            "NRSEGZ is not supported for " +
-            std::to_string(nact) +
-            " active phases"
-        };
-    }
 
     Opm::RestartIO::InteHEAD::WellSegDims
     getWellSegDims(const int              num_water_tracer,
@@ -401,7 +383,7 @@ namespace {
             wsd.maxSegmentsPerWell(),
             wsd.maxLateralBranchesPerWell(),
             22,           // Number of entries per segment in ISEG (2017.2)
-            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)
         };
     }

src/opm/output/eclipse/InteHEAD.cpp

@@ -27,6 +27,7 @@
 #include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/parser/eclipse/Units/UnitSystem.hpp>
+#include <opm/parser/eclipse/EclipseState/Runspec.hpp>
 
 #include <opm/common/utility/TimeService.hpp>
 
@@ -822,6 +823,24 @@ Opm::RestartIO::InteHEAD::netBalanceData(const NetBalanceDims& nwbaldim)
     return *this;
 }
 
+int Opm::RestartIO::InteHEAD::numRsegElem(const ::Opm::Phases& phase)
+{
+    const auto nact = phase.active(::Opm::Phase::OIL)
+        + phase.active(::Opm::Phase::GAS)
+        + phase.active(::Opm::Phase::WATER);
+
+    switch (nact) {
+    case 1: return 126;
+    case 2: return 134;
+    case 3: return 146;
+    }
+
+    throw std::invalid_argument {
+        "NRSEGZ is not supported for " +
+            std::to_string(nact) +
+            " active phases"
+            };
+}
 
 // =====================================================================