diff --git a/test/porousmediumflow/1p/implicit/CMakeLists.txt b/test/porousmediumflow/1p/implicit/CMakeLists.txt
index 170567869d4c069f4565b957e23a8cf8252e1eb2..4b66f73a77d21bb734e83b1dac4c71370215bb91 100644
--- a/test/porousmediumflow/1p/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/1p/implicit/CMakeLists.txt
@@ -8,3 +8,4 @@ add_subdirectory(fracture2d3d)
add_subdirectory(isothermal)
add_subdirectory(nonisothermal)
add_subdirectory(network1d3d)
+add_subdirectory(rootbenchmark)
diff --git a/test/porousmediumflow/1p/implicit/rootbenchmark/CMakeLists.txt b/test/porousmediumflow/1p/implicit/rootbenchmark/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e37556813e76b233b7cec4d62d36101448cbf3e7
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/rootbenchmark/CMakeLists.txt
@@ -0,0 +1,15 @@
+dune_symlink_to_source_files(FILES params.input convergencetest.py)
+
+dumux_add_test(NAME test_1p_rootbenchmark_tpfa
+ LABELS porousmediumflow 1p
+ SOURCES main.cc
+ COMPILE_DEFINITIONS TYPETAG=RootBenchmarkCCTpfa
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_1p_rootbenchmark_tpfa -Problem.EnableGravity false)
+
+dumux_add_test(NAME test_1p_rootbenchmark_tpfa_gravity
+ TARGET test_1p_rootbenchmark_tpfa
+ LABELS porousmediumflow 1p
+ COMPILE_DEFINITIONS TYPETAG=RootBenchmarkCCTpfa
+ COMMAND ./convergencetest.py
+ CMD_ARGS test_1p_rootbenchmark_tpfa test_1p_rootbenchmark_tpfa_gravity.log)
diff --git a/test/porousmediumflow/1p/implicit/rootbenchmark/convergencetest.py b/test/porousmediumflow/1p/implicit/rootbenchmark/convergencetest.py
new file mode 100755
index 0000000000000000000000000000000000000000..b103ea1476e26ff6ebd2d74abeb81b1ff94fbd67
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/rootbenchmark/convergencetest.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python3
+
+from math import *
+import subprocess
+import sys
+
+if len(sys.argv) < 2:
+ sys.stderr.write('Please provide at least a single argument <testname> to the script\n')
+ sys.exit(1)
+
+testname = str(sys.argv[1])
+if len(sys.argv) > 2:
+ if sys.argv[2].endswith(".log"):
+ logfile = str(sys.argv[2])
+ testargs = [str(i) for i in sys.argv][3:]
+ else:
+ logfile = testname + ".log"
+ testargs = [str(i) for i in sys.argv][2:]
+
+# remove the old log file
+subprocess.run(['rm', logfile])
+print("Removed old log file ({})!".format(logfile))
+
+
+# do the runs with different refinement
+for i in [0, 1, 2, 3, 4, 5]:
+ subprocess.run(['./' + testname] + testargs + ['-Problem.Name', logfile.rstrip(".log")] + ['-Grid.Refinement', str(i)])
+
+# check the rates and append them to the log file
+logfile = open(logfile, "r+")
+
+error = []
+hmax = []
+for line in logfile:
+ line = line.strip("\n")
+ line = line.strip(r"\[ConvergenceTest\]")
+ line = line.split()
+ error.append(float(line[2]))
+ hmax.append(float(line[5]))
+
+results = []
+logfile.truncate(0)
+logfile.write("n\thmax\t\terror\t\trate\n")
+logfile.write("-"*50 + "\n")
+for i in range(len(error)-1):
+ if isnan(error[i]) or isinf(error[i]):
+ continue
+ if not (error[i] < 1e-12 or error[i+1] < 1e-12):
+ rate = (log(error[i])-log(error[i+1]))/(log(hmax[i])-log(hmax[i+1]))
+ message = "{}\t{:0.4e}\t{:0.4e}\t{:0.4e}\n".format(i, hmax[i], error[i], rate)
+ logfile.write(message)
+ results.append(rate)
+ else:
+ logfile.write("error: exact solution!?")
+i = len(error)-1
+message = "{}\t{:0.4e}\t{:0.4e}\n\n\n".format(i, hmax[i], error[i])
+logfile.write(message)
+
+logfile.close()
+print("\nComputed the following convergence rates for {}:\n".format(testname))
+subprocess.call(['cat', testname + '.log'])
+
+# check the rates, we expect rates around 2
+for r in results:
+ if int(round(r)) != 2:
+ sys.stderr.write("*"*70 + "\n" + "The convergence rates were not close enough to 2! Test failed.\n" + "*"*70 + "\n")
+ sys.exit(1)
+
+sys.exit(0)
diff --git a/test/porousmediumflow/1p/implicit/rootbenchmark/main.cc b/test/porousmediumflow/1p/implicit/rootbenchmark/main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..34c53eebee7da27a4aef8772493c881774de1113
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/rootbenchmark/main.cc
@@ -0,0 +1,136 @@
+// -*- 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 OnePTests
+ * \brief Root benchmark case Schnepf et al 2020 M3.1 doi: 10.3389/fpls.2020.00316
+ */
+
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/parallel/mpihelper.hh>
+
+#include <dumux/common/properties.hh>
+
+#include <dumux/linear/seqsolverbackend.hh>
+#include <dumux/linear/pdesolver.hh>
+
+#include <dumux/assembly/fvassembler.hh>
+
+#if HAVE_GNUPLOT
+#include <dumux/io/gnuplotinterface.hh>
+#endif
+#include <dumux/io/vtkoutputmodule.hh>
+#include <dumux/io/grid/gridmanager_yasp.hh>
+
+#include "properties.hh"
+
+int main(int argc, char** argv)
+{
+ using namespace Dumux;
+
+ // define the type tag for this problem
+ using TypeTag = Properties::TTag::TYPETAG;
+
+ // initialize MPI, finalize is done automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // 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<TypeTag, Properties::Grid>> gridManager;
+ gridManager.init();
+
+ // we compute on the leaf grid view
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // create the finite volume grid geometry
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
+ gridGeometry->update();
+
+ // the problem (initial and boundary conditions)
+ using Problem = GetPropType<TypeTag, Properties::Problem>;
+ auto problem = std::make_shared<Problem>(gridGeometry);
+
+ // the solution vector
+ using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
+ SolutionVector x(gridGeometry->numDofs());
+
+ // the grid variables
+ using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
+ auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
+ gridVariables->init(x);
+
+ // intialize the vtk output module
+ VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
+ vtkWriter.addVolumeVariable([](const auto& v) { return v.pressure(); }, "p (Pa)");
+ vtkWriter.addVolumeVariable([](const auto& v) { return 100*(v.pressure() - 1e5)/(1000*9.81); }, "psi (cm)");
+ vtkWriter.addField(problem->analyticalPressure(), "p_exact");
+ vtkWriter.addField(problem->analyticalHead(), "psi_exact");
+ vtkWriter.write(0.0);
+
+ // the system assembler
+ using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
+ auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
+
+ // the linear solver
+ using LinearSolver = UMFPackBackend;
+ auto linearSolver = std::make_shared<LinearSolver>();
+
+ // the system solver
+ using Solver = Dumux::LinearPDESolver<Assembler, LinearSolver>;
+ Solver solver(assembler, linearSolver);
+
+ // solve & write output
+ solver.solve(x);
+ vtkWriter.write(1.0);
+
+ // write more output for benchmark plot
+ problem->outputL2Norm(x);
+
+#if HAVE_GNUPLOT
+ // plot with gnuplot
+ bool enablePlot = getParam<bool>("Problem.EnablePlot", false);
+ if (enablePlot)
+ {
+ GnuplotInterface<double> gnuplot;
+ std::vector<double> psi(x.size()), z(x.size());
+ for (const auto& element : elements(leafGridView))
+ {
+ const auto eIdx = gridGeometry->elementMapper().index(element);
+ z[eIdx] = element.geometry().center()[2];
+ psi[eIdx] = 100*(x[eIdx][0] - 1.0e5)/(1000*9.81);
+ }
+ gnuplot.addDataSetToPlot(psi, z, "psi.dat", "w l t 'DuMux'");
+ gnuplot.addDataSetToPlot(problem->analyticalHead(), z, "psi_exact.dat", "w p t 'exact'");
+ gnuplot.plot();
+ }
+#endif
+
+ // print parameters (used and unused)
+ if (mpiHelper.rank() == 0)
+ Parameters::print();
+
+ return 0;
+
+} // end main
diff --git a/test/porousmediumflow/1p/implicit/rootbenchmark/params.input b/test/porousmediumflow/1p/implicit/rootbenchmark/params.input
new file mode 100644
index 0000000000000000000000000000000000000000..af72dc6d06303bf09f5074fd0fc6d139d0253e0f
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/rootbenchmark/params.input
@@ -0,0 +1,20 @@
+[Grid]
+LowerLeft = -0.5
+UpperRight = 0.0
+Cells = 10
+
+[Problem]
+Name = test_1p_rootbenchmark
+EnableGravity = true
+RootCollarPressure = 1900.0 # Pa (psi = -1000cm, p = psi*0.01*9.81*1000 + 1e5)
+SoilPressure = 80380.0 # Pa (psi = -200cm, p = psi*0.01*9.81*1000 + 1e5)
+Radius = 2e-3 # 2mm
+AxialConductivity = 4.32e-2 # cm^3/day (contains factor rho*g)
+RadialConductivity = 1.728e-4 # 1/day (contains factor rho*g)
+
+[Component]
+LiquidDensity = 1000
+LiquidKinematicViscosity = 1.0e-6
+
+[Assembly.NumericDifference]
+PriVarMagnitude = 1e5
diff --git a/test/porousmediumflow/1p/implicit/rootbenchmark/problem.hh b/test/porousmediumflow/1p/implicit/rootbenchmark/problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..0b748697ef5db13e1aa78619bcd4bfc5dbd1d7de
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/rootbenchmark/problem.hh
@@ -0,0 +1,240 @@
+// -*- 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 OnePTests
+ * \brief Root benchmark case Schnepf et al 2020 M3.1 doi: 10.3389/fpls.2020.00316
+ */
+
+#ifndef DUMUX_ROOTBENCHMARK_TEST_PROBLEM_HH
+#define DUMUX_ROOTBENCHMARK_TEST_PROBLEM_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dumux/common/properties.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/common/boundarytypes.hh>
+
+#include <dumux/io/format.hh>
+
+#include <dumux/porousmediumflow/problem.hh>
+
+#include <dumux/material/components/constant.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup OnePTests
+ * \brief Root benchmark case Schnepf et al 2020 M3.1 doi: 10.3389/fpls.2020.00316
+ */
+template <class TypeTag>
+class RootBenchmarkProblem : public PorousMediumFlowProblem<TypeTag>
+{
+ using ParentType = PorousMediumFlowProblem<TypeTag>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using GridView = typename GridGeometry::GridView;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+ static constexpr int dimWorld = GridView::dimensionworld;
+
+ using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+ using BoundaryTypes = Dumux::BoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
+ using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
+
+public:
+ RootBenchmarkProblem(std::shared_ptr<const GridGeometry> gridGeometry)
+ : ParentType(gridGeometry)
+ {
+ name_ = getParam<std::string>("Problem.Name");
+ rootCollarPressure_ = getParam<double>("Problem.RootCollarPressure");
+ soilPressure_ = getParam<double>("Problem.SoilPressure");
+ radius_ = getParam<double>("Problem.Radius");
+ density_ = getParam<double>("Component.LiquidDensity");
+ Kr_ = this->spatialParams().radialConductivity();
+ Kx_ = this->spatialParams().axialConductivity();
+ const auto enableGravity = getParam<bool>("Problem.EnableGravity");
+ const Scalar gravity = enableGravity ? 9.81 : 0.0;
+
+ // cf. Schnepf et al 2020 doi: 10.3389/fpls.2020.00316
+ // Equation is 2*pi*R*Kr(ps - p0) = d/ds(-Kx (dp/ds + rho*g))
+ // Analytical solution is p = ps + A*exp(Cs) + B*exp(-Cs)
+ // where C = sqrt(2*pi*R*Kr/Kx)
+ // Two boundary conditions are needed to determine A and B
+ // Dirichlet at root collar (s=0, p=p0)
+ // (1) A + B = p0 - ps
+ // Neumann no-flow at tip (s=-L, dp/ds + rho*g = 0)
+ // (2) A*C*exp(-CL) - B*C*exp(CL) = -rho*g
+ // Assemble (1) and (2) and solve for A and B
+ Dune::FieldMatrix<double, 2, 2> M;
+ Dune::FieldVector<double, 2> b;
+ using std::exp; using std::sqrt;
+ const auto C = sqrt(2.0*M_PI*radius_*Kr_/Kx_); C_ = C;
+ const auto L = this->gridGeometry().bBoxMax()[dimWorld-1] - this->gridGeometry().bBoxMin()[dimWorld-1];
+ M[0][0] = 1.0; M[0][1] = 1.0; b[0] = rootCollarPressure_ - soilPressure_;
+ M[1][0] = C*exp(-C*L); M[1][1] = -C*exp(C*L); b[1] = -density_*gravity;
+ M.solve(d_, b);
+
+ std::cout << "Computed analytical solution with:\n";
+ std::cout << Fmt::format("A or d1 = {}, B or d2 = {}, c = {}, √c = {}\n", d_[0], d_[1], C*C, C);
+ std::cout << Fmt::format("Kx = {}, Kr = {}\n", Kx_, Kr_);
+
+ const auto psiCollar = 100*(rootCollarPressure_ - 1e5)/(1000*9.81);
+ const auto psiSoil = 100*(soilPressure_ - 1e5)/(1000*9.81);
+ const auto psi0 = 100*(analyticalSolution(gridGeometry->bBoxMax()) - 1e5)/(1000*9.81);
+ const auto psiTip = 100*(analyticalSolution(gridGeometry->bBoxMin()) - 1e5)/(1000*9.81);
+ std::cout << Fmt::format("ψ(0) = {}, ψ(-0.5) = {}, ψ0 = {}, ψs = {} cm\n", psi0, psiTip, psiCollar, psiSoil);
+
+ analyticalPressures_.resize(gridGeometry->gridView().size(0));
+ analyticalHead_.resize(gridGeometry->gridView().size(0));
+ for (const auto& element : elements(gridGeometry->gridView()))
+ {
+ const auto eIdx = gridGeometry->elementMapper().index(element);
+ const auto& globalPos = element.geometry().center();
+ analyticalPressures_[eIdx] = analyticalSolution(globalPos);
+ analyticalHead_[eIdx] = 100*(analyticalPressures_[eIdx] - 1e5)/(1000*9.81);
+ }
+ }
+
+ /*!
+ * \brief The problem name.
+ */
+ const std::string& name() const
+ { return name_; }
+
+ /*!
+ * \brief Returns the temperature within the domain in [K].
+ */
+ Scalar temperature() const
+ { return 273.15 + 10.0; }
+
+ /*!
+ * \brief Returns how much the domain is extruded at a given sub-control volume.
+ * Assume circular cross section with given radius
+ */
+ Scalar extrusionFactorAtPos(const GlobalPosition& globalPos) const
+ { return M_PI*radius_*radius_; }
+
+ /*!
+ * \brief Specifies which kind of boundary condition should be
+ * used for which equation on a given boundary control volume.
+ */
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ {
+ BoundaryTypes bcTypes;
+ if (globalPos[dimWorld-1] > this->gridGeometry().bBoxMax()[dimWorld-1] - eps_)
+ bcTypes.setAllDirichlet(); // Dirichlet at root collar
+ else
+ bcTypes.setAllNeumann(); // No-flow at root tips
+ return bcTypes;
+ }
+
+ /*!
+ * \brief Evaluates the boundary conditions for a Dirichlet control volume.
+ */
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ { return PrimaryVariables(rootCollarPressure_); }
+
+ /*!
+ * \brief Evaluates the source term for all phases within a given sub-control volume.
+ * Positive values mean that the conserved quantity is created, negative ones mean that it vanishes.
+ * E.g. for the mass balance that would be a mass rate in \f$ [ kg / (m^3 \cdot s)] \f$.
+ */
+ template<class ElementVolumeVariables>
+ NumEqVector source(const Element& element,
+ const FVElementGeometry& fvGeometry,
+ const ElementVolumeVariables& elemVolVars,
+ const SubControlVolume& scv) const
+ {
+ NumEqVector values(0.0);
+ const auto& volVars = elemVolVars[scv];
+ // source term with static soil pressure
+ values[0] = density_*2.0*M_PI*radius_*Kr_*(soilPressure_ - volVars.pressure());
+ // make this volume-specific source
+ values[0] /= volVars.extrusionFactor();
+ return values;
+ }
+
+ template<class SolutionVector>
+ void outputL2Norm(const SolutionVector& solution) const
+ {
+ // calculate the discrete L2-Norm and maximum element size
+ Scalar l2Norm = 0.0;
+ Scalar l2NormNormalization = 0.0;
+ Scalar hMax = 0.0;
+
+ const auto& gg = this->gridGeometry();
+ for (const auto& element : elements(gg.gridView()))
+ {
+ const auto eIdx = gg.elementMapper().index(element);
+ const auto geo = element.geometry();
+ const auto globalPos = geo.center();
+ const auto pe = analyticalSolution(globalPos);
+ const auto p = solution[eIdx][0];
+ const auto dV = geo.volume();
+
+ l2Norm += (p - pe)*(p - pe)*dV;
+ l2NormNormalization += pe*pe*dV;
+ hMax = std::max(dV, hMax);
+ }
+
+ l2Norm = std::sqrt(l2Norm/l2NormNormalization);
+
+ // write the norm into a log file
+ std::ofstream logFile(this->name() + ".log", std::ios::app);
+ logFile << "[ConvergenceTest] L2-norm(pressure) = " << l2Norm << " hMax = " << hMax << std::endl;
+ }
+
+ Scalar analyticalSolution(const GlobalPosition& globalPos) const
+ {
+ using std::exp;
+ const auto s = globalPos[dimWorld-1];
+ return soilPressure_ + d_[0]*exp(s*C_) + d_[1]*exp(-s*C_);
+ }
+
+ const std::vector<Scalar>& analyticalPressure() const
+ { return analyticalPressures_; }
+
+ const std::vector<Scalar>& analyticalHead() const
+ { return analyticalHead_; }
+
+private:
+ Scalar rootCollarPressure_, soilPressure_;
+ Scalar radius_;
+ Scalar density_;
+ Scalar Kr_, Kx_;
+
+ static constexpr Scalar eps_ = 1e-8;
+ std::string name_;
+
+ // analytical solution
+ Dune::FieldVector<Scalar, 2> d_;
+ Scalar C_;
+
+ std::vector<Scalar> analyticalPressures_;
+ std::vector<Scalar> analyticalHead_;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/1p/implicit/rootbenchmark/properties.hh b/test/porousmediumflow/1p/implicit/rootbenchmark/properties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..00760a570fbfc0e34dd138cf933b0a07142c68df
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/rootbenchmark/properties.hh
@@ -0,0 +1,80 @@
+// -*- 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 OnePTests
+ * \brief Root benchmark
+ */
+
+#ifndef DUMUX_ONEP_ROOT_BENCHMARK_PROPERTIES_HH
+#define DUMUX_ONEP_ROOT_BENCHMARK_PROPERTIES_HH
+
+#include <dune/grid/yaspgrid.hh>
+
+#include <dumux/discretization/cctpfa.hh>
+
+#include <dumux/porousmediumflow/1p/model.hh>
+
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidsystems/1pliquid.hh>
+
+#include "problem.hh"
+#include "spatialparams.hh"
+
+namespace Dumux::Properties {
+
+// Create new type tags
+namespace TTag {
+struct RootBenchmark { using InheritsFrom = std::tuple<OneP>; };
+struct RootBenchmarkCCTpfa { using InheritsFrom = std::tuple<RootBenchmark, CCTpfaModel>; };
+} // end namespace TTag
+
+// Set the grid type
+template<class TypeTag>
+struct Grid<TypeTag, TTag::RootBenchmark>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<Scalar, 1>>;
+};
+
+// Set the problem property
+template<class TypeTag>
+struct Problem<TypeTag, TTag::RootBenchmark>
+{ using type = RootBenchmarkProblem<TypeTag>; };
+
+// Set the spatial parameters
+template<class TypeTag>
+struct SpatialParams<TypeTag, TTag::RootBenchmark>
+{
+ using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = RootBenchmarkSpatialParams<GridGeometry, Scalar>;
+};
+
+// the fluid system
+template<class TypeTag>
+struct FluidSystem<TypeTag, TTag::RootBenchmark>
+{
+ using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+ using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
+};
+
+} // end namespace Dumux::Properties
+
+#endif
diff --git a/test/porousmediumflow/1p/implicit/rootbenchmark/spatialparams.hh b/test/porousmediumflow/1p/implicit/rootbenchmark/spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d6be25d3a7b6cad63bf2383c675c0ee5d696c9d3
--- /dev/null
+++ b/test/porousmediumflow/1p/implicit/rootbenchmark/spatialparams.hh
@@ -0,0 +1,108 @@
+// -*- 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 OnePTests
+ * \brief Root benchmark
+ */
+
+#ifndef DUMUX_ONEP_ROOTBENCHMARK_SPATIALPARAMS_HH
+#define DUMUX_ONEP_ROOTBENCHMARK_SPATIALPARAMS_HH
+
+#include <dumux/porousmediumflow/properties.hh>
+#include <dumux/material/spatialparams/fv1p.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup OnePTests
+ * \brief Root benchmark spatial params
+ */
+template<class GridGeometry, class Scalar>
+class RootBenchmarkSpatialParams
+: public FVSpatialParamsOneP<GridGeometry, Scalar,
+ RootBenchmarkSpatialParams<GridGeometry, Scalar>>
+{
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename FVElementGeometry::SubControlVolume;
+ using GridView = typename GridGeometry::GridView;
+ using Element = typename GridView::template Codim<0>::Entity;
+ using ParentType = FVSpatialParamsOneP<GridGeometry, Scalar,
+ RootBenchmarkSpatialParams<GridGeometry, Scalar>>;
+
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ // export permeability type
+ using PermeabilityType = Scalar;
+
+ RootBenchmarkSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
+ : ParentType(gridGeometry)
+ {
+ radius_ = getParam<double>("Problem.Radius");
+ const auto rho = getParam<double>("Component.LiquidDensity");
+ const auto mu = getParam<double>("Component.LiquidKinematicViscosity")*rho;
+ Kx_ = getParam<double>("Problem.AxialConductivity")/(86400.0*1000*9.81)*1e-6; // cm^3/day -> m^4/(Pa*s)
+ Kr_ = getParam<double>("Problem.RadialConductivity")/(86400.0*1000*9.81); // 1/day -> m/(Pa*s)
+ permeability_ = Kx_*mu/(M_PI*radius_*radius_);
+ }
+
+ /*!
+ * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$.
+ */
+ template<class ElementSolution>
+ PermeabilityType permeability(const Element& element,
+ const SubControlVolume& scv,
+ const ElementSolution& elemSol) const
+ {
+ return permeability_;
+ }
+
+ /*!
+ * \brief Returns the root radius in \f$[m]\f$.
+ */
+ Scalar radius() const
+ { return radius_; }
+
+ /*!
+ * \brief The axial conductivity \f$[m^4 Pa^{-1} s^{-1}]\f$.
+ */
+ Scalar axialConductivity() const
+ { return Kx_; }
+
+ /*!
+ * \brief The radial conductivity \f$[m Pa^{-1} s^{-1}]\f$.
+ */
+ Scalar radialConductivity() const
+ { return Kr_; }
+
+ /*!
+ * \brief The porosity \f$\mathrm{[-]}\f$.
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return 1.0; }
+
+private:
+ Scalar radius_, permeability_, Kx_, Kr_;
+ static constexpr Scalar eps_ = 1e-8;
+};
+
+} // end namespace Dumux
+
+#endif