GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/examples/porenetwork_upscaling/main.cc
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 62 65 95.4%
Functions: 3 3 100.0%
Branches: 136 322 42.2%
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 use a single-phase
9 // Pore-Network-Model to evaluate the upscaled Darcy permeability of a given network.
10 // [[content]]
11 // ### Includes
12 // [[details]] includes
13 // [[codeblock]]
14 #include <config.h>
15
16 #include <iostream>
17
18 #include <algorithm>
19
20 #include <dune/common/float_cmp.hh> // for floating point comparison
21 #include <dune/common/exceptions.hh>
22
23 #include <dumux/common/properties.hh> // for GetPropType
24 #include <dumux/common/parameters.hh> // for getParam
25 #include <dumux/common/initialize.hh>
26
27 #include <dumux/linear/istlsolvers.hh>
28 #include <dumux/linear/linearsolvertraits.hh>
29 #include <dumux/linear/linearalgebratraits.hh>
30 #include <dumux/linear/pdesolver.hh>
31 #include <dumux/nonlinear/newtonsolver.hh>
32 #include <dumux/assembly/fvassembler.hh>
33
34 #include <dumux/io/vtkoutputmodule.hh>
35 #include <dumux/porenetwork/common/pnmvtkoutputmodule.hh>
36
37 #include <dumux/io/grid/gridmanager_yasp.hh>
38 #include <dumux/io/grid/porenetwork/gridmanager.hh> // for pore-network grid
39 #include <dumux/porenetwork/common/boundaryflux.hh> // for getting the total mass flux leaving the network
40
41 #include "upscalinghelper.hh"
42 #include "properties.hh"
43 // [[/codeblock]]
44 // [[/details]]
45 //
46 // ### The driver function
47 //
48 // The template argument `TypeTag` determines if we run the example assuming
49 // a creeping flow regime or not. Which regime is selected is set with the parameter
50 // `Problem.AssumeCreepingFlow` in the input file.
51 namespace Dumux {
52
53 template<class TypeTag>
54 4 void runExample()
55 {
56 // ### Create the grid and the grid geometry
57 // [[codeblock]]
58 // The grid manager can be used to create a grid from the input file
59 using GridManager = PoreNetwork::GridManager<3>;
60 4 GridManager gridManager;
61
5/12
✓ 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 9 not taken.
✓ Branch 10 taken 2 times.
✓ Branch 12 taken 2 times.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
 8 gridManager.init();
62
63 // We compute on the leaf grid view.
64
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 4 const auto& leafGridView = gridManager.grid().leafGridView();
65
66 // instantiate the grid geometry
67
1/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 8 auto gridData = gridManager.getGridData();
68 using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>;
69
3/8
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 12 auto gridGeometry = std::make_shared<GridGeometry>(leafGridView, *gridData);
70 // [[/codeblock]]
71
72 // ### Initialize the problem and grid variables
73 // [[codeblock]]
74 using Problem = GetPropType<TypeTag, Properties::Problem>;
75
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 8 auto problem = std::make_shared<Problem>(gridGeometry);
76
77 // instantiate and initialize the discrete and exact solution vectors
78 using SolutionVector = GetPropType<TypeTag, Properties::SolutionVector>;
79
4/10
✓ 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 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
 16 SolutionVector x(gridGeometry->numDofs()); // zero-initializes
80
81 // instantiate and initialize the grid variables
82 using GridVariables = GetPropType<TypeTag, Properties::GridVariables>;
83
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 8 auto gridVariables = std::make_shared<GridVariables>(problem, gridGeometry);
84
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 8 gridVariables->init(x);
85 // [[/codeblock]]
86
87 // ### Instantiate the solver
88 // We use the `NewtonSolver` class, which is instantiated on the basis
89 // of an assembler and a linear solver. When the `solve()` function of the
90 // Newton solver instance is called, it uses the assembler and
91 // linear solver to assemble and solve the linear systems in each Newton step.
92 // [[codeblock]]
93 using Assembler = FVAssembler<TypeTag, DiffMethod::numeric>;
94
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 8 auto assembler = std::make_shared<Assembler>(problem, gridGeometry, gridVariables);
95
96 using LinearSolver = Dumux::UMFPackIstlSolver<SeqLinearSolverTraits, LinearAlgebraTraitsFromAssembler<Assembler>>;
97
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 8 auto linearSolver = std::make_shared<LinearSolver>();
98
99 using NewtonSolver = NewtonSolver<Assembler, LinearSolver>;
100
12/28
✓ 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 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✓ Branch 20 taken 2 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 2 times.
✗ Branch 24 not taken.
✓ Branch 25 taken 2 times.
✗ Branch 26 not taken.
✓ Branch 27 taken 2 times.
✓ Branch 29 taken 2 times.
✗ Branch 30 not taken.
✓ Branch 32 taken 2 times.
✗ Branch 33 not taken.
✗ Branch 34 not taken.
✗ Branch 35 not taken.
✗ Branch 36 not taken.
✗ Branch 37 not taken.
 20 NewtonSolver nonLinearSolver(assembler, linearSolver);
101 // [[/codeblock]]
102
103 // ### Initialize VTK output
104 using VtkOutputFields = GetPropType<TypeTag, Properties::IOFields>;
105 using VtkWriter = PoreNetwork::VtkOutputModule<GridVariables, GetPropType<TypeTag, Properties::FluxVariables>, SolutionVector>;
106
7/14
✓ 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 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 16 taken 2 times.
✗ Branch 17 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
 16 VtkWriter vtkWriter(*gridVariables, x, problem->name());
107
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 VtkOutputFields::initOutputModule(vtkWriter);
108 // specify the field type explicitly since it may not be possible
109 // to deduce this from the vector size in a pore network
110
7/16
✓ 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 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✗ Branch 15 not taken.
✓ Branch 16 taken 2 times.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 8 vtkWriter.addField(gridGeometry->poreVolume(), "poreVolume", Vtk::FieldType::vertex);
111
7/16
✓ 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 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 12 vtkWriter.addField(gridGeometry->throatShapeFactor(), "throatShapeFactor", Vtk::FieldType::element);
112
7/16
✓ 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 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✓ Branch 18 taken 2 times.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
 12 vtkWriter.addField(gridGeometry->throatCrossSectionalArea(), "throatCrossSectionalArea", Vtk::FieldType::element);
113
114 // ### Prepare the upscaling procedure.
115 // Set up a helper class to determine the total mass flux leaving the network
116
4/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 taken 2 times.
✗ Branch 11 not taken.
 20 const auto boundaryFlux = PoreNetwork::BoundaryFlux(*gridVariables, assembler->localResidual(), x);
117
118 // ### The actual upscaling procedure
119 // #### Instantiate the upscaling helper
120 // [[codeblock]]
121 using Scalar = GetPropType<TypeTag, Properties::Scalar>;
122 using UpscalingHelper = UpscalingHelper<Scalar>;
123 4 UpscalingHelper upscalingHelper;
124 // [[/codeblock]]
125
126 // Set the side lengths used for applying the pressure gradient and calculating the REV outflow area.
127 // One can either specify these values manually (usually more accurate) or let the UpscalingHelper struct
128 // determine it automatically based on the network's bounding box.
129 // [[codeblock]]
130 2 const auto sideLengths = [&]()
131 {
132 using GlobalPosition = typename GridGeometry::GlobalCoordinate;
133 10 if (hasParam("Problem.SideLength"))
134 ✗ return getParam<GlobalPosition>("Problem.SideLength");
135 else
136 4 return upscalingHelper.getSideLengths(*gridGeometry);
137
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 }();
138
139 // pass the side lengths to the problem
140
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 8 problem->setSideLengths(sideLengths);
141 // [[/codeblock]]
142
143 // Get the maximum and minimum pressure gradient and the population of sample points specified in the input file
144 // [[codeblock]]
145
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 const Scalar minPressureGradient = getParam<Scalar>("Problem.MinimumPressureGradient", 1e1);
146
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 const Scalar maxPressureGradient = getParam<Scalar>("Problem.MaximumPressureGradient", 1e10);
147
148
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 if (!(minPressureGradient < maxPressureGradient))
149 ✗ DUNE_THROW(Dune::InvalidStateException, "Maximum pressure gradient must be greater than minimum pressure gradient");
150
151
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 const int numberOfSamples = getParam<int>("Problem.NumberOfPressureGradients", 1);
152 // [[/codeblock]]
153
154 // Iterate over all directions specified before, apply several pressure gradients, calculate the mass flux
155 // and finally determine the the upscaled properties.
156 // [[codeblock]]
157
5/10
✓ 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 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 1 times.
✗ Branch 13 not taken.
 18 const auto directions = getParam<std::vector<std::size_t>>("Problem.Directions", std::vector<std::size_t>{0, 1, 2});
158
3/8
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 4 upscalingHelper.setDirections(directions);
159
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
 16 for (int dimIdx : directions)
160 {
161 // set the direction in which the pressure gradient will be applied
162 8 problem->setDirection(dimIdx);
163
164
2/2
✓ Branch 0 taken 20 times.
✓ Branch 1 taken 2 times.
 44 for (int i = 0; i < numberOfSamples; i++)
165 {
166 // reset the solution
167 40 x = 0;
168
169 // set the pressure gradient to be applied
170
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 18 times.
 40 Scalar pressureGradient = maxPressureGradient * std::exp(i + 1 - numberOfSamples);
171
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 18 times.
 40 if (i == 0)
172
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 4 pressureGradient = std::min(minPressureGradient, pressureGradient);
173
174
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 80 problem->setPressureGradient(pressureGradient);
175
176 // solve problem
177
1/2
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
 40 nonLinearSolver.solve(x);
178
179 // set the sample points
180
7/16
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 20 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 20 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 20 times.
✗ Branch 14 not taken.
✓ Branch 17 taken 20 times.
✗ Branch 18 not taken.
✓ Branch 20 taken 20 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
 200 const Scalar totalFluidMassFlux = boundaryFlux.getFlux(std::vector<int>{ problem->outletPoreLabel() })[0];
181
2/4
✓ Branch 1 taken 20 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 20 times.
✗ Branch 5 not taken.
 80 upscalingHelper.setDataPoints(*problem, totalFluidMassFlux);
182 }
183
184 // write a vtu file for the given direction for the last sample
185
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 vtkWriter.write(dimIdx);
186 }
187
188 // calculate and report the upscaled properties
189 4 constexpr bool isCreepingFlow = std::is_same_v<TypeTag, Properties::TTag::PNMUpscalingCreepingFlow>;
190
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 8 upscalingHelper.calculateUpscaledProperties(*problem, isCreepingFlow);
191
1/2
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
 4 upscalingHelper.report(isCreepingFlow);
192
193 // compare the Darcy permeability with reference data if provided in input file and report in case of inconsistency
194
5/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✗ Branch 13 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✗ Branch 16 not taken.
 4 static const auto referenceData = getParam<std::vector<Scalar>>("Problem.ReferencePermeability", std::vector<Scalar>{});
195
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
 8 if (!referenceData.empty())
196
3/8
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 8 upscalingHelper.compareWithReference(referenceData);
197
198 // plot the results just for non-creeping flow
199 // creeping flow would just result in a straight line (permeability is independent of the pressure gradient)
200
2/6
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
 4 bool plotConductivity = getParam<bool>("Output.PlotConductivity", true);
201
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 2 if (!isCreepingFlow && plotConductivity)
202 ✗ upscalingHelper.plot();
203 4 };
204
205 } // end namespace Dumux
206 // [[/codeblock]]
207
208 // ### The main function
209 // [[details]] main
210 // [[codeblock]]
211 2 int main(int argc, char** argv)
212 {
213 2 using namespace Dumux;
214
215 // Initialize MPI+X environment
216 2 Dumux::initialize(argc, argv);
217
218 // We parse the command line arguments.
219
9/24
✓ 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 taken 2 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 2 times.
✗ Branch 14 not taken.
✓ Branch 15 taken 2 times.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✓ Branch 18 taken 2 times.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✓ Branch 21 taken 2 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.
 8 Parameters::init(argc, argv);
220
221 // Convenience alias for the type tag of the problem.
222 2 using CreepingFlowTypeTag = Properties::TTag::PNMUpscalingCreepingFlow;
223 2 using NonCreepingFlowTypeTag = Properties::TTag::PNMUpscalingNonCreepingFlow;
224 // // [[/codeblock]]
225
226 // user decides whether creeping flow or non-creeping flow should be run
227
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 1 times.
 2 if (getParam<bool>("Problem.AssumeCreepingFlow", false))
228 1 runExample<CreepingFlowTypeTag>();
229 else
230 1 runExample<NonCreepingFlowTypeTag>();
231
232 // program end, return with 0 exit code (success)
233 2 return 0;
234 }
235 // [[/codeblock]]
236 // [[/details]]
237 // [[/content]]
238