| + +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` + + | +
+ |
+
+
+ 
+ 
+
+
+
+ 
+ 
+
+
+
+
+```
+
+We want to use `YaspGrid`, an implementation of the dune grid interface for structured grids:
+
+```cpp
+#include
+```
+
+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
+```
+
+The fluid properties are specified in the following headers (we use a liquid with constant properties as the fluid phase):
+
+```cpp
+#include
+#include
+```
+
+We include the problem header used for this simulation.
+
+```cpp
+#include "problem.hh"
+```
+
+
+
+### 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; };
+} // 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
+struct FluidSystem
+{
+ using Scalar = GetPropType;
+ using type = FluidSystems::OnePLiquid >;
+};
+
+// This sets the grid type used for the simulation. Here, we use a structured 2D grid.
+template
+struct Grid { using type = Dune::YaspGrid<2>; };
+
+// This sets our problem class (see problem.hh) containing initial and boundary conditions.
+template
+struct Problem { using type = Dumux::LidDrivenCavityExampleProblem ; };
+```
+
+We also set some properties related to memory management
+throughout the simulation.
+
+struct EnableGridGeometryCache { static constexpr bool value = true; };
+//This enables grid wide caching for the flux variables.
+template
+struct EnableGridFluxVariablesCache { static constexpr bool value = true; };
+// This enables grid-wide caching for the finite volume grid geometry
+template
+struct EnableGridVolumeVariablesCache { static constexpr bool value = true; };
+} // end namespace Dumux::Properties
+```
+
+
+
+
+
+
+
+## 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.
+
+
+Click to hide/show the file documentation (or inspect the [source code](../properties.hh))
+ + +### Includes +Click to show 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. + +```cpp +#includeClick to show 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. + + +```cpp +// This enables grid-wide caching of the volume variables. +template
+
+#include
+```
+
+Include the `NavierStokesProblem` class, the base
+class from which we will derive.
+
+```cpp
+#include
+```
+
+Include the `NavierStokesBoundaryTypes` class which specifies the boundary types set in this problem.
+
+```cpp
+#include
+```
+
+### 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 LidDrivenCavityExampleProblem : public NavierStokesProblem
+{
+ using ParentType = NavierStokesProblem;
+
+ using BoundaryTypes = Dumux::NavierStokesBoundaryTypes::numEq()>;
+ using GridGeometry = GetPropType;
+ using FVElementGeometry = typename GridGeometry::LocalView;
+ using SubControlVolume = typename GridGeometry::SubControlVolume;
+ using Indices = typename GetPropType::Indices;
+ using NumEqVector = GetPropType;
+ using PrimaryVariables = GetPropType;
+ using Scalar = GetPropType;
+
+ 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 gridGeometry)
+ : ParentType(gridGeometry)
+ {
+ lidVelocity_ = getParam("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
+```
+
+
+
+
+
+
+## The main file (`main.cc`)
+
+Click to hide/show the file documentation (or inspect the [source code](../problem.hh))
+ + +### Include files + + +```cpp +#include
+
+```
+
+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
+#include
+```
+
+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
+#include
+#include
+```
+
+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
+```
+
+This class contains functionality for VTK output for models using the staggered finite volume scheme.
+
+```cpp
+#include
+```
+
+We include the problem header used for this simulation.
+
+```cpp
+#include "properties.hh"
+```
+
+
+
+The following function writes the velocities and coordinates at x = 0.5 and y = 0.5 into a log file.
+
+```cpp
+template
+void writeSteadyVelocityAndCoordinates(const Problem& problem, const SolutionVector &sol, const GridGeometry gridGeometry)
+{
+ std::ofstream logFilevx(problem->name() + "_vx.log"), logFilevy(problem->name() + "_vy.log");
+ logFilevx << "y vx\n";
+ logFilevy << "x vy\n";
+
+ static constexpr double eps_ = 1.0e-7;
+ for (const auto& element : elements(gridGeometry->gridView()))
+ {
+ auto fvGeometry = localView(*gridGeometry);
+ fvGeometry.bind(element);
+ for (const auto& scvf : scvfs(fvGeometry))
+ {
+ if (!scvf.boundary() && scvf.insideScvIdx() > scvf.outsideScvIdx())
+ {
+ const auto& globalPos = scvf.ipGlobal();
+ const auto velocity = sol[gridGeometry->faceIdx()][scvf.dofIndex()][0];
+
+ if (std::abs(globalPos[0]-0.5) < eps_)
+ logFilevx << globalPos[1] << " " << velocity << "\n";
+ else if (std::abs(globalPos[1]-0.5) < eps_)
+ logFilevy << globalPos[0] << " " << velocity << "\n";
+ }
+ }
+ }
+}
+```
+
+### 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)
+{
+ using namespace Dumux;
+
+ // The Dune MPIHelper must be instantiated for each program using Dune, it is finalized automatically on exit
+ const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
+
+ // 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`.
+
+```cpp
+ using TypeTag = Properties::TTag::LidDrivenCavityExample;
+```
+
+#### 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> 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;
+ auto gridGeometry = std::make_shared(leafGridView);
+ gridGeometry->update();
+```
+
+We now instantiate the problem, in which we define the boundary and initial conditions.
+
+```cpp
+ using Problem = GetPropType;
+ auto problem = std::make_shared(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;
+ SolutionVector x;
+ x[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
+ x[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
+ problem->applyInitialSolution(x);
+ auto xOld = x;
+```
+
+We use the initial solution vector to intialize the `gridVariables`.
+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;
+ auto gridVariables = std::make_shared(problem, gridGeometry);
+ gridVariables->init(x);
+```
+
+We get some time loop parameters from the input file
+and instantiate the time loop
+
+```cpp
+ using Scalar = GetPropType;
+ const auto tEnd = getParam("TimeLoop.TEnd");
+ const auto maxDt = getParam("TimeLoop.MaxTimeStepSize");
+ const auto dt = getParam("TimeLoop.DtInitial");
+
+ auto timeLoop = std::make_shared>(0, dt, tEnd);
+ timeLoop->setMaxTimeStepSize(maxDt);
+```
+
+We then initialize the predefined model-specific output vtk output.
+
+```cpp
+ using IOFields = GetPropType;
+ StaggeredVtkOutputModule vtkWriter(*gridVariables, x, problem->name());
+ IOFields::initOutputModule(vtkWriter); // Add model specific output fields
+ vtkWriter.write(0.0);
+```
+
+To solve the non-linear problem at hand, we use the `NewtonSolver`,
+which we have to tell how to assemble and solve the system in each
+iteration. Here, we use the direct linear solver UMFPack.
+
+```cpp
+ using Assembler = StaggeredFVAssembler;
+ auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld);
+
+ using LinearSolver = Dumux::UMFPackBackend;
+ auto linearSolver = std::make_shared();
+
+ using NewtonSolver = Dumux::NewtonSolver;
+ NewtonSolver nonLinearSolver(assembler, linearSolver);
+```
+
+##### The time loop
+In each time step, we solve the non-linear system of equations, write
+the current solution into .vtk files and prepare for the next time step.
+
+```cpp
+ timeLoop->start(); do
+ {
+ // We solve the non-linear system with time step control.
+ nonLinearSolver.solve(x, *timeLoop);
+
+ // We make the new solution the old solution.
+ xOld = x;
+ gridVariables->advanceTimeStep();
+
+ // We advance to the time loop to the next step.
+ timeLoop->advanceTimeStep();
+
+ // We write vtk output for each time step.
+ vtkWriter.write(timeLoop->time());
+
+ // We report statistics of this time step.
+ timeLoop->reportTimeStep();
+
+ // We set a new dt as suggested by the newton solver for the next time step.
+ timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
+
+ } while (!timeLoop->finished());
+```
+
+We write the velocities and coordinates at x = 0.5 and y = 0.5 into a file
+
+```cpp
+ writeSteadyVelocityAndCoordinates(problem, x, gridGeometry);
+```
+
+The following piece of code prints a final status report of the time loop
+before the program is terminated.
+
+```cpp
+ timeLoop->finalize(leafGridView.comm());
+
+ // print used and unused parameters
+ if (mpiHelper.rank() == 0)
+ Parameters::print();
+
+ return 0;
+} // end main
+```
+
+
+
+
+
+| [:arrow_left: Back to the main documentation](../README.md) | [:arrow_right: Continue with part 2](postprocessing.md) |
+|---|---:|
+
diff --git a/examples/liddrivencavity/doc/problem_intro.md b/examples/liddrivencavity/doc/problem_intro.md
new file mode 100644
index 0000000000000000000000000000000000000000..6ff6d9b269baf68ba6cd3aa4f6e5909f5d5d2ae4
--- /dev/null
+++ b/examples/liddrivencavity/doc/problem_intro.md
@@ -0,0 +1,5 @@
+# Part 1: Implementation
+
+The implementation of simulation setup and main flow is structured as follows:
+
+[[_TOC_]]
diff --git a/examples/liddrivencavity/img/lidverification.png b/examples/liddrivencavity/img/lidverification.png
new file mode 100644
index 0000000000000000000000000000000000000000..1293cf14dc10fcaaee46a172b3c986f03498dfbe
Binary files /dev/null and b/examples/liddrivencavity/img/lidverification.png differ
diff --git a/examples/liddrivencavity/img/result.svg b/examples/liddrivencavity/img/result.svg
new file mode 100644
index 0000000000000000000000000000000000000000..3640f745a4e57539a42c93343506f282e615851b
--- /dev/null
+++ b/examples/liddrivencavity/img/result.svg
@@ -0,0 +1,367 @@
+
+
diff --git a/examples/liddrivencavity/img/setup.png b/examples/liddrivencavity/img/setup.png
new file mode 100644
index 0000000000000000000000000000000000000000..220b9ee0877bea3ca48f0f4c77e93c36fb0feddb
Binary files /dev/null and b/examples/liddrivencavity/img/setup.png differ
diff --git a/examples/liddrivencavity/main.cc b/examples/liddrivencavity/main.cc
new file mode 100644
index 0000000000000000000000000000000000000000..d8c88134d8293980af9fce9ff33c91e1a519e09b
--- /dev/null
+++ b/examples/liddrivencavity/main.cc
@@ -0,0 +1,223 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see