Batchwise tracer calculations.

ebos/eclgenerictracermodel.cc

@@ -243,6 +243,43 @@ linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
     return result.converged;
 }
 
+template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
+bool  EclGenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
+linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::vector<TracerVector>& b)
+{
+#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
+    Dune::FMatrixPrecision<Scalar>::set_singular_limit(1.e-30);
+    Dune::FMatrixPrecision<Scalar>::set_absolute_limit(1.e-30);
+#endif
+    Scalar tolerance = 1e-2;
+    int maxIter = 100;
+
+    int verbosity = 0;
+    using TracerSolver = Dune::BiCGSTABSolver<TracerVector>;
+    using TracerOperator = Dune::MatrixAdapter<TracerMatrix,TracerVector,TracerVector>;
+    using TracerScalarProduct = Dune::SeqScalarProduct<TracerVector>;
+    using TracerPreconditioner = Dune::SeqILU< TracerMatrix,TracerVector,TracerVector>;
+
+    TracerOperator tracerOperator(M);
+    TracerScalarProduct tracerScalarProduct;
+    TracerPreconditioner tracerPreconditioner(M, 0, 1); // results in ILU0
+
+    TracerSolver solver (tracerOperator, tracerScalarProduct,
+                         tracerPreconditioner, tolerance, maxIter,
+                         verbosity);
+
+    bool converged = true;
+    for (int nrhs =0; nrhs < b.size(); ++nrhs) {
+        x[nrhs] = 0.0;
+        Dune::InverseOperatorResult result;
+        solver.apply(x[nrhs], b[nrhs], result);
+        converged = (converged && result.converged);
+    }
+
+    // return the result of the solver
+    return converged;
+}
+
 #if HAVE_DUNE_FEM
 template class EclGenericTracerModel<Dune::CpGrid,
                                      Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,

ebos/eclgenerictracermodel.hh

@@ -50,6 +50,8 @@ class EclGenericTracerModel {
 public:
     using TracerMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 1, 1>>;
     using TracerVector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
+    using DualTracerMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 2, 2>>;
+    using DualTracerVector = Dune::BlockVector<Dune::FieldVector<Scalar,2>>;
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
 
     /*!
@@ -68,6 +70,12 @@ public:
      */
     Scalar tracerConcentration(int tracerIdx, int globalDofIdx) const;
 
+    /*!
+    * \brief Return well tracer rates
+    */
+    const std::map<std::pair<std::string, std::string>, double>&
+    getWellTracerRates() const {return wellTracerRate_;}
+
 protected:
     EclGenericTracerModel(const GridView& gridView,
                           const EclipseState& eclState,
@@ -85,6 +93,8 @@ protected:
 
     bool linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b);
 
+    bool linearSolveBatchwise_(const TracerMatrix& M, std::vector<TracerVector>& x, std::vector<TracerVector>& b);
+
     const GridView& gridView_;
     const EclipseState& eclState_;
     const CartesianIndexMapper& cartMapper_;
@@ -98,6 +108,10 @@ protected:
     TracerVector tracerResidual_;
     std::vector<int> cartToGlobal_;
     std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> storageOfTimeIndex1_;
+
+    // <wellName, tracerIdx> -> wellRate
+    std::map<std::pair<std::string, std::string>, double> wellTracerRate_;
+
 };
 
 } // namespace Opm

ebos/ecltracermodel.hh

@@ -77,6 +77,9 @@ class EclTracerModel : public EclGenericTracerModel<GetPropType<TypeTag, Propert
 
     using TracerEvaluation = DenseAd::Evaluation<Scalar,1>;
 
+    using TracerVector = typename BaseType::TracerVector;
+    using DualTracerVector = typename BaseType::DualTracerVector;
+
     enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
     enum { numPhases = FluidSystem::numPhases };
     enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
@@ -90,6 +93,11 @@ public:
                    simulator.vanguard().cartesianIndexMapper(),
                    simulator.model().dofMapper())
         , simulator_(simulator)
