Extend tracer model to solution tracers.

Solve an extended linear system with free and solution tracers with mass transfer coupling term.

opm/simulators/flow/EclWriter.hpp

@@ -500,8 +500,10 @@ public:
 
         {
             const auto& tracers = simulator_.vanguard().eclState().tracer();
-            for (const auto& tracer : tracers)
+            for (const auto& tracer : tracers) {
                 solutionKeys.emplace_back(tracer.fname(), UnitSystem::measure::identity, true);
+                solutionKeys.emplace_back(tracer.sname(), UnitSystem::measure::identity, true);
+            }
         }
 
         // The episodeIndex is rewined one back before beginRestart is called
@@ -521,17 +523,20 @@ public:
             auto restartValues = loadParallelRestart(this->eclIO_.get(), actionState, summaryState, solutionKeys, extraKeys,
                                                      gridView.grid().comm());
             for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
-                unsigned globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
+                unsigned globalIdx = this->collectToIORank_.localIdxToGlobalIdx(elemIdx);
-                outputModule_->setRestart(restartValues.solution, elemIdx, globalIdx);
+                eclOutputModule_->setRestart(restartValues.solution, elemIdx, globalIdx);
             }
 
             auto& tracer_model = simulator_.problem().tracerModel();
             for (int tracer_index = 0; tracer_index < tracer_model.numTracers(); tracer_index++) {
-                const auto& tracer_name = tracer_model.fname(tracer_index);
+                const auto& free_tracer_name = tracer_model.fname(tracer_index);
-                const auto& tracer_solution = restartValues.solution.template data<double>(tracer_name);
+                const auto& free_tracer_solution = restartValues.solution.template data<double>(free_tracer_name);
+                const auto& sol_tracer_name = tracer_model.sname(tracer_index);
+                const auto& sol_tracer_solution = restartValues.solution.template data<double>(sol_tracer_name);
                 for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
-                    unsigned globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
+                    unsigned globalIdx = this->collectToIORank_.localIdxToGlobalIdx(elemIdx);
-                    tracer_model.setTracerConcentration(tracer_index, globalIdx, tracer_solution[globalIdx]);
+                    tracer_model.setFreeTracerConcentration(tracer_index, globalIdx, free_tracer_solution[globalIdx]);
+                    tracer_model.setSolTracerConcentration(tracer_index, globalIdx, sol_tracer_solution[globalIdx]);
                 }
             }
 

opm/simulators/flow/GenericOutputBlackoilModule.cpp

@@ -700,21 +700,37 @@ assignToSolution(data::Solution& sol)
     }
 
     // Tracers
-    if (! this->tracerConcentrations_.empty()) {
+    if (! this->freeTracerConcentrations_.empty()) {
         const auto& tracers = this->eclState_.tracer();
         for (auto tracerIdx = 0*tracers.size();
              tracerIdx < tracers.size(); ++tracerIdx)
         {
             sol.insert(tracers[tracerIdx].fname(),
                        UnitSystem::measure::identity,
-                       std::move(tracerConcentrations_[tracerIdx]),
+                       std::move(freeTracerConcentrations_[tracerIdx]),
                        data::TargetType::RESTART_TRACER_SOLUTION);
         }
 
-        // Put tracerConcentrations container into a valid state.  Otherwise
+        // Put freeTracerConcentrations container into a valid state.  Otherwise
         // we'll move from vectors that have already been moved from if we
         // get here and it's not a restart step.
-        this->tracerConcentrations_.clear();
+        this->freeTracerConcentrations_.clear();
+    }
+    if (! this->solTracerConcentrations_.empty()) {
+        const auto& tracers = this->eclState_.tracer();
+        for (auto tracerIdx = 0*tracers.size();
+             tracerIdx < tracers.size(); ++tracerIdx)
+        {
+            sol.insert(tracers[tracerIdx].sname(),
+                       UnitSystem::measure::identity,
+                       std::move(solTracerConcentrations_[tracerIdx]),
+                       data::TargetType::RESTART_TRACER_SOLUTION);
+        }
+
+        // Put solTracerConcentrations container into a valid state.  Otherwise
+        // we'll move from vectors that have already been moved from if we
+        // get here and it's not a restart step.
+        this->solTracerConcentrations_.clear();
     }
 }
 
@@ -1304,10 +1320,15 @@ doAllocBuffers(const unsigned bufferSize,
 
     // tracers
     if (numTracers > 0) {
-        tracerConcentrations_.resize(numTracers);
+        freeTracerConcentrations_.resize(numTracers);
         for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx)
         {
-            tracerConcentrations_[tracerIdx].resize(bufferSize, 0.0);
+            freeTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0);
+        }
+        solTracerConcentrations_.resize(numTracers);
+        for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx)
+        {
+            solTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0);
         }
     }
 

opm/simulators/flow/GenericOutputBlackoilModule.hpp

@@ -504,7 +504,8 @@ protected:
     std::array<ScalarBuffer, numPhases> viscosity_;
     std::array<ScalarBuffer, numPhases> relativePermeability_;
 
-    std::vector<ScalarBuffer> tracerConcentrations_;
+    std::vector<ScalarBuffer> freeTracerConcentrations_;
+    std::vector<ScalarBuffer> solTracerConcentrations_;
 
     std::array<ScalarBuffer, numPhases> residual_;
 

opm/simulators/flow/GenericTracerModel.hpp

