README.md 30.8 KB
Newer Older
1
2
<!-- Important: This file has been automatically generated by generate_example_docs.py. Do not edit this file directly! -->

Martin Utz's avatar
Martin Utz committed
3
# Shallow water flow with bottom friction
Martin Utz's avatar
Martin Utz committed
4
5
6
7
This example shows how the shallow water flow model can be
applied to simulate steady subcritical flow including
bottom friction (bed shear stress).

8
__You will learn how to__
Martin Utz's avatar
Martin Utz committed
9

10
11
* solve a shallow water flow problem including bottom friction
* compute and output (VTK) an analytical reference solution
Martin Utz's avatar
Martin Utz committed
12

13
__Result__. The numerical and analytical solutions for the problem will look like this:
14
15
16

![Result Logo](img/result.png)

17
__Table of contents__. This description is structured as follows:
18

19
[[_TOC_]]
20
21
22
23
24

## Mathematical model
The 2D shallow water equations (SWEs) are given by

```math
Martin Utz's avatar
Martin Utz committed
25
26
27
\frac{\partial \mathbf{U}}{\partial t} +
\frac{\partial \mathbf{F}}{\partial x} +
\frac{\partial \mathbf{G}}{\partial y} - \mathbf{S_b} - \mathbf{S_f} = 0
28
```
Martin Utz's avatar
Martin Utz committed
29

30
where $`\mathbf{U}`$, $`\mathbf{F}`$ and $`\mathbf{G}`$ defined as
Martin Utz's avatar
Martin Utz committed
31

32
```math
Martin Utz's avatar
Martin Utz committed
33
34
35
\mathbf{U} = \begin{bmatrix} h \\ uh \\ vh \end{bmatrix},
\mathbf{F} = \begin{bmatrix} hu \\ hu^2  + \frac{1}{2} gh^2 \\ huv \end{bmatrix},
\mathbf{G} = \begin{bmatrix} hv \\ huv \\ hv^2  + \frac{1}{2} gh^2 \end{bmatrix}
36
```
Martin Utz's avatar
Martin Utz committed
37

38
39
$`Z`$ is the bedSurface, $`h`$ the water depth, $`u`$ the velocity in
x-direction and $`v`$ the velocity in y-direction, $`g`$ is the constant of gravity.
Martin Utz's avatar
Martin Utz committed
40

41
42
43
44
The source terms for the bed friction $`\mathbf{S_b}`$ and bed slope
$`\mathbf{S_f}`$ are given as

```math
Martin Utz's avatar
Martin Utz committed
45
46
47
\mathbf{S_b} = \begin{bmatrix} 0 \\ -gh \frac{\partial z}{\partial x}
               \\ -gh \frac{\partial z}{\partial y}\end{bmatrix},
\mathbf{S_f} = \begin{bmatrix} 0 \\ -ghS_{fx} \\ -ghS_{fy}\end{bmatrix}.
48
```
Martin Utz's avatar
Martin Utz committed
49

50
For this example, a cell-centered finite volume method (`cctpfa`) is applied to solve the SWEs
Martin Utz's avatar
Martin Utz committed
51
52
in combination with a fully-implicit time discretization. For cases where no sharp fronts or
traveling waves occur it is possible to apply time steps larger than CFL number = 1 to reduce
53
the computation time. Even if a steady state solution is considered, an implicit time stepping method
Martin Utz's avatar
Martin Utz committed
54
55
is applied.

Martin Utz's avatar
Martin Utz committed
56
## Problem set-up
57

Martin Utz's avatar
Martin Utz committed
58
The model domain is given by a rough channel with a slope of 0.001.
Martin Utz's avatar
Martin Utz committed
59
The domain is 500 meters long and 10 meters wide.
60
![Domain](img/domain.png).
Martin Utz's avatar
Martin Utz committed
61

Martin Utz's avatar
Martin Utz committed
62
63
64
65
66
Bottom friction is considered by applying
the friction law of Manning (Manning n = 0.025). At the lateral sides no friction is considered and  a
no-flow no slip boundary condition is applied. This is the default boundary condition for the shallow water model.

At the left border a discharge boundary condition
67
68
is applied as inflow boundary condition with q = -1.0 ($`m^2 s^{-1}`$). At the right border a water fixed depth boundary condition
is applied for the outflow. Normal flow is assumed, therefore the water depth at the right border is calculated after
Martin Utz's avatar
Martin Utz committed
69
70
the of Gaukler-Manning-Strickler equation:

71
72
73
```math
v_m = n^{-1} R_{hy}^{2/3} I_s^{1/2}
```
Martin Utz's avatar
Martin Utz committed
74

75
Where the mean velocity $`v_m`$ is given as $`v_m = \frac{q}{h}`$,
Martin Utz's avatar
Martin Utz committed
76
77
$`n`$ is the friction value after Manning. $`R_{hy}`$ the hydraulic radius, which is assumed to be equal to
the water depth. $`I_s`$ is the bed slope and $`q`$ the unity inflow discharge
Martin Utz's avatar
Martin Utz committed
78
79

The water depth h can be calculated as
80
81
82
```math
h = \left(\frac{n q}{\sqrt{I_s}} \right)^{3/5}
```
Martin Utz's avatar
Martin Utz committed
83
84