+        , wat_(TracerBatch<TracerVector>(waterPhaseIdx))
+        , oil_(TracerBatch<TracerVector>(oilPhaseIdx))
+        , gas_(TracerBatch<TracerVector>(gasPhaseIdx))
+        , oilS_(TracerBatch<DualTracerVector>(oilPhaseIdx))
+        , gasS_(TracerBatch<DualTracerVector>(gasPhaseIdx))
     { }
 
 
@@ -101,28 +109,8 @@ public:
         bool enabled = EWOMS_GET_PARAM(TypeTag, bool, EnableTracerModel);
         this->doInit(enabled, simulator_.model().numGridDof(),
                      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);
-       }
+        prepareTracerBatches();
     }
 
     void beginTimeStep()
@@ -130,21 +118,9 @@ public:
         if (this->numTracers()==0)
             return;
 
-        this->tracerConcentrationInitial_ = this->tracerConcentration_;
-
-        // compute storageCache
-        ElementContext elemCtx(simulator_);
-        auto elemIt = simulator_.gridView().template begin</*codim=*/0>();
-        auto elemEndIt = simulator_.gridView().template end</*codim=*/0>();
-        for (; elemIt != elemEndIt; ++ elemIt) {
-            elemCtx.updateAll(*elemIt);
-            int globalDofIdx = elemCtx.globalSpaceIndex(0, 0);
-            for (int tracerIdx = 0; tracerIdx < this->numTracers(); ++ tracerIdx){
-                Scalar storageOfTimeIndex1;
-                computeStorage_(storageOfTimeIndex1, elemCtx, 0, /*timIdx=*/0, tracerIdx);
-                this->storageOfTimeIndex1_[tracerIdx][globalDofIdx] = storageOfTimeIndex1;
-            }
-        }
+        updateStorageCache(wat_);
+        updateStorageCache(oil_);
+        updateStorageCache(gas_);
     }
 
     /*!
@@ -155,19 +131,9 @@ public:
         if (this->numTracers()==0)
             return;
 
-        for (int tracerIdx = 0; tracerIdx < this->numTracers(); ++ tracerIdx){
-
-            typename BaseType::TracerVector dx(this->tracerResidual_.size());
-            // Newton step (currently the system is linear, converge in one iteration)
-            for (int iter = 0; iter < 5; ++ iter){
-                linearize_(tracerIdx);
-                this->linearSolve_(*this->tracerMatrix_, dx, this->tracerResidual_);
-                this->tracerConcentration_[tracerIdx] -= dx;
-
-                if (dx.two_norm()<1e-2)
-                    break;
-            }
-        }
+        advanceTracerFields(wat_);
+        advanceTracerFields(oil_);
+        advanceTracerFields(gas_);
     }
 
     /*!
@@ -189,40 +155,58 @@ public:
     { /* not implemented */ }
 
 protected:
