Commit 8732166a authored by Martin Utz's avatar Martin Utz Committed by Timo Koch
Browse files

[examples][shallowwater] Improve example structure

parent d76f4101
......@@ -29,7 +29,7 @@ You learn how to
### [:open_file_folder: Example 3: Shallow water model](shallowwaterfriction/README.md)
The shallow water flow model is applied to simulate steady subcritical flow in a river including a bottom friction model.
The shallow water flow model is applied to simulate steady subcritical flow in a channel including a bottom friction model.
You learn how to
* solve a shallow water flow problem including bottom friction
......
{
"README.md" : [
"doc/intro.md",
"doc/_intro.md"
],
"doc/swe.md" : [
"doc/swe_intro.md",
"properties.hh",
"spatialparams.hh",
"problem.hh",
"spatialparams.hh"
],
"doc/main.md" : [
"doc/main_intro.md",
"main.cc"
]
],
"navigation" : {
"mainpage" : "README.md",
"subpages" : [
"doc/swe.md",
"doc/main.md"
],
"subtitles" : [
"Shallow water flow simulation setup",
"Main program flow"
]
}
}
This diff is collapsed.
# Shallow water flow with bottom friction
This example shows how the shallow water flow model can be
applied to simulate steady subcritical flow including
bottom friction (bed shear stress).
In this example, the shallow water flow model is applied to simulate
a steady subcritical flow including bottom friction (bed shear stress).
__You will learn how to__
* solve a shallow water flow problem including bottom friction
* compute and output (VTK) an analytical reference solution
__Result__. The numerical and analytical solutions for the problem will look like this:
__Result__. The numerical and analytical solutions for the free surface will look like this:
![Result Logo](img/result.png)
<figure>
<center>
<img src="img/swe_result.png" alt="Shallow water result" width="60%"/>
<figcaption> <b> Fig.1 </b> - Setup and result for the shallow water problem with bottom friction.</figcaption>
</center>
</figure>
__Table of contents__. This description is structured as follows:
[[_TOC_]]
## Mathematical model
The 2D shallow water equations (SWEs) are given by
## Problem set-up
### Model domain
The model domain is given by a rough channel with a slope of 0.001.
The domain is 500 meters long and 5 meters wide.
The bottom altitude is 10 m at the inflow and hence 9.5 m at the outflow.
Bottom friction is considered by applying
Manning's law ($`n`$ = 0.025).
### Boundary conditions
At the lateral sides a no-flow boundary condition is applied. Also no friction is
considered there and therefore a no slip boundary
condition is applied. These are the default boundary condition for the shallow
water model. At the left border a discharge boundary condition
is applied as inflow boundary condition with $`q = -1.0 m^2 s^{-1}`$.
At the right border a fixed water depth boundary condition
is applied for the outflow. Normal flow is assumed, therefore the water
depth at the right border is calculated using the equation
of Gauckler, Manning and Strickler.
### Initial conditons
The initial water depth is set to 1 m, which is slightly higher than the normal flow
water depth (0.87 m). Therefore, we expect a decreasing
water level during the simulation until the normal flow condition is reached in
the entire model domain. The inital velocity is set to zero.
## Model description
As mentioned above, this examples uses the shallow water equations (SWEs) to solve the problem.
These are a depth averaged simplification of the Navier-Stokes equations. To calculate the
bottom friction Manning's law is used. An alternative is Nikuradse's law, which is also implemented
in DuMu<sup>x</sup>.
### Shallow water model
The shallow water equations are given as:
```math
\frac{\partial \mathbf{U}}{\partial t} +
......@@ -33,54 +68,61 @@ where $`\mathbf{U}`$, $`\mathbf{F}`$ and $`\mathbf{G}`$ defined as
\mathbf{G} = \begin{bmatrix} hv \\ huv \\ hv^2 + \frac{1}{2} gh^2 \end{bmatrix}
```
$`Z`$ is the bedSurface, $`h`$ the water depth, $`u`$ the velocity in
x-direction and $`v`$ the velocity in y-direction, $`g`$ is the constant of gravity.
$`h`$ the water depth, $`u`$ the velocity in x-direction and $`v`$ the velocity in y-direction,
$`g`$ is the constant of gravity.
The source terms for the bed friction $`\mathbf{S_b}`$ and bed slope
The source terms for the bed slope $`\mathbf{S_b}`$ and friction
$`\mathbf{S_f}`$ are given as
```math
\mathbf{S_b} = \begin{bmatrix} 0 \\ -gh \frac{\partial z}{\partial x}
\\ -gh \frac{\partial z}{\partial y}\end{bmatrix},
\mathbf{S_f} = \begin{bmatrix} 0 \\ -ghS_{fx} \\ -ghS_{fy}\end{bmatrix}.
\mathbf{S_f} = \begin{bmatrix} 0 \\ghS_{fx} \\ghS_{fy}\end{bmatrix}.
```
For this example, a cell-centered finite volume method (`cctpfa`) is applied to solve the SWEs
in combination with a fully-implicit time discretization. For cases where no sharp fronts or
traveling waves occur it is possible to apply time steps larger than CFL number = 1 to reduce
the computation time. Even if a steady state solution is considered, an implicit time stepping method
is applied.
with the bedSurface $`z`$. $`S_{fx}`$ and $`S_{fy}`$ are the bed shear stess
components in x- and y-direction, which are calculated by Manning's law.
## Problem set-up
### Mannings law
The empirical Manning model specifies the bed shear stress by the following equations:
The model domain is given by a rough channel with a slope of 0.001.
The domain is 500 meters long and 10 meters wide.
![Domain](img/domain.png).
```math
S_{fx} = \frac{n^2u}{R_{hy}^{4/3}} \sqrt(u^2 + v^2),
Bottom friction is considered by applying
the friction law of Manning (Manning n = 0.025). At the lateral sides no friction is considered and a
no-flow no slip boundary condition is applied. This is the default boundary condition for the shallow water model.
S_{fy} = \frac{n^2v}{R_{hy}^{4/3}} \sqrt(u^2 + v^2)
```
$`n`$ is Manning's friction value and $`R_{hy}`$ is the hydraulic radius,
which is assumed to be equal to the water depth $`h`$.
At the left border a discharge boundary condition
is applied as inflow boundary condition with q = -1.0 ($`m^2 s^{-1}`$). At the right border a water fixed depth boundary condition
is applied for the outflow. Normal flow is assumed, therefore the water depth at the right border is calculated after
the of Gaukler-Manning-Strickler equation:
### Analytical solution
Since normal flow conditions are assumed, the analytic solution is calculated using the equation
of Gauckler, Manning and Strickler:
```math
v_m = n^{-1} R_{hy}^{2/3} I_s^{1/2}
```
Where the mean velocity $`v_m`$ is given as $`v_m = \frac{q}{h}`$,
$`n`$ is the friction value after Manning. $`R_{hy}`$ the hydraulic radius, which is assumed to be equal to
the water depth. $`I_s`$ is the bed slope and $`q`$ the unity inflow discharge
Where the mean velocity $`v_m`$ is given as
```math
v_m = \frac{q}{h}
```
$`I_s`$ is the bed slope and $`q`$ the unity inflow discharge.
Hence, the water depth $`h`$ can be calculated by
The water depth h can be calculated as
```math
h = \left(\frac{n q}{\sqrt{I_s}} \right)^{3/5}
```
The formula of Gaukler Manning and Strickler is also used to calculate the analytic solution. All parameters
for the simulation are given in the file `params.input`.
### Discretisation
For this example, a cell-centered finite volume method (cctpfa) is applied to solve the SWEs
in combination with a fully-implicit time discretization. For cases where no sharp fronts or
traveling waves occur it is possible to apply time steps larger than CFL number = 1 to reduce
the computation time. Even if a steady state solution is considered, an implicit time stepping method
is applied.
# Implementation
......
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](swe.md) |
|---|---:|
# Part 2: Main program flow
We want to solve a shallow water flow problem in a rough
channel to obtain the water table and
compare it to the analytic solution. This is done with the
`main()` function
of the program which is defined in the file `main.cc` described below.
The code documentation is structured as follows:
[[_TOC_]]
## The main file `main.cc`
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
### Included header files
<details><summary> Click to show includes</summary>
These are DUNE helper classes related to parallel computations, time measurements and file I/O
```cpp
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
```
The following headers include functionality related to property definition or retrieval, as well as
the retrieval of input parameters specified in the input file or via the command line.
```cpp
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
```
The following files contains the available linear solver backends, the non linear Newton Solver
and the assembler for the linear systems arising from finite volume discretizations
(box-scheme, tpfa-approximation, mpfa-approximation).
```cpp
#include <dumux/linear/amgbackend.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/assembly/fvassembler.hh>
```
The following class provides a convenient way of writing of dumux simulation results to VTK format.
```cpp
#include <dumux/io/vtkoutputmodule.hh>
```
The gridmanager constructs a grid from the information in the input or grid file.
Many different Dune grid implementations are supported, of which a list can be found
in `gridmanager.hh`.
```cpp
#include <dumux/io/grid/gridmanager.hh>
```
We include the header file specifing the properties of this example
```cpp
#include "properties.hh"
```
</details>
### The main function
We will now discuss the main program flow implemented within the `main` function.
At the beginning of each program using Dune, an instance of `Dune::MPIHelper` has to
be created. Moreover, we parse the run-time arguments from the command line and the
input file:
```cpp
int main(int argc, char** argv) try
{
using namespace Dumux;
// The Dune MPIHelper must be instantiated for each program using Dune
Dune::MPIHelper::instance(argc, argv);
// We parse command line arguments and input file
Parameters::init(argc, argv);
```
We define a convenience alias for the type tag of the problem. The type
tag contains all the properties that are needed to define the model and the problem
setup. Throughout the main file, we will obtain types defined for these type tag
using the property system, i.e. with `GetPropType`.
```cpp
using TypeTag = Properties::TTag::RoughChannel;
```
#### Step 1: Create the grid
The `GridManager` class creates the grid from information given in the input file.
This can either be a grid file or, in the case of structured grids, one can specify the coordinates
of the corners of the grid and the number of cells to be used to discretize each spatial direction.
```cpp
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
// We compute on the leaf grid view
const auto& leafGridView = gridManager.grid().leafGridView();
```
#### Step 2: Solving the shallow water problem
First, a finite volume grid geometry is constructed from the grid that was created above.
This builds the sub-control volumes (scv) and sub-control volume faces (scvf) for each element
of the grid partition.
```cpp
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
gridGeometry->update();
```
We now instantiate the problem, in which we define the boundary and initial conditions.
```cpp
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
```
We initialize the solution vector. The shallow water problem is transient,
therefore the initial solution defined in the problem is applied to the
solution vector. On the basis of this solution, we initialize then the grid variables.
```cpp
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x(gridGeometry->numDofs());
problem->applyInitialSolution(x);
auto xOld = x;
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
```
Let us now instantiate the time loop. Therefore, we read in some time loop parameters from the input file.
The parameter `tEnd` defines the duration of the simulation, `dt` the initial time step size and `maxDt` the maximal time step size.
Moreover, we define the end of the simulation `tEnd` as check point in the time loop at which we will write the solution to vtk files.
```cpp
using Scalar = GetPropType<TypeTag, Properties::Scalar>; // type for scalar values
const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
auto dt = getParam<Scalar>("TimeLoop.DtInitial");
// We instantiate time loop.
auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
timeLoop->setMaxTimeStepSize(maxDt);
timeLoop->setCheckPoint(tEnd);
```
We initialize the assembler with a time loop for the transient problem.
Within the time loop, we will use this assembler in each time step to assemble the linear system.
```cpp
using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
```
We initialize the linear solver.
```cpp
using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
```
We initialize the non-linear solver.
```cpp
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
```
The following lines of code initialize the vtk output module, add the velocity output facility
and write out the initial solution. At each checkpoint, we will use the output module to write
the solution of a time step into a corresponding vtk file.
```cpp
VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables,x, problem->name());
// add model-specific output fields to the writer
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
IOFields::initOutputModule(vtkWriter);
// We add the analytical solution ("exactWaterDepth" and "exactVelocityX") to the predefined specific output.
vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth");
vtkWriter.addField(problem->getExactVelocityX(), "exactVelocityX");
// We calculate the analytic solution.
problem->analyticalSolution();
// write initial solution (including the above calculated analytical solution.
vtkWriter.write(0.0);
```
##### The time loop
We start the time loop and solve a new time step as long as `tEnd` is not reached. In every time step,
the problem is assembled and solved, the solution is updated, and when a checkpoint is reached the solution
is written to a new vtk file. In addition, statistics related to CPU time, the current simulation time
and the time step sizes used is printed to the terminal.
```cpp
timeLoop->start(); do
{
// First we define the old solution as the solution of the previous time step for storage evaluations.
assembler->setPreviousSolution(xOld);
// We solve the non-linear system with time step control, using Newthon's method.
nonLinearSolver.solve(x,*timeLoop);
// We make the new solution the old solution.
xOld = x;
gridVariables->advanceTimeStep();
// We advance to the time loop to the next step.
timeLoop->advanceTimeStep();
// We write vtk output, if we reached the check point (end of time loop)
if (timeLoop->isCheckPoint())
vtkWriter.write(timeLoop->time());
// We report statistics of this time step.
timeLoop->reportTimeStep();
// We set new dt as suggested by newton controller for the next time step.
timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
} while (!timeLoop->finished());
```
The following piece of code prints a final status report of the time loop
before the program is terminated.
```cpp
timeLoop->finalize(leafGridView.comm());
return 0;
}
```
#### Exception handling
In this part of the main file we catch and print possible exceptions that could
occur during the simulation.
<details><summary> Click to show exception handler</summary>
```cpp
// errors related to run-time parameters
catch (const Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
// errors related to the parsing of Dune grid files
catch (const Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
"). Most likely, the DGF file name is wrong "
"or the DGF file is corrupted, "
"e.g. missing hash at end of file or wrong number (dimensions) of entries."
<< " ---> Abort!" << std::endl;
return 2;
}
// generic error handling with Dune::Exception
catch (const Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
}
// other exceptions
catch (...)
{
std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
return 4;
}
```
</details>
</details>
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](swe.md) |
|---|---:|
# Part 2: Main program flow
We want to solve a shallow water flow problem in a rough
channel to obtain the water table and
compare it to the analytic solution. This is done with the
`main()` function
of the program which is defined in the file `main.cc` described below.
The code documentation is structured as follows:
[[_TOC_]]
This diff is collapsed.
# Part 1: Implementation of the shallow water flow simulation setup
The shallow water flow setup, including the bottom friction,
is implemented in the files `properties.hh`,
`problem.hh` and `spatialparams.hh`. In the first of these files, a new
type tag is declared for this problem. This allows the specialization
of DuMu<sup>x</sup> `properties` for this type tag, which can be used to customize compile-time
settings for the simulation. Two exemplary `properties`, that are mandatory to be
specialized, are `Problem` and `SpatialParams`. With the first, one sets the
`Problem` class to be used, in which users can define initial and boundary conditions.
Similarly, in the `SpatialParams` class one implements the parameter distributions
(e.g. friction value) that should be used by the model.
The documentation provided in the sequel is structured as follows:
[[_TOC_]]
......@@ -17,138 +17,168 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
// ## The file `main.cc`
// ## The main file `main.cc`
// [[content]]
//
//
// This is the main file for the shallow water example. Here we can see the programme sequence and how the system is solved using newton's method.
// ### Includes
// ### Included header files
// [[details]] includes
// [[exclude]]
// Some generic includes.
#include <config.h>
// Standard header file for C++, to get time and date information.
#include <ctime>
// Standard header file for C++, for in- and output.
#include <iostream>
// Dumux is based on DUNE, the Distributed and Unified Numerics Environment, which provides several grid managers and linear solvers. So we need some includes from that.
// [[/exclude]]
//
// These are DUNE helper classes related to parallel computations, time measurements and file I/O
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/timer.hh>
#include <dune/grid/io/file/dgfparser/dgfexception.hh>
#include <dune/grid/io/file/vtk.hh>
// We need the following class to simplify the writing of dumux simulation data to VTK format.
#include <dumux/io/vtkoutputmodule.hh>
// In Dumux a property system is used to specify the model. For this, different properties are defined containing type definitions, values and methods. All properties are declared in the file properties.hh.
// The following headers include functionality related to property definition or retrieval, as well as
// the retrieval of input parameters specified in the input file or via the command line.
#include <dumux/common/properties.hh>
// The following file contains the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.
#include <dumux/common/parameters.hh>
// The file dumuxmessage.hh contains the class defining the start and end message of the simulation.
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
// The gridmanager constructs a grid from the information in the input or grid file. There is a specification for the different supported grid managers.
#include <dumux/io/grid/gridmanager.hh>
// We include the linear solver to be used to solve the linear system
// The following files contains the available linear solver backends, the non linear Newton Solver
// and the assembler for the linear systems arising from finite volume discretizations
// (box-scheme, tpfa-approximation, mpfa-approximation).
#include <dumux/linear/amgbackend.hh>
#include <dumux/linear/linearsolvertraits.hh>
// We include the nonlinear newtons method
#include <dumux/nonlinear/newtonsolver.hh>
// Further we include assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation)
#include <dumux/assembly/fvassembler.hh>
// We include the properties
// The following class provides a convenient way of writing of dumux simulation results to VTK format.
#include <dumux/io/vtkoutputmodule.hh>
// The gridmanager constructs a grid from the information in the input or grid file.
// Many different Dune grid implementations are supported, of which a list can be found
// in `gridmanager.hh`.
#include <dumux/io/grid/gridmanager.hh>
// We include the header file specifing the properties of this example
#include "properties.hh"
// [[/details]]
//
// ### Beginning of the main function
// ### The main function
// We will now discuss the main program flow implemented within the `main` function.
// At the beginning of each program using Dune, an instance of `Dune::MPIHelper` has to
// be created. Moreover, we parse the run-time arguments from the command line and the
// input file:
// [[codeblock]]