GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/examples/freeflowchannel/main.cc
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 64 77 83.1%
Functions: 1 1 100.0%
Branches: 125 298 41.9%
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 <ctime>
17 #include <iostream>
18 // [[/exclude]]
19
20 // These are DUNE helper classes related to parallel computations and file I/O.
21 #include <dune/common/parallel/mpihelper.hh>
22 #include <dune/grid/io/file/dgfparser/dgfexception.hh>
23
24 // The following headers include functionality related to property definition or retrieval, as well as
25 // the retrieval of input parameters specified in the input file or via the command line.
26 #include <dumux/common/properties.hh>
27 #include <dumux/common/parameters.hh>
28 #include <dumux/common/initialize.hh>
29
30 // The following files contain the non-linear Newton solver, the available linear solver backends and the assembler for the linear
31 // systems arising from the staggered-grid discretization.
32 #include <dumux/nonlinear/newtonsolver.hh>
33 #include <dumux/linear/istlsolvers.hh>
34 #include <dumux/linear/linearalgebratraits.hh>
35 #include <dumux/linear/linearsolvertraits.hh>
36 #include <dumux/assembly/staggeredfvassembler.hh>
37 #include <dumux/assembly/diffmethod.hh> // analytic or numeric differentiation
38
39 // The following class provides a convenient way of writing of dumux simulation results to VTK format.
40 #include <dumux/io/staggeredvtkoutputmodule.hh>
41 // The gridmanager constructs a grid from the information in the input or grid file.
42 // Many different Dune grid implementations are supported, of which a list can be found
43 // in `gridmanager.hh`.
44 #include <dumux/io/grid/gridmanager.hh>
45
46 // This class contains functionality for additional flux output.
47 #include <dumux/freeflow/navierstokes/staggered/fluxoversurface.hh>
48
49
50 // In this header, a `TypeTag` is defined, which collects
51 // the properties that are required for the simulation.
52 // It also contains the actual problem with initial and boundary conditions.
53 // For detailed information, please have a look
54 // at the documentation provided therein.
55 #include "properties.hh"
56 // [[/details]]
57 //
58
59 // ### The main function
60 // We will now discuss the main program flow implemented within the `main` function.
61 // At the beginning of each program using Dune, an instance of `Dune::MPIHelper` has to
62 // be created. Moreover, we parse the run-time arguments from the command line and the
63 // input file:
64 // [[codeblock]]
65 1 int main(int argc, char** argv) try
66 {
67 1 using namespace Dumux;
68
69 // maybe initialize MPI and/or multithreading backend
70
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Dumux::initialize(argc, argv);
71
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const auto& mpiHelper = Dune::MPIHelper::instance();
72
73 // parse command line arguments and input file
74
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);
75 // [[/codeblock]]
76
77 // We define a convenience alias for the type tag of the problem (`Properties::TTag::ChannelExample`). The type
78 // tag contains all the properties that are needed to define the model and the problem
79 // setup. Throughout the main file, we will obtain types defined for this type tag
80 // using the property system, i.e. with `GetPropType`.
81 1 using TypeTag = Properties::TTag::ChannelExample;
82
83 // #### Step 1: Create the grid
84 // The `GridManager` class creates the grid from information given in the input file.
85 // This can either be a grid file, or in the case of structured grids, one can specify the coordinates
86 // of the corners of the grid and the number of cells to be used to discretize each spatial direction.
87 // [[codeblock]]
88
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
89
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();
90
91 // We compute on the leaf grid view.
92
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();
93 // [[/codeblock]]
94
95 // #### Step 2: Setting up and solving the problem
96 // First, a finite volume grid geometry is constructed from the grid that was created above.
97 // This builds the sub-control volumes (`scv`) and sub-control volume faces (`scvf`) for each element
98 // of the grid partition.
99 1 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
100
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);
101
102 // We now instantiate the `problem`, in which we define the boundary and initial conditions.
103 1 using Problem = GetPropType<TypeTag, Properties::Problem>;
104
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);
105
106 // We set a solution vector `x` which consist of two parts: one part (indexed by `cellCenterIdx`)
107 // is for the pressure degrees of freedom (`dofs`) living in grid cell centers. Another part
108 // (indexed by `faceIdx`) is for degrees of freedom defining the normal velocities on grid cell faces.
109 // We initialize the solution vector by what was defined as the initial solution of the the problem.
110 1 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
111
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 2 SolutionVector x;
112
4/8
✓ 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.
 3 x[GridGeometry::cellCenterIdx()].resize(gridGeometry->numCellCenterDofs());
