From 1ae9db04d6a3b82af5ef0ff838200f821511bba0 Mon Sep 17 00:00:00 2001
From: Felix Weinhardt <felix.weinhardt@gmx.de>
Date: Wed, 24 Apr 2019 12:19:37 +0200
Subject: [PATCH] added 1ptracer to examples
---
examples/1ptracer/CMakeLists.txt | 13 +
examples/1ptracer/README.md | 28 ++
examples/1ptracer/main.cc | 328 ++++++++++++++++++++++
examples/1ptracer/params.input | 21 ++
examples/1ptracer/problem_1p.hh | 165 +++++++++++
examples/1ptracer/problem_tracer.hh | 231 +++++++++++++++
examples/1ptracer/spatialparams_1p.hh | 131 +++++++++
examples/1ptracer/spatialparams_tracer.hh | 114 ++++++++
examples/CMakeLists.txt | 1 +
9 files changed, 1032 insertions(+)
create mode 100644 examples/1ptracer/CMakeLists.txt
create mode 100644 examples/1ptracer/README.md
create mode 100644 examples/1ptracer/main.cc
create mode 100644 examples/1ptracer/params.input
create mode 100644 examples/1ptracer/problem_1p.hh
create mode 100644 examples/1ptracer/problem_tracer.hh
create mode 100644 examples/1ptracer/spatialparams_1p.hh
create mode 100644 examples/1ptracer/spatialparams_tracer.hh
diff --git a/examples/1ptracer/CMakeLists.txt b/examples/1ptracer/CMakeLists.txt
new file mode 100644
index 0000000000..845a699be3
--- /dev/null
+++ b/examples/1ptracer/CMakeLists.txt
@@ -0,0 +1,13 @@
+dune_symlink_to_source_files(FILES "params.input")
+
+dune_add_test(NAME example_1ptracer
+ 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_1ptracer_transport-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/test_1ptracer-00010.vtu
+ ${CMAKE_SOURCE_DIR}/test/references/test_1ptracer_pressure-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/1p.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/example_1ptracer params.input -Problem.Name example_1ptracer")
diff --git a/examples/1ptracer/README.md b/examples/1ptracer/README.md
new file mode 100644
index 0000000000..82211db997
--- /dev/null
+++ b/examples/1ptracer/README.md
@@ -0,0 +1,28 @@
+This tutorial was copied from /home/felix/Dokumente/Aktuelles_U/19_03_27_DumuxTag/dumux/dumux/test/porousmediumflow/tracer/1ptracer.
+
+# One-phase flow with with random permeability distribution (and a tracer model)
+
+## Problem set-up
+This example combines a stationary One-phase flow in a porous medium with randomly generated permeabilities with a tracer model. In the first step, the groundwater-velocity is evaluated under stationary conditions. Based on the velocity field the tracer model is running instationary.
+
+## Random permeability generation
+### spatialparams_1p.hh
+The follwing code can be found in lines 64-72.
+
+```C++
+// generate random fields
+ std::mt19937 rand(0);
+ std::lognormal_distribution<Scalar> K(std::log(permeability_), std::log(permeability_)*0.1);
+ std::lognormal_distribution<Scalar> KLens(std::log(permeabilityLens_), std::log(permeabilityLens_)*0.1);
+ for (const auto& element : elements(fvGridGeometry->gridView()))
+ {
+ const auto eIdx = fvGridGeometry->elementMapper().index(element);
+ const auto globalPos = element.geometry().center();
+ K_[eIdx] = isInLens_(globalPos) ? KLens(rand) : K(rand);
+ }
+```
+Two lognormal distributions are defined: K and KLens with different means and standard deviations. In this case the mean is defined as the log of the permeability given from the input file and the standard deviation is 10% of the mean.
+Within a loop through all elements the randomly generated permeabilities are assigned according to their position in the domain (inside or outside the lense).
+
+## Tracer model
+...
\ No newline at end of file
diff --git a/examples/1ptracer/main.cc b/examples/1ptracer/main.cc
new file mode 100644
index 0000000000..caae7310a4
--- /dev/null
+++ b/examples/1ptracer/main.cc
@@ -0,0 +1,328 @@
+// -*- 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 Test for the tracer CC model.
+ */
+
+#include <config.h>
+
+#include "problem_1p.hh"
+#include "problem_tracer.hh"
+
+#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 <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/dumuxmessage.hh>
+
+#include <dumux/linear/seqsolverbackend.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+#include <dumux/assembly/diffmethod.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 OnePTypeTag = Properties::TTag::IncompressibleTest;
+ using TracerTypeTag = Properties::TTag::TracerTestCC;
+
+ //! 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 the command line arguments and input file
+ ////////////////////////////////////////////////////////////
+
+ //! parse command line arguments
+ Parameters::init(argc, argv);
+
+ //////////////////////////////////////////////////////////////////////
+ // try to create a grid (from the given grid file or the input file)
+ /////////////////////////////////////////////////////////////////////
+
+ // only create the grid once using the 1p type tag
+ GridManager<GetPropType<OnePTypeTag, Properties::Grid>> gridManager;
+ gridManager.init();
+
+ //! we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ ////////////////////////////////////////////////////////////
+ // setup & solve 1p problem on this grid
+ ////////////////////////////////////////////////////////////
+
+ //! create the finite volume grid geometry
+ using FVGridGeometry = GetPropType<OnePTypeTag, Properties::FVGridGeometry>;
+ auto fvGridGeometry = std::make_shared<FVGridGeometry>(leafGridView);
+ fvGridGeometry->update();
+
+ //! the problem (boundary conditions)
+ using OnePProblem = GetPropType<OnePTypeTag, Properties::Problem>;
+ auto problemOneP = std::make_shared<OnePProblem>(fvGridGeometry);
+
+ //! the solution vector
+ using JacobianMatrix = GetPropType<OnePTypeTag, Properties::JacobianMatrix>;
+ using SolutionVector = GetPropType<OnePTypeTag, Properties::SolutionVector>;
+ SolutionVector p(leafGridView.size(0));
+
+ //! the linear system
+ auto A = std::make_shared<JacobianMatrix>();
+ auto r = std::make_shared<SolutionVector>();
+
+ //! the grid variables
+ using OnePGridVariables = GetPropType<OnePTypeTag, Properties::GridVariables>;
+ auto onePGridVariables = std::make_shared<OnePGridVariables>(problemOneP, fvGridGeometry);
+ onePGridVariables->init(p);
+
+ //! the assembler
+ using OnePAssembler = FVAssembler<OnePTypeTag, DiffMethod::analytic>;
+ auto assemblerOneP = std::make_shared<OnePAssembler>(problemOneP, fvGridGeometry, onePGridVariables);
+ assemblerOneP->setLinearSystem(A, r);
+
+ Dune::Timer timer;
+ // assemble the local jacobian and the residual
+ Dune::Timer assemblyTimer; std::cout << "Assembling linear system ..." << std::flush;
+ assemblerOneP->assembleJacobianAndResidual(p);
+ assemblyTimer.stop(); std::cout << " took " << assemblyTimer.elapsed() << " seconds." << std::endl;
+
+ // we solve Ax = -r
+ (*r) *= -1.0;
+
+ //! solve the 1p problem
+ using LinearSolver = UMFPackBackend;
+ Dune::Timer solverTimer; std::cout << "Solving linear system ..." << std::flush;
+ auto linearSolver = std::make_shared<LinearSolver>();
+ linearSolver->solve(*A, p, *r);
+ solverTimer.stop(); std::cout << " took " << solverTimer.elapsed() << " seconds." << std::endl;
+
+ //! update the grid variables
+ Dune::Timer updateTimer; std::cout << "Updating variables ..." << std::flush;
+ onePGridVariables->update(p);
+ updateTimer.elapsed(); std::cout << " took " << updateTimer.elapsed() << std::endl;
+
+ //! write output to vtk
+ using GridView = GetPropType<OnePTypeTag, Properties::GridView>;
+ Dune::VTKWriter<GridView> onepWriter(leafGridView);
+ onepWriter.addCellData(p, "p");
+ const auto& k = problemOneP->spatialParams().getKField();
+ onepWriter.addCellData(k, "permeability");
+ onepWriter.write("1p");
+
+ timer.stop();
+
+ const auto& comm = Dune::MPIHelper::getCollectiveCommunication();
+ std::cout << "Simulation took " << timer.elapsed() << " seconds on "
+ << comm.size() << " processes.\n"
+ << "The cumulative CPU time was " << timer.elapsed()*comm.size() << " seconds.\n";
+
+ ////////////////////////////////////////////////////////////
+ // compute volume fluxes for the tracer model
+ ////////////////////////////////////////////////////////////
+ using Scalar = GetPropType<OnePTypeTag, Properties::Scalar>;
+ std::vector<Scalar> volumeFlux(fvGridGeometry->numScvf(), 0.0);
+
+ using FluxVariables = GetPropType<OnePTypeTag, Properties::FluxVariables>;
+ auto upwindTerm = [](const auto& volVars) { return volVars.mobility(0); };
+ for (const auto& element : elements(leafGridView))
+ {
+ auto fvGeometry = localView(*fvGridGeometry);
+ fvGeometry.bind(element);
+
+ auto elemVolVars = localView(onePGridVariables->curGridVolVars());
+ elemVolVars.bind(element, fvGeometry, p);
+
+ auto elemFluxVars = localView(onePGridVariables->gridFluxVarsCache());
+ elemFluxVars.bind(element, fvGeometry, elemVolVars);
+
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ const auto idx = scvf.index();
+
+ if (!scvf.boundary())
+ {
+ FluxVariables fluxVars;
+ fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
+ volumeFlux[idx] = fluxVars.advectiveFlux(0, upwindTerm);
+ }
+ else
+ {
+ const auto bcTypes = problemOneP->boundaryTypes(element, scvf);
+ if (bcTypes.hasOnlyDirichlet())
+ {
+ FluxVariables fluxVars;
+ fluxVars.init(*problemOneP, element, fvGeometry, elemVolVars, scvf, elemFluxVars);
+ volumeFlux[idx] = fluxVars.advectiveFlux(0, upwindTerm);
+ }
+ }
+ }
+ }
+
+ ////////////////////////////////////////////////////////////
+ // setup & solve tracer problem on the same grid
+ ////////////////////////////////////////////////////////////
+
+ //! the problem (initial and boundary conditions)
+ using TracerProblem = GetPropType<TracerTypeTag, Properties::Problem>;
+ auto tracerProblem = std::make_shared<TracerProblem>(fvGridGeometry);
+
+ // set the flux from the 1p problem
+ tracerProblem->spatialParams().setVolumeFlux(volumeFlux);
+
+ //! the solution vector
+ SolutionVector x(leafGridView.size(0));
+ tracerProblem->applyInitialSolution(x);
+ auto xOld = x;
+
+ //! the grid variables
+ using GridVariables = GetPropType<TracerTypeTag, Properties::GridVariables>;
+ auto gridVariables = std::make_shared<GridVariables>(tracerProblem, fvGridGeometry);
+ gridVariables->init(x);
+
+ //! get some time loop parameters
+ const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
+ auto dt = getParam<Scalar>("TimeLoop.DtInitial");
+ const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
+
+ //! instantiate time loop
+ auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ //! the assembler with time loop for instationary problem
+ using TracerAssembler = FVAssembler<TracerTypeTag, DiffMethod::analytic, /*implicit=*/false>;
+ auto assembler = std::make_shared<TracerAssembler>(tracerProblem, fvGridGeometry, gridVariables, timeLoop);
+ assembler->setLinearSystem(A, r);
+
+ //! intialize the vtk output module
+ VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, tracerProblem->name());
+ using IOFields = GetPropType<TracerTypeTag, Properties::IOFields>;
+ IOFields::initOutputModule(vtkWriter); // Add model specific output fields
+ using VelocityOutput = GetPropType<TracerTypeTag, Properties::VelocityOutput>;
+ vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
+ vtkWriter.write(0.0);
+
+ /////////////////////////////////////////////////////////////////////////////////////////////////
+ // run instationary non-linear simulation
+ /////////////////////////////////////////////////////////////////////////////////////////////////
+
+ //! set some check points for the time loop
+ timeLoop->setPeriodicCheckPoint(tEnd/10.0);
+
+ //! start the time loop
+ timeLoop->start(); do
+ {
+ // set previous solution for storage evaluations
+ assembler->setPreviousSolution(xOld);
+
+ Dune::Timer assembleTimer;
+ assembler->assembleJacobianAndResidual(x);
+ assembleTimer.stop();
+
+ // solve the linear system A(xOld-xNew) = r
+ Dune::Timer solveTimer;
+ SolutionVector xDelta(x);
+ linearSolver->solve(*A, xDelta, *r);
+ solveTimer.stop();
+
+ // update solution and grid variables
+ updateTimer.reset();
+ x -= xDelta;
+ gridVariables->update(x);
+ updateTimer.stop();
+
+ // statistics
+ 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;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // write vtk output on check points
+ if (timeLoop->isCheckPoint())
+ vtkWriter.write(timeLoop->time());
+
+ // report statistics of this time step
+ timeLoop->reportTimeStep();
+
+ // set new dt
+ timeLoop->setTimeStepSize(dt);
+
+ } while (!timeLoop->finished());
+
+ timeLoop->finalize(leafGridView.comm());
+
+ ////////////////////////////////////////////////////////////
+ // finalize, print dumux message to say goodbye
+ ////////////////////////////////////////////////////////////
+
+ //! print dumux end message
+ if (mpiHelper.rank() == 0)
+ DumuxMessage::print(/*firstCall=*/false);
+
+ return 0;
+
+}
+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;
+}
diff --git a/examples/1ptracer/params.input b/examples/1ptracer/params.input
new file mode 100644
index 0000000000..fb7cadbe32
--- /dev/null
+++ b/examples/1ptracer/params.input
@@ -0,0 +1,21 @@
+[Grid]
+UpperRight = 1 1
+Cells = 50 50
+
+[SpatialParams]
+LensLowerLeft = 0.2 0.2
+LensUpperRight = 0.8 0.8
+
+Permeability = 1e-10 # [m^2]
+PermeabilityLens = 1e-11 # [m^2]
+
+[TimeLoop]
+DtInitial = 10 # [s]
+MaxTimeStepSize = 10
+TEnd = 5000 # [s]
+
+[Problem]
+Name = tracer # name passed to the output routines
+
+[Vtk]
+AddVelocity = true
diff --git a/examples/1ptracer/problem_1p.hh b/examples/1ptracer/problem_1p.hh
new file mode 100644
index 0000000000..150c0b869d
--- /dev/null
+++ b/examples/1ptracer/problem_1p.hh
@@ -0,0 +1,165 @@
+// -*- 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 The properties for the incompressible test
+ */
+
+#ifndef DUMUX_ONEP_TRACER_TEST_PROBLEM_HH
+#define DUMUX_ONEP_TRACER_TEST_PROBLEM_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cctpfa.hh>
+#include <dumux/porousmediumflow/1p/model.hh>
+#include <dumux/porousmediumflow/problem.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh>
+
+#include "spatialparams_1p.hh"
+
+namespace Dumux {
+
+/*!
+ * \ingroup TracerTests
+ * \brief The properties for the incompressible test
+ */
+// forward declarations
+template<class TypeTag>
+class OnePTestProblem;
+
+namespace Properties {
+// Create new type tags
+namespace TTag {
+struct IncompressibleTest { using InheritsFrom = std::tuple<OneP, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::IncompressibleTest> { using type = Dune::YaspGrid<2>; };
+
+// Set the problem type
+template<class TypeTag>
+struct Problem<TypeTag, TTag::IncompressibleTest> { using type = OnePTestProblem<TypeTag>; };
+
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::IncompressibleTest>
+{
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = OnePTestSpatialParams<FVGridGeometry, Scalar>;
+};
+
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::IncompressibleTest> { using type = OnePIncompressibleLocalResidual<TypeTag>; };
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::IncompressibleTest>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Components::SimpleH2O<Scalar> >;
+};
+
+// Enable caching
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::IncompressibleTest> { static constexpr bool value = true; };
+} // end namespace Properties
+
+template<class TypeTag>
+class OnePTestProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using GridView = GetPropType<TypeTag, Properties::GridView>;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using FVElementGeometry = typename GetPropType<TypeTag, Properties::FVGridGeometry>::LocalView;
+ using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
+
+ static constexpr int dimWorld = GridView::dimensionworld;
+
+public:
+ OnePTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry) {}
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ *
+ * \param element The finite element
+ * \param scvf The sub-control volume face
+ */
+ BoundaryTypes boundaryTypes(const Element &element,
+ const SubControlVolumeFace &scvf) const
+ {
+ BoundaryTypes values;
+ const auto globalPos = scvf.ipGlobal();
+
+ Scalar eps = 1.0e-6;
+ if (globalPos[dimWorld-1] < eps || globalPos[dimWorld-1] > this->fvGridGeometry().bBoxMax()[dimWorld-1] - eps)
+ values.setAllDirichlet();
+ else
+ values.setAllNeumann();
+
+ return values;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet control volume.
+ *
+ * \param element The finite element
+ * \param scvf The sub-control volume face
+ *
+ * For this method, the \a values parameter stores primary variables.
+ */
+ PrimaryVariables dirichlet(const Element &element,
+ const SubControlVolumeFace &scvf) const
+ {
+ const auto& pos = scvf.ipGlobal();
+ PrimaryVariables values(0);
+ values[0] = 1.0e+5*(1.1 - pos[dimWorld-1]*0.1);
+ return values;
+ }
+
+ /*!
+ * \brief Returns the temperature \f$\mathrm{[K]}\f$ for an isothermal problem.
+ *
+ * This is not specific to the discretization. By default it just
+ * throws an exception so it must be overloaded by the problem if
+ * no energy equation is used.
+ */
+ Scalar temperature() const
+ {
+ return 283.15; // 10°C
+ }
+
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/examples/1ptracer/problem_tracer.hh b/examples/1ptracer/problem_tracer.hh
new file mode 100644
index 0000000000..929566b151
--- /dev/null
+++ b/examples/1ptracer/problem_tracer.hh
@@ -0,0 +1,231 @@
+// -*- 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 Definition of a problem, for the tracer problem:
+ * A rotating velocity field mixes a tracer band in a porous groundwater reservoir.
+ */
+
+#ifndef DUMUX_TRACER_TEST_PROBLEM_HH
+#define DUMUX_TRACER_TEST_PROBLEM_HH
+
+#include <dune/grid/yaspgrid.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 Definition of a problem, for the tracer problem:
+ * A rotating velocity field mixes a tracer band in a porous groundwater reservoir.
+ */
+template <class TypeTag>
+class TracerTestProblem;
+
+namespace Properties {
+// Create new type tags
+namespace TTag {
+struct TracerTest { using InheritsFrom = std::tuple<Tracer>; };
+struct TracerTestCC { using InheritsFrom = std::tuple<TracerTest, CCTpfaModel>; };
+} // end namespace TTag
+
+// enable caching
+template<class TypeTag>
+struct EnableGridVolumeVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableGridFluxVariablesCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; };
+template<class TypeTag>
+struct EnableFVGridGeometryCache<TypeTag, TTag::TracerTest> { static constexpr bool value = true; };
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::TracerTest> { using type = Dune::YaspGrid<2>; };
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::TracerTest> { using type = TracerTestProblem<TypeTag>; };
+
+// Set the spatial parameters
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::TracerTest>
+{
+ using FVGridGeometry = GetPropType<TypeTag, Properties::FVGridGeometry>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = TracerTestSpatialParams<FVGridGeometry, Scalar>;
+};
+
+// Define whether mole(true) or mass (false) fractions are used
+template<class TypeTag>
+struct UseMoles<TypeTag, TTag::TracerTest> { static constexpr bool value = false; };
+template<class TypeTag>
+struct SolutionDependentMolecularDiffusion<TypeTag, TTag::TracerTestCC> { 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); }
+
+ //! Human readable phase name (index phaseIdx) (for velocity vtk output)
+ static std::string phaseName(int phaseIdx = 0)
+ { return "Groundwater"; }
+
+ //! 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)
+ { return 0.0; }
+};
+
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::TracerTest> { 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_boxtracer -ParameterFile ./test_boxtracer.input</tt> or
+ * <tt>./test_cctracer -ParameterFile ./test_cctracer.input</tt>
+ */
+template <class TypeTag>
+class TracerTestProblem : 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:
+ TracerTestProblem(std::shared_ptr<const FVGridGeometry> fvGridGeom)
+ : ParentType(fvGridGeom)
+ {
+ // 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';
+ }
+
+ /*!
+ * \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;
+ values.setAllNeumann();
+ return values;
+ }
+ // \}
+
+ /*!
+ * \name Volume terms
+ */
+ // \{
+
+ /*!
+ * \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 (globalPos[1] < 0.1 + eps_)
+ {
+ 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;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/examples/1ptracer/spatialparams_1p.hh b/examples/1ptracer/spatialparams_1p.hh
new file mode 100644
index 0000000000..b5ba7b3e42
--- /dev/null
+++ b/examples/1ptracer/spatialparams_1p.hh
@@ -0,0 +1,131 @@
+// -*- 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 The spatial parameters for the incompressible test
+ */
+
+#ifndef DUMUX_INCOMPRESSIBLE_ONEP_TEST_SPATIAL_PARAMS_HH
+#define DUMUX_INCOMPRESSIBLE_ONEP_TEST_SPATIAL_PARAMS_HH
+
+#include <random>
+#include <dumux/porousmediumflow/properties.hh>
+#include <dumux/material/spatialparams/fv1p.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup TracerTests
+ * \brief The spatial parameters for the incompressible test
+ */
+template<class FVGridGeometry, class Scalar>
+class OnePTestSpatialParams
+: public FVSpatialParamsOneP<FVGridGeometry, Scalar,
+ OnePTestSpatialParams<FVGridGeometry, Scalar>>
+{
+ using GridView = typename FVGridGeometry::GridView;
+ using FVElementGeometry = typename FVGridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ParentType = FVSpatialParamsOneP<FVGridGeometry, Scalar,
+ OnePTestSpatialParams<FVGridGeometry, Scalar>>;
+
+ static constexpr int dimWorld = GridView::dimensionworld;
+ using GlobalPosition = typename SubControlVolume::GlobalPosition;
+
+public:
+ using PermeabilityType = Scalar;
+ OnePTestSpatialParams(std::shared_ptr<const FVGridGeometry> fvGridGeometry)
+ : ParentType(fvGridGeometry), K_(fvGridGeometry->gridView().size(0), 0.0)
+ {
+ permeability_ = getParam<Scalar>("SpatialParams.Permeability");
+ permeabilityLens_ = getParam<Scalar>("SpatialParams.PermeabilityLens");
+
+ lensLowerLeft_ = getParam<GlobalPosition>("SpatialParams.LensLowerLeft");
+ lensUpperRight_ =getParam<GlobalPosition>("SpatialParams.LensUpperRight");
+
+ // generate random fields
+ std::mt19937 rand(0);
+ std::lognormal_distribution<Scalar> K(std::log(permeability_), std::log(permeability_)*0.1);
+ std::lognormal_distribution<Scalar> KLens(std::log(permeabilityLens_), std::log(permeabilityLens_)*0.1);
+ for (const auto& element : elements(fvGridGeometry->gridView()))
+ {
+ const auto eIdx = fvGridGeometry->elementMapper().index(element);
+ const auto globalPos = element.geometry().center();
+ K_[eIdx] = isInLens_(globalPos) ? KLens(rand) : K(rand);
+ }
+ }
+
+ /*!
+ * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+ *
+ * \param element The element
+ * \param scv The sub-control volume
+ * \param elemSol The element solution vector
+ * \return the intrinsic permeability
+ */
+ template<class ElementSolution>
+ const PermeabilityType& permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ return K_[scv.dofIndex()];
+
+ // if (isInLens_(scv.dofPosition()))
+ // return permeabilityLens_;
+ // else
+ // return permeability_;
+ }
+
+ /*!
+ * \brief Function for defining the porosity which is possibly solution dependent.
+ *
+ * \return the porosity
+ */
+ Scalar porosityAtPos(const GlobalPosition &globalPos) const
+ { return 0.2; }
+
+ //! Reference to the k field
+ const std::vector<Scalar>& getKField() const
+ { return K_; }
+
+private:
+ bool isInLens_(const GlobalPosition &globalPos) const
+ {
+ for (int i = 0; i < dimWorld; ++i) {
+ if (globalPos[i] < lensLowerLeft_[i] + eps_ || globalPos[i] > lensUpperRight_[i] - eps_)
+ return false;
+ }
+ return true;
+ }
+
+ GlobalPosition lensLowerLeft_;
+ GlobalPosition lensUpperRight_;
+
+ Scalar permeability_, permeabilityLens_;
+
+ std::vector<Scalar> K_;
+
+ static constexpr Scalar eps_ = 1.5e-7;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/examples/1ptracer/spatialparams_tracer.hh b/examples/1ptracer/spatialparams_tracer.hh
new file mode 100644
index 0000000000..08760fed12
--- /dev/null
+++ b/examples/1ptracer/spatialparams_tracer.hh
@@ -0,0 +1,114 @@
+// -*- 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 Definition of the spatial parameters for the tracer problem
+ */
+
+#ifndef DUMUX_TRACER_TEST_SPATIAL_PARAMS_HH
+#define DUMUX_TRACER_TEST_SPATIAL_PARAMS_HH
+
+#include <dumux/porousmediumflow/properties.hh>
+#include <dumux/material/spatialparams/fv1p.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup TracerTests
+ * \brief Definition of the spatial parameters for the tracer problem
+ */
+template<class FVGridGeometry, class Scalar>
+class TracerTestSpatialParams
+: public FVSpatialParamsOneP<FVGridGeometry, Scalar,
+ TracerTestSpatialParams<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 = FVSpatialParamsOneP<FVGridGeometry, Scalar,
+ TracerTestSpatialParams<FVGridGeometry, Scalar>>;
+
+ static const int dimWorld = GridView::dimensionworld;
+ using GlobalPosition = typename Dune::FieldVector<Scalar, dimWorld>;
+
+public:
+
+ TracerTestSpatialParams(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.2; }
+
+ /*!
+ * \brief Defines the dispersivity.
+ *
+ * \param element The finite element
+ * \param scv The sub-control volume
+ * \param elemSol The solution for all dofs of the element
+ */
+ template<class ElementSolution>
+ Scalar dispersivity(const Element &element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ { return 0; }
+
+ //! 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 1000; }
+
+ //! Fluid molar mass
+ Scalar fluidMolarMass(const Element &element,
+ const SubControlVolume& scv) const
+ { return 18.0; }
+
+ Scalar fluidMolarMass(const GlobalPosition &globalPos) const
+ { return 18.0; }
+
+ //! 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; }
+
+private:
+ std::vector<Scalar> volumeFlux_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index 5f0d34d1ee..3a825767d0 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1 +1,2 @@
add_subdirectory(2pinfiltration)
+add_subdirectory(1ptracer)
--
GitLab