GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/examples/1protationsymmetry/main.cc
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 66 76 86.8%
Functions: 2 3 66.7%
Branches: 116 280 41.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 // ## The main program (`main.cc`)
8 // This file contains the main program flow. In this example, we solve a stationary
9 // and rotationally symmetric single-phase problem for a sequence of refined grids
10 // and compute the convergence rates.
11 // [[content]]
12 // ### Includes
13 // [[details]] includes
14 // [[codeblock]]
15 #include <config.h>
16
17 #include <iostream>
18 #include <vector>
19
20 #include <dumux/common/initialize.hh>
21 #include <dumux/common/properties.hh> // for GetPropType
22 #include <dumux/common/parameters.hh> // for getParam
23 #include <dumux/common/integrate.hh> // for integrateL2Error
24
25 #include <dumux/linear/istlsolvers.hh>
26 #include <dumux/linear/linearsolvertraits.hh>
27 #include <dumux/linear/linearalgebratraits.hh>
28 #include <dumux/linear/pdesolver.hh>
29 #include <dumux/assembly/fvassembler.hh>
30 #include <dumux/assembly/diffmethod.hh>
31
32 #include <dumux/io/vtkoutputmodule.hh>
33 #include <dumux/io/grid/gridmanager_yasp.hh>
34
35 #include "properties.hh"
36 // [[/codeblock]]
37 // [[/details]]
38 //
39 // ### Beginning of the main function
40 // [[codeblock]]
41 1 int main(int argc, char** argv) try
42 {
43 1 using namespace Dumux;
44
45 // maybe initialize MPI and/or multithreading backend
46
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 Dumux::initialize(argc, argv);
47
48 // We parse the command line arguments.
49
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);
50
51 // Convenience alias for the type tag of the problem.
52 1 using TypeTag = Properties::TTag::OnePRotSym;
53 // [[/codeblock]]
54
55 // ### Create the grid and the grid geometry
56 // [[codeblock]]
57 // The grid manager can be used to create a grid from the input file
58 1 using Grid = GetPropType<TypeTag, Properties::Grid>;
59
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 2 GridManager<Grid> gridManager;
60
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();
61
62 // We compute on the leaf grid view.
63
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();
64
65 // instantiate the grid geometry
66 1 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
67
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);
68 // [[/codeblock]]
69
70 // ### Initialize the problem and grid variables
71 // [[codeblock]]
72 1 using Problem = GetPropType<TypeTag, Properties::Problem>;
73
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);
74
75 // We define a function to update the discrete analytical solution vector
76 // using the exactSolution() function in the problem
77 11 const auto updateAnalyticalSolution = [&](auto& pExact, auto& vExact)
78 {
79 30 pExact.resize(gridGeometry->numDofs());
80 40 vExact.resize(gridGeometry->elementMapper().size());
81
82
2/12
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3105 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✓ Branch 6 taken 3100 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
 12430 for (const auto& element : elements(gridGeometry->gridView()))
83 {
84
2/4
✓ Branch 1 taken 3100 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3100 times.
✗ Branch 5 not taken.
 12400 auto fvGeometry = localView(*gridGeometry);
85
1/2
✓ Branch 1 taken 3100 times.
✗ Branch 2 not taken.
 6200 fvGeometry.bindElement(element);
86
4/4
✓ Branch 0 taken 6200 times.
✓ Branch 1 taken 3100 times.
✓ Branch 2 taken 6200 times.
✓ Branch 3 taken 3100 times.
 37200 for (auto&& scv : scvs(fvGeometry))
87 {
88 37200 pExact[scv.dofIndex()] = problem->exactSolution(scv.dofPosition());
89 }
90 18600 const auto eIdx = gridGeometry->elementMapper().index(element);
91 12400 vExact[eIdx] = problem->exactVelocity(element.geometry().center());
92 }
93 11 };
94
95 // instantiate and initialize the discrete and exact solution vectors
96 1 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
97
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 SolutionVector p(gridGeometry->numDofs());
98
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 SolutionVector pExact;
99
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 std::vector<double> vExact;
100
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 updateAnalyticalSolution(pExact, vExact);
101
102 // instantiate and initialize the grid variables
103 1 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
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 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(p);
106 // [[/codeblock]]
107
108 // ### Initialize VTK output
109
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 VtkOutputModule<GridVariables, SolutionVector> vtkWriter(*gridVariables, p, problem->name());
110
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 GetPropType<TypeTag, Properties::IOFields>::initOutputModule(vtkWriter);
111
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 vtkWriter.addField(pExact, "pExact"); // add the exact solution to the output fields
112
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 vtkWriter.addField(vExact, "vExact"); // add the exact velocity to the output fields
113
114 1 using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
115
4/8
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 1 times.
 3 vtkWriter.addVelocityOutput(std::make_shared<VelocityOutput>(*gridVariables));
