Merge branch 'feature/2p-tracer-v2' into 'master'

Feature/2p tracer

Closes #498

See merge request !1526

dumux/porousmediumflow/tracer/localresidual.hh

@@ -85,7 +85,8 @@ public:
             for (int compIdx = 0; compIdx < numComponents; ++compIdx)
                 storage[compIdx] += volVars.porosity()
                                     * volVars.molarDensity(phaseIdx)
-                                    * volVars.moleFraction(phaseIdx, compIdx);
+                                    * volVars.moleFraction(phaseIdx, compIdx)
+                                    * volVars.saturation(phaseIdx);
         }
         // formulation with mass balances
         else
@@ -93,7 +94,9 @@ public:
             for (int compIdx = 0; compIdx < numComponents; ++compIdx)
                 storage[compIdx] += volVars.porosity()
                                     * volVars.density(phaseIdx)
-                                    * volVars.massFraction(phaseIdx, compIdx);
+                                    * volVars.massFraction(phaseIdx, compIdx)
+                                    * volVars.saturation(phaseIdx);
+
         }
 
         return storage;

dumux/porousmediumflow/tracer/volumevariables.hh

@@ -21,17 +21,30 @@
  * \ingroup TracerModel
  * \brief Quantities required by the tracer model in a control volume.
  */
-
 #ifndef DUMUX_TRACER_VOLUME_VARIABLES_HH
 #define DUMUX_TRACER_VOLUME_VARIABLES_HH
 
+#include <type_traits>
 #include <array>
