Add ecl tracer model

reads tracer input from deck, solves tracer equation fully implicit as a post processing step in endTimeStep

tested on a simple modified SPE1CASE1 deck and compared with eclipse

TODO: restart and parallel

ebos/ecloutputblackoilmodule.hh

@@ -321,6 +321,16 @@ public:
             dewPointPressure_.resize(bufferSize, 0.0);
         }
 
+        // tracers
+        const int numTracers = simulator_.problem().tracerModel().numTracers();
+        if (numTracers > 0){
+            tracerConcentrations_.resize(numTracers);
+            for (int tracerIdx = 0; tracerIdx < numTracers; ++tracerIdx)
+            {
+                tracerConcentrations_[tracerIdx].resize(bufferSize, 0.0);
+            }
+        }
+
         //Warn for any unhandled keyword
         if (log) {
             for (auto& keyValue: rstKeywords) {
@@ -341,6 +351,7 @@ public:
             saturatedOilFormationVolumeFactor_.resize(bufferSize, 0.0);
         if (false)
             oilSaturationPressure_.resize(bufferSize, 0.0);
+
     }
 
     /*!
@@ -581,6 +592,17 @@ public:
             if (gasConnectionSaturations_.count(cartesianIdx) > 0) {
                 gasConnectionSaturations_[cartesianIdx] = Opm::getValue(fs.saturation(gasPhaseIdx));
             }
+
+            // tracers
+            const auto& tracerModel = simulator_.problem().tracerModel();
+            if (tracerConcentrations_.size()>0) {
+                for (int tracerIdx = 0; tracerIdx < tracerModel.numTracers(); tracerIdx++){
+                    if (tracerConcentrations_[tracerIdx].size() == 0)
+                        continue;
+
+                    tracerConcentrations_[tracerIdx][globalDofIdx] = tracerModel.tracerConcentration(tracerIdx, globalDofIdx);
+                }
+            }
         }
     }
 
@@ -814,6 +836,15 @@ public:
                            Opm::data::TargetType::SUMMARY);
             }
         }
+
+        // tracers
+        const auto& tracerModel = simulator_.problem().tracerModel();
+        if (tracerConcentrations_.size()>0) {
+            for (int tracerIdx = 0; tracerIdx<tracerModel.numTracers(); tracerIdx++){
+                std::string tmp = tracerModel.tracerName(tracerIdx) + "F";
+                sol.insert (tmp, Opm::UnitSystem::measure::identity, std::move(tracerConcentrations_[tracerIdx]), Opm::data::TargetType::RESTART_SOLUTION);
+            }
+        }
     }
 
     // write Fluid In Place to output log
@@ -912,7 +943,7 @@ public:
 
     }
 
-    void setRestart(const Opm::data::Solution& sol, unsigned elemIdx, unsigned globalDofIndex)
+    void setRestart(const Opm::data::Solution& sol, unsigned elemIdx, unsigned globalDofIndex) 
     {
         Scalar so = 1.0;
         if (saturation_[waterPhaseIdx].size() > 0 && sol.has("SWAT")) {
@@ -1356,6 +1387,7 @@ private:
     std::map<size_t, Scalar> oilConnectionPressures_;
     std::map<size_t, Scalar> waterConnectionSaturations_;
     std::map<size_t, Scalar> gasConnectionSaturations_;
+    std::vector<ScalarBuffer> tracerConcentrations_;
 };
 } // namespace Ewoms
 

ebos/eclpeacemanwell.hh

@@ -487,6 +487,15 @@ public:
         matrix.setBlock(wellGlobalDofIdx, wellGlobalDofIdx, diagBlock);
     }
 
+    Scalar volumetricSurfaceRateForConnection(int globalDofIdx, int phaseIdx) const {
+        const DofVariables& dofVars = *dofVariables_.at(globalDofIdx);
+        std::array<Scalar, numPhases> volumetricReservoirRates;
+        computeVolumetricDofRates_(volumetricReservoirRates, actualBottomHolePressure_, dofVars);
+        std::array<Scalar, numPhases> volumetricSurfaceRates;
+        computeSurfaceRates_(volumetricSurfaceRates, volumetricReservoirRates, dofVars);
+        return volumetricSurfaceRates[phaseIdx];
+    }
+
 
     // reset the well to the initial state, i.e. remove all degrees of freedom...
     void clear()

ebos/eclproblem.hh

@@ -60,6 +60,8 @@
 #include "ecldummygradientcalculator.hh"
 #include "eclfluxmodule.hh"
 #include "eclbaseaquifermodel.hh"
+#include "ecltracermodel.hh"
+#include "vtkecltracermodule.hh"
 
 #include <ewoms/common/pffgridvector.hh>
 #include <ewoms/models/blackoil/blackoilmodel.hh>
@@ -108,7 +110,7 @@ BEGIN_PROPERTIES
 #if EBOS_USE_ALUGRID
 NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclAluGridVanguard, EclOutputBlackOil));
 #else
-NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclCpGridVanguard, EclOutputBlackOil));
+NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclCpGridVanguard, EclOutputBlackOil, VtkEclTracer));
 //NEW_TYPE_TAG(EclBaseProblem, INHERITS_FROM(EclPolyhedralGridVanguard, EclOutputBlackOil));
 #endif
 
@@ -369,6 +371,8 @@ class EclProblem : public GET_PROP_TYPE(TypeTag, BaseProblem)
 
     typedef EclWriter<TypeTag> EclWriterType;
 
+    typedef EclTracerModel<TypeTag> TracerModel;
+
     typedef typename GridView::template Codim<0>::Iterator ElementIterator;
 
     struct RockParams {
@@ -384,6 +388,7 @@ public:
     {
         ParentType::registerParameters();
         EclWriterType::registerParameters();
+        VtkEclTracerModule<TypeTag>::registerParameters();
 
         EWOMS_REGISTER_PARAM(TypeTag, bool, EnableWriteAllSolutions,
                              "Write all solutions to disk instead of only the ones for the "
@@ -492,7 +497,9 @@ public:
         , wellModel_(simulator)
         , aquiferModel_(simulator)
         , pffDofData_(simulator.gridView(), this->elementMapper())
+        , tracerModel_(simulator)
     {
+        this->model().addOutputModule(new VtkEclTracerModule<TypeTag>(simulator));
         // Tell the black-oil extensions to initialize their internal data structures
         const auto& vanguard = simulator.vanguard();
         SolventModule::initFromDeck(vanguard.deck(), vanguard.eclState());
@@ -579,6 +586,8 @@ public:
             eclWriter_->writeInit();
             this->simulator().vanguard().releaseGlobalTransmissibilities();
         }
+
+        tracerModel_.init();
     }
 
     void prefetch(const Element& elem) const
@@ -722,6 +731,8 @@ public:
         }
 
         aquiferModel_.beginTimeStep();
+        tracerModel_.beginTimeStep();
+
     }
 
     /*!
@@ -765,6 +776,8 @@ public:
             wellModel_.endTimeStep();
 
         aquiferModel_.endTimeStep();
+        tracerModel_.endTimeStep();
+
 
         // we no longer need the initial soluiton
         if (this->simulator().episodeIndex() == 0 && !initialFluidStates_.empty())  {
@@ -950,6 +963,9 @@ public:
     EclThresholdPressure<TypeTag>& thresholdPressure()
     { return thresholdPressures_; }
 
+    const EclTracerModel<TypeTag>& tracerModel() const
+    { return tracerModel_; }
+
     /*!
      * \copydoc FvBaseMultiPhaseProblem::porosity
      */
@@ -1779,6 +1795,10 @@ private:
             if (enablePolymer)
                  polymerConcentration_[elemIdx] = eclWriter_->eclOutputModule().getPolymerConcentration(elemIdx);
         }