The formula of Gaukler Manning and Strickler is also used to calculate the analytic solution. All parameters
85
86
87
for the simulation are given in the file `params.input`.

# Implementation
Martin Utz's avatar
Martin Utz committed
88

89
90
91
92
93
94
95
96
97
98
99
100
## Folder layout and files

```
└── shallowwaterfriction/
    ├── CMakeLists.txt          -> build system file
    ├── main.cc                 -> main program flow
    ├── params.input            -> runtime parameters
    ├── properties.hh           -> compile time configuration
    ├── problem.hh              -> boundary & initial conditions
    └── spatialparams.hh        -> spatial parameter fields
```

Martin Utz's avatar
Martin Utz committed
101

102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
## The file `properties.hh`

The header includes will be mentioned in the text below.
<details><summary>Click to show the header includes</summary>

```cpp
#include <dune/grid/yaspgrid.hh>

#include <dumux/common/properties.hh>
#include <dumux/discretization/cctpfa.hh>
#include <dumux/freeflow/shallowwater/model.hh>

#include "spatialparams.hh"
#include "problem.hh"
```

</details>

Let's define the properties for our simulation

```cpp
namespace Dumux::Properties {
```

First, a so-called TypeTag is created. Properties are traits specialized for this TypeTag (a simple `struct`).
The properties of two other TypeTags are inherited by adding the alias `InheritsFrom`.
Here, properties from the shallow water model (`TTag::ShallowWater`) and the
cell-centered finite volume scheme with two-point-flux approximation (`TTag::CCTpfaModel`)
are inherited. These other TypeTag definitions can be found in the included
headers `dumux/freeflow/shallowwater/model.hh` and `dumux/discretization/cctpfa.hh`.

```cpp
namespace TTag {
struct RoughChannel { using InheritsFrom = std::tuple<ShallowWater, CCTpfaModel>; };
}
```

We use a structured Cartesian grid with tensor product structure.
`Dune::YaspGrid` (Yet Another Structure Parallel Grid) is defined in `dune/grid/yaspgrid.hh`
in the Dune module `dune-grid`.

```cpp
template<class TypeTag>
struct Grid<TypeTag, TTag::RoughChannel>
{ using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates<GetPropType<TypeTag, Properties::Scalar>, 2> >; };
```

Next, we specialize the properties `Problem` and `SpatialParams` for our new TypeTag and
set the type to our problem and spatial parameter classes implemented
in `problem.hh` and `spatialparams.hh`.

```cpp
template<class TypeTag>
struct Problem<TypeTag, TTag::RoughChannel>
{ using type = Dumux::RoughChannelProblem<TypeTag>; };

template<class TypeTag>
struct SpatialParams<TypeTag, TTag::RoughChannel>
{
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
    using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;

    using type = RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>;
};
```

Finally, we enable caching for the grid geometry. The cache
stores values that were already calculated for later usage.
This makes the simulation run faster but it uses more memory.

```cpp
template<class TypeTag>
struct EnableGridGeometryCache<TypeTag, TTag::RoughChannel>
{ static constexpr bool value = true; };

} // end namespace Dumux::Properties
```




Martin Utz's avatar
Martin Utz committed
185
186
187
## The file `spatialparams.hh`


Martin Utz's avatar
Martin Utz committed
188
We include the basic spatial parameters for finite volumes file from which we will inherit
189

Martin Utz's avatar
Martin Utz committed
190
191
192
```cpp
#include <dumux/material/spatialparams/fv.hh>
```
193

Martin Utz's avatar
Martin Utz committed
194
The parameters header is needed to retrieve run-time parameters.
195

Martin Utz's avatar
Martin Utz committed
196
197
198
```cpp
#include <dumux/common/parameters.hh>
```
199

Martin Utz's avatar
Martin Utz committed
200
We include all friction laws, between we can choose for the calculation of the bottom friction source.
201

Martin Utz's avatar
Martin Utz committed
202
203
204
205
206
```cpp
#include <dumux/material/fluidmatrixinteractions/frictionlaws/frictionlaw.hh>
#include <dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh>
#include <dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh>
```
207