@@ -52,8 +52,9 @@ class Well;
 template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar>
 class GenericTracerModel {
 public:
-    using TracerMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 1, 1>>;
+    using TracerVectorSingle = Dune::BlockVector<Dune::FieldVector<Scalar, 1>>;
-    using TracerVector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
+    using TracerMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 2, 2>>;
+    using TracerVector = Dune::BlockVector<Dune::FieldVector<Scalar, 2>>;
     using CartesianIndexMapper = Dune::CartesianIndexMapper<Grid>;
     static constexpr int dimWorld = Grid::dimensionworld;
     /*!
@@ -66,13 +67,16 @@ public:
      */
     const std::string& name(int tracerIdx) const;
     std::string fname(int tracerIdx) const;
+    std::string sname(int tracerIdx) const;
 
 
     /*!
      * \brief Return the tracer concentration for tracer index and global DofIdx
      */
-    Scalar tracerConcentration(int tracerIdx, int globalDofIdx) const;
+    Scalar freeTracerConcentration(int tracerIdx, int globalDofIdx) const;
-    void setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value);
+    Scalar solTracerConcentration(int tracerIdx, int globalDofIdx) const;
+    void setFreeTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value);
+    void setSolTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value);
 
     /*!
     * \brief Return well tracer rates
@@ -117,9 +121,10 @@ protected:
     const DofMapper& dofMapper_;
 
     std::vector<int> tracerPhaseIdx_;
-    std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> tracerConcentration_;
+    std::vector<TracerVector> tracerConcentration_;
     std::unique_ptr<TracerMatrix> tracerMatrix_;
-    std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> storageOfTimeIndex1_;
+    std::vector<TracerVectorSingle> freeTracerConcentration_;
+    std::vector<TracerVectorSingle> solTracerConcentration_;
 
     // <wellName, tracerIdx> -> wellRate
     std::map<std::pair<std::string, std::string>, Scalar> wellTracerRate_;

opm/simulators/flow/GenericTracerModel_impl.hpp

@@ -114,19 +114,36 @@ GenericTracerModel(const GridView& gridView,
 
 template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
 Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
-tracerConcentration(int tracerIdx, int globalDofIdx) const
+freeTracerConcentration(int tracerIdx, int globalDofIdx) const
 {
-    if (tracerConcentration_.empty())
+    if (freeTracerConcentration_.empty())
         return 0.0;
 
-    return tracerConcentration_[tracerIdx][globalDofIdx];
+    return freeTracerConcentration_[tracerIdx][globalDofIdx];
+}
+
+template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
+Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
+solTracerConcentration(int tracerIdx, int globalDofIdx) const
+{
+    if (solTracerConcentration_.empty())
+        return 0.0;
+
+    return solTracerConcentration_[tracerIdx][globalDofIdx];
 }
 
 template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
 void GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
-setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value)
+setFreeTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value)
 {
-    this->tracerConcentration_[tracerIdx][globalDofIdx] = value;
+    this->freeTracerConcentration_[tracerIdx][globalDofIdx] = value;
+}
+
+template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
+void GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
+setSolTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value)
+{
+    this->solTracerConcentration_[tracerIdx][globalDofIdx] = value;
 }
 
 template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
@@ -143,6 +160,13 @@ fname(int tracerIdx) const
     return this->eclState_.tracer()[tracerIdx].fname();
 }
 
+template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
+std::string GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
+sname(int tracerIdx) const
+{
+    return this->eclState_.tracer()[tracerIdx].sname();
+}
+
 template<class Grid,class GridView, class DofMapper, class Stencil, class Scalar>
 Scalar GenericTracerModel<Grid,GridView,DofMapper,Stencil,Scalar>::
 currentConcentration_(const Well& eclWell, const std::string& name) const
@@ -170,7 +194,8 @@ doInit(bool rst, std::size_t numGridDof,
     // retrieve the number of tracers from the deck
     const std::size_t numTracers = tracers.size();
     tracerConcentration_.resize(numTracers);
-    storageOfTimeIndex1_.resize(numTracers);
+    freeTracerConcentration_.resize(numTracers);
+    solTracerConcentration_.resize(numTracers);
 
     // the phase where the tracer is
     tracerPhaseIdx_.resize(numTracers);
@@ -185,34 +210,71 @@ doInit(bool rst, std::size_t numGridDof,
             tracerPhaseIdx_[tracerIdx] = gasPhaseIdx;
 
         tracerConcentration_[tracerIdx].resize(numGridDof);
-        storageOfTimeIndex1_[tracerIdx].resize(numGridDof);
+        tracerConcentration_[tracerIdx].resize(numGridDof);
+        freeTracerConcentration_[tracerIdx].resize(numGridDof);
+        solTracerConcentration_[tracerIdx].resize(numGridDof);
 
         if (rst)
             continue;
 
 
-        //TBLK keyword
+        // TBLKF keyword
         if (tracer.free_concentration.has_value()){
             const auto& free_concentration = tracer.free_concentration.value();
             int tblkDatasize = free_concentration.size();
             if (tblkDatasize < cartMapper_.cartesianSize()){
-                throw std::runtime_error("Wrong size of TBLK for" + tracer.name);
+                throw std::runtime_error("Wrong size of TBLKF for" + tracer.name);
             }
             for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
                 int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
-                tracerConcentration_[tracerIdx][globalDofIdx] = free_concentration[cartDofIdx];
+                tracerConcentration_[tracerIdx][globalDofIdx][0] = free_concentration[cartDofIdx];
+                freeTracerConcentration_[tracerIdx][globalDofIdx] = free_concentration[cartDofIdx];
             }
         }
-        //TVDPF keyword
+        // TVDPF keyword
         else if (tracer.free_tvdp.has_value()) {
             const auto& free_tvdp = tracer.free_tvdp.value();
             for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
-                tracerConcentration_[tracerIdx][globalDofIdx] =
+                tracerConcentration_[tracerIdx][globalDofIdx][0] =
+                    free_tvdp.evaluate("TRACER_CONCENTRATION",
+                                       centroids_(globalDofIdx)[2]);
+                freeTracerConcentration_[tracerIdx][globalDofIdx] =
                     free_tvdp.evaluate("TRACER_CONCENTRATION",
                                        centroids_(globalDofIdx)[2]);
             }
-        } else
+        } 
-            throw std::logic_error(fmt::format("Can not initialize tracer: {}", tracer.name));
+        else {
+            throw std::logic_error(fmt::format("Can not initialize free tracer concentration: {}", tracer.name));
+        }
+
+        // TBLKS keyword
+        if (tracer.solution_concentration.has_value()){
+            const auto& solution_concentration = tracer.solution_concentration.value();
+            int tblkDatasize = solution_concentration.size();
+            if (tblkDatasize < cartMapper_.cartesianSize()){
+                throw std::runtime_error("Wrong size of TBLKS for" + tracer.name);
+            }
+            for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
+                int cartDofIdx = cartMapper_.cartesianIndex(globalDofIdx);
+                tracerConcentration_[tracerIdx][globalDofIdx][1] = solution_concentration[cartDofIdx];
+                solTracerConcentration_[tracerIdx][globalDofIdx] = solution_concentration[cartDofIdx];
+            }
+        }
+        // TVDPS keyword
+        else if (tracer.solution_tvdp.has_value()) {
+            const auto& solution_tvdp = tracer.solution_tvdp.value();
+            for (std::size_t globalDofIdx = 0; globalDofIdx < numGridDof; ++globalDofIdx) {
+                tracerConcentration_[tracerIdx][globalDofIdx][1] =
+                    solution_tvdp.evaluate("TRACER_CONCENTRATION",
+                                           centroids_(globalDofIdx)[2]);
+                solTracerConcentration_[tracerIdx][globalDofIdx] =
+                    solution_tvdp.evaluate("TRACER_CONCENTRATION",
+                                           centroids_(globalDofIdx)[2]);
+            }
+        } 
+        else {
+            throw std::logic_error(fmt::format("Can not initialize solution tracer concentration: {}", tracer.name));
+        }
     }
 
     // allocate matrix for storing the Jacobian of the tracer residual
@@ -237,8 +299,9 @@ doInit(bool rst, std::size_t numGridDof,
     }
 
     // allocate space for the rows of the matrix
-    for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx)
+    for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
         tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
+    }
     tracerMatrix_->endrowsizes();
 
     // fill the rows with indices. each degree of freedom talks to
@@ -247,8 +310,9 @@ doInit(bool rst, std::size_t numGridDof,
     for (unsigned dofIdx = 0; dofIdx < numGridDof; ++ dofIdx) {
         typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
         typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
-        for (; nIt != nEndIt; ++nIt)
+        for (; nIt != nEndIt; ++nIt) {
             tracerMatrix_->addindex(dofIdx, *nIt);
+        }
     }
     tracerMatrix_->endindices();
 }

opm/simulators/flow/OutputBlackoilModule.hpp

@@ -661,14 +661,23 @@ public:
 
             // tracers
             const auto& tracerModel = simulator_.problem().tracerModel();
-            if (! this->tracerConcentrations_.empty()) {
+            if (! this->freeTracerConcentrations_.empty()) {
                 for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) {
-                    if (this->tracerConcentrations_[tracerIdx].empty()) {
+                    if (this->freeTracerConcentrations_[tracerIdx].empty()) {
                         continue;
                     }
-
+                    this->freeTracerConcentrations_[tracerIdx][globalDofIdx] =
-                    this->tracerConcentrations_[tracerIdx][globalDofIdx] =
+                        tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx);
-                        tracerModel.tracerConcentration(tracerIdx, globalDofIdx);
+                }
+            }
+            if (! this->solTracerConcentrations_.empty()) {
+                for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); ++tracerIdx) {
+                    if (this->solTracerConcentrations_[tracerIdx].empty()) {
+                        continue;
+                    }
+                    this->solTracerConcentrations_[tracerIdx][globalDofIdx] =
+                        tracerModel.solTracerConcentration(tracerIdx, globalDofIdx);
+                    
                 }
             }
 

opm/simulators/flow/TracerModel.hpp

@@ -113,7 +113,7 @@ public:
          some phase index stuff. If this is a normal run the initial tracer
          concentrations will be assigned from the TBLK or TVDPF keywords.
 
-      2. [Restart only:] The tracer concenntration are read from the restart
+      2. [Restart only:] The tracer concentration are read from the restart
          file and the concentrations are applied with repeated calls to the
          setTracerConcentration() method. This is currently done in the
          eclwriter::beginRestart() method.
@@ -144,9 +144,6 @@ public:
                 if (! FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)){
                     throw std::runtime_error("Oil tracer specified for non-oil fluid system:" + this->name(tracerIdx));
                 }
-                if (FluidSystem::enableVaporizedOil()) {
-                    throw std::runtime_error("Oil tracer in combination with kw VAPOIL is not supported: " + this->name(tracerIdx));
-                }
 
                 oil_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
             }
@@ -154,12 +151,13 @@ public:
                 if (! FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)){
                     throw std::runtime_error("Gas tracer specified for non-gas fluid system:" + this->name(tracerIdx));
                 }
-                if (FluidSystem::enableDissolvedGas()) {
-                    throw std::runtime_error("Gas tracer in combination with kw DISGAS is not supported: " + this->name(tracerIdx));
-                }
 
                 gas_.addTracer(tracerIdx, this->tracerConcentration_[tracerIdx]);
             }
+
+            // resize free and solution volume storages
+            fVol1_[this->tracerPhaseIdx_[tracerIdx]].resize(this->freeTracerConcentration_[tracerIdx].size());
+            sVol1_[this->tracerPhaseIdx_[tracerIdx]].resize(this->solTracerConcentration_[tracerIdx].size());
         }
 
         // will be valid after we move out of tracerMatrix_
@@ -220,39 +218,63 @@ public:
 
 protected:
 
-    // evaluate water storage volume(s) in a single cell
+    // compute volume associated with free concentration
-    template <class LhsEval>
+    Scalar computeFreeVolume_(const int tracerPhaseIdx,
-    void computeVolume_(LhsEval& freeVolume,
+                              unsigned globalDofIdx,
-                        const int tracerPhaseIdx,
+                              unsigned timeIdx)
-                        const ElementContext& elemCtx,
-                        unsigned scvIdx,
-                        unsigned timeIdx)
     {
-        const auto& intQuants = elemCtx.intensiveQuantities(scvIdx, timeIdx);
+        const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx);
         const auto& fs = intQuants.fluidState();
+
         Scalar phaseVolume =
             decay<Scalar>(fs.saturation(tracerPhaseIdx))
             *decay<Scalar>(fs.invB(tracerPhaseIdx))
             *decay<Scalar>(intQuants.porosity());
 
-        // avoid singular matrix if no water is present.
+        return max(phaseVolume, 1e-10);
-        phaseVolume = max(phaseVolume, 1e-10);
+    }
 
-        if (std::is_same<LhsEval, Scalar>::value)
+    // compute volume associated with solution concentration
-            freeVolume = phaseVolume;
+    Scalar computeSolutionVolume_(const int tracerPhaseIdx,
-        else
+                                  unsigned globalDofIdx,
-            freeVolume = phaseVolume * variable<LhsEval>(1.0, 0);
+                                  unsigned timeIdx)
+    {
+        const auto& intQuants = simulator_.model().intensiveQuantities(globalDofIdx, timeIdx);
+        const auto& fs = intQuants.fluidState();
+
+        Scalar phaseVolume;
+
+        // vaporized oil
+        if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
+            phaseVolume =
+                decay<Scalar>(fs.saturation(FluidSystem::gasPhaseIdx))
+                * decay<Scalar>(fs.invB(FluidSystem::gasPhaseIdx))
+                * decay<Scalar>(fs.Rv())
+                * decay<Scalar>(intQuants.porosity());
+        }
+
+        // dissolved gas
+        else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
+            phaseVolume =
+                decay<Scalar>(fs.saturation(FluidSystem::oilPhaseIdx))
+                * decay<Scalar>(fs.invB(FluidSystem::oilPhaseIdx))
+                * decay<Scalar>(fs.Rs())
+                * decay<Scalar>(intQuants.porosity());
+        }
+        else {
+            phaseVolume = 0.0;
+        }
+
+        return max(phaseVolume, 1e-10);
 
     }
 
-    // evaluate the flux(es) over one face
+    void computeFreeFlux_(TracerEvaluation & freeFlux,
-    void computeFlux_(TracerEvaluation & freeFlux,
+                          bool & isUp,
-                      bool & isUpFree,
+                          const int tracerPhaseIdx,
-                      const int tracerPhaseIdx,
+                          const ElementContext& elemCtx,
-                      const ElementContext& elemCtx,
+                          unsigned scvfIdx,
-                      unsigned scvfIdx,
+                          unsigned timeIdx)
-                      unsigned timeIdx)
-
     {
         const auto& stencil = elemCtx.stencil(timeIdx);
         const auto& scvf = stencil.interiorFace(scvfIdx);
@@ -265,17 +287,72 @@ protected:
         const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
         const auto& fs = intQuants.fluidState();
 
-        Scalar A = scvf.area();
+        Scalar v = 
-        Scalar v = decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx));
+                decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx)) 
-        Scalar b = decay<Scalar>(fs.invB(tracerPhaseIdx));
+                * decay<Scalar>(fs.invB(tracerPhaseIdx));
 
+        Scalar A = scvf.area();
         if (inIdx == upIdx) {
-            freeFlux = A*v*b*variable<TracerEvaluation>(1.0, 0);
+            freeFlux = A*v*variable<TracerEvaluation>(1.0, 0);
-            isUpFree = true;
+            isUp = true;
         }
         else {
-            freeFlux = A*v*b*1.0;
+            freeFlux = A*v;
-            isUpFree = false;
+            isUp = false;
+        }
+    }
+
+    void computeSolFlux_(TracerEvaluation & solFlux,
+                         bool & isUp,
+                         const int tracerPhaseIdx,
+                         const ElementContext& elemCtx,
+                         unsigned scvfIdx,
+                         unsigned timeIdx)
+    {
+        const auto& stencil = elemCtx.stencil(timeIdx);
+        const auto& scvf = stencil.interiorFace(scvfIdx);
+
+        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+        unsigned inIdx = extQuants.interiorIndex();
+
+        Scalar v;
+        unsigned upIdx;
+
+        // vaporized oil
+        if (tracerPhaseIdx == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
+            upIdx = extQuants.upstreamIndex(FluidSystem::gasPhaseIdx);
+
+            const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
+            const auto& fs = intQuants.fluidState();
+            v =
+                decay<Scalar>(fs.invB(FluidSystem::gasPhaseIdx))
+                * decay<Scalar>(extQuants.volumeFlux(FluidSystem::gasPhaseIdx))
+                * decay<Scalar>(fs.Rv());
+        }
+        // dissolved gas
+        else if (tracerPhaseIdx == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
+            upIdx = extQuants.upstreamIndex(FluidSystem::oilPhaseIdx);
+
+            const auto& intQuants = elemCtx.intensiveQuantities(upIdx, timeIdx);
+            const auto& fs = intQuants.fluidState();
+            v =
+                decay<Scalar>(fs.invB(FluidSystem::oilPhaseIdx))
+                * decay<Scalar>(extQuants.volumeFlux(FluidSystem::oilPhaseIdx))
+                * decay<Scalar>(fs.Rs());
+        }
+        else {
+            upIdx = 0;
+            v = 0.0;
+        }
+
+        Scalar A = scvf.area();
+        if (inIdx == upIdx) {
+            solFlux = A*v*variable<TracerEvaluation>(1.0, 0);
+            isUp = true;
+        }
+        else {
+            solFlux = A*v;
+            isUp = false;
         }
     }
 
@@ -290,29 +367,47 @@ protected:
     {
         if (tr.numTracer() == 0)
             return;
-
+        
         std::vector<Scalar> storageOfTimeIndex1(tr.numTracer());
+        std::vector<Scalar> fStorageOfTimeIndex1(tr.numTracer());
+        std::vector<Scalar> sStorageOfTimeIndex1(tr.numTracer());
         if (elemCtx.enableStorageCache()) {
             for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
-                storageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I];
+                fStorageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I][0];
+                sStorageOfTimeIndex1[tIdx] = tr.storageOfTimeIndex1_[tIdx][I][1];
             }
         }
         else {
-            Scalar fVolume1;
+            // ///
-            computeVolume_(fVolume1, tr.phaseIdx_, elemCtx, 0, /*timeIdx=*/1);
+            // OBS: below code will give runtime error since we cannot access intensive quantities at timeIdx=1
+            // //
+            Scalar fVolume1 = computeFreeVolume_(tr.phaseIdx_, I1, 1);
+            Scalar sVolume1 = computeSolutionVolume_(tr.phaseIdx_, I1, 1);
             for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
