Timo Koch committed Mar 31, 2020 1 2   Martin Utz committed Oct 08, 2019 3 # Shallow water flow with bottom friction  Martin Utz committed Oct 08, 2019 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).  Timo Koch committed Apr 02, 2020 8 __You will learn how to__  Martin Utz committed Oct 08, 2019 9   Timo Koch committed Apr 02, 2020 10 11 * solve a shallow water flow problem including bottom friction * compute and output (VTK) an analytical reference solution  Martin Utz committed Oct 08, 2019 12   Timo Koch committed Apr 02, 2020 13 __Result__. The numerical and analytical solutions for the problem will look like this:  Timo Koch committed Apr 02, 2020 14 15 16  ![Result Logo](img/result.png)  Timo Koch committed Apr 02, 2020 17 __Table of contents__. This description is structured as follows:  Timo Koch committed Apr 02, 2020 18   Timo Koch committed Apr 02, 2020 19 [[_TOC_]]  Timo Koch committed Apr 02, 2020 20 21 22 23 24  ## Mathematical model The 2D shallow water equations (SWEs) are given by math  Martin Utz committed Oct 08, 2019 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  Timo Koch committed Apr 02, 2020 28   Martin Utz committed Oct 08, 2019 29   Timo Koch committed Apr 02, 2020 30 where $\mathbf{U}$, $\mathbf{F}$ and $\mathbf{G}$ defined as  Martin Utz committed Oct 08, 2019 31   Timo Koch committed Apr 02, 2020 32 math  Martin Utz committed Oct 08, 2019 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}  Timo Koch committed Apr 02, 2020 36   Martin Utz committed Oct 08, 2019 37   Timo Koch committed Apr 02, 2020 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 committed Oct 08, 2019 40   Timo Koch committed Apr 02, 2020 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 committed Oct 08, 2019 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}.  Timo Koch committed Apr 02, 2020 48   Martin Utz committed Oct 08, 2019 49   Timo Koch committed Apr 02, 2020 50 For this example, a cell-centered finite volume method (cctpfa) is applied to solve the SWEs  Martin Utz committed Oct 08, 2019 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 committed Oct 08, 2019 54 55 is applied.  Martin Utz committed Oct 08, 2019 56 ## Problem set-up  Timo Koch committed Apr 02, 2020 57   Martin Utz committed Oct 08, 2019 58 The model domain is given by a rough channel with a slope of 0.001.  Martin Utz committed Oct 08, 2019 59 The domain is 500 meters long and 10 meters wide.  60 ![Domain](img/domain.png).  Martin Utz committed Oct 08, 2019 61   Martin Utz committed Oct 08, 2019 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 committed Oct 08, 2019 69 70 the of Gaukler-Manning-Strickler equation:  Timo Koch committed Apr 02, 2020 71 72 73 math v_m = n^{-1} R_{hy}^{2/3} I_s^{1/2}   Martin Utz committed Oct 08, 2019 74   Timo Koch committed Apr 02, 2020 75 Where the mean velocity $v_m$ is given as $v_m = \frac{q}{h}$,  Martin Utz committed Oct 08, 2019 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 committed Oct 08, 2019 78 79  The water depth h can be calculated as  Timo Koch committed Apr 02, 2020 80 81 82 math h = \left(\frac{n q}{\sqrt{I_s}} \right)^{3/5}   Martin Utz committed Oct 08, 2019 83 84  The formula of Gaukler Manning and Strickler is also used to calculate the analytic solution. All parameters  Timo Koch committed Apr 02, 2020 85 86 87 for the simulation are given in the file params.input. # Implementation  Martin Utz committed Oct 08, 2019 88   Timo Koch committed Apr 02, 2020 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 committed Oct 08, 2019 101   Timo Koch committed Apr 02, 2020 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.
Click to show the header includes cpp #include #include #include #include #include "spatialparams.hh" #include "problem.hh" 
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; }; }  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 struct Grid { using type = Dune::YaspGrid<2, Dune::TensorProductCoordinates, 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 struct Problem { using type = Dumux::RoughChannelProblem; }; template struct SpatialParams { using GridGeometry = GetPropType; using Scalar = GetPropType; using ElementVolumeVariables = typename GetPropType::LocalView; using VolumeVariables = typename ElementVolumeVariables::VolumeVariables; using type = RoughChannelSpatialParams; };  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 struct EnableGridGeometryCache { static constexpr bool value = true; }; } // end namespace Dumux::Properties   Martin Utz committed Oct 08, 2019 185 186 187 ## The file spatialparams.hh  Martin Utz committed Oct 08, 2019 188 We include the basic spatial parameters for finite volumes file from which we will inherit  Timo Koch committed Mar 30, 2020 189   Martin Utz committed Oct 08, 2019 190 191 192 cpp #include   Timo Koch committed Mar 30, 2020 193   Martin Utz committed Oct 08, 2019 194 The parameters header is needed to retrieve run-time parameters.  Timo Koch committed Mar 30, 2020 195   Martin Utz committed Oct 08, 2019 196 197 198 cpp #include   Timo Koch committed Mar 30, 2020 199   Martin Utz committed Oct 08, 2019 200 We include all friction laws, between we can choose for the calculation of the bottom friction source.  Timo Koch committed Mar 30, 2020 201   Martin Utz committed Oct 08, 2019 202 203 204 205 206 cpp #include #include #include   Timo Koch committed Mar 30, 2020 207   Martin Utz committed Oct 08, 2019 208 We enter the namespace Dumux. All Dumux functions and classes are in a namespace Dumux, to make sure they dont clash with symbols from other libraries you may want to use in conjunction with Dumux.  Timo Koch committed Mar 30, 2020 209   Martin Utz committed Oct 08, 2019 210 211 212 cpp namespace Dumux {   Timo Koch committed Mar 30, 2020 213   Martin Utz committed Oct 08, 2019 214 In the RoughChannelSpatialParams class we define all functions needed to describe the spatial distributed parameters.  Timo Koch committed Mar 30, 2020 215   Martin Utz committed Oct 08, 2019 216 cpp  Simon Scholz committed Nov 22, 2019 217 template  Martin Utz committed Oct 08, 2019 218 class RoughChannelSpatialParams  Simon Scholz committed Nov 22, 2019 219 220 : public FVSpatialParams>  Martin Utz committed Oct 08, 2019 221 222 {   Timo Koch committed Mar 30, 2020 223   Martin Utz committed Oct 08, 2019 224 We introduce using declarations that are derived from the property system which we need in this class  Timo Koch committed Mar 30, 2020 225   Martin Utz committed Oct 08, 2019 226 cpp  Simon Scholz committed Nov 22, 2019 227 228 229 230  using ThisType = RoughChannelSpatialParams; using ParentType = FVSpatialParams; using GridView = typename GridGeometry::GridView; using FVElementGeometry = typename GridGeometry::LocalView;  Martin Utz committed Oct 08, 2019 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:   Timo Koch committed Mar 30, 2020 237   Martin Utz committed Oct 08, 2019 238 In the constructor be read some values from the params.input and initialize the friciton law.  Timo Koch committed Mar 30, 2020 239   Martin Utz committed Oct 08, 2019 240 cpp  Simon Scholz committed Nov 22, 2019 241 242  RoughChannelSpatialParams(std::shared_ptr gridGeometry) : ParentType(gridGeometry)  Martin Utz committed Oct 08, 2019 243 244 245 246 247 248 249  { gravity_ = getParam("Problem.Gravity"); bedSlope_ = getParam("Problem.BedSlope"); frictionLawType_ = getParam("Problem.FrictionLaw"); initFrictionLaw(); }   Timo Koch committed Mar 30, 2020 250   Martin Utz committed Oct 08, 2019 251 We initialize the friction law based on the law specified in params.input.  Timo Koch committed Mar 30, 2020 252   Martin Utz committed Oct 08, 2019 253 254 255 256 257 258 259 260 cpp void initFrictionLaw() { if (frictionLawType_ == "Manning") { Scalar manningN = getParam("Problem.ManningN"); frictionLaw_ = std::make_unique>(gravity_, manningN); }  Simon Scholz committed Nov 22, 2019 261  else if (frictionLawType_ == "Nikuradse")  Martin Utz committed Oct 08, 2019 262 263 264 265 266 267  { Scalar ks = getParam("Problem.Ks"); frictionLaw_ = std::make_unique>(ks); } else {  Martin Utz committed Oct 08, 2019 268  std::cout<<"The FrictionLaw in params.input is unknown. Valid entries are Manning and Nikuradse!"<& frictionLaw(const Element& element, const SubControlVolume& scv) const { return *frictionLaw_; }   Timo Koch committed Mar 30, 2020 300   Martin Utz committed Oct 08, 2019 301 Define the bed surface based on the bedSlope_.  Timo Koch committed Mar 30, 2020 302   Martin Utz committed Oct 08, 2019 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_; };   Timo Koch committed Mar 30, 2020 317   Martin Utz committed Oct 08, 2019 318 end of namespace Dumux.  Timo Koch committed Mar 30, 2020 319   Martin Utz committed Oct 08, 2019 320 321 322 323 324 325 cpp }   Timo Koch committed Mar 30, 2020 326   Martin Utz committed Oct 08, 2019 327 ## The file problem.hh  Timo Koch committed Apr 02, 2020 328 We start with includes  Timo Koch committed Mar 30, 2020 329   Martin Utz committed Oct 08, 2019 330 331 cpp #include  Timo Koch committed Apr 02, 2020 332 #include  Martin Utz committed Oct 08, 2019 333 334 335 #include #include   Timo Koch committed Mar 30, 2020 336   Martin Utz committed Oct 08, 2019 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.  Timo Koch committed Mar 30, 2020 339   Martin Utz committed Oct 08, 2019 340 cpp  Timo Koch committed Apr 02, 2020 341 342 namespace Dumux {  Martin Utz committed Oct 08, 2019 343 344 345 346 template class RoughChannelProblem : public ShallowWaterProblem {   Timo Koch committed Mar 30, 2020 347   Martin Utz committed Oct 08, 2019 348 We use convenient declarations that we derive from the property system.  Timo Koch committed Mar 30, 2020 349   Martin Utz committed Oct 08, 2019 350 351 352 353 354 355 cpp using ParentType = ShallowWaterProblem; using PrimaryVariables = GetPropType; using BoundaryTypes = GetPropType; using Scalar = GetPropType; using Indices = typename GetPropType::Indices;  Simon Scholz committed Nov 22, 2019 356  using GridGeometry = GetPropType;  Martin Utz committed Oct 08, 2019 357 358  using NeumannFluxes = GetPropType; using ElementVolumeVariables = typename GetPropType::LocalView;  Katharina Heck committed Oct 08, 2019 359 360  using GridVariables = GetPropType; using ElementFluxVariablesCache = typename GridVariables::GridFluxVariablesCache::LocalView;  Martin Utz committed Oct 08, 2019 361  using VolumeVariables = typename ElementVolumeVariables::VolumeVariables;  Simon Scholz committed Nov 22, 2019 362  using FVElementGeometry = typename GetPropType::LocalView;  Martin Utz committed Oct 08, 2019 363  using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace;  364  using GridView = typename GetPropType::GridView;  Martin Utz committed Oct 08, 2019 365 366 367 368 369 370 371  using Element = typename GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; using NumEqVector = GetPropType; using SubControlVolume = typename FVElementGeometry::SubControlVolume; public:   Timo Koch committed Mar 30, 2020 372   Martin Utz committed Oct 08, 2019 373 This is the constructor of our problem class.  Timo Koch committed Mar 30, 2020 374   Martin Utz committed Oct 08, 2019 375 cpp  Simon Scholz committed Nov 22, 2019 376 377  RoughChannelProblem(std::shared_ptr gridGeometry) : ParentType(gridGeometry)  Martin Utz committed Oct 08, 2019 378 379  {   Timo Koch committed Mar 30, 2020 380   Martin Utz committed Oct 08, 2019 381 We read the parameters from the params.input file.  Timo Koch committed Mar 30, 2020 382   Martin Utz committed Oct 08, 2019 383 384 385 386 387 388 cpp name_ = getParam("Problem.Name"); constManningN_ = getParam("Problem.ManningN"); bedSlope_ = getParam("Problem.BedSlope"); discharge_ = getParam("Problem.Discharge");   Timo Koch committed Mar 30, 2020 389   390 We calculate the outflow boundary condition using the Gauckler-Manning-Strickler formula.  Timo Koch committed Mar 30, 2020 391   Martin Utz committed Oct 08, 2019 392 393 394 cpp hBoundary_ = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_);   Timo Koch committed Mar 30, 2020 395   Martin Utz committed Oct 08, 2019 396 We initialize the analytic solution to a verctor of the appropriate size filled with zeros.  Timo Koch committed Mar 30, 2020 397   Martin Utz committed Oct 08, 2019 398 cpp  Simon Scholz committed Nov 22, 2019 399 400  exactWaterDepth_.resize(gridGeometry->numDofs(), 0.0); exactVelocityX_.resize(gridGeometry->numDofs(), 0.0);  Martin Utz committed Oct 08, 2019 401 402  }   Timo Koch committed Mar 30, 2020 403   Martin Utz committed Oct 08, 2019 404 Get the analytical water depth  Timo Koch committed Mar 30, 2020 405   Martin Utz committed Oct 08, 2019 406 407 408 409 410 411 cpp const std::vector& getExactWaterDepth() { return exactWaterDepth_; }   Timo Koch committed Mar 30, 2020 412   Martin Utz committed Oct 08, 2019 413 Get the analytical velocity  Timo Koch committed Mar 30, 2020 414   Martin Utz committed Oct 08, 2019 415 416 417 418 419 420 cpp const std::vector& getExactVelocityX() { return exactVelocityX_; }   Timo Koch committed Mar 30, 2020 421   422 Get the water depth with Gauckler-Manning-Strickler  Timo Koch committed Mar 30, 2020 423   Martin Utz committed Oct 08, 2019 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); }   Timo Koch committed Mar 30, 2020 434   Martin Utz committed Oct 08, 2019 435 Get the analytical solution  Timo Koch committed Mar 30, 2020 436   Martin Utz committed Oct 08, 2019 437 438 439 440 441 cpp void analyticalSolution() { using std::abs;  Simon Scholz committed Nov 22, 2019 442  for (const auto& element : elements(this->gridGeometry().gridView()))  Martin Utz committed Oct 08, 2019 443 444 445 446  { const Scalar h = this->gauklerManningStrickler(discharge_,constManningN_,bedSlope_); const Scalar u = abs(discharge_)/h;  Simon Scholz committed Nov 22, 2019 447  const auto eIdx = this->gridGeometry().elementMapper().index(element);  Martin Utz committed Oct 08, 2019 448 449 450 451 452  exactWaterDepth_[eIdx] = h; exactVelocityX_[eIdx] = u; } }   Timo Koch committed Mar 30, 2020 453   Martin Utz committed Oct 08, 2019 454 Get the problem name. It is used as a prefix for files generated by the simulation.  Timo Koch committed Mar 30, 2020 455   Martin Utz committed Oct 08, 2019 456 457 458 459 460 461 cpp const std::string& name() const { return name_; }   Timo Koch committed Mar 30, 2020 462   Martin Utz committed Oct 08, 2019 463 Get the source term.  Timo Koch committed Mar 30, 2020 464   Martin Utz committed Oct 08, 2019 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);   Timo Koch committed Mar 30, 2020 474   Martin Utz committed Oct 08, 2019 475 In this model the bottom friction is the only source.  Timo Koch committed Mar 30, 2020 476   Martin Utz committed Oct 08, 2019 477 478 479 480 481 482 cpp source += bottomFrictionSource(element, fvGeometry, elemVolVars, scv); return source; }   Timo Koch committed Mar 30, 2020 483   Martin Utz committed Oct 08, 2019 484 Get the source term due to bottom friction.  Timo Koch committed Mar 30, 2020 485   Martin Utz committed Oct 08, 2019 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];   Timo Koch committed Mar 30, 2020 495   Martin Utz committed Oct 08, 2019 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).  Timo Koch committed Mar 30, 2020 497   Martin Utz committed Oct 08, 2019 498 499 500 cpp Dune::FieldVector bottomShearStress = this->spatialParams().frictionLaw(element, scv).shearStress(volVars);   Timo Koch committed Mar 30, 2020 501   Martin Utz committed Oct 08, 2019 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.  Timo Koch committed Mar 30, 2020 503   Martin Utz committed Oct 08, 2019 504 505 506 507 508 509 510 511 cpp bottomFrictionSource[0] = 0.0; bottomFrictionSource[1] = bottomShearStress[0]; bottomFrictionSource[2] = bottomShearStress[1]; return bottomFrictionSource; }   Timo Koch committed Mar 30, 2020 512   Martin Utz committed Oct 08, 2019 513 We specify the boundary condition type.  Timo Koch committed Mar 30, 2020 514   Martin Utz committed Oct 08, 2019 515 516 517 518 519 cpp BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const { BoundaryTypes bcTypes;   Timo Koch committed Mar 30, 2020 520   Martin Utz committed Oct 08, 2019 521 Since we use a weak imposition all boundary conditions are of Neumann type.  Timo Koch committed Mar 30, 2020 522   Martin Utz committed Oct 08, 2019 523 524 525 526 527 cpp bcTypes.setAllNeumann(); return bcTypes; }   Timo Koch committed Mar 30, 2020 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 committed Oct 08, 2019 535 536 537 538 cpp NeumannFluxes neumann(const Element& element, const FVElementGeometry& fvGeometry, const ElementVolumeVariables& elemVolVars,  Katharina Heck committed Oct 08, 2019 539  const ElementFluxVariablesCache& elemFluxVarsCache,  Martin Utz committed Oct 08, 2019 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 boundaryStateVariables;   Timo Koch committed Mar 30, 2020 550   Martin Utz committed Oct 08, 2019 551 Calculate the rieman invariants for imposed discharge at the left side.  Timo Koch committed Mar 30, 2020 552   Martin Utz committed Oct 08, 2019 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); }   Timo Koch committed Mar 30, 2020 564   Martin Utz committed Oct 08, 2019 565 Calculate the rieman invariants for impose water depth at the right side.  Timo Koch committed Mar 30, 2020 566   Martin Utz committed Oct 08, 2019 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); }   Timo Koch committed Mar 30, 2020 578   Martin Utz committed Oct 08, 2019 579 Calculate the rieman invarianty for the no-flow boundary.  Timo Koch committed Mar 30, 2020 580   Martin Utz committed Oct 08, 2019 581 582 583 584 585 586 587 588 cpp else { boundaryStateVariables[0] = insideVolVars.waterDepth(); boundaryStateVariables[1] = -insideVolVars.velocity(0); boundaryStateVariables[2] = -insideVolVars.velocity(1); }   Timo Koch committed Mar 30, 2020 589   Martin Utz committed Oct 08, 2019 590 We calculate the boundary fluxes based on a Riemann problem.  Timo Koch committed Mar 30, 2020 591   Martin Utz committed Oct 08, 2019 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; }   Timo Koch committed Mar 30, 2020 611   Martin Utz committed Oct 08, 2019 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.  Timo Koch committed Mar 30, 2020 613   Martin Utz committed Oct 08, 2019 614 615 616 617 618 cpp PrimaryVariables initialAtPos(const GlobalPosition &globalPos) const { PrimaryVariables values(0.0);   Timo Koch committed Mar 30, 2020 619   Martin Utz committed Oct 08, 2019 620 We set the initial water depth to one meter.  Timo Koch committed Mar 30, 2020 621   Martin Utz committed Oct 08, 2019 622 623 624 cpp values[0] = 1.0;   Timo Koch committed Mar 30, 2020 625   Martin Utz committed Oct 08, 2019 626 We set the x-component of the initial velocity to zero.  Timo Koch committed Mar 30, 2020 627   Martin Utz committed Oct 08, 2019 628 629 630 cpp values[1] = 0.0;   Timo Koch committed Mar 30, 2020 631   Martin Utz committed Oct 08, 2019 632 We set the y-component of the initial velocity to zero.  Timo Koch committed Mar 30, 2020 633   Martin Utz committed Oct 08, 2019 634 635 636 637 638 639 cpp values[2] = 0.0; return values; };   Timo Koch committed Mar 30, 2020 640   Martin Utz committed Oct 08, 2019 641 \}  Timo Koch committed Mar 30, 2020 642   Martin Utz committed Oct 08, 2019 643 644 645 646 cpp private:   Timo Koch committed Mar 30, 2020 647   648 We declare the private variables of the problem. They are initialized in the problems constructor.  Martin Utz committed Oct 08, 2019 649 We declare the variable for the analytic solution.  Timo Koch committed Mar 30, 2020 650   Martin Utz committed Oct 08, 2019 651 652 653 654 cpp std::vector exactWaterDepth_; std::vector exactVelocityX_;   Timo Koch committed Mar 30, 2020 655   Martin Utz committed Oct 08, 2019 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.  Timo Koch committed Mar 30, 2020 657   Martin Utz committed Oct 08, 2019 658 659 660 cpp Scalar constManningN_;   Timo Koch committed Mar 30, 2020 661   Martin Utz committed Oct 08, 2019 662 The constant bed slope.  Timo Koch committed Mar 30, 2020 663   Martin Utz committed Oct 08, 2019 664 665 666 cpp Scalar bedSlope_;   Timo Koch committed Mar 30, 2020 667   Martin Utz committed Oct 08, 2019 668 The water depth at the outflow boundary.  Timo Koch committed Mar 30, 2020 669   Martin Utz committed Oct 08, 2019 670 671 672 cpp Scalar hBoundary_;   Timo Koch committed Mar 30, 2020 673   Martin Utz committed Oct 08, 2019 674 The discharge at the inflow boundary.  Timo Koch committed Mar 30, 2020 675   Martin Utz committed Oct 08, 2019 676 677 678 cpp Scalar discharge_;   Timo Koch committed Mar 30, 2020 679   Martin Utz committed Oct 08, 2019 680 eps is used as a small value for the definition of the boundry conditions  Timo Koch committed Mar 30, 2020 681   Martin Utz committed Oct 08, 2019 682 683 684 685 cpp static constexpr Scalar eps_ = 1.0e-6; std::string name_; };  Timo Koch committed Mar 30, 2020 686   Timo Koch committed Apr 02, 2020 687 } // end namespace Dumux  Martin Utz committed Oct 08, 2019 688 689 690 691   Timo Koch committed Mar 30, 2020 692   Martin Utz committed Oct 08, 2019 693 694 ## The file main.cc  Timo Koch committed Mar 30, 2020 695   Martin Utz committed Oct 08, 2019 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  Timo Koch committed Mar 30, 2020 698   Martin Utz committed Oct 08, 2019 699 700 701 cpp #include   Timo Koch committed Mar 30, 2020 702   Martin Utz committed Oct 08, 2019 703 Standard header file for C++, to get time and date information.  Timo Koch committed Mar 30, 2020 704   Martin Utz committed Oct 08, 2019 705 706 707 cpp #include   Timo Koch committed Mar 30, 2020 708   Martin Utz committed Oct 08, 2019 709 Standard header file for C++, for in- and output.  Timo Koch committed Mar 30, 2020 710   Martin Utz committed Oct 08, 2019 711 712 713 cpp #include   Timo Koch committed Mar 30, 2020 714   Martin Utz committed Oct 08, 2019 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.  Timo Koch committed Mar 30, 2020 716   Martin Utz committed Oct 08, 2019 717 718 719 720 721 722 cpp #include #include #include #include   Timo Koch committed Mar 30, 2020 723   Martin Utz committed Oct 08, 2019 724 We need the following class to simplify the writing of dumux simulation data to VTK format.  Timo Koch committed Mar 30, 2020 725   Martin Utz committed Oct 08, 2019 726 727 728 cpp #include   Timo Koch committed Mar 30, 2020 729   Martin Utz committed Oct 08, 2019 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.  Timo Koch committed Mar 30, 2020 731   Martin Utz committed Oct 08, 2019 732 733 734 cpp #include   Timo Koch committed Mar 30, 2020 735   Martin Utz committed Oct 08, 2019 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.  Timo Koch committed Mar 30, 2020 737   Martin Utz committed Oct 08, 2019 738 739 740 cpp #include   Timo Koch committed Mar 30, 2020 741   Martin Utz committed Oct 08, 2019 742 The file dumuxmessage.hh contains the class defining the start and end message of the simulation.  Timo Koch committed Mar 30, 2020 743   Martin Utz committed Oct 08, 2019 744 745 746 747 cpp #include #include   Timo Koch committed Mar 30, 2020 748   Martin Utz committed Oct 08, 2019 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.  Timo Koch committed Mar 30, 2020 750   Martin Utz committed Oct 08, 2019 751 752 753 cpp #include   Timo Koch committed Mar 30, 2020 754   Martin Utz committed Oct 08, 2019 755 We include the linear solver to be used to solve the linear system  Timo Koch committed Mar 30, 2020 756   Martin Utz committed Oct 08, 2019 757 758 cpp #include  Timo Koch committed Mar 17, 2020 759 #include  Martin Utz committed Oct 08, 2019 760   Timo Koch committed Mar 30, 2020 761   Martin Utz committed Oct 08, 2019 762 We include the nonlinear newtons method  Timo Koch committed Mar 30, 2020 763   Martin Utz committed Oct 08, 2019 764 765 766 cpp #include   Timo Koch committed Mar 30, 2020 767   Martin Utz committed Oct 08, 2019 768 Further we include assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation)  Timo Koch committed Mar 30, 2020 769   Martin Utz committed Oct 08, 2019 770 771 772 cpp #include   Timo Koch committed Mar 30, 2020 773   Timo Koch committed Apr 02, 2020 774 We include the properties  Timo Koch committed Mar 30, 2020 775   Martin Utz committed Oct 08, 2019 776 cpp  Timo Koch committed Apr 02, 2020 777 #include "properties.hh"  Martin Utz committed Oct 08, 2019 778   Timo Koch committed Mar 30, 2020 779   Martin Utz committed Oct 08, 2019 780 ### Beginning of the main function  Timo Koch committed Mar 30, 2020 781   Martin Utz committed Oct 08, 2019 782 783 784 785 786 cpp int main(int argc, char** argv) try { using namespace Dumux;   Timo Koch committed Mar 30, 2020 787   Martin Utz committed Oct 08, 2019 788 We define the type tag for this problem  Timo Koch committed Mar 30, 2020 789   Martin Utz committed Oct 08, 2019 790 791 792 cpp using TypeTag = Properties::TTag::RoughChannel;   Timo Koch committed Mar 30, 2020 793   Martin Utz committed Oct 08, 2019 794 We initialize MPI, finalize is done automatically on exit  Timo Koch committed Mar 30, 2020 795   Martin Utz committed Oct 08, 2019 796 797 798 cpp const auto& mpiHelper = Dune::MPIHelper::instance(argc, argv);   Timo Koch committed Mar 30, 2020 799   Martin Utz committed Oct 08, 2019 800 We print dumux start message  Timo Koch committed Mar 30, 2020 801   Martin Utz committed Oct 08, 2019 802 803 804 805 cpp if (mpiHelper.rank() == 0) DumuxMessage::print(/*firstCall=*/true);   Timo Koch committed Mar 30, 2020 806   Martin Utz committed Oct 08, 2019 807 We parse command line arguments and input file  Timo Koch committed Mar 30, 2020 808   Martin Utz committed Oct 08, 2019 809 810 811 cpp Parameters::init(argc, argv);   Timo Koch committed Mar 30, 2020 812   Martin Utz committed Oct 08, 2019 813 814 ### Create the grid A gridmanager tries to create the grid either from a grid file or the input file.  Timo Koch committed Mar 30, 2020 815   Martin Utz committed Oct 08, 2019 816 817 818 819 cpp GridManager> gridManager; gridManager.init();   Timo Koch committed Mar 30, 2020 820   Martin Utz committed Oct 08, 2019 821 We compute on the leaf grid view  Timo Koch committed Mar 30, 2020 822   Martin Utz committed Oct 08, 2019 823 824 825 cpp const auto& leafGridView = gridManager.grid().leafGridView();   Timo Koch committed Mar 30, 2020 826   Martin Utz committed Oct 08, 2019 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.  Timo Koch committed Mar 30, 2020 831   Martin Utz committed Oct 08, 2019 832 cpp  Simon Scholz committed Nov 22, 2019 833 834 835  using GridGeometry = GetPropType; auto gridGeometry = std::make_shared(leafGridView); gridGeometry->update();  Martin Utz committed Oct 08, 2019 836   Timo Koch committed Mar 30, 2020 837   Martin Utz committed Oct 08, 2019 838 In the problem, we define the boundary and initial conditions.  Timo Koch committed Mar 30, 2020 839   Martin Utz committed Oct 08, 2019 840 841 cpp using Problem = GetPropType;  Simon Scholz committed Nov 22, 2019 842  auto problem = std::make_shared(gridGeometry);  Martin Utz committed Oct 08, 2019 843   Timo Koch committed Mar 30, 2020 844   Martin Utz committed Oct 08, 2019 845 We initialize the solution vector  Timo Koch committed Mar 30, 2020 846   Martin Utz committed Oct 08, 2019 847 848 cpp using SolutionVector = GetPropType;  Simon Scholz committed Nov 22, 2019 849  SolutionVector x(gridGeometry->numDofs());  Martin Utz committed Oct 08, 2019 850 851 852  problem->applyInitialSolution(x); auto xOld = x;   Timo Koch committed Mar 30, 2020 853   Martin Utz committed Oct 08, 2019 854 And then use the solutionvector to intialize the gridVariables.  Timo Koch committed Mar 30, 2020 855   Martin Utz committed Oct 08, 2019 856 857 cpp using GridVariables = GetPropType;  Simon Scholz committed Nov 22, 2019 858  auto gridVariables = std::make_shared(problem, gridGeometry);  Martin Utz committed Oct 08, 2019 859 860  gridVariables->init(x);   Timo Koch committed Mar 30, 2020 861   Martin Utz committed Oct 08, 2019 862 We get some time loop parameters from the input file.  Timo Koch committed Mar 30, 2020 863   Martin Utz committed Oct 08, 2019 864 865 866 867 868 869 cpp using Scalar = GetPropType; const auto tEnd = getParam("TimeLoop.TEnd"); const auto maxDt = getParam("TimeLoop.MaxTimeStepSize"); auto dt = getParam("TimeLoop.DtInitial");   Timo Koch committed Mar 30, 2020 870   Martin Utz committed Oct 08, 2019 871 We intialize the vtk output module. Each model has a predefined model specific output with relevant parameters for that model.  Timo Koch committed Mar 30, 2020 872   Martin Utz committed Oct 08, 2019 873 874 875 876 cpp using IOFields = GetPropType; VtkOutputModule vtkWriter(*gridVariables,x, problem->name());   Timo Koch committed Mar 30, 2020 877   Martin Utz committed Oct 08, 2019 878 We add the analytical solution ("exactWaterDepth" and "exactVelocityX") to the predefined specific output.  Timo Koch committed Mar 30, 2020 879   Martin Utz committed Oct 08, 2019 880 881 882 883 cpp vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth"); vtkWriter.addField(problem->getExactVelocityX(), "exactVelocityX");   Timo Koch committed Mar 30, 2020 884   Martin Utz committed Oct 08, 2019 885 We calculate the analytic solution.  Timo Koch committed Mar 30, 2020 886   Martin Utz committed Oct 08, 2019 887 888 889 890 891 cpp problem->analyticalSolution(); IOFields::initOutputModule(vtkWriter); vtkWriter.write(0.0);   Timo Koch committed Mar 30, 2020 892   Martin Utz committed Oct 08, 2019 893 We instantiate time loop.  Timo Koch committed Mar 30, 2020 894   Martin Utz committed Oct 08, 2019 895 896 897 898 cpp auto timeLoop = std::make_shared>(0, dt, tEnd); timeLoop->setMaxTimeStepSize(maxDt);   Timo Koch committed Mar 30, 2020 899   Martin Utz committed Oct 08, 2019 900 we set the assembler with the time loop because we have an instationary problem.  Timo Koch committed Mar 30, 2020 901   Martin Utz committed Oct 08, 2019 902 903 cpp using Assembler = FVAssembler;  Timo Koch committed Apr 14, 2020 904  auto assembler = std::make_shared(problem, gridGeometry, gridVariables, timeLoop, xOld);  Martin Utz committed Oct 08, 2019 905   Timo Koch committed Mar 30, 2020 906   Martin Utz committed Oct 08, 2019 907 We set the linear solver.  Timo Koch committed Mar 30, 2020 908   Martin Utz committed Oct 08, 2019 909 cpp  Timo Koch committed Mar 17, 2020 910  using LinearSolver = AMGBiCGSTABBackend>;  Simon Scholz committed Nov 22, 2019 911  auto linearSolver = std::make_shared(leafGridView, gridGeometry->dofMapper());  Martin Utz committed Oct 08, 2019 912   Timo Koch committed Mar 30, 2020 913   Martin Utz committed Oct 08, 2019 914 Additionaly, we set the non-linear solver.  Timo Koch committed Mar 30, 2020 915   Martin Utz committed Oct 08, 2019 916 917 918 919 cpp using NewtonSolver = Dumux::NewtonSolver; NewtonSolver nonLinearSolver(assembler, linearSolver);   Timo Koch committed Mar 30, 2020 920   Martin Utz committed Oct 08, 2019 921 We set some check point at the end of the time loop. The check point is used to trigger the vtk output.  Timo Koch committed Mar 30, 2020 922   Martin Utz committed Oct 08, 2019 923 924 925 cpp timeLoop->setCheckPoint(tEnd);   Timo Koch committed Mar 30, 2020 926   Martin Utz committed Oct 08, 2019 927 We start the time loop.  Timo Koch committed Mar 30, 2020 928   Martin Utz committed Oct 08, 2019 929 930 931 932 cpp timeLoop->start(); do {   Timo Koch committed Mar 30, 2020 933   Martin Utz committed Oct 08, 2019 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.  Timo Koch committed Mar 30, 2020 935   Martin Utz committed Oct 08, 2019 936 937 938 cpp assembler->setPreviousSolution(xOld);   Timo Koch committed Mar 30, 2020 939   Martin Utz committed Oct 08, 2019 940 We solve the non-linear system with time step control.  Timo Koch committed Mar 30, 2020 941   Martin Utz committed Oct 08, 2019 942 943 944 cpp nonLinearSolver.solve(x,*timeLoop);   Timo Koch committed Mar 30, 2020 945   Martin Utz committed Oct 08, 2019 946 We make the new solution the old solution.  Timo Koch committed Mar 30, 2020 947   Martin Utz committed Oct 08, 2019 948 949 950 951 cpp xOld = x; gridVariables->advanceTimeStep();   Timo Koch committed Mar 30, 2020 952   Martin Utz committed Oct 08, 2019 953 We advance to the time loop to the next step.  Timo Koch committed Mar 30, 2020 954   Martin Utz committed Oct 08, 2019 955 956 957 cpp timeLoop->advanceTimeStep();   Timo Koch committed Mar 30, 2020 958   Martin Utz committed Oct 08, 2019 959 We write vtk output, if we reached the check point (end of time loop)  Timo Koch committed Mar 30, 2020 960   Martin Utz committed Oct 08, 2019 961 962 963 964 cpp if (timeLoop->isCheckPoint()) vtkWriter.write(timeLoop->time());   Timo Koch committed Mar 30, 2020 965   Martin Utz committed Oct 08, 2019 966 We report statistics of this time step.  Timo Koch committed Mar 30, 2020 967   Martin Utz committed Oct 08, 2019 968 969 970 cpp timeLoop->reportTimeStep();   Timo Koch committed Mar 30, 2020 971   Martin Utz committed Oct 08, 2019 972 We set new dt as suggested by newton controller for the next time step.  Timo Koch committed Mar 30, 2020 973   Martin Utz committed Oct 08, 2019 974 975 976 977 978 979 980 981 cpp timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize())); } while (!timeLoop->finished()); timeLoop->finalize(leafGridView.comm());   Timo Koch committed Mar 30, 2020 982   Martin Utz committed Oct 08, 2019 983 984 ### Final Output We print dumux end message.  Timo Koch committed Mar 30, 2020 985   Martin Utz committed Oct 08, 2019 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; }