Martin Utz's avatar
Martin Utz committed
208
We enter the namespace Dumux. All Dumux functions and classes are in a namespace Dumux, to make sure they don`t clash with symbols from other libraries you may want to use in conjunction with Dumux.
209

Martin Utz's avatar
Martin Utz committed
210
211
212
```cpp
namespace Dumux {
```
213

Martin Utz's avatar
Martin Utz committed
214
In the RoughChannelSpatialParams class we define all functions needed to describe the spatial distributed parameters.
215

Martin Utz's avatar
Martin Utz committed
216
```cpp
217
template<class GridGeometry, class Scalar, class VolumeVariables>
Martin Utz's avatar
Martin Utz committed
218
class RoughChannelSpatialParams
219
220
: public FVSpatialParams<GridGeometry, Scalar,
                         RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>>
Martin Utz's avatar
Martin Utz committed
221
222
{
```
223

Martin Utz's avatar
Martin Utz committed
224
We introduce using declarations that are derived from the property system which we need in this class
225

Martin Utz's avatar
Martin Utz committed
226
```cpp
227
228
229
230
    using ThisType = RoughChannelSpatialParams<GridGeometry, Scalar, VolumeVariables>;
    using ParentType = FVSpatialParams<GridGeometry, Scalar, ThisType>;
    using GridView = typename GridGeometry::GridView;
    using FVElementGeometry = typename GridGeometry::LocalView;
Martin Utz's avatar
Martin Utz committed
231
232
233
234
235
236
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;
    using Element = typename GridView::template Codim<0>::Entity;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;

public:
```
237

Martin Utz's avatar
Martin Utz committed
238
In the constructor be read some values from the `params.input` and initialize the friciton law.
239

Martin Utz's avatar
Martin Utz committed
240
```cpp
241
242
    RoughChannelSpatialParams(std::shared_ptr<const GridGeometry> gridGeometry)
    : ParentType(gridGeometry)
Martin Utz's avatar
Martin Utz committed
243
244
245
246
247
248
249
    {
        gravity_ = getParam<Scalar>("Problem.Gravity");
        bedSlope_ = getParam<Scalar>("Problem.BedSlope");
        frictionLawType_ = getParam<std::string>("Problem.FrictionLaw");
        initFrictionLaw();
    }
```
250

Martin Utz's avatar
Martin Utz committed
251
We initialize the friction law based on the law specified in `params.input`.
252

Martin Utz's avatar
Martin Utz committed
253
254
255
256
257
258
259
260
```cpp
    void initFrictionLaw()
    {
      if (frictionLawType_ == "Manning")
      {
          Scalar manningN = getParam<Scalar>("Problem.ManningN");
          frictionLaw_ = std::make_unique<FrictionLawManning<VolumeVariables>>(gravity_, manningN);
      }
261
      else if (frictionLawType_ == "Nikuradse")
Martin Utz's avatar
Martin Utz committed
262
263
264
265
266
267
      {
          Scalar ks = getParam<Scalar>("Problem.Ks");
          frictionLaw_ = std::make_unique<FrictionLawNikuradse<VolumeVariables>>(ks);
      }
      else
      {
Martin Utz's avatar
Martin Utz committed
268
          std::cout<<"The FrictionLaw in params.input is unknown. Valid entries are `Manning` and `Nikuradse`!"<<std::endl;
Martin Utz's avatar
Martin Utz committed
269
270
271
      }
    }
```
272

Martin Utz's avatar
Martin Utz committed
273
Use this function, if you want to vary the value for the gravity.
274

Martin Utz's avatar
Martin Utz committed
275
276
277
278
279
280
```cpp
    Scalar gravity(const GlobalPosition& globalPos) const
    {
        return gravity_;
    }
```
281

Martin Utz's avatar
Martin Utz committed
282
Use this function for a constant gravity.
283

Martin Utz's avatar
Martin Utz committed
284
285
286
287
288
289
```cpp
    Scalar gravity() const
    {
        return gravity_;
    }
```
290

Martin Utz's avatar
Martin Utz committed
291
This function returns an object of the friction law class, which is initialized with the appropriate friction values. If you want to use different friciton values or laws, you have to use a vector of unique_ptr for `frictionLaw_` and pick the right friction law instances via the `element` argument.
292

Martin Utz's avatar
Martin Utz committed
293
294
295
296
297
298
299
```cpp
    const FrictionLaw<VolumeVariables>& frictionLaw(const Element& element,
                                                    const SubControlVolume& scv) const
    {
        return *frictionLaw_;
    }
```
300

Martin Utz's avatar
Martin Utz committed
301
Define the bed surface based on the `bedSlope_`.
302

Martin Utz's avatar
Martin Utz committed
303
304
305
306
307
308
309
310
311
312
313
314
315
316
```cpp
    Scalar bedSurface(const Element& element,
                      const SubControlVolume& scv) const
    {
        return 10.0 - element.geometry().center()[0] * bedSlope_;
    }

private:
    Scalar gravity_;
    Scalar bedSlope_;
    std::string frictionLawType_;
    std::unique_ptr<FrictionLaw<VolumeVariables>> frictionLaw_;
};
```
317

Martin Utz's avatar
Martin Utz committed
318
end of namespace Dumux.
319

Martin Utz's avatar
Martin Utz committed
320
321
322
323
324
325
```cpp
}
```



326

Martin Utz's avatar
Martin Utz committed
327
## The file `problem.hh`
328
We start with includes
329

Martin Utz's avatar
Martin Utz committed
330
331
```cpp
#include <dumux/common/parameters.hh>
332
#include <dumux/common/properties.hh>
Martin Utz's avatar
Martin Utz committed
333
334
335
#include <dumux/freeflow/shallowwater/problem.hh>
#include <dumux/freeflow/shallowwater/boundaryfluxes.hh>
```
336

Martin Utz's avatar
Martin Utz committed
337
338
We enter the problem class where all necessary boundary conditions and initial conditions are set for our simulation.
As this is a shallow water problem, we inherit from the basic ShallowWaterProblem.
339

Martin Utz's avatar
Martin Utz committed
340
```cpp
341
342
namespace Dumux {

Martin Utz's avatar
Martin Utz committed
343
344
345
346
template <class TypeTag>
class RoughChannelProblem : public ShallowWaterProblem<TypeTag>
{
```
347

Martin Utz's avatar
Martin Utz committed
348
We use convenient declarations that we derive from the property system.
349

Martin Utz's avatar
Martin Utz committed
350
351
352
353
354
355
```cpp
    using ParentType = ShallowWaterProblem<TypeTag>;
    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
    using BoundaryTypes = GetPropType<TypeTag, Properties::BoundaryTypes>;
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    using Indices = typename GetPropType<TypeTag, Properties::ModelTraits>::Indices;
356
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
Martin Utz's avatar
Martin Utz committed
357
358
    using NeumannFluxes = GetPropType<TypeTag, Properties::NumEqVector>;
    using ElementVolumeVariables = typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView;
359
360
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
    using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;
Martin Utz's avatar
Martin Utz committed
361
    using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;
362
    using FVElementGeometry = typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView;
Martin Utz's avatar
Martin Utz committed
363
    using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;
364
    using GridView = typename GetPropType<TypeTag, Properties::GridGeometry>::GridView;
Martin Utz's avatar
Martin Utz committed
365
366
367
368
369
370
371
    using Element = typename GridView::template Codim<0>::Entity;
    using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
    using NumEqVector = GetPropType<TypeTag, Properties::NumEqVector>;
    using SubControlVolume = typename FVElementGeometry::SubControlVolume;

public:
```
372

Martin Utz's avatar
Martin Utz committed
373
This is the constructor of our problem class.
374

Martin Utz's avatar
Martin Utz committed
375
```cpp
376
377
    RoughChannelProblem(std::shared_ptr<const GridGeometry> gridGeometry)
    : ParentType(gridGeometry)
Martin Utz's avatar
Martin Utz committed
378
379
    {
```
380

Martin Utz's avatar
Martin Utz committed
381
We read the parameters from the params.input file.
382

Martin Utz's avatar
Martin Utz committed
383
384
385
386
387
388
```cpp
        name_ = getParam<std::string>("Problem.Name");
        constManningN_ = getParam<Scalar>("Problem.ManningN");
        bedSlope_ = getParam<Scalar>("Problem.BedSlope");
        discharge_ = getParam<Scalar>("Problem.Discharge");
```
389

390
We calculate the outflow boundary condition using the Gauckler-Manning-Strickler formula.
391

Martin Utz's avatar
Martin Utz committed
392
393
394
```cpp
        hBoundary_ = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);
```
395

Martin Utz's avatar
Martin Utz committed
396
We initialize the analytic solution to a verctor of the appropriate size filled with zeros.
397

Martin Utz's avatar
Martin Utz committed
398
```cpp
399
400
        exactWaterDepth_.resize(gridGeometry->numDofs(), 0.0);
        exactVelocityX_.resize(gridGeometry->numDofs(), 0.0);
Martin Utz's avatar
Martin Utz committed
401
402
    }
```
403

Martin Utz's avatar
Martin Utz committed
404
Get the analytical water depth
405

Martin Utz's avatar
Martin Utz committed
406
407
408
409
410
411
```cpp
    const std::vector<Scalar>& getExactWaterDepth()
    {
        return exactWaterDepth_;
    }
```
412

Martin Utz's avatar
Martin Utz committed
413
Get the analytical velocity
414

Martin Utz's avatar
Martin Utz committed
415
416
417
418
419
420
```cpp
    const std::vector<Scalar>& getExactVelocityX()
    {
        return exactVelocityX_;
    }
```
421

422
Get the water depth with Gauckler-Manning-Strickler
423

Martin Utz's avatar
Martin Utz committed
424
425
426
427
428
429
430
431
432
433
```cpp
    Scalar gauklerManningStrickler(Scalar discharge, Scalar manningN, Scalar bedSlope)
    {
        using std::pow;
        using std::abs;
        using std::sqrt;

        return pow(abs(discharge)*manningN/sqrt(bedSlope), 0.6);
    }
```
434

Martin Utz's avatar
Martin Utz committed
435
Get the analytical solution
436

Martin Utz's avatar
Martin Utz committed
437
438
439
440
441
```cpp
    void analyticalSolution()
    {
        using std::abs;

442
        for (const auto& element : elements(this->gridGeometry().gridView()))
Martin Utz's avatar
Martin Utz committed
443
444
445
446
        {
            const Scalar h = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);
            const Scalar u = abs(discharge_)/h;

447
            const auto eIdx = this->gridGeometry().elementMapper().index(element);
Martin Utz's avatar
Martin Utz committed
448
449
450
451
452
            exactWaterDepth_[eIdx] = h;
            exactVelocityX_[eIdx] = u;
        }
    }