113
4/8
✓ 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.
 3 x[GridGeometry::faceIdx()].resize(gridGeometry->numFaceDofs());
114
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 problem->applyInitialSolution(x);
115
116 // The grid variables are used store variables (primary and secondary variables) on sub-control volumes and faces (volume and flux variables).
117 1 using GridVariables =GetPropType<TypeTag, Properties::GridVariables>;
118
1/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2 auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
119
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 gridVariables->init(x);
120
121 // We then initialize the predefined model-specific output VTK output.
122 1 using IOFields = GetPropType<TypeTag, Properties::IOFields>;
123
8/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 16 taken 1 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 1 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 1 times.
✗ Branch 22 not taken.
 4 StaggeredVtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, x, problem->name());
124
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 IOFields::initOutputModule(vtkWriter); // Add model specific output fields
125
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vtkWriter.write(0.0);
126
127 // [[details]] calculation of surface fluxes
128 // We set up two surfaces over which fluxes are calculated.
129 // We determine the extend $`[xMin,xMax] \times [yMin,yMax]`$ of the physical domain.
130 // The first surface (added by the first call of `addSurface`) shall be placed at the middle of the channel.
131 // If we have an odd number of cells in x-direction, there would not be any cell faces
132 // at the position of the surface (which is required for the flux calculation).
133 // In this case, we add half a cell-width to the x-position in order to make sure that
134 // the cell faces lie on the surface. This assumes a regular cartesian grid.
135 // The second surface (second call of `addSurface`) is placed at the outlet of the channel.
136 1 FluxOverSurface<GridVariables,
137 SolutionVector,
138 GetPropType<TypeTag, Properties::ModelTraits>,
139
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 3 GetPropType<TypeTag, Properties::LocalResidual>> flux(*gridVariables, x);
140
141 1 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
142
143
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 const Scalar xMin = gridGeometry->bBoxMin()[0];
144
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 const Scalar xMax = gridGeometry->bBoxMax()[0];
145
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 const Scalar yMin = gridGeometry->bBoxMin()[1];
146
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 const Scalar yMax = gridGeometry->bBoxMax()[1];
147
148 1 const Scalar planePosMiddleX = xMin + 0.5*(xMax - xMin);
149
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 int numCellsX = getParam<std::vector<int>>("Grid.Cells")[0];
150
151
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const unsigned int refinement = getParam<unsigned int>("Grid.Refinement", 0);
152 1 numCellsX *= (1<<refinement);
153
154
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 const Scalar offsetX = (numCellsX % 2 == 0) ? 0.0 : 0.5*((xMax - xMin) / numCellsX);
155
156 1 using GridView = typename GridGeometry::GridView;
157 1 using Element = typename GridView::template Codim<0>::Entity;
158 1 using GlobalPosition = typename Element::Geometry::GlobalCoordinate;
159
160
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const auto p0middle = GlobalPosition{planePosMiddleX + offsetX, yMin};
161
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const auto p1middle = GlobalPosition{planePosMiddleX + offsetX, yMax};
162
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 flux.addSurface("middle", p0middle, p1middle);
163
164
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const auto p0outlet = GlobalPosition{xMax, yMin};
165
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const auto p1outlet = GlobalPosition{xMax, yMax};
166
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 flux.addSurface("outlet", p0outlet, p1outlet);
167 // [[/details]]
168
169 // We create and initialize the `assembler` for the stationary problem.
170 // This is where the Jacobian matrix for the Newton solver is assembled.
171 1 using Assembler = StaggeredFVAssembler<TypeTag, DiffMethod::numeric>;
172
1/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 2 auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
173
174 // We use UMFPack as direct linear solver within each Newton iteration.
175 1 using LinearSolver = Dumux::UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
176
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 linearSolver = std::make_shared<LinearSolver>();
177
178 // This example considers a linear problem (incompressible Stokes flow), therefore
179 // the non-linear Newton solver is not really necessary.
180 // For sake of generality, we nevertheless use it here such that the example can be easily
181 // changed to a non-linear problem by switching on the inertia terms in the input file or by choosing a compressible fluid.
182 // In the following piece of code we instantiate the non-linear newton solver and let it solve
183 // the problem.
184 // [[codeblock]]
185 // alias for and instantiation of the newton solver
186 1 using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
187
12/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 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 1 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 1 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 1 times.
✓ Branch 29 taken 1 times.
✗ Branch 30 not taken.
✓ Branch 31 taken 1 times.
✗ Branch 32 not taken.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
 5 NewtonSolver nonLinearSolver(assembler, linearSolver);