-                storageOfTimeIndex1[tIdx] = fVolume1*tr.concentrationInitial_[tIdx][I1];
+                fStorageOfTimeIndex1[tIdx] = fVolume1 * tr.concentration_[tIdx][I1][0];
+                sStorageOfTimeIndex1[tIdx] = sVolume1 * tr.concentration_[tIdx][I1][1];
             }
         }
 
-        TracerEvaluation fVolume;
+        TracerEvaluation fVol = computeFreeVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
-        computeVolume_(fVolume, tr.phaseIdx_, elemCtx, 0, /*timeIdx=*/0);
+        TracerEvaluation sVol = computeSolutionVolume_(tr.phaseIdx_, I, 0) * variable<TracerEvaluation>(1.0, 0);
         for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
-            Scalar storageOfTimeIndex0 = fVolume.value()*tr.concentration_[tIdx][I];
+            // Free part
-            Scalar localStorage = (storageOfTimeIndex0 - storageOfTimeIndex1[tIdx]) * scvVolume/dt;
+            Scalar fStorageOfTimeIndex0 = fVol.value() * tr.concentration_[tIdx][I][0];
-            tr.residual_[tIdx][I][0] += localStorage; //residual + flux
+            Scalar fLocalStorage = (fStorageOfTimeIndex0 - fStorageOfTimeIndex1[tIdx]) * scvVolume/dt;
+            tr.residual_[tIdx][I][0] += fLocalStorage; //residual + flux
+
+            // Solution part
+            Scalar sStorageOfTimeIndex0 = sVol.value() * tr.concentration_[tIdx][I][1];
+            Scalar sLocalStorage = (sStorageOfTimeIndex0 - sStorageOfTimeIndex1[tIdx]) * scvVolume/dt;
+            tr.residual_[tIdx][I][1] += sLocalStorage; //residual + flux
+
         }
