From babac4bbfbb2d941c132818a43400aaf89c60a18 Mon Sep 17 00:00:00 2001 From: Timo Koch Date: Wed, 4 May 2022 19:55:23 +0200 Subject: [PATCH] [benchmark][richards] Add benchmark case on annulus with analytic solution --- test/porousmediumflow/richards/CMakeLists.txt | 1 + .../richards/annulus/CMakeLists.txt | 25 ++ .../richards/annulus/README.md | 29 ++ .../porousmediumflow/richards/annulus/main.cc | 181 ++++++++++ .../richards/annulus/main_instationary.cc | 204 +++++++++++ .../richards/annulus/params.input | 62 ++++ .../richards/annulus/problem.hh | 323 ++++++++++++++++++ .../richards/annulus/properties.hh | 79 +++++ .../richards/annulus/spatialparams.hh | 95 ++++++ 9 files changed, 999 insertions(+) create mode 100644 test/porousmediumflow/richards/annulus/CMakeLists.txt create mode 100644 test/porousmediumflow/richards/annulus/README.md create mode 100644 test/porousmediumflow/richards/annulus/main.cc create mode 100644 test/porousmediumflow/richards/annulus/main_instationary.cc create mode 100644 test/porousmediumflow/richards/annulus/params.input create mode 100644 test/porousmediumflow/richards/annulus/problem.hh create mode 100644 test/porousmediumflow/richards/annulus/properties.hh create mode 100644 test/porousmediumflow/richards/annulus/spatialparams.hh diff --git a/test/porousmediumflow/richards/CMakeLists.txt b/test/porousmediumflow/richards/CMakeLists.txt index 6597b11e79..44b49bdddd 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 0000000000..85be915d7b --- /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 0000000000..e9ac5cafaf --- /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 0000000000..4ca5cc906b --- /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 0000000000..d319b11a73 --- /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 0000000000..afe39435cb --- /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 0000000000..4943b6c5f9 --- /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 0000000000..3d1fc01ac1 --- /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 0000000000..d81487ce2a --- /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 -- GitLab