-    // evaluate storage term for all tracers in a single cell
-    template <class LhsEval>
-    void computeStorage_(LhsEval& tracerStorage,
-                         const ElementContext& elemCtx,
-                         unsigned scvIdx,
-                         unsigned timeIdx,
-                         const int tracerIdx)
-    {
-        int globalDofIdx = elemCtx.globalSpaceIndex(scvIdx, timeIdx);
 
+    // evaluate water storage volume(s) in a single cell
+    template <class LhsEval>
+    void computeVolume_(LhsEval& freeVolume,
+                        LhsEval& dissolvedVolume,
+                        const int tracerPhaseIdx,
+                        const ElementContext& elemCtx,
+                        unsigned scvIdx,
+                        unsigned timeIdx)
+    {
         const auto& intQuants = elemCtx.intensiveQuantities(scvIdx, timeIdx);
         const auto& fs = intQuants.fluidState();
         Scalar phaseVolume =
-            decay<Scalar>(fs.saturation(this->tracerPhaseIdx_[tracerIdx]))
-            *decay<Scalar>(fs.invB(this->tracerPhaseIdx_[tracerIdx]))
+            decay<Scalar>(fs.saturation(tracerPhaseIdx))
+            *decay<Scalar>(fs.invB(tracerPhaseIdx))
             *decay<Scalar>(intQuants.porosity());
 
         // avoid singular matrix if no water is present.
         phaseVolume = max(phaseVolume, 1e-10);
 
         if (std::is_same<LhsEval, Scalar>::value)
-            tracerStorage = phaseVolume * this->tracerConcentrationInitial_[tracerIdx][globalDofIdx];
+            freeVolume = phaseVolume;
         else
-            tracerStorage =
-                    phaseVolume
-                    * variable<LhsEval>(this->tracerConcentration_[tracerIdx][globalDofIdx][0], 0);
+            freeVolume = phaseVolume * variable<LhsEval>(1.0, 0);
+
+        if (FluidSystem::enableDissolvedGas() && tracerPhaseIdx == gasPhaseIdx) {
+            dissolvedVolume =
+                decay<Scalar>(fs.saturation(oilPhaseIdx))
+                *decay<Scalar>(fs.invB(oilPhaseIdx))
+                *decay<Scalar>(fs.Rs())
+                *decay<Scalar>(intQuants.porosity());
+        }
+
+        if (FluidSystem::enableVaporizedOil() && tracerPhaseIdx == oilPhaseIdx) {
+            dissolvedVolume =
+                decay<Scalar>(fs.saturation(gasPhaseIdx))
+                *decay<Scalar>(fs.invB(gasPhaseIdx))
+                *decay<Scalar>(fs.Rv())
+                *decay<Scalar>(intQuants.porosity());
+        }
+
     }
 
-    // evaluate the tracerflux over one face
-    void computeFlux_(TracerEvaluation & tracerFlux,
+    // evaluate the flux(es) over one face
+    void computeFlux_(TracerEvaluation & freeFlux,
+                      TracerEvaluation & dissolvedFlux,
+                      bool & isUpFree,
+                      bool & isUpDissolved,
+                      const int tracerPhaseIdx,
                       const ElementContext& elemCtx,
                       unsigned scvfIdx,
-                      unsigned timeIdx,
-                      const int tracerIdx)
+                      unsigned timeIdx)
 
     {
         const auto& stencil = elemCtx.stencil(timeIdx);
@@ -231,33 +215,66 @@ protected:
         const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
         unsigned inIdx = extQuants.interiorIndex();
 
-        const int tracerPhaseIdx = this->tracerPhaseIdx_[tracerIdx];
-
         unsigned upIdx = extQuants.upstreamIndex(tracerPhaseIdx);
-        int globalUpIdx = elemCtx.globalSpaceIndex(upIdx, timeIdx);
 
         const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
         const auto& fs = intQuants.fluidState();
 
         Scalar A = scvf.area();
         Scalar v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx));
-        Scalar b = decay<Scalar>(fs.invB(this->tracerPhaseIdx_[tracerIdx]));
-        Scalar c = this->tracerConcentration_[tracerIdx][globalUpIdx];
+        Scalar b = decay<Scalar>(fs.invB(tracerPhaseIdx));
 