-        (*tr.mat)[I][I][0][0] += fVolume.derivative(0) * scvVolume/dt;
+
+        // Derivative matrix
+        (*tr.mat)[I][I][0][0] += fVol.derivative(0) * scvVolume/dt;
+        (*tr.mat)[I][I][1][1] += sVol.derivative(0) * scvVolume/dt;
+
     }
 
     template<class TrRe>
@@ -325,16 +420,28 @@ protected:
         if (tr.numTracer() == 0)
             return;
 
-        TracerEvaluation flux;
+        TracerEvaluation fFlux;
+        TracerEvaluation sFlux;
         bool isUpF;
-        computeFlux_(flux, isUpF, tr.phaseIdx_, elemCtx, scvfIdx, 0);
+        bool isUpS;
-        int globalUpIdx = isUpF ? I : J;
+        computeFreeFlux_(fFlux, isUpF, tr.phaseIdx_, elemCtx, scvfIdx, 0);
+        computeSolFlux_(sFlux, isUpS, tr.phaseIdx_, elemCtx, scvfIdx, 0);
+        int fGlobalUpIdx = isUpF ? I : J;
+        int sGlobalUpIdx = isUpS ? I : J;
         for (int tIdx =0; tIdx < tr.numTracer(); ++tIdx) {
-            tr.residual_[tIdx][I][0] += flux.value()*tr.concentration_[tIdx][globalUpIdx]; //residual + flux
+            // Free and solution fluxes
+            tr.residual_[tIdx][I][0] += fFlux.value()*tr.concentration_[tIdx][fGlobalUpIdx][0]; //residual + flux
+            tr.residual_[tIdx][I][1] += sFlux.value()*tr.concentration_[tIdx][sGlobalUpIdx][1]; //residual + flux
         }
