GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/examples/shallowwaterfriction/main.cc
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 51 64 79.7%
Functions: 1 1 100.0%
Branches: 111 262 42.4%
Line Branch Exec Source
1 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2 // vi: set et ts=4 sw=4 sts=4:
3 //
4 // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7
8 // ## The main file `main.cc`
9 // [[content]]
10 //
11 // ### Included header files
12 // [[details]] includes
13 // [[exclude]]
14 // Some generic includes.
15 #include <config.h>
16 #include <iostream>
17 // [[/exclude]]
18 //
19 // DUNE helper class for MPI
20 #include <dune/common/parallel/mpihelper.hh>
21
22 // The following headers include functionality related to property definition or retrieval, as well as
23 // the retrieval of input parameters specified in the input file or via the command line.
24 #include <dumux/common/properties.hh>
25 #include <dumux/common/parameters.hh>
26 #include <dumux/common/initialize.hh>
27
28 // The following files contains the available linear solver backends, the non linear Newton Solver
29 // and the assembler for the linear systems arising from finite volume discretizations
30 // (box-scheme, tpfa-approximation, mpfa-approximation).
31 #include <dumux/linear/istlsolvers.hh>
32 #include <dumux/linear/linearsolvertraits.hh>
33 #include <dumux/linear/linearalgebratraits.hh>
34 #include <dumux/nonlinear/newtonsolver.hh>
35 #include <dumux/assembly/fvassembler.hh>
36
37 // The following class provides a convenient way of writing of dumux simulation results to VTK format.
38 #include <dumux/io/vtkoutputmodule.hh>
39 // The gridmanager constructs a grid from the information in the input or grid file.
40 // Many different Dune grid implementations are supported, of which a list can be found
41 // in `gridmanager.hh`.
42 #include <dumux/io/grid/gridmanager_yasp.hh>
43
44 // We include the header file specifying the properties of this example
45 #include "properties.hh"
46 // [[/details]]
47 //
48
49 // ### The main function
50 // We will now discuss the main program flow implemented within the `main` function.
51 // At the beginning of each program we initialize (e.g. potential parallel backends like MPI).
52 // Moreover, we parse the run-time arguments from the command line and the
53 // input file:
54 // [[codeblock]]
55 1 int main(int argc, char** argv) try
56 {
57 1 using namespace Dumux;
58
59 // maybe initialize MPI and/or multithreading backend
60
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Dumux::initialize(argc, argv);
61
62 // We parse command line arguments and input file
63
9/28
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 1 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 1 times.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✗ Branch 30 not taken.
✗ Branch 31 not taken.
✗ Branch 32 not taken.
 4 Parameters::init(argc, argv);
64 // [[/codeblock]]
65
66 // We define a convenience alias for the type tag of the problem. The type
67 // tag contains all the properties that are needed to define the model and the problem
68 // setup. Throughout the main file, we will obtain types defined for these type tag
69 // using the property system, i.e. with `GetPropType`.
70 1 using TypeTag = Properties::TTag::RoughChannel;
71
72 // #### Step 1: Create the grid
73 // The `GridManager` class creates the grid from information given in the input file.
74 // This can either be a grid file or, in the case of structured grids, one can specify the coordinates
75 // of the corners of the grid and the number of cells to be used to discretize each spatial direction.
76 //[[codeblock]]
77
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
78
5/12
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 2 gridManager.init();
79
80 // We compute on the leaf grid view
81
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 1 const auto& leafGridView = gridManager.grid().leafGridView();
82 //[[/codeblock]]
83
84 // #### Step 2: Solving the shallow water problem
85 // First, a finite volume grid geometry is constructed from the grid that was created above.
86 // This builds the sub-control volumes (scv) and sub-control volume faces (scvf) for each element
87 // of the grid partition.
88 1 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
89
1/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2 auto gridGeometry = std::make_shared<GridGeometry>(leafGridView);
90
91 // We now instantiate the problem, in which we define the boundary and initial conditions.
92 1 using Problem = GetPropType<TypeTag, Properties::Problem>;
93
2/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 2 auto problem = std::make_shared<Problem>(gridGeometry);
94
95 // We initialize the solution vector. The shallow water problem is transient,
96 // therefore the initial solution defined in the problem is applied to the
97 // solution vector. On the basis of this solution, we initialize then the grid variables.
98 1 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
99
4/10
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 3 SolutionVector x(gridGeometry->numDofs());
100
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 problem->applyInitialSolution(x);
101
2/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 2 auto xOld = x;
102
103 1 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
104
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 2 auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
105
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 gridVariables->init(x);
106
107 // Let us now instantiate the time loop. Therefore, we read in some time loop parameters from the input file.
108 // The parameter `tEnd` defines the duration of the simulation, `dt` the initial time step size and `maxDt` the maximal time step size.
109 // Moreover, we define the end of the simulation `tEnd` as check point in the time loop at which we will write the solution to vtk files.
110 // [[codeblock]]
111 1 using Scalar = GetPropType<TypeTag, Properties::Scalar>; // type for scalar values
112
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
113
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
114
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 auto dt = getParam<Scalar>("TimeLoop.DtInitial");
115
116 // We instantiate time loop.
117
2/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 2 auto timeLoop = std::make_shared<CheckPointTimeLoop<Scalar>>(0.0, dt, tEnd);
118
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 timeLoop->setMaxTimeStepSize(maxDt);
119
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 timeLoop->setCheckPoint(tEnd);
120 // [[/codeblock]]
121
122 // We initialize the assembler with a time loop for the transient problem.
123 // Within the time loop, we will use this assembler in each time step to assemble the linear system.
124 1 using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
125
2/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 2 auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables, timeLoop, xOld);
126
127 // We initialize the linear solver.
128 1 using LinearSolver = AMGBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>,
129 LinearAlgebraTraitsFromAssembler<Assembler>>;
130
4/10
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 4 auto linearSolver = std::make_shared<LinearSolver>(leafGridView, gridGeometry->dofMapper());
131
132 // We initialize the non-linear solver.
133 1 using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
134
8/20
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 22 taken 1 times.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 27 not taken.
 5 NewtonSolver nonLinearSolver(assembler, linearSolver);