-        if (inIdx == upIdx)
-            tracerFlux = A*v*b*variable<TracerEvaluation>(c, 0);
-        else
-            tracerFlux = A*v*b*c;
+        if (inIdx == upIdx) {
+            freeFlux = A*v*b*variable<TracerEvaluation>(1.0, 0);
+            isUpFree = true;
+        }
+        else {
+            freeFlux = A*v*b*1.0;
+            isUpFree = false;
+        }
+
+        if (FluidSystem::enableDissolvedGas() && tracerPhaseIdx == gasPhaseIdx) {
+            v = decay<Scalar>(extQuants.volumeFlux(oilPhaseIdx));
+            b = decay<Scalar>(fs.invB(oilPhaseIdx))*decay<Scalar>(fs.Rs());
+            upIdx = extQuants.upstreamIndex(oilPhaseIdx);
+            if (inIdx == upIdx) {
+                dissolvedFlux += A*v*b*variable<TracerEvaluation>(1.0, 0);
+                isUpDissolved = true;
+            }
+            else {
+                dissolvedFlux += A*v*b*1.0;
+                isUpDissolved = false;
+            }
+        }
+
+        if (FluidSystem::enableVaporizedOil() && tracerPhaseIdx == oilPhaseIdx) {
+            v = decay<Scalar>(extQuants.volumeFlux(gasPhaseIdx));
+            b = decay<Scalar>(fs.invB(gasPhaseIdx))*decay<Scalar>(fs.Rv());
+            upIdx = extQuants.upstreamIndex(gasPhaseIdx);
+            if (inIdx == upIdx) {
+                dissolvedFlux += A*v*b*variable<TracerEvaluation>(1.0, 0);
+                isUpDissolved = true;
+            }
+            else {
+                dissolvedFlux += A*v*b*1.0;
+                isUpDissolved = false;
+            }
+        }
 
     }
 
-    void linearize_(int tracerIdx)
+    template <class TrRe>
+    void linearizeTracer_(TrRe & tr)
     {
+        if (tr.numTracer() == 0)
+            return;
+
         (*this->tracerMatrix_) = 0.0;
-        this->tracerResidual_ = 0.0;
+        for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx)
+            tr.residual_[tIdx] = 0.0;
 
         size_t numGridDof =  simulator_.model().numGridDof();
-        std::vector<double> volumes(numGridDof, 0.0);
+        //std::vector<double> volumes(numGridDof, 0.0);
         ElementContext elemCtx(simulator_);
         auto elemIt = simulator_.gridView().template begin</*codim=*/0>();
         auto elemEndIt = simulator_.gridView().template end</*codim=*/0>();
@@ -275,44 +292,69 @@ protected:
             Scalar dt = elemCtx.simulator().timeStepSize();
 
             size_t I = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timIdx=*/0);
-            volumes[I] = scvVolume;
-            TracerEvaluation localStorage;
-            TracerEvaluation storageOfTimeIndex0;
-            Scalar storageOfTimeIndex1;
-            computeStorage_(storageOfTimeIndex0, elemCtx, 0, /*timIdx=*/0, tracerIdx);
-            if (elemCtx.enableStorageCache())
-                storageOfTimeIndex1 = this->storageOfTimeIndex1_[tracerIdx][I];
-            else
-                computeStorage_(storageOfTimeIndex1, elemCtx, 0, /*timIdx=*/1, tracerIdx);
+            size_t I1 = elemCtx.globalSpaceIndex(/*dofIdx=*/ 0, /*timIdx=*/1);
+            //volumes[I] = scvVolume;
+
+            Scalar storageOfTimeIndex1[tr.numTracer()];
+            if (elemCtx.enableStorageCache()) {
+                for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+                    storageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I];
+                }
+            }
+            else {
+                Scalar fVolume1, sVolume1;
+                computeVolume_(fVolume1, sVolume1, tr.phaseIdx_, elemCtx, 0, /*timIdx=*/1);
+                for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+                    storageOfTimeIndex1[tIdx] = fVolume1*tr.concentrationInitial_[tIdx][I1];
+                }
+            }
+
+            TracerEvaluation fVolume, sVolume;
+            computeVolume_(fVolume, sVolume, tr.phaseIdx_, elemCtx, 0, /*timIdx=*/0);
+            for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+                Scalar storageOfTimeIndex0 = fVolume.value()*tr.concentration_[tIdx][I];
+                Scalar localStorage = (storageOfTimeIndex0 - storageOfTimeIndex1[tIdx]) * scvVolume/dt;
+                tr.residual_[tIdx][I][0] += localStorage; //residual + flux
+            }
+            (*this->tracerMatrix_)[I][I][0][0] += fVolume.derivative(0) * scvVolume/dt;
 
