GCC Code Coverage Report


Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/examples/2pinfiltration/main.cc
Date: 2024-05-04 19:09:25
Exec Total Coverage
Lines: 77 86 89.5%
Functions: 3 3 100.0%
Branches: 127 291 43.6%

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 file `main.cc`
8 // [[content]]
9 //
10 // This is the main file for the 2pinfiltration example. Here we can see the programme sequence and how the system is solved using Newton's method
11 // ### Included header files
12 // [[details]] includes
13 // [[codeblock]]
14 #include <config.h>
15
16 // Standard header file for C++, for in- and output.
17 #include <iostream>
18 // [[/codeblock]]
19
20 // A Dune helper class for enabling parallel simulations with MPI
21 #include <dune/common/parallel/mpihelper.hh>
22
23 // 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. Additionally, we include the parameter class, which manages the definition of input parameters by a default value, the inputfile or the command line.
24 #include <dumux/common/properties.hh>
25 #include <dumux/common/parameters.hh>
26 #include <dumux/common/initialize.hh>
27
28 //We include the linear solver to be used to solve the linear system and the nonlinear Newton's method
29 #include <dumux/linear/istlsolvers.hh>
30 #include <dumux/linear/linearsolvertraits.hh>
31 #include <dumux/linear/linearalgebratraits.hh>
32 #include <dumux/nonlinear/newtonsolver.hh>
33
34 // Further, we include assembler, which assembles the linear systems for finite volume schemes (box-scheme, tpfa-approximation, mpfa-approximation) and a file that defines the different differentiation methods used to compute the derivatives of the residual
35 #include <dumux/assembly/fvassembler.hh>
36 #include <dumux/assembly/diffmethod.hh>
37
38 // We need the following class to simplify the writing of dumux simulation data to VTK format.
39 #include <dumux/io/vtkoutputmodule.hh>
40 // The gridmanager constructs a grid from the information in the input or grid file.
41 #include <dumux/io/grid/gridmanager_alu.hh>
42
43 //We include several files which are needed for the adaptive grid
44 #include <dumux/adaptive/adapt.hh>
45 #include <dumux/adaptive/markelements.hh>
46 #include <dumux/adaptive/initializationindicator.hh>
47 #include <dumux/porousmediumflow/2p/griddatatransfer.hh>
48 #include <dumux/porousmediumflow/2p/gridadaptindicator.hh>
49
50 // Finally, we include the properties which configure the simulation
51 #include "properties.hh"
52 // [[/details]]
53 //
54
55 // ### The main function
56 // At the beginning of the simulation, we create an alias for our problem type tag, for which we have
57 // defined the properties of our simulation (see `properties.hh`). Additionally, we have to create an
58 // instance of `Dune::MPIHelper` and parse the run-time arguments:
59 // [[codeblock]]
60 1 int main(int argc, char** argv) try
61 {
62 1 using namespace Dumux;
63
64 // we define the type tag for this problem
65 1 using TypeTag = Properties::TTag::PointSourceExample;
66
67 // (maybe) initialize MPI and/or multithreading backend
68
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 Dumux::initialize(argc, argv);
69
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 const auto& mpiHelper = Dune::MPIHelper::instance();
70
71 // Construct parameter tree from command line arguments (and input file if given in the arguments)
72
9/27
✓ 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.
4 Parameters::init(argc, argv);
73 // [[/codeblock]]
74
75 // #### Create the grid
76 // The `gridManager` creates a grid from our parameter choices in the input file (which could be a grid
77 // file or specified domain dimensions and number of cells, as done in this example).
78 // [[codeblock]]
79
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
2 GridManager<GetPropType<TypeTag, Properties::Grid>> gridManager;
80
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();
81
82 // The instationary non-linear problem is run on this grid.
83 //
84 // we compute on the leaf grid view
85
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();
86 // [[/codeblock]]
87
88 // #### Set-up of the problem
89 // We build the finite volume geometry, which allows us to iterate over subcontrolvolumes (scv) and
90 // subcontrolvolume faces (scvf) embedded in the elements of the grid partition.
91 1 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
92
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);
93
94 // In the problem, we define the boundary and initial conditions and compute the point sources.
95 // The `computePointSourceMap` function is inherited from `FVProblem` (see `dumux/common/fvproblem.hh`).
96 // It calls the `addPointSources` method, which is empty per default, but can be overloaded in user problems.
97 // As we have seen, in this example we do specify a point source (see `problem.hh`).
98 // [[codeblock]]
99 1 using Problem = GetPropType<TypeTag, Properties::Problem>;
100
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);
101
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
2 problem->computePointSourceMap();
102 // [[/codeblock]]
103
104 // We initialize the solution vector (all primary variables on the grid) and then use it to initialize the
105 // `gridVariables`, which provide access to both primary & secondary variables (per sub-control volume).
106 1 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
107
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 x;
108
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
2 problem->applyInitialSolution(x);
109
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;
110
111 1 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
112
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);
113
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
2 gridVariables->init(x);
114
115 // ##### Grid adaption
116 // The following piece of code shows how we instantiate an indicator class which determines where to adapt the grid,
117 // and a data transfer class that maps a solution to an adapted grid.
118 // [[codeblock]]
119 1 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
120
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 const Scalar refineTol = getParam<Scalar>("Adaptive.RefineTolerance");
121
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 const Scalar coarsenTol = getParam<Scalar>("Adaptive.CoarsenTolerance");
122 // We use an indicator for a two-phase flow problem that is saturation-dependent.
123 // It is defined in the file `dumux/porousmediumflow/2p/gridadaptindicator.hh.`
124 // and allows to set the minimum and maximum allowed refinement levels via the input parameters.
125
5/12
✓ 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.
✓ Branch 14 taken 1 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
4 TwoPGridAdaptIndicator<TypeTag> indicator(gridGeometry);
126 // The data transfer performs the transfer of data on a grid from before to after adaptation
127 // and is defined in the file `dumux/porousmediumflow/2p/griddatatransfer.hh`.
128 // Its main functions are to store and reconstruct the primary variables.
129
3/10
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 1 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 10 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
4 TwoPGridDataTransfer<TypeTag> dataTransfer(problem, gridGeometry, gridVariables, x);
130 // [[/codeblock]]
131
132 // Convenience functions to perform a grid adaptation given a pre-calculated indicator.
133 // [[codeblock]]
134 17 const auto adaptGridAndVariables = [&] (auto& indicator) {
135 // We mark and adapt the elements according to the indicator.
136 16 bool wasAdapted = false;
137
2/2
✓ Branch 2 taken 7 times.
✓ Branch 3 taken 1 times.
16 if (markElements(gridManager.grid(), indicator))
138 14 wasAdapted = adapt(gridManager.grid(), dataTransfer);
139
140 // In case of any adapted elements, the gridvariables and the pointsourcemap are updated.
141
1/2
✓ Branch 0 taken 7 times.
✗ Branch 1 not taken.
14 if (wasAdapted)
142 {
143 // We overwrite the old solution with the new (resized & interpolated) one
144 14 xOld = x;
145 // We initialize the secondary variables to the new (and "new old") solution
146 28 gridVariables->updateAfterGridAdaption(x);
147 // we update the point source map after adaption
148 28 problem->computePointSourceMap();
149 }
150 17 };
151 // [[/codeblock]]
152
153 // We do initial refinements around sources/BCs, for which we use the `GridAdaptInitializationIndicator` defined in `dumux/adaptive/initializationindicator.hh`.
154 // Afterwards, depending on the initial conditions, another grid adaptation using our `indicator` above might be necessary, which uses
155 // `Adaptive.RefineTolerance` and `Adaptive.CoarsenTolerance` for this step.
156
3/10
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 1 times.
✗ Branch 9 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
4 GridAdaptInitializationIndicator<TypeTag> initIndicator(problem, gridGeometry, gridVariables);
157
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 const auto maxLevel = getParam<std::size_t>("Adaptive.MaxLevel", 0);
158
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 times.
3 for (std::size_t i = 0; i < maxLevel; ++i)
159 {
160
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
2 initIndicator.calculate(x);
161
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
2 adaptGridAndVariables(initIndicator);
162 }
163
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 indicator.calculate(x, refineTol, coarsenTol);
164
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 adaptGridAndVariables(indicator);
165 // [[/codeblock]]
166
167 // #### Solving the problem
168
169 // We get some time loop parameters from the input file params.input
170 1 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
171
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
172
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
173
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 const auto dt = getParam<Scalar>("TimeLoop.DtInitial");
174
175 // and initialize the vtkoutput. Each model has a predefined model specific output with relevant parameters for that model.
176 1 using IOFields = GetPropType<TypeTag, Properties::IOFields>;
177
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());
178 1 using VelocityOutput = GetPropType<TypeTag, Properties::VelocityOutput>;
179
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));
180
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 IOFields::initOutputModule(vtkWriter); // Add model specific output fields
181
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 vtkWriter.write(0.0);
182
183 // We instantiate the time loop
184
1/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
2 auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
185
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
2 timeLoop->setMaxTimeStepSize(maxDt);
186
187 // and set the assembler with the time loop because we have an instationary problem
188 1 using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
189
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);
190
191 // We set the linear solver and the non-linear solver
192 1 using LinearSolver = AMGBiCGSTABIstlSolver<LinearSolverTraits<GridGeometry>,
193 LinearAlgebraTraitsFromAssembler<Assembler>>;
194
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());
195
196 1 using NewtonSolver = Dumux::NewtonSolver<Assembler, LinearSolver>;
197
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 18 taken 1 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 1 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 24 not taken.
✗ Branch 25 not taken.
5 NewtonSolver nonLinearSolver(assembler, linearSolver);
198
199 // ##### The time loop
200 // We start the time loop. In each time step before we start calculating a new solution we check if we have to refine the mesh again based on the solution of the previous time step. If the grid is adapted we update the gridvariables and the pointsourcemap. Afterwards, the solution of that time step is calculated.
201 // [[codeblock]]
202
2/4
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
2 timeLoop->start(); do
203 {
204 // We only want to refine/coarsen after first time step is finished, not before.
205 // The initial refinement was already done before the start of the time loop.
206 // This means we only refine when the time is greater than 0. Note that we now also
207 // have to update the assembler, since the sizes of the residual vector and jacobian matrix change.
208
6/6
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 1 times.
18 if (timeLoop->time() > 0)
209 {
210
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
5 indicator.calculate(x, refineTol, coarsenTol);
211
1/2
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
5 adaptGridAndVariables(indicator);
212
2/4
✓ Branch 1 taken 5 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 5 times.
✗ Branch 5 not taken.
10 assembler->updateAfterGridAdaption();
213 }
214
215 // We solve the non-linear system with time step control.
216
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
12 nonLinearSolver.solve(x, *timeLoop);
217
218 // We make the new solution the old solution.
219
1/2
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
6 xOld = x;
220
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
12 gridVariables->advanceTimeStep();
221
222 // We advance to the time loop to the next step.
223
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
12 timeLoop->advanceTimeStep();
224
225 // We write vtk output for each time step
226
3/6
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 6 times.
✗ Branch 8 not taken.
18 vtkWriter.write(timeLoop->time());
227
228 // We report statistics of this time step
229
2/4
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
12 timeLoop->reportTimeStep();
230
231 // We set a new dt as suggested by the newton solver for the next time step
232
5/10
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 6 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 6 times.
✓ Branch 9 taken 6 times.
✗ Branch 10 not taken.
30 timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
233
234
4/6
✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 6 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 5 times.
✓ Branch 7 taken 1 times.
12 } while (!timeLoop->finished());
235 // [[/codeblock]]
236
237 // The following piece of code prints a final status report of the time loop
238 // before the program is terminated and we print he dumux end message
239 // [[codeblock]]
240
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
2 timeLoop->finalize(leafGridView.comm());
241
242 // print dumux end message
243
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)
244
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
1 Parameters::print();
245
246 1 return 0;
247 } // end main
248 //[[/codeblock]]
249 // ### Exception handling
250 // In this part of the main file we catch and print possible exceptions that could
251 // occur during the simulation.
252 // [[details]] exception handler
253 // [[codeblock]]
254 // errors related to run-time parameters
255 catch (const Dumux::ParameterException &e)
256 {
257 std::cerr << std::endl << e << " ---> Abort!" << std::endl;
258 return 1;
259 }
260 catch (const Dune::Exception &e)
261 {
262 std::cerr << "Dune reported error: " << e << " ---> Abort!" << std::endl;
263 return 3;
264 }
265 catch (...)
266 {
267 std::cerr << "Unknown exception thrown! ---> Abort!" << std::endl;
268 return 4;
269 }
270 // [[/codeblock]]
271 // [[/details]]
272 // [[/content]]
273