-        if (isUpF) {
+
-            (*tr.mat)[J][I][0][0] = -flux.derivative(0);
+        // Derivative matrix
-            (*tr.mat)[I][I][0][0] += flux.derivative(0);
+        if (isUpF){
+            (*tr.mat)[J][I][0][0] = -fFlux.derivative(0);
+            (*tr.mat)[I][I][0][0] += fFlux.derivative(0);
+        }
+        if (isUpS) {
+            (*tr.mat)[J][I][1][1] = -sFlux.derivative(0);
+            (*tr.mat)[I][I][1][1] += sFlux.derivative(0);
         }
     }
 
@@ -355,24 +462,87 @@ protected:
             wtracer[tIdx] = this->currentConcentration_(eclWell, this->name(tr.idx_[tIdx]));
         }
 
-        for (auto& perfData : well.perforationData()) {
+        std::size_t well_index = simulator_.problem().wellModel().wellState().index(well.name()).value();
-            const auto I = perfData.cell_index;
+        const auto& ws = simulator_.problem().wellModel().wellState().well(well_index);
+        for (std::size_t i = 0; i < ws.perf_data.size(); ++i) {
+            const auto I = ws.perf_data.cell_index[i];
             Scalar rate = well.volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
-            if (rate > 0) {
+
+            Scalar rate_s;
+            if (tr.phaseIdx_ == FluidSystem::oilPhaseIdx && FluidSystem::enableVaporizedOil()) {
+                rate_s = ws.perf_data.phase_mixing_rates[i][ws.dissolved_gas];
+            }
+            else if (tr.phaseIdx_ == FluidSystem::gasPhaseIdx && FluidSystem::enableDissolvedGas()) {
+                rate_s = ws.perf_data.phase_mixing_rates[i][ws.vaporized_oil];
+            }
+            else {
+                rate_s = 0.0;
+            }
+
+            Scalar rate_f = rate - rate_s;
+
+            if (rate_f > 0) {
                 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
-                    tr.residual_[tIdx][I][0] -= rate*wtracer[tIdx];
+                    // Injection of free tracer only
+                    tr.residual_[tIdx][I][0] -= rate_f*wtracer[tIdx];
+
                     // Store _injector_ tracer rate for reporting
-                    this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) += rate*wtracer[tIdx];
+                    this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) += rate_f*wtracer[tIdx];
                 }
             }
