Commit babac4bb authored by Timo Koch's avatar Timo Koch Committed by Timo Koch
Browse files

[benchmark][richards] Add benchmark case on annulus with analytic solution

parent 7c0e9f09
Pipeline #17666 passed with stages
add_subdirectory(analytical)
add_subdirectory(annulus)
add_subdirectory(benchmarks)
add_subdirectory(lens)
add_subdirectory(nonisothermal)
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)
# 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%.
// -*- 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/>. *
*****************************************************************************/
#include <config.h>
#include <iostream>
#include <dune/common/exceptions.hh>
#include <dune/common/timer.hh>
#include <dumux/common/initialize.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/integrate.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/io/format.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
#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<TypeTag, Properties::Grid>;
GridManager<Grid> gridManager;
gridManager.init();
// We compute on the leaf grid view.
const auto& leafGridView = gridManager.grid().leafGridView();
// instantiate the grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
// Initialize the problem and grid variables
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(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<TypeTag, Properties::SolutionVector>;
SolutionVector p(gridGeometry->numDofs()); problem->applyInitialSolution(p);
SolutionVector pExact; updateAnalyticalSolution(pExact);
// instantiate and initialize the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(p);
// initialize VTK output
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, p, problem->name());
GetPropType<TypeTag, Properties::IOFields>::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<TypeTag, DiffMethod::analytic>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
using LinearSolver = ILU0BiCGSTABBackend;
auto linearSolver = std::make_shared<LinearSolver>();
NewtonSolver<Assembler, LinearSolver> solver(assembler, linearSolver);
// Solution of the problem and error computation
solver.solve(p);
vtkWriter.write(0.0);
const auto initialPressure = problem->headToPressure(getParam<double>("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<int>("Grid.RefinementSteps");
std::vector<double> 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;
}
// -*- 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/>. *
*****************************************************************************/
#include <config.h>
#include <iostream>
#include <dune/common/exceptions.hh>
#include <dune/common/timer.hh>
#include <dumux/common/initialize.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/integrate.hh>
#include <dumux/common/timeloop.hh>
#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/io/format.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
#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<TypeTag, Properties::Grid>;
GridManager<Grid> gridManager;
gridManager.init();
// We compute on the leaf grid view.
const auto& leafGridView = gridManager.grid().leafGridView();
// instantiate the grid geometry
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
// Initialize the problem and grid variables
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(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<TypeTag, Properties::SolutionVector>;
SolutionVector p(gridGeometry->numDofs()); problem->applyInitialSolution(p);
const auto initialPressure = problem->headToPressure(getParam<double>("Problem.InitialPressureHeadInCm", -100));
SolutionVector pAnalyticApproximation; updateAnalyticalSolution(pAnalyticApproximation, problem->pressureToPsi(initialPressure));
auto pOld = p;
// instantiate and initialize the grid variables
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(p);
auto gridVariablesAnalytic = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariablesAnalytic->init(pAnalyticApproximation);
const auto initialSaturation = problem->saturation(initialPressure);
// get some time loop parameters
const auto tEnd = getParam<double>("TimeLoop.TEnd", 24*60*60*40);
const auto maxDt = getParam<double>("TimeLoop.MaxTimeStepSize", 0.5*24*60*60);
auto dt = getParam<double>("TimeLoop.DtInitial", 0.5*24*60*60);
// instantiate time loop
auto timeLoop = std::make_shared<TimeLoop<double>>(0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
std::vector<double> storageDerivative(leafGridView.size(Grid::dimension), 0.0);
// initialize VTK output
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, p, problem->name());
GetPropType<TypeTag, Properties::IOFields>::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<TypeTag, DiffMethod::analytic>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, pOld);
using LinearSolver = ILU0BiCGSTABBackend;
auto linearSolver = std::make_shared<LinearSolver>();
NewtonSolver<Assembler, LinearSolver> 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;
}
[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
// -*- 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/>. *
*****************************************************************************/
#ifndef DUMUX_RICHARDS_ANNULUS_PROBLEM_HH
#define DUMUX_RICHARDS_ANNULUS_PROBLEM_HH
#include <iostream>
#include <cmath>
#include <dumux/common/boundarytypes.hh>
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
#include <dumux/common/numeqvector.hh>
#include <dumux/common/integrate.hh>
#include <dumux/nonlinear/findscalarroot.hh>
#include <dumux/discretization/extrusion.hh>
#include <dumux/porousmediumflow/problem.hh>
namespace Dumux {
template<class TypeTag, class FluidSystem>
class RichardsAnnulusProblem : public PorousMediumFlowProblem<TypeTag>
{
using ParentType = PorousMediumFlowProblem<TypeTag>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using NumEqVector = Dumux::NumEqVector<PrimaryVariables>;
using BoundaryTypes = Dumux::BoundaryTypes<PrimaryVariables::size()>;
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
RichardsAnnulusProblem(std::shared_ptr<const GridGeometry> 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<Scalar>(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<Scalar>("Problem.InnerFlux", 0.1); // positive means outflow (cm/day)
innerFlowRate_ = fi*2*M_PI*innerRadius_*0.01/(24*60*60);
const auto fo = getParam<Scalar>("Problem.OuterFlux", 0.0); // positive means outflow (cm/day)
outerFlowRate_ = fo*2*M_PI*outerRadius_*0.01/(24*