From 8ac386b6597e98eb2626243a6a5a67cbdc02daff Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Sat, 18 Apr 2020 10:06:12 +0200
Subject: [PATCH] Merge branch 'feature/improve-ff-example' into 'master'
Improve free-flow example
Closes #855
See merge request dumux-repositories/dumux!1976
(cherry picked from commit 2c74b446d04242ae019ebdaece1641fbb2183b59)
---
examples/freeflowchannel/.doc_config | 2 +-
examples/freeflowchannel/README.md | 523 ++++++++++++++-----------
examples/freeflowchannel/doc/_intro.md | 65 +++
examples/freeflowchannel/doc/intro.md | 45 ---
examples/freeflowchannel/main.cc | 208 +++++-----
examples/freeflowchannel/params.input | 2 +-
examples/freeflowchannel/problem.hh | 160 ++++----
examples/freeflowchannel/properties.hh | 81 +++-
8 files changed, 591 insertions(+), 495 deletions(-)
create mode 100644 examples/freeflowchannel/doc/_intro.md
delete mode 100644 examples/freeflowchannel/doc/intro.md
diff --git a/examples/freeflowchannel/.doc_config b/examples/freeflowchannel/.doc_config
index ad433016f2..6994794918 100644
--- a/examples/freeflowchannel/.doc_config
+++ b/examples/freeflowchannel/.doc_config
@@ -1,6 +1,6 @@
{
"README.md" : [
- "doc/intro.md",
+ "doc/_intro.md",
"properties.hh",
"problem.hh",
"main.cc"
diff --git a/examples/freeflowchannel/README.md b/examples/freeflowchannel/README.md
index d7221b3877..8a767204bf 100644
--- a/examples/freeflowchannel/README.md
+++ b/examples/freeflowchannel/README.md
@@ -1,10 +1,10 @@
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->
-# Freeflow through a channel
+# Free flow through a channel
__You learn how to__
-* solve a free flow channel problem
+* solve a free-flow channel problem
* set outflow boundary conditions in the free-flow context
__Results__. In this example we will obtain the following stationary velocity profile:
@@ -16,22 +16,42 @@ __Table of contents__. This description is structured as follows:
[[_TOC_]]
## Mathematical model
-The Stokes model without gravitation and without sources or sinks for a stationary, incompressible, laminar, single phase, one-component, isothermal ($`T=10^\circ C`$) flow is considered assuming a Newtonian fluid of constant density $` \varrho = 1~\frac{\text{kg}}{\text{m}^3} `$ and constant kinematic viscosity $` \nu = 1~\frac{\text{m}^2}{\text{s}} `$. The momentum balance
+In this example, the Stokes model for stationary and incompressible single phase flow is considered.
+Thus, the momentum balance equations
+
```math
- \nabla\cdot\left(\mu\left(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{\text{T}}\right)\right)+ \nabla p = 0
```
-with density $`\varrho`$, velocity $`\boldsymbol{u}`$, dynamic viscosity $`\mu=\varrho\nu`$ and pressure $`p`$ and the mass balance
+
+and the mass balance
+
```math
-\nabla \cdot \left(\varrho\boldsymbol{u}\right) =0
+\nabla \cdot \left(\boldsymbol{u}\right) =0
```
-are discretized using a staggered-grid finite-volume scheme as spatial discretization with pressures and velocity components as primary variables. For details on the discretization scheme, have a look at the Dumux [handbook](https://dumux.org/handbook).
-## Problem set-up
-This example contains a stationary free flow of a fluid through two parallel solid plates in two dimensions from left to right. The figure below shows the simulation set-up. The fluid flows into the system at the left with a constant velocity of $` v = 1~\frac{\text{m}}{\text{s}} `$. The inflow velocity profile is a block profile. Due to the no-slip, no-flow boundary conditions at the top and bottom plate, the velocity profile gradually assumes a parabolic shape along the channel. At the outlet, the pressure is fixed and a zero velocity gradient in x-direction is assumed. The physical domain, which is modeled is the rectangular domain $`x\in[0,10],~y\in[0,1]`$.
+are solved, where $`\varrho`$ and $`\mu`$ are the density and viscosity of the fluid,
+$`\boldsymbol{u}`$ is the fluid velocity and $`p`$ is the pressure. Here, we use constant fluid
+properties with $`\varrho = 1~\frac{\text{kg}}{\text{m}^3}`$ and $`\mu = 1~\text{Pa}\text{s}`$.
+Furthermore, isothermal conditions with a homogeneous temperature distribution of $`T=10^\circ C`$ are assumed.
-
+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,
+have a look at the Dumux [handbook](https://dumux.org/handbook).
-In the following, we take a close look at the files containing the set-up: At first, boundary conditions are set in `problem.hh` for the Navier-Stokes model. Afterwards, we show the different steps for solving the model in the source file `main.cc`.
+## Problem set-up
+This example considers stationary flow of a fluid between two parallel solid plates in two dimensions.
+Flow is enforced from left to right by prescribing an inflow velocity of $` v = 1~\frac{\text{m}}{\text{s}} `$
+on the left boundary, while a fixed pressure of $`p = 1.1 \text{bar}`$ and a zero velocity gradient
+in x-direction are prescribed on the right boundary. On the top and bottom boundaries, no-slip
+conditions are applied, which cause a parabolic velocity profile to develop along the channel.
+Take a look at Figure 1 for an illustration of the domain and the boundary conditions.
+
+<figure>
+ <center>
+ <img src="img/setup.png" alt="Free-flow setup" width="80%"/>
+ <figcaption> <b> Fig.1 </b> - Setup for the free flow problem.</figcaption>
+ </center>
+</figure>
# Implementation
@@ -47,95 +67,163 @@ In the following, we take a close look at the files containing the set-up: At fi
```
-## The file `properties.hh`
-In the following, we set the properties for our simulation
-(click [here](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/blob/master/slides/dumux-course-properties.pdf) for DuMux course slides on the property system).
-We start with includes
-<details><summary>Click to show the header includes</summary>
+## 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>
+```
-#include <dumux/common/properties.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>
-#include <dumux/freeflow/navierstokes/model.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/fluidsystems/1pliquid.hh>
#include <dumux/material/components/constant.hh>
+```
+We include the problem header used for this simulation.
+
+```cpp
#include "problem.hh"
```
</details>
-Then we start setting the properties.
-1. For every test problem, a new `TypeTag` has to be created, which is done within the namespace `TTag` (subnamespace of `Properties`). It inherits from the Navier-Stokes flow model and the staggered-grid discretization scheme.
-2. The grid is chosen to be a two-dimensional Cartesian structured grid (`YaspGrid`).
-3. We set the `FluidSystem` to be a one-phase liquid with a single component. The class `Component::Constant` refers to a component with constant fluid properties (density, viscosity, ...) that can be set via the input file in the group `[0.Component]` where the number is the identifier given as template argument to the class template `Component::Constant`.
-4. The problem class `ChannelExampleProblem`, which is forward declared before we enter `namespace Dumux` and defined later in this file, is defined to be the problem used in this test problem (charaterized by the TypeTag `ChannelExample`). The fluid system, which contains information about the properties such as density, viscosity or diffusion coefficient of the fluid we're simulating, is set to a constant one phase liquid.
-5. We enable caching for the following classes (which stores values that were already calculated for later usage and thus results in higher memory usage but improved CPU speed): the grid volume variables, the grid flux variables, the finite volume grid geometry.
+### Type tag definition
+
+We define a type tag for our simulation with the name `ChannelExample`
+and inherit the properties specialized for the type tags `NavierStokes` and `StaggeredFreeFlowModel`.
+This way, most of the properties required for Navier-Stokes single-phase flow simulations
+using the staggered-grid scheme are conveniently specialized for our new type tag.
+However, some properties depend on user choices and no meaningful default value can be set.
+Those properties will be addressed later in this file.
+Please note that, in this example, we actually want to solve the Stokes instead of the
+Navier-Stokes equations. This can be achieved at runtime by setting the parameter
+`Problem.EnableInertiaTerms = false`. Have a look at the input file `params.input`
+to see how this is done in this example.
```cpp
+// We enter the namespace Dumux::Properties in order to import the entire Dumux namespace for general use:
namespace Dumux::Properties {
+// declaration of the `ChannelExample` type tag for the single-phase flow problem
namespace TTag {
struct ChannelExample { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
} // namespace TTag
+```
+
+### Property specializations
+
+In the following piece of code, mandatory properties for which no meaningful
+default can be set, are specialized for our type tag `ChannelExample`.
+```cpp
+// This sets the grid type used for the simulation. Here, we use a structured 2D grid.
template<class TypeTag>
struct Grid<TypeTag, TTag::ChannelExample> { 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::ChannelExample> { using type = Dumux::ChannelExampleProblem<TypeTag> ; };
+// 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::ChannelExample>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
};
+```
+
+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 EnableGridVolumeVariablesCache<TypeTag, TTag::ChannelExample> { static constexpr bool value = true; };
-
+//This enables grid wide caching for the flux variables.
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::ChannelExample> { static constexpr bool value = true; };
-
+// This enables grid-wide caching for the finite volume grid geometry
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::ChannelExample> { static constexpr bool value = true; };
} // end namespace Dumux::Properties
```
+</details>
+</details>
-## The file `problem.hh`
+## 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>
-### The problem class
-We enter the problem class where all necessary initial and boundary conditions are set for our simulation.
-As this is a Stokes problem, we inherit from the basic <code>NavierStokesProblem</code>.
-<details><summary>Toggle to show code:</summary>
+### Include files
+
+The only include we need here is the `NavierStokesProblem` class, the base
+class from which we will derive.
```cpp
-#include <dumux/common/properties.hh>
#include <dumux/freeflow/navierstokes/problem.hh>
+```
+### The problem class
+We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
+As we are solving a problem related to free flow, we inherit from the base class `NavierStokesProblem`.
+
+```cpp
namespace Dumux {
template <class TypeTag>
class ChannelExampleProblem : public NavierStokesProblem<TypeTag>
{
-```
-
-</details>
-
-We use convenient declarations that we derive from the property system.
-<details>
-<summary>Toggle to show code (convenient declarations)</summary>
-
-
-```cpp
+ // A few convenience aliases used throughout this class.
using ParentType = NavierStokesProblem<TypeTag>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
@@ -150,18 +238,9 @@ We use convenient declarations that we derive from the property system.
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
-```
-
-</details>
-
-There follows the constructor of our problem class:
-Within the constructor, we set the inlet velocity to a run-time specified value.
-As no run-time value is specified, we set the outlet pressure to 1.1e5 Pa.
-<details>
-<summary>Toggle to expand code (constructor)</summary>
-
-
-```cpp
+ // This is the constructor of our problem class:
+ // Within the constructor, we set the inlet velocity to a run-time specified value.
+ // If no run-time value is specified, we set the outlet pressure to 1.1e5 Pa.
ChannelExampleProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
@@ -170,39 +249,34 @@ As no run-time value is specified, we set the outlet pressure to 1.1e5 Pa.
}
```
-</details>
-
-Now, we define the type of initial and boundary conditions depending on location.
-Two types of boundary conditions can be specified: Dirichlet and Neumann. On a Dirichlet boundary,
-the values of the primary variables need to be fixed.
-On a Neumann boundary condition, values for derivatives need to be fixed.
-When Dirichlet conditions are set for the pressure, the derivative of the velocity
-vector with respect to the direction normal to the boundary is automatically set to
-zero. This boundary condition is called in-/outflow boundary condition in Dumux.
-In the following we specify Dirichlet boundaries for velocity on the left of our domain
-if isInlet_ is true, Dirichlet boundaries for pressure on the right of our domain
-if isOutlet_ is true and specify Dirichlet boundaries for velocity on the top and bottom
-of our domain else.
-<details>
-<summary>Toggle to expand code (<code>boundaryTypesAtPos</code>)</summary>
-
+#### 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 boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
BoundaryTypes values;
- if(isInlet_(globalPos))
+ if (isInlet_(globalPos))
{
+ // We specify Dirichlet boundary conditions for the velocity on the left of our domain
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
}
- else if(isOutlet_(globalPos))
+ else if (isOutlet_(globalPos))
{
+ // We fix the pressure on the right side of the domain
values.setDirichlet(Indices::pressureIdx);
}
else
{
+ // We specify Dirichlet boundary conditions for the velocity on the remaining boundaries (lower and upper wall)
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
}
@@ -211,42 +285,31 @@ of our domain else.
}
```
-</details>
-
-Second, we specify the values for the Dirichlet boundaries. We need to fix the values of our primary variables.
-To ensure a no-slip boundary condition at the top and bottom of the channel, the Dirichlet velocity
-in x-direction is set to zero if not at the inlet.
-<details>
-<summary>Toggle to expand code (<code>dirichletAtPos</code>)</summary>
-
+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 dirichletAtPos(const GlobalPosition& globalPos) const
{
+ // Use the initial values as default Dirichlet values
PrimaryVariables values = initialAtPos(globalPos);
- if(!isInlet_(globalPos))
- {
+ // Set a no-slip condition at the top and bottom wall of the channel
+ if (!isInlet_(globalPos))
values[Indices::velocityXIdx] = 0.0;
- }
return values;
}
```
-</details>
-
-We specify the values for the initial conditions.
-We assign constant values for pressure and velocity components.
-<details>
-<summary>Toggle to expand code (<code>initialAtPos</code>)</summary>
-
+The following function defines the initial conditions.
```cpp
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
PrimaryVariables values;
+ // Set the pressure and velocity values
values[Indices::pressureIdx] = outletPressure_;
values[Indices::velocityXIdx] = inletVelocity_;
values[Indices::velocityYIdx] = 0.0;
@@ -255,219 +318,201 @@ We assign constant values for pressure and velocity components.
}
```
-</details>
-
+#### Temperature distribution
We need to specify a constant temperature for our isothermal problem.
-We set it to 10°C.
-<details>
-<summary>Toggle to expand code (<code>temperature</code>)</summary>
-
+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; }
-private:
```
-</details>
-
-The inlet is at the left side of the physical domain.
-<details>
-<summary>Toggle to expand code (<code>isInlet_</code>)</summary>
-
+The inlet is on the left side of the physical domain.
```cpp
+private:
bool isInlet_(const GlobalPosition& globalPos) const
- {
- return globalPos[0] < eps_;
- }
+ { return globalPos[0] < eps_; }
```
-</details>
-
-The outlet is at the right side of the physical domain.
-<details>
-<summary>Toggle to expand code (<code>isOutlet_</code>)</summary>
-
+The outlet is on the right side of the physical domain.
```cpp
bool isOutlet_(const GlobalPosition& globalPos) const
- {
- return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_;
- }
+ { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
```
-</details>
-
Finally, private variables are declared:
-<details>
-<summary>Toggle to expand code (private variables)</summary>
-
```cpp
- static constexpr Scalar eps_=1e-6;
+ static constexpr Scalar eps_ = 1e-6;
Scalar inletVelocity_;
Scalar outletPressure_;
-};
-}
-```
-
-</details>
-
-
-
-
-## The file `main.cc`
+}; // end class definition of ChannelExampleProblem
+} // end namespace Dumux
+```
-We look now at the main file for the channel problem.
-### Includes
-Necessary files are included.
-<details>
-<summary>Toggle to expand details</summary>
-First, the Dune configuration file is include, the standard header file for C++ to get time and date information
-and another standard header for in- and output.
+</details>
-<details>
-<summary>Toggle to expand code (includes of problem file and of standard headers)</summary>
-```cpp
-#include <config.h>
+## The main file (`main.cc`)
-#include <ctime>
-#include <iostream>
-```
+<details open>
+<summary><b>Click to hide/show the file documentation</b> (or inspect the [source code](main.cc))</summary>
-</details>
-
-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.
-<details>
-<summary>Toggle to expand code (dune includes)</summary>
+### Included header files
+<details><summary> Click to show includes</summary>
+These are DUNE helper classes related to parallel computations 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>
-#include <dune/istl/io.hh>
```
-</details>
-
-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 file `parameters.hh` contains the parameter class, which manages the definition of input parameters by a default
-value, the inputfile or the command line.
-The file `dumuxmessage.hh` contains the class defining the start and end message of the simulation.
-
-The file `seqsolverbackend.hh` contains the class, which defines the sequential linear solver backends.
-The nonlinear Newton's method is included, as well as the assembler, which assembles the linear systems for staggered-grid finite volume schemes.
-The containing class in the file `diffmethod.hh` defines the different differentiation methods used to compute the derivatives of the residual.
-
-We need the class in `staggeredvtkoutputmodule.hh` to simplify the writing of dumux simulation data to VTK format.
-The gridmanager constructs a grid from the information in the input or grid file. There is a specification for the
-different supported grid managers.
-
-The following class contains functionality for additional flux output to the console.
-<details>
-<summary>Toggle to expand code (dumux includes)</summary>
-
+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>
-#include <dumux/common/dumuxmessage.hh>
-#include <dune/common/deprecated.hh>
+```
-#include <dumux/linear/seqsolverbackend.hh>
+The following files contain the non-linear Newton solver, the available linear solver backends and the assembler for the linear
+systems arising from the staggered-grid discretization.
+
+```cpp
#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/staggeredfvassembler.hh>
-#include <dumux/assembly/diffmethod.hh>
+#include <dumux/assembly/diffmethod.hh> // analytic or numeric differentiation
+```
+
+The following class provides a convenient way of writing of dumux simulation results to VTK format.
+```cpp
#include <dumux/io/staggeredvtkoutputmodule.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>
+```
+This class contains functionality for additional flux output.
+
+```cpp
#include <dumux/freeflow/navierstokes/staggered/fluxoversurface.hh>
+```
+In this header, a `TypeTag` is defined, which collects
+the properties that are required for the simulation.
+It also contains the actual problem with initial and boundary conditions.
+For detailed information, please have a look
+at the documentation provided therein.
+
+```cpp
#include "properties.hh"
```
-</details>
</details>
-### Beginning of the main function
-We begin the main function by making the type tag `ChannelExample`, that we defined in `problem.hh` for this test problem available here.
-Then we initializing the message passing interface (MPI), even if we do not plan to run the application in parallel. Finalizing of the MPI is done automatically on exit.
-We continue by printing the dumux start message and parsing the command line arguments and runtimeparameters from the input file in the init 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:
```cpp
int main(int argc, char** argv) try
{
using namespace Dumux;
- using TypeTag = Properties::TTag::ChannelExample;
-
+ // The Dune MPIHelper must be instantiated for each program using Dune
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
- if (mpiHelper.rank() == 0)
- DumuxMessage::print(/*firstCall=*/true);
-
+ // 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 this type tag
+using the property system, i.e. with `GetPropType`.
-### Set-up and solving of the problem
-
-A gridmanager tries to create the grid either from a grid file or the input file. You can learn more about grids in
-Dumux in the [grid exercise](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/-/tree/master/exercises/exercise-grids).
-Then, we compute the leaf grid view, which is the overlap of all grid levels in hierarchical grids.
-We create (`std::make_shared`) and initialize (`update`) the finite volume geometry to build up the subcontrolvolumes
-and subcontrolvolume faces for each element of the grid partition.
-In the problem, we define the boundary and initial conditions.
-We set a solution vector which has a part (indexed by `cellCenterIdx`) for degrees of freedom (`dofs`) living
-in grid cell centers - pressures - and a part (indexed by `faceIdx`) for degrees of freedom livin on grid cell faces.
-We initialize the solution vector by what was defined as the initial solution in `problem.hh` (`applyInitialSolution`)
-and then use the solution vector to intialize the `gridVariables`. Grid variables in general contain the s
-primary variables (velocities, pressures) as well as secondary variables (density, viscosity, ...).
-We then initialize the vtkoutput. Each model has a predefined model-specific output with relevant parameters
-for that model. Here, it is pressure, velocity, density and process rank (relevant in the case of parallelisation).
+```cpp
+ using TypeTag = Properties::TTag::ChannelExample;
+```
+#### 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: Setting up and solving the 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 set a solution vector which consist of two parts: one part (indexed by `cellCenterIdx`)
+is for the pressure degrees of freedom (`dofs`) living in grid cell centers. Another part
+(indexed by `faceIdx`) is for degrees of freedom defining the normal velocities on grid cell faces.
+We initialize the solution vector by what was defined as the initial solution of the the problem.
+```cpp
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x;
x[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
x[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
problem->applyInitialSolution(x);
+```
+
+The grid variables are used store variables (primary and secondary variables) on sub-control volumes and faces (volume and flux variables).
+```cpp
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
+```
+
+We then initialize the predefined model-specific output vtk output.
+```cpp
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
vtkWriter.write(0.0);
```
-
+<details><summary> Click to show calculation of surface fluxes</summary>
We set up two surfaces over which fluxes are calculated.
We determine the extensions [xMin,xMax]x[yMin,yMax] of the physical domain.
The first surface (added by the first call of addSurface) shall be placed at the middle of the channel.
@@ -477,7 +522,6 @@ In this case, we add half a cell-width to the x-position in order to make sure t
the cell faces lie on the surface. This assumes a regular cartesian grid.
The second surface (second call of addSurface) is placed at the outlet of the channel.
-
```cpp
FluxOverSurface<GridVariables,
SolutionVector,
@@ -514,47 +558,54 @@ The second surface (second call of addSurface) is placed at the outlet of the ch
flux.addSurface("outlet", p0outlet, p1outlet);
```
-
-The incompressible Stokes equation depends only linearly on the velocity, so we could use a linear solver to solve the problem.
-Here, we use the show the more general case which would also work for incompressible fluids or the
-Navier-Stokes equation. We use non-linear Newton solver for the solution.
-In each step of the Newton solver, we assemble the respective linear system, including the jacobian matrix and the residual by the
-`assembler`. The linear systems are here solved by the direct linear solver `UMFPack`.
-As a postprocessing, we calculate mass and volume fluxes over the planes specified above.
-<details>
-<summary>Toggle to expand code (assembly, solution process, postprocessing)</summary>
-
+</details>
+We create and initialize the assembler for the stationary problem.
+This is where the Jacobian matrix for the Newton solver is assembled.
```cpp
using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
+```
+
+We use UMFPack as direct linear solver within each Newton iteration.
+```cpp
using LinearSolver = Dumux::UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
+```
+
+This example considers a linear problem (incompressible Stokes flow), therefore
+the non-linear Newton solver is not really necessary.
+For sake of generality, we nevertheless use it here such that the example can be easily
+changed to a non-linear problem by switching on the inertia terms in the input file or by choosing a compressible fluid.
+In the following piece of code we instantiate the non-linear newton solver and let it solve
+the problem.
+```cpp
+ // alias for and instantiation of the newton solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
+ // Solve the (potentially non-linear) system.
nonLinearSolver.solve(x);
+```
+
+In the following we calculate mass and volume fluxes over the planes specified above
+(you have to click to unfold the code showing how to set up the surface fluxes).
+```cpp
flux.calculateMassOrMoleFluxes();
flux.calculateVolumeFluxes();
```
-</details>
-
-### Final Output
-We write the vtk output and print the mass/energy/volume fluxes over the planes.
-We conclude by printing the dumux end message. After the end of the main function,
-possibly catched error messages are printed.
-<details>
-<summary>Toggle to expand code (final output)</summary>
-
+#### Final Output
+We write the VTK output and print the mass/energy/volume fluxes over the planes.
+We conclude by printing the dumux end message.
```cpp
vtkWriter.write(1.0);
- if(GetPropType<TypeTag, Properties::ModelTraits>::enableEnergyBalance())
+ if (GetPropType<TypeTag, Properties::ModelTraits>::enableEnergyBalance())
{
std::cout << "mass / energy flux at middle is: " << flux.netFlux("middle") << std::endl;
std::cout << "mass / energy flux at outlet is: " << flux.netFlux("outlet") << std::endl;
@@ -569,18 +620,25 @@ possibly catched error messages are printed.
std::cout << "volume flux at outlet is: " << flux.netFlux("outlet")[0] << std::endl;
if (mpiHelper.rank() == 0)
- {
Parameters::print();
- DumuxMessage::print(/*firstCall=*/false);
- }
return 0;
} // end main
+```
+
+#### 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 (Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
+// errors related to the parsing of Dune grid files
catch (Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
@@ -590,11 +648,13 @@ catch (Dune::DGFException & e)
<< " ---> Abort!" << std::endl;
return 2;
}
+// generic error handling with Dune::Exception
catch (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;
@@ -604,4 +664,5 @@ catch (...)
</details>
+</details>
diff --git a/examples/freeflowchannel/doc/_intro.md b/examples/freeflowchannel/doc/_intro.md
new file mode 100644
index 0000000000..4007f34ec7
--- /dev/null
+++ b/examples/freeflowchannel/doc/_intro.md
@@ -0,0 +1,65 @@
+# Free flow through a channel
+
+__You learn how to__
+
+* solve a free-flow channel problem
+* set outflow boundary conditions in the free-flow context
+
+__Results__. In this example we will obtain the following stationary velocity profile:
+
+
+
+__Table of contents__. This description is structured as follows:
+
+[[_TOC_]]
+
+## Mathematical model
+In this example, the Stokes model for stationary and incompressible single phase flow is considered.
+Thus, the momentum balance equations
+
+```math
+- \nabla\cdot\left(\mu\left(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{\text{T}}\right)\right)+ \nabla p = 0
+```
+
+and the mass balance
+
+```math
+\nabla \cdot \left(\boldsymbol{u}\right) =0
+```
+
+are solved, where $`\varrho`$ and $`\mu`$ are the density and viscosity of the fluid,
+$`\boldsymbol{u}`$ is the fluid velocity and $`p`$ is the pressure. Here, we use constant fluid
+properties with $`\varrho = 1~\frac{\text{kg}}{\text{m}^3}`$ and $`\mu = 1~\text{Pa}\text{s}`$.
+Furthermore, isothermal conditions with a homogeneous temperature distribution of $`T=10^\circ C`$ are assumed.
+
+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,
+have a look at the Dumux [handbook](https://dumux.org/handbook).
+
+## Problem set-up
+This example considers stationary flow of a fluid between two parallel solid plates in two dimensions.
+Flow is enforced from left to right by prescribing an inflow velocity of $` v = 1~\frac{\text{m}}{\text{s}} `$
+on the left boundary, while a fixed pressure of $`p = 1.1 \text{bar}`$ and a zero velocity gradient
+in x-direction are prescribed on the right boundary. On the top and bottom boundaries, no-slip
+conditions are applied, which cause a parabolic velocity profile to develop along the channel.
+Take a look at Figure 1 for an illustration of the domain and the boundary conditions.
+
+<figure>
+ <center>
+ <img src="img/setup.png" alt="Free-flow setup" width="80%"/>
+ <figcaption> <b> Fig.1 </b> - Setup for the free flow problem.</figcaption>
+ </center>
+</figure>
+
+# Implementation
+
+## Folder layout and files
+
+```
+└── freeflowchannel/
+ ├── CMakeLists.txt -> build system file
+ ├── main.cc -> main program flow
+ ├── params.input -> runtime parameters
+ ├── properties.hh -> compile time configuration
+ └── problem.hh -> boundary & initial conditions
+```
diff --git a/examples/freeflowchannel/doc/intro.md b/examples/freeflowchannel/doc/intro.md
deleted file mode 100644
index baec2112c9..0000000000
--- a/examples/freeflowchannel/doc/intro.md
+++ /dev/null
@@ -1,45 +0,0 @@
-# Freeflow through a channel
-
-__You learn how to__
-
-* solve a free flow channel problem
-* set outflow boundary conditions in the free-flow context
-
-__Results__. In this example we will obtain the following stationary velocity profile:
-
-
-
-__Table of contents__. This description is structured as follows:
-
-[[_TOC_]]
-
-## Mathematical model
-The Stokes model without gravitation and without sources or sinks for a stationary, incompressible, laminar, single phase, one-component, isothermal ($`T=10^\circ C`$) flow is considered assuming a Newtonian fluid of constant density $` \varrho = 1~\frac{\text{kg}}{\text{m}^3} `$ and constant kinematic viscosity $` \nu = 1~\frac{\text{m}^2}{\text{s}} `$. The momentum balance
-```math
-- \nabla\cdot\left(\mu\left(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{\text{T}}\right)\right)+ \nabla p = 0
-```
-with density $`\varrho`$, velocity $`\boldsymbol{u}`$, dynamic viscosity $`\mu=\varrho\nu`$ and pressure $`p`$ and the mass balance
-```math
-\nabla \cdot \left(\varrho\boldsymbol{u}\right) =0
-```
-are discretized using a staggered-grid finite-volume scheme as spatial discretization with pressures and velocity components as primary variables. For details on the discretization scheme, have a look at the Dumux [handbook](https://dumux.org/handbook).
-
-## Problem set-up
-This example contains a stationary free flow of a fluid through two parallel solid plates in two dimensions from left to right. The figure below shows the simulation set-up. The fluid flows into the system at the left with a constant velocity of $` v = 1~\frac{\text{m}}{\text{s}} `$. The inflow velocity profile is a block profile. Due to the no-slip, no-flow boundary conditions at the top and bottom plate, the velocity profile gradually assumes a parabolic shape along the channel. At the outlet, the pressure is fixed and a zero velocity gradient in x-direction is assumed. The physical domain, which is modeled is the rectangular domain $`x\in[0,10],~y\in[0,1]`$.
-
-
-
-In the following, we take a close look at the files containing the set-up: At first, boundary conditions are set in `problem.hh` for the Navier-Stokes model. Afterwards, we show the different steps for solving the model in the source file `main.cc`.
-
-# Implementation
-
-## Folder layout and files
-
-```
-└── freeflowchannel/
- ├── CMakeLists.txt -> build system file
- ├── main.cc -> main program flow
- ├── params.input -> runtime parameters
- ├── properties.hh -> compile time configuration
- └── problem.hh -> boundary & initial conditions
-```
diff --git a/examples/freeflowchannel/main.cc b/examples/freeflowchannel/main.cc
index d42b694d71..ee74271e2f 100644
--- a/examples/freeflowchannel/main.cc
+++ b/examples/freeflowchannel/main.cc
@@ -17,138 +17,123 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
-// ## The file `main.cc`
-//
-//
-// We look now at the main file for the channel problem.
-// ### Includes
-// Necessary files are included.
-//<details>
-// <summary>Toggle to expand details</summary>
-//
-// First, the Dune configuration file is include, the standard header file for C++ to get time and date information
-// and another standard header for in- and output.
-//
-//<details>
-// <summary>Toggle to expand code (includes of problem file and of standard headers)</summary>
-//
+ // ## The main file (`main.cc`)
+ // [[content]]
+ //
+ // ### Included header files
+ // [[details]] includes
+ // [[exclude]]
+ // Some generic includes.
#include <config.h>
-
#include <ctime>
#include <iostream>
-// </details>
-//
-// 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.
-//<details>
-// <summary>Toggle to expand code (dune includes)</summary>
-//
+// [[/exclude]]
+
+// These are DUNE helper classes related to parallel computations 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>
-#include <dune/istl/io.hh>
-// </details>
-//
-// 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 file `parameters.hh` contains the parameter class, which manages the definition of input parameters by a default
-// value, the inputfile or the command line.
-// The file `dumuxmessage.hh` contains the class defining the start and end message of the simulation.
-//
-// The file `seqsolverbackend.hh` contains the class, which defines the sequential linear solver backends.
-// The nonlinear Newton's method is included, as well as the assembler, which assembles the linear systems for staggered-grid finite volume schemes.
-// The containing class in the file `diffmethod.hh` defines the different differentiation methods used to compute the derivatives of the residual.
-//
-// We need the class in `staggeredvtkoutputmodule.hh` to simplify the writing of dumux simulation data to VTK format.
-// The gridmanager constructs a grid from the information in the input or grid file. There is a specification for the
-// different supported grid managers.
-//
-// The following class contains functionality for additional flux output to the console.
-//<details>
-// <summary>Toggle to expand code (dumux includes)</summary>
-//
+
+// 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>
#include <dumux/common/parameters.hh>
-#include <dumux/common/dumuxmessage.hh>
-#include <dune/common/deprecated.hh>
-#include <dumux/linear/seqsolverbackend.hh>
+// The following files contain the non-linear Newton solver, the available linear solver backends and the assembler for the linear
+// systems arising from the staggered-grid discretization.
#include <dumux/nonlinear/newtonsolver.hh>
+#include <dumux/linear/seqsolverbackend.hh>
#include <dumux/assembly/staggeredfvassembler.hh>
-#include <dumux/assembly/diffmethod.hh>
+#include <dumux/assembly/diffmethod.hh> // analytic or numeric differentiation
+// The following class provides a convenient way of writing of dumux simulation results to VTK format.
#include <dumux/io/staggeredvtkoutputmodule.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>
+// This class contains functionality for additional flux output.
#include <dumux/freeflow/navierstokes/staggered/fluxoversurface.hh>
+
+// In this header, a `TypeTag` is defined, which collects
+// the properties that are required for the simulation.
+// It also contains the actual problem with initial and boundary conditions.
+// For detailed information, please have a look
+// at the documentation provided therein.
#include "properties.hh"
-// </details>
-// </details>
+// [[/details]]
//
-// ### Beginning of the main function
-// We begin the main function by making the type tag `ChannelExample`, that we defined in `problem.hh` for this test problem available here.
-// Then we initializing the message passing interface (MPI), even if we do not plan to run the application in parallel. Finalizing of the MPI is done automatically on exit.
-// We continue by printing the dumux start message and parsing the command line arguments and runtimeparameters from the input file in the init 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]]
int main(int argc, char** argv) try
{
using namespace Dumux;
- using TypeTag = Properties::TTag::ChannelExample;
-
+ // The Dune MPIHelper must be instantiated for each program using Dune
const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
- if (mpiHelper.rank() == 0)
- DumuxMessage::print(/*firstCall=*/true);
-
+ // parse command line arguments and input file
Parameters::init(argc, argv);
- //
- // ### Set-up and solving of the problem
- //
- // A gridmanager tries to create the grid either from a grid file or the input file. You can learn more about grids in
- // Dumux in the [grid exercise](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/-/tree/master/exercises/exercise-grids).
- // Then, we compute the leaf grid view, which is the overlap of all grid levels in hierarchical grids.
- // We create (`std::make_shared`) and initialize (`update`) the finite volume geometry to build up the subcontrolvolumes
- // and subcontrolvolume faces for each element of the grid partition.
- // In the problem, we define the boundary and initial conditions.
- // We set a solution vector which has a part (indexed by `cellCenterIdx`) for degrees of freedom (`dofs`) living
- // in grid cell centers - pressures - and a part (indexed by `faceIdx`) for degrees of freedom livin on grid cell faces.
- // We initialize the solution vector by what was defined as the initial solution in `problem.hh` (`applyInitialSolution`)
- // and then use the solution vector to intialize the `gridVariables`. Grid variables in general contain the s
- // primary variables (velocities, pressures) as well as secondary variables (density, viscosity, ...).
- // We then initialize the vtkoutput. Each model has a predefined model-specific output with relevant parameters
- // for that model. Here, it is pressure, velocity, density and process rank (relevant in the case of parallelisation).
- //
+ // [[/codeblock]]
+
+ // 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 this type tag
+ // using the property system, i.e. with `GetPropType`.
+ using TypeTag = Properties::TTag::ChannelExample;
+
+ // #### 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.
+ // [[codeblock]]
GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
gridManager.init();
+ // We compute on the leaf grid view.
const auto& leafGridView = gridManager.grid().leafGridView();
+ // [[/codeblock]]
+ // #### Step 2: Setting up and solving the 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.
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.
using Problem = GetPropType<TypeTag, Properties::Problem>;
auto problem = std::make_shared<Problem>(gridGeometry);
+ // We set a solution vector which consist of two parts: one part (indexed by `cellCenterIdx`)
+ // is for the pressure degrees of freedom (`dofs`) living in grid cell centers. Another part
+ // (indexed by `faceIdx`) is for degrees of freedom defining the normal velocities on grid cell faces.
+ // We initialize the solution vector by what was defined as the initial solution of the the problem.
using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
SolutionVector x;
x[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
x[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
problem->applyInitialSolution(x);
+ // The grid variables are used store variables (primary and secondary variables) on sub-control volumes and faces (volume and flux variables).
using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
gridVariables->init(x);
+ // We then initialize the predefined model-specific output vtk output.
using IOFields = GetPropType<TypeTag, Properties::IOFields>;
StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
IOFields::initOutputModule(vtkWriter); // Add model specific output fields
vtkWriter.write(0.0);
- //
+
+ // [[details]] calculation of surface fluxes
// We set up two surfaces over which fluxes are calculated.
// We determine the extensions [xMin,xMax]x[yMin,yMax] of the physical domain.
// The first surface (added by the first call of addSurface) shall be placed at the middle of the channel.
@@ -157,7 +142,6 @@ int main(int argc, char** argv) try
// In this case, we add half a cell-width to the x-position in order to make sure that
// the cell faces lie on the surface. This assumes a regular cartesian grid.
// The second surface (second call of addSurface) is placed at the outlet of the channel.
- //
FluxOverSurface<GridVariables,
SolutionVector,
GetPropType<TypeTag, Properties::ModelTraits>,
@@ -191,41 +175,43 @@ int main(int argc, char** argv) try
const auto p0outlet = GlobalPosition{xMax, yMin};
const auto p1outlet = GlobalPosition{xMax, yMax};
flux.addSurface("outlet", p0outlet, p1outlet);
- //
- // The incompressible Stokes equation depends only linearly on the velocity, so we could use a linear solver to solve the problem.
- // Here, we use the show the more general case which would also work for incompressible fluids or the
- // Navier-Stokes equation. We use non-linear Newton solver for the solution.
- // In each step of the Newton solver, we assemble the respective linear system, including the jacobian matrix and the residual by the
- // `assembler`. The linear systems are here solved by the direct linear solver `UMFPack`.
- // As a postprocessing, we calculate mass and volume fluxes over the planes specified above.
- //<details>
- // <summary>Toggle to expand code (assembly, solution process, postprocessing)</summary>
- //
+ // [[/details]]
+
+ // We create and initialize the assembler for the stationary problem.
+ // This is where the Jacobian matrix for the Newton solver is assembled.
using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
+ // We use UMFPack as direct linear solver within each Newton iteration.
using LinearSolver = Dumux::UMFPackBackend;
auto linearSolver = std::make_shared<LinearSolver>();
+ // This example considers a linear problem (incompressible Stokes flow), therefore
+ // the non-linear Newton solver is not really necessary.
+ // For sake of generality, we nevertheless use it here such that the example can be easily
+ // changed to a non-linear problem by switching on the inertia terms in the input file or by choosing a compressible fluid.
+ // In the following piece of code we instantiate the non-linear newton solver and let it solve
+ // the problem.
+ // [[codeblock]]
+ // alias for and instantiation of the newton solver
using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
NewtonSolver nonLinearSolver(assembler, linearSolver);
+ // Solve the (potentially non-linear) system.
nonLinearSolver.solve(x);
-
+ // [[/codeblock]]
+ // In the following we calculate mass and volume fluxes over the planes specified above
+ // (you have to click to unfold the code showing how to set up the surface fluxes).
flux.calculateMassOrMoleFluxes();
flux.calculateVolumeFluxes();
- // </details>
- //
- // ### Final Output
- // We write the vtk output and print the mass/energy/volume fluxes over the planes.
- // We conclude by printing the dumux end message. After the end of the main function,
- // possibly catched error messages are printed.
- //<details>
- // <summary>Toggle to expand code (final output)</summary>
- //
+
+ // #### Final Output
+ // We write the VTK output and print the mass/energy/volume fluxes over the planes.
+ // We conclude by printing the dumux end message.
+ // [[codeblock]]
vtkWriter.write(1.0);
- if(GetPropType<TypeTag, Properties::ModelTraits>::enableEnergyBalance())
+ if (GetPropType<TypeTag, Properties::ModelTraits>::enableEnergyBalance())
{
std::cout << "mass / energy flux at middle is: " << flux.netFlux("middle") << std::endl;
std::cout << "mass / energy flux at outlet is: " << flux.netFlux("outlet") << std::endl;
@@ -240,18 +226,23 @@ int main(int argc, char** argv) try
std::cout << "volume flux at outlet is: " << flux.netFlux("outlet")[0] << std::endl;
if (mpiHelper.rank() == 0)
- {
Parameters::print();
- DumuxMessage::print(/*firstCall=*/false);
- }
return 0;
} // end main
+// [[/codeblock]]
+// #### Exception handling
+// In this part of the main file we catch and print possible exceptions that could
+// occur during the simulation.
+// [[details]] exception handler
+// [[codeblock]]
+// errors related to run-time parameters
catch (Dumux::ParameterException &e)
{
std::cerr << std::endl << e << " ---> Abort!" << std::endl;
return 1;
}
+// errors related to the parsing of Dune grid files
catch (Dune::DGFException & e)
{
std::cerr << "DGF exception thrown (" << e <<
@@ -261,15 +252,18 @@ catch (Dune::DGFException & e)
<< " ---> Abort!" << std::endl;
return 2;
}
+// generic error handling with Dune::Exception
catch (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>
-//
+// [[/codeblock]]
+// [[/details]]
+// [[/content]]
diff --git a/examples/freeflowchannel/params.input b/examples/freeflowchannel/params.input
index 4f1d82a9bc..623d2f8e56 100644
--- a/examples/freeflowchannel/params.input
+++ b/examples/freeflowchannel/params.input
@@ -6,7 +6,7 @@ Cells = 100 50 # [m]
Name = example_ff_channel # name passed to the output routines
InletVelocity = 1 # [m/s]
EnableGravity = false
-EnableInertiaTerms = false # Stokes, not Navier Stokes, in this example
+EnableInertiaTerms = false # only consider the Stokes equations in this example
[Component]
LiquidDensity = 1.0 # [kg/m^3]
diff --git a/examples/freeflowchannel/problem.hh b/examples/freeflowchannel/problem.hh
index 927dcf023b..2707c73767 100644
--- a/examples/freeflowchannel/problem.hh
+++ b/examples/freeflowchannel/problem.hh
@@ -20,28 +20,29 @@
#ifndef DUMUX_EXAMPLES_FREEFLOW_CHANNEL_PROBLEM_HH
#define DUMUX_EXAMPLES_FREEFLOW_CHANNEL_PROBLEM_HH
-// ## The file `problem.hh`
+// ## 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.
//
-// ### The problem class
-// We enter the problem class where all necessary initial and boundary conditions are set for our simulation.
+// [[content]]
+//
+// ### Include files
//
-// As this is a Stokes problem, we inherit from the basic <code>NavierStokesProblem</code>.
-// <details><summary>Toggle to show code:</summary>
-#include <dumux/common/properties.hh>
+// The only include we need here is the `NavierStokesProblem` class, the base
+// class from which we will derive.
#include <dumux/freeflow/navierstokes/problem.hh>
+// ### The problem class
+// We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
+// As we are solving a problem related to free flow, we inherit from the base class `NavierStokesProblem`.
+// [[codeblock]]
namespace Dumux {
template <class TypeTag>
class ChannelExampleProblem : public NavierStokesProblem<TypeTag>
{
- // </details>
- //
- // We use convenient declarations that we derive from the property system.
- //<details>
- // <summary>Toggle to show code (convenient declarations)</summary>
- //
+ // A few convenience aliases used throughout this class.
using ParentType = NavierStokesProblem<TypeTag>;
using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
@@ -56,134 +57,109 @@ class ChannelExampleProblem : public NavierStokesProblem<TypeTag>
using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
public:
- // </details>
- //
- // There follows the constructor of our problem class:
+ // This is the constructor of our problem class:
// Within the constructor, we set the inlet velocity to a run-time specified value.
- // As no run-time value is specified, we set the outlet pressure to 1.1e5 Pa.
- //<details>
- // <summary>Toggle to expand code (constructor)</summary>
- //
+ // If no run-time value is specified, we set the outlet pressure to 1.1e5 Pa.
ChannelExampleProblem(std::shared_ptr<const GridGeometry> gridGeometry)
: ParentType(gridGeometry)
{
inletVelocity_ = getParam<Scalar>("Problem.InletVelocity");
outletPressure_ = getParam<Scalar>("Problem.OutletPressure", 1.1e5);
}
- // </details>
- //
- // Now, we define the type of initial and boundary conditions depending on location.
- // Two types of boundary conditions can be specified: Dirichlet and Neumann. On a Dirichlet boundary,
- // the values of the primary variables need to be fixed.
- // On a Neumann boundary condition, values for derivatives need to be fixed.
- // When Dirichlet conditions are set for the pressure, the derivative of the velocity
- // vector with respect to the direction normal to the boundary is automatically set to
- // zero. This boundary condition is called in-/outflow boundary condition in Dumux.
- // In the following we specify Dirichlet boundaries for velocity on the left of our domain
- // if isInlet_ is true, Dirichlet boundaries for pressure on the right of our domain
- // if isOutlet_ is true and specify Dirichlet boundaries for velocity on the top and bottom
- // of our domain else.
- //<details>
- // <summary>Toggle to expand code (<code>boundaryTypesAtPos</code>)</summary>
- //
- BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
+ // [[/codeblock]]
+
+ // #### 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.
+ // [[codeblock]]
+ BoundaryTypes boundaryTypesAtPos(const GlobalPosition& globalPos) const
{
BoundaryTypes values;
- if(isInlet_(globalPos))
+ if (isInlet_(globalPos))
{
+ // We specify Dirichlet boundary conditions for the velocity on the left of our domain
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
}
- else if(isOutlet_(globalPos))
+ else if (isOutlet_(globalPos))
{
+ // We fix the pressure on the right side of the domain
values.setDirichlet(Indices::pressureIdx);
}
else
{
+ // We specify Dirichlet boundary conditions for the velocity on the remaining boundaries (lower and upper wall)
values.setDirichlet(Indices::velocityXIdx);
values.setDirichlet(Indices::velocityYIdx);
}
return values;
}
- // </details>
- //
- // Second, we specify the values for the Dirichlet boundaries. We need to fix the values of our primary variables.
- // To ensure a no-slip boundary condition at the top and bottom of the channel, the Dirichlet velocity
- // in x-direction is set to zero if not at the inlet.
- //<details>
- // <summary>Toggle to expand code (<code>dirichletAtPos</code>)</summary>
- //
- PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const
+ // [[/codeblock]]
+
+ // The following function specifies the __values on Dirichlet boundaries__.
+ // We need to define values for the primary variables (velocity and pressure).
+ // [[codeblock]]
+ PrimaryVariables dirichletAtPos(const GlobalPosition& globalPos) const
{
+ // Use the initial values as default Dirichlet values
PrimaryVariables values = initialAtPos(globalPos);
- if(!isInlet_(globalPos))
- {
+ // Set a no-slip condition at the top and bottom wall of the channel
+ if (!isInlet_(globalPos))
values[Indices::velocityXIdx] = 0.0;
- }
return values;
}
- // </details>
- //
- // We specify the values for the initial conditions.
- // We assign constant values for pressure and velocity components.
- //<details>
- // <summary>Toggle to expand code (<code>initialAtPos</code>)</summary>
- //
- PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
+ // [[/codeblock]]
+
+ // The following function defines the initial conditions.
+ // [[codeblock]]
+ PrimaryVariables initialAtPos(const GlobalPosition& globalPos) const
{
PrimaryVariables values;
+ // Set the pressure and velocity values
values[Indices::pressureIdx] = outletPressure_;
values[Indices::velocityXIdx] = inletVelocity_;
values[Indices::velocityYIdx] = 0.0;
return values;
}
- // </details>
- //
+ // [[/codeblock]]
+
+ // #### Temperature distribution
// We need to specify a constant temperature for our isothermal problem.
- // We set it to 10°C.
- //<details>
- // <summary>Toggle to expand code (<code>temperature</code>)</summary>
- //
+ // Fluid properties that depend on temperature will be calculated with this value.
+ // This would be important if another fluidsystem was used.
Scalar temperature() const
{ return 273.15 + 10; }
+
+// The inlet is on the left side of the physical domain.
+// [[codeblock]]
private:
- // </details>
- //
- // The inlet is at the left side of the physical domain.
- //<details>
- // <summary>Toggle to expand code (<code>isInlet_</code>)</summary>
- //
bool isInlet_(const GlobalPosition& globalPos) const
- {
- return globalPos[0] < eps_;
- }
- // </details>
- //
- // The outlet is at the right side of the physical domain.
- //<details>
- // <summary>Toggle to expand code (<code>isOutlet_</code>)</summary>
- //
+ { return globalPos[0] < eps_; }
+ // [[/codeblock]]
+
+ // The outlet is on the right side of the physical domain.
bool isOutlet_(const GlobalPosition& globalPos) const
- {
- return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_;
- }
- // </details>
- //
+ { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; }
+
// Finally, private variables are declared:
- //<details>
- // <summary>Toggle to expand code (private variables)</summary>
- //
- static constexpr Scalar eps_=1e-6;
+ // [[codeblock]]
+ static constexpr Scalar eps_ = 1e-6;
Scalar inletVelocity_;
Scalar outletPressure_;
-};
-}
+
+}; // end class definition of ChannelExampleProblem
+} // end namespace Dumux
+// [[/codeblock]]
+// [[/content]]
#endif
-// </details>
-//
diff --git a/examples/freeflowchannel/properties.hh b/examples/freeflowchannel/properties.hh
index 06dc602b5c..56e717db39 100644
--- a/examples/freeflowchannel/properties.hh
+++ b/examples/freeflowchannel/properties.hh
@@ -20,57 +20,102 @@
#ifndef DUMUX_EXAMPLES_FREEFLOW_CHANNEL_PROPERTIES_HH
#define DUMUX_EXAMPLES_FREEFLOW_CHANNEL_PROPERTIES_HH
-// ## The file `properties.hh`
-// In the following, we set the properties for our simulation
-// (click [here](https://git.iws.uni-stuttgart.de/dumux-repositories/dumux-course/blob/master/slides/dumux-course-properties.pdf) for DuMux course slides on the property system).
-// We start with includes
-// <details><summary>Click to show the header includes</summary>
-#include <dune/grid/yaspgrid.hh>
+// ## 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.
+//
+// [[content]]
+//
+// ### Includes
+// [[details]] includes
+//
+// 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.
+#include <dumux/freeflow/navierstokes/model.hh>
-#include <dumux/common/properties.hh>
+// We want to use `YaspGrid`, an implementation of the dune grid interface for structured grids:
+#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:
#include <dumux/discretization/staggered/freeflow/properties.hh>
-#include <dumux/freeflow/navierstokes/model.hh>
+// The fluid properties are specified in the following headers (we use a liquid with constant properties as the fluid phase):
#include <dumux/material/fluidsystems/1pliquid.hh>
#include <dumux/material/components/constant.hh>
+// We include the problem header used for this simulation.
#include "problem.hh"
-// </details>
+// [[/details]]
//
-// Then we start setting the properties.
-// 1. For every test problem, a new `TypeTag` has to be created, which is done within the namespace `TTag` (subnamespace of `Properties`). It inherits from the Navier-Stokes flow model and the staggered-grid discretization scheme.
-// 2. The grid is chosen to be a two-dimensional Cartesian structured grid (`YaspGrid`).
-// 3. We set the `FluidSystem` to be a one-phase liquid with a single component. The class `Component::Constant` refers to a component with constant fluid properties (density, viscosity, ...) that can be set via the input file in the group `[0.Component]` where the number is the identifier given as template argument to the class template `Component::Constant`.
-// 4. The problem class `ChannelExampleProblem`, which is forward declared before we enter `namespace Dumux` and defined later in this file, is defined to be the problem used in this test problem (charaterized by the TypeTag `ChannelExample`). The fluid system, which contains information about the properties such as density, viscosity or diffusion coefficient of the fluid we're simulating, is set to a constant one phase liquid.
-// 5. We enable caching for the following classes (which stores values that were already calculated for later usage and thus results in higher memory usage but improved CPU speed): the grid volume variables, the grid flux variables, the finite volume grid geometry.
+// ### Type tag definition
+//
+// We define a type tag for our simulation with the name `ChannelExample`
+// and inherit the properties specialized for the type tags `NavierStokes` and `StaggeredFreeFlowModel`.
+// This way, most of the properties required for Navier-Stokes single-phase flow simulations
+// using the staggered-grid scheme are conveniently specialized for our new type tag.
+// However, some properties depend on user choices and no meaningful default value can be set.
+// Those properties will be addressed later in this file.
+// Please note that, in this example, we actually want to solve the Stokes instead of the
+// Navier-Stokes equations. This can be achieved at runtime by setting the parameter
+// `Problem.EnableInertiaTerms = false`. Have a look at the input file `params.input`
+// to see how this is done in this example.
// [[codeblock]]
+// We enter the namespace Dumux::Properties in order to import the entire Dumux namespace for general use:
namespace Dumux::Properties {
+// declaration of the `ChannelExample` type tag for the single-phase flow problem
namespace TTag {
struct ChannelExample { using InheritsFrom = std::tuple<NavierStokes, StaggeredFreeFlowModel>; };
} // namespace TTag
+// [[/codeblock]]
+// ### Property specializations
+//
+// In the following piece of code, mandatory properties for which no meaningful
+// default can be set, are specialized for our type tag `ChannelExample`.
+// [[codeblock]]
+// This sets the grid type used for the simulation. Here, we use a structured 2D grid.
template<class TypeTag>
struct Grid<TypeTag, TTag::ChannelExample> { 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::ChannelExample> { using type = Dumux::ChannelExampleProblem<TypeTag> ; };
+// 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::ChannelExample>
{
using Scalar = GetPropType<TypeTag, Properties::Scalar>;
using type = FluidSystems::OnePLiquid<Scalar, Components::Constant<1, Scalar> >;
};
+// [[/codeblock]]
+// We also set some properties related to memory management
+// throughout the simulation.
+// [[details]] caching properties
+//
+// 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.
+//
+// [[codeblock]]
+// This enables grid-wide caching of the volume variables.
template<class TypeTag>
struct EnableGridVolumeVariablesCache<TypeTag, TTag::ChannelExample> { static constexpr bool value = true; };
-
+//This enables grid wide caching for the flux variables.
template<class TypeTag>
struct EnableGridFluxVariablesCache<TypeTag, TTag::ChannelExample> { static constexpr bool value = true; };
-
+// This enables grid-wide caching for the finite volume grid geometry
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::ChannelExample> { static constexpr bool value = true; };
} // end namespace Dumux::Properties
// [[/codeblock]]
-
+// [[/details]]
+// [[/content]]
#endif
--
GitLab