116
117 // ### Instantiate the solver
118 // We use the `LinearPDESolver` class, which is instantiated on the basis
119 // of an assembler and a linear solver. When the `solve` function of the
120 // `LinearPDESolver` is called, it uses the assembler and linear
121 // solver classes to assemble and solve the linear system around the provided
122 // solution and stores the result therein.
123 // [[codeblock]]
124 1 using Assembler = FVAssembler<TypeTag, DiffMethod::analytic>;
125
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);
126
127 1 using LinearSolver = ILUBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>,
128 LinearAlgebraTraitsFromAssembler<Assembler>>;
129
6/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 15 taken 1 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 18 not taken.
 6 auto linearSolver = std::make_shared<LinearSolver>(gridGeometry->gridView(), gridGeometry->dofMapper());
130
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 21 taken 1 times.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
 5 LinearPDESolver<Assembler, LinearSolver> solver(assembler, linearSolver);
131
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 solver.setVerbosity(0); // suppress output during solve()
132 // [[/codeblock]]
133
134 // ### Solution of the problem and error computation
135 // The problem is solved by calling `solve` on the instance of `LinearPDESolver`
136 // that we have created above. In the following piece of code, we solve the
137 // problem on the initial refinement and compute the corresponding L2 error.
138 // For a convenient way of computing the L2 error, the function `integrateL2Error`
139 // can be used.
140 // [[codeblock]]
141
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 solver.solve(p);
142
143 // container to store the L2 errors for the different refinements
144
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 const int numRefinements = getParam<int>("Grid.RefinementSteps");
145
2/6
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
 3 std::vector<double> l2Errors(numRefinements);
146
147 // use third order error integration
148 1 constexpr int orderQuadratureRule = 3;
149
150 // compute initial L2 error
151
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 l2Errors[0] = integrateL2Error(*gridGeometry, p, pExact, orderQuadratureRule);
152 // [[/codeblock]]
153
154 // This procedure is now repeated for the number of refinements as specified
155 // in the input file.
156 // [[codeblock]]
157
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 1 times.
 5 for (int stepIdx = 1; stepIdx < numRefinements; stepIdx++)
158 {
159 // Globally refine the grid once
160
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 4 gridManager.grid().globalRefine(1);
161
162 // update the grid geometry, the grid variables and
163 // the solution vectors now that the grid has been refined
164
4/8
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 4 times.
✗ Branch 11 not taken.
 8 gridGeometry->update(gridManager.grid().leafGridView());
165
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 8 gridVariables->updateAfterGridAdaption(p);
166
167 // this recreates the linear system, i.e. the sizes of
168 // the right hand side vector and the Jacobian matrix,
169 // and its sparsity pattern.
170
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 8 assembler->updateAfterGridAdaption();
171
172
3/6
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 4 times.
✗ Branch 8 not taken.
 12 p.resize(gridGeometry->numDofs());
173
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 updateAnalyticalSolution(pExact, vExact);
174
175 // solve problem on refined grid
176
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 solver.solve(p);
177 // [[/codeblock]]
178 // #### Post-processing and output
179 // At the end of each refinement step, the convergence
180 // rate is printed to the terminal.
181 // [[codeblock]]
182 // Calculate the L2 error using the numerical solution
183
2/4
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 4 times.
✗ Branch 5 not taken.
 8 l2Errors[stepIdx] = integrateL2Error(*gridGeometry, p, pExact, orderQuadratureRule);
184
185 // Print the error and convergence rate
186 4 const auto rate = std::log(l2Errors[stepIdx]/l2Errors[stepIdx-1])/std::log(0.5);
187 8 const auto numDofs = gridGeometry->numDofs();
188 12 std::cout << std::setprecision(8) << std::scientific
189
5/10
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 4 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 4 times.
✗ Branch 9 not taken.
✓ Branch 12 taken 4 times.
✗ Branch 13 not taken.
✓ Branch 15 taken 4 times.
✗ Branch 16 not taken.
 24 << "-- L2 error for " << std::setw(5) << numDofs << " dofs: " << l2Errors[stepIdx]
190
1/2
✓ Branch 2 taken 4 times.
✗ Branch 3 not taken.
 8 << ", rate: " << rate
191
1/2
✓ Branch 1 taken 4 times.
✗ Branch 2 not taken.
 4 << std::endl;
192 }
193 // [[/codeblock]]
194
195 // After the last refinement, we write the solution to VTK file format on the
196 // finest grid and exit the main function.
197 // [[codeblock]]
198
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 vtkWriter.write(0.0);
199
200 // program end, return with 0 exit code (success)
201
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 return 0;
202 }
203 // [[/codeblock]]
204 // ### Exception handling
205 // In this part of the main file we catch and print possible exceptions that could
206 // occur during the simulation.
207 // [[details]] error handler
208 ✗ catch (const Dumux::ParameterException &e)
209 {
210 ✗ std::cerr << std::endl << e << " ---> Abort!" << std::endl;
211 ✗ return 1;
212 }
213 ✗ catch (const Dune::DGFException & e)
214 {
215 ✗ std::cerr << "DGF exception thrown (" << e <<
216 "). Most likely, the DGF file name is wrong "
217 "or the DGF file is corrupted, "
218 "e.g. missing hash at end of file or wrong number (dimensions) of entries."
219 ✗ << " ---> Abort!" << std::endl;
220 ✗ return 2;
221 }
222 ✗ catch (const Dune::Exception &e)
223 {
224 ✗ std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
225 ✗ return 3;
226 }
227 // [[/details]]
228 // [[/content]]
229