188
189 // Solve the (potentially non-linear) system.
190
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 nonLinearSolver.solve(x);
191 // [[/codeblock]]
192 // In the following we calculate mass and volume fluxes over the planes specified above
193 // (you have to click to unfold the code showing how to set up the surface fluxes).
194
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 flux.calculateMassOrMoleFluxes();
195
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 flux.calculateVolumeFluxes();
196
197 // #### Final Output
198 // We write the VTK output and print the mass/energy/volume fluxes over the planes.
199 // We conclude by printing the dumux end message.
200 // [[codeblock]]
201
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vtkWriter.write(1.0);
202
203 1 if (GetPropType<TypeTag, Properties::ModelTraits>::enableEnergyBalance())
204 {
205 std::cout << "mass / energy flux at middle is: " << flux.netFlux("middle") << std::endl;
206 std::cout << "mass / energy flux at outlet is: " << flux.netFlux("outlet") << std::endl;
207 }
208 else
209 {
210
6/14
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 3 std::cout << "mass flux at middle is: " << flux.netFlux("middle") << std::endl;
211
6/14
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 3 std::cout << "mass flux at outlet is: " << flux.netFlux("outlet") << std::endl;
212 }
213
214
6/14
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 3 std::cout << "volume flux at middle is: " << flux.netFlux("middle")[0] << std::endl;
215
7/16
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✓ Branch 11 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✓ Branch 17 taken 1 times.
✓ Branch 18 taken 1 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 3 std::cout << "volume flux at outlet is: " << flux.netFlux("outlet")[0] << std::endl;
216
217
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 if (mpiHelper.rank() == 0)
218
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Parameters::print();
219
220 1 return 0;
221 } // end main
222 // [[/codeblock]]
223 // #### Exception handling
224 // In this part of the main file we catch and print possible exceptions that could
225 // occur during the simulation.
226 // [[details]] exception handler
227 // [[codeblock]]
228 // errors related to run-time parameters
229 ✗ catch (Dumux::ParameterException &e)
230 {
231 ✗ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
232 ✗ return 1;
233 }
234 // errors related to the parsing of Dune grid files
235 ✗ catch (Dune::DGFException & e)
236 {
237 ✗ std::cerr << "DGF exception thrown (" << e <<
238 "). Most likely, the DGF file name is wrong "
239 "or the DGF file is corrupted, "
240 "e.g. missing hash at end of file or wrong number (dimensions) of entries."
241 ✗ << " ---> Abort!" << std::endl;
242 ✗ return 2;
243 }
244 // generic error handling with Dune::Exception
245 ✗ catch (Dune::Exception &e)
246 {
247 ✗ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
248 ✗ return 3;
249 }
250 // other exceptions
251 ✗ catch (...)
252 {
253 ✗ std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
254 ✗ return 4;
255 }
256 // [[/codeblock]]
257 // [[/details]]
258 // [[/content]]
259