```
453

Martin Utz's avatar
Martin Utz committed
454
Get the problem name. It is used as a prefix for files generated by the simulation.
455

Martin Utz's avatar
Martin Utz committed
456
457
458
459
460
461
```cpp
    const std::string& name() const
    {
        return name_;
    }
```
462

Martin Utz's avatar
Martin Utz committed
463
Get the source term.
464

Martin Utz's avatar
Martin Utz committed
465
466
467
468
469
470
471
472
473
```cpp
     NumEqVector source(const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
                        const SubControlVolume &scv) const
    {

        NumEqVector source (0.0);
```
474

Martin Utz's avatar
Martin Utz committed
475
In this model the bottom friction is the only source.
476

Martin Utz's avatar
Martin Utz committed
477
478
479
480
481
482
```cpp
        source += bottomFrictionSource(element, fvGeometry, elemVolVars, scv);

        return source;
    }
```
483

Martin Utz's avatar
Martin Utz committed
484
Get the source term due to bottom friction.
485

Martin Utz's avatar
Martin Utz committed
486
487
488
489
490
491
492
493
494
```cpp
     NumEqVector bottomFrictionSource(const Element& element,
                                      const FVElementGeometry& fvGeometry,
                                      const ElementVolumeVariables& elemVolVars,
                                      const SubControlVolume &scv) const
     {
        NumEqVector bottomFrictionSource(0.0);
        const auto& volVars = elemVolVars[scv];
```
495

Martin Utz's avatar
Martin Utz committed
496
For the calculation of the source term due to bottom friction the two-dimensional bottom shear stess vector is needed. This is the force per area, which works between the flow and the bed. It is calculated within the `FrictionLaw`, which is a spatialParameter. In this model the `FrictionLawManning` is used (see `params.input`).
497

Martin Utz's avatar
Martin Utz committed
498
499
500
```cpp
        Dune::FieldVector<Scalar, 2> bottomShearStress = this->spatialParams().frictionLaw(element, scv).shearStress(volVars);
```
501

Martin Utz's avatar
Martin Utz committed
502
The bottom shear stress causes a pure loss of momentum. Thus the first entry of the `bottomFrictionSource`, which is related to the mass balance equation is zero. The second entry of the `bottomFricitonSource` corresponds to the momentum equation in x-direction and is therefore equal to the first, the x-component, of the `bottomShearStress`. Accordingly the third entry of the `bottomFrictionSource` is equal to the second component of the `bottomShearStress`.
503

Martin Utz's avatar
Martin Utz committed
504
505
506
507
508
509
510
511
```cpp
        bottomFrictionSource[0] = 0.0;
        bottomFrictionSource[1] = bottomShearStress[0];
        bottomFrictionSource[2] = bottomShearStress[1];

        return bottomFrictionSource;
     }
