diff --git a/test/porousmediumflow/richards/CMakeLists.txt b/test/porousmediumflow/richards/CMakeLists.txt
index 6597b11e79290f63aa891e2973d4a37aecf3a090..44b49bdddda7e463638175001890790478032a4a 100644
--- a/test/porousmediumflow/richards/CMakeLists.txt
+++ b/test/porousmediumflow/richards/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory(analytical)
+add_subdirectory(annulus)
add_subdirectory(benchmarks)
add_subdirectory(lens)
add_subdirectory(nonisothermal)
diff --git a/test/porousmediumflow/richards/annulus/CMakeLists.txt b/test/porousmediumflow/richards/annulus/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..85be915d7ba95c394ce0697f6848dc53e1cd83bf
--- /dev/null
+++ b/test/porousmediumflow/richards/annulus/CMakeLists.txt
@@ -0,0 +1,25 @@
+dune_symlink_to_source_files(FILES "params.input")
+
+dumux_add_test(NAME test_richards_annulus
+ LABELS porousmediumflow richards
+ SOURCES main.cc
+ COMMAND ${CMAKE_CURRENT_BINARY_DIR}/test_richards_annulus
+ CMD_ARGS params.input -Problem.Name test_richards_annulus -SpatialParams.SoilType Clay)
+
+dumux_add_test(NAME test_richards_annulus_loam
+ LABELS porousmediumflow richards
+ TARGET test_richards_annulus
+ COMMAND ${CMAKE_CURRENT_BINARY_DIR}/test_richards_annulus
+ CMD_ARGS params.input -Problem.Name test_richards_annulus -SpatialParams.SoilType Loam)
+
+dumux_add_test(NAME test_richards_annulus_instationary
+ LABELS porousmediumflow richards
+ SOURCES main_instationary.cc
+ COMMAND ${CMAKE_CURRENT_BINARY_DIR}/test_richards_annulus_instationary
+ CMD_ARGS params.input -Problem.Name test_richards_annulus_instationary -SpatialParams.SoilType Clay)
+
+dumux_add_test(NAME test_richards_annulus_instationary_loam
+ LABELS porousmediumflow richards
+ TARGET test_richards_annulus_instationary
+ COMMAND ${CMAKE_CURRENT_BINARY_DIR}/test_richards_annulus_instationary
+ CMD_ARGS params.input -Problem.Name test_richards_annulus_instationary -SpatialParams.SoilType Loam)
diff --git a/test/porousmediumflow/richards/annulus/README.md b/test/porousmediumflow/richards/annulus/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..e9ac5cafafbbd6def47a700d6521409b436b14ab
--- /dev/null
+++ b/test/porousmediumflow/richards/annulus/README.md
@@ -0,0 +1,29 @@
+# Richards test on annulus
+
+This may corresponds to a rhizosphere model, where we have water uptake
+by a root on the inner boundary. By transforming the Richards equation
+using a Kirchhoff transformation, the flux term becomes linear and
+we can compute an analytic solution for the stationary equation.
+Under some assumptions on the storage term, we can also compute
+approximative solutions to the instationary equation. Here, we
+test convergence against the analytical solution by implementing
+storage approximation as volumetric source terms (possibly nonlinear).
+
+The benchmark case is described in https://doi.org/10.3389/fpls.2020.00316
+See problem documentation for the analytic solution and how it's derived.
+
+# Tests
+
+There are two (four) tests implemented.
+
+* `test_richards_annulus` is a convergence
+test against the stationary analytical solution when the root-soil interface
+pressure reached the wilting point pressure (-15000 cm pressure head).
+The test runs for loam and clay and makes sure the convergence rate is better than 1.95.
+
+* `test_richards_annulus_instationary` is a time-dependent test. At each time step,
+the analytical solution corresponding to the steady-rate approximation with root-soil
+pressure given by the numerical solution is computed. This solution corresponds to
+a runtime (since the uptake rate is constant). The numerical runtime and the analytically
+approximated runtime are compared. The analytic approximation is very accurate and the test
+makes sure the difference is below 5%.
diff --git a/test/porousmediumflow/richards/annulus/main.cc b/test/porousmediumflow/richards/annulus/main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..4ca5cc906b77ec47759ee94fbf105f602fba9bd1
--- /dev/null
+++ b/test/porousmediumflow/richards/annulus/main.cc
@@ -0,0 +1,181 @@
+// -*- 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 . *
+ *****************************************************************************/
+
+#include
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include "properties.hh"
+
+int main(int argc, char** argv)
+{
+ using namespace Dumux;
+
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+
+ // We parse the command line arguments.
+ Parameters::init(argc, argv);
+
+ // Initialize timer
+ Dune::Timer timer;
+
+ // Convenience alias for the type tag of the problem.
+ using TypeTag = Properties::TTag::RichardsAnnulus;
+
+ // The grid manager can be used to create a grid from the input file
+ using Grid = GetPropType;
+ GridManager gridManager;
+ gridManager.init();
+
+ // We compute on the leaf grid view.
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // instantiate the grid geometry
+ using GridGeometry = GetPropType;
+ auto gridGeometry = std::make_shared(leafGridView);
+
+ // Initialize the problem and grid variables
+ using Problem = GetPropType;
+ auto problem = std::make_shared(gridGeometry);
+
+ // We define a function to update the discrete analytical solution vector
+ // using the exactSolution() function in the problem
+ const auto updateAnalyticalSolution = [&](auto& pExact)
+ {
+ pExact.resize(gridGeometry->numDofs());
+ for (const auto& element : elements(gridGeometry->gridView()))
+ {
+ auto fvGeometry = localView(*gridGeometry);
+ fvGeometry.bindElement(element);
+ for (auto&& scv : scvs(fvGeometry))
+ pExact[scv.dofIndex()] = problem->exactSolution(scv.dofPosition());
+ }
+ };
+
+ // instantiate and initialize the discrete and exact solution vectors
+ using SolutionVector = GetPropType;
+ SolutionVector p(gridGeometry->numDofs()); problem->applyInitialSolution(p);
+ SolutionVector pExact; updateAnalyticalSolution(pExact);
+
+ // instantiate and initialize the grid variables
+ using GridVariables = GetPropType;
+ auto gridVariables = std::make_shared(problem, gridGeometry);
+ gridVariables->init(p);
+
+ // initialize VTK output
+ VtkOutputModule vtkWriter(*gridVariables, p, problem->name());
+ GetPropType::initOutputModule(vtkWriter);
+ vtkWriter.addField(pExact, "pExact"); // add the exact solution to the output fields
+ vtkWriter.addVolumeVariable([&](const auto& v){
+ return (v.pressure(0) - problem->nonwettingReferencePressure())*100/(v.density(0)*9.81); },
+ "head in cm"
+ );
+ vtkWriter.write(-1.0);
+
+ // instantiate the solver
+ using Assembler = FVAssembler;
+ auto assembler = std::make_shared(problem, gridGeometry, gridVariables);
+
+ using LinearSolver = ILU0BiCGSTABBackend;
+ auto linearSolver = std::make_shared();
+ NewtonSolver solver(assembler, linearSolver);
+
+ // Solution of the problem and error computation
+ solver.solve(p);
+ vtkWriter.write(0.0);
+
+ const auto initialPressure = problem->headToPressure(getParam("Problem.InitialPressureHeadInCm", -100));
+ const auto initialSaturation = problem->saturation(initialPressure);
+ std::cout << "Initial uniform profile with saturation: " << initialSaturation << std::endl;
+ std::cout << "Time until this profile is reached: "
+ << problem->computeTime(p, *gridVariables, initialSaturation)/(24*60*60)
+ << " days" << std::endl;
+
+ // container to store the L2 errors for the different refinements
+ const int numRefinements = getParam("Grid.RefinementSteps");
+ std::vector l2Errors(numRefinements);
+
+ // use third order error integration
+ constexpr int orderQuadratureRule = 2;
+
+ // compute initial L2 error
+ l2Errors[0] = integrateL2Error(*gridGeometry, p, pExact, orderQuadratureRule);
+ for (int stepIdx = 1; stepIdx < numRefinements; stepIdx++)
+ {
+ // Globally refine the grid once
+ gridManager.grid().globalRefine(1);
+
+ // update the grid geometry, the grid variables and
+ // the solution vectors now that the grid has been refined
+ gridGeometry->update(gridManager.grid().leafGridView());
+ p.resize(gridGeometry->numDofs());
+ problem->applyInitialSolution(p);
+ updateAnalyticalSolution(pExact);
+ gridVariables->updateAfterGridAdaption(p);
+
+ // this recreates the linear system, i.e. the sizes of
+ // the right hand side vector and the Jacobian matrix,
+ // and its sparsity pattern.
+ assembler->updateAfterGridAdaption();
+
+ // solve problem on refined grid
+ solver.solve(p);
+ vtkWriter.write(stepIdx);
+
+ std::cout << "Time until this profile is reached: "
+ << problem->computeTime(p, *gridVariables, initialSaturation)/(24*60*60)
+ << " days" << std::endl;
+
+ // Calculate the L2 error using the numerical solution
+ l2Errors[stepIdx] = integrateL2Error(*gridGeometry, p, pExact, orderQuadratureRule);
+
+ // Print the error and convergence rate
+ const auto rate = std::log(l2Errors[stepIdx]/l2Errors[stepIdx-1])/std::log(0.5);
+ const auto numDofs = gridGeometry->numDofs();
+ std::cout << Fmt::format(
+ "L2 error: {}, rate: {} ({:>5} dofs)\n",
+ l2Errors[stepIdx], rate, numDofs
+ );
+
+ if (rate < 1.95)
+ DUNE_THROW(Dune::Exception, "Convergence rate too small (< 1.95): " << rate);
+ }
+
+ std::cout << "Simulation took " << timer.elapsed() << " second." << std::endl;
+
+ return 0;
+}
diff --git a/test/porousmediumflow/richards/annulus/main_instationary.cc b/test/porousmediumflow/richards/annulus/main_instationary.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d319b11a7349a9029748e185b84e9dcf2b354013
--- /dev/null
+++ b/test/porousmediumflow/richards/annulus/main_instationary.cc
@@ -0,0 +1,204 @@
+// -*- 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 . *
+ *****************************************************************************/
+
+#include
+
+#include
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include "properties.hh"
+
+int main(int argc, char** argv)
+{
+ using namespace Dumux;
+
+ // maybe initialize MPI and/or multithreading backend
+ Dumux::initialize(argc, argv);
+
+ // We parse the command line arguments.
+ Parameters::init(argc, argv);
+
+ // Initialize timer
+ Dune::Timer timer;
+
+ // Convenience alias for the type tag of the problem.
+ using TypeTag = Properties::TTag::RichardsAnnulus;
+
+ // The grid manager can be used to create a grid from the input file
+ using Grid = GetPropType;
+ GridManager gridManager;
+ gridManager.init();
+
+ // We compute on the leaf grid view.
+ const auto& leafGridView = gridManager.grid().leafGridView();
+
+ // instantiate the grid geometry
+ using GridGeometry = GetPropType;
+ auto gridGeometry = std::make_shared(leafGridView);
+
+ // Initialize the problem and grid variables
+ using Problem = GetPropType;
+ auto problem = std::make_shared(gridGeometry);
+ problem->disableSource();
+ problem->disableStationaryInitialSolution();
+ problem->enableInnerNeumannBC();
+
+ // We define a function to update the discrete analytical solution vector
+ // using the exactSolution() function in the problem
+ // Note that this refers to the exact solution in the stationary case
+ // In the instationary case, this analytical solution is just an approximation
+ const auto updateAnalyticalSolution = [&](auto& pAnalyticApproximation, const double innerPsi)
+ {
+ pAnalyticApproximation.resize(gridGeometry->numDofs());
+ for (const auto& element : elements(gridGeometry->gridView()))
+ {
+ auto fvGeometry = localView(*gridGeometry);
+ fvGeometry.bindElement(element);
+ for (auto&& scv : scvs(fvGeometry))
+ pAnalyticApproximation[scv.dofIndex()] = problem->exactSolution(scv.dofPosition(), innerPsi);
+ }
+ };
+
+ // instantiate and initialize the discrete and exact solution vectors
+ using SolutionVector = GetPropType;
+ SolutionVector p(gridGeometry->numDofs()); problem->applyInitialSolution(p);
+ const auto initialPressure = problem->headToPressure(getParam("Problem.InitialPressureHeadInCm", -100));
+ SolutionVector pAnalyticApproximation; updateAnalyticalSolution(pAnalyticApproximation, problem->pressureToPsi(initialPressure));
+ auto pOld = p;
+
+ // instantiate and initialize the grid variables
+ using GridVariables = GetPropType;
+ auto gridVariables = std::make_shared(problem, gridGeometry);
+ gridVariables->init(p);
+
+ auto gridVariablesAnalytic = std::make_shared(problem, gridGeometry);
+ gridVariablesAnalytic->init(pAnalyticApproximation);
+ const auto initialSaturation = problem->saturation(initialPressure);
+
+ // get some time loop parameters
+ const auto tEnd = getParam("TimeLoop.TEnd", 24*60*60*40);
+ const auto maxDt = getParam("TimeLoop.MaxTimeStepSize", 0.5*24*60*60);
+ auto dt = getParam("TimeLoop.DtInitial", 0.5*24*60*60);
+
+ // instantiate time loop
+ auto timeLoop = std::make_shared>(0.0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+
+ std::vector storageDerivative(leafGridView.size(Grid::dimension), 0.0);
+
+ // initialize VTK output
+ VtkOutputModule vtkWriter(*gridVariables, p, problem->name());
+ GetPropType::initOutputModule(vtkWriter);
+ vtkWriter.addField(pAnalyticApproximation, "pAnalyticApproximation");
+ vtkWriter.addField(storageDerivative, "dS/dt");
+ vtkWriter.addVolumeVariable([&](const auto& v){
+ return (v.pressure(0) - problem->nonwettingReferencePressure())*100/(v.density(0)*9.81); },
+ "head in cm"
+ );
+ vtkWriter.write(0.0);
+
+ // instantiate the solver
+ using Assembler = FVAssembler;
+ auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, pOld);
+
+ using LinearSolver = ILU0BiCGSTABBackend;
+ auto linearSolver = std::make_shared();
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+
+ // time loop
+ bool finished = false;
+ const auto wiltingPressure = problem->headToPressure(-15000);
+ timeLoop->start(); do
+ {
+ // solve the non-linear system with time step control
+ nonLinearSolver.solve(p, *timeLoop);
+
+ auto interfacePressure = p[0][0];
+ while (interfacePressure < 1.001*wiltingPressure || interfacePressure > 1.1e5)
+ {
+ p = pOld;
+ assembler->resetTimeStep(p);
+ const auto dt = timeLoop->timeStepSize() * 0.5;
+ std::cout << Fmt::format("Retry with time step size: {} s", dt) << std::endl;
+ timeLoop->setTimeStepSize(dt);
+ nonLinearSolver.solve(p, *timeLoop);
+ interfacePressure = p[0][0];
+ }
+
+ std::cout << Fmt::format("Root-soil interface pressure: {} Pa", interfacePressure) << std::endl;
+ problem->updateStorageDerivative(p, pOld, *gridVariables, timeLoop->timeStepSize(), storageDerivative);
+
+ // make the new solution the old solution
+ pOld = p;
+ gridVariables->advanceTimeStep();
+
+ // advance to the time loop to the next step
+ timeLoop->advanceTimeStep();
+
+ // compute and compare with analytic approximation
+ updateAnalyticalSolution(pAnalyticApproximation, problem->pressureToPsi(interfacePressure));
+ gridVariablesAnalytic->update(pAnalyticApproximation);
+ const auto timePrediction = problem->computeTime(pAnalyticApproximation, *gridVariablesAnalytic, initialSaturation);
+ std::cout << Fmt::format("It took {} days until this pressure was reached\n", timeLoop->time()/86400.0);
+ std::cout << Fmt::format("The analytical approximation predicted {} days\n", timePrediction/86400.0);
+
+ // write vtk output
+ vtkWriter.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()));
+
+ // check if we are done
+ if (interfacePressure <= wiltingPressure)
+ finished = true;
+
+ finished |= timeLoop->finished();
+
+ } while (!finished);
+
+ const auto timePrediction = problem->computeTime(pAnalyticApproximation, *gridVariablesAnalytic, initialSaturation);
+ if (auto d = std::abs(timeLoop->time()-timePrediction)/timePrediction; d > 0.05)
+ DUNE_THROW(Dune::Exception, "Deviation between estimate and numerical solution too big: " << d*100 << "%");
+
+ timeLoop->finalize(leafGridView.comm());
+ nonLinearSolver.report();
+ Parameters::print();
+
+ return 0;
+}
diff --git a/test/porousmediumflow/richards/annulus/params.input b/test/porousmediumflow/richards/annulus/params.input
new file mode 100644
index 0000000000000000000000000000000000000000..afe39435cb62d9923d0778ed5d65261ef1c8994f
--- /dev/null
+++ b/test/porousmediumflow/richards/annulus/params.input
@@ -0,0 +1,62 @@
+[Problem]
+EnableGravity = false
+Name = test_richards_annulus
+InnerFlux = 0.1 # 0.05 # cm/day
+InitialPressureHeadInCm = -100.0
+EnableOuterNeumann = true
+EnableInnerNeumann = false
+
+[Grid]
+Cells0 = 300
+Grading0 = 0.95
+Positions0 = 0.02e-2 0.6e-2
+RefinementSteps = 3
+
+[Flux]
+UpwindWeight = 0.5
+
+[Newton]
+MaxRelativeShift = 1e-7
+
+[SpatialParams]
+SoilType = Clay # Loam # Sand
+Temperature = 283.15
+
+[SpatialParams.Clay]
+Name = Clay
+Swr = 0.25
+Snr = 0.0
+VanGenuchtenPcLowSweThreshold = 0.0
+VanGenuchtenPcHighSweThreshold = 1.0
+VanGenuchtenKrwHighSweThreshold = 1.0
+VanGenuchtenAlpha = 1.01936799e-4
+VanGenuchtenN = 1.1
+VanGenuchtenL = 0.5
+Porosity = 0.40
+Permeability = 1.1798241e-13
+
+[SpatialParams.Loam]
+Name = Loam
+Swr = 0.18604651162
+Snr = 0.0
+VanGenuchtenPcLowSweThreshold = 0.0
+VanGenuchtenPcHighSweThreshold = 1.0
+VanGenuchtenKrwHighSweThreshold = 1.0
+VanGenuchtenAlpha = 4.07747197e-4
+VanGenuchtenN = 1.6
+VanGenuchtenL = 0.5
+Porosity = 0.43
+Permeability = 5.8991203e-13
+
+[SpatialParams.Sand]
+Name = Sand
+Swr = 0.10465116279
+Snr = 0.0
+VanGenuchtenPcLowSweThreshold = 0.0
+VanGenuchtenPcHighSweThreshold = 1.0
+VanGenuchtenKrwHighSweThreshold = 1.0
+VanGenuchtenAlpha = 1.52905199e-3
+VanGenuchtenN = 3.0
+VanGenuchtenL = 0.5
+Porosity = 0.43
+Permeability = 1.1798241e-11
diff --git a/test/porousmediumflow/richards/annulus/problem.hh b/test/porousmediumflow/richards/annulus/problem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..4943b6c5f9e30ff783ab76785082368c411e97cb
--- /dev/null
+++ b/test/porousmediumflow/richards/annulus/problem.hh
@@ -0,0 +1,323 @@
+// -*- 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 . *
+ *****************************************************************************/
+
+#ifndef DUMUX_RICHARDS_ANNULUS_PROBLEM_HH
+#define DUMUX_RICHARDS_ANNULUS_PROBLEM_HH
+
+#include
+#include
+
+#include
+#include
+#include
+#include
+
+#include
+#include
+#include
+
+#include
+
+namespace Dumux {
+
+template
+class RichardsAnnulusProblem : public PorousMediumFlowProblem
+{
+ using ParentType = PorousMediumFlowProblem;
+
+ using GridGeometry = GetPropType;
+ using Indices = typename GetPropType::Indices;
+ using Scalar = GetPropType;
+ using PrimaryVariables = GetPropType;
+
+ using NumEqVector = Dumux::NumEqVector;
+ using BoundaryTypes = Dumux::BoundaryTypes;
+ using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
+ using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
+
+public:
+ RichardsAnnulusProblem(std::shared_ptr gridGeometry)
+ : ParentType(gridGeometry)
+ {
+ // The default parameters correspond to the benchmark cases C1.1
+ // of Schnepf et al (2020) https://doi.org/10.3389/fpls.2020.00316
+
+ // material parameters
+ k_ = getParam(this->spatialParams().paramGroup() + ".Permeability");
+ static_assert(!FluidSystem::isCompressible(0), "Assumes incompressible fluid");
+ rho_ = FluidSystem::H2O::liquidDensity(0, 0);
+ static_assert(FluidSystem::viscosityIsConstant(0), "Assumes constant viscosity fluid");
+ mu_ = FluidSystem::H2O::liquidViscosity(0, 0);
+
+ // geometry
+ innerRadius_ = gridGeometry->bBoxMin()[0];
+ outerRadius_ = gridGeometry->bBoxMax()[0];
+
+ // boundary conditions (volumetric flow rate in m^2/s, pressure/psi in Pa)
+ const auto fi = getParam("Problem.InnerFlux", 0.1); // positive means outflow (cm/day)
+ innerFlowRate_ = fi*2*M_PI*innerRadius_*0.01/(24*60*60);
+ const auto fo = getParam("Problem.OuterFlux", 0.0); // positive means outflow (cm/day)
+ outerFlowRate_ = fo*2*M_PI*outerRadius_*0.01/(24*60*60);
+
+ const auto innerPressureHead = getParam("Problem.InnerPressureHeadInCm", -15000);
+ const auto innerPressure = headToPressure(innerPressureHead);
+ innerPsi_ = applyKirchhoffTransformation_(innerPressure);
+ source_ = (innerFlowRate_ + outerFlowRate_)/(M_PI*(innerRadius_*innerRadius_ - outerRadius_*outerRadius_));
+
+ // boundary condition types
+ enableOuterNeumann_ = getParam("Problem.EnableOuterNeumann", true);
+ enableInnerNeumann_ = getParam("Problem.EnableInnerNeumann", false);
+
+ std::cout << "Inner flow rate: " << innerFlowRate_ << " m^2/s" << std::endl;
+ std::cout << "Outer flow rate: " << outerFlowRate_ << " m^2/s" << std::endl;
+ std::cout << "Critical inner pressure: " << innerPressure << " Pa "
+ << "(-> head: " << innerPressureHead << " cm)" << std::endl;
+
+ // initial condition
+ const auto initialHead = getParam("Problem.InitialPressureHeadInCm", -100);
+ initialPressure_ = headToPressure(initialHead);
+ }
+
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
+ {
+ BoundaryTypes values;
+ if (enableOuterNeumann_ && globalPos[0] > this->gridGeometry().bBoxMax()[0] - 1e-10)
+ values.setAllNeumann();
+ else if (enableInnerNeumann_ && globalPos[0] < this->gridGeometry().bBoxMin()[0] + 1e-10)
+ values.setAllNeumann();
+ else
+ values.setAllDirichlet();
+ return values;
+ }
+
+ // negative means extraction
+ NumEqVector sourceAtPos(const GlobalPosition& globalPos) const
+ {
+ if (enableSource_)
+ return { -rho_*source_ };
+ else
+ return { 0.0 };
+ }
+
+ // positive means outflow
+ NumEqVector neumannAtPos(const GlobalPosition& globalPos) const
+ {
+ if (globalPos[0] > this->gridGeometry().bBoxMax()[0] - 1e-10)
+ return { rho_*outerFlowRate_/(2*M_PI*outerRadius_) };
+ else
+ return { rho_*innerFlowRate_/(2*M_PI*innerRadius_) };
+ }
+
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
+ { return exactSolution(globalPos); }
+
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
+ {
+ if (useStationaryAnalyticInitialSolution_)
+ return exactSolution(globalPos);
+ else
+ {
+ PrimaryVariables values(initialPressure_);
+ values.setState(Indices::bothPhases);
+ return values;
+ }
+ }
+
+ /*! \brief Analytical solution
+
+ Analytical solution of stationary incompressible
+ Richards equations with source term B:
+
+ ψ = -µ/k A/(2π) ln(r) + r^2/4 µ/k B + C
+
+ Three constants to fix with three boundary conditions:
+ - flux at inner boundary (r=r_i): q_i = 2πr_i k/µ ∂ψ/∂r
+ - flux at outer boundary (r=r_o): q_o = -2πr_o k/µ ∂ψ/∂r
+ - ψ = T(p) at inner boundary -> ψ = ψ_ri
+
+ Yields the solution in terms of q_i, q_o, ψ_ri:
+
+ ψ = ψ_ri - µ/k A/(2π) ln(r/r_i) + (r^2 - r_i^2)/4 µ/k B
+ where
+ A = - q_i + πr_i^2 B
+ B = (q_i + q_o)/(π(r_i^2 - r_o^2)) [source]
+ */
+ PrimaryVariables exactSolution(const GlobalPosition& globalPos, const Scalar innerPsi) const
+ {
+ const auto r = globalPos[0];
+ const auto ri2 = innerRadius_*innerRadius_;
+ const auto A = -innerFlowRate_ + M_PI*ri2*source_;
+ const auto psi = innerPsi
+ - A*mu_/(2*M_PI*k_)*std::log(r/innerRadius_)
+ + source_*0.25*(r*r - ri2)*mu_/k_;
+
+ PrimaryVariables priVars(applyInverseKirchhoffTransformation_(psi));
+ priVars.setState(Indices::bothPhases);
+ return priVars;
+ }
+
+ PrimaryVariables exactSolution(const GlobalPosition& globalPos) const
+ { return exactSolution(globalPos, innerPsi_); }
+
+ Scalar nonwettingReferencePressure() const
+ { return 1.0e5; };
+
+ template
+ Scalar computeTime(const SolutionVector& p, const GridVariables& vars, const Scalar initialSaturation) const
+ {
+ const auto& gg = this->gridGeometry();
+ auto fvGeometry = localView(gg);
+ auto elemVolVars = localView(vars.curGridVolVars());
+ Scalar totalWaterVolume = 0.0;
+ Scalar refWaterVolume = 0.0;
+ using Extrusion = Extrusion_t;
+ for (const auto& element : elements(gg.gridView()))
+ {
+ fvGeometry.bindElement(element);
+ elemVolVars.bindElement(element, fvGeometry, p);
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ totalWaterVolume += Extrusion::volume(scv)*volVars.porosity()*volVars.saturation(0);
+ refWaterVolume += Extrusion::volume(scv)*volVars.porosity()*initialSaturation;
+ }
+ }
+
+ return (refWaterVolume - totalWaterVolume)/(innerFlowRate_ + outerFlowRate_);
+ }
+
+ template
+ void updateStorageDerivative(const SolutionVector& p, const SolutionVector& pOld,
+ const GridVariables& vars, const Scalar dt,
+ std::vector& storageDerivative) const
+ {
+ const auto& gg = this->gridGeometry();
+ auto fvGeometry = localView(gg);
+ auto elemVolVars = localView(vars.curGridVolVars());
+ auto oldElemVolVars = localView(vars.prevGridVolVars());
+ storageDerivative.assign(storageDerivative.size(), 0.0);
+ auto volumes = storageDerivative;
+ using Extrusion = Extrusion_t;
+ for (const auto& element : elements(gg.gridView()))
+ {
+ fvGeometry.bindElement(element);
+ elemVolVars.bindElement(element, fvGeometry, p);
+ oldElemVolVars.bindElement(element, fvGeometry, pOld);
+ for (const auto& scv : scvs(fvGeometry))
+ {
+ const auto& volVars = elemVolVars[scv];
+ const auto& oldVolVars = oldElemVolVars[scv];
+ storageDerivative[scv.dofIndex()] += Extrusion::volume(scv)*(volVars.saturation(0) - oldVolVars.saturation(0))/dt;
+ volumes[scv.dofIndex()] += Extrusion::volume(scv);
+ }
+ }
+
+ for (int i = 0; i < storageDerivative.size(); ++i)
+ storageDerivative[i] /= volumes[i];
+ }
+
+ // return saturation provided pressure in Pa
+ const Scalar saturation(Scalar const pw) const
+ {
+ const auto& pcKrSwCurve = this->spatialParams().pcKrSwCurve();
+ const Scalar pc = this->nonwettingReferencePressure() - pw;
+ return pcKrSwCurve.sw(pc);
+ }
+
+ /*!
+ * \brief convert pressure head in cm to pressure in Pa
+ */
+ Scalar headToPressure(const Scalar head) const
+ {
+ return this->nonwettingReferencePressure() + head / 100.0 * rho_ * 9.81;
+ }
+
+ /*!
+ * \brief convert pressure in Pa to potential psi in Pa (Kirchhoff transformation)
+ */
+ Scalar pressureToPsi(const Scalar pw) const
+ {
+ return applyKirchhoffTransformation_(pw);
+ }
+
+ /*!
+ * \brief Disable the source term (for the instationary problem)
+ */
+ void disableSource()
+ { enableSource_ = false; }
+
+ /*!
+ * \brief Disable using the stationary initial solution (for the instationary problem)
+ */
+ void disableStationaryInitialSolution()
+ { useStationaryAnalyticInitialSolution_ = false; }
+
+ /*!
+ * \brief Enable using Neumann boundary conditions on the inner boundary
+ */
+ void enableInnerNeumannBC()
+ { enableInnerNeumann_ = true; }
+
+private:
+
+ /*!
+ * \brief relative wetting phase permeability
+ */
+ Scalar kr_(Scalar pw) const
+ {
+ const auto& pcKrSwCurve = this->spatialParams().pcKrSwCurve();
+ const Scalar pc = this->nonwettingReferencePressure() - pw;
+ const Scalar sw = pcKrSwCurve.sw(pc);
+ return pcKrSwCurve.krw(sw);
+ }
+
+ /*!
+ * \brief transformed variable p->psi
+ */
+ Scalar applyKirchhoffTransformation_(Scalar pw) const
+ {
+ return integrateScalarFunction(
+ [&](const double p){ return kr_(p); }, -1.5e6, pw, 1e-20
+ );
+ }
+
+ /*!
+ * \brief inverse transformation psi->p
+ */
+ Scalar applyInverseKirchhoffTransformation_(Scalar psi) const
+ {
+ const auto residual = [&, psi](const auto& p){
+ return psi - applyKirchhoffTransformation_(p);
+ };
+
+ return findScalarRootBrent(-1.5e6, 1.0e6, residual, 1e-14, 200000);
+ }
+
+ Scalar k_, mu_, rho_;
+ Scalar innerRadius_, outerRadius_, source_;
+ Scalar innerFlowRate_, outerFlowRate_, innerPsi_;
+ Scalar initialPressure_;
+ bool enableOuterNeumann_, enableInnerNeumann_;
+ bool useStationaryAnalyticInitialSolution_ = true;
+ bool enableSource_ = true;
+};
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/richards/annulus/properties.hh b/test/porousmediumflow/richards/annulus/properties.hh
new file mode 100644
index 0000000000000000000000000000000000000000..3d1fc01ac1281083aaddaf0453b5d40d1f28c08c
--- /dev/null
+++ b/test/porousmediumflow/richards/annulus/properties.hh
@@ -0,0 +1,79 @@
+// -*- 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 . *
+ *****************************************************************************/
+
+#ifndef DUMUX_RICHARDS_ANNULUS_PROPERTIES_HH
+#define DUMUX_RICHARDS_ANNULUS_PROPERTIES_HH
+
+#include
+
+#include
+#include
+#include
+
+#include "problem.hh"
+#include "spatialparams.hh"
+
+namespace Dumux::Properties {
+
+namespace TTag {
+struct RichardsAnnulus { using InheritsFrom = std::tuple; };
+} // end namespace TTag
+
+template
+struct Grid
+{ using type = Dune::YaspGrid<1, Dune::TensorProductCoordinates, 1>>; };
+
+template
+struct Problem
+{
+ using type = RichardsAnnulusProblem<
+ TypeTag, GetPropType
+ >;
+};
+
+template
+struct SpatialParams
+{
+ using type = RichardsAnnulusSpatialParams<
+ GetPropType, GetPropType
+ >;
+};
+
+template
+struct GridGeometry
+{
+private:
+ static constexpr bool enableCache = getPropValue();
+ using Scalar = GetPropType;
+ using GridView = typename GetPropType::LeafGridView;
+
+ // We take the default traits as basis and exchange the extrusion type
+ // The first axis (x-axis) is the radial axis, hence the zero. That means we rotate about the second axis (y-axis).
+ struct GGTraits : public BoxDefaultGridGeometryTraits
+ { using Extrusion = RotationalExtrusion<0>; };
+
+public:
+ // Pass the above traits to the box grid geometry such that it uses the
+ // rotation-symmetric sub-control volumes and faces.
+ using type = BoxFVGridGeometry;
+};
+
+} // end namespace Dumux::Properties
+
+#endif
diff --git a/test/porousmediumflow/richards/annulus/spatialparams.hh b/test/porousmediumflow/richards/annulus/spatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..d81487ce2abd352f4562d031072f987c9c589c04
--- /dev/null
+++ b/test/porousmediumflow/richards/annulus/spatialparams.hh
@@ -0,0 +1,95 @@
+// -*- 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 . *
+ *****************************************************************************/
+
+#ifndef DUMUX_RICHARDS_ANNULUS_SPATIAL_PARAMETERS_HH
+#define DUMUX_RICHARDS_ANNULUS_SPATIAL_PARAMETERS_HH
+
+#include
+#include
+
+#include
+#include
+
+namespace Dumux {
+
+/*!
+ * \ingroup RichardsTests
+ * \brief Spatial parameters for the Richards benchmarks
+ */
+template
+class RichardsAnnulusSpatialParams
+: public FVPorousMediumFlowSpatialParamsMP>
+{
+ using ThisType = RichardsAnnulusSpatialParams;
+ using ParentType = FVPorousMediumFlowSpatialParamsMP;
+ using GridView = typename GridGeometry::GridView;
+ using GlobalPosition = typename GridView::template Codim<0>::Geometry::GlobalCoordinate;
+ using PcKrSwCurve = FluidMatrix::VanGenuchtenDefault;
+
+public:
+ // export permeability type
+ using PermeabilityType = Scalar;
+
+ RichardsAnnulusSpatialParams(std::shared_ptr gridGeometry)
+ : ParentType(gridGeometry)
+ , paramGroup_("SpatialParams." + getParam("SpatialParams.SoilType"))
+ , pcKrSwCurve_(paramGroup_)
+ , permeability_(getParam(paramGroup_ + ".Permeability"))
+ , porosity_(getParam(paramGroup_ + ".Porosity"))
+ {
+ std::cout << "Using " << getParam(paramGroup_ + ".Name") << " parameters" << std::endl;
+ }
+
+ /*!
+ * \brief Returns the intrinsic permeability tensor [m^2] at a given location
+ */
+ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const
+ { return permeability_; }
+
+ /*!
+ * \brief Returns the porosity [-] at a given location
+ */
+ Scalar porosityAtPos(const GlobalPosition& globalPos) const
+ { return porosity_; }
+
+ /*!
+ * \brief Returns the fluid-matrix interaction law for the sub-control volume
+ */
+ auto fluidMatrixInteractionAtPos(const GlobalPosition& globalPos) const
+ { return makeFluidMatrixInteraction(pcKrSwCurve_); }
+
+ /*!
+ * \brief material law for analytic solution
+ */
+ const PcKrSwCurve& pcKrSwCurve() const { return pcKrSwCurve_; }
+
+ /*!
+ * \brief the parameter group selecting the soil type
+ */
+ const std::string& paramGroup() const { return paramGroup_; }
+
+private:
+ std::string paramGroup_;
+ const PcKrSwCurve pcKrSwCurve_;
+ const Scalar permeability_, porosity_;
+};
+
+} // end namespace Dumux
+
+#endif