+#include <dune/common/std/type_traits.hh>
 
 #include <dumux/porousmediumflow/volumevariables.hh>
 #include <dumux/material/solidstates/updatesolidvolumefractions.hh>
 
 namespace Dumux {
 
+namespace Detail {
+// helper structs and functions detecting if the user-defined spatial params class
+// has user-specified functions saturation() for multi-phase tracer.
+template <typename T, typename ...Ts>
+using SaturationDetector = decltype(std::declval<T>().spatialParams().saturation(std::declval<Ts>()...));
+
+template<class T, typename ...Args>
+static constexpr bool hasSaturation()
+{ return Dune::Std::is_detected<SaturationDetector, T, Args...>::value; }
+
+} // end namespace Detail
+
 /*!
  * \ingroup TracerModel
  * \brief Contains the quantities which are constant within a
@@ -74,6 +87,7 @@ public:
         // the spatial params special to the tracer model
         fluidDensity_ = problem.spatialParams().fluidDensity(element, scv);
         fluidMolarMass_ = problem.spatialParams().fluidMolarMass(element, scv);
+        fluidSaturation_ = saturation_(problem, element, scv);
 
         for (int compIdx = 0; compIdx < ParentType::numFluidComponents(); ++compIdx)
         {
@@ -87,9 +101,9 @@ public:
      *
      * We always forward to the fluid state with the phaseIdx property (see class description).
      *
-     * \param pIdx TODO docme!
+     * \param phaseIdx The phase index
      */
-    Scalar density(int pIdx = 0) const
+    Scalar density(int phaseIdx = 0) const
     { return fluidDensity_; }
 
     /*!
@@ -102,12 +116,14 @@ public:
      * \brief Returns the saturation.
      *
      * This method is here for compatibility reasons with other models. The saturation
-     * is always 1.0 in a one-phasic context.
-     *
-     * \param pIdx The phase index
+     * is always 1.0 in a one-phasic context, if two-phases or richards are considered,
+     * the spatialParams serve as way to pass the saturation from the main-file to the
+     * volVars and then to the localresidual for the tracer model.
+
+     * \param phaseIdx The phase index
      */
-    Scalar saturation(int pIdx = 0) const
-    { return 1.0; }
+    Scalar saturation(int phaseIdx = 0) const
+    { return fluidSaturation_ ; }
 
     /*!
      * \brief Returns the mobility.
@@ -115,53 +131,53 @@ public:
      * This method is here for compatibility reasons with other models. The mobility is always 1
      * for one-phasic models where the velocity field is given
      *
-     * \param pIdx The phase index
+     * \param phaseIdx The phase index
      */
-    Scalar mobility(int pIdx = 0) const
+    Scalar mobility(int phaseIdx = 0) const
     { return 1.0; }
 
     /*!
      * \brief Returns the molar density \f$\mathrm{[mol/m^3]}\f$ the of the fluid phase.
      *
-     * \param pIdx The phase index
+     * \param phaseIdx The phase index
      */
-    Scalar molarDensity(int pIdx = 0) const
+    Scalar molarDensity(int phaseIdx = 0) const
     { return fluidDensity_/fluidMolarMass_; }
 
     /*!
      * \brief Returns the mole fraction \f$\mathrm{[mol/mol]}\f$ of a component in the phase.
      *
-     * \param pIdx The phase index
+     * \param phaseIdx The phase index
      * \param compIdx The index of the component
      */
-    Scalar moleFraction(int pIdx, int compIdx) const
+    Scalar moleFraction(int phaseIdx, int compIdx) const
     { return useMoles ? moleOrMassFraction_[compIdx] : moleOrMassFraction_[compIdx]/FluidSystem::molarMass(compIdx)*fluidMolarMass_; }
 
     /*!
      * \brief Returns the mass fraction \f$\mathrm{[kg/kg]}\f$ of a component in the phase.
      *
-     * \param pIdx The phase index
+     * \param phaseIdx The phase index
      * \param compIdx The index of the component
      */
-    Scalar massFraction(int pIdx, int compIdx) const
+    Scalar massFraction(int phaseIdx, int compIdx) const
     { return useMoles ? moleOrMassFraction_[compIdx]*FluidSystem::molarMass(compIdx)/fluidMolarMass_ : moleOrMassFraction_[compIdx]; }
 
     /*!
      * \brief Returns the concentration \f$\mathrm{[mol/m^3]}\f$  of a component in the phase.
      *
-     * \param pIdx The phase index
+     * \param phaseIdx The phase index
      * \param compIdx The index of the component
      */
-    Scalar molarity(int pIdx, int compIdx) const
-    { return moleFraction(pIdx, compIdx)*molarDensity(); }
+    Scalar molarity(int phaseIdx, int compIdx) const
+    { return moleFraction(phaseIdx, compIdx)*molarDensity(); }
 
     /*!
      * \brief Returns the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ in the fluid.
      *
-     * \param pIdx The phase index
+     * \param phaseIdx The phase index
      * \param compIdx The index of the component
      */
-    Scalar diffusionCoefficient(int pIdx, int compIdx) const
+    Scalar diffusionCoefficient(int phaseIdx, int compIdx) const
     { return diffCoeff_[compIdx]; }
 
     // /*!
@@ -180,6 +196,37 @@ public:
 protected:
     SolidState solidState_;
     Scalar fluidDensity_, fluidMolarMass_;
+    Scalar fluidSaturation_ = 1.0;
+    /*!
+     * \brief Gets the saturation in an scv.
+     *
+     * \param problem the problem to solve
+     * \param element the element (codim-0-entity) the scv belongs to
+     * \param scv the sub control volume
+     * \note this gets selected if the user uses the multiphase tracer
+     */
+     template<class Problem, class Element, class Scv,
+              std::enable_if_t<Detail::hasSaturation<Problem, Element, Scv>(), int> = 0>
+     Scalar saturation_(const Problem& problem,
+                        const Element& element,
+                        const Scv& scv)
+     { return problem.spatialParams().saturation(element, scv); }
+
+    /*!
+     * \brief Gets the saturation in an scv.
+     *
+     * \param problem the problem to solve
+     * \param element the element (codim-0-entity) the scv belongs to
+     * \param scv the sub control volume
+     * \note this gets selected if the user a single phase tracer
+     */
+    template<class Problem, class Element, class Scv,
+             std::enable_if_t<!Detail::hasSaturation<Problem, Element, Scv>(), int> = 0>
+    Scalar saturation_(const Problem& problem,
+                       const Element &element,
+                       const Scv &scv)
+    { return 1.0; }
+
     // DispersivityType dispersivity_;
     std::array<Scalar, ParentType::numFluidComponents()> diffCoeff_;
     std::array<Scalar, ParentType::numFluidComponents()> moleOrMassFraction_;

test/porousmediumflow/2p/implicit/incompressible/spatialparams.hh

@@ -78,8 +78,8 @@ public:
         outerMaterialParams_.setVgAlpha(0.0037);
         outerMaterialParams_.setVgn(4.7);
 
-        lensK_ = 9.05e-12;
-        outerK_ = 4.6e-10;
+        lensK_ = getParam<Scalar>("SpatialParams.lensK", 9.05e-12);
+        outerK_ = getParam<Scalar>("SpatialParams.outerK", 4.6e-10);
     }
 
     /*!

test/porousmediumflow/tracer/2ptracer/CMakeLists.txt

+dune_symlink_to_source_files(FILES "params.input")
+
+dune_add_test(NAME test_2ptracer_lens_tpfa
+              LABELS porousmediumflow tracer
+              SOURCES main.cc
+              CMAKE_GUARD HAVE_UMFPACK
+              COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+              CMD_ARGS --script fuzzy
+                       --files ${CMAKE_SOURCE_DIR}/test/references/test_2ptracer_lens_tpfa_transport-reference.vtu
+                               ${CMAKE_CURRENT_BINARY_DIR}/test_2ptracer_lens_tpfa-00050.vtu
+                               ${CMAKE_SOURCE_DIR}/test/references/test_2ptracer_lens_tpfa_2p-reference.vtu
+                               ${CMAKE_CURRENT_BINARY_DIR}/test_2ptracer_lens_tpfa_2p-00050.vtu
+                       --command "${CMAKE_CURRENT_BINARY_DIR}/test_2ptracer_lens_tpfa params.input -Problem.Name test_2ptracer_lens_tpfa")

test/porousmediumflow/tracer/2ptracer/main.cc

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program 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.                                     *
+ *                                                                           *
+ *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TracerTests
+ * \brief Test for the 2p tracer CC model
+ */
+#include <config.h>
+
+#include <ctime>
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/timer.hh>
+#include <dune/grid/io/file/dgfparser/dgfexception.hh>
+#include <dune/grid/io/file/vtk.hh>
+
+#include <test/porousmediumflow/2p/implicit/incompressible/problem.hh>
+#include "problem_tracer.hh"
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+
+#include <dumux/linear/amgbackend.hh>
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/nonlinear/newtonsolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.hh>
+
+#include <dumux/discretization/method.hh>
+
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager.hh>
+
+int main(int argc, char** argv) try
+{
+    using namespace Dumux;
+
+    //! define the type tags for this problem
+    using TwoPTypeTag = Properties::TTag::TwoPIncompressibleTpfa;
+    using TracerTypeTag = Properties::TTag::TwoPTracerTestTpfa;
+    //! initialize MPI, finalize is done automatically on exit
+    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+    //! print dumux start message
+    if (mpiHelper.rank() == 0)
+        DumuxMessage::print(/*firstCall=*/true);
+
+    // parse command line arguments and input file
+    Parameters::init(argc, argv);
+
+    // try to create a grid (from the given grid file or the input file)
+    GridManager<GetPropType<TwoPTypeTag, Properties::Grid>> gridManager;
+    gridManager.init();
+
+    // get some time loop parameters
+    using Scalar =  GetPropType<TwoPTypeTag, Properties::Scalar>;
+    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+    auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+
+    // check if we are about to restart a previously interrupted simulation
+    Scalar restartTime = getParam<Scalar>("Restart.Time", 0);
+
+    // instantiate time loop
+    auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(restartTime, dt, tEnd);
+    timeLoop->setMaxTimeStepSize(maxDt);
+
+    ////////////////////////////////////////////////////////////
+    // set 2p Problem
+    ////////////////////////////////////////////////////////////
+
+    // we compute on the leaf grid view
+    const auto& leafGridView = gridManager.grid().leafGridView();
+
+    // create the finite volume grid geometry
+    using FVGridGeometry = GetPropType<TwoPTypeTag, Properties::FVGridGeometry>;
+    auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+    fvGridGeometry->update();
+
+    // the problem (initial and boundary conditions)
+    using TwoPProblem = GetPropType<TwoPTypeTag, Properties::Problem>;
+    auto twoPProblem = std::make_shared<TwoPProblem>(fvGridGeometry);
+
+    // the solution vector
+    using TwoPSolutionVector = GetPropType<TwoPTypeTag, Properties::SolutionVector>;
+    TwoPSolutionVector p(fvGridGeometry->numDofs());
+    twoPProblem->applyInitialSolution(p);
+    auto pOld = p;
+
+    // maybe update the interface parameters
+    if (ENABLEINTERFACESOLVER)
+        twoPProblem->spatialParams().updateMaterialInterfaceParams(p);
+
+    // the grid variables
+    using TwoPGridVariables = GetPropType<TwoPTypeTag, Properties::GridVariables>;
+    auto twoPGridVariables = std::make_shared<TwoPGridVariables>(twoPProblem, fvGridGeometry);
+    twoPGridVariables->init(p);
+
+    // intialize the vtk output module
+    using TwoPIOFields = GetPropType<TwoPTypeTag, Properties::IOFields>;
+
+    // use non-conforming output for the test with interface solver
+    const auto ncOutput = getParam<bool>("Problem.UseNonConformingOutput", false);
+    VtkOutputModule<TwoPGridVariables, TwoPSolutionVector> twoPVtkWriter(*twoPGridVariables, p, twoPProblem->name()+"_2p", "",
+                                                                        (ncOutput ? Dune::VTK::nonconforming : Dune::VTK::conforming));
+
+    TwoPIOFields::initOutputModule(twoPVtkWriter); //!< Add model specific output fields
+    twoPVtkWriter.write(0.0);
+
+    // the assembler with time loop for instationary problem
+    using TwoPAssembler = FVAssembler<TwoPTypeTag, DiffMethod::numeric>;
+    auto twoPAssembler = std::make_shared<TwoPAssembler>(twoPProblem, fvGridGeometry, twoPGridVariables, timeLoop);
+
+    // the linear solver
+    using TwoPLinearSolver = AMGBackend<TwoPTypeTag>;
+    auto twoPLinearSolver = std::make_shared<TwoPLinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+    // the non-linear solver
+    using NewtonSolver = Dumux::NewtonSolver<TwoPAssembler, TwoPLinearSolver>;
+    NewtonSolver nonLinearSolver(twoPAssembler, twoPLinearSolver);
+
+    ////////////////////////////////////////////////////////////
+    // set tracer Problem
+    ////////////////////////////////////////////////////////////
+
+    //! the problem (initial and boundary conditions)
+    using TracerProblem = GetPropType<TracerTypeTag, Properties::Problem>;
+    auto tracerProblem = std::make_shared<TracerProblem>(fvGridGeometry);
+
+    //! the solution vector
+    using TracerSolutionVector = GetPropType<TracerTypeTag, Properties::SolutionVector>;
+    TracerSolutionVector x(leafGridView.size(0));
+    tracerProblem->applyInitialSolution(x);
+    auto xOld = x;
+
+    //! initialize the flux, density and saturation vectors
+    std::vector<Scalar> volumeFlux_(fvGridGeometry->numScvf(), 0.0);
+    std::vector<Scalar> density_(fvGridGeometry->numScv(), 0.0);
+    std::vector<Scalar> saturation_(fvGridGeometry->numScv(), 0.0);
+
+    //! the grid variables
+    using TracerGridVariables = GetPropType<TracerTypeTag, Properties::GridVariables>;
+    auto tracerGridVariables = std::make_shared<TracerGridVariables>(tracerProblem, fvGridGeometry);
+    tracerGridVariables->init(x);
+
+    // the linear solver
+    using TracerLinearSolver = AMGBackend<TracerTypeTag>;
+    auto tracerLinearSolver = std::make_shared<TracerLinearSolver>(leafGridView, fvGridGeometry->dofMapper());
+
+     //! the linear system
+    using JacobianMatrix = GetPropType<TracerTypeTag, Properties::JacobianMatrix>;
+    auto A = std::make_shared<JacobianMatrix>();
+    auto r = std::make_shared<TracerSolutionVector>();
+
+    //! the assembler with time loop for instationary problem
+    using TracerAssembler = FVAssembler<TracerTypeTag, DiffMethod::analytic, /*implicit=*/false>;
+    auto tracerAssembler = std::make_shared<TracerAssembler>(tracerProblem, fvGridGeometry, tracerGridVariables, timeLoop);
+    tracerAssembler->setLinearSystem(A, r);
+
+    // set the flux, density and saturation from the 2p problem
+    tracerProblem->spatialParams().setVolumeFlux(volumeFlux_);
+    tracerProblem->spatialParams().setDensity(density_);
+    tracerProblem->spatialParams().setSaturation(saturation_);
+
+    //! initialize the vtk output module
+    VtkOutputModule<TracerGridVariables, TracerSolutionVector> vtkWriter(*tracerGridVariables, x, tracerProblem->name());
+    using TracerIOFields = GetPropType<TracerTypeTag, Properties::IOFields>;
+    TracerIOFields::initOutputModule(vtkWriter); //!< Add model specific output fields
+    using VelocityOutput = GetPropType<TracerTypeTag, Properties::VelocityOutput>;
+    vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*tracerGridVariables));
+    vtkWriter.write(0.0);
+
+    //! set some check points for the time loop
+    timeLoop->setPeriodicCheckPoint(tEnd/50.0);
+
+    ////////////////////////////////////////////////////////////
+    // run instationary non-linear problem
+    ////////////////////////////////////////////////////////////
+
+    // time loop
+    timeLoop->start(); do
+    {
+        // set previous solution for storage evaluations
+        twoPAssembler->setPreviousSolution(pOld);
+
+        // solve the non-linear system with time step control
+        nonLinearSolver.solve(p, *timeLoop);
+
+        // make the new solution the old solution
+        pOld = p;
+        twoPGridVariables->advanceTimeStep();
+
+        // write vtk output
+        twoPVtkWriter.write(timeLoop->time());
+
+        // report statistics of this time step
+        timeLoop->reportTimeStep();
+
+        // set new dt as suggested by the Newton solver
+        timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+
+        // loop over elements to compute fluxes, saturations, densities for tracer
+        using FluxVariables =  typename GET_PROP_TYPE(TwoPTypeTag, FluxVariables);
+        auto upwindTerm = [](const auto& volVars) { return volVars.mobility(0); };
+        for (const auto& element : elements(leafGridView))
+        {
+            ////////////////////////////////////////////////////////////
+            // compute volume fluxes for the tracer model
+            ///////////////////////////////////////////////////////////
+
+            auto fvGeometry = localView(*fvGridGeometry);
+            fvGeometry.bind(element);
+
+            auto elemVolVars = localView(twoPGridVariables->curGridVolVars());
+            elemVolVars.bind(element, fvGeometry, p);
+
+            auto elemFluxVars = localView(twoPGridVariables->gridFluxVarsCache());
+            elemFluxVars.bind(element, fvGeometry, elemVolVars);
+
+            for (const auto& scvf : scvfs(fvGeometry))
+            {
+                const auto idx = scvf.index();
+
+                if (!scvf.boundary())
+                {
+                    FluxVariables fluxVars;
+                    fluxVars.init(*twoPProblem, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
+                            volumeFlux_[idx] = fluxVars.advectiveFlux(0, upwindTerm);
+                }
+                else
+                {
+                    const auto bcTypes = twoPProblem->boundaryTypes(element, scvf);
+                    if (bcTypes.hasOnlyDirichlet())
+                    {
+                        FluxVariables fluxVars;
+                        fluxVars.init(*twoPProblem, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
+                                volumeFlux_[idx] = fluxVars.advectiveFlux(0, upwindTerm);
+                    }
+                }
+            }
+
+            ////////////////////////////////////////////////////////////
+            // compute densities and saturations for the tracer model
+            ///////////////////////////////////////////////////////////
+            for (const auto& scv : scvs(fvGeometry))
+            {
+                const auto& volVars = elemVolVars[scv];
+                const auto idx = scv.dofIndex();
+
+                density_[idx] = volVars.density(0);
+                saturation_[idx] = volVars.saturation(0);
+            }
+        }
+
+        ////////////////////////////////////////////////////////////
+        // solve tracer problem on the same grid
+        ////////////////////////////////////////////////////////////
+
+        // set the flux from the 2p problem
+        tracerProblem->spatialParams().setVolumeFlux(volumeFlux_);
+        tracerProblem->spatialParams().setDensity(density_);
+        tracerProblem->spatialParams().setSaturation(saturation_);
+
+        // set previous solution for storage evaluations
+        tracerAssembler->setPreviousSolution(xOld);
+
+        Dune::Timer tracerAssembleTimer;
+        tracerAssembler->assembleJacobianAndResidual(x);
+        tracerAssembleTimer.stop();
+
+        // solve the linear system A(xOld-xNew) = r
+        Dune::Timer solveTimer;
+        TracerSolutionVector xDelta(x);
+        tracerLinearSolver->solve(*A, xDelta, *r);
+        solveTimer.stop();
+
+        // update solution and grid variables
+        Dune::Timer updateTimer;
+        updateTimer.reset();
+        x -= xDelta;
+        tracerGridVariables->update(x);
+        updateTimer.stop();
+
+        // statistics
+        Dune::Timer assembleTimer;
+        const auto elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
+        std::cout << "Assemble/solve/update time: "
+                  <<  assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
+                  <<  solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/"
+                  <<  updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
+                  <<  std::endl;
+
+        // make the new solution the old solution
+        xOld = x;
+        tracerGridVariables->advanceTimeStep();
+
+        // advance the time loop to the next step
+        timeLoop->advanceTimeStep();
+
+        // write vtk output
+        vtkWriter.write(timeLoop->time());
+
+    } while (!timeLoop->finished());
+
+    timeLoop->finalize(leafGridView.comm());
+
+    ////////////////////////////////////////////////////////////
+    // finalize, print dumux message to say goodbye
+    ////////////////////////////////////////////////////////////
+
+    //! print dumux end message
+    if (mpiHelper.rank() == 0)
+    {
+        Parameters::print();
+        DumuxMessage::print(/*firstCall=*/false);
+    }
+
+    return 0;
+} // end main
+
+catch (Dumux::ParameterException &e)
+{
+    std::cerr << std::endl << e << " ---> Abort!" << std::endl;
+    return 1;
+}
+catch (Dune::DGFException & e)
+{
+    std::cerr << "DGF exception thrown (" << e <<
+                 "). Most likely, the DGF file name is wrong "
+                 "or the DGF file is corrupted, "
+                 "e.g. missing hash at end of file or wrong number (dimensions) of entries."
+                 << " ---> Abort!" << std::endl;
+    return 2;
+}
+catch (Dune::Exception &e)
+{
+    std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
+    return 3;
+}
+catch (...)
+{
+    std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
+    return 4;
+}

test/porousmediumflow/tracer/2ptracer/params.input

+[TimeLoop]
+DtInitial = 10 # [s]
+TEnd = 450 # [s]
+MaxTimeStepSize = 10
+
+[Grid]
+UpperRight = 6 4
+Cells = 48 32
+
+[SpatialParams]
+LensLowerLeft = 1.0 2.0 # [m] coordinates of the lower left lens corner
+LensUpperRight = 4.0 3.0 # [m] coordinates of the upper right lens corner
+outerK = 1e-09 # [m^2]
+lensK = 1e-12 # [m^2]
+
+[Problem]
+Name = test_2ptracer_lens_tpfa
+EnableGravity = true
+BinaryDiffusionCoefficient = 0.0 # only advective transport
+
+[Newton]
+EnablePartialReassembly = true
+
+[Vtk]
+AddVelocity = true

test/porousmediumflow/tracer/2ptracer/problem_tracer.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program 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.                                     *
+ *                                                                           *
+ *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TracerTests
+ * \brief  A 2p problem with multiple tracer bands in a porous groundwater reservoir with a lens
+ */
+#ifndef DUMUX_TWOP_TRACER_TEST_PROBLEM_HH
+#define DUMUX_TWOP_TRACER_TEST_PROBLEM_HH
+
+#include <dumux/discretization/cctpfa.hh>
+#include <dumux/porousmediumflow/tracer/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/material/fluidsystems/base.hh>
+
+#include "spatialparams_tracer.hh"
+
+namespace Dumux {
+/*!
+ * \ingroup TracerTests
+ * \brief A 2p problem with multiple tracer bands in a porous groundwater reservoir with a lens
+ */
+template <class TypeTag>
+class TwoPTracerTestProblem;
+
+namespace Properties {
+//Create new type tags
+namespace TTag {
+struct TwoPTracerTest { using InheritsFrom = std::tuple<Tracer>; };
+struct TwoPTracerTestTpfa { using InheritsFrom = std::tuple<TwoPTracerTest, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::TwoPTracerTest> { using type = Dune::YaspGrid<2>; };
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::TwoPTracerTest> { using type = TwoPTracerTestProblem<TypeTag>; };
+
+// Set the spatial parameters
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::TwoPTracerTest>
+{
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using type = TwoPTracerTestSpatialParams<FVGridGeometry, Scalar>;
+};
+
+// Define whether mole(true) or mass (false) fractions are used
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::TwoPTracerTest> { static constexpr bool value = false; };
+template<class TypeTag>
+struct SolutionDependentMolecularDiffusion<TypeTag, TTag::TwoPTracerTestTpfa> { static constexpr bool value = false; };
+
+//! A simple fluid system with one tracer component
+template<class TypeTag>
+class TracerFluidSystem : public FluidSystems::Base<GetPropType<TypeTag, Properties::Scalar>,
+                                                                TracerFluidSystem<TypeTag>>
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+
+public:
+    //! If the fluid system only contains tracer components
+    static constexpr bool isTracerFluidSystem()
+    { return true; }
+
+    //! No component is the main component
+    static constexpr int getMainComponent(int phaseIdx)
+    { return -1; }
+
+    //! The number of components
+    static constexpr int numComponents = 1;
+
+    //! Human readable component name (index compIdx) (for vtk output)
+    static std::string componentName(int compIdx)
+    { return "tracer_" + std::to_string(compIdx); }
+
+    //! Molar mass in kg/mol of the component with index compIdx
+    static Scalar molarMass(unsigned int compIdx)
+    { return 0.300; }
+
+    //! Binary diffusion coefficient
+    //! (might depend on spatial parameters like pressure / temperature)
+    static Scalar binaryDiffusionCoefficient(unsigned int compIdx,
+                                             const Problem& problem,
+                                             const Element& element,
+                                             const SubControlVolume& scv)
+    {
+        static const Scalar D = getParam<Scalar>("Problem.BinaryDiffusionCoefficient");
+        return D;
+    }
+};
+
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::TwoPTracerTest> { using type = TracerFluidSystem<TypeTag>; };
+
+} // end namespace Properties
+
+/*!
+ * \ingroup TracerTests
+ *
+ * \brief Definition of a problem, for the tracer problem:
+ * A lens of contaminant tracer is diluted by diffusion and a base groundwater flow
+ *
+ * This problem uses the \ref TracerModel model.
+ *
+ * To run the simulation execute the following line in shell:
+ * <tt>./test_2ptracer -ParameterFile ./params.input</tt> or
+ * <tt>./test_2ptracer -ParameterFile ./params.input</tt>
+ */
+template <class TypeTag>
+class TwoPTracerTestProblem : public PorousMediumFlowProblem<TypeTag>
+{
+    using ParentType = PorousMediumFlowProblem<TypeTag>;
+
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using SpatialParams = GetPropType<TypeTag, Properties::SpatialParams>;
+
+    //! property that defines whether mole or mass fractions are used
+    static constexpr bool useMoles = getPropValue<TypeTag, Properties::UseMoles>();
+
+    using Element = typename FVGridGeometry::GridView::template Codim<0>::Entity;
+    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+    TwoPTracerTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry)
+    {
+        // stating in the console whether mole or mass fractions are used
+        if(useMoles)
+            std::cout<<"problem uses mole fractions" << '\n';
+        else
+            std::cout<<"problem uses mass fractions" << '\n';
+
+        stripeWidth_ = getParam<Scalar>("Problem.StripeWidth", 0.125);
+    }
+
+    /*!
+     * \name Boundary conditions
+     */
+    // \{
+
+    /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary segment.
+     *
+     * \param globalPos The position for which the bc type should be evaluated
+     */
+    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+    {
+        BoundaryTypes values;
+        if (onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
+            values.setAllDirichlet();
+        else
+            values.setAllNeumann();
+        return values;
+    }
+    // \}
+
+    /*!
+     * \name Volume terms
+     */
+    // \{
+
+    /*!
+     * \brief Evaluates the boundary conditions for a Dirichlet boundary segment.
+     *
+     * \param globalPos The global position
+     */
+    PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+    {
+        PrimaryVariables initialValues(0.0);
+        return initialValues;
+    }
+
+    /*!
+     * \brief Evaluates the initial value for a control volume.
+     *
+     * \param globalPos The position for which the initial condition should be evaluated
+     *
+     * For this method, the \a values parameter stores primary
+     * variables.
+     */
+    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+    {
+        PrimaryVariables initialValues(0.0);
+
+        if (onStripe1_(globalPos) || onStripe2_(globalPos) || onStripe3_(globalPos))
+        {
+            if (useMoles)
+                initialValues = 1e-9;
+            else
+                initialValues = 1e-9*FluidSystem::molarMass(0)/this->spatialParams().fluidMolarMass(globalPos);
+        }
+        return initialValues;
+    }
+
+private:
+    static constexpr Scalar eps_ = 1e-6;
+    Scalar stripeWidth_;
+
+    bool onUpperBoundary_(const GlobalPosition &globalPos) const
+    {
+        return globalPos[1] > this->fvGridGeometry().bBoxMax()[1] - 0.1 - eps_;
+    }
+
+    bool onLeftBoundary_(const GlobalPosition &globalPos) const
+    {
+        return globalPos[0] < this->fvGridGeometry().bBoxMin()[0] + eps_;
+    }
+
+    bool onRightBoundary_(const GlobalPosition &globalPos) const
+    {
+        return globalPos[0] > this->fvGridGeometry().bBoxMax()[0] - eps_;
+    }
+
+    bool onStripe1_(const GlobalPosition &globalPos) const
+    {
+       const auto xMax = this->fvGridGeometry().bBoxMax()[0];
+       return  (
+           ( (xMax /4.0 - stripeWidth_*0.5) < globalPos[0] + eps_ ) &&
+           ( (xMax/4.0 + stripeWidth_*0.5) > globalPos[0] + eps_  )
+       );
+    }
+
+    bool onStripe2_(const GlobalPosition &globalPos) const
+    {
+        const auto xMax = this->fvGridGeometry().bBoxMax()[0];
+        return  (
+            ( (2.0 * xMax /4.0 - stripeWidth_*0.5) < globalPos[0] + eps_ ) &&
+            ( (2.0 * xMax/4.0 + stripeWidth_*0.5) > globalPos[0] + eps_ )
+        );
+    }
+
+    bool onStripe3_(const GlobalPosition &globalPos) const
+    {
+        const auto xMax = this->fvGridGeometry().bBoxMax()[0];
+        return  (
+            ( (3.0 * xMax /4.0 - stripeWidth_*0.5) < globalPos[0] + eps_ ) &&
+            ( (3.0 * xMax/4.0 + stripeWidth_*0.5) > globalPos[0] + eps_ )
+        );
+    }
+
+};
+
+} //end namespace Dumux
+
+#endif

test/porousmediumflow/tracer/2ptracer/spatialparams_tracer.hh

+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program 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.                                     *
+ *                                                                           *
+ *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup TracerTests
+ * \brief Definition of the spatial parameters for the tracer problem
+ */
+#ifndef DUMUX_TWOP_TRACER_TEST_SPATIAL_PARAMS_HH
+#define DUMUX_TWOP_TRACER_TEST_SPATIAL_PARAMS_HH
+
+#include <dumux/porousmediumflow/properties.hh>
+#include <dumux/material/spatialparams/fv.hh>
+#include <iostream>
+
+namespace Dumux {
+
+/*!
+ * \ingroup TracerTests
+ * \brief Definition of the spatial parameters for the tracer problem
+ */
+template<class FVGridGeometry, class Scalar>
+class TwoPTracerTestSpatialParams
+: public FVSpatialParams<FVGridGeometry, Scalar,
+                         TwoPTracerTestSpatialParams<FVGridGeometry, Scalar>>
+{
+    using GridView = typename FVGridGeometry::GridView;
+    using FVElementGeometry = typename FVGridGeometry::LocalView;
+    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+    using Element = typename GridView::template Codim<0>::Entity;
+    using ParentType = FVSpatialParams<FVGridGeometry, Scalar,
+                                       TwoPTracerTestSpatialParams<FVGridGeometry, Scalar>>;
+
+    static const int dimWorld = GridView::dimensionworld;
+    using GlobalPosition = typename Dune::FieldVector<Scalar, dimWorld>;
+
+
+public:
+
+    TwoPTracerTestSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+    : ParentType(fvGridGeometry) {}
+
+    /*!
+     * \brief Defines the porosity \f$\mathrm{[-]}\f$.
+     *
+     * \param globalPos The global position
+     */
+    Scalar porosityAtPos(const GlobalPosition& globalPos) const
+    { return 0.4; }
+
+    //! Fluid properties that are spatial parameters in the tracer model
+    //! They can possible vary with space but are usually constants
+
+    //! Fluid density
+    Scalar fluidDensity(const Element &element,
+                        const SubControlVolume& scv) const
+    { return density_[this->fvGridGeometry().elementMapper().index(element)];; }
+
+    void setDensity(const std::vector<Scalar>& d)
+    { density_ = d; }
+
+    //! fluid molar mass
+    Scalar fluidMolarMass(const Element &element,
+                          const SubControlVolume& scv) const
+    { return 0.018; }
+
+    Scalar fluidMolarMass(const GlobalPosition &globalPos) const
+    { return 0.018; }
+
+    //! Velocity field
+    template<class ElementVolumeVariables>
+    Scalar volumeFlux(const Element &element,
+                      const FVElementGeometry& fvGeometry,
+                      const ElementVolumeVariables& elemVolVars,
+                      const SubControlVolumeFace& scvf) const
+    { return volumeFlux_[scvf.index()]; }
+
+    void setVolumeFlux(const std::vector<Scalar>& f)
+    { volumeFlux_ = f; }
+
+    //! saturation from twoPProblem
+    Scalar saturation(const Element &element,
+                      const SubControlVolume& scv) const
+    { return saturation_[scv.dofIndex()]; }
+
+    void setSaturation(const std::vector<Scalar>& s)
+    { saturation_ = s; }
+
+private:
+    std::vector<Scalar> volumeFlux_;
+    std::vector<Scalar> density_;
+    std::vector<Scalar> saturation_;
+};
+
+} // end namespace Dumux
+
+#endif

test/porousmediumflow/tracer/CMakeLists.txt

-add_subdirectory(constvel)
 add_subdirectory(1ptracer)
+add_subdirectory(2ptracer)
+add_subdirectory(constvel)
 add_subdirectory(multicomp)