135
136 // The following lines of code initialize the vtk output module, add the velocity output facility
137 // and write out the initial solution. At each checkpoint, we will use the output module to write
138 // the solution of a time step into a corresponding vtk file.
139 // [[codeblock]]
140
7/14
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
 3 VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables,x, problem->name());
141
142 // add model-specific output fields to the writer
143 1 using IOFields = GetPropType<TypeTag, Properties::IOFields>;
144
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 IOFields::initOutputModule(vtkWriter);
145
146 // We add the analytical solution ("exactWaterDepth" and "exactVelocityX") to the predefined specific output.
147
7/16
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 2 vtkWriter.addField(problem->getExactWaterDepth(), "exactWaterDepth");
148
7/16
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 1 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 1 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 1 times.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 2 vtkWriter.addField(problem->getExactVelocityX(), "exactVelocityX");
149
150 // We calculate the analytic solution.
151
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 problem->analyticalSolution();
152
153 // write initial solution (including the above calculated analytical solution.
154
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vtkWriter.write(0.0);
155 // [[/codeblock]]
156
157 // ##### The time loop
158 // We start the time loop and solve a new time step as long as `tEnd` is not reached. In every time step,
159 // the problem is assembled and solved, the solution is updated, and when a checkpoint is reached the solution
160 // is written to a new vtk file. In addition, statistics related to CPU time, the current simulation time
161 // and the time step sizes used is printed to the terminal.
162 // [[codeblock]]
163
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 timeLoop->start(); do
164 {
165 // We solve the non-linear system with time step control, using Newthon's method.
166
2/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
 48 nonLinearSolver.solve(x,*timeLoop);
167
168 // We make the new solution the old solution.
169
1/2
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
 24 xOld = x;
170
2/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
 48 gridVariables->advanceTimeStep();
171
172 // We advance to the time loop to the next step.
173
2/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
 48 timeLoop->advanceTimeStep();
174
175 // We write vtk output, if we reached the check point (end of time loop)
176
4/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 23 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 23 times.
 48 if (timeLoop->isCheckPoint())
177
3/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
 3 vtkWriter.write(timeLoop->time());
178
179 // We report statistics of this time step.
180
2/4
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
 48 timeLoop->reportTimeStep();
181
182 // We set new dt as suggested by newton controller for the next time step.
183
5/10
✗ Branch 0 not taken.
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 24 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 24 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 24 times.
✓ Branch 9 taken 24 times.
✗ Branch 10 not taken.
 120 timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
184
185
186
4/6
✓ Branch 1 taken 24 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 24 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 23 times.
✓ Branch 7 taken 1 times.
 48 } while (!timeLoop->finished());
187 // [[/codeblock]]
188
189 // The following piece of code prints a final status report of the time loop
190 // before the program is terminated.
191
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 2 timeLoop->finalize(leafGridView.comm());
192
193 1 return 0;
194 }
195 // #### Exception handling
196 // In this part of the main file we catch and print possible exceptions that could
197 // occur during the simulation.
198 // [[details]] exception handler
199 // [[codeblock]]
200 // errors related to run-time parameters
201 ✗ catch (const Dumux::ParameterException &e)
202 {
203 ✗ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
204 ✗ return 1;
205 }
206 // errors related to the parsing of Dune grid files
207 ✗ catch (const Dune::DGFException & e)
208 {
209 ✗ std::cerr << "DGF exception thrown (" << e <<
210 "). Most likely, the DGF file name is wrong "
211 "or the DGF file is corrupted, "
212 "e.g. missing hash at end of file or wrong number (dimensions) of entries."
213 ✗ << " ---> Abort!" << std::endl;
214 ✗ return 2;
215 }
216 // generic error handling with Dune::Exception
217 ✗ catch (const Dune::Exception &e)
218 {
219 ✗ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
220 ✗ return 3;
221 }
222 // other exceptions
223 ✗ catch (...)
224 {
225 ✗ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
226 ✗ return 4;
227 }
228 // [[/codeblock]]
229 // [[/details]]
230 // [[/content]]
231