```
512

Martin Utz's avatar
Martin Utz committed
513
We specify the boundary condition type.
514

Martin Utz's avatar
Martin Utz committed
515
516
517
518
519
```cpp
    BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const
    {
        BoundaryTypes bcTypes;
```
520

Martin Utz's avatar
Martin Utz committed
521
Since we use a weak imposition all boundary conditions are of Neumann type.
522

Martin Utz's avatar
Martin Utz committed
523
524
525
526
527
```cpp
        bcTypes.setAllNeumann();
        return bcTypes;
    }
```
528
529
530
531
532
533
534

We specify the neumann boundary. Due to the weak imposition we calculate the flux at the
boundary, with a  Rieman solver. For this the state of a virtual cell outside of the boundary
is needed (`boundaryStateVariables`), wich is calculated with the Riemann invariants
(see Yoon and Kang, Finite Volume Model for Two-Dimensional Shallow Water Flows on Unstructured Grids).
The calculation of the Riemann invariants differ depending on the type of the boundary (h, q or no-flow boundary).

Martin Utz's avatar
Martin Utz committed
535
536
537
538
```cpp
    NeumannFluxes neumann(const Element& element,
                          const FVElementGeometry& fvGeometry,
                          const ElementVolumeVariables& elemVolVars,
539
                          const ElementFluxVariablesCache& elemFluxVarsCache,
Martin Utz's avatar
Martin Utz committed
540
541
542
543
544
545
546
547
548
549
                          const SubControlVolumeFace& scvf) const
    {
        NeumannFluxes values(0.0);

        const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
        const auto& insideVolVars = elemVolVars[insideScv];
        const auto& nxy = scvf.unitOuterNormal();
        const auto gravity = this->spatialParams().gravity(scvf.center());
        std::array<Scalar, 3> boundaryStateVariables;
```
550

Martin Utz's avatar
Martin Utz committed
551
Calculate the rieman invariants for imposed discharge at the left side.
552

Martin Utz's avatar
Martin Utz committed
553
554
555
556
557
558
559
560
561
562
563
```cpp
        if (scvf.center()[0] < 0.0 + eps_)
        {
            boundaryStateVariables = ShallowWater::fixedDischargeBoundary(discharge_,
                                                                          insideVolVars.waterDepth(),
                                                                          insideVolVars.velocity(0),
                                                                          insideVolVars.velocity(1),
                                                                          gravity,
                                                                          nxy);
        }
```
564

Martin Utz's avatar
Martin Utz committed
565
Calculate the rieman invariants for impose water depth at the right side.
566

Martin Utz's avatar
Martin Utz committed
567
568
569
570
571
572
573
574
575
576
577
```cpp
        else if (scvf.center()[0] > 100.0 - eps_)
        {
            boundaryStateVariables =  ShallowWater::fixedWaterDepthBoundary(hBoundary_,
                                                                            insideVolVars.waterDepth(),
                                                                            insideVolVars.velocity(0),
                                                                            insideVolVars.velocity(1),
                                                                            gravity,
                                                                            nxy);
        }
```
578

Martin Utz's avatar
Martin Utz committed
579
Calculate the rieman invarianty for the no-flow boundary.
580

Martin Utz's avatar
Martin Utz committed
581
582
583
584
585
586
587
588
```cpp
        else
        {
            boundaryStateVariables[0] = insideVolVars.waterDepth();
            boundaryStateVariables[1] = -insideVolVars.velocity(0);
            boundaryStateVariables[2] = -insideVolVars.velocity(1);
        }
```
589

Martin Utz's avatar
Martin Utz committed
590
We calculate the boundary fluxes based on a Riemann problem.
591

Martin Utz's avatar
Martin Utz committed
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
```cpp
        auto riemannFlux = ShallowWater::riemannProblem(insideVolVars.waterDepth(),
                                                        boundaryStateVariables[0],
                                                        insideVolVars.velocity(0),
                                                        boundaryStateVariables[1],
                                                        insideVolVars.velocity(1),
                                                        boundaryStateVariables[2],
                                                        insideVolVars.bedSurface(),
                                                        insideVolVars.bedSurface(),
                                                        gravity,
                                                        nxy);

        values[Indices::massBalanceIdx] = riemannFlux[0];
        values[Indices::velocityXIdx]   = riemannFlux[1];
        values[Indices::velocityYIdx]   = riemannFlux[2];

        return values;
    }