+
+        if (tracerModel().numTracers() > 0)
+            std::cout << "Warning: Restart is not implemented for the tracer model, it will initialize with initial tracer concentration" << std::endl;
+
     }
 
     void readExplicitInitialCondition_()
@@ -2128,10 +2148,11 @@ private:
 
     EclWellModel wellModel_;
     EclAquiferModel aquiferModel_;
-
     std::unique_ptr<EclWriterType> eclWriter_;
 
     PffGridVector<GridView, Stencil, PffDofData_, DofMapper> pffDofData_;
+    TracerModel tracerModel_;
+
 
 };
 

ebos/ecltracermodel.hh

@@ -0,0 +1,506 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/**
+ * \file
+ *
+ * \copydoc Ewoms::EclTracerModel
+ */
+#ifndef EWOMS_ECL_TRACER_MODEL_HH
+#define EWOMS_ECL_TRACER_MODEL_HH
+#include "tracervdtable.hh"
+#include <dune/istl/operators.hh>
+#include <dune/istl/solvers.hh>
+#include <dune/istl/preconditioners.hh>
+#include <string>
+#include <vector>
+#include <iostream>
+
+
+namespace Ewoms {
+
+/*!
+ * \ingroup EclBlackOilSimulator
+ *
+ * \brief A class which handles tracers as specified in ecl deck
+ */
+template <class TypeTag>
+class EclTracerModel
+{
+    typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Stencil) Stencil;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+    typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+
+    typedef Opm::DenseAd::Evaluation<Scalar,1> TracerEvaluation;
+
+    enum { numEq = GET_PROP_VALUE(TypeTag, NumEq) };
+    enum { numPhases = FluidSystem::numPhases };
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
+    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+
+    typedef Dune::BCRSMatrix<Dune::FieldMatrix<Scalar, 1, 1>> TracerMatrix;
+    typedef Dune::BlockVector< Dune::FieldVector<Scalar,1>> TracerVector;
+
+public:
+    EclTracerModel(Simulator& simulator)
+        : simulator_(simulator)
+    { }
+
+    /*!
+     * \brief Initialize all internal data structures needed by the tracer module
+     */
+    void init()
+    {
+        const Opm::Deck& deck = simulator_.vanguard().deck();
+
+        if (!deck.hasKeyword("TRACERS"))
+            return; // tracer treatment is supposed to be disabled
+
+        if (!deck.hasKeyword("TRACER")){
+            throw std::runtime_error("the deck does not contain the TRACER keyword");
+        }
+        if (simulator_.gridView().comm().size() > 1) {
+            tracerNames_.resize(0);
+            std::cout << "Warning: Tracer model is not compatible with mpi run"  << std::endl;
+            return;
+        }
+
+        //how many tracers?
+        const int numTracers = deck.getKeyword("TRACER").size();
+        tracerNames_.resize(numTracers);
+        tracerConcentration_.resize(numTracers);
+        wTracer_.resize(numTracers,0.0);
+        storageOfTimeIndex1_.resize(numTracers);
+        size_t numAllDof =  simulator_.model().numTotalDof();
+
+        // the phase where the tracer is
+        tracerPhaseIdx_.resize(numTracers);
+        for (int tracerIdx = 0;  tracerIdx < numTracers; ++tracerIdx) {
+            const auto& tracerRecord = deck.getKeyword("TRACER").getRecord(tracerIdx);
+            tracerNames_[tracerIdx] = tracerRecord.getItem("NAME").template get<std::string>(0);
+            const std::string& fluidName = tracerRecord.getItem("FLUID").template get<std::string>(0);
+            if (fluidName == "WAT")
+                tracerPhaseIdx_[tracerIdx] = waterPhaseIdx;
+            else if (fluidName == "OIL")
+                tracerPhaseIdx_[tracerIdx] = oilPhaseIdx;
+            else if (fluidName == "GAS")
+                tracerPhaseIdx_[tracerIdx] = gasPhaseIdx;
+            else
+                throw std::invalid_argument("Tracer: invalid fluid name "
+                                            +fluidName+" for "+tracerNames_[tracerIdx]);
+
+            tracerConcentration_[tracerIdx].resize(numAllDof);
+            storageOfTimeIndex1_[tracerIdx].resize(numAllDof);
+            std::string tmp = "TVDPF" +tracerNames_[tracerIdx];
+
+
+            //TBLK keyword
+            if (deck.hasKeyword("TBLKF" +tracerNames_[tracerIdx])){
+                const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper();
+                const auto& tblkData =
+                        deck.getKeyword("TBLKF"
+                                        +tracerNames_
+                                        [tracerIdx]).getRecord(0).getItem(0).getSIDoubleData();
+                int tblkDatasize = tblkData.size();
+                if (tblkDatasize < simulator_.vanguard().cartesianSize()){
+                    throw std::runtime_error("Uninitialized tracer concentration (TBLKF) for tracer "
+                                             + tracerName(tracerIdx));
+                }
+                for (size_t globalDofIdx = 0; globalDofIdx < numAllDof; ++globalDofIdx){
+                    int cartDofIdx = cartMapper.cartesianIndex(globalDofIdx);
+                    tracerConcentration_[tracerIdx][globalDofIdx] = tblkData[cartDofIdx];
+                }
+            }
+            //TVDPF keyword
+            else if (deck.hasKeyword(tmp)){
+                TracerVdTable dtable(deck.getKeyword(tmp).getRecord(0).getItem(0));
+                const auto& eclGrid = simulator_.vanguard().eclState().getInputGrid();
+                const auto& cartMapper = simulator_.vanguard().cartesianIndexMapper();
+
+                for (size_t globalDofIdx = 0; globalDofIdx < numAllDof; ++globalDofIdx){
+                    int cartDofIdx = cartMapper.cartesianIndex(globalDofIdx);
+                    const auto& center = eclGrid.getCellCenter(cartDofIdx);
+                    tracerConcentration_[tracerIdx][globalDofIdx]
+                            = dtable.evaluate("TRACER_CONCENTRATION", center[2]);
+                }
+            }
+            else {
+                throw std::runtime_error("Uninitialized tracer concentration for tracer "
+                                         + tracerName(tracerIdx));
+            }
+
+        }
+
+        //initialize tracerConcentration
+        tracerConcentration0_ = tracerConcentration_;
+
+
+        //tracerResidual
+        tracerResidual_.resize(numAllDof);
+        // allocate raw matrix
+        tracerMatrix_ = new TracerMatrix(numAllDof, numAllDof, TracerMatrix::random);
+
+        Stencil stencil(simulator_.gridView(), simulator_.model().dofMapper() );
+
+        // for the main model, find out the global indices of the neighboring degrees of
+        // freedom of each primary degree of freedom
+        typedef std::set<unsigned> NeighborSet;
+        std::vector<NeighborSet> neighbors(numAllDof);
+
+        ElementIterator elemIt = simulator_.gridView().template begin<0>();
+        const ElementIterator elemEndIt = simulator_.gridView().template end<0>();
+        for (; elemIt != elemEndIt; ++elemIt) {
+            const Element& elem = *elemIt;
+            stencil.update(elem);
+
+            for (unsigned primaryDofIdx = 0; primaryDofIdx < stencil.numPrimaryDof(); ++primaryDofIdx)
+            {
+                unsigned myIdx = stencil.globalSpaceIndex(primaryDofIdx);
+
+                for (unsigned dofIdx = 0; dofIdx < stencil.numDof(); ++dofIdx) {
+                    unsigned neighborIdx = stencil.globalSpaceIndex(dofIdx);
+                    neighbors[myIdx].insert(neighborIdx);
+                }
+            }
+        }
+
+        // allocate space for the rows of the matrix
+        for (unsigned dofIdx = 0; dofIdx < numAllDof; ++ dofIdx)
+            tracerMatrix_->setrowsize(dofIdx, neighbors[dofIdx].size());
+        tracerMatrix_->endrowsizes();
+
+        // fill the rows with indices. each degree of freedom talks to
+        // all of its neighbors. (it also talks to itself since
+        // degrees of freedom are sometimes quite egocentric.)
+        for (unsigned dofIdx = 0; dofIdx < numAllDof; ++ dofIdx) {
+            typename NeighborSet::iterator nIt = neighbors[dofIdx].begin();
+            typename NeighborSet::iterator nEndIt = neighbors[dofIdx].end();
+            for (; nIt != nEndIt; ++nIt)
+                tracerMatrix_->addindex(dofIdx, *nIt);
+        }
+        tracerMatrix_->endindices();
+
+        const int sizeCartGrid = simulator_.vanguard().cartesianSize();
+        cartToGlobal_.resize(sizeCartGrid);
+        for (unsigned i = 0; i < numAllDof; ++i) {
+            int cartIdx = simulator_.vanguard().cartesianIndex(i);
+            cartToGlobal_[cartIdx] = i;
+        }
+
+    }
+
+    /*!
+     * \brief Return the number of tracers considered by the tracerModel.
+     */
+    int numTracers() const
+    { return tracerNames_.size(); }
+
+    /*!
+     * \brief Return the tracer name
+     */
+    const std::string& tracerName(int tracerIdx) const    
+    {
+        if (numTracers()==0)
+            throw std::logic_error("This method should never be called when there are no tracers in the model");
+
+        return tracerNames_[tracerIdx];
+    }
+
+    /*!
+     * \brief Return the tracer concentration for tracer index and global DofIdx
+     */
+    Scalar tracerConcentration(int tracerIdx, int globalDofIdx) const {
+        if (numTracers()==0)
+            return 0.0;
+
+        return tracerConcentration_[tracerIdx][globalDofIdx];
+    }
+
+    void beginTimeStep()
+    {
+
+        if (numTracers()==0)
+            return;
+
+        tracerConcentration0_ = tracerConcentration_;
+
+        if (!EWOMS_GET_PARAM(TypeTag, bool, EnableStorageCache))
+            return;
+
+        // 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 < numTracers(); ++ tracerIdx){
+                Scalar storageOfTimeIndex1;
+                computeStorage_(storageOfTimeIndex1, elemCtx, 0, /*timIdx=*/0, tracerIdx);
+                storageOfTimeIndex1_[tracerIdx][globalDofIdx] = storageOfTimeIndex1;
+            }
+        }
+    }
+
+    /*!
+     * \brief Informs the tracer model that a time step has just been finished.
+     */
+    void endTimeStep()
+    {
+        if (numTracers()==0)
+            return;
+
+        for (int tracerIdx = 0; tracerIdx < numTracers(); ++ tracerIdx){
+
+            TracerVector dx(tracerResidual_.size());
+            // Newton step (currently the system is linear, converge in one iteration)
+            for (int iter = 0; iter < 5; ++ iter){
+                linearize_(tracerIdx);
+                linearSolve_(*tracerMatrix_, dx, tracerResidual_);
+                tracerConcentration_[tracerIdx] -= dx;
+
+                if (dx.two_norm()<1e-2)
+                    break;
+            }
+        }
+    }
+
+    /*!
+     * \brief This method writes the complete state of all tracer
+     *        to the hard disk.
+     */
+    template <class Restarter>
+    void serialize(Restarter& res OPM_UNUSED)
+    {
+        /* not implemented */
+    }
+
+    /*!
+     * \brief This method restores the complete state of the tracer
+     *        from disk.
+     *
+     * It is the inverse of the serialize() method.
+     */
+    template <class Restarter>
+    void deserialize(Restarter& res OPM_UNUSED)
+    {
+        /* 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);
+
+        const auto& intQuants = elemCtx.intensiveQuantities(scvIdx, timeIdx);
+        const auto& fs = intQuants.fluidState();
+        Scalar phaseVolume = Opm::decay<Scalar>(fs.saturation(tracerPhaseIdx_[tracerIdx]))
+                *Opm::decay<Scalar>(fs.invB(tracerPhaseIdx_[tracerIdx]))
+                *Opm::decay<Scalar>(intQuants.porosity());
+
+        // avoid singular matrix if no water is present.
+        phaseVolume = Opm::max(phaseVolume, 1e-10);
+
+        if (std::is_same<LhsEval, Scalar>::value)
+            tracerStorage = phaseVolume * tracerConcentration0_[tracerIdx][globalDofIdx];
+        else
+            tracerStorage =
+                    phaseVolume
+                    * Opm::variable<LhsEval>(tracerConcentration_[tracerIdx][globalDofIdx][0], 0);
+    }
+
+    // evaluate the tracerflux over one face
+    void computeFlux_(TracerEvaluation & tracerFlux,
+                      const ElementContext& elemCtx,
+                      unsigned scvfIdx,
+                      unsigned timeIdx,
+                      const int tracerIdx)
+
+    {
+        const auto& stencil = elemCtx.stencil(timeIdx);
+        const auto& scvf = stencil.interiorFace(scvfIdx);
+
+        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+        unsigned inIdx = extQuants.interiorIndex();
+
+        const int tracerPhaseIdx = 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 = Opm::decay<Scalar>(extQuants.volumeFlux(tracerPhaseIdx));
+        Scalar b = Opm::decay<Scalar>(fs.invB(tracerPhaseIdx_[tracerIdx]));
+        Scalar c = tracerConcentration_[tracerIdx][globalUpIdx];
+
+        if (inIdx == upIdx)
+            tracerFlux = A*v*b*Opm::variable<TracerEvaluation>(c, 0);
+        else
+            tracerFlux = A*v*b*c;
+
+    }
+
+    bool linearSolve_(const TracerMatrix& M, TracerVector& x, TracerVector& b)
+    {
+#if ! DUNE_VERSION_NEWER(DUNE_COMMON, 2,7)
+        Dune::FMatrixPrecision<LinearSolverScalar>::set_singular_limit(1.e-30);
+        Dune::FMatrixPrecision<LinearSolverScalar>::set_absolute_limit(1.e-30);
+#endif
+        x = 0.0;
+        Scalar tolerance = 1e-2;
+        int maxIter = 100;
+
+        int verbosity = 0;
+        typedef Dune::BiCGSTABSolver<TracerVector> TracerSolver;
+        typedef Dune::MatrixAdapter<TracerMatrix, TracerVector , TracerVector > TracerOperator;
+        typedef Dune::SeqScalarProduct< TracerVector > TracerScalarProduct ;
+        typedef Dune::SeqILU< TracerMatrix, TracerVector, TracerVector  > TracerPreconditioner;
+
+        TracerOperator tracerOperator(M);
+        TracerScalarProduct tracerScalarProduct;
+        TracerPreconditioner tracerPreconditioner(M, 1);
+
+        TracerSolver solver (tracerOperator, tracerScalarProduct,
+                             tracerPreconditioner, tolerance, maxIter,
+                             verbosity);
+
+        Dune::InverseOperatorResult result;
+        solver.apply(x, b, result);
+
+        // return the result of the solver
+        return result.converged;
+    }
+
+    void linearize_(int tracerIdx)
+    {
+        (*tracerMatrix_) = 0.0;
+        tracerResidual_ = 0.0;
+
+        size_t numAllDof =  simulator_.model().numTotalDof();
+        std::vector<double> volumes(numAllDof, 0.0);
+        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);
+
+            Scalar extrusionFactor =
+                    elemCtx.intensiveQuantities(/*dofIdx=*/ 0, /*timeIdx=*/0).extrusionFactor();
+            Opm::Valgrind::CheckDefined(extrusionFactor);
+            assert(Opm::isfinite(extrusionFactor));
+            assert(extrusionFactor > 0.0);
+            Scalar scvVolume =
+                    elemCtx.stencil(/*timeIdx=*/0).subControlVolume(/*dofIdx=*/ 0).volume()
+                    * extrusionFactor;
+            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 = storageOfTimeIndex1_[tracerIdx][I];
+            else
+                computeStorage_(storageOfTimeIndex1, elemCtx, 0, /*timIdx=*/1, tracerIdx);
+
+            localStorage = (storageOfTimeIndex0 - storageOfTimeIndex1) * scvVolume/dt;
+            tracerResidual_[I][0] += localStorage.value(); //residual + flux
+            (*tracerMatrix_)[I][I][0][0] = localStorage.derivative(0);
+            size_t numInteriorFaces = elemCtx.numInteriorFaces(/*timIdx=*/0);
+            for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
+                TracerEvaluation flux;
+                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);
+                tracerResidual_[I][0] += flux.value(); //residual + flux
+                (*tracerMatrix_)[J][I][0][0] = -flux.derivative(0);
+                (*tracerMatrix_)[I][J][0][0] = flux.derivative(0);
+            }
+
+        }
+        // Wells
+        const int episodeIdx = simulator_.episodeIndex();
+        const auto& wells = simulator_.vanguard().schedule().getWells(episodeIdx);
+        for (auto well : wells) {
+
+            if (well->getStatus(episodeIdx) == Opm::WellCommon::SHUT)
+                continue;
+
+            const double wtracer = well->getTracerProperties(episodeIdx).getConcentration(tracerNames_[tracerIdx]);
+            std::array<int, 3> cartesianCoordinate;
+            for( auto& connection : well->getConnections(episodeIdx)) {
+
+                if (connection.state() == Opm::WellCompletion::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 = cartToGlobal_[cartIdx];
+                Scalar rate = simulator_.problem().wellModel().well(well->name())->volumetricSurfaceRateForConnection(I, tracerPhaseIdx_[tracerIdx]);
+                if (rate > 0)
+                    tracerResidual_[I][0] -= rate*wtracer;
+                else if (rate < 0)
+                    tracerResidual_[I][0] -= 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>>> tracerConcentration0_;
+    TracerMatrix *tracerMatrix_;
+    TracerVector tracerResidual_;
+    std::vector<Scalar> wTracer_;
+    std::vector<int> cartToGlobal_;
+    std::vector<Dune::BlockVector<Dune::FieldVector<Scalar, 1>>> storageOfTimeIndex1_;
+
+};
+} // namespace Ewoms
+
+#endif

ebos/tracervdtable.hh

@@ -0,0 +1,54 @@
+/*
+  Copyright (C) 2018 by NORCE
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#ifndef EWOMS_TRACER_VD_TABLE_HH
+#define	EWOMS_TRACER_VD_TABLE_HH
+#include <opm/parser/eclipse/EclipseState/Tables/SimpleTable.hpp>
+
+namespace Ewoms {
+
+
+    /*!
+     * \brief A class that contains tracer concentration vs depth table
+     */
+    class TracerVdTable : public Opm::SimpleTable {
+    public:
+        TracerVdTable( const Opm::DeckItem& item ) {
+            this->m_schema.addColumn( Opm::ColumnSchema( "DEPTH" , Opm::Table::STRICTLY_INCREASING , Opm::Table::DEFAULT_NONE ));
+            this->m_schema.addColumn( Opm::ColumnSchema( "TRACER_CONCENTRATION" , Opm::Table::RANDOM , Opm::Table::DEFAULT_NONE ));
+
+            Opm::SimpleTable::init(item);
+        }
+        /*!
+         * \brief Return the depth column
+         */
+        const Opm::TableColumn& getDepthColumn() const {
+            return Opm::SimpleTable::getColumn(0);
+        }
+        /*!
+         * \brief Return the tracer concentration column
+         */
+        const Opm::TableColumn& getTracerConcentration() const {
+            return Opm::SimpleTable::getColumn(1);
+        }
+    };
+}
+
+
+#endif // TRACERVDTABLE_HH
+