-            localStorage = (storageOfTimeIndex0 - storageOfTimeIndex1) * scvVolume/dt;
-            this->tracerResidual_[I][0] += localStorage.value(); //residual + flux
-            (*this->tracerMatrix_)[I][I][0][0] = localStorage.derivative(0);
             size_t numInteriorFaces = elemCtx.numInteriorFaces(/*timIdx=*/0);
             for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
-                TracerEvaluation flux;
+                TracerEvaluation flux, sFlux;
                 const auto& face = elemCtx.stencil(0).interiorFace(scvfIdx);
                 unsigned j = face.exteriorIndex();
                 unsigned J = elemCtx.globalSpaceIndex(/*dofIdx=*/ j, /*timIdx=*/0);
-                computeFlux_(flux, elemCtx, scvfIdx, 0, tracerIdx);
-                this->tracerResidual_[I][0] += flux.value(); //residual + flux
-                (*this->tracerMatrix_)[J][I][0][0] = -flux.derivative(0);
-                (*this->tracerMatrix_)[I][J][0][0] = flux.derivative(0);
+                bool isUpF, isUpS;
+                computeFlux_(flux, sFlux, isUpF, isUpS, tr.phaseIdx_, elemCtx, scvfIdx, 0);
+                int globalUpIdx = isUpF ? I : J;
+                for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+                    tr.residual_[tIdx][I][0] += flux.value()*tr.concentration_[tIdx][globalUpIdx]; //residual + flux
+                }
+                if (isUpF) {
+                    (*this->tracerMatrix_)[J][I][0][0] = -flux.derivative(0);
+                    (*this->tracerMatrix_)[I][I][0][0] += flux.derivative(0);
+                }
             }
 
         }
 
-        // Wells
+        // Wells  terms
         const int episodeIdx = simulator_.episodeIndex();
         const auto& wells = simulator_.vanguard().schedule().getWells(episodeIdx);
         for (const auto& well : wells) {
 
-            wellTracerRate_[std::make_pair(well.name(),tracerIdx)] = 0.0;
+            for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+                this->wellTracerRate_[std::make_pair(well.name(),this->tracerNames_[tr.idx_[tIdx]])] = 0.0;
+            }
 
             if (well.getStatus() == Well::Status::SHUT)
                 continue;
 
-            const double wtracer = well.getTracerProperties().getConcentration(this->tracerNames_[tracerIdx]);
+            double wtracer[tr.numTracer()];
+            for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+                wtracer[tIdx] = well.getTracerProperties().getConcentration(this->tracerNames_[tr.idx_[tIdx]]);
+            }
+
             std::array<int, 3> cartesianCoordinate;
             for (auto& connection : well.getConnections()) {
 
@@ -323,34 +365,162 @@ protected:
                 cartesianCoordinate[1] = connection.getJ();
                 cartesianCoordinate[2] = connection.getK();
                 const size_t cartIdx = simulator_.vanguard().cartesianIndex(cartesianCoordinate);
-                const int I = cartToGlobal_[cartIdx];
-                Scalar rate = simulator_.problem().wellModel().well(well.name())->volumetricSurfaceRateForConnection(I, tracerPhaseIdx_[tracerIdx]);
+                const int I = this->cartToGlobal_[cartIdx];
+                Scalar rate = simulator_.problem().wellModel().well(well.name())->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
                 if (rate > 0) {
-                    tracerResidual_[I][0] -= rate*wtracer;
-                    wellTracerRate_.at(std::make_pair(well.name(),tracerIdx)) += rate*wtracer;
+                    for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+                        tr.residual_[tIdx][I][0] -= rate*wtracer[tIdx];
+                        // Store _injector_ tracer rate for reporting
+                        this->wellTracerRate_.at(std::make_pair(well.name(),this->tracerNames_[tr.idx_[tIdx]])) += rate*wtracer[tIdx];
+                    }
                 }
                 else if (rate < 0) {
-                    tracerResidual_[I][0] -= rate*tracerConcentration_[tracerIdx][I];
-                    wellTracerRate_.at(std::make_pair(well.name(),tracerIdx)) += rate*tracerConcentration_[tracerIdx][I];
+                    for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+                        tr.residual_[tIdx][I][0] -= rate*tr.concentration_[tIdx][I];
+                    }
+                    (*this->tracerMatrix_)[I][I][0][0] -= rate*variable<TracerEvaluation>(1.0, 0).derivative(0);
                 }
             }
         }
     }
 
