Commit 34bbeecc authored by Timo Koch's avatar Timo Koch Committed by Dennis Gläser
Browse files

[examples][rotsym] add structural elements and compile docu

parent 297e729b
{
"README.md" : [
"doc/_intro.md"
],
"doc/problem.md" : [
"doc/problem_intro.md",
"properties.hh",
"problem.hh",
"spatialparams.hh"
],
"doc/main.md" : [
"doc/main_intro.md",
"main.cc"
],
"doc/paraview.md" : [
"doc/paraview_doc.md"
],
"navigation" : {
"mainpage" : "README.md",
"subpages" : [
"doc/problem.md",
"doc/main.md",
"doc/paraview.md"
],
"subtitles" : [
"Rotation-symmetric one-phase flow simulation setup",
"Main program flow",
"Post-processing with ParaView"
]
}
}
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
# Rotation-symmetric pressure distribution
__In this example, you will learn how to__
* solve a rotation-symmetric problem one-dimensionally
* perform a convergence test against an analytical solution
* apply the `Rotational Extrusion` filters in [ParaView](https://www.paraview.org/) for a two-dimensional visualization of the one-dimensional results
__Result__. With the `Rotational Extrusion` and the `Warp By Scalar` filters in [ParaView](https://www.paraview.org/),
the pressure distribution of this example looks as shown in the following picture:
<figure>
<center>
<img src="img/result.png" alt="Rotation-symmetric pressure distribution" width="60%"/>
<figcaption> <b> Fig.1 </b> - Rotation-symmetric pressure distribution on a disc (warped to 3D). </figcaption>
</center>
</figure>
__Table of contents__. This description is structured as follows:
[[_TOC_]]
## Problem setup
We consider a single-phase problem that leads to a rotation-symmetric pressure distribution.
The following figure illustrates the setup:
<figure>
<center>
<img src="img/setup.svg" alt="Rotation-symmetric setup" width="60%"/>
<figcaption> <b> Fig.2 </b> - Setup for the rotation-symmetric problem. The pressure boundary conditions are shown by the colored lines and the simulation domain is depicted in grey.</figcaption>
</center>
</figure>
This could, for example, represent a cross section of an injection/extraction well in a homogeneous
and isotropic porous medium, where the well with radius $`r_1`$ is cut out and the
injection/extraction pressure $`p_1`$ is prescribed as a Dirichlet boundary condition. At the outer
radius $`r_2`$, we set the pressure $`p_2`$. In the polar coordinates $`r`$ and $`\varphi`$, the
solution to this problem is independent of the angular coordinate $`\varphi`$ and can be reduced to
a one-dimensional problem in the radial coordinate $`r`$. Therefore, in this example, we want to
solve the problem on a one-dimensional computational domain as illustrated by the orange line in
the above figure.
## Mathematical model
In this example we are using the single-phase model of DuMu<sup>x</sup>, which considers Darcy's law to relate
the Darcy velocity $`\textbf u`$ to gradients of the pressure $`p`$. In the case of rotational
symmetry, the mass balance equation for the fluid phase can be transformed using polar coordinates:
```math
-\frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\varrho k}{\mu} \frac{\partial p}{\partial r} \right) = 0,
```
where we identify the Darcy velocity in radial direction $`u_r = -\frac{k}{\mu} \frac{\partial p}{\partial r}`$,
and where $`k`$ is the permeability of the porous medium, $`\mu`$ is the dynamic viscosity of the
fluid, and $`\varrho`$ is the fluid density.
## Discretization
We employ a finite-volume scheme to spatially discretize the mass balance equation shown above.
Let us consider a discretization of the one-dimensional domain into control volumes
$`K_i = \left[ r_i, r_{i+1} \right]`$. The discrete equation describing mass conservation inside a control volume
$`K_i`$ is obtained by integration and reads:
```math
- 2 \pi r_{i+1} \left( \varrho u_r \right)_{r_{i+1}}
+ 2 \pi r_i \left( \varrho u_r \right)_{r_i}
= 0.
```
For this type of equation, the implementation of the finite-volume schemes in DuMu<sup>x</sup> is based on
the general form:
```math
\sum_{\sigma \in \mathcal{S}_K} | \sigma | \left( \varrho \textbf u \cdot \textbf n \right)_\sigma = 0,
```
where $`\sigma`$ are the faces of the control volume and where the notation
$`( \cdot )_\sigma`$ was used to denote quantities evaluated for a face $`\sigma`$.
The area of a face is denoted with $`| \sigma |`$. Thus, comparing the two equations
we identify $`| \sigma | = 2 \pi r_\sigma`$ for the case of rotational symmetry
on a disc. Here, $`r_\sigma`$ refers to the radius at which the face is situated
in the one-dimensional discretization.
In DuMu<sup>x</sup>, this is realized in the classes `RotationSymmetricSubControlVolume` and
`RotationSymmetricSubControlVolumeFace`, which implement one-dimensional control
volumes and faces, that take into account the extrusion about the rotation axes
of symmetry in the computations of volumes and areas. This will be discussed in part 1
of the documentation.
# Implementation & Post processing
## Part 1: Rotation-symmetric one-phase flow simulation setup
| [:arrow_right: Click to continue with part 1 of the documentation](doc/problem.md) |
|---:|
## Part 2: Main program flow
| [:arrow_right: Click to continue with part 2 of the documentation](doc/main.md) |
|---:|
## Part 3: Post-processing with ParaView
| [:arrow_right: Click to continue with part 3 of the documentation](doc/paraview.md) |
|---:|
\ No newline at end of file
<!-- 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](problem.md) | [:arrow_right: Continue with part 3](paraview.md) |
|---|---|---:|
# Part 2: Main program flow
We want to solve a rotational symmetric Laplace problem on a disc and
conduct a grid convergence study against an analytical solution.
The main program flow is implemented in file `main.cc` described below.
The code documentation is structured as follows:
[[_TOC_]]
## The main program (`main.cc`)
This file contains the main program flow. In this example, we solve a stationary
and rotationally symmetric single-phase problem for a sequence of refined grids
and compute the convergence rates.
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../main.cc))</summary>
### Includes
<details><summary> Click to show includes</summary>
```cpp
#include <config.h>
#include <iostream>
#include <dune/common/parallel/mpihelper.hh>
#include <dumux/common/properties.hh> // for GetPropType
#include <dumux/common/parameters.hh> // for getParam
#include <dumux/common/integrate.hh> // for integrateL2Error
#include <dumux/linear/seqsolverbackend.hh> // for UMFPackBackend
#include <dumux/linear/pdesolver.hh> // for LinearPDESolver
#include <dumux/assembly/fvassembler.hh>
#include <dumux/assembly/diffmethod.hh>
#include <dumux/io/vtkoutputmodule.hh>
#include <dumux/io/grid/gridmanager_yasp.hh>
#include "properties.hh"
```
</details>
### Beginning of the main function
```cpp
int main(int argc, char** argv) try
{
using namespace Dumux;
// We initialize MPI. Finalization is done automatically on exit.
Dune::MPIHelper::instance(argc, argv);
// We parse the command line arguments.
Parameters::init(argc, argv);
// Convenience alias for the type tag of the problem.
using TypeTag = Properties::TTag::OnePRotSym;
```
### Create the grid and the grid geometry
```cpp
// 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);
gridGeometry->update();
```
### Initialize the problem and grid variables
```cpp
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());
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
```cpp
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
```
### Instantiate the solver
We use the `LinearPDESolver` class, which is instantiated on the basis
of an assembler and a linear solver. When the `solve` function of the
`LinearPDESolver` is called, it uses the assembler and linear
solver classes to assemble and solve the linear system around the provided
solution and stores the result therein.
```cpp
using Assembler = FVAssembler<TypeTag, DiffMethod::analytic>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
using LinearSolver = UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
LinearPDESolver<Assembler, LinearSolver> solver(assembler, linearSolver);
solver.setVerbose(false); // suppress output during solve()
```
### Solution of the problem and error computation
The problem is solved by calling `solve` on the instance of `LinearPDESolver`
that we have created above. In the following piece of code, we solve the
problem on the initial refinement and compute the corresponding L2 error.
For a convenient way of computing the L2 error, the function `integrateL2Error`
can be used.
```cpp
solver.solve(p);
// 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 = 3;
// compute initial L2 error
l2Errors[0] = integrateL2Error(*gridGeometry, p, pExact, orderQuadratureRule);
```
This procedure is now repeated for the number of refinements as specified
in the input file.
```cpp
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();
gridVariables->updateAfterGridAdaption(p);
p.resize(gridGeometry->numDofs());
updateAnalyticalSolution(pExact);
// this recreates the linear system, i.e. the sizes of
// the right hand side vector and the Jacobian matrix,
// and its sparsity pattern.
assembler->setLinearSystem();
// solve problem on refined grid
solver.solve(p);
```
#### Post-processing and output
At the end of each refinement step, the convergence
rate is printed to the terminal.
```cpp
// 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 << std::setprecision(8) << std::scientific
<< "-- L2 error for " << std::setw(5) << numDofs << " dofs: " << l2Errors[stepIdx]
<< ", rate: " << rate
<< std::endl;
}
```
After the last refinement, we write the solution to VTK file format on the
finest grid and exit the main function.
```cpp
vtkWriter.write(0.0);
// program end, return with 0 exit code (success)
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 error handler</summary>
```cpp
catch (const Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
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;
}
catch (const Dune::Exception &e)
{
std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
return 3;
}
```
</details>
</details>
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](problem.md) | [:arrow_right: Continue with part 3](paraview.md) |
|---|---|---:|
# Part 2: Main program flow
We want to solve a rotational symmetric Laplace problem on a disc and
conduct a grid convergence study against an analytical solution.
The main program flow is implemented in file `main.cc` described below.
The code documentation is structured as follows:
[[_TOC_]]
<!-- 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 2](main.md) |
|---|---:|
# Part 3: Post-processing with ParaView
The result file `example_rotationsymmetry.pvd` can be opened with the software [ParaView](https://www.paraview.org/).
To obtain a visualisation as shown in the introduction of this documented example, after loading
the result file, choose `Filters`>`Alphabetical`>`Rotational Extrusion`.
You might have to reset your view and switch to `3D`. Then apply `Filters`>`Alphabetical`>`Warp By Scalar`.
The result should look like this:
<figure>
<center>
<img src="../img/result.png" alt="Rotation-symmetric pressure distribution" width="60%"/>
<figcaption> <b> Fig.1 - </b>Rotation-symmetric pressure distribution on a disc (warped to 3D)</figcaption>
</center>
</figure>
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 2](main.md) |
|---|---:|
# Part 3: Post-processing with ParaView
The result file `example_rotationsymmetry.pvd` can be opened with the software [ParaView](https://www.paraview.org/).
To obtain a visualisation as shown in the introduction of this documented example, after loading
the result file, choose `Filters`>`Alphabetical`>`Rotational Extrusion`.
You might have to reset your view and switch to `3D`. Then apply `Filters`>`Alphabetical`>`Warp By Scalar`.
The result should look like this:
<figure>
<center>
<img src="../img/result.png" alt="Rotation-symmetric pressure distribution" width="60%"/>
<figcaption> <b> Fig.1 - </b>Rotation-symmetric pressure distribution on a disc (warped to 3D)</figcaption>
</center>
</figure>
<!-- 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_right: Continue with part 2](main.md) |
|---|---:|
# Part 1: Simulation setup
The code documentation is structured as follows:
[[_TOC_]]
## Compile-time settings (`properties.hh`)
This file defines the `TypeTag` used for the simulation in this example, for
which we specialize a number of compile-time `properties`.
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../properties.hh))</summary>
### Includes
<details><summary> Click to show includes</summary>
```cpp
#include <dune/grid/yaspgrid.hh> // for `Dune::YaspGrid`
#include <dumux/discretization/box.hh> // for `TTag::BoxModel`
```
The `OneP` type tag specializes most of the `properties` required for single-
phase flow simulations in DuMu<sup>x</sup>. We will use this in the following to inherit the
respective properties, and subsequently specialize those properties for our
type tag, which we want to modify or for which no meaningful default can be set.
```cpp
#include <dumux/porousmediumflow/1p/model.hh> // for `TTag::OneP`
```
The local residual for incompressible flow is included.
The one-phase flow model (included above) uses a default implementation of the
local residual for single-phase flow. However, in this example we are using an
incompressible fluid phase. Therefore, we are including the specialized local
residual which contains functionality to analytically compute the entries of
the Jacobian matrix. We will use this in the main file.
```cpp
#include <dumux/porousmediumflow/1p/incompressiblelocalresidual.hh>
```
We will use a single liquid phase consisting of a component with constant fluid properties.
```cpp
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
```
As mentioned at the beginning of the documentation of this example, DuMu<sup>x</sup>
provides specialized implementations of control volumes and faces for
rotation-symmetric problems. These take care of adjusting volume and area
computations to account for the extrusion about the symmetry axes.
These implementations are exported by the `RotationSymmetricGridGeometryTraits`.
```cpp
#include <dumux/discretization/rotationsymmetricgridgeometrytraits.hh>
```
The classes that define the problem and parameters used in this simulation
```cpp
#include "problem.hh"
#include "spatialparams.hh"
```
</details>
### `TypeTag` definition
A `TypeTag` for our simulation is defined, which inherits properties from the
single-phase flow model and the box scheme.
```cpp
namespace Dumux::Properties {
namespace TTag {
struct OnePRotSym { using InheritsFrom = std::tuple<OneP, BoxModel>; };
}
```
### Property specializations
In the following piece of code, mandatory `properties` for which no meaningful
default can be set, are specialized for our type tag `OnePRotSym`.
```cpp
// We use a structured 1D grid with an offset. This allows us to define the
// computational domain to be between the radii $`r_1`$ and $`r_2`$ as illustrated
// in the beginning of the documentation of this example
template<class TypeTag>
struct Grid<TypeTag, TTag::OnePRotSym>
{ using type = Dune::YaspGrid<1, Dune::EquidistantOffsetCoordinates<double, 1>>; };
// The problem class specifying initial and boundary conditions:
template<class TypeTag>
struct Problem<TypeTag, TTag::OnePRotSym>
{ using type = RotSymExampleProblem<TypeTag>; };
// Our spatial parameters class defining the permeability and porosity of the porous medium:
template<class TypeTag>
struct SpatialParams<TypeTag, TTag::OnePRotSym>
{
private:
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = RotSymExampleSpatialParams<GridGeometry, Scalar>;
};
// We use a single liquid phase consisting of a component with constant fluid properties.
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::OnePRotSym>
{
private:
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
public:
using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
};
```
As mentioned before, DuMu<sup>x</sup> provides specialized implementations of sub-control
volumes and faces for rotation-symmetric problems, which are exported by the
`RotationSymmetricGridGeometryTraits`.
Here, we pass these traits to the grid geometry of the box scheme (the scheme
that we use here) and specialize the `GridGeometry` property accordingly.
```cpp
template<class TypeTag>
struct GridGeometry<TypeTag, TTag::OnePRotSym>
{
private:
static constexpr bool enableCache = getPropValue<TypeTag, Properties::EnableGridGeometryCache>();
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using GridView = typename GetPropType<TypeTag, Properties::Grid>::LeafGridView;
// The default traits for box grid geometries
using DefaultTraits = BoxDefaultGridGeometryTraits<GridView>;
// On the basis of the default traits, define the traits for rotational symmetry.
// These will export the corresponding rotation-symmetric sub-control volumes and faces.
using GGTraits = RotationSymmetricGridGeometryTraits<DefaultTraits, RotationPolicy::disc>;
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<Scalar, GridView, enableCache, GGTraits>;
};
```
Moreover, here we use a local residual specialized for incompressible flow
that contains functionality related to analytic differentiation.