```
611

Martin Utz's avatar
Martin Utz committed
612
We set the initial conditions. In this example constant initial conditions are used. Therefore the argument `globalPos` is not needed. If you want to impose spatial variable initial conditions, you have to use the `globalPos`.
613

Martin Utz's avatar
Martin Utz committed
614
615
616
617
618
```cpp
    PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const
    {
        PrimaryVariables values(0.0);
```
619

Martin Utz's avatar
Martin Utz committed
620
We set the initial water depth to one meter.
621

Martin Utz's avatar
Martin Utz committed
622
623
624
```cpp
        values[0] = 1.0;
```
625

Martin Utz's avatar
Martin Utz committed
626
We set the x-component of the initial velocity to zero.
627

Martin Utz's avatar
Martin Utz committed
628
629
630
```cpp
        values[1] = 0.0;
```
631

Martin Utz's avatar
Martin Utz committed
632
We set the y-component of the initial velocity to zero.
633

Martin Utz's avatar
Martin Utz committed
634
635
636
637
638
639
```cpp
        values[2] = 0.0;

        return values;
    };
```
640

Martin Utz's avatar
Martin Utz committed
641
\}
642

Martin Utz's avatar
Martin Utz committed
643
644
645
646
```cpp

private:
```
647

648
We declare the private variables of the problem. They are initialized in the problems constructor.
Martin Utz's avatar
Martin Utz committed
649
We declare the variable for the analytic solution.
650

Martin Utz's avatar
Martin Utz committed
651
652
653
654
```cpp
    std::vector<Scalar> exactWaterDepth_;
    std::vector<Scalar> exactVelocityX_;
```
655

Martin Utz's avatar
Martin Utz committed
656
constant friction value. An analytic solution is only available for const friction. If you want to run the simulation with a non constant friciton value (specified in the spatialParams) you have to remove the analytic solution.
657

Martin Utz's avatar
Martin Utz committed
658
659
660
```cpp
    Scalar constManningN_;
```
661

Martin Utz's avatar
Martin Utz committed
662
The constant bed slope.
663

Martin Utz's avatar
Martin Utz committed
664
665
666
```cpp
    Scalar bedSlope_;
```
667

Martin Utz's avatar
Martin Utz committed
668
The water depth at the outflow boundary.
669

Martin Utz's avatar
Martin Utz committed
670
671
672
```cpp
    Scalar hBoundary_;
```
673

Martin Utz's avatar
Martin Utz committed
674
The discharge at the inflow boundary.
675

Martin Utz's avatar
Martin Utz committed
676
677
678
```cpp
    Scalar discharge_;
```
679

Martin Utz's avatar
Martin Utz committed
680
eps is used as a small value for the definition of the boundry conditions
681

Martin Utz's avatar
Martin Utz committed
682
683
684
685
```cpp
    static constexpr Scalar eps_ = 1.0e-6;
    std::string name_;
};
686