+    template <class TrRe>
+    void updateStorageCache(TrRe & tr)
+    {
+        if (tr.numTracer() == 0)
+            return;
+
+        tr.concentrationInitial_ = tr.concentration_;
+
+        ElementContext elemCtx(simulator_);
+        auto elemIt = simulator_.gridView().template begin</*codim=*/0>();
+        auto elemEndIt = simulator_.gridView().template end</*codim=*/0>();
+        for (; elemIt != elemEndIt; ++ elemIt) {
+            elemCtx.updateAll(*elemIt);
+            int globalDofIdx = elemCtx.globalSpaceIndex(0, /*timIdx=*/0);
+            Scalar fVolume, sVolume;
+            computeVolume_(fVolume, sVolume, tr.phaseIdx_, elemCtx, 0, /*timIdx=*/0);
+            for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+                tr.storageOfTimeIndex1_[tIdx][globalDofIdx] = fVolume*tr.concentrationInitial_[tIdx][globalDofIdx];
+            }
+        }
+    }
+
+    template <class TrRe>
+    void advanceTracerFields(TrRe & tr)
+    {
+        if (tr.numTracer() == 0)
+            return;
+
+        std::vector<TracerVector> dx(tr.concentration_);
+        for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx)
+            dx[tIdx] = 0.0;
+
+        linearizeTracer_(tr);
+
+        bool converged = this->linearSolveBatchwise_(*this->tracerMatrix_, dx, tr.residual_);
+        if (!converged)
+            std::cout << "### Tracer model: Warning, linear solver did not converge. ###" << std::endl;
+
+        for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+            tr.concentration_[tIdx] -= dx[tIdx];
+            this->tracerConcentration_[tr.idx_[tIdx]] = tr.concentration_[tIdx];
+        }
+
+        // Store _producer_ tracer rate for reporting
+        const int episodeIdx = simulator_.episodeIndex();
+        const auto& wells = simulator_.vanguard().schedule().getWells(episodeIdx);
+        for (const auto& well : wells) {
+
+            if (well.getStatus() == Well::Status::SHUT)
+                continue;
+
+            std::array<int, 3> cartesianCoordinate;
+            for (auto& connection : well.getConnections()) {
+
+                if (connection.state() == Connection::State::SHUT)
+                    continue;
+
+                cartesianCoordinate[0] = connection.getI();
+                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, tr.phaseIdx_);
+                if (rate < 0) { //Injection rates already reported during assembly
+                    for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
+                        this->wellTracerRate_.at(std::make_pair(well.name(),this->tracerNames_[tr.idx_[tIdx]])) += rate*tr.concentration_[tIdx][I];
+                    }
+                }
+            }
+        }
+    }
+
+    void prepareTracerBatches()
+    {
+        for (int tracerIdx=0; tracerIdx<this->tracerPhaseIdx_.size(); ++tracerIdx) {
+            if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::waterPhaseIdx) {
+                if (! FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)){
+                    throw std::runtime_error("Water tracer specified for non-water fluid system:" + this->tracerNames_[tracerIdx]);
+                }
+
+                wat_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
+            }
+            else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::oilPhaseIdx) {
+                if (! FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
+                    throw std::runtime_error("Oil tracer specified for non-oil fluid system:" + this->tracerNames_[tracerIdx]);
+                }
+                if (FluidSystem::enableVaporizedOil()) {
+                    throw std::runtime_error("Oil tracer in combination with kw VAPOIL is not supported: " + this->tracerNames_[tracerIdx]);
+                }
+
+                oil_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
+            }
+            else if (this->tracerPhaseIdx_[tracerIdx] == FluidSystem::gasPhaseIdx) {
+                if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
+                    throw std::runtime_error("Gas tracer specified for non-gas fluid system:" + this->tracerNames_[tracerIdx]);
+                }
+                if (FluidSystem::enableDissolvedGas()) {
+                    throw std::runtime_error("Gas tracer in combination with kw DISGAS is not supported: " + this->tracerNames_[tracerIdx]);
+                }
+
+                gas_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
+            }
+        }
+    }
+
     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_;
