Support restart for TEMP (and gas tracers)

opm/input/eclipse/EclipseState/Runspec.cpp

@@ -598,6 +598,12 @@ Tracers Tracers::serializationTestObject() {
 int Tracers::water_tracers() const {
     return this->m_water_tracers;
 }
+int Tracers::gas_tracers() const {
+    return this->m_gas_tracers;
+}
+int Tracers::oil_tracers() const {
+    return this->m_oil_tracers;
+}
 
 
 Tracers::Tracers(const Deck& deck) {

opm/input/eclipse/EclipseState/Runspec.hpp

@@ -449,6 +449,8 @@ class Tracers {
 
     explicit Tracers(const Deck& );
     int water_tracers() const;
+    int gas_tracers() const;
+    int oil_tracers() const;
 
     template<class Serializer>
     void serializeOp(Serializer& serializer) {

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

@@ -465,12 +465,31 @@ Well::Well(const RestartIO::RstWell& rst_well,
 
         this->updateInjection(std::move(i));
 
+        const auto isTemp = (rst_well.inj_temperature < RestartIO::RstWell::UNDEFINED_VALUE);
+        std::size_t tracer_conc_index = 0;
+        if (isTemp) {
+            this->well_inj_temperature = rst_well.inj_temperature;
+            ++tracer_conc_index;
+        }
+
         if (!rst_well.tracer_concentration_injection.empty()) {
             auto tracer = std::make_shared<WellTracerProperties>(this->getTracerProperties());
-            for (std::size_t tracer_index = 0; tracer_index < tracer_config.size(); tracer_index++) {
+            for (std::size_t tracer_index = 0; tracer_index < tracer_config.size(); ++tracer_index) {
                 const auto& tname = tracer_config[tracer_index].name;
-                const auto concentration = rst_well.tracer_concentration_injection[tracer_index];
-                tracer->setConcentration(tname, concentration);
+                const auto phase = tracer_config[tracer_index].phase;
+                if (phase == Phase::WATER) {
+                    const auto concentration = rst_well.tracer_concentration_injection[tracer_conc_index];
+                    tracer->setConcentration(tname, concentration);
+                } else {
+                    const auto free_conc = rst_well.tracer_concentration_injection[tracer_conc_index];
+                    const auto sol_conc = rst_well.tracer_concentration_injection[++tracer_conc_index];
+                    if (WellType::gas_injector(this->wtype.ecl_wtype()) || WellType::oil_injector(this->wtype.ecl_wtype())) {
+                        tracer->setConcentration(tname, free_conc);
+                        if (sol_conc > 0.0) {
+                            OpmLog::warning(fmt::format("Well {}: Restoring a non-zero solution concentration of tracer {} is not yet supported.", rst_well.name, tname));
+                        }
+                    }
+                }
             }
             this->updateTracer(tracer);
         }

opm/io/eclipse/rst/well.cpp

@@ -181,14 +181,21 @@ Opm::RestartIO::RstWell::RstWell(const UnitSystem&  unit_system,
     water_void_rate(     unit_system.to_si(M::liquid_surface_volume, xwel[VI::XWell::WatVoidPrRate])),
     gas_void_rate(       unit_system.to_si(M::gas_surface_volume,    xwel[VI::XWell::GasVoidPrRate]))
 {
-    // For E100 it appears that +1 instead of -1 is written for group_controllable_flag when the 
-    // group control is aactive, so using this to correct active_control (where ind.ctrl. is written)
+    // For E100 it appears that +1 instead of -1 is written for group_controllable_flag when the
+    // group control is active, so using this to correct active_control (where ind.ctrl. is written)
     if (group_controllable_flag > 0) {
         active_control = VI::IWell::Value::WellCtrlMode::Group;
     }
 
-    for (std::size_t tracer_index = 0;
-         tracer_index < static_cast<std::size_t>(header.runspec.tracers().water_tracers());
+    // If the TEMP option is active, this is the first tracer
+    std::size_t tracer_index = 0;
+    const bool isTemp = header.runspec.temp();
+    if (isTemp) {
+        this->inj_temperature = unit_system.to_si(M::temperature, swel[VI::SWell::TracerOffset + tracer_index++]);
+    }
+    const auto& tracers = header.runspec.tracers();
+    for (const auto num_tracer_injconcs = tracers.water_tracers() + 2*(tracers.gas_tracers() + tracers.oil_tracers()) + tracer_index;
+         tracer_index < static_cast<std::size_t>(num_tracer_injconcs);
          ++tracer_index)
     {
         this->tracer_concentration_injection.push_back(swel[VI::SWell::TracerOffset + tracer_index]);

opm/io/eclipse/rst/well.hpp

@@ -129,6 +129,7 @@ struct RstWell
     float dfac_corr_coeff_a{};
     float dfac_corr_exponent_b{};
     float dfac_corr_exponent_c{};
+    float inj_temperature{UNDEFINED_VALUE};
     std::vector<float> tracer_concentration_injection;
 
     double oil_rate;

opm/output/eclipse/AggregateWellData.cpp

@@ -1080,14 +1080,22 @@ namespace {
         void assignTracerData(const Opm::TracerConfig& tracers,
                               const Opm::SummaryState& smry,
                               const std::string&       wname,
-                              SWellArray&              sWell)
+                              SWellArray&              sWell,
+                              const bool               isTemp = false)
         {
             auto output_index = static_cast<std::size_t>(VI::SWell::index::TracerOffset);
+            // Temperature tracer is first, if present
+            if (isTemp) sWell[output_index++] = smry.get_well_var(wname, "WTICHEA", 0.0);
 
             for (const auto& tracer : tracers) {
                 if (tracer.phase == Opm::Phase::WATER) {
                     sWell[output_index++] =
                         smry.get_well_var(wname, fmt::format("WTIC{}", tracer.name), 0.0);
+                } else {
+                    sWell[output_index++] =
+                        smry.get_well_var(wname, fmt::format("WTICF{}", tracer.name), 0.0);
+                    sWell[output_index++] =
+                        smry.get_well_var(wname, fmt::format("WTICS{}", tracer.name), 0.0);
                 }
             }
         }
@@ -1133,7 +1141,8 @@ namespace {
             assignDFactorCorrelation(well, units, sWell);
             assignEconomicLimits(well, swprop, sWell);
             assignWellTest(well.name(), sched, wtest_state, sim_step, swprop, sWell);
-            assignTracerData(tracers, smry, well.name(), sWell);
+
+            assignTracerData(tracers, smry, well.name(), sWell, sched.runspec().temp());
             assignBhpVfpAdjustment(well, swprop, sWell);
         }
     } // SWell
@@ -1340,70 +1349,90 @@ namespace {
                               const Opm::Tracers& tracer_dims,
                               const Opm::SummaryState& smry,
                               const Opm::Well& well,
-                              XWellArray& xWell)
+                              XWellArray& xWell,
+                              const bool isTemp)
         {
             if (tracers.empty() || tracer_dims.water_tracers() == 0)
                 return;
 
             using Ix = ::Opm::RestartIO::Helpers::VectorItems::XWell::index;
             std::fill(xWell.begin() + Ix::TracerOffset, xWell.end(), 0);
 
+            // For each vector in rate, prod_total, inj_total, inj_conc, prod_conc:
+            //    TEMP (if present), then each tracer in definition order (TRACER) [free then solution conc for HC tracers]
 
-            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;
+            auto output_index = static_cast<std::size_t>(Ix::TracerOffset);
+            // Rates
+            if (well.isInjector()) {
+                if (isTemp) xWell[output_index++] = -smry.get_well_var(well.name(), "WTIRHEA", 0);
+                for (const auto& tracer : tracers) {
+                    if (tracer.phase == Opm::Phase::WATER) {
+                        xWell[output_index++] = -smry.get_well_var(well.name(), fmt::format("WTIR{}", tracer.name), 0);
+                    } else {
+                        xWell[output_index++] = -smry.get_well_var(well.name(), fmt::format("WTIRF{}", tracer.name), 0);
+                        xWell[output_index++] = -smry.get_well_var(well.name(), fmt::format("WTIRS{}", tracer.name), 0);
+                    }
+                }
+            } else {
+                if (isTemp) xWell[output_index++] = smry.get_well_var(well.name(), "WTPRHEA", 0);
+                for (const auto& tracer : tracers) {
+                    if (tracer.phase == Opm::Phase::WATER) {
+                        xWell[output_index++] = smry.get_well_var(well.name(), fmt::format("WTPR{}", tracer.name), 0);
+                    } else {
+                        xWell[output_index++] = smry.get_well_var(well.name(), fmt::format("WTPRF{}", tracer.name), 0);
+                        xWell[output_index++] = smry.get_well_var(well.name(), fmt::format("WTPRS{}", tracer.name), 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_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;
+            // Production totals
+            if (isTemp) xWell[output_index++] = smry.get_well_var(well.name(), "WTPTHEA", 0);
+            for (const auto& tracer : tracers) {
+                if (tracer.phase == Opm::Phase::WATER) {
+                    xWell[output_index++] = smry.get_well_var(well.name(), fmt::format("WTPT{}", tracer.name), 0);
+                } else {
+                    xWell[output_index++] = smry.get_well_var(well.name(), fmt::format("WTPTF{}", tracer.name), 0);
+                    xWell[output_index++] = smry.get_well_var(well.name(), fmt::format("WTPTS{}", tracer.name), 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 + 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;
+            // Injection totals
+            if (isTemp) xWell[output_index++] = smry.get_well_var(well.name(), "WTITHEA", 0);
+            for (const auto& tracer : tracers) {
+                if (tracer.phase == Opm::Phase::WATER) {
+                    xWell[output_index++] = smry.get_well_var(well.name(), fmt::format("WTIT{}", tracer.name), 0);
+                } else {
+                    xWell[output_index++] = smry.get_well_var(well.name(), fmt::format("WTITF{}", tracer.name), 0);
+                    xWell[output_index++] = smry.get_well_var(well.name(), fmt::format("WTITS{}", tracer.name), 0);
                 }
             }
 
-            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;
+            // inj and prod conc
+            const bool injector = well.isInjector();
+            for (std::size_t n=0; n < 2; ++n) {
+                if (isTemp) {
+                        xWell[output_index++] = injector ? smry.get_well_var(well.name(), "WTICHEA", 0) : smry.get_well_var(well.name(), "WTPCHEA", 0);
+                }
+                for (const auto& tracer : tracers) {
+                    if (tracer.phase == Opm::Phase::WATER) {
+                        xWell[output_index++] = injector ? smry.get_well_var(well.name(), fmt::format("WTIC{}", tracer.name), 0) : smry.get_well_var(well.name(), fmt::format("WTPC{}", tracer.name), 0);
+                    } else {
+                        xWell[output_index++] = injector ? smry.get_well_var(well.name(), fmt::format("WTICF{}", tracer.name), 0) : smry.get_well_var(well.name(), fmt::format("WTPCF{}", tracer.name), 0);
+                        xWell[output_index++] = injector ? smry.get_well_var(well.name(), fmt::format("WTICS{}", tracer.name), 0) : smry.get_well_var(well.name(), fmt::format("WTPCS{}", tracer.name), 0);
+                    }
                 }
             }
 
-            std::size_t output_index = Ix::TracerOffset + 5*tracer_dims.water_tracers();
-            xWell[output_index] = 0;
-            xWell[output_index + 1] = 0;
+            xWell[output_index++] = 0;
+            xWell[output_index++] = 0;
         }
 
         template <class XWellArray>
         void dynamicContrib(const ::Opm::Well&         well,
                             const Opm::TracerConfig&   tracers,
                             const Opm::Tracers&        tracer_dims,
                             const ::Opm::SummaryState& smry,
+                            const bool                 isTemp,
                             XWellArray&                xWell)
         {
             if (well.isProducer()) {
@@ -1433,7 +1462,7 @@ namespace {
                 }
             }
             assignCumulatives(well.name(), smry, xWell);
-            assignTracerData(tracers, tracer_dims, smry, well, xWell);
+            assignTracerData(tracers, tracer_dims, smry, well, xWell, isTemp);
         }
     } // XWell
 
@@ -1598,6 +1627,6 @@ captureDynamicWellData(const Opm::Schedule&       sched,
     {
         auto xwell = this->xWell_[wellID];
 
-        XWell::dynamicContrib(well, tracers, sched.runspec().tracers(), smry, xwell);
+        XWell::dynamicContrib(well, tracers, sched.runspec().tracers(), smry, sched.runspec().temp(), xwell);
     });
 }

opm/output/eclipse/CreateInteHead.cpp

@@ -22,6 +22,7 @@
 
 #include <opm/output/eclipse/InteHEAD.hpp>
 #include <opm/output/eclipse/VectorItems/intehead.hpp>
+#include <opm/output/eclipse/VectorItems/well.hpp>
 
 #include <opm/input/eclipse/EclipseState/Aquifer/AquiferConfig.hpp>
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
@@ -265,7 +266,7 @@ namespace {
     std::array<int, 4>
     getNGRPZ(const int             grpsz,
              const int             ngrp,
-             const int             num_water_tracer,
+             const int             num_tracers,
              const ::Opm::Runspec& rspec)
     {
         const auto& wd = rspec.wellDimensions();
@@ -275,7 +276,7 @@ namespace {
 
         const int nigrpz = 97 + std::max(nwgmax, ngmax);
         const int nsgrpz = 112;
-        const int nxgrpz = 180 + 4*num_water_tracer;
+        const int nxgrpz = 180 + 4*num_tracers + 1;  // +1 since late 2022
         const int nzgrpz = 5;
 
         return {{
@@ -360,7 +361,7 @@ namespace {
 
 
     Opm::RestartIO::InteHEAD::WellSegDims
-    getWellSegDims(const int              num_water_tracer,
+    getWellSegDims(const int              num_tracers,
                    const ::Opm::Runspec&  rspec,
                    const ::Opm::Schedule& sched,
                    const std::size_t      report_step,
@@ -379,7 +380,7 @@ namespace {
             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
+               + 8*num_tracers, // Number of entries per segment in RSEG
             10            // Number of entries per segment in ILBR (2017.2)
         };
     }
@@ -573,8 +574,11 @@ createInteHead(const EclipseState& es,
     const auto& tdim  = es.getTableManager();
     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 = num_water_tracer*5 + 2 * (num_water_tracer > 0);
+    const auto& tracers = es.runspec().tracers();
+    // TEMP is a tracer, oil&gas tracers have both free and solution parts.
+    const auto num_tracers = tracers.water_tracers() + tracers.oil_tracers() + tracers.gas_tracers() + (es.runspec().temp() ? 1 : 0);
+    const auto num_tracer_comps = num_tracers + tracers.oil_tracers() + tracers.gas_tracers();
+    int nxwelz_tracer_shift = num_tracer_comps*5 + 2 * (num_tracers > 0);
 
     const auto ih = InteHEAD{}
         .dimensions         (grid.getNXYZ())
@@ -588,18 +592,18 @@ createInteHead(const EclipseState& es,
              // across a range of reference cases, but are not guaranteed to be
              // universally valid.
         .drsdt(sched, lookup_step)
-        .params_NWELZ       (155 + num_water_tracer, 122 + 2*num_water_tracer, 130 + nxwelz_tracer_shift, 3) // n{isxz}welz: number of data elements per well in {ISXZ}WELL
-        .params_NCON        (25, 41, 58 + 5*num_water_tracer)       // n{isx}conz: number of data elements per completion in ICON
-        .params_GRPZ        (getNGRPZ(nwgmax, ngmax, num_water_tracer, rspec))
+        .params_NWELZ       (155 + num_tracers, 122 + 2 * num_tracer_comps, VectorItems::XWell::index::TracerOffset + nxwelz_tracer_shift, 3) // n{isxz}welz: number of data elements per well in {ISXZ}WELL
+        .params_NCON        (25, 41, 58 + 5*num_tracer_comps)       // n{isx}conz: number of data elements per completion in ICON
+        .params_GRPZ        (getNGRPZ(nwgmax, ngmax, num_tracer_comps, rspec))
         .aquiferDimensions  (inferAquiferDimensions(es, sched[lookup_step]))
         .stepParam          (num_solver_steps, report_step)
         .tuningParam        (getTuningPars(sched[lookup_step].tuning()))
         .liftOptParam       (getLiftOptPar(sched, report_step, lookup_step))
-        .wellSegDimensions  (getWellSegDims(num_water_tracer, rspec, sched, report_step, lookup_step))
+        .wellSegDimensions  (getWellSegDims(num_tracer_comps, rspec, sched, report_step, lookup_step))
         .regionDimensions   (getRegDims(tdim, rdim))
         .ngroups            ({ ngmax })
         .params_NGCTRL      (GroupControl(sched, report_step, lookup_step))
-        .variousParam       (201802, 100)  // Output should be compatible with Eclipse 100, 2017.02 version.
+        .variousParam       (201802, 100, num_tracer_comps)  // Output should be compatible with Eclipse 100, 2022.04 version.
         .udqParam_1         (getUdqParam(rspec, sched, report_step, lookup_step))
         .actionParam        (getActionParam(rspec, acts, report_step))
         .variousUDQ_ACTIONXParam()

opm/output/eclipse/CreateLogiHead.cpp

@@ -54,7 +54,7 @@ createLogiHead(const EclipseState& es)
           tabMgr.getPvcdoTable().empty());
 
     const auto lh = LogiHEAD{}
-        .variousParam(false, false, wsd.maxSegmentedWells(), hystPar.active())
+        .variousParam(false, false, wsd.maxSegmentedWells(), hystPar.active(), rspec.temp())
         .pvtModel(pvt)
         .network(rspec.networkDimensions().maxNONodes())
         ;