687
} // end namespace Dumux
Martin Utz's avatar
Martin Utz committed
688
689
690
691
```



692

Martin Utz's avatar
Martin Utz committed
693
694
## The file `main.cc`

695

Martin Utz's avatar
Martin Utz committed
696
697
This is the main file for the shallow water example. Here we can see the programme sequence and how the system is solved using newton's method.
### Includes
698

Martin Utz's avatar
Martin Utz committed
699
700
701
```cpp
#include <config.h>
```
702

Martin Utz's avatar
Martin Utz committed
703
Standard header file for C++, to get time and date information.
704

Martin Utz's avatar
Martin Utz committed
705
706
707
```cpp
#include <ctime>
```
708

Martin Utz's avatar
Martin Utz committed
709
Standard header file for C++, for in- and output.
710

Martin Utz's avatar
Martin Utz committed
711
712
713
```cpp
#include <iostream>
```
714

Martin Utz's avatar
Martin Utz committed
715
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.
716

Martin Utz's avatar
Martin Utz committed
717
718
719
720
721
722
```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>
```
723

Martin Utz's avatar
Martin Utz committed
724
We need the following class to simplify the writing of dumux simulation data to VTK format.
725

Martin Utz's avatar
Martin Utz committed
726
727
728
```cpp
#include <dumux/io/vtkoutputmodule.hh>
```
729

Martin Utz's avatar
Martin Utz committed
730
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.
731

Martin Utz's avatar
Martin Utz committed
732
733
734
```cpp
#include <dumux/common/properties.hh>
```
735

Martin Utz's avatar
Martin Utz committed
736
The following file contains the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.
737

Martin Utz's avatar
Martin Utz committed
738
739
740
```cpp
#include <dumux/common/parameters.hh>
```
741

Martin Utz's avatar
Martin Utz committed
742
The file dumuxmessage.hh contains the class defining the start and end message of the simulation.
743

Martin Utz's avatar
Martin Utz committed
744
745
746
747
```cpp
#include <dumux/common/dumuxmessage.hh>
#include <dumux/common/defaultusagemessage.hh>
```
748

Martin Utz's avatar
Martin Utz committed
749
The gridmanager constructs a grid from the information in the input or grid file. There is a specification for the different supported grid managers.
750

Martin Utz's avatar
Martin Utz committed
751
752
753
```cpp
#include <dumux/io/grid/gridmanager.hh>
```
754

Martin Utz's avatar
Martin Utz committed
755
We include the linear solver to be used to solve the linear system
756

Martin Utz's avatar
Martin Utz committed
757
758
```cpp
#include <dumux/linear/amgbackend.hh>
Timo Koch's avatar
Timo Koch committed
759
#include <dumux/linear/linearsolvertraits.hh>
Martin Utz's avatar
Martin Utz committed
760
```
761

Martin Utz's avatar
Martin Utz committed
762
We include the nonlinear newtons method
763

Martin Utz's avatar
Martin Utz committed
764
765
766
```cpp
#include <dumux/nonlinear/newtonsolver.hh>
```
767

Martin Utz's avatar
Martin Utz committed
768
Further we include assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation)
769

Martin Utz's avatar
Martin Utz committed
770
771
772
```cpp
#include <dumux/assembly/fvassembler.hh>
```
773

774
We include the properties
775

Martin Utz's avatar
Martin Utz committed
776
```cpp
777
#include "properties.hh"
Martin Utz's avatar
Martin Utz committed
778
```
779

Martin Utz's avatar
Martin Utz committed
780
### Beginning of the main function
781

Martin Utz's avatar
Martin Utz committed
782
783
784
785
786
```cpp
int main(int argc, char** argv) try
{
    using namespace Dumux;
```
787

Martin Utz's avatar
Martin Utz committed
788
We define the type tag for this problem
789

Martin Utz's avatar
Martin Utz committed
790
791
792
```cpp
    using TypeTag = Properties::TTag::RoughChannel;
```
793

Martin Utz's avatar
Martin Utz committed
794
We initialize MPI, finalize is done automatically on exit
795

Martin Utz's avatar
Martin Utz committed
796
797
798
```cpp
    const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);
```
799

Martin Utz's avatar
Martin Utz committed
800
We print dumux start message
801

Martin Utz's avatar
Martin Utz committed
802
803
804
805
```cpp
    if (mpiHelper.rank() == 0)
        DumuxMessage::print(/*firstCall=*/true);
```
806

Martin Utz's avatar
Martin Utz committed
807
We parse command line arguments and input file
808

Martin Utz's avatar
Martin Utz committed
809
810
811
```cpp
    Parameters::init(argc, argv);
```
812

Martin Utz's avatar
Martin Utz committed
813
814
### Create the grid
A gridmanager tries to create the grid either from a grid file or the input file.
815

Martin Utz's avatar
Martin Utz committed
816
817
818
819
```cpp
    GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
    gridManager.init();
```
820

Martin Utz's avatar
Martin Utz committed
821
We compute on the leaf grid view
822

Martin Utz's avatar
Martin Utz committed
823
824
825
```cpp
    const auto& leafGridView = gridManager.grid().leafGridView();
```
826

Martin Utz's avatar
Martin Utz committed
827
828
829
830
### Setup and solving of the problem
#### Setup
We create and initialize the finite volume grid geometry, the problem, the linear system, including the jacobian matrix, the residual and the solution vector and the gridvariables.
We need the finite volume geometry to build up the subcontrolvolumes (scv) and subcontrolvolume faces (scvf) for each element of the grid partition.
831

Martin Utz's avatar
Martin Utz committed
832
```cpp
833
834
835
    using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
    auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
    gridGeometry->update();
Martin Utz's avatar
Martin Utz committed
836
```
837

Martin Utz's avatar
Martin Utz committed
838
In the problem, we define the boundary and initial conditions.
839

Martin Utz's avatar
Martin Utz committed
840
841
```cpp
    using Problem = GetPropType<TypeTag, Properties::Problem>;
842
    auto problem = std::make_shared<Problem>(gridGeometry);
Martin Utz's avatar
Martin Utz committed
843
```
844

Martin Utz's avatar
Martin Utz committed
845
We initialize the solution vector
846

Martin Utz's avatar
Martin Utz committed
847
848
```cpp
    using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
849
    SolutionVector x(gridGeometry->numDofs());
Martin Utz's avatar
Martin Utz committed
850
851
852
    problem->applyInitialSolution(x);
    auto xOld = x;
```
853

Martin Utz's avatar
Martin Utz committed
854
And then use the solutionvector to intialize the gridVariables.
855

Martin Utz's avatar
Martin Utz committed
856
857
```cpp
    using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
858
    auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
Martin Utz's avatar
Martin Utz committed
859
860
    gridVariables->init(x);