ebos/vtkecltracermodule.hh

@@ -0,0 +1,170 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ * \copydoc Ewoms::VtkEclTracerModule
+ */
+#ifndef EWOMS_VTK_ECL_TRACER_MODULE_HH
+#define EWOMS_VTK_ECL_TRACER_MODULE_HH
+
+
+#include <ewoms/io/vtkmultiwriter.hh>
+#include <ewoms/io/baseoutputmodule.hh>
+
+#include <ewoms/common/propertysystem.hh>
+#include <ewoms/common/parametersystem.hh>
+#include <ewoms/models/blackoil/blackoilproperties.hh>
+
+
+#include <dune/common/fvector.hh>
+
+#include <cstdio>
+
+BEGIN_PROPERTIES
+
+// create new type tag for the VTK tracer output
+NEW_TYPE_TAG(VtkEclTracer);
+
+// create the property tags needed for the tracer model
+NEW_PROP_TAG(EnableVtkOutput);
+NEW_PROP_TAG(VtkOutputFormat);
+NEW_PROP_TAG(VtkWriteEclTracerConcentration);
+
+// set default values for what quantities to output
+SET_BOOL_PROP(VtkEclTracer, VtkWriteEclTracerConcentration, false);
+END_PROPERTIES
+
+namespace Ewoms {
+    /*!
+ * \ingroup Vtk
+ *
+ * \brief VTK output module for the tracer model's parameters.
+ */
+    template <class TypeTag>
+    class VtkEclTracerModule : public BaseOutputModule<TypeTag>
+    {
+        typedef BaseOutputModule<TypeTag> ParentType;
+
+        typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+        typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+        typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
+        typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
+
+        typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+        typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+
+        static const int vtkFormat = GET_PROP_VALUE(TypeTag, VtkOutputFormat);
+        typedef Ewoms::VtkMultiWriter<GridView, vtkFormat> VtkMultiWriter;
+
+
+        typedef typename ParentType::ScalarBuffer ScalarBuffer;
+
+    public:
+        VtkEclTracerModule(const Simulator& simulator)
+            : ParentType(simulator)
+        { }
+
+        /*!
+     * \brief Register all run-time parameters for the tracer VTK output
+     * module.
+     */
+        static void registerParameters()
+        {
+            EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteEclTracerConcentration,
+                                 "Include the tracer concentration "
+                                 "in the VTK output files");
+        }
+
+        /*!
+     * \brief Allocate memory for the scalar fields we would like to
+     *        write to the VTK file.
+     */
+        void allocBuffers()
+        {
+            if (eclTracerConcentrationOutput_()){
+                const auto& tracerModel = this->simulator_.problem().tracerModel();
+                eclTracerConcentration_.resize(tracerModel.numTracers());
+                for(size_t tracerIdx=0; tracerIdx<eclTracerConcentration_.size();++tracerIdx){
+
+                    this->resizeScalarBuffer_(eclTracerConcentration_[tracerIdx]);
+                }
+            }
+
+        }
+
+        /*!
+     * \brief Modify the internal buffers according to the intensive quantities relevant for
+     *        an element
+     */
+        void processElement(const ElementContext& elemCtx)
+        {
+            if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput))
+                return;
+
+            const auto& tracerModel = elemCtx.problem().tracerModel();
+
+            for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
+                unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
+
+                if (eclTracerConcentrationOutput_()){
+                    for(size_t tracerIdx=0; tracerIdx<eclTracerConcentration_.size();++tracerIdx){
+                        eclTracerConcentration_[tracerIdx][globalDofIdx] = tracerModel.tracerConcentration(tracerIdx, globalDofIdx);
+                    }
+                }
+            }
+        }
+
+        /*!
+     * \brief Add all buffers to the VTK output writer.
+     */
+        void commitBuffers(BaseOutputWriter& baseWriter)
+        {
+            VtkMultiWriter *vtkWriter = dynamic_cast<VtkMultiWriter*>(&baseWriter);
+            if (!vtkWriter)
+                return;
+
+            if (eclTracerConcentrationOutput_()){
+                const auto& tracerModel = this->simulator_.problem().tracerModel();
+                for(size_t tracerIdx=0; tracerIdx<eclTracerConcentration_.size();++tracerIdx){
+                    const std::string tmp = "tracerConcentration_" + tracerModel.tracerName(tracerIdx);
+                    this->commitScalarBuffer_(baseWriter,tmp.c_str(), eclTracerConcentration_[tracerIdx]);
+                }
+            }
+
+
+
+        }
+
+    private:
+        static bool eclTracerConcentrationOutput_()
+        {
+            static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteEclTracerConcentration);
+            return val;
+        }
+
+
+        std::vector<ScalarBuffer> eclTracerConcentration_;
+    };
+} // namespace Ewoms
+
+#endif