Partitioning tracers

opm/simulators/flow/EclWriter.hpp

@@ -500,8 +500,12 @@ class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>
 
         {
             const auto& tracers = simulator_.vanguard().eclState().tracer();
-            for (const auto& tracer : tracers)
+            for (const auto& tracer : tracers) {
+                bool enableSolTracer = (tracer.phase == Phase::GAS && FluidSystem::enableDissolvedGas()) ||
+                                       (tracer.phase == Phase::OIL && FluidSystem::enableVaporizedOil());
                 solutionKeys.emplace_back(tracer.fname(), UnitSystem::measure::identity, true);
+                solutionKeys.emplace_back(tracer.sname(), UnitSystem::measure::identity, enableSolTracer);
+            }
         }
 
         // The episodeIndex is rewined one back before beginRestart is called
@@ -527,11 +531,30 @@ class EclWriter : public EclGenericWriter<GetPropType<TypeTag, Properties::Grid>
 
             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& tracer_solution = restartValues.solution.template data<double>(tracer_name);
+                // Free tracers
+                const auto& free_tracer_name = tracer_model.fname(tracer_index);
+                const auto& free_tracer_solution = restartValues.solution.template data<double>(free_tracer_name);
                 for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
                     unsigned globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
-                    tracer_model.setTracerConcentration(tracer_index, globalIdx, tracer_solution[globalIdx]);
+                    tracer_model.setFreeTracerConcentration(tracer_index, globalIdx, free_tracer_solution[globalIdx]);
+                }
+                // Solution tracer (only if DISGAS/VAPOIL are active for gas/oil tracers)
+                if ((tracer_model.phase(tracer_index) == Phase::GAS && FluidSystem::enableDissolvedGas()) ||
+                    (tracer_model.phase(tracer_index) == Phase::OIL && FluidSystem::enableVaporizedOil())) {
+                        tracer_model.setEnableSolTracers(tracer_index, true);
+                        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);
+                            tracer_model.setSolTracerConcentration(tracer_index, globalIdx, sol_tracer_solution[globalIdx]);
+                        }
+                }
+                else {
+                    tracer_model.setEnableSolTracers(tracer_index, false);
+                    for (unsigned elemIdx = 0; elemIdx < numElements; ++elemIdx) {
+                            unsigned globalIdx = this->collectOnIORank_.localIdxToGlobalIdx(elemIdx);
+                            tracer_model.setSolTracerConcentration(tracer_index, globalIdx, 0.0);
+                        }
                 }
             }
 

opm/simulators/flow/GenericOutputBlackoilModule.cpp

@@ -700,21 +700,40 @@ 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)
+        {
+            if (solTracerConcentrations_[tracerIdx].empty())
+                continue;
+
+            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();
     }
 }
 
@@ -838,6 +857,7 @@ doAllocBuffers(const unsigned bufferSize,
                const bool     vapparsActive,
                const bool     enableHysteresis,
                const unsigned numTracers,
+               const std::vector<bool>& enableSolTracers,
                const unsigned numOutputNnc)
 {
     // Output RESTART_OPM_EXTENDED only when explicitly requested by user.
@@ -1304,10 +1324,16 @@ doAllocBuffers(const unsigned bufferSize,
 
     // tracers
     if (numTracers > 0) {
-        tracerConcentrations_.resize(numTracers);
+        freeTracerConcentrations_.resize(numTracers);
+        for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx)
+        {
+            freeTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0);
+        }
+        solTracerConcentrations_.resize(numTracers);
         for (unsigned tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx)
         {
-            tracerConcentrations_[tracerIdx].resize(bufferSize, 0.0);
+            if (enableSolTracers[tracerIdx])
+                solTracerConcentrations_[tracerIdx].resize(bufferSize, 0.0);
         }
     }
 

opm/simulators/flow/GenericOutputBlackoilModule.hpp