```
861

Martin Utz's avatar
Martin Utz committed
862
We get some time loop parameters from the input file.
863

Martin Utz's avatar
Martin Utz committed
864
865
866
867
868
869
```cpp
    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
    const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
    const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
    auto dt = getParam<Scalar>("TimeLoop.DtInitial");
```
870

Martin Utz's avatar
Martin Utz committed
871
We intialize the vtk output module. Each model has a predefined model specific output with relevant parameters for that model.
872

Martin Utz's avatar
Martin Utz committed
873
874
875
876
```cpp
    using IOFields = GetPropType<TypeTag, Properties::IOFields>;
    VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables,x, problem->name());
```
877

Martin Utz's avatar
Martin Utz committed
878
We add the analytical solution ("exactWaterDepth" and "exactVelocityX") to the predefined specific output.
879

Martin Utz's avatar
Martin Utz committed
880
881
882
883
```cpp
    vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth");
    vtkWriter.addField(problem->getExactVelocityX(), "exactVelocityX");
```
884

Martin Utz's avatar
Martin Utz committed
885
We calculate the analytic solution.
886

Martin Utz's avatar
Martin Utz committed
887
888
889
890
891
```cpp
    problem->analyticalSolution();
    IOFields::initOutputModule(vtkWriter);
    vtkWriter.write(0.0);
```
892

Martin Utz's avatar
Martin Utz committed
893
We instantiate time loop.
894

Martin Utz's avatar
Martin Utz committed
895
896
897
898
```cpp
    auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0, dt, tEnd);
    timeLoop->setMaxTimeStepSize(maxDt);
```
899

Martin Utz's avatar
Martin Utz committed
900
we set the assembler with the time loop because we have an instationary problem.
901

Martin Utz's avatar
Martin Utz committed
902
903
```cpp
    using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
904
    auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
Martin Utz's avatar
Martin Utz committed
905
```
906

Martin Utz's avatar
Martin Utz committed
907
We set the linear solver.
908

Martin Utz's avatar
Martin Utz committed
909
```cpp
Timo Koch's avatar
Timo Koch committed
910
    using LinearSolver = AMGBiCGSTABBackend<LinearSolverTraits<GridGeometry>>;
911
    auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
Martin Utz's avatar
Martin Utz committed
912
```
913

Martin Utz's avatar
Martin Utz committed
914
Additionaly, we set the non-linear solver.
915

Martin Utz's avatar
Martin Utz committed
916
917
918
919
```cpp
    using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
    NewtonSolver nonLinearSolver(assembler, linearSolver);
```
920

Martin Utz's avatar
Martin Utz committed
921
We set some check point at the end of the time loop. The check point is used to trigger the vtk output.
922

Martin Utz's avatar
Martin Utz committed
923
924
925
```cpp
    timeLoop->setCheckPoint(tEnd);
```
926

Martin Utz's avatar
Martin Utz committed
927
We start the time loop.
928

Martin Utz's avatar
Martin Utz committed
929
930
931
932
```cpp
    timeLoop->start(); do
    {
```
933

Martin Utz's avatar
Martin Utz committed
934
We start to calculate the new solution of that time step. First we define the old solution as the solution of the previous time step for storage evaluations.
935

Martin Utz's avatar
Martin Utz committed
936
937
938
```cpp
        assembler->setPreviousSolution(xOld);
```
939

Martin Utz's avatar
Martin Utz committed
940
We solve the non-linear system with time step control.
941

Martin Utz's avatar
Martin Utz committed
942
943
944
```cpp
        nonLinearSolver.solve(x,*timeLoop);
```
945

Martin Utz's avatar
Martin Utz committed
946
We make the new solution the old solution.
947

Martin Utz's avatar
Martin Utz committed
948
949
950
951
```cpp
        xOld = x;
        gridVariables->advanceTimeStep();
```
952

Martin Utz's avatar
Martin Utz committed
953
We advance to the time loop to the next step.
954

Martin Utz's avatar
Martin Utz committed
955
956
957
```cpp
        timeLoop->advanceTimeStep();
```
958

Martin Utz's avatar
Martin Utz committed
959
We write vtk output, if we reached the check point (end of time loop)
960

Martin Utz's avatar
Martin Utz committed
961
962
963
964
```cpp
        if (timeLoop->isCheckPoint())
            vtkWriter.write(timeLoop->time());
```
965

Martin Utz's avatar
Martin Utz committed
966
We report statistics of this time step.
967

Martin Utz's avatar
Martin Utz committed
968
969
970
```cpp
        timeLoop->reportTimeStep();
```
971

Martin Utz's avatar
Martin Utz committed
972
We set new dt as suggested by newton controller for the next time step.
973

Martin Utz's avatar
Martin Utz committed
974
975
976
977
978
979
980
981
```cpp
        timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));


    } while (!timeLoop->finished());

    timeLoop->finalize(leafGridView.comm());
```
982

Martin Utz's avatar
Martin Utz committed
983
984
### Final Output
We print dumux end message.
985

Martin Utz's avatar
Martin Utz committed
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
```cpp
    if (mpiHelper.rank() == 0)
    {
        Parameters::print();
        DumuxMessage::print(/*firstCall=*/false);
    }

    return 0;
} // end main

catch (const Dumux::ParameterException &e)
{
    std::cerr << std::endl << e << " ---> Abort!" << std::endl;
    return 1;
}
For faster browsing, not all history is shown. View entire blame