+    template <typename TV>
+    struct TracerBatch {
+      std::vector<int> idx_;
+      const int phaseIdx_;
+      std::vector<TV> concentrationInitial_;
+      std::vector<TV> concentration_;
+      std::vector<TV> storageOfTimeIndex1_;
+      std::vector<TV> residual_;
 
-    // <wellName, tracerIdx> -> wellRate
-    std::map<std::pair<std::string, int>, double> wellTracerRate_;
+      TracerBatch(int phaseIdx) : phaseIdx_(phaseIdx) {}
 
+      const int numTracer() const {return idx_.size(); }
+
+      void addTracer(const int idx, const TV & concentration)
+      {
+          int numGridDof = concentration.size();
+          idx_.emplace_back(idx);
+          concentrationInitial_.emplace_back(concentration);
+          concentration_.emplace_back(concentration);
+          storageOfTimeIndex1_.emplace_back(numGridDof);
+          residual_.emplace_back(numGridDof);
+      }
+    };
+
+    TracerBatch<TracerVector> wat_;
+    TracerBatch<TracerVector> oil_;
+    TracerBatch<TracerVector> gas_;
+
+    TracerBatch<DualTracerVector> oilS_;
+    TracerBatch<DualTracerVector> gasS_;
 };
 
 } // namespace Opm

ebos/eclwriter.hh

@@ -170,8 +170,6 @@ 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);
 

opm/simulators/wells/BlackoilWellModel.hpp

@@ -63,6 +63,7 @@
 #include <opm/simulators/wells/WellProdIndexCalculator.hpp>
 #include <opm/simulators/wells/ParallelWellInfo.hpp>
 #include <opm/simulators/timestepping/gatherConvergenceReport.hpp>
+#include <ebos/eclgenerictracermodel.hh>
 #include <dune/common/fmatrix.hh>
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/matrixmatrix.hh>
@@ -218,6 +219,8 @@ namespace Opm {
                                                           return this->wasDynamicallyShutThisTimeStep(well_ndex);
                                                       });
 
+                this->assignWellTracerRates(wsrpt, ebosSimulator_.problem().tracerModel().getWellTracerRates());
+
                 this->assignWellGuideRates(wsrpt);
                 this->assignShutConnections(wsrpt, this->reportStepIndex());
 
@@ -377,7 +380,10 @@ namespace Opm {
                            const int pvtreg,
                            std::vector<double>& resv_coeff) override;
 
-             void computeWellTemperature();
+            void computeWellTemperature();
+
+            void assignWellTracerRates(data::Wells& wsrpt,
+                                       const std::map<std::pair<std::string, std::string>, double>& wellTracerRates) const;
 
         private:
             BlackoilWellModel(Simulator& ebosSimulator, const PhaseUsage& pu);

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -1561,4 +1561,33 @@ namespace Opm {
             this->wellState().update_temperature(wellID, weighted_temperature/total_weight);
         }
     }
+
+
+
+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    assignWellTracerRates(data::Wells& wsrpt,
+                          const std::map<std::pair<std::string, std::string>, double>& wellTracerRates) const
+    {
+
+        if (wellTracerRates.empty())
+            return; // no tracers
+
+        for (const auto& wTR : wellTracerRates) {
+            std::string wellName = wTR.first.first;
+            auto xwPos = wsrpt.find(wellName);
+            if (xwPos == wsrpt.end()) { // No well results.
+                continue;
+            }
+            std::string tracerName = wTR.first.second;
+            double rate = wTR.second;
+            xwPos->second.rates.set(data::Rates::opt::tracer, rate, tracerName);
+        }
+    }
+
+
+
+
+
 } // namespace Opm