-            else if (rate < 0) {
+            else if (rate_f < 0) {
                 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
-                    tr.residual_[tIdx][I][0] -= rate*tr.concentration_[tIdx][I];
+                    // Free and solution tracer production
+                    tr.residual_[tIdx][I][0] -= rate_f * tr.concentration_[tIdx][I][0];
+                    tr.residual_[tIdx][I][1] -= rate_s * tr.concentration_[tIdx][I][1];
+
+                    // Store _producer_ tracer rate for reporting
+                    this->wellTracerRate_.at(std::make_pair(eclWell.name(),this->name(tr.idx_[tIdx]))) += 
+                        rate_f * tr.concentration_[tIdx][I][0] + rate_s * tr.concentration_[tIdx][I][1];
                 }
-                (*tr.mat)[I][I][0][0] -= rate*variable<TracerEvaluation>(1.0, 0).derivative(0);
+
+                // Derivative matrix for producer
+                (*tr.mat)[I][I][0][0] -= rate_f * variable<TracerEvaluation>(1.0, 0).derivative(0);
+                (*tr.mat)[I][I][1][1] -= rate_s * variable<TracerEvaluation>(1.0, 0).derivative(0);
             }
         }
     }
+    
+    template<class TrRe>
+    void assembleTracerEquationSource(TrRe& tr,
+                                      const Scalar scvVolume,
+                                      const Scalar dt,
+                                      unsigned I)
+    {
+        if (tr.numTracer() == 0)
+            return;
+
+        // Diff dissolved phase volume
+        Scalar sVol0 = computeSolutionVolume_(tr.phaseIdx_, I, 0) * scvVolume;
+        Scalar dsVol = sVol0 - sVol1_[tr.phaseIdx_][I];
+
+        // Source term determined by sign of dsVol: if dsVol > 0 then ms -> mf, else mf -> ms
+        for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
+            if (dsVol >= 0) {
+                tr.residual_[tIdx][I][0] += (dsVol / dt) * tr.concentration_[tIdx][I][0];
+                tr.residual_[tIdx][I][1] -= (dsVol / dt) * tr.concentration_[tIdx][I][0];
+            }
+            else {
+                tr.residual_[tIdx][I][0] += (dsVol / dt) * tr.concentration_[tIdx][I][1];
+                tr.residual_[tIdx][I][1] -= (dsVol / dt) * tr.concentration_[tIdx][I][1];
+            }
+        }
+
+        // Derivative matrix
+        if (dsVol >= 0) {
+            (*tr.mat)[I][I][0][0] += (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
+            (*tr.mat)[I][I][1][0] -= (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
+        }
+        else {
+            (*tr.mat)[I][I][1][1] += (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
+            (*tr.mat)[I][I][0][1] -= (dsVol / dt) * variable<TracerEvaluation>(1.0, 0).derivative(0);
+        }
+    }                                      
 
     void assembleTracerEquations_()
     {
@@ -385,8 +555,9 @@ protected:
         for (auto& tr : tbatch) {
             if (tr.numTracer() != 0) {
                 (*tr.mat) = 0.0;
-                for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx)
+                for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
                     tr.residual_[tIdx] = 0.0;
+                }
             }
         }
 
@@ -400,8 +571,10 @@ protected:
             {
                 // Dirichlet boundary conditions needed for the parallel matrix
                 for (auto& tr : tbatch) {
-                    if (tr.numTracer() != 0)
+                    if (tr.numTracer() != 0) {
                         (*tr.mat)[I][I][0][0] = 1.;
+                        (*tr.mat)[I][I][1][1] = 1.;
+                    }
                 }
                 continue;
             }
@@ -433,6 +606,12 @@ protected:
                     this->assembleTracerEquationFlux(tr, elemCtx, scvfIdx, I, J);
                 }
             }
