Assign tracer data to XWEL

opm/output/eclipse/AggregateWellData.hpp

@@ -60,6 +60,7 @@ namespace Opm { namespace RestartIO { namespace Helpers {
                                      const std::vector<int>& 	 inteHead);
 
         void captureDynamicWellData(const Opm::Schedule&        sched,
+                                    const TracerConfig&       tracer,
                                     const std::size_t           sim_step,
                                     const Opm::data::Wells&     xw,
                                     const Opm::SummaryState&    smry);

opm/output/eclipse/VectorItems/well.hpp

@@ -247,6 +247,8 @@ namespace Opm { namespace RestartIO { namespace Helpers { namespace VectorItems
 
             WatVoidPrRate = 122, // Well's voidage production rate
             GasVoidPrRate = 123, // Well's voidage production rate
+
+            TracerOffset  = 130, // Tracer data start at this index
         };
     } // XWell
 

opm/parser/eclipse/EclipseState/TracerConfig.hpp

@@ -103,6 +103,7 @@ public:
     static TracerConfig serializeObject();
 
     size_t size() const;
+    bool empty() const;
 
     const std::vector<TracerEntry>::const_iterator begin() const;
     const std::vector<TracerEntry>::const_iterator end() const;

src/opm/output/eclipse/AggregateWellData.cpp

@@ -911,7 +911,75 @@ namespace {
         }
 
         template <class XWellArray>
-        void dynamicContrib(const ::Opm::Well&        well,
+        void assignTracerData(const Opm::TracerConfig& tracers,
+                              const Opm::Tracers& tracer_dims,
+                              const Opm::SummaryState& smry,
+                              const Opm::Well& well,
+                              XWellArray& xWell)
+        {
+            if (tracers.empty())
+                return;
+
+            using Ix = ::Opm::RestartIO::Helpers::VectorItems::XWell::index;
+            std::fill(xWell.begin() + Ix::TracerOffset, xWell.end(), 0);
+
+
+            for (std::size_t tracer_index=0; tracer_index < tracers.size(); tracer_index++) {
+                const auto& tracer = tracers[tracer_index];
+                std::size_t output_index = Ix::TracerOffset + tracer_index;
+                if (well.isInjector()) {
+                    const auto& wtir = smry.get_well_var(well.name(), fmt::format("WTIR{}", tracer.name), 0);
+                    xWell[output_index] = -wtir;
+                } else {
+                    const auto& wtpr = smry.get_well_var(well.name(), fmt::format("WTPR{}", tracer.name), 0);
+                    xWell[output_index] = wtpr;
+                }
+                output_index++;
+            }
+
+
+            for (std::size_t tracer_index=0; tracer_index < tracers.size(); tracer_index++) {
+                const auto& tracer = tracers[tracer_index];
+                std::size_t output_index = Ix::TracerOffset + tracer_dims.water_tracers() + tracer_index;
+                if (well.isProducer()) {
+                    const auto& wtpr = smry.get_well_var(well.name(), fmt::format("WTPT{}", tracer.name), 0);
+                    xWell[output_index] = wtpr;
+                }
+                output_index++;
+            }
+
+            for (std::size_t tracer_index=0; tracer_index < tracers.size(); tracer_index++) {
+                const auto& tracer = tracers[tracer_index];
+                std::size_t output_index = Ix::TracerOffset + 2*tracer_dims.water_tracers() + tracer_index;
+                if (well.isInjector()) {
+                    const auto& wtir = smry.get_well_var(well.name(), fmt::format("WTIT{}", tracer.name), 0);
+                    xWell[output_index] = -wtir;
+                }
+            }
+
+            for (std::size_t n=0; n < 2; n++) {
+                for (std::size_t tracer_index=0; tracer_index < tracers.size(); tracer_index++) {
+                    const auto& tracer = tracers[tracer_index];
+                    std::size_t output_index = Ix::TracerOffset + (3 + n)*tracer_dims.water_tracers() + tracer_index;
+                    const auto& wtic = smry.get_well_var(well.name(), fmt::format("WTIC{}", tracer.name), 0);
+                    const auto& wtpc = smry.get_well_var(well.name(), fmt::format("WTPC{}", tracer.name), 0);
+
+                    if (std::abs(wtic) > 0)
+                        xWell[output_index] = wtic;
+                    else
+                        xWell[output_index] = wtpc;
+                }
+            }
+
+            std::size_t output_index = Ix::TracerOffset + 5*tracer_dims.water_tracers();
+            xWell[output_index] = 0;
+            xWell[output_index + 1] = 0;
+        }
+
+        template <class XWellArray>
+        void dynamicContrib(const ::Opm::Well&         well,
+                            const Opm::TracerConfig&   tracers,
+                            const Opm::Tracers&        tracer_dims,
                             const ::Opm::SummaryState& smry,
                             XWellArray&                xWell)
         {
@@ -942,6 +1010,7 @@ namespace {
                 }
             }
             assignCumulatives(well.name(), smry, xWell);
+            assignTracerData(tracers, tracer_dims, smry, well, xWell);
         }
     } // XWell
 
@@ -1064,6 +1133,7 @@ captureDeclaredWellData(const Schedule&   sched,
 void
 Opm::RestartIO::Helpers::AggregateWellData::
 captureDynamicWellData(const Opm::Schedule&        sched,
+                       const TracerConfig&         tracers,
                        const std::size_t           sim_step,
                        const Opm::data::Wells& xw,
                        const ::Opm::SummaryState&  smry)
@@ -1089,11 +1159,11 @@ captureDynamicWellData(const Opm::Schedule&        sched,
     });
 
     // Dynamic contributions to XWEL array.
-    wellLoop(wells, sched, sim_step, [this, &smry]
+    wellLoop(wells, sched, sim_step, [this, &sched, &tracers, &smry]
         (const Well& well, const std::size_t wellID) -> void
     {
         auto xwell = this->xWell_[wellID];
 
-        XWell::dynamicContrib(well, smry, xwell);
+        XWell::dynamicContrib(well, tracers, sched.runspec().tracers(), smry, xwell);
     });
 }

src/opm/output/eclipse/CreateInteHead.cpp

@@ -577,11 +577,7 @@ createInteHead(const EclipseState& es,
     const auto& rdim  = tdim.getRegdims();
     const auto& rckcfg = es.getSimulationConfig().rock_config();
     auto num_water_tracer = es.runspec().tracers().water_tracers();
-    int nxwelz_tracer_shift;
+    int nxwelz_tracer_shift = num_water_tracer*5 + 2 * (num_water_tracer > 0);
-    if (num_water_tracer > 1)
-        nxwelz_tracer_shift = 7 + (num_water_tracer - 1)*5;
-    else
-        nxwelz_tracer_shift = 7*num_water_tracer;
 
     const auto ih = InteHEAD{}
         .dimensions         (grid.getNXYZ())

src/opm/output/eclipse/RestartIO.cpp

@@ -392,7 +392,7 @@ namespace {
     {
         auto wellData = Helpers::AggregateWellData(ih);
         wellData.captureDeclaredWellData(schedule, tracers, sim_step, action_state, wtest_state, sumState, ih);
-        wellData.captureDynamicWellData(schedule, sim_step, wells, sumState);
+        wellData.captureDynamicWellData(schedule, tracers, sim_step, wells, sumState);
 
         rstFile.write("IWEL", wellData.getIWell());
         rstFile.write("SWEL", wellData.getSWell());

src/opm/parser/eclipse/EclipseState/TracerConfig.cpp

@@ -168,6 +168,11 @@ size_t TracerConfig::size() const {
     return this->tracers.size();
 }
 
+bool TracerConfig::empty() const {
+    return this->tracers.empty();
+}
+
+
 const std::vector<TracerConfig::TracerEntry>::const_iterator TracerConfig::begin() const {
     return this->tracers.begin();
 }

tests/test_AggregateWellData.cpp

@@ -941,7 +941,7 @@ BOOST_AUTO_TEST_CASE (Dynamic_Well_Data_Step1)
     auto awd = Opm::RestartIO::Helpers::AggregateWellData{ih.value};
     Opm::WellTestState wtest_state;
 
-    awd.captureDynamicWellData(simCase.sched, rptStep, xw, smry);
+    awd.captureDynamicWellData(simCase.sched, simCase.es.tracer(), rptStep, xw, smry);
 
     // IWEL (OP_1)
     {
@@ -1085,7 +1085,7 @@ BOOST_AUTO_TEST_CASE (Dynamic_Well_Data_Step2)
     auto awd = Opm::RestartIO::Helpers::AggregateWellData{ih.value};
     Opm::WellTestState wtest_state;
 
-    awd.captureDynamicWellData(simCase.sched, rptStep, xw, smry);
+    awd.captureDynamicWellData(simCase.sched, simCase.es.tracer(), rptStep, xw, smry);
 
     // IWEL (OP_1) -- closed producer
     {
@@ -1279,7 +1279,7 @@ BOOST_AUTO_TEST_CASE(WELL_POD) {
 
     auto wellData = Opm::RestartIO::Helpers::AggregateWellData(ih);
     wellData.captureDeclaredWellData(simCase.sched, simCase.es.tracer(), sim_step, action_state, wtest_state, sumState, ih);
-    wellData.captureDynamicWellData(simCase.sched, sim_step, xw , sumState);
+    wellData.captureDynamicWellData(simCase.sched, simCase.es.tracer(), sim_step, xw , sumState);
 
     auto connectionData = Opm::RestartIO::Helpers::AggregateConnectionData(ih);
     connectionData.captureDeclaredConnData(simCase.sched, simCase.grid, units, xw , sumState, sim_step);

tests/test_rst.cpp

@@ -298,7 +298,7 @@ BOOST_AUTO_TEST_CASE(State_test) {
 
     auto wellData = Opm::RestartIO::Helpers::AggregateWellData(ih);
     wellData.captureDeclaredWellData(simCase.sched, simCase.es.tracer(), sim_step, action_state, wtest_state, sumState, ih);
-    wellData.captureDynamicWellData(simCase.sched, sim_step, {} , sumState);
+    wellData.captureDynamicWellData(simCase.sched, simCase.es.tracer(), sim_step, {} , sumState);
 
     auto connectionData = Opm::RestartIO::Helpers::AggregateConnectionData(ih);
     connectionData.captureDeclaredConnData(simCase.sched, simCase.grid, units, {} , sumState, sim_step);