GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/common/math.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 137 151 90.7%
Functions: 20 23 87.0%
Branches: 157 209 75.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-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
5 // SPDX-License-Identifier: GPL-3.0-or-later
6 //
7 /*!
8 * \file
9 * \ingroup Core
10 * \brief Define some often used mathematical functions
11 */
12 #ifndef DUMUX_MATH_HH
13 #define DUMUX_MATH_HH
14
15 #include <algorithm>
16 #include <numeric>
17 #include <cmath>
18 #include <utility>
19
20 #include <dune/common/typetraits.hh>
21 #include <dune/common/fvector.hh>
22 #include <dune/common/fmatrix.hh>
23 #include <dune/common/dynmatrix.hh>
24 #include <dune/common/float_cmp.hh>
25
26 namespace Dumux {
27
28 /*!
29 * \ingroup Core
30 * \brief Calculate the (weighted) arithmetic mean of two scalar values.
31 *
32 * \param x The first input value
33 * \param y The second input value
34 * \param wx The first weight
35 * \param wy The second weight
36 */
37 template <class Scalar>
38 constexpr Scalar arithmeticMean(Scalar x, Scalar y, Scalar wx = 1.0, Scalar wy = 1.0) noexcept
39 {
40 static_assert(Dune::IsNumber<Scalar>::value, "The arguments x, y, wx, and wy have to be numbers!");
41
42
4/4
✓ Branch 0 taken 125500978 times.
✓ Branch 1 taken 608176 times.
✓ Branch 2 taken 110665516 times.
✓ Branch 3 taken 566378 times.
 236785052 if (x*y <= 0)
43 return 0;
44 236166494 return (x * wx + y * wy)/(wx + wy);
45 }
46
47 /*!
48 * \ingroup Core
49 * \brief Calculate the (weighted) harmonic mean of two scalar values.
50 *
51 * \param x The first input value
52 * \param y The second input value
53 * \param wx The first weight
54 * \param wy The second weight
55 */
56 template <class Scalar>
57 constexpr Scalar harmonicMean(Scalar x, Scalar y, Scalar wx = 1.0, Scalar wy = 1.0) noexcept
58 {
59 static_assert(Dune::IsNumber<Scalar>::value, "The arguments x, y, wx, and wy have to be numbers!");
60
61
14/18
✓ Branch 0 taken 1032633713 times.
✓ Branch 1 taken 1010558 times.
✓ Branch 2 taken 678253082 times.
✓ Branch 3 taken 1036604 times.
✓ Branch 4 taken 183573648 times.
✓ Branch 5 taken 2020608 times.
✓ Branch 6 taken 274976364 times.
✓ Branch 7 taken 10000 times.
✓ Branch 8 taken 22986628 times.
✓ Branch 9 taken 10000 times.
✓ Branch 10 taken 11335942 times.
✗ Branch 11 not taken.
✓ Branch 12 taken 4286880 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 2373056 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 2362080 times.
✗ Branch 17 not taken.
 2216869163 if (x*y <= 0)
62 return 0;
63 2212781393 return (wx + wy) * x * y / (wy * x + wx * y);
64 }
65
66 /*!
67 * \ingroup Core
68 * \brief Calculate the geometric mean of two scalar values.
69 *
70 * \param x The first input value
71 * \param y The second input value
72 * \note as std::sqrt is not constexpr this function is not constexpr
73 */
74 template <class Scalar>
75 Scalar geometricMean(Scalar x, Scalar y) noexcept
76 {
77 static_assert(Dune::IsNumber<Scalar>::value, "The arguments x and y have to be numbers!");
78
79 if (x*y <= 0)
80 return 0;
81 using std::sqrt;
82 return sqrt(x*y)*sign(x);
83 }
84
85 /*!
86 * \ingroup Core
87 * \brief Calculate the harmonic mean of a fixed-size matrix.
88 *
89 * This is done by calculating the harmonic mean for each entry
90 * individually. The result is stored the first argument.
91 *
92 * \param K The matrix for storing the result
93 * \param Ki The first input matrix
94 * \param Kj The second input matrix
95 */
96 template <class Scalar, int m, int n>
97 void harmonicMeanMatrix(Dune::FieldMatrix<Scalar, m, n> &K,
98 const Dune::FieldMatrix<Scalar, m, n> &Ki,
99 const Dune::FieldMatrix<Scalar, m, n> &Kj)
100 {
101 for (int rowIdx=0; rowIdx < m; rowIdx++){
102 for (int colIdx=0; colIdx< n; colIdx++){
103 if (Dune::FloatCmp::ne<Scalar>(Ki[rowIdx][colIdx], Kj[rowIdx][colIdx])) {
104 K[rowIdx][colIdx] =
105 harmonicMean(Ki[rowIdx][colIdx],
106 Kj[rowIdx][colIdx]);
107 }
108 else
109 K[rowIdx][colIdx] = Ki[rowIdx][colIdx];
110 }
111 }
112 }
113
114 /*!
115 * \ingroup Core
116 * \brief A smoothed minimum function (using cubic interpolation)
117 * \note The returned value is guaranteed to be smaller or equal to the minimum of the two input values.
118 * \note The smoothing kicks in if the difference between the two input values is smaller than the smoothing parameter.
119 * \param a The first input value
120 * \param b The second input value
121 * \param k The smoothing parameter (k > 0)
122 */
123 template<class Scalar>
124 Scalar smoothMin(const Scalar a, const Scalar b, const Scalar k)
125 {
126 using std::max; using std::min; using std::abs;
127
8/8
✓ Branch 0 taken 1681 times.
✓ Branch 1 taken 150 times.
✓ Branch 2 taken 1681 times.
✓ Branch 3 taken 150 times.
✓ Branch 4 taken 199 times.
✓ Branch 5 taken 18 times.
✓ Branch 6 taken 199 times.
✓ Branch 7 taken 18 times.
 4098 const auto h = max(k-abs(a-b), 0.0 )/k;
128
4/4
✓ Branch 0 taken 1110 times.
✓ Branch 1 taken 721 times.
✓ Branch 2 taken 164 times.
✓ Branch 3 taken 53 times.
 3324 return min(a, b) - h*h*h*k*(1.0/6.0);
129 }
130
131 /*!
132 * \ingroup Core
133 * \brief A smoothed maximum function (using cubic interpolation)
134 *
135 * \param a The first input value
136 * \param b The second input value
137 * \param k The smoothing parameter (k > 0)
138 */
139 template<class Scalar>
140 Scalar smoothMax(const Scalar a, const Scalar b, const Scalar k)
141 4 { return -smoothMin(-a, -b, k); }
142
143 /*!
144 * \ingroup Core
145 * \brief Invert a linear polynomial analytically
146 *
147 * The polynomial is defined as
148 * \f[ p(x) = a\; x + b \f]
149 *
150 * This method Returns the number of solutions which are in the real
151 * numbers, i.e. 1 except if the slope of the line is 0.
152 *
153 * \param sol Container into which the solutions are written
154 * \param a The coefficiont for the linear term
155 * \param b The coefficiont for the constant term
156 */
157 template <class Scalar, class SolContainer>
158 int invertLinearPolynomial(SolContainer &sol,
159 Scalar a,
160 Scalar b)
161 {
162 600000 if (a == 0.0)
163 return 0;
164
165 600000 sol[0] = -b/a;
166 return 1;
167 }
168
169 /*!
170 * \ingroup Core
171 * \brief Invert a quadratic polynomial analytically
172 *
173 * The polynomial is defined as
174 * \f[ p(x) = a\; x^2 + + b\;x + c \f]
175 *
176 * This method returns the number of solutions which are in the real
177 * numbers. The "sol" argument contains the real roots of the parabola
178 * in order with the smallest root first.
179 *
180 * \param sol Container into which the solutions are written
181 * \param a The coefficiont for the quadratic term
182 * \param b The coefficiont for the linear term
183 * \param c The coefficiont for the constant term
184 */
185 template <class Scalar, class SolContainer>
186 1200000 int invertQuadraticPolynomial(SolContainer &sol,
187 Scalar a,
188 Scalar b,
189 Scalar c)
190 {
191 // check for a line
192
2/2
✓ Branch 0 taken 600000 times.
✓ Branch 1 taken 600000 times.
 1200000 if (a == 0.0)
193
1/2
✓ Branch 0 taken 600000 times.
✗ Branch 1 not taken.
 600000 return invertLinearPolynomial(sol, b, c);
194
195 // discriminant
196 600000 Scalar Delta = b*b - 4*a*c;
197
1/2
✓ Branch 0 taken 600000 times.
✗ Branch 1 not taken.
 600000 if (Delta < 0)
198 return 0; // no real roots
199
200 using std::sqrt;
201 600000 Delta = sqrt(Delta);
202 600000 sol[0] = (- b + Delta)/(2*a);
203 600000 sol[1] = (- b - Delta)/(2*a);
204
205 // sort the result
206
2/2
✓ Branch 0 taken 400000 times.
✓ Branch 1 taken 200000 times.
 600000 if (sol[0] > sol[1])
207 {
208 using std::swap;
209 400000 swap(sol[0], sol[1]);
210 }
211 return 2; // two real roots
212 }
213
214 //! \cond false
215 template <class Scalar, class SolContainer>
216 3801593 void invertCubicPolynomialPostProcess_(SolContainer &sol,
217 int numSol,
218 Scalar a, Scalar b, Scalar c, Scalar d,
219 std::size_t numIterations = 1)
220 {
221 28525529 const auto eval = [&](auto x){ return d + x*(c + x*(b + x*a)); };
222 20213642 const auto evalDeriv = [&](auto x){ return c + x*(2*b + x*3*a); };
223
224 // do numIterations Newton iterations on the analytic solution
225 // and update result if the precision is increased
226
2/2
✓ Branch 0 taken 8311887 times.
✓ Branch 1 taken 3801593 times.
 12113480 for (int i = 0; i < numSol; ++i)
227 {
228 8311887 Scalar x = sol[i];
229 16623774 Scalar fCurrent = eval(x);
230 8311887 const Scalar fOld = fCurrent;
231
2/2
✓ Branch 0 taken 16412049 times.
✓ Branch 1 taken 8311887 times.
 24723936 for (int j = 0; j < numIterations; ++j)
232 {
233 32824098 const Scalar fPrime = evalDeriv(x);
234
1/2
✓ Branch 0 taken 16412049 times.
✗ Branch 1 not taken.
 16412049 if (fPrime == 0.0)
235 break;
236 16412049 x -= fCurrent/fPrime;
237 16412049 fCurrent = eval(x);
238 }
239
240 using std::abs;
241
6/6
✓ Branch 0 taken 5208773 times.
✓ Branch 1 taken 3103114 times.
✓ Branch 2 taken 5208773 times.
✓ Branch 3 taken 3103114 times.
✓ Branch 4 taken 5208773 times.
✓ Branch 5 taken 3103114 times.
 24935661 if (abs(fCurrent) < abs(fOld))
242 5208773 sol[i] = x;
243 }
244 3801593 }
245 //! \endcond
246
247 /*!
248 * \ingroup Core
249 * \brief Invert a cubic polynomial analytically
250 *
251 * The polynomial is defined as
252 * \f[ p(x) = a\; x^3 + + b\;x^3 + c\;x + d \f]
253 *
254 * This method returns the number of solutions which are in the real
255 * numbers. The "sol" argument contains the real roots of the cubic
256 * polynomial in order with the smallest root first.
257 *
258 * \note The closer the roots are to each other the less
259 * precise the inversion becomes. Increase number of post-processing iterations for improved results.
260 *
261 * \param sol Container into which the solutions are written
262 * \param a The coefficient for the cubic term
263 * \param b The coefficient for the quadratic term
264 * \param c The coefficient for the linear term
265 * \param d The coefficient for the constant term
266 * \param numPostProcessIterations The number of iterations to increase precision of the analytical result
267 */
268 template <class Scalar, class SolContainer>
269 3801593 int invertCubicPolynomial(SolContainer *sol,
270 Scalar a, Scalar b, Scalar c, Scalar d,
271 std::size_t numPostProcessIterations = 1)
272 {
273 // reduces to a quadratic polynomial
274
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3801593 times.
 3801593 if (a == 0)
275 1200000 return invertQuadraticPolynomial(sol, b, c, d);
276
277 // normalize the polynomial
278 3801593 b /= a;
279 3801593 c /= a;
280 3801593 d /= a;
281 3801593 a = 1;
282
283 // get rid of the quadratic term by substituting x = t - b/3
284 3801593 Scalar p = c - b*b/3;
285 3801593 Scalar q = d + (2*b*b*b - 9*b*c)/27;
286
287 // now we are at the form t^3 + p*t + q = 0. First we handle some
288 // special cases to avoid divisions by zero later...
289
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3801593 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 3801593 if (p == 0.0 && q == 0.0) {
290 // t^3 = 0, i.e. triple root at t = 0
291 ✗ sol[0] = sol[1] = sol[2] = 0.0 - b/3;
292 ✗ return 3;
293 }
294
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 3801593 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 3801593 else if (p == 0.0 && q != 0.0) {
295 // t^3 + q = 0,
296 //
297 // i. e. single real root at t=curt(q)
298 using std::cbrt;
299 ✗ Scalar t = cbrt(q);
300 ✗ sol[0] = t - b/3;
301
302 ✗ return 1;
303 }
304
2/4
✓ Branch 0 taken 3801593 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3801593 times.
 3801593 else if (p != 0.0 && q == 0.0) {
305 // t^3 + p*t = 0 = t*(t^2 + p),
306 //
307 // i. e. roots at t = 0, t^2 + p = 0
308 ✗ if (p > 0) {
309 ✗ sol[0] = 0.0 - b/3;
310 ✗ return 1; // only a single real root at t=0
311 }
312
313 // two additional real roots at t = sqrt(-p) and t = -sqrt(-p)
314 using std::sqrt;
315 ✗ sol[0] = -sqrt(-p) - b/3;
316 ✗ sol[1] = 0.0 - b/3;
317 ✗ sol[2] = sqrt(-p) - b/3;
318
319 ✗ return 3;
320 }
321
322 // At this point
323 //
324 // t^3 + p*t + q = 0
325 //
326 // with p != 0 and q != 0 holds. Introducing the variables u and v
327 // with the properties
328 //
329 // u + v = t and 3*u*v + p = 0
330 //
331 // leads to
332 //
333 // u^3 + v^3 + q = 0 .
334 //
335 // multiplying both sides with u^3 and taking advantage of the
336 // fact that u*v = -p/3 leads to
337 //
338 // u^6 + q*u^3 - p^3/27 = 0
339 //
340 // Now, substituting u^3 = w yields
341 //
342 // w^2 + q*w - p^3/27 = 0
343 //
344 // This is a quadratic equation with the solutions
345 //
346 // w = -q/2 + sqrt(q^2/4 + p^3/27) and
347 // w = -q/2 - sqrt(q^2/4 + p^3/27)
348 //
349 // Since w is equivalent to u^3 it is sufficient to only look at
350 // one of the two cases. Then, there are still 2 cases: positive
351 // and negative discriminant.
352 3801593 Scalar wDisc = q*q/4 + p*p*p/27;
353
2/2
✓ Branch 0 taken 1546446 times.
✓ Branch 1 taken 2255147 times.
 3801593 if (wDisc >= 0) { // the positive discriminant case:
354 // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
355 using std::cbrt;
356 using std::sqrt;
357
358 // Choose the root that is safe against the loss of precision
359 // that can cause wDisc to be equal to q*q/4 despite p != 0.
360 // We do not want u to be zero in that case. Mathematically,
361 // we are happy with either root.
362 1546446 const Scalar u = [&]{
363
2/2
✓ Branch 0 taken 1504327 times.
✓ Branch 1 taken 42119 times.
 1546446 return q > 0 ? cbrt(-0.5*q - sqrt(wDisc)) : cbrt(-0.5*q + sqrt(wDisc));
364 1546446 }();
365 // at this point, u != 0 since p^3 = 0 is necessary in order
366 // for u = 0 to hold, so
367 1546446 sol[0] = u - p/(3*u) - b/3;
368 // does not produce a division by zero. the remaining two
369 // roots of u are rotated by +- 2/3*pi in the complex plane
370 // and thus not considered here
371 1546446 invertCubicPolynomialPostProcess_(sol, 1, a, b, c, d, numPostProcessIterations);
372 1546446 return 1;
373 }
374 else { // the negative discriminant case:
375 2255147 Scalar uCubedRe = - q/2;
376 using std::sqrt;
377 2255147 Scalar uCubedIm = sqrt(-wDisc);
378 // calculate the cube root of - q/2 + sqrt(q^2/4 + p^3/27)
379 using std::cbrt;
380 2255147 Scalar uAbs = cbrt(sqrt(uCubedRe*uCubedRe + uCubedIm*uCubedIm));
381 using std::atan2;
382 2255147 Scalar phi = atan2(uCubedIm, uCubedRe)/3;
383
384 // calculate the length and the angle of the primitive root
385
386 // with the definitions from above it follows that
387 //
388 // x = u - p/(3*u) - b/3
389 //
390 // where x and u are complex numbers. Rewritten in polar form
391 // this is equivalent to
392 //
393 // x = |u|*e^(i*phi) - p*e^(-i*phi)/(3*|u|) - b/3 .
394 //
395 // Factoring out the e^ terms and subtracting the additional
396 // terms, yields
397 //
398 // x = (e^(i*phi) + e^(-i*phi))*(|u| - p/(3*|u|)) - y - b/3
399 //
400 // with
401 //
402 // y = - |u|*e^(-i*phi) + p*e^(i*phi)/(3*|u|) .
403 //
404 // The crucial observation is the fact that y is the conjugate
405 // of - x + b/3. This means that after taking advantage of the
406 // relation
407 //
408 // e^(i*phi) + e^(-i*phi) = 2*cos(phi)
409 //
410 // the equation
411 //
412 // x = 2*cos(phi)*(|u| - p / (3*|u|)) - conj(x) - 2*b/3
413 //
414 // holds. Since |u|, p, b and cos(phi) are real numbers, it
415 // follows that Im(x) = - Im(x) and thus Im(x) = 0. This
416 // implies
417 //
418 // Re(x) = x = cos(phi)*(|u| - p / (3*|u|)) - b/3 .
419 //
420 // Considering the fact that u is a cubic root, we have three
421 // values for phi which differ by 2/3*pi. This allows to
422 // calculate the three real roots of the polynomial:
423
2/2
✓ Branch 0 taken 6765441 times.
✓ Branch 1 taken 2255147 times.
 9020588 for (int i = 0; i < 3; ++i) {
424 using std::cos;
425 6765441 sol[i] = cos(phi)*(uAbs - p/(3*uAbs)) - b/3;
426 6765441 phi += 2*M_PI/3;
427 }
428
429 // post process the obtained solution to increase numerical
430 // precision
431 2255147 invertCubicPolynomialPostProcess_(sol, 3, a, b, c, d, numPostProcessIterations);
432
433 // sort the result
434 using std::sort;
435 4510294 sort(sol, sol + 3);
436
437 2255147 return 3;
438 }
439
440 // NOT REACHABLE!
441 return 0;
442 }
443
444 /*!
445 * \ingroup Core
446 * \brief Comparison of two position vectors
447 *
448 * Compares an current position vector with a reference vector, and returns true
449 * if the position vector is larger.
450 * "Larger" in this case means that all the entries of each spatial dimension are
451 * larger compared to the reference vector.
452 *
453 * \param pos Vector holding the current Position that is to be checked
454 * \param smallerVec Reference vector, holding the minimum values for comparison.
455 */
456 template <class Scalar, int dim>
457 bool isLarger(const Dune::FieldVector<Scalar, dim> &pos,
458 const Dune::FieldVector<Scalar, dim> &smallerVec)
459 {
460 for (int i=0; i < dim; i++)
461 {
462 if (pos[i]<= smallerVec[i])
463 {
464 return false;
465 }
466 }
467 return true;
468 }
469
470 /*!
471 * \ingroup Core
472 * \brief Comparison of two position vectors
473 *
474 * Compares an current position vector with a reference vector, and returns true
475 * if the position vector is smaller.
476 * "Smaller" in this case means that all the entries of each spatial dimension are
477 * smaller in comparison with the reference vector.
478 *
479 * \param pos Vector holding the current Position that is to be checked
480 * \param largerVec Reference vector, holding the maximum values for comparison.
481 */
482 template <class Scalar, int dim>
483 bool isSmaller(const Dune::FieldVector<Scalar, dim> &pos,
484 const Dune::FieldVector<Scalar, dim> &largerVec)
485 {
486 for (int i=0; i < dim; i++)
487 {
488 if (pos[i]>= largerVec[i])
489 {
490 return false;
491 }
492 }
493 return true;
494 }
495
496 /*!
497 * \ingroup Core
498 * \brief Comparison of three position vectors
499 *
500 * Compares an current position vector with two reference vector, and returns true
501 * if the position vector lies in between them.
502 * "Between" in this case means that all the entries of each spatial dimension are
503 * smaller in comparison with the larger reference vector as well as larger compared
504 * to the smaller reference.
505 * This is comfortable to cheack weather the current position is located inside or
506 * outside of a lense with different properties.
507 *
508 * \param pos Vector holding the current Position that is to be checked
509 * \param smallerVec Reference vector, holding the minimum values for comparison.
510 * \param largerVec Reference vector, holding the maximum values for comparison.
511 */
512 template <class Scalar, int dim>
513 bool isBetween(const Dune::FieldVector<Scalar, dim> &pos,
514 const Dune::FieldVector<Scalar, dim> &smallerVec,
515 const Dune::FieldVector<Scalar, dim> &largerVec)
516 {
517 if (isLarger(pos, smallerVec) && isSmaller(pos, largerVec))
518 {
519 return true;
520 }
521 else
522 return false;
523 }
524
525 //! forward declaration of the linear interpolation policy (default)
526 namespace InterpolationPolicy { struct Linear; }
527
528 /*!
529 * \ingroup Core
530 * \brief a generic function to interpolate given a set of parameters and an interpolation point
531 * \param params the parameters used for interpolation (depends on the policy used)
532 * \param ip the interpolation point
533 */
534 template <class Policy = InterpolationPolicy::Linear, class Scalar, class ... Parameter>
535 Scalar interpolate(Scalar ip, Parameter&& ... params)
536
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
 51291200 { return Policy::interpolate(ip, std::forward<Parameter>(params) ...); }
537
538 /*!
539 * \ingroup Core
540 * \brief Interpolation policies
541 */
542 namespace InterpolationPolicy {
543
544 /*!
545 * \ingroup Core
546 * \brief interpolate linearly between two given values
547 */
548 struct Linear
549 {
550 /*!
551 * \brief interpolate linearly between two given values
552 * \param ip the interpolation point in [0,1]
553 * \param params array with the lower and upper bound
554 */
555 template<class Scalar>
556 static constexpr Scalar interpolate(Scalar ip, const std::array<Scalar, 2>& params)
557 {
558 48851177 return params[0]*(1.0 - ip) + params[1]*ip;
559 }
560 };
561
562 /*!
563 * \ingroup Core
564 * \brief interpolate linearly in a piecewise linear function (tabularized function)
565 */
566 struct LinearTable
567 {
568 /*!
569 * \brief interpolate linearly in a piecewise linear function (tabularized function)
570 * \param ip the interpolation point
571 * \param range positions of values
572 * \param values values to interpolate from
573 * \note if the interpolation point is out of bounds this will return the bounds
574 */
575 template<class Scalar, class RandomAccessContainer0, class RandomAccessContainer1>
576 51292377 static constexpr Scalar interpolate(Scalar ip, const RandomAccessContainer0& range, const RandomAccessContainer1& values)
577 {
578 // check bounds
579
6/6
✓ Branch 0 taken 2440317 times.
✓ Branch 1 taken 48852060 times.
✓ Branch 2 taken 2440317 times.
✓ Branch 3 taken 48852060 times.
✓ Branch 4 taken 2440317 times.
✓ Branch 5 taken 48852060 times.
 153877131 if (ip > range[range.size()-1]) return values[values.size()-1];
580
2/2
✓ Branch 0 taken 882 times.
✓ Branch 1 taken 48851178 times.
 48852060 if (ip < range[0]) return values[0];
581
582 // if we are within bounds find the index of the lower bound
583
4/4
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 48851177 times.
✓ Branch 2 taken 1 times.
✓ Branch 3 taken 48851177 times.
 97702356 const auto lookUpIndex = std::distance(range.begin(), std::lower_bound(range.begin(), range.end(), ip));
584
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 48851177 times.
 48851178 if (lookUpIndex == 0)
585 1 return values[0];
586
587 97702354 const auto ipLinear = (ip - range[lookUpIndex-1])/(range[lookUpIndex] - range[lookUpIndex-1]);
588 195404708 return Dumux::interpolate<Linear>(ipLinear, std::array<Scalar, 2>{{values[lookUpIndex-1], values[lookUpIndex]}});
589 }
590
591 template<class Scalar, class RandomAccessContainer>
592 static constexpr Scalar interpolate(Scalar ip, const std::pair<RandomAccessContainer, RandomAccessContainer>& table)
593 {
594 2352 const auto& [range, values] = table;
595 1176 return interpolate(ip, range, values);
596 }
597 };
598
599 } // end namespace InterpolationPolicy
600
601 /*!
602 * \ingroup Core
603 * \brief Generates linearly spaced vectors
604 *
605 * \param begin The first value in the vector
606 * \param end The last value in the vector
607 * \param samples The size of the vector
608 * \param endPoint if the range is including the interval's end point or not
609 */
610 template <class Scalar>
611 50 std::vector<Scalar> linspace(const Scalar begin, const Scalar end,
612 std::size_t samples,
613 bool endPoint = true)
614 {
615 using std::max;
616
1/2
✓ Branch 0 taken 50 times.
✗ Branch 1 not taken.
 50 samples = max(std::size_t{2}, samples); // only makes sense for 2 or more samples
617
2/2
✓ Branch 0 taken 49 times.
✓ Branch 1 taken 1 times.
 50 const Scalar divisor = endPoint ? samples-1 : samples;
618 50 const Scalar delta = (end-begin)/divisor;
619 100 std::vector<Scalar> vec(samples);
620
2/2
✓ Branch 0 taken 8048755 times.
✓ Branch 1 taken 50 times.
 8048805 for (std::size_t i = 0; i < samples; ++i)
621 16097510 vec[i] = begin + i*delta;
622 50 return vec;
623 }
624
625
626 /*!
627 * \ingroup Core
628 * \brief Evaluates the Antoine equation used to calculate the vapour
629 * pressure of various liquids.
630 *
631 * See http://en.wikipedia.org/wiki/Antoine_equation
632 *
633 * \param temperature The temperature [K] of the fluid
634 * \param A The first coefficient for the Antoine equation
635 * \param B The first coefficient for the Antoine equation
636 * \param C The first coefficient for the Antoine equation
637 */
638 template <class Scalar>
639 Scalar antoine(Scalar temperature,
640 Scalar A,
641 Scalar B,
642 Scalar C)
643 {
644 const Scalar ln10 = 2.3025850929940459;
645 using std::exp;
646 return exp(ln10*(A - B/(C + temperature)));
647 }
648
649 /*!
650 * \ingroup Core
651 * \brief Sign or signum function.
652 *
653 * Returns 1 for a positive argument.
654 * Returns -1 for a negative argument.
655 * Returns 0 if the argument is zero.
656 */
657 template<class ValueType>
658 ✗ constexpr int sign(const ValueType& value) noexcept
659 {
660
13/13
✓ Branch 0 taken 26 times.
✓ Branch 1 taken 158776510 times.
✓ Branch 2 taken 101242356 times.
✓ Branch 3 taken 43613912 times.
✓ Branch 4 taken 124074537 times.
✓ Branch 5 taken 69538 times.
✓ Branch 6 taken 146894738 times.
✓ Branch 7 taken 252987 times.
✓ Branch 8 taken 56753133 times.
✓ Branch 9 taken 1152762 times.
✓ Branch 10 taken 1835 times.
✓ Branch 11 taken 1003572 times.
✓ Branch 12 taken 15 times.
 1169128681 return (ValueType(0) < value) - (value < ValueType(0));
661 }
662
663 /*!
664 * \ingroup Core
665 * \brief Cross product of two vectors in three-dimensional Euclidean space
666 *
667 * \param vec1 The first vector
668 * \param vec2 The second vector
669 */
670 template <class Scalar>
671 39345780 Dune::FieldVector<Scalar, 3> crossProduct(const Dune::FieldVector<Scalar, 3> &vec1,
672 const Dune::FieldVector<Scalar, 3> &vec2)
673 {
674 157383120 return {vec1[1]*vec2[2]-vec1[2]*vec2[1],
675 118037340 vec1[2]*vec2[0]-vec1[0]*vec2[2],
676 236074680 vec1[0]*vec2[1]-vec1[1]*vec2[0]};
677 }
678
679 /*!
680 * \ingroup Core
681 * \brief Cross product of two vectors in two-dimensional Euclidean space retuning scalar
682 *
683 * \param vec1 The first vector
684 * \param vec2 The second vector
685 */
686 template <class Scalar>
687 Scalar crossProduct(const Dune::FieldVector<Scalar, 2> &vec1,
688 const Dune::FieldVector<Scalar, 2> &vec2)
689
20/20
✓ Branch 0 taken 13628205 times.
✓ Branch 1 taken 21035116 times.
✓ Branch 2 taken 13628205 times.
✓ Branch 3 taken 21035116 times.
✓ Branch 4 taken 13628205 times.
✓ Branch 5 taken 21035116 times.
✓ Branch 6 taken 13628205 times.
✓ Branch 7 taken 21035116 times.
✓ Branch 8 taken 13628205 times.
✓ Branch 9 taken 21035116 times.
✓ Branch 10 taken 200 times.
✓ Branch 11 taken 200 times.
✓ Branch 12 taken 200 times.
✓ Branch 13 taken 200 times.
✓ Branch 14 taken 200 times.
✓ Branch 15 taken 200 times.
✓ Branch 16 taken 200 times.
✓ Branch 17 taken 200 times.
✓ Branch 18 taken 200 times.
✓ Branch 19 taken 200 times.
 361779670 { return vec1[0]*vec2[1]-vec1[1]*vec2[0]; }
690
691 /*!
692 * \ingroup Core
693 * \brief Triple product of three vectors in three-dimensional Euclidean space retuning scalar
694 *
695 * \param vec1 The first vector
696 * \param vec2 The second vector
697 * \param vec3 The third vector
698 */
699 template <class Scalar>
700 Scalar tripleProduct(const Dune::FieldVector<Scalar, 3> &vec1,
701 const Dune::FieldVector<Scalar, 3> &vec2,
702 const Dune::FieldVector<Scalar, 3> &vec3)
703 10272607 { return crossProduct<Scalar>(vec1, vec2)*vec3; }
704
705 /*!
706 * \ingroup Core
707 * \brief Transpose a FieldMatrix
708 *
709 * \param M The matrix to be transposed
710 */
711 template <class Scalar, int m, int n>
712 Dune::FieldMatrix<Scalar, n, m> getTransposed(const Dune::FieldMatrix<Scalar, m, n>& M)
713 {
714 2185722 Dune::FieldMatrix<Scalar, n, m> T;
715
6/6
✓ Branch 0 taken 1753604 times.
✓ Branch 1 taken 876802 times.
✓ Branch 2 taken 38758 times.
✓ Branch 3 taken 12920 times.
✓ Branch 4 taken 2592000 times.
✓ Branch 5 taken 1296000 times.
 6570084 for (std::size_t i = 0; i < m; ++i)
716
6/6
✓ Branch 0 taken 3507208 times.
✓ Branch 1 taken 1753604 times.
✓ Branch 2 taken 116270 times.
✓ Branch 3 taken 38758 times.
✓ Branch 4 taken 5184000 times.
✓ Branch 5 taken 2592000 times.
 13191840 for (std::size_t j = 0; j < n; ++j)
717 44037390 T[j][i] = M[i][j];
718
719 return T;
720 }
721
722 /*!
723 * \ingroup Core
724 * \brief Transpose a DynamicMatrix
725 *
726 * \param M The matrix to be transposed
727 */
728 template <class Scalar>
729 Dune::DynamicMatrix<Scalar> getTransposed(const Dune::DynamicMatrix<Scalar>& M)
730 {
731 std::size_t rows_T = M.M();
732 std::size_t cols_T = M.N();
733
734 Dune::DynamicMatrix<Scalar> M_T(rows_T, cols_T, 0.0);
735
736 for (std::size_t i = 0; i < rows_T; ++i)
737 for (std::size_t j = 0; j < cols_T; ++j)
738 M_T[i][j] = M[j][i];
739
740 return M_T;
741 }
742
743 /*!
744 * \ingroup Core
745 * \brief Multiply two dynamic matrices
746 *
747 * \param M1 The first dynamic matrix
748 * \param M2 The second dynamic matrix (to be multiplied to M1 from the right side)
749 */
750 template <class Scalar>
751 52461048 Dune::DynamicMatrix<Scalar> multiplyMatrices(const Dune::DynamicMatrix<Scalar> &M1,
752 const Dune::DynamicMatrix<Scalar> &M2)
753 {
754 using size_type = typename Dune::DynamicMatrix<Scalar>::size_type;
755 52461048 const size_type rows = M1.N();
756 52461048 const size_type cols = M2.M();
757
758 DUNE_ASSERT_BOUNDS(M1.M() == M2.N());
759
760 52461048 Dune::DynamicMatrix<Scalar> result(rows, cols, 0.0);
761
2/2
✓ Branch 1 taken 207548397 times.
✓ Branch 2 taken 52461048 times.
 260009445 for (size_type i = 0; i < rows; i++)
762
2/2
✓ Branch 0 taken 822834623 times.
✓ Branch 1 taken 207548397 times.
 1030383020 for (size_type j = 0; j < cols; j++)
763
2/2
✓ Branch 1 taken 3395029373 times.
✓ Branch 2 taken 822834623 times.
 4217863996 for (size_type k = 0; k < M1.M(); k++)
764 23765205611 result[i][j] += M1[i][k]*M2[k][j];
765
766 52461048 return result;
767 }
768
769 /*!
770 * \ingroup Core
771 * \brief Multiply two field matrices
772 *
773 * \param M1 The first field matrix
774 * \param M2 The second field matrix (to be multiplied to M1 from the right side)
775 */
776 template <class Scalar, int rows1, int cols1, int cols2>
777 706929 Dune::FieldMatrix<Scalar, rows1, cols2> multiplyMatrices(const Dune::FieldMatrix<Scalar, rows1, cols1> &M1,
778 const Dune::FieldMatrix<Scalar, cols1, cols2> &M2)
779 {
780 using size_type = typename Dune::FieldMatrix<Scalar, rows1, cols2>::size_type;
781
782 706929 Dune::FieldMatrix<Scalar, rows1, cols2> result(0.0);
783
2/2
✓ Branch 0 taken 2827716 times.
✓ Branch 1 taken 706929 times.
 3534645 for (size_type i = 0; i < rows1; i++)
784
2/2
✓ Branch 0 taken 11310864 times.
✓ Branch 1 taken 2827716 times.
 14138580 for (size_type j = 0; j < cols2; j++)
785
2/2
✓ Branch 0 taken 45243456 times.
✓ Branch 1 taken 11310864 times.
 56554320 for (size_type k = 0; k < cols1; k++)
786 316704192 result[i][j] += M1[i][k]*M2[k][j];
787
788 706929 return result;
789 }
790
791
792 /*!
793 * \ingroup Core
794 * \brief Trace of a dense matrix
795 *
796 * \param M The dense matrix
797 */
798 template <class MatrixType>
799 typename Dune::DenseMatrix<MatrixType>::field_type
800 trace(const Dune::DenseMatrix<MatrixType>& M)
801 {
802 1 const auto rows = M.N();
803 DUNE_ASSERT_BOUNDS(rows == M.M()); // rows == cols
804
805 using MatType = Dune::DenseMatrix<MatrixType>;
806 1 typename MatType::field_type trace = 0.0;
807
808
4/4
✓ Branch 0 taken 1204579 times.
✓ Branch 1 taken 3585729 times.
✓ Branch 2 taken 651 times.
✓ Branch 3 taken 1297 times.
 4792255 for (typename MatType::size_type i = 0; i < rows; ++i)
809 10761090 trace += M[i][i];
810
811 return trace;
812 }
813
814 /*!
815 * \ingroup Core
816 * \brief Returns the result of the projection of
817 * a vector v with a Matrix M.
818 *
819 * Note: We use DenseVector and DenseMatrix here so that
820 * it can be used with the statically and dynamically
821 * allocated Dune Vectors/Matrices. Size mismatch
822 * assertions are done in the respective Dune classes.
823 *
824 * \param M The matrix
825 * \param v The vector
826 */
827 template <class MAT, class V>
828 typename Dune::DenseVector<V>::derived_type
829 85257460 mv(const Dune::DenseMatrix<MAT>& M,
830 const Dune::DenseVector<V>& v)
831 {
832
2/4
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
 93716269 typename Dune::DenseVector<V>::derived_type res(v);
833 93716262 M.mv(v, res);
834 85257460 return res;
835 }
836
837 /*!
838 * \ingroup Core
839 * \brief Returns the result of a vector v multiplied by a scalar m.
840 *
841 * Note: We use DenseVector and DenseMatrix here so that
842 * it can be used with the statically and dynamically
843 * allocated Dune Vectors/Matrices. Size mismatch
844 * assertions are done in the respective Dune classes.
845 *
846 * \note We need the enable_if to make sure that only Scalars
847 * fit here. Matrix types are forced to use the above
848 * mv(DenseMatrix, DenseVector) instead.
849 *
850 * \param m The scale factor
851 * \param v The vector
852 */
853
854 template <class FieldScalar, class V>
855 typename std::enable_if_t<Dune::IsNumber<FieldScalar>::value,
856 typename Dune::DenseVector<V>::derived_type>
857 1764 mv(const FieldScalar m, const Dune::DenseVector<V>& v)
858 {
859 1766 typename Dune::DenseVector<V>::derived_type res(v);
860 1764 res *= m;
861 1764 return res;
862 }
863
864 /*!
865 * \ingroup Core
866 * \brief Evaluates the scalar product of a vector v2, projected by
867 * a matrix M, with a vector v1.
868 *
869 * Note: We use DenseVector and DenseMatrix here so that
870 * it can be used with the statically and dynamically
871 * allocated Dune Vectors/Matrices. Size mismatch
872 * assertions are done in the respective Dune classes.
873 *
874 * \param v1 The first vector
875 * \param M The matrix
876 * \param v2 The second vector
877 */
878 template <class V1, class MAT, class V2>
879 typename Dune::DenseMatrix<MAT>::value_type
880 6 vtmv(const Dune::DenseVector<V1>& v1,
881 const Dune::DenseMatrix<MAT>& M,
882 const Dune::DenseVector<V2>& v2)
883 {
884
5/10
✓ Branch 1 taken 3988534 times.
✓ Branch 2 taken 83478 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1 times.
✗ Branch 12 not taken.
✓ Branch 13 taken 1 times.
 104267742 return v1*mv(M, v2);
885 }
886
887 /*!
888 * \ingroup Core
889 * \brief Evaluates the scalar product of a vector v2, scaled by
890 * a scalar m, with a vector v1.
891 *
892 * Note: We use DenseVector and DenseMatrix here so that
893 * it can be used with the statically and dynamically
894 * allocated Dune Vectors/Matrices. Size mismatch
895 * assertions are done in the respective Dune classes.
896 *
897 * \note We need the enable_if to make sure that only Scalars
898 * fit here. Matrix types are forced to use the above
899 * vtmv(DenseVector, DenseMatrix, DenseVector) instead.
900 *
901 * \param v1 The first vector
902 * \param m The scale factor
903 * \param v2 The second vector
904 */
905
906 template <class V1, class FieldScalar, class V2>
907 typename std::enable_if_t<Dune::IsNumber<FieldScalar>::value, FieldScalar>
908 vtmv(const Dune::DenseVector<V1>& v1,
909 const FieldScalar m,
910 const Dune::DenseVector<V2>& v2)
911 {
912
8/12
✓ Branch 0 taken 466732665 times.
✓ Branch 1 taken 7888094 times.
✓ Branch 2 taken 154976125 times.
✓ Branch 3 taken 1679216 times.
✓ Branch 4 taken 875901 times.
✓ Branch 5 taken 1083504 times.
✓ Branch 6 taken 770556 times.
✓ Branch 7 taken 49896 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
 2759015103 return m*(v1*v2);
913 }
914
915 /*!
916 * \ingroup Core
917 * \brief Returns the [intercept, slope] of the regression line
918 * fitted to a set of (x, y) data points.
919 *
920 * Note: We use least-square regression method to find
921 * the regression line.
922 *
923 * \param x x-values of the data set
924 * \param y y-values of the data set
925 */
926 template <class Scalar>
927 3 std::array<Scalar, 2> linearRegression(const std::vector<Scalar>& x,
928 const std::vector<Scalar>& y)
929 {
930
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3 times.
 9 if (x.size() != y.size())
931 ✗ DUNE_THROW(Dune::InvalidStateException, "x and y array must have the same length.");
932
933 3 const Scalar averageX = std::accumulate(x.begin(), x.end(), 0.0)/x.size();
934 9 const Scalar averageY = std::accumulate(y.begin(), y.end(), 0.0)/y.size();
935
936 // calculate temporary variables necessary for slope computation
937 12 const Scalar numerator = std::inner_product(
938 x.begin(), x.end(), y.begin(), 0.0, std::plus<Scalar>(),
939 208 [&](auto xx, auto yy) { return (xx - averageX) * (yy - averageY); }
940 );
941 3 const Scalar denominator = std::inner_product(
942 x.begin(), x.end(), x.begin(), 0.0, std::plus<Scalar>(),
943 208 [&](auto xx, auto yy) { return (xx - averageX) * (yy - averageX); }
944 );
945
946 // compute slope and intercept of the regression line
947 3 const Scalar slope = numerator / denominator;
948 3 const Scalar intercept = averageY - slope * averageX;
949
950 3 return {intercept, slope};
951 }
952
953 } // end namespace Dumux
954
955 #endif
956