+            
+            // Source terms (mass transfer between free and solution tracer)
+            for (auto& tr : tbatch) {
+                this->assembleTracerEquationSource(tr, scvVolume, dt, I);
+            }
+
         }
 
         // Well terms
@@ -457,22 +636,29 @@ protected:
     void updateStorageCache()
     {
         for (auto& tr : tbatch) {
-            if (tr.numTracer() != 0)
+            if (tr.numTracer() != 0) {
                 tr.concentrationInitial_ = tr.concentration_;
+            }
         }
 
         ElementContext elemCtx(simulator_);
         for (const auto& elem : elements(simulator_.gridView())) {
             elemCtx.updatePrimaryStencil(elem);
             elemCtx.updatePrimaryIntensiveQuantities(/*timeIdx=*/0);
+            Scalar extrusionFactor = elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
+            Scalar scvVolume = elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume() * extrusionFactor;
             int globalDofIdx = elemCtx.globalSpaceIndex(0, /*timeIdx=*/0);
             for (auto& tr : tbatch) {
                 if (tr.numTracer() == 0)
                     continue;
-                Scalar fVolume;
+                
-                computeVolume_(fVolume, tr.phaseIdx_, elemCtx, 0, /*timeIdx=*/0);
+                Scalar fVol1 = computeFreeVolume_(tr.phaseIdx_, globalDofIdx, 0);
+                Scalar sVol1 = computeSolutionVolume_(tr.phaseIdx_, globalDofIdx, 0);
+                fVol1_[tr.phaseIdx_][globalDofIdx] = fVol1 * scvVolume;
+                sVol1_[tr.phaseIdx_][globalDofIdx] = sVol1 * scvVolume;
                 for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
-                    tr.storageOfTimeIndex1_[tIdx][globalDofIdx] = fVolume*tr.concentrationInitial_[tIdx][globalDofIdx];
+                    tr.storageOfTimeIndex1_[tIdx][globalDofIdx][0] = fVol1 * tr.concentrationInitial_[tIdx][globalDofIdx][0];
+                    tr.storageOfTimeIndex1_[tIdx][globalDofIdx][1] = sVol1 * tr.concentrationInitial_[tIdx][globalDofIdx][1];
                 }
             }
         }
@@ -489,8 +675,9 @@ protected:
             // Note that we solve for a concentration update (compared to previous time step)
             // Confer also assembleTracerEquations_(...) above.
             std::vector<TracerVector> dx(tr.concentration_);
-            for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx)
+            for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
                 dx[tIdx] = 0.0;
+            }
 
             bool converged = this->linearSolveBatchwise_(*tr.mat, dx, tr.residual_);
             if (!converged) {
@@ -498,9 +685,14 @@ protected:
             }
 
             for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
+                // New concetration
                 tr.concentration_[tIdx] -= dx[tIdx];
-                // Tracer concentrations for restart report
+                
-                this->tracerConcentration_[tr.idx_[tIdx]] = tr.concentration_[tIdx];
+                // Partition concentration into free and solution tracers for output
+                for (std::size_t globalDofIdx = 0; globalDofIdx < tr.concentration_[tIdx].size(); ++globalDofIdx) {
+                    this->freeTracerConcentration_[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][0];
+                    this->solTracerConcentration_[tr.idx_[tIdx]][globalDofIdx] = tr.concentration_[tIdx][globalDofIdx][1];
+                }
             }
 
             // Store _producer_ tracer rate for reporting
@@ -518,9 +710,6 @@ protected:
                     Scalar rate = wellPtr->volumetricSurfaceRateForConnection(I, tr.phaseIdx_);
                     if (rate < 0) {
                         rateWellNeg += rate;
-                        for (int tIdx = 0; tIdx < tr.numTracer(); ++tIdx) {
-                            this->wellTracerRate_.at(std::make_pair(well.name(),this->name(tr.idx_[tIdx]))) += rate*tr.concentration_[tIdx][I];
-                        }
                     }
                     else {
                         rateWellPos += rate;
@@ -548,6 +737,7 @@ protected:
         }
     }
 