@@ -340,6 +340,7 @@ class GenericOutputBlackoilModule {
                         const bool vapparsActive,
                         const bool enableHysteresis,
                         unsigned numTracers,
+                        const std::vector<bool>& enableSolTracers,
                         unsigned numOutputNnc);
 
     void makeRegionSum(Inplace& inplace,
@@ -504,7 +505,8 @@ class GenericOutputBlackoilModule {
     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.cpp

@@ -24,6 +24,8 @@
 
 #include <opm/simulators/flow/GenericTracerModel_impl.hpp>
 
+#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
+
 #if HAVE_DUNE_FEM
 #include <dune/common/version.hh>
 #include <dune/fem/gridpart/adaptiveleafgridpart.hh>
@@ -39,6 +41,7 @@ template class GenericTracerModel<Dune::CpGrid,
                                   Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,
                                   Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>>,
                                   Opm::EcfvStencil<double,Dune::GridView<Dune::DefaultLeafGridViewTraits<Dune::CpGrid>>,false,false>,
+                                  BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
                                   double>;
 
 #if HAVE_DUNE_FEM
@@ -50,12 +53,14 @@ template class GenericTracerModel<Dune::CpGrid,
                                   GV,
                                   Dune::MultipleCodimMultipleGeomTypeMapper<GV>,
                                   EcfvStencil<double, GV, false, false>,
+                                  BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
                                   double>;
 #else
 template class GenericTracerModel<Dune::CpGrid,
                                   Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,
                                   Dune::MultipleCodimMultipleGeomTypeMapper<Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>>,
                                   EcfvStencil<double,Dune::GridView<Dune::Fem::GridPart2GridViewTraits<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false>>>,false,false>,
+                                  BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
                                   double>;
 template class GenericTracerModel<Dune::CpGrid,
                                   Dune::Fem::GridPart2GridViewImpl<Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, (Dune::PartitionIteratorType)4, false> >,
@@ -65,6 +70,7 @@ template class GenericTracerModel<Dune::CpGrid,
                                   EcfvStencil<double, Dune::Fem::GridPart2GridViewImpl<
                                                                Dune::Fem::AdaptiveLeafGridPart<Dune::CpGrid, Dune::PartitionIteratorType(4), false> >,
                                                    false, false>,
+                                  BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,
                                   double>;
 #endif
 #endif // HAVE_DUNE_FEM

opm/simulators/flow/GenericTracerModel.hpp

@@ -36,6 +36,8 @@
 
 #include <opm/simulators/linalg/matrixblock.hh>
 
+#include <opm/input/eclipse/EclipseState/Phase.hpp>
+
 #include <array>
 #include <cstddef>
 #include <functional>
@@ -49,11 +51,12 @@ namespace Opm {
 class EclipseState;
 class Well;
 
-template<class Grid, class GridView, class DofMapper, class Stencil, class Scalar>
+template<class Grid, class GridView, class DofMapper, class Stencil, class FluidSystem, class Scalar>
 class GenericTracerModel {
 public:
-    using TracerMatrix = Dune::BCRSMatrix<Opm::MatrixBlock<Scalar, 1, 1>>;
-    using TracerVector = Dune::BlockVector<Dune::FieldVector<Scalar,1>>;
+    using TracerVectorSingle = 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,25 +69,44 @@ class GenericTracerModel {
      */
     const std::string& name(int tracerIdx) const;
     std::string fname(int tracerIdx) const;
+    std::string sname(int tracerIdx) const;
+    std::string wellfname(int tracerIdx) const;
+    std::string wellsname(int tracerIdx) const;
 
+    Phase phase(int tracerIdx) const;
+    const std::vector<bool>& enableSolTracers() const;
 
     /*!
      * \brief Return the tracer concentration for tracer index and global DofIdx
      */
-    Scalar tracerConcentration(int tracerIdx, int globalDofIdx) const;
-    void setTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value);
+    Scalar freeTracerConcentration(int tracerIdx, int globalDofIdx) const;
+    Scalar solTracerConcentration(int tracerIdx, int globalDofIdx) const;
+    void setFreeTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value);
+    void setSolTracerConcentration(int tracerIdx, int globalDofIdx, Scalar value);
+    void setEnableSolTracers(int tracerIdx, bool enableSolTracer);
 
     /*!
     * \brief Return well tracer rates
     */
     const std::map<std::pair<std::string, std::string>, Scalar>&
     getWellTracerRates() const {return wellTracerRate_;}
+    const std::map<std::pair<std::string, std::string>, Scalar>&
+    getWellFreeTracerRates() const {return wellFreeTracerRate_;}
+    const std::map<std::pair<std::string, std::string>, Scalar>&
+    getWellSolTracerRates() const {return wellSolTracerRate_;}
+    const std::map<std::tuple<std::string, std::string, std::size_t>, Scalar>&
+    getMswTracerRates() const {return mSwTracerRate_;}
 
     template<class Serializer>
     void serializeOp(Serializer& serializer)
     {
         serializer(tracerConcentration_);
+        serializer(freeTracerConcentration_);
+        serializer(solTracerConcentration_);
         serializer(wellTracerRate_);
+        serializer(wellFreeTracerRate_);
+        serializer(wellSolTracerRate_);
+        serializer(mSwTracerRate_);
     }
 
 protected:
@@ -117,12 +139,20 @@ class GenericTracerModel {
     const DofMapper& dofMapper_;
 
     std::vector<int> tracerPhaseIdx_;
-    std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> tracerConcentration_;
+    std::vector<bool> enableSolTracers_;
+    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
+    // <wellName, tracerName> -> wellRate
     std::map<std::pair<std::string, std::string>, Scalar> wellTracerRate_;
+    std::map<std::pair<std::string, std::string>, Scalar> wellFreeTracerRate_;
+    std::map<std::pair<std::string, std::string>, Scalar> wellSolTracerRate_;
+
+    // <wellName, tracerName, segNum> -> wellRate
+    std::map<std::tuple<std::string, std::string, std::size_t>, Scalar> mSwTracerRate_;
+
     /// \brief Function returning the cell centers
     std::function<std::array<double,dimWorld>(int)> centroids_;
 };