Commit 44d5851c authored by Hanchuan Wu's avatar Hanchuan Wu Committed by Dennis Gläser
Browse files

[examples] Add lid driven cavity example

Upgrade the tests into a documented example with
verfication benchmark against published reference data.
parent cd8f5a65
......@@ -4,3 +4,4 @@ add_subdirectory(biomineralization)
add_subdirectory(shallowwaterfriction)
add_subdirectory(freeflowchannel)
add_subdirectory(1protationsymmetry)
add_subdirectory(liddrivencavity)
......@@ -112,3 +112,20 @@ You learn how to
<figure><img src="biomineralization/img/pore_scale_w_processes_named.png" alt="biomin result"/></figure></td>
</a></td>
</tr></table>
### [:open_file_folder: Example 7: Lid-driven cavity](liddrivencavity/README.md)
<table><tr><td>
In this example, we simulate laminar incompressible flow in a cavity with the Navier-Stokes equations.
You learn how to
* solve a single-phase Navier-Stokes flow problem
* compare the results of Stokes flow (Re = 1) and Navier-Stokes flow (Re = 1000)
* compare the numerical results with the reference data using the plotting library `matplotlib`
</td>
<td width="20%"><a href="liddrivencavity/README.md">
<figure><img src="liddrivencavity/img/setup.png" alt="liddriven result"/></figure></td>
</a></td>
</tr></table>
{
"README.md" : [
"doc/_intro.md"
],
"doc/problem.md" : [
"doc/problem_intro.md",
"properties.hh",
"problem.hh",
"main.cc"
],
"doc/postprocessing.md" : [
"doc/postprocessing_text.md"
],
"navigation" : {
"mainpage" : "README.md",
"subpages" : [
"doc/problem.md",
"doc/postprocessing.md"
],
"subtitles" : [
"Implementation",
"Post-processing"
]
}
}
add_subdirectory(reference_data)
dumux_add_test(NAME example_ff_liddrivencavity
SOURCES main.cc
LABELS freeflow navierstokes example
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_closedsystem_ldc_re1-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/example_ff_liddrivencavity_re1-00002.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/example_ff_liddrivencavity params_re1.input
-Problem.Name example_ff_liddrivencavity_re1
-Grid.Cells \"64 64\"")
dumux_add_test(NAME example_ff_liddrivencavity_re1000
TARGET example_ff_liddrivencavity
LABELS freeflow navierstokes example
CMAKE_GUARD HAVE_UMFPACK
COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
CMD_ARGS --script fuzzy
--files ${CMAKE_SOURCE_DIR}/test/references/test_ff_navierstokes_closedsystem_ldc_re1000-reference.vtu
${CMAKE_CURRENT_BINARY_DIR}/example_ff_liddrivencavity_re1000-00009.vtu
--command "${CMAKE_CURRENT_BINARY_DIR}/example_ff_liddrivencavity params_re1000.input
-Problem.Name example_ff_liddrivencavity_re1000
-Grid.Cells \"64 64\" -TimeLoop.TEnd 50")
dune_symlink_to_source_files(FILES "params_re1.input" "params_re1000.input" "run_and_plot.py")
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
# Shear-driven cavity flow
We use the Navier-Stokes equations to simulate laminar incompressible flow in a
cavity whose lid moves with a constant velocity u = 1 m/s.
We will verify the numerical model by comparing the simulation results with reference data published in [Ghia et al. (1982)](https://doi.org/10.1016/0021-9991(82)90058-4) and [Jurjević (1999)](https://doi.org/10.1002/(SICI)1097-0363(19991015)31:3<601::AID-FLD892>3.0.CO;2-Z).
__Results__. After simulating a few time steps, we will obtain the following velocity field for Reynolds number Re = 1 and Re = 1000:
<figure>
<center>
<img src="img/result.svg" alt="Numerical results" width="50%"/>
<figcaption> <b> Fig.1 </b> - Steady velocity field for Stokes (left) and Navier-Stokes flow for the lid-driven cavity problem.</figcaption>
</center>
</figure>
Our numerical results agree well with the reference data:
<figure>
<center>
<img src="img/lidverification.png" alt="Lid-driven cavity verification" width="90%"/>
</center>
<figcaption> <b> Fig.2 </b> - Horizontal and vertical velocity profiles at x = 0.5 m and y = 0.5 m for Re = 1 (left) and Re = 1000 (right). exp: experimental data; num: numerical data.</figcaption>
</figure>
__In this example, you will__
* solve a single-phase Navier-Stokes flow problem
* see the differences between Stokes flow (Re = 1) and Navier-Stokes flow (Re = 1000)
* compare the numerical results with the reference data using the plotting library `matplotlib`
__Table of contents__. This description is structured as follows:
[[_TOC_]]
## Problem setup
Flow in a cavity with the dimensions 1 m × 1 m is considered, where the top lid is
moving at a constant speed of 1 m/s to the right.
The following figure illustrates the setup:
<figure>
<center>
<img src="img/setup.png" alt="Lid-driven cavity setup" width="35%"/>
<figcaption> <b> Fig.3 </b> - Setup for the lid-driven cavity problem.</figcaption>
</center>
</figure>
Two different flow regimes at Re = 1 ($`\nu = 1`$ $`\text{m}^2/\text{s}`$) and Re = 1000 ($`\nu = 1000`$ $`\text{m}^2/\text{s}`$) are simulated, where the Reynolds number is defined with respect to the cavity’s side length.
## Mathematical & numerical models
Mass and momentum balance are given by
```math
\nabla \cdot \bold{v} =0,
```
```math
\frac{(\partial\rho\bold{v})}{\partial t} + \nabla \cdot (\rho\bold{v}\bold{v}^{\text{T}}) =\nabla\cdot\left[\mu\left(\nabla\bold{v}+\nabla\bold{v}^{\text{T}}\right)\right]- \nabla p,
```
where $`\bold{v}`$ and p are the velocity and pressure of the fluid (primary variables). $`\rho`$ and $`\mu=\rho\nu`$ are the mass density and dynamic viscosity (fluid properties).
All equations are discretized with the staggered-grid finite-volume scheme as spatial discretization with pressures and velocity components as primary variables. For details on the discretization scheme, we refer to the DuMu<sup>x</sup> [handbook](https://dumux.org/docs/handbook/master/dumux-handbook.pdf).
# Implementation & Postprocessing
## Part 1: Implementation
| [:arrow_right: Click to continue with part 1 of the documentation](doc/problem.md) |
|---:|
## Part 2: Post-processing
| [:arrow_right: Click to continue with part 2 of the documentation](doc/postprocessing.md) |
|---:|
# Shear-driven cavity flow
We use the Navier-Stokes equations to simulate laminar incompressible flow in a
cavity whose lid moves with a constant velocity u = 1 m/s.
We will verify the numerical model by comparing the simulation results with reference data published in [Ghia et al. (1982)](https://doi.org/10.1016/0021-9991(82)90058-4) and [Jurjević (1999)](https://doi.org/10.1002/(SICI)1097-0363(19991015)31:3<601::AID-FLD892>3.0.CO;2-Z).
__Results__. After simulating a few time steps, we will obtain the following velocity field for Reynolds number Re = 1 and Re = 1000:
<figure>
<center>
<img src="img/result.svg" alt="Numerical results" width="50%"/>
<figcaption> <b> Fig.1 </b> - Steady velocity field for Stokes (left) and Navier-Stokes flow for the lid-driven cavity problem.</figcaption>
</center>
</figure>
Our numerical results agree well with the reference data:
<figure>
<center>
<img src="img/lidverification.png" alt="Lid-driven cavity verification" width="90%"/>
</center>
<figcaption> <b> Fig.2 </b> - Horizontal and vertical velocity profiles at x = 0.5 m and y = 0.5 m for Re = 1 (left) and Re = 1000 (right). exp: experimental data; num: numerical data.</figcaption>
</figure>
__In this example, you will__
* solve a single-phase Navier-Stokes flow problem
* see the differences between Stokes flow (Re = 1) and Navier-Stokes flow (Re = 1000)
* compare the numerical results with the reference data using the plotting library `matplotlib`
__Table of contents__. This description is structured as follows:
[[_TOC_]]
## Problem setup
Flow in a cavity with the dimensions 1 m × 1 m is considered, where the top lid is
moving at a constant speed of 1 m/s to the right.
The following figure illustrates the setup:
<figure>
<center>
<img src="img/setup.png" alt="Lid-driven cavity setup" width="35%"/>
<figcaption> <b> Fig.3 </b> - Setup for the lid-driven cavity problem.</figcaption>
</center>
</figure>
Two different flow regimes at Re = 1 ($`\nu = 1`$ $`\text{m}^2/\text{s}`$) and Re = 1000 ($`\nu = 1000`$ $`\text{m}^2/\text{s}`$) are simulated, where the Reynolds number is defined with respect to the cavity’s side length.
## Mathematical & numerical models
Mass and momentum balance are given by
```math
\nabla \cdot \bold{v} =0,
```
```math
\frac{(\partial\rho\bold{v})}{\partial t} + \nabla \cdot (\rho\bold{v}\bold{v}^{\text{T}}) =\nabla\cdot\left[\mu\left(\nabla\bold{v}+\nabla\bold{v}^{\text{T}}\right)\right]- \nabla p,
```
where $`\bold{v}`$ and p are the velocity and pressure of the fluid (primary variables). $`\rho`$ and $`\mu=\rho\nu`$ are the mass density and dynamic viscosity (fluid properties).
All equations are discretized with the staggered-grid finite-volume scheme as spatial discretization with pressures and velocity components as primary variables. For details on the discretization scheme, we refer to the DuMu<sup>x</sup> [handbook](https://dumux.org/docs/handbook/master/dumux-handbook.pdf).
# Implementation & Postprocessing
<!-- 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) |
|---|---:|
# Part 2: Post processing
In this part we first visualize our simulation result and then a verification and validation of the Navier-Stokes model implemented in DuMu<sup>x</sup> is conducted.
## Step 1. Visualize results with Paraview
The result files (`.vtu` and `.pvd` files) can be opened with the software [ParaView](https://www.paraview.org/)
To obtain a visualization as shown in the introduction of this documented example, after loading the result file(s), choose `Filters`>`Common`>`Stream Tracer`.
For Re = 1 and Re = 1000, the result should look like this:
<figure>
<center>
<img src="../img/result.svg" alt="Lid-driven cavity setup" width="50%"/>
<figcaption> <b> Fig.1 </b> - Steady velocity field for Stokes (left) and Navier-Stokes flow for the lid-driven cavity problem.</figcaption>
</center>
</figure>
## Step 2. Compare our data with reference
The verification and validation are essential to guarantee the accuracy and credibility of the numerical models.
The velocity components for the velocity components at x = 0.5m and y = 0.5m are obtained as we run the test cases. We compare our results with the reference data reconstructed from [Ghia et al.](https://doi.org/10.1016/0021-9991(82)90058-4) and [Jurjević](https://doi.org/10.1002/(SICI)1097-0363(19991015)31:3<601::AID-FLD892>3.0.CO;2-Z).
For convenience, we placed the reference data in the folder named `reference_data`. For instance, the files `ghia_x.csv` and `ghia_y.csv` represent the reference vertical velocity component $`v_x`$ at x = 0.5 m and $`v_y`$ at y = 0.5 m for the scenario Re = 1000, respectively. The other files in the folder `reference_data` represent the numerical and experimental data for the scenario Re = 1.
Assuming that you have `python3` with `matplotlib` installed on your device, this comparison process can be done via the script `run_and_plot.py`. Type `python3 run_and_plot.py` in the command and you should see the following plot:
<figure>
<center>
<img src="../img/lidverification.png" alt="Lid-driven cavity verification" width="90%"/>
</center>
<figcaption> <b> Fig.2 </b> - Horizontal and vertical velocity profiles at x = 0.5 m and y = 0.5 m for Re = 1 (left) and Re = 1000 (right).</figcaption>
</figure>
It can be seen that the numerical results calculated by DuMu<sup>x</sup> are in good agreement with the reference data.
| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_left: Go back to part 1](problem.md) |
|---|---:|
# Part 2: Post processing
In this part we first visualize our simulation result and then a verification and validation of the Navier-Stokes model implemented in DuMu<sup>x</sup> is conducted.
## Step 1. Visualize results with Paraview
The result files (`.vtu` and `.pvd` files) can be opened with the software [ParaView](https://www.paraview.org/)
To obtain a visualization as shown in the introduction of this documented example, after loading the result file(s), choose `Filters`>`Common`>`Stream Tracer`.
For Re = 1 and Re = 1000, the result should look like this:
<figure>
<center>
<img src="../img/result.svg" alt="Lid-driven cavity setup" width="50%"/>
<figcaption> <b> Fig.1 </b> - Steady velocity field for Stokes (left) and Navier-Stokes flow for the lid-driven cavity problem.</figcaption>
</center>
</figure>
## Step 2. Compare our data with reference
The verification and validation are essential to guarantee the accuracy and credibility of the numerical models.
The velocity components for the velocity components at x = 0.5m and y = 0.5m are obtained as we run the test cases. We compare our results with the reference data reconstructed from [Ghia et al.](https://doi.org/10.1016/0021-9991(82)90058-4) and [Jurjević](https://doi.org/10.1002/(SICI)1097-0363(19991015)31:3<601::AID-FLD892>3.0.CO;2-Z).
For convenience, we placed the reference data in the folder named `reference_data`. For instance, the files `ghia_x.csv` and `ghia_y.csv` represent the reference vertical velocity component $`v_x`$ at x = 0.5 m and $`v_y`$ at y = 0.5 m for the scenario Re = 1000, respectively. The other files in the folder `reference_data` represent the numerical and experimental data for the scenario Re = 1.
Assuming that you have `python3` with `matplotlib` installed on your device, this comparison process can be done via the script `run_and_plot.py`. Type `python3 run_and_plot.py` in the command and you should see the following plot:
<figure>
<center>
<img src="../img/lidverification.png" alt="Lid-driven cavity verification" width="90%"/>
</center>
<figcaption> <b> Fig.2 </b> - Horizontal and vertical velocity profiles at x = 0.5 m and y = 0.5 m for Re = 1 (left) and Re = 1000 (right).</figcaption>
</figure>
It can be seen that the numerical results calculated by DuMu<sup>x</sup> are in good agreement with the reference data.
<!-- 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](postprocessing.md) |
|---|---:|
# Part 1: Implementation
The implementation of simulation setup and main flow is structured as follows:
[[_TOC_]]
## Compile-time settings (`properties.hh`)
In this file, the type tag used for this simulation is defined,
for which we then specialize properties (compile time options) to the needs of the desired setup.
<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>
The `NavierStokes` type tag specializes most of the properties required for Navier-
Stokes single-phase flow simulations in DuMuX. 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/freeflow/navierstokes/model.hh>
```
We want to use `YaspGrid`, an implementation of the dune grid interface for structured grids:
```cpp
#include <dune/grid/yaspgrid.hh>
```
In this example, we want to discretize the equations with the staggered-grid
scheme which is so far the only available option for free-flow models in DuMux:
```cpp
#include <dumux/discretization/staggered/freeflow/properties.hh>
```
The fluid properties are specified in the following headers (we use a liquid with constant properties as the fluid phase):
```cpp
#include <dumux/material/components/constant.hh>
#include <dumux/material/fluidsystems/1pliquid.hh>
```
We include the problem header used for this simulation.
```cpp
#include "problem.hh"
```
</details>
### Type tag definition
We define a type tag for our simulation with the name `LidDrivenCavityExample`
and inherit the properties specialized for the type tags `NavierStokes` and `StaggeredFreeFlowModel`.
```cpp
namespace Dumux::Properties {
// We define the `LidDrivenCavityExample` type tag and let it inherit from the single-phase `NavierStokes`
// tag (model) and the `StaggeredFreeFlowModel` (discretization scheme).
namespace TTag {
struct LidDrivenCavityExample { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
} // end namespace TTag
```
### Property specializations
In the following piece of code, mandatory properties for which no meaningful
default exist are specialized for our type tag `LidDrivenCavityExample`.
```cpp
// This sets the fluid system to be used. Here, we use a liquid with constant properties as fluid phase.
template<class TypeTag>
struct FluidSystem<TypeTag, TTag::LidDrivenCavityExample>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
};
// This sets the grid type used for the simulation. Here, we use a structured 2D grid.
template<class TypeTag>
struct Grid<TypeTag, TTag::LidDrivenCavityExample> { using type = Dune::YaspGrid<2>; };
// This sets our problem class (see problem.hh) containing initial and boundary conditions.
template<class TypeTag>
struct Problem<TypeTag, TTag::LidDrivenCavityExample> { using type = Dumux::LidDrivenCavityExampleProblem<TypeTag> ; };
```
We also set some properties related to memory management
throughout the simulation.
<details><summary> Click to show caching properties</summary>
In Dumux, one has the option to activate/deactivate the grid-wide caching of
geometries and variables. If active, the CPU time can be significantly reduced
as less dynamic memory allocation procedures are necessary. Per default, grid-wide
caching is disabled to ensure minimal memory requirements, however, in this example we
want to active all available caches, which significantly increases the memory
demand but makes the simulation faster.
```cpp
// This enables grid-wide caching of the volume variables.
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::LidDrivenCavityExample> { static constexpr bool value = true; };
//This enables grid wide caching for the flux variables.
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::LidDrivenCavityExample> { static constexpr bool value = true; };
// This enables grid-wide caching for the finite volume grid geometry
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::LidDrivenCavityExample> { static constexpr bool value = true; };
} // end namespace Dumux::Properties
```
</details>
</details>
## Initial and boundary conditions (`problem.hh`)
This file contains the __problem class__ which defines the initial and boundary
conditions for the Navier-Stokes single-phase flow simulation.
<details open>
<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](../problem.hh))</summary>
### Include files
```cpp
#include <dumux/common/properties.hh>
#include <dumux/common/parameters.hh>
```
Include the `NavierStokesProblem` class, the base
class from which we will derive.
```cpp
#include <dumux/freeflow/navierstokes/problem.hh>
```
Include the `NavierStokesBoundaryTypes` class which specifies the boundary types set in this problem.
```cpp
#include <dumux/freeflow/navierstokes/boundarytypes.hh>
```
### The problem class
As we are solving a problem related to free flow, we create a new class called `LidDrivenCavityExampleProblem`
and let it inherit from the class `NavierStokesProblem`.
```cpp
namespace Dumux {
template <class TypeTag>
class LidDrivenCavityExampleProblem : public NavierStokesProblem<TypeTag>
{
using ParentType = NavierStokesProblem<TypeTag>;
using BoundaryTypes = Dumux::NavierStokesBoundaryTypes<GetPropType<TypeTag, Properties::ModelTraits>::numEq()>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
using FVElementGeometry = typename GridGeometry::LocalView;
using SubControlVolume = typename GridGeometry::SubControlVolume;
using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using Element = typename GridGeometry::GridView::template Codim<0>::Entity;
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
// Within the constructor, we set the lid velocity to a run-time specified value.
LidDrivenCavityExampleProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
lidVelocity_ = getParam<Scalar>("Problem.LidVelocity");
}
```
#### Temperature distribution
We need to specify a constant temperature for our isothermal problem.
Fluid properties that depend on temperature will be calculated with this value.
This would be important if another fluidsystem was used.
```cpp
Scalar temperature() const
{ return 273.15 + 10; } // 10°C
```
#### Boundary conditions
With the following function we define the __type of boundary conditions__ depending on the location.
Three types of boundary conditions can be specified: Dirichlet, Neumann or outflow boundary conditions. On
Dirichlet boundaries, the values of the primary variables need to be fixed. On a Neumann boundaries,
values for derivatives need to be fixed. Outflow conditions set a gradient of zero in normal direction towards the boundary
for the respective primary variables (excluding pressure).
When Dirichlet conditions are set for the pressure, the velocity gradient
with respect to the direction normal to the boundary is automatically set to zero.
```cpp
BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
{
BoundaryTypes values;
// We set Dirichlet values for the velocity at each boundary
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
return values;
}
```
We define a function for setting a fixed Dirichlet pressure value at a given internal cell.
This is required for having a defined pressure level in our closed system domain.
```cpp
bool isDirichletCell(const Element& element,
const FVElementGeometry& fvGeometry,
const SubControlVolume& scv,
int pvIdx) const
{
auto isLowerLeftCell = [&](const SubControlVolume& scv)
{ return scv.dofIndex() == 0; };
```
We set a fixed pressure in one cell
```cpp
return (isLowerLeftCell(scv) && pvIdx == Indices::pressureIdx);
}
```
The following function specifies the __values on Dirichlet boundaries__.
We need to define values for the primary variables (velocity and pressure).
```cpp
PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values;
values[Indices::pressureIdx] = 1.1e+5;
values[Indices::velocityXIdx] = 0.0;
values[Indices::velocityYIdx] = 0.0;
// We set the no slip-condition at the top, that means the fluid has the same velocity as the lid
if (globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_)
values[Indices::velocityXIdx] = lidVelocity_;
return values;
}
```
The following function defines the initial conditions.
```cpp
PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
{
PrimaryVariables values;
values[Indices::pressureIdx] = 1.0e+5;
values[Indices::velocityXIdx] = 0.0;
values[Indices::velocityYIdx] = 0.0;
return values;
}
```
the data members of the problem class
```cpp
private:
static constexpr Scalar eps_ = 1e-6;
Scalar lidVelocity_;
};
} // end namespace Dumux
```
</details>
## 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 is DUNE helper class related to parallel computation
```cpp