+    
 
     Simulator& simulator_;
 
@@ -561,13 +751,13 @@ protected:
 
     template <typename TV>
     struct TracerBatch {
-      std::vector<int> idx_;
+        std::vector<int> idx_;
-      const int phaseIdx_;
+        const int phaseIdx_;
-      std::vector<TV> concentrationInitial_;
+        std::vector<TV> concentrationInitial_;
-      std::vector<TV> concentration_;
+        std::vector<TV> concentration_;
-      std::vector<TV> storageOfTimeIndex1_;
+        std::vector<TV> storageOfTimeIndex1_;
-      std::vector<TV> residual_;
+        std::vector<TV> residual_;
-      std::unique_ptr<TracerMatrix> mat;
+        std::unique_ptr<TracerMatrix> mat;
 
       bool operator==(const TracerBatch& rhs) const
       {
@@ -600,12 +790,12 @@ protected:
 
       void addTracer(const int idx, const TV & concentration)
       {
-          int numGridDof = concentration.size();
+            int numGridDof = concentration.size();
-          idx_.emplace_back(idx);
+            idx_.emplace_back(idx);
-          concentrationInitial_.emplace_back(concentration);
+            concentrationInitial_.emplace_back(concentration);
-          concentration_.emplace_back(concentration);
+            concentration_.emplace_back(concentration);
-          storageOfTimeIndex1_.emplace_back(numGridDof);
+            residual_.emplace_back(numGridDof);
-          residual_.emplace_back(numGridDof);
+            storageOfTimeIndex1_.emplace_back(numGridDof);
       }
     };
 
@@ -613,6 +803,8 @@ protected:
     TracerBatch<TracerVector>& wat_;
     TracerBatch<TracerVector>& oil_;
     TracerBatch<TracerVector>& gas_;
+    std::array<std::vector<Scalar>, 3> fVol1_;
+    std::array<std::vector<Scalar>, 3> sVol1_;
 };
 
 } // namespace Opm

opm/simulators/flow/VtkTracerModule.hpp

@@ -107,10 +107,12 @@ namespace Opm {
         {
             if (eclTracerConcentrationOutput_()){
                 const auto& tracerModel = this->simulator_.problem().tracerModel();
-                eclTracerConcentration_.resize(tracerModel.numTracers());
+                eclFreeTracerConcentration_.resize(tracerModel.numTracers());
-                for (std::size_t tracerIdx = 0; tracerIdx < eclTracerConcentration_.size(); ++tracerIdx) {
+                eclSolTracerConcentration_.resize(tracerModel.numTracers());
+                for (std::size_t tracerIdx = 0; tracerIdx < eclFreeTracerConcentration_.size(); ++tracerIdx) {
 
-                    this->resizeScalarBuffer_(eclTracerConcentration_[tracerIdx]);
+                    this->resizeScalarBuffer_(eclFreeTracerConcentration_[tracerIdx]);
+                    this->resizeScalarBuffer_(eclSolTracerConcentration_[tracerIdx]);
                 }
             }
 
@@ -131,8 +133,9 @@ namespace Opm {
                 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
 
                 if (eclTracerConcentrationOutput_()){
-                    for (std::size_t tracerIdx  = 0; tracerIdx < eclTracerConcentration_.size(); ++tracerIdx) {
+                    for (std::size_t tracerIdx  = 0; tracerIdx < eclFreeTracerConcentration_.size(); ++tracerIdx) {
-                        eclTracerConcentration_[tracerIdx][globalDofIdx] = tracerModel.tracerConcentration(tracerIdx, globalDofIdx);
+                        eclFreeTracerConcentration_[tracerIdx][globalDofIdx] = tracerModel.freeTracerConcentration(tracerIdx, globalDofIdx);
+                        eclSolTracerConcentration_[tracerIdx][globalDofIdx] = tracerModel.solTracerConcentration(tracerIdx, globalDofIdx);
                     }
                 }
             }
@@ -149,9 +152,11 @@ namespace Opm {
 
             if (eclTracerConcentrationOutput_()){
                 const auto& tracerModel = this->simulator_.problem().tracerModel();
-                for (std::size_t tracerIdx = 0; tracerIdx < eclTracerConcentration_.size(); ++tracerIdx) {
+                for (std::size_t tracerIdx = 0; tracerIdx < eclFreeTracerConcentration_.size(); ++tracerIdx) {
-                    const std::string tmp = "tracerConcentration_" + tracerModel.name(tracerIdx);
+                    const std::string tmp = "freeTracerConcentration_" + tracerModel.name(tracerIdx);
-                    this->commitScalarBuffer_(baseWriter,tmp.c_str(), eclTracerConcentration_[tracerIdx]);
+                    this->commitScalarBuffer_(baseWriter,tmp.c_str(), eclFreeTracerConcentration_[tracerIdx]);
+                    const std::string tmp2 = "solTracerConcentration_" + tracerModel.name(tracerIdx);
+                    this->commitScalarBuffer_(baseWriter,tmp2.c_str(), eclSolTracerConcentration_[tracerIdx]);
                 }
             }
 
@@ -167,7 +172,8 @@ namespace Opm {
         }
 
 
-        std::vector<ScalarBuffer> eclTracerConcentration_;
+        std::vector<ScalarBuffer> eclFreeTracerConcentration_;
+        std::vector<ScalarBuffer> eclSolTracerConcentration_;
     };
 } // namespace Opm