GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: dumux/examples/liddrivencavity/main.cc
Date: 2025-04-12 19:19:20
Exec Total Coverage
Lines: 92 92 100.0%
Functions: 2 2 100.0%
Branches: 92 180 51.1%
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-FileCopyrightText: 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 #include <dune/common/timer.hh>
18 #include <dumux/common/dumuxmessage.hh>
19 // [[/exclude]]
20
21 // These is DUNE helper class related to parallel computation
22 #include <dune/common/parallel/mpihelper.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/parameters.hh>
27 #include <dumux/common/properties.hh>
28 #include <dumux/common/initialize.hh>
29
30 // The following files contain the multi-domain Newton solver, the available linear solver backends and the assembler for the linear
31 // systems arising from the staggered-grid discretization.
32 #include <dumux/linear/istlsolvers.hh>
33 #include <dumux/linear/linearalgebratraits.hh>
34 #include <dumux/linear/linearsolvertraits.hh>
35 #include <dumux/multidomain/fvassembler.hh>
36 #include <dumux/multidomain/traits.hh>
37 #include <dumux/multidomain/newtonsolver.hh>
38
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 // This class contains functionality for VTK output for models using the staggered finite volume scheme.
45 #include <dumux/io/vtkoutputmodule.hh>
46 #include <dumux/freeflow/navierstokes/velocityoutput.hh>
47
48 // We include the problem header used for this simulation.
49 #include "properties.hh"
50 // [[/details]]
51
52 // The following function writes the velocities and coordinates at x = 0.5 and y = 0.5 into a log file.
53 // [[codeblock]]
54 template<class Problem, class SolutionVector>
55 2 void writeSteadyVelocityAndCoordinates(const Problem& problem, const SolutionVector& sol)
56 {
57 2 const auto& gridGeometry = problem.gridGeometry();
58
3/6
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 2 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
 4 std::ofstream logFilevx(problem.name() + "_vx.log"), logFilevy(problem.name() + "_vy.log");
59
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 logFilevx << "y vx\n";
60
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 logFilevy << "x vy\n";
61
62 static constexpr double eps_ = 1.0e-7;
63
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 2 std::vector<int> dofHandled(gridGeometry.numDofs(), false);
64
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 8194 times.
✗ Branch 5 not taken.
 24578 for (const auto& element : elements(gridGeometry.gridView()))
65 {
66 8192 auto fvGeometry = localView(gridGeometry);
67 8192 fvGeometry.bind(element);
68
4/4
✓ Branch 0 taken 16128 times.
✓ Branch 1 taken 16640 times.
✓ Branch 2 taken 32768 times.
✓ Branch 3 taken 8192 times.
 40960 for (const auto& scv : scvs(fvGeometry))
69 {
70
2/2
✓ Branch 0 taken 16128 times.
✓ Branch 1 taken 16640 times.
 32768 if (dofHandled[scv.dofIndex()])
71 16128 continue;
72
73
2/2
✓ Branch 0 taken 16128 times.
✓ Branch 1 taken 512 times.
 16640 if (!scv.boundary())
74 {
75 16128 const auto& globalPos = scv.dofPosition();
76
2/2
✓ Branch 0 taken 128 times.
✓ Branch 1 taken 16000 times.
 16128 const auto velocity = sol[scv.dofIndex()][0];
77
2/2
✓ Branch 0 taken 128 times.
✓ Branch 1 taken 16000 times.
 16128 dofHandled[scv.dofIndex()] = true;
78
79
2/2
✓ Branch 0 taken 128 times.
✓ Branch 1 taken 16000 times.
 16128 if (std::abs(globalPos[0]-0.5) < eps_)
80
4/8
✓ Branch 1 taken 128 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 128 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 128 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 128 times.
✗ Branch 11 not taken.
 128 logFilevx << globalPos[1] << " " << velocity << "\n";
81
2/2
✓ Branch 0 taken 128 times.
✓ Branch 1 taken 15872 times.
 16000 else if (std::abs(globalPos[1]-0.5) < eps_)
82
4/8
✓ Branch 1 taken 128 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 128 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 128 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 128 times.
✗ Branch 11 not taken.
 128 logFilevy << globalPos[0] << " " << velocity << "\n";
83 }
84 }
85 }
86 2 }
87 // [[/codeblock]]
88
89 // ### The main function
90 // We will now discuss the main program flow implemented within the `main` function.
91 // At the beginning of each program using Dune, an instance of `Dune::MPIHelper` has to
92 // be created. Moreover, we parse the run-time arguments from the command line and the
93 // input file:
94 // [[codeblock]]
95 2 int main(int argc, char** argv)
96 {
97 2 using namespace Dumux;
98
99 // maybe initialize MPI and/or multithreading backend
100 2 Dumux::initialize(argc, argv);
101 2 const auto& mpiHelper = Dune::MPIHelper::instance();
102
103 // parse command line arguments and input file
104
3/8
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 4 Parameters::init(argc, argv);
105 // [[/codeblock]]
106
107 // We define a convenience alias for the type tags of the problems. The type
108 // tag contains all the properties that are needed to define the model and the problem
109 // setup. Throughout the main file, we will obtain types defined each type tag
110 // using the property system, i.e. with `GetPropType`.
111 2 using MomentumTypeTag = Properties::TTag::LidDrivenCavityExampleMomentum;
112 2 using MassTypeTag = Properties::TTag::LidDrivenCavityExampleMass;
113
114 // #### Step 1: Create the grid
115 // The `GridManager` class creates the grid from information given in the input file.
116 // This can either be a grid file, or in the case of structured grids, one can specify the coordinates
117 // of the corners of the grid and the number of cells to be used to discretize each spatial direction.
118 // [[codeblock]]
119 2 GridManager<GetPropType<MomentumTypeTag, Properties::Grid>> gridManager;
120
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 gridManager.init();
121
122 // We compute on the leaf grid view.
123
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 2 const auto& leafGridView = gridManager.grid().leafGridView();
124 // [[/codeblock]]
125
126 // #### Step 2: Setting up and solving the problem
127 // First, a finite volume grid geometry is constructed from the grid that was created above.
128 // This builds the sub-control volumes (scv) and sub-control volume faces (scvf) for each element
129 // of the grid partition.
130 // This is done for both the momentum and mass grid geometries
131 2 using MomentumGridGeometry = GetPropType<MomentumTypeTag, Properties::GridGeometry>;
132
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto momentumGridGeometry = std::make_shared<MomentumGridGeometry>(leafGridView);
133 2 using MassGridGeometry = GetPropType<MassTypeTag, Properties::GridGeometry>;
134
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto massGridGeometry = std::make_shared<MassGridGeometry>(leafGridView);
135
136 // We introduce the multidomain coupling manager, which will couple the mass and the momentum problems
137 // We can obtain the type from either the `MomentumTypeTag` or the `MassTypeTag` because they are mutually coupled with the same manager
138 2 using CouplingManager = GetPropType<MomentumTypeTag, Properties::CouplingManager>;
139
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto couplingManager = std::make_shared<CouplingManager>();
140
141 // We now instantiate the problems, in which we define the boundary and initial conditions.
142 2 using MomentumProblem = GetPropType<MomentumTypeTag, Properties::Problem>;
143
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto momentumProblem = std::make_shared<MomentumProblem>(momentumGridGeometry, couplingManager);
144 2 using MassProblem = GetPropType<MassTypeTag, Properties::Problem>;
145
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto massProblem = std::make_shared<MassProblem>(massGridGeometry, couplingManager);
146
147 // We set a solution vector which consist of two parts: one part (indexed by `massIdx`)
148 // is for the pressure degrees of freedom (`dofs`) living in grid cell centers. Another part
149 // (indexed by `momentumIdx`) is for degrees of freedom defining the normal velocities on grid cell faces.
150 // We initialize the solution vector by what was defined as the initial solution of the the problem.
151 2 constexpr auto momentumIdx = CouplingManager::freeFlowMomentumIndex;
152 2 constexpr auto massIdx = CouplingManager::freeFlowMassIndex;
153 2 using Traits = MultiDomainTraits<MomentumTypeTag, MassTypeTag>;
154 2 using SolutionVector = typename Traits::SolutionVector;
155
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 SolutionVector x;
156
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 momentumProblem->applyInitialSolution(x[momentumIdx]);
157
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 massProblem->applyInitialSolution(x[massIdx]);
158
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto xOld = x;
159
160 // We use the initial solution vector to create the `gridVariables`.
161 // The grid variables are used store variables (primary and secondary variables) on sub-control volumes and faces (volume and flux variables).
162 2 using MomentumGridVariables = GetPropType<MomentumTypeTag, Properties::GridVariables>;
163
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto momentumGridVariables = std::make_shared<MomentumGridVariables>(momentumProblem, momentumGridGeometry);
164 2 using MassGridVariables = GetPropType<MassTypeTag, Properties::GridVariables>;
165
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto massGridVariables = std::make_shared<MassGridVariables>(massProblem, massGridGeometry);
166
167 // using the problems and the grid variables, the coupling manager and the grid variables are initialized with the initial solution.
168 // The grid variables have to be initialized _after_ the coupling manager.
169 // This is because they require the correct coupling context between the mass and momentum model to be initialized.
170
2/6
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
 8 couplingManager->init(momentumProblem, massProblem, std::make_tuple(momentumGridVariables, massGridVariables), x, xOld);
171
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 momentumGridVariables->init(x[momentumIdx]);
172
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 massGridVariables->init(x[massIdx]);
173
174 // We get some time loop parameters from the input file
175 // and instantiate the time loop
176 2 using Scalar = typename Traits::Scalar;
177
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 const auto tEnd = getParam<Scalar>("TimeLoop.TEnd");
178
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 const auto maxDt = getParam<Scalar>("TimeLoop.MaxTimeStepSize");
179
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 const auto dt = getParam<Scalar>("TimeLoop.DtInitial");
180
181
1/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 2 auto timeLoop = std::make_shared<TimeLoop<Scalar>>(0, dt, tEnd);
182
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 timeLoop->setMaxTimeStepSize(maxDt);
183
184 // We then initialize the predefined model-specific output vtk output.
185 2 using IOFields = GetPropType<MassTypeTag, Properties::IOFields>;
186
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
 2 VtkOutputModule vtkWriter(*massGridVariables, x[massIdx], massProblem->name());
187
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 IOFields::initOutputModule(vtkWriter); // Add model specific output fields
188
3/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
 4 vtkWriter.addVelocityOutput(std::make_shared<NavierStokesVelocityOutput<MassGridVariables>>());
189
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 vtkWriter.write(0.0);
190
191 // To solve the non-linear problem at hand, we use the `NewtonSolver`,
192 // which we have to tell how to assemble and solve the system in each
193 // iteration. Here, we use the direct linear solver UMFPack.
194 2 using Assembler = MultiDomainFVAssembler<Traits, CouplingManager, DiffMethod::numeric>;
195
1/2
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 6 auto assembler = std::make_shared<Assembler>(std::make_tuple(momentumProblem, massProblem),
196 4 std::make_tuple(momentumGridGeometry, massGridGeometry),
197 2 std::make_tuple(momentumGridVariables, massGridVariables),
198 2 couplingManager, timeLoop, xOld);
199 // the linear solver
200 2 using LinearSolver = Dumux::UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
201
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 auto linearSolver = std::make_shared<LinearSolver>();
202
203 // the non-linear solver
204 2 using NewtonSolver = MultiDomainNewtonSolver<Assembler, LinearSolver, CouplingManager>;
205
4/14
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 7 taken 2 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 2 times.
✗ Branch 12 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
 10 NewtonSolver nonLinearSolver(assembler, linearSolver, couplingManager);
206
207 // ##### The time loop
208 // In each time step, we solve the non-linear system of equations, write
209 // the current solution into VTK files and prepare for the next time step.
210 // [[codeblock]]
211
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 timeLoop->start(); do
212 {
213 // We solve the non-linear system with time step control.
214
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 nonLinearSolver.solve(x, *timeLoop);
215
216 // We make the new solution the old solution.
217
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 xOld = x;
218
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 massGridVariables->advanceTimeStep();
219
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 momentumGridVariables->advanceTimeStep();
220
221 // We advance to the time loop to the next step.
222
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 timeLoop->advanceTimeStep();
223
224 // We write vtk output for each time step.
225
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 vtkWriter.write(timeLoop->time());
226
227 // We report statistics of this time step.
228
1/2
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
 10 timeLoop->reportTimeStep();
229
230 // We set a new dt as suggested by the newton solver for the next time step.
231
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✓ Branch 3 taken 10 times.
✗ Branch 4 not taken.
 20 timeLoop->setTimeStepSize(nonLinearSolver.suggestTimeStepSize(timeLoop->timeStepSize()));
232
233
3/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8 times.
✓ Branch 4 taken 2 times.
 10 } while (!timeLoop->finished());
234 // [[/codeblock]]
235
236 // We write the velocities and coordinates at x = 0.5 and y = 0.5 into a file.
237
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 writeSteadyVelocityAndCoordinates(*momentumProblem, x[momentumIdx]);
238
239 // The following piece of code prints a final status report of the time loop
240 // before the program is terminated.
241 // [[codeblock]]
242
2/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 2 timeLoop->finalize(leafGridView.comm());
243
244 // print used and unused parameters
245
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if (mpiHelper.rank() == 0)
246
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 2 Parameters::print();
247
248 2 return 0;
249 22 } // end main
250 // [[/codeblock]]
251 // [[/content]]
252