From 6d76f6c55720385100f917d58e65d9094232a946 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ove=20S=C3=A6vareid?= <ovsa@norceresearch.no>
Date: Tue, 11 May 2021 13:01:45 +0200
Subject: [PATCH] Tracer summary curves, collect well rates.

---
 ebos/ecltracermodel.hh | 53 +++++++++++++++++++++++++++++++++++++-----
 ebos/eclwriter.hh      |  2 ++
 2 files changed, 49 insertions(+), 6 deletions(-)

diff --git a/ebos/ecltracermodel.hh b/ebos/ecltracermodel.hh
index c77f280bc..eb910eedf 100644
--- a/ebos/ecltracermodel.hh
+++ b/ebos/ecltracermodel.hh
@@ -103,6 +103,28 @@ public:
                      gasPhaseIdx, oilPhaseIdx, waterPhaseIdx);
     }
 
+    /*!
+     * \brief Collect well tracer rates for reporting
+     */
+    void addTracerRatesToWells(Opm::data::WellRates& wellDatas) const
+    {
+
+       if (numTracers()==0)
+          return; // no tracers today
+
+       for (const auto& wTR : wellTracerRate_) {
+          std::string wellName = wTR.first.first;
+          std::string tracerName = tracerNames_[wTR.first.second];
+          double rate = wTR.second;
+          if (!wellDatas.count(wellName)) {
+                Opm::data::Well wellData;
+                wellDatas.emplace(std::make_pair(wellName, wellData));
+          }
+          Opm::data::Well& wellData = wellDatas.at(wellName);
+          wellData.rates.set(Opm::data::Rates::opt::tracer, rate, tracerName);
+       }
+    }
+
     void beginTimeStep()
     {
         if (this->numTracers()==0)
@@ -285,6 +307,8 @@ protected:
         const auto& wells = simulator_.vanguard().schedule().getWells(episodeIdx);
         for (const auto& well : wells) {
 
+            wellTracerRate_[std::make_pair(well.name(),tracerIdx)] = 0.0;
+
             if (well.getStatus() == Well::Status::SHUT)
                 continue;
 
@@ -299,17 +323,34 @@ protected:
                 cartesianCoordinate[1] = connection.getJ();
                 cartesianCoordinate[2] = connection.getK();
                 const size_t cartIdx = simulator_.vanguard().cartesianIndex(cartesianCoordinate);
-                const int I = this->cartToGlobal_[cartIdx];
-                Scalar rate = simulator_.problem().wellModel().well(well.name())->volumetricSurfaceRateForConnection(I, this->tracerPhaseIdx_[tracerIdx]);
-                if (rate > 0)
-                    this->tracerResidual_[I][0] -= rate*wtracer;
-                else if (rate < 0)
-                    this->tracerResidual_[I][0] -= rate*this->tracerConcentration_[tracerIdx][I];
+                const int I = cartToGlobal_[cartIdx];
+                Scalar rate = simulator_.problem().wellModel().well(well.name())->volumetricSurfaceRateForConnection(I, tracerPhaseIdx_[tracerIdx]);
+                if (rate > 0) {
+                    tracerResidual_[I][0] -= rate*wtracer;
+                    wellTracerRate_.at(std::make_pair(well.name(),tracerIdx)) += rate*wtracer;
+                }
+                else if (rate < 0) {
+                    tracerResidual_[I][0] -= rate*tracerConcentration_[tracerIdx][I];
+                    wellTracerRate_.at(std::make_pair(well.name(),tracerIdx)) += rate*tracerConcentration_[tracerIdx][I];
+                }
             }
         }
     }
 
     Simulator& simulator_;
+
+    std::vector<std::string> tracerNames_;
+    std::vector<int> tracerPhaseIdx_;
+    std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> tracerConcentration_;
+    std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> tracerConcentrationInitial_;
+    TracerMatrix *tracerMatrix_;
+    TracerVector tracerResidual_;
+    std::vector<int> cartToGlobal_;
+    std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> storageOfTimeIndex1_;
+
+    // <wellName, tracerIdx> -> wellRate
+    std::map<std::pair<std::string, int>, double> wellTracerRate_;
+
 };
 
 } // namespace Opm
diff --git a/ebos/eclwriter.hh b/ebos/eclwriter.hh
index 7d370307a..0ce6aad99 100644
--- a/ebos/eclwriter.hh
+++ b/ebos/eclwriter.hh
@@ -170,6 +170,8 @@ public:
             simulator_.vanguard().externalSetupTime();
 
         const auto localWellData            = simulator_.problem().wellModel().wellData();
+        simulator_.problem().tracerModel().addTracerRatesToWells(const_cast<Opm::data::WellRates&>(localWellData));
+
         const auto localGroupAndNetworkData = simulator_.problem().wellModel()
             .groupAndNetworkData(reportStepNum);