GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/examples/porenetwork_upscaling/upscalinghelper.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 68 117 58.1%
Functions: 8 11 72.7%
Branches: 100 402 24.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 #ifndef DUMUX_PNM_ONEP_PERMEABILITY_UPSCALING_HELPER_HH
9 #define DUMUX_PNM_ONEP_PERMEABILITY_UPSCALING_HELPER_HH
10
11 // ## Upscaling helper struct (`upscalinghelper.hh`)
12 //
13 // This file contains the __upscaling helper struct__ which considers the volume flux leaving
14 // the pore network in flow direction in order to find the upscaled Darcy permeability.
15 // [[content]]
16 #include <iostream>
17 #include <ostream>
18 #include <iomanip>
19 #include <numeric>
20 #include <functional>
21
22 #include <dune/common/float_cmp.hh>
23
24 #include <dumux/io/gnuplotinterface.hh>
25 #include <dumux/io/format.hh>
26 #include <dumux/common/parameters.hh>
27
28 namespace Dumux {
29
30 template<class Scalar>
31 class UpscalingHelper
32 {
33 public:
34 // ### Set data points to calculate intrinsic permeability and Forchheimer coefficient
35 // This function first evaluates the mass flux leaving the network in the direction of the applied pressure gradient.
36 // Afterwards, the mass flux is converted into an volume flux which is used to calculate the apparent velocity.
37 // Then apparent permeability of the network is computed and stored for furthure calculations.
38 // [[codeblock]]
39 template <class Problem>
40 40 void setDataPoints(const Problem &problem, const Scalar totalMassFlux)
41 {
42 // get the domain side lengths from the problem
43 40 auto sideLengths = problem.sideLengths();
44
45
46 // get the applied pressure gradient
47 40 const auto pressureGradient = problem.pressureGradient();
48 40 const auto pressureDrop = pressureGradient * sideLengths[problem.direction()];
49
50 // get the fluid properties
51 40 const auto liquidDensity = problem.liquidDensity();
52 40 const auto liquidDynamicViscosity = problem.liquidDynamicViscosity();
53
54 // convert mass to volume flux
55 40 const auto volumeFlux = totalMassFlux / liquidDensity;
56 ;
57
58 // calculate apparent velocity
59 40 sideLengths[problem.direction()] = 1.0;
60 80 const auto outflowArea = std::accumulate(sideLengths.begin(), sideLengths.end(), 1.0, std::multiplies<Scalar>());
61 40 const auto vApparent = volumeFlux / outflowArea;
62
63 // compute apparent permeability
64 40 const auto KApparent = vApparent / pressureGradient * liquidDynamicViscosity;
65 // calculate rho v / mu, called inertia to viscous ratio in the rest of the code
66 40 const auto inertiaToViscousRatio = liquidDensity * vApparent / liquidDynamicViscosity;
67
68 // store the required data for further calculations
69 80 totalPressureDrop_[problem.direction()].push_back(pressureDrop);
70 80 apparentVelocity_[problem.direction()].push_back(vApparent);
71 80 apparentPermeability_[problem.direction()].push_back(KApparent);
72 80 inertiaToViscousRatio_[problem.direction()].push_back(inertiaToViscousRatio);
73 40 }
74 // [[/codeblock]]
75
76 // ### Calculate intrinsic permeability and Forchheimer coefficient.
77 // This function first calculate intrinsic permeability and Forchheimer coefficient using linear least squares regression method
78 // and reports them. It also plot the apparent permeability of the porous medium versus Forchheimer number/pressure gradient in each
79 // simulation.
80 // [[codeblock]]
81 template <class Problem>
82 4 void calculateUpscaledProperties(const Problem &problem, bool isCreepingFlow)
83 {
84 4 const auto sideLengths = problem.sideLengths();
85 4 const auto liquidDynamicViscosity = problem.liquidDynamicViscosity();
86
87
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
 16 for (const auto dirIdx : directions_)
88 {
89 // determine Darcy permeability as the maximum permeability of the domain
90
7/12
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 1 times.
✓ Branch 11 taken 1 times.
 24 darcyPermeability_[dirIdx] = *max_element(apparentPermeability_[dirIdx].begin(), apparentPermeability_[dirIdx].end());
91
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 4 if (!isCreepingFlow)
92 {
93
6/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 1 times.
 26 for (int i = 0; i < totalPressureDrop_[dirIdx].size(); i++)
94 {
95 // calculate the Darcy pressure drop.
96
8/8
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 8 times.
✓ Branch 5 taken 2 times.
✓ Branch 6 taken 8 times.
✓ Branch 7 taken 2 times.
 80 const Scalar darcyPressureDrop = liquidDynamicViscosity * apparentVelocity_[dirIdx][i] * sideLengths[dirIdx] / darcyPermeability_[dirIdx];
97
98 // calculate the ratio of Dracy to total pressure drop
99
4/4
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 8 times.
✓ Branch 3 taken 2 times.
 40 const Scalar pressureDropRatio = darcyPressureDrop / totalPressureDrop_[dirIdx][i];
100
101 // set sample points for upscaling of Forchheimer parameters.
102 // first, check the permability ratio to see if the flow regime is Forchheimer.
103
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 2 times.
 20 if (pressureDropRatio < 0.99)
104 {
105 64 samplePointsX_[dirIdx].push_back(inertiaToViscousRatio_[dirIdx][i]);
106 64 samplePointsY_[dirIdx].push_back(1 / apparentPermeability_[dirIdx][i]);
107 }
108 }
109 // determine regression line and accordingly the Forchheimer permeability and the Forchheimer coefficient
110 6 const auto [intercept, slope] = linearRegression(samplePointsX_[dirIdx], samplePointsY_[dirIdx]);
111 2 forchheimerPermeability_[dirIdx] = 1.0 / intercept;
112 2 forchheimerCoefficient_[dirIdx] = slope;
113 2 writePlotDataToFile(dirIdx);
114 }
115 }
116 4 }
117 // [[/codeblock]]
118
119 // ### Determine the domain's side lengths
120 //
121 // We determine the domain side length by using the bounding box of the network
122 // [[codeblock]]
123 template<class GridGeometry>
124 ✗ auto getSideLengths(const GridGeometry& gridGeometry)
125 {
126 using GlobalPosition = typename GridGeometry::GlobalCoordinate;
127 ✗ GlobalPosition result(0.0);
128
129 ✗ std::cout << "Automatically determining side lengths of REV based on bounding box of pore network" << std::endl;
130 ✗ for (int dimIdx = 0; dimIdx < GridGeometry::GridView::dimensionworld; ++dimIdx)
131 ✗ result[dimIdx] = gridGeometry.bBoxMax()[dimIdx] - gridGeometry.bBoxMin()[dimIdx];
132
133 ✗ return result;
134 }
135 // [[/codeblock]]
136
137 // ### Plot the data using Gnuplot
138 //
139 // [[codeblock]]
140 void plot()
141 {
142 // plot permeability ratio vs. Forchheimer number
143 ✗ plotPermeabilityratioVsForchheimerNumber_();
144
145 // plot inverse of apparent permability vs. rho v / mu
146 ✗ plotInversePrmeabilityVsInertiaToViscousRatio_();
147 }
148 // [[/codeblock]]
149 //
150 // ### Save the relevant data for plot
151 // [[codeblock]]
152 void writePlotDataToFile(std::size_t dirIdx)
153 {
154 // write permeability ratio vs. Forchheimer number
155 1 writePermeabilityratioVsForchheimerNumber_(dirIdx);
156
157 // write inverse of apparent permability vs. rho v / mu
158 1 writeInversePrmeabilityVsInertiaToViscousRatio_(dirIdx);
159 }
160 // [[/codeblock]]
161 //
162 // ### Report the upscaled data
163 // [[codeblock]]
164 2 void report(bool isCreepingFlow)
165 {
166 // Report the results for each direction
167
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
 10 for (const auto dirIdx : directions_)
168 {
169 2 std::cout << Fmt::format("\n{:#>{}}\n\n", "", 40)
170 4 << Fmt::format("{}-direction:\n", dirName_[dirIdx])
171
7/20
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
✓ Branch 14 taken 2 times.
✗ Branch 15 not taken.
✓ Branch 17 taken 2 times.
✗ Branch 18 not taken.
✓ Branch 19 taken 2 times.
✗ Branch 20 not taken.
✗ Branch 21 not taken.
✓ Branch 22 taken 2 times.
✓ Branch 23 taken 2 times.
✗ 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.
 6 << Fmt::format("-- Darcy (intrinsic) permeability = {:.3e} m^2\n", darcyPermeability_[dirIdx]);
172
173 // Report non-creeping flow upscaled properties
174
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1 times.
 2 if (!isCreepingFlow)
175 {
176
2/6
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 2 std::cout << Fmt::format("-- Forchheimer permeability = {:.3e} m^2\n", forchheimerPermeability_[dirIdx]);
177
2/6
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
 3 std::cout << Fmt::format("-- Forchheimer coefficient = {:.3e} m^-1\n", forchheimerCoefficient_[dirIdx]);
178 }
179
180
3/8
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 2 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 4 std::cout << Fmt::format("\n{:#>{}}\n", "", 40) << std::endl;
181 }
182 2 }
183 // [[/codeblock]]
184 //
185 // ### Compare with reference data provided in input file
186 // [[codeblock]]
187 2 void compareWithReference(std::vector<Scalar> referenceData)
188 {
189
4/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
 8 for (const auto dirIdx : directions_)
190 {
191
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 const auto K = darcyPermeability_[dirIdx];
192
3/6
✓ 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.
 2 static const Scalar eps = getParam<Scalar>("Problem.TestEpsilon", 1e-3);
193
3/6
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 2 times.
 6 if (Dune::FloatCmp::ne<Scalar>(K, referenceData[dirIdx], eps))
194 {
195 ✗ std::cerr << "Calculated permeability of " << K << " in "
196 ✗ <<dirName_[dirIdx]<<"-direction does not match with reference value of "
197 ✗ << referenceData[dirIdx] << std::endl;
198 }
199 }
200 2 }
201 // [[/codeblock]]
202 //
203 // ### Set the directions that need to be considered
204 //
205 // [[codeblock]]
206 void setDirections(std::vector<std::size_t> directions)
207 {
208
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 directions_ = directions;
209 }
210 // [[/codeblock]]
211
212 // ### Save the relevant data for plot of permeability ratio vs. Forchheimer number
213 // [[codeblock]]
214 private:
215 1 void writePermeabilityratioVsForchheimerNumber_(std::size_t dirIdx)
216 {
217 // open a logfile
218
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 3 std::ofstream logfile(dirName_[dirIdx] + "-dir-PermeabilityratioVsForchheimerNumber.dat");
219
220 // save the data needed to be plotted in logfile
221
6/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 1 times.
 31 for (int i = 0; i < apparentPermeability_[dirIdx].size(); i++)
222 {
223 // compute the Forchheimer number
224
4/8
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 10 times.
✗ Branch 11 not taken.
 40 const Scalar forchheimerNumber = darcyPermeability_[dirIdx] * forchheimerCoefficient_[dirIdx] * inertiaToViscousRatio_[dirIdx][i];
225 // ratio between apparrent permeability and darcy permeability
226
3/6
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 10 times.
✗ Branch 8 not taken.
 30 const Scalar permeabilityRatio = apparentPermeability_[dirIdx][i] / darcyPermeability_[dirIdx];
227
228
3/6
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 10 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 10 times.
✗ Branch 9 not taken.
 20 logfile << forchheimerNumber << " " << permeabilityRatio << std::endl;
229 }
230 1 }
231 // [[/codeblock]]
232
233 // ### Save the relevant data for plot of inverse of apparent permability vs. rho v / mu
234 // [[codeblock]]
235 1 void writeInversePrmeabilityVsInertiaToViscousRatio_(std::size_t dirIdx)
236 {
237 // open a logfile and write inverese of apparent permeability given by the model vs. inertial to viscous ratio (rho v / mu)
238
1/2
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
 3 std::ofstream logfile(dirName_[dirIdx] + "-dir-InversePrmeabilityVsInertiaToViscousRatio.dat");
239
240 // save the data needed to be plotted in logfile
241
6/6
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 10 times.
✓ Branch 3 taken 1 times.
✓ Branch 4 taken 10 times.
✓ Branch 5 taken 1 times.
 31 for (int i = 0; i < apparentPermeability_[dirIdx].size(); i++)
242 {
243
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
 20 const Scalar inertiaToViscousRatio = inertiaToViscousRatio_[dirIdx][i];
244
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
 20 const Scalar inverseAppPermeability = 1 / apparentPermeability_[dirIdx][i];
245
246 // compute inverse of apparent permeability using the Forchheimer permeability and coefficient
247
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
 20 const Scalar inverseAppPermeabilityForchheimer = 1 / forchheimerPermeability_[dirIdx] + inertiaToViscousRatio * forchheimerCoefficient_[dirIdx];
248
249
4/8
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 10 times.
✗ Branch 6 not taken.
✓ Branch 9 taken 10 times.
✗ Branch 10 not taken.
✓ Branch 12 taken 10 times.
✗ Branch 13 not taken.
 30 logfile << inertiaToViscousRatio << " " << 1e-12 * inverseAppPermeability << " " << 1e-12 * inverseAppPermeabilityForchheimer << std::endl;
250 }
251 1 }
252 // [[/codeblock]]
253 //
254 // ### Plot permeability ratio vs. Forchheimer number using Gnuplot
255 //
256 // [[codeblock]]
257 ✗ void plotPermeabilityratioVsForchheimerNumber_()
258 {
259 // using gnuplot interface
260 ✗ Dumux::GnuplotInterface<Scalar> gnuplot(true);
261 ✗ gnuplot.setOpenPlotWindow(true);
262
263 ✗ for (const auto dirIdx : directions_)
264 {
265 ✗ gnuplot.resetAll();
266 ✗ std::string title{}, option{};
267
268 // add the data in each direction for plot
269 ✗ gnuplot.addFileToPlot(dirName_[dirIdx] + "-dir-PermeabilityratioVsForchheimerNumber.dat", "notitle with lines");
270 // set the properties of lines to be plotted
271 ✗ option += "set linetype 1 linecolor 1 linewidth 7\n";
272 // report the darcy permeability in each direction as the title of the plot
273 ✗ title += Fmt::format("{}-direction, Darcy permeability= {:.3e} m^2 ", dirName_[dirIdx], darcyPermeability_[dirIdx]);
274
275 ✗ option += "set title \"" + title + "\"\n";
276 ✗ option += "set logscale x""\n";
277 ✗ option += "set format x '10^{%L}'""\n";
278
279 ✗ gnuplot.setXlabel("Forchheimer Number [-]");
280 ✗ gnuplot.setYlabel("Apparent permeability / Darcy permeability [-]");
281 ✗ gnuplot.setOption(option);
282 ✗ gnuplot.plot("permeability_ratio_versus_forchheimer_number");
283 }
284 ✗ }
285 // [[/codeblock]]
286 //
287 // ### Plot inverse of apparent permability vs. rho v / mu using Gnuplot
288 //
289 // [[codeblock]]
290 ✗ void plotInversePrmeabilityVsInertiaToViscousRatio_()
291 {
292 // using gnuplot interface
293 ✗ Dumux::GnuplotInterface<Scalar> gnuplot(true);
294 ✗ gnuplot.setOpenPlotWindow(true);
295
296 ✗ for (const auto dirIdx : directions_)
297 {
298 ✗ gnuplot.resetAll();
299 ✗ std::string title{}, option{};
300 ✗ std::string legend0 = "u 1:2 title \"Network model\" with lines";
301 // add the data in each direction for plot, first set of data
302 ✗ gnuplot.addFileToPlot(dirName_[dirIdx] + "-dir-InversePrmeabilityVsInertiaToViscousRatio.dat", legend0);
303
304 // set the properties of lines to be plotted
305 ✗ option += "set linetype 1 linecolor 1 linewidth 5\n";
306
307 ✗ std::string legend1 = "u 1:3 title \"Forchheimer equation\" with lines";
308 // add the data in each direction for plot, second set of data
309 ✗ gnuplot.addFileToPlot(dirName_[dirIdx] + "-dir-InversePrmeabilityVsInertiaToViscousRatio.dat", legend1);
310
311 // set the properties of lines to be plotted
312 ✗ option += "set linetype 2 linecolor 2 linewidth 5\n";
313
314 // report the darcy permeability in each direction as the title of the plot
315 ✗ title += Fmt::format("{}-direction, Darcy permeability= {:.3e} m^2 ", dirName_[dirIdx], darcyPermeability_[dirIdx]);
316
317 ✗ option += "set title \"" + title + "\"\n";
318 ✗ option += "set logscale x""\n";
319 ✗ option += "set format x '10^{%L}'""\n";
320
321 ✗ gnuplot.setXlabel("{/Symbol r} v / {/Symbol m} [1/m]");
322 ✗ gnuplot.setYlabel("1/ Apparent permeability [1/m^2] x 1e12");
323 ✗ gnuplot.setOption(option);
324 ✗ gnuplot.plot("inverse_apppermeability_versus_rhovmu-1");
325 }
326 ✗ }
327 // [[/codeblock]]
328 // [[details]] private data members
329 // [[codeblock]]
330 std::array<std::vector<Scalar>, 3> samplePointsX_;
331 std::array<std::vector<Scalar>, 3> samplePointsY_;
332 std::array<std::vector<Scalar>, 3>totalPressureDrop_;
333 std::array<std::vector<Scalar>, 3> apparentVelocity_;
334 std::array<std::vector<Scalar>, 3> apparentPermeability_;
335 std::array<std::vector<Scalar>, 3> inertiaToViscousRatio_;
336 std::array<Scalar, 3> darcyPermeability_;
337 std::array<Scalar, 3> forchheimerPermeability_;
338 std::array<Scalar, 3> forchheimerCoefficient_;
339 const std::array<std::string, 3> dirName_ = {"X", "Y", "Z"};
340 std::vector<std::size_t> directions_;
341 };
342
343 } // end namespace Dumux
344 // [[/codeblock]]
345 // [[/details]]
346 // [[/content]]
347 #endif
348