]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add results section
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 30 Dec 2019 21:16:21 +0000 (22:16 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 24 Mar 2020 06:59:18 +0000 (07:59 +0100)
examples/step-67/doc/results.dox

index 207317d6cecebe195498a44dd8bc267955ccc985..7a8f6d9289a9fc8a3d7504cb5cb52fc460d8f13e 100644 (file)
@@ -1,4 +1,765 @@
 <h1>Results</h1>
 
-Running the program ...
+<h3>Program output</h3>
 
+Running the program with the default settings on a machine with 40 processes
+produces the following output:
+@code
+Running with 40 MPI processes
+Vectorization over 8 doubles = 512 bits (AVX512), VECTORIZATION_LEVEL=3
+Number of degrees of freedom: 147,456 ( = 4 [vars] x 1,024 [cells] x 36 [dofs/cell/var] )
+Time step size: 0.00689325, minimal h: 0.3125, initial transport scaling: 0.102759
+
+Time:       0 , dt:   0.0069 , error rho:   2.76e-07 , rho * u:  1.259e-06 , energy: 2.987e-06
+Time:    1.01 , dt:   0.0069 , error rho:   1.37e-06 , rho * u:  2.252e-06 , energy: 4.153e-06
+Time:    2.01 , dt:   0.0069 , error rho:  1.561e-06 , rho * u:   2.43e-06 , energy: 4.493e-06
+Time:    3.01 , dt:   0.0069 , error rho:  1.714e-06 , rho * u:  2.591e-06 , energy: 4.762e-06
+Time:    4.01 , dt:   0.0069 , error rho:  1.843e-06 , rho * u:  2.625e-06 , energy: 4.985e-06
+Time:    5.01 , dt:   0.0069 , error rho:  1.496e-06 , rho * u:  1.961e-06 , energy: 4.142e-06
+Time:       6 , dt:   0.0083 , error rho:  1.007e-06 , rho * u:  7.119e-07 , energy: 2.972e-06
+Time:       7 , dt:   0.0095 , error rho:  9.096e-07 , rho * u:  3.786e-07 , energy: 2.626e-06
+Time:       8 , dt:   0.0096 , error rho:  8.439e-07 , rho * u:  3.338e-07 , energy:  2.43e-06
+Time:       9 , dt:   0.0096 , error rho:  7.822e-07 , rho * u:  2.984e-07 , energy: 2.248e-06
+Time:      10 , dt:   0.0096 , error rho:  7.231e-07 , rho * u:  2.666e-07 , energy: 2.074e-06
+
++-------------------------------------------+------------------+------------+------------------+
+| Total wallclock time elapsed              |     2.249s    30 |     2.249s |     2.249s     8 |
+|                                           |                  |                               |
+| Section                       | no. calls |   min time  rank |   avg time |   max time  rank |
++-------------------------------------------+------------------+------------+------------------+
+| compute errors                |        11 |  0.008066s    13 |   0.00952s |   0.01041s    20 |
+| compute transport speed       |       258 |   0.01012s    13 |   0.05392s |   0.08574s    25 |
+| output                        |        11 |    0.9597s    13 |    0.9613s |    0.9623s     6 |
+| rk time step total            |      1283 |    0.9827s    25 |     1.015s |      1.06s    13 |
+| rk_stage - integrals L_h      |      6415 |    0.8803s    26 |    0.9198s |    0.9619s    14 |
+| rk_stage - inv mass + vec upd |      6415 |   0.05677s    15 |   0.06487s |   0.07597s    13 |
++-------------------------------------------+------------------+------------+------------------+
+@endcode
+
+The program output shows that all errors are small. This is due to the fact
+that we use a relatively fine mesh of $32^2$ cells with polynomials of degree
+5 for a solution that is smooth. An interesting pattern shows for the time
+step size: whereas it is 0.0069 up to time 5, it increases to 0.0096 for later
+times. The step size increases once the vortex with some motion on top of the
+speed of sound (and thus faster propagation) leaves the computational domain
+between times 5 and 6.5. Our time step formula recognizes this and only
+detects the acoustic limit in the last part of the simulation.
+
+The summary of wall clock times shows that 1283 time steps have been performed
+in 1.02 seconds (looking at the average time among all MPI processes), while
+the output of 11 files has taken additional 0.96 seconds. Broken down per time
+step and into the five Runge--Kutta stages, the compute time per evaluation is
+0.16 milliseconds. This high performance is typical of matrix-free evaluators
+and a reason why explicit time integration is very competitive against
+implicit solvers, especially for large-scale simulations. The breakdown of
+computational times at the end of the program run shows that the evaluation of
+integrals in $\mathcal L_h$ contributes with around 0.92 seconds and the
+application of the inverse mass matrix with 0.06 seconds. Furthermore, the
+estimation of the transport speed for the time step size computation
+contributes with another 0.05 seconds of compute time.
+
+If we use three more levels of global refinement and 9.4 million DoFs in total,
+the final statistics are as follows (for the modified Lax--Friedrichs flux,
+$p=5$, and the same system of 40 cores of dual-socket Intel Xeon Gold 6230):
+@code
++-------------------------------------------+------------------+------------+------------------+
+| Total wallclock time elapsed              |     244.9s    12 |     244.9s |     244.9s    34 |
+|                                           |                  |                               |
+| Section                       | no. calls |   min time  rank |   avg time |   max time  rank |
++-------------------------------------------+------------------+------------+------------------+
+| compute errors                |        11 |    0.4239s    12 |    0.4318s |    0.4408s     9 |
+| compute transport speed       |      2053 |     3.962s    12 |     6.727s |     10.12s     7 |
+| output                        |        11 |     30.35s    12 |     30.36s |     30.37s     9 |
+| rk time step total            |     10258 |     201.7s     7 |     205.1s |     207.8s    12 |
+| rk_stage - integrals L_h      |     51290 |     121.3s     6 |     126.6s |     136.3s    16 |
+| rk_stage - inv mass + vec upd |     51290 |     66.19s    16 |     77.52s |     81.84s    10 |
++-------------------------------------------+------------------+------------+------------------+
+@endcode
+
+Per time step, the solver now takes 0.02 seconds, about 25 times as long as
+for the small problem with 147k unknowns. Given that the problem involves 64
+times as many unknowns, the increase in computing time is not
+surprising. Since we also do 8 times as many time steps, the compute time
+should in theory increase by a factor of 512. The actual increase is 205 s /
+1.02 s = 202. This is because the small problem size cannot fully utilize the
+40 cores due to communication overhead. This becomes clear if we look into the
+details of the operations done per time step. The evaluation of the
+differential operator $\mathcal L_h$ with nearest neighbor communication goes
+from 0.92 seconds to 127 seconds, i.e., it increases with a factor of 138. On
+the other hand, the cost for application of the inverse mass matrix and the
+vector updates, which do not need to communicate between the MPI processes at
+all, has increased by a factor of 1195. The increase is more than the
+theoretical factor of 512 because the operation is limited by the bandwidth
+from RAM memory for the larger size while for the smaller size, all vectors
+fit into the caches of the CPU. The numbers show that the mass matrix
+evaluation and vector update part consume almost 40% of the time spent by the
+Runge--Kutta stages -- despite using a low-storage Runge--Kutta integrator and
+merging of vector operations! And despite using over-integration for the
+$\mathcal L_h$ operator. For simpler differential operators and more expensive
+time integrators, the proportion spent in the mass matrix and vector update
+part can also reach 70%. If we compute a throughput number in terms of DoFs
+processed per second and Runge--Kutta stage, we obtain @f[ \text{throughput} =
+\frac{n_\mathrm{time steps} n_\mathrm{stages}
+n_\mathrm{dofs}}{t_\mathrm{compute}} = \frac{10258 \cdot 5 \cdot
+9.4\,\text{MDoFs}}{205s} = 2360\, \text{MDoFs/s} @f] This throughput number is
+very high, given that simply copying one vector to another one runs at around
+10,000 MDoFs/s.
+
+If we go to the next-larger size with 37.7 million DoFs, the overall
+simulation time is 2196 seconds, with 1978 seconds spent in the time
+stepping. The increase in run time is a factor of 9.3 for the L_h operator
+(1179 versus 127 seconds) and a factor of 10.3 for the inverse mass matrix and
+vector updates (797 vs 77.5 seconds). The reason for this non-optimal increase
+in run time can be traced back to cache effects on the given hardware (with 40
+MB of L2 cache and 55 MB of L3 cache): While not all of the relevant data fits
+into caches for 9.4 million DoFs (one vector takes 75 MB and we have three
+vectors plus some additional data in MatrixFree), there is capacity for almost
+half of it nonetheless. Given that modern caches are more sophisticated than
+the naive least-recently-used strategy (where we would have little re-use as
+the data is used in a streaming-like fashion), we can assume that a sizeable
+fraction of data can indeed be delivered from caches for the 9.4 million DoFs
+case. For the larger case, even with optimal caching less than 10 percent of
+data would fit into caches, with an associated loss in performance.
+
+<h3>Convergence rates for the analytical test case</h3>
+
+For the modified Lax--Friedrichs flux and measuring the error in the momentum
+variable, we obtain the following convergence table (the rates are very
+similar for the density and energy variables):
+
+<table align="center" class="doxtable">
+  <tr>
+    <th>&nbsp;</th>
+    <th colspan="3"><i>p</i>=2</th>
+    <th colspan="3"><i>p</i>=3</th>
+    <th colspan="3"><i>p</i>=5</th>
+  </tr>
+  <tr>
+    <th>n_cells</th>
+    <th>n_dofs</th>
+    <th>error mom</th>
+    <th>rate</th>
+    <th>n_dofs</th>
+    <th>error mom</th>
+    <th>rate</th>
+    <th>n_dofs</th>
+    <th>error mom</th>
+    <th>rate</th>
+  </tr>
+  <tr>
+    <td align="right">16</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td align="center">2,304</td>
+    <td align="center">1.373e-01</td>
+    <td>&nbsp;</td>
+  </tr>
+  <tr>
+    <td align="right">64</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td align="center">4,096</td>
+    <td align="center">9.130e-02</td>
+    <td>&nbsp;</td>
+    <td align="center">9,216</td>
+    <td align="center">8.899e-03</td>
+    <td>3.94</td>
+  </tr>
+  <tr>
+    <td align="right">256</td>
+    <td align="center">9,216</td>
+    <td align="center">5.577e-02</td>
+    <td>&nbsp;</td>
+    <td align="center">16,384</td>
+    <td align="center">7.381e-03</td>
+    <td>3.64</td>
+    <td align="center">36,864</td>
+    <td align="center">2.082e-04</td>
+    <td>5.42</td>
+  </tr>
+  <tr>
+    <td align="right">1024</td>
+    <td align="center">36,864</td>
+    <td align="center">4.724e-03</td>
+    <td>3.56</td>
+    <td align="center">65,536</td>
+    <td align="center">3.072e-04</td>
+    <td>4.59</td>
+    <td align="center">147,456</td>
+    <td align="center">2.625e-06</td>
+    <td>6.31</td>
+  </tr>
+  <tr>
+    <td align="right">4096</td>
+    <td align="center">147,456</td>
+    <td align="center">6.205e-04</td>
+    <td>2.92</td>
+    <td align="center">262,144</td>
+    <td align="center">1.880e-05</td>
+    <td>4.03</td>
+    <td align="center">589,824</td>
+    <td align="center">3.268e-08</td>
+    <td>6.33</td>
+  </tr>
+  <tr>
+    <td align="right">16,384</td>
+    <td align="center">589,824</td>
+    <td align="center">8.279e-05</td>
+    <td>2.91</td>
+    <td align="center">1,048,576</td>
+    <td align="center">1.224e-06</td>
+    <td>3.94</td>
+    <td align="center">2,359,296</td>
+    <td align="center">9.252e-10</td>
+    <td>5.14</td>
+  </tr>
+  <tr>
+    <td align="right">65,536</td>
+    <td align="center">2,359,296</td>
+    <td align="center">1.105e-05</td>
+    <td>2.91</td>
+    <td align="center">4,194,304</td>
+    <td align="center">7.871e-08</td>
+    <td>3.96</td>
+    <td align="center">9,437,184</td>
+    <td align="center">1.369e-10</td>
+    <td>2.77</td>
+  </tr>
+  <tr>
+    <td align="right">262,144</td>
+    <td align="center">9,437,184</td>
+    <td align="center">1.615e-06</td>
+    <td>2.77</td>
+    <td align="center">16,777,216</td>
+    <td align="center">4.961e-09</td>
+    <td>3.99</td>
+    <td align="center">37,748,736</td>
+    <td align="center">7.091e-11</td>
+    <td>0.95</td>
+  </tr>
+</table>
+
+If we switch to the Harten-Lax-van Leer flux, the results are as follows:
+<table align="center" class="doxtable">
+  <tr>
+    <th>&nbsp;</th>
+    <th colspan="3"><i>p</i>=2</th>
+    <th colspan="3"><i>p</i>=3</th>
+    <th colspan="3"><i>p</i>=5</th>
+  </tr>
+  <tr>
+    <th>n_cells</th>
+    <th>n_dofs</th>
+    <th>error mom</th>
+    <th>rate</th>
+    <th>n_dofs</th>
+    <th>error mom</th>
+    <th>rate</th>
+    <th>n_dofs</th>
+    <th>error mom</th>
+    <th>rate</th>
+  </tr>
+  <tr>
+    <td align="right">16</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td align="center">2,304</td>
+    <td align="center">1.339e-01</td>
+    <td>&nbsp;</td>
+  </tr>
+  <tr>
+    <td align="right">64</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td>&nbsp;</td>
+    <td align="center">4,096</td>
+    <td align="center">9.037e-02</td>
+    <td>&nbsp;</td>
+    <td align="center">9,216</td>
+    <td align="center">8.849e-03</td>
+    <td>3.92</td>
+  </tr>
+  <tr>
+    <td align="right">256</td>
+    <td align="center">9,216</td>
+    <td align="center">4.204e-02</td>
+    <td>&nbsp;</td>
+    <td align="center">16,384</td>
+    <td align="center">9.143e-03</td>
+    <td>3.31</td>
+    <td align="center">36,864</td>
+    <td align="center">2.501e-04</td>
+    <td>5.14</td>
+  </tr>
+  <tr>
+    <td align="right">1024</td>
+    <td align="center">36,864</td>
+    <td align="center">4.913e-03</td>
+    <td>3.09</td>
+    <td align="center">65,536</td>
+    <td align="center">3.257e-04</td>
+    <td>4.81</td>
+    <td align="center">147,456</td>
+    <td align="center">3.260e-06</td>
+    <td>6.26</td>
+  </tr>
+  <tr>
+    <td align="right">4096</td>
+    <td align="center">147,456</td>
+    <td align="center">7.862e-04</td>
+    <td>2.64</td>
+    <td align="center">262,144</td>
+    <td align="center">1.588e-05</td>
+    <td>4.36</td>
+    <td align="center">589,824</td>
+    <td align="center">2.953e-08</td>
+    <td>6.79</td>
+  </tr>
+  <tr>
+    <td align="right">16,384</td>
+    <td align="center">589,824</td>
+    <td align="center">1.137e-04</td>
+    <td>2.79</td>
+    <td align="center">1,048,576</td>
+    <td align="center">9.400e-07</td>
+    <td>4.08</td>
+    <td align="center">2,359,296</td>
+    <td align="center">4.286e-10</td>
+    <td>6.11</td>
+  </tr>
+  <tr>
+    <td align="right">65,536</td>
+    <td align="center">2,359,296</td>
+    <td align="center">1.476e-05</td>
+    <td>2.95</td>
+    <td align="center">4,194,304</td>
+    <td align="center">5.799e-08</td>
+    <td>4.02</td>
+    <td align="center">9,437,184</td>
+    <td align="center">2.789e-11</td>
+    <td>3.94</td>
+  </tr>
+  <tr>
+    <td align="right">262,144</td>
+    <td align="center">9,437,184</td>
+    <td align="center">2.038e-06</td>
+    <td>2.86</td>
+    <td align="center">16,777,216</td>
+    <td align="center">3.609e-09</td>
+    <td>4.01</td>
+    <td align="center">37,748,736</td>
+    <td align="center">5.730e-11</td>
+    <td>-1.04</td>
+  </tr>
+</table>
+
+The tables show that we get optimal $\mathcal O\left(h^{p+1}\right)$
+convergence rates for both numerical fluxes. The errors are slighly smaller
+for the Lax--Friedrichs flux for $p=2$, but the picture is reversed for
+$p=3$. For $p=5$, we reach the roundoff accuracy of $10^{-11}$ with both
+fluxes on the finest grids. Also note that the errors are absolute with a
+domain length of $10^2$, so relative errors are below $10^{-12}$. The HLL flux
+is somewhat better for the highest degree, which is due to a slight inaccuracy
+of the Lax--Friedrichs flux: The Lax--Friedrichs flux sets a Dirichlet
+condition on the solution that leaves the domain, which results in a small
+artificial reflection, which is accentuated for the Lax--Friedrichs
+flux. Apart from that, we see that the influence of the numerical flux is
+minor, as the polynomial part inside elements is the main driver of the
+accucary. The limited influence of the flux also has consequences when trying
+to approach more challenging setups with the higher-order DG setup: Taking for
+example the parameters and grid of step-33, we get oscillations (which in turn
+make density negative and make the solution explode) with both fluxes once the
+high-mass part comes near the boundary, as opposed to the low-order finite
+volume case ($p=0$). Thus, any case that leads to shocks in the solution
+necessitates some form of limiting or artificial dissipation. For another
+alternative, see the step-69 tutorial program.
+
+<h3>Results for flow in channel around cylinder in 2D</h3>
+
+For the test case of the flow around a cylinder in a channel, we need to
+change the first code line to `constexpr unsigned int testcase = 1;`.  This
+test case starts with a background field of a constant velocity of Ma=0.31 and
+density around an obstacle in the form of a cylinder. Since we impose a
+no-penetration condition on the cylinder walls, the flow has to rearrange,
+which creates a big sound wave. The following pictures show the pressure at
+times 0.1, 0.25, 0.5, and 1.0 (top left to bottom right) for the 2D case with
+5 levels of global refinement. We clearly see the discontinuity that
+propagates slowly in the upstream direction and more quickly in downstream
+direction in the first snapshot at time 0.1. At time 0.25, the sound wave has
+reached the top and bottom walls and reflected back to the interior. From the
+different position of the reflected waves from lower and upper walls we can
+see the slight asymmetry of the Sch&auml;fer-Turek test case represented by
+GridGenerator::channel_with_cylinder() with somewhat more space above the
+cylinder compared to below. At later times, the picture is more chaotic with
+many sound waves all over the place.
+
+<table align="center" class="doxtable" style="width:85%">
+  <tr>
+    <td>
+        <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_010.png" alt="" width="100%">
+    </td>
+    <td>
+        <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_025.png" alt="" width="100%">
+    </td>
+  </tr>
+  <tr>
+    <td>
+        <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_050.png" alt="" width="100%">
+    </td>
+    <td>
+        <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_100.png" alt="" width="100%">
+    </td>
+  </tr>
+</table>
+
+The next picture shows an elevation plot of the pressure at time 1.0 looking
+from the channel inlet towards the outlet -- here, we can see the large number
+of reflections. In the figure, two types of waves are visible. The
+larger-amplitude waves correspond to various reflections that happened as the
+initial discontinuity hit the walls, whereas the small-amplitude waves of the
+size similar to the elements correspond to numerical artifacts. They have the
+origin in the finite resolution of the scheme and appear as the discontinuity
+travels through elements with high-order polynomials. This effect can be cured
+by increasing resolution. Apart from this effect, the rich wave structure is
+the result of the transport accuracy of the high-order DG method.
+
+<img src="https://www.dealii.org/images/steps/developer/step-67.pressure_elevated.jpg" alt="" width="40%">
+
+With 2 levels of global refinement, the mesh and its partitioning on 40 MPI
+processes looks as follows:
+
+<img src="https://www.dealii.org/images/steps/developer/step-67.grid-owner.png" alt="" width="70%">
+
+When we run the code with 4 levels of global refinements on 40 cores, we get
+the following output:
+@code
+Running with 40 MPI processes
+Vectorization over 8 doubles = 512 bits (AVX512), VECTORIZATION_LEVEL=3
+Number of degrees of freedom: 3,686,400 ( = 4 [vars] x 25,600 [cells] x 36 [dofs/cell/var] )
+Time step size: 7.39876e-05, minimal h: 0.001875, initial transport scaling: 0.00110294
+
+Time:       0 , dt:  7.4e-05 , norm rho:   4.17e-16 , rho * u:  1.629e-16 , energy: 1.381e-15
+Time:    0.05 , dt:  6.3e-05 , norm rho:    0.02075 , rho * u:    0.03801 , energy:   0.08772
+Time:     0.1 , dt:  5.9e-05 , norm rho:    0.02211 , rho * u:    0.04515 , energy:   0.08953
+Time:    0.15 , dt:  5.7e-05 , norm rho:    0.02261 , rho * u:    0.04592 , energy:   0.08967
+Time:     0.2 , dt:  5.8e-05 , norm rho:    0.02058 , rho * u:    0.04361 , energy:   0.08222
+Time:    0.25 , dt:  5.9e-05 , norm rho:    0.01695 , rho * u:    0.04203 , energy:   0.06873
+Time:     0.3 , dt:  5.9e-05 , norm rho:    0.01653 , rho * u:     0.0401 , energy:   0.06604
+Time:    0.35 , dt:  5.7e-05 , norm rho:    0.01774 , rho * u:    0.04264 , energy:    0.0706
+
+...
+
+Time:    1.95 , dt:  5.8e-05 , norm rho:    0.01488 , rho * u:    0.03923 , energy:   0.05185
+Time:       2 , dt:  5.7e-05 , norm rho:    0.01432 , rho * u:    0.03969 , energy:   0.04889
+
++-------------------------------------------+------------------+------------+------------------+
+| Total wallclock time elapsed              |     273.6s    13 |     273.6s |     273.6s     0 |
+|                                           |                  |                               |
+| Section                       | no. calls |   min time  rank |   avg time |   max time  rank |
++-------------------------------------------+------------------+------------+------------------+
+| compute errors                |        41 |   0.01112s    35 |    0.0672s |    0.1337s     0 |
+| compute transport speed       |      6914 |     5.422s    35 |     15.96s |     29.99s     1 |
+| output                        |        41 |     37.24s    35 |      37.3s |     37.37s     0 |
+| rk time step total            |     34564 |     205.4s     1 |     219.5s |     230.1s    35 |
+| rk_stage - integrals L_h      |    172820 |     153.6s     1 |     164.9s |     175.6s    27 |
+| rk_stage - inv mass + vec upd |    172820 |     47.13s    13 |     53.09s |     64.05s    33 |
++-------------------------------------------+------------------+------------+------------------+
+@endcode
+
+We note that the norms we print for the various quantities are the deviations
+$\rho'$, $(\rho u)'$, and $E'$ against the background field which is the
+initial condition. The distribution of run time is overall similar as in the
+previous test case. The only slight difference is the larger proportion of
+time spent in L_h as compared to the inverse mass matrix and vector
+updates. This is because the geometry is deformed and the matrix-free
+framework needs to load additional arrays for the geometry from memory that
+are compressed in the affine mesh case.
+
+Increasing the number of global refinements to 5, the output becomes:
+@code
+Running with 40 MPI processes
+Vectorization over 8 doubles = 512 bits (AVX512), VECTORIZATION_LEVEL=3
+Number of degrees of freedom: 14,745,600 ( = 4 [vars] x 102,400 [cells] x 36 [dofs/cell/var] )
+
+...
+
++-------------------------------------------+------------------+------------+------------------+
+| Total wallclock time elapsed              |      2693s    32 |      2693s |      2693s    23 |
+|                                           |                  |                               |
+| Section                       | no. calls |   min time  rank |   avg time |   max time  rank |
++-------------------------------------------+------------------+------------+------------------+
+| compute errors                |        41 |   0.04537s    32 |     0.173s |    0.3489s     0 |
+| compute transport speed       |     13858 |     40.75s    32 |     85.99s |     149.8s     0 |
+| output                        |        41 |     153.8s    32 |     153.9s |     154.1s     0 |
+| rk time step total            |     69284 |      2386s     0 |      2450s |      2496s    32 |
+| rk_stage - integrals L_h      |    346420 |      1365s    32 |      1574s |      1718s    19 |
+| rk_stage - inv mass + vec upd |    346420 |     722.5s    10 |     870.7s |      1125s    32 |
++-------------------------------------------+------------------+------------+------------------+
+@endcode
+
+The effect on performance is similar as for the analytical test case -- in
+theory, computation times should increase by a factor of 8, but we actually
+see an increase by a factor of 11 for the time steps (219.5 seconds versus
+2450 seconds). This can be traced back to caches, with the small case mostly
+fitting in caches. An interesting effect, typical of programs with a mix of
+local communication (integrals L_h) and global communication (computation of
+transport speed) with some load imbalance, can be observed by looking at the
+MPI rank that measure the minimal and maximal time of different phases,
+respectively. Rank 0 reports the fastest throughput for the "rk time step
+total" part. At the same time, it appears to be slowest for the "compute
+transport speed" part, more almost a factor of 2 slower than the
+average. Since the latter involves global communication, we can attribute the
+slowness in this part to the fact that the local Runge--Kutta stages have
+advanced more quickly on this rank and need to wait until the other processors
+catch up. At this point, one can wonder about the reason for this imbalance:
+The number of cells is almost the same on all MPI processes due to the default
+weights. However, the matrix-free framework is faster on affine and Cartesian
+cells located towards the outlet of the channel, to which the lower MPI ranks
+are assigned. On the other hand, rank 32, which reports the highest run time
+for the Runga--Kutta stages, owns the curved cells near the cylinder, for
+which no data compression is possible. To improve throughput, we could assign
+different weights to different cell types, or even measure the run time for a
+few time steps and try to rebalance then.
+
+The throughput per Runge--Kutta stage can be computed to 2085 MDoFs/s for the
+14.7 million DoFs test case over 346k Runge--Kutta stages, slightly slower
+than the Cartesian mesh throughput of 2360 MDoFs/s reported above.
+
+Finally, if we add one additional refinement, we record the following output:
+@code
+Running with 40 MPI processes
+Vectorization over 8 doubles = 512 bits (AVX512), VECTORIZATION_LEVEL=3
+Number of degrees of freedom: 58,982,400 ( = 4 [vars] x 409,600 [cells] x 36 [dofs/cell/var] )
+
+...
+
+Time:    1.95 , dt:  1.4e-05 , norm rho:    0.01488 , rho * u:    0.03923 , energy:   0.05183
+Time:       2 , dt:  1.4e-05 , norm rho:    0.01431 , rho * u:    0.03969 , energy:   0.04887
+
++-------------------------------------------+------------------+------------+------------------+
+| Total wallclock time elapsed              | 2.166e+04s    26 | 2.166e+04s | 2.166e+04s    24 |
+|                                           |                  |                               |
+| Section                       | no. calls |   min time  rank |   avg time |   max time  rank |
++-------------------------------------------+------------------+------------+------------------+
+| compute errors                |        41 |    0.1758s    30 |     0.672s |     1.376s     1 |
+| compute transport speed       |     27748 |     321.3s    34 |     678.8s |      1202s     1 |
+| output                        |        41 |     616.3s    32 |     616.4s |     616.4s    34 |
+| rk time step total            |    138733 | 1.983e+04s     1 | 2.036e+04s | 2.072e+04s    34 |
+| rk_stage - integrals L_h      |    693665 | 1.052e+04s    32 | 1.248e+04s | 1.387e+04s    19 |
+| rk_stage - inv mass + vec upd |    693665 |      6404s    10 |      7868s | 1.018e+04s    32 |
++-------------------------------------------+------------------+------------+------------------+
+@endcode
+
+The "rk time stop total" part corresponds to a throughput of 2010 MDoFs/s. The
+overall run time to perform 139k time steps is 20k seconds (5.7 hours) or 7
+time steps per second. More throughput can be achieved by adding more cores to
+the computation.
+
+<h3>Results for flow in channel around cylinder in 3D</h3>
+
+Switching the channel test case to 3D with 3 global refinements, the output is
+@code
+Running with 40 MPI processes
+Vectorization over 8 doubles = 512 bits (AVX512), VECTORIZATION_LEVEL=3
+Number of degrees of freedom: 221,184,000 ( = 5 [vars] x 204,800 [cells] x 216 [dofs/cell/var] )
+
+...
+
+Time:    1.95 , dt:  0.00011 , norm rho:    0.01131 , rho * u:    0.03056 , energy:   0.04091
+Time:       2 , dt:  0.00011 , norm rho:     0.0119 , rho * u:    0.03142 , energy:   0.04425
+
++-------------------------------------------+------------------+------------+------------------+
+| Total wallclock time elapsed              | 1.734e+04s     4 | 1.734e+04s | 1.734e+04s    38 |
+|                                           |                  |                               |
+| Section                       | no. calls |   min time  rank |   avg time |   max time  rank |
++-------------------------------------------+------------------+------------+------------------+
+| compute errors                |        41 |    0.6551s    34 |     3.216s |     7.281s     0 |
+| compute transport speed       |      3546 |       160s    34 |     393.2s |     776.9s     0 |
+| output                        |        41 |      1350s    34 |      1353s |      1357s     0 |
+| rk time step total            |     17723 | 1.519e+04s     0 | 1.558e+04s | 1.582e+04s    34 |
+| rk_stage - integrals L_h      |     88615 | 1.005e+04s    32 | 1.126e+04s |  1.23e+04s    11 |
+| rk_stage - inv mass + vec upd |     88615 |      3056s    11 |      4322s |      5759s    32 |
++-------------------------------------------+------------------+------------+------------------+
+@endcode
+
+The physics are similar to the 2D case, with a slight motion in the z
+direction due to the gravitational force. The throughput per Runge--Kutta
+stage in this case is
+@f[
+\text{throughput} = \frac{n_\mathrm{time steps} n_\mathrm{stages}
+n_\mathrm{dofs}}{t_\mathrm{compute}} =
+\frac{17723 \cdot 5 \cdot 221.2\,\text{M}}{15580s} = 1258\, \text{MDoFs/s}.
+@f]
+
+The throughput is lower than in 2D because the computation of the $L_h$ term
+is more expensive. This is due to over-integration with `degree+2` points and
+the larger fraction of face integrals (worse volume-to-surface ratio) with
+more expensive flux computations. If we only consider the inverse mass matrix
+and vector update part, we record a throughput of 4857 MDoFs/s for the 2D case
+of the isentropic vortex with 37.7 million unknowns, whereas the 3D case
+runs with 4535 MDoFs/s. The performance is similar because both cases are in
+fact limited by the memory bandwidth.
+
+If we go to four levels of global refinement, we need to increase the number
+of processes to fit everything in memory -- the computation needs around 350
+GB of RAM memory in this case. Also, the time it takes to complete 35k time
+steps becomes more tolerable by adding additional resources. We therefore use
+6 nodes with 40 cores each, resulting in a computation with 240 MPI processes:
+@code
+Running with 240 MPI processes
+Vectorization over 8 doubles = 512 bits (AVX512), VECTORIZATION_LEVEL=3
+Number of degrees of freedom: 1,769,472,000 ( = 5 [vars] x 1,638,400 [cells] x 216 [dofs/cell/var] )
+
+...
+
+Time:    1.95 , dt:  5.6e-05 , norm rho:    0.01129 , rho * u:     0.0306 , energy:   0.04086
+Time:       2 , dt:  5.6e-05 , norm rho:    0.01189 , rho * u:    0.03145 , energy:   0.04417
+
++-------------------------------------------+------------------+------------+------------------+
+| Total wallclock time elapsed              | 5.396e+04s   151 | 5.396e+04s | 5.396e+04s     0 |
+|                                           |                  |                               |
+| Section                       | no. calls |   min time  rank |   avg time |   max time  rank |
++-------------------------------------------+------------------+------------+------------------+
+| compute errors                |        41 |     2.632s   178 |     7.221s |     16.56s     0 |
+| compute transport speed       |      7072 |       714s   193 |      1553s |      3351s     0 |
+| output                        |     35350 |      8065s   176 |      8070s |      8079s     0 |
+| rk time step total            |     35350 |  4.25e+04s     0 |  4.43e+04s | 4.515e+04s   193 |
+| rk_stage - integrals L_h      |    176750 | 2.936e+04s   134 | 3.222e+04s |  3.67e+04s    99 |
+| rk_stage - inv mass + vec upd |    176750 |      7004s    99 | 1.207e+04s |  1.55e+04s   132 |
++-------------------------------------------+------------------+------------+------------------+
+@endcode
+
+<h3>Possibilities for extensions</h3>
+
+The code presented here straight-forwardly extends to adaptive meshes, given
+appropriate indicators for setting the refinement flags. Large-scale
+adaptivity of a similar solver in the context of the acoustic wave equation
+has been derived by the <a href="https://github.com/kronbichler/exwave">exwave
+project</a>. However, in the present context the effect of adaptivity is often
+limited to early times and effects close to the origin of sound waves, as the
+waves eventually reflect and diffract. This leads to steep gradients all over
+the place, similar to turbulent flow.
+
+Another topic that we did not discuss in the results section is a comparison
+of different time integration schemes. The program provides four variants of
+low-storage Runga--Kutta integrators that each have slightly different
+accuracy and stability behavior. Among the schemes implemented here, the
+higher-order ones provide additional accuracy but come with slightly lower
+efficiency in terms of step size per stage before they violate the CFL
+condition. An interesting extension would be to compare the low-storage
+variants proposed here with standard Runge--Kutta integrators or to use vector
+operations that are run separate from the mass matrix operation and compare
+performance.
+
+<h4>More advanced numerical flux functions and skew-symmetric formulations</h4>
+
+As mentioned in the introduction, the modified Lax--Friedrichs flux and the
+HLL flux employed in this program are only two variants of a large body of
+numerical fluxes available in the literature on the Euler equations. One
+example is the HLLC flux (Harten-Lax-van Leer-Contact) flux which adds the
+effect of rarefaction waves missing in the HLL flux, or the Roe flux. As
+mentioned in the introduction, the effect of numerical fluxes on high-order DG
+schemes is debatable.
+
+A related improvement to increase the stability of the solver is to also
+consider the spatial integral terms. A shortcoming in the rather naive
+implementation used above is the fact that the energy conservation of the
+original Euler equations (in the absence of shocks) only holds up to a
+discretization error. If the solution is under-resolved, the discretization
+error can give rise to an increase in the numerical energy and eventually
+render the discretization unstable. This is because of the inexact numerical
+integration of the terms in the Euler equations, which both contain rational
+nonlinearities and higher-degree content from curved cells. A way out of this
+dilemma are so-called skew-symmetric formulations, see e.g. <a
+href=https://doi.org/10.1137%2F120890144>Gassner (2013)</a> for a simple
+variant. Skew symmetry means that switching the role of the solution
+$\mathbf{w}$ and test functions $\mathbf{v}$ in the weak form produces the
+exact negative of the original quantity, apart from some boundary terms. In
+the discrete setting, the challenge is to keep this skew symmetry also when
+the integrals are only computed approximately (in the continuous case,
+skew-symmetry is a consequence of integration by parts). Skew-symmetric
+numerical schemes balance spatial derivatives in the conservative form
+$(\nabla \mathbf v, \mathbf{F}(\mathbf w))_{K}$ with contributions in the
+convective form $(\mathbf v, \tilde{\mathbf{F}}(\mathbf w)\nabla
+\mathbf{w})_{K}$ for some $\tilde{\mathbf{F}}$. The precise terms depend on
+the equation and the integration formula, and can in some cases by understood
+by special skew-symmetric finite difference schemes.
+
+<h4>Equipping the code for supersonic calculations</h4>
+
+As mentioned in the introduction, the solution to the Euler equations develops
+shocks as the Mach number increases, which require additional mechanisms to
+stabilize the scheme, e.g. in the form of limiters. The main challenge besides
+actually implementing the limiter or artificial viscosity approach would be to
+load-balance the computations, as the additional computations involved for
+limiting the oscillations in troubled cells would make them heavier than the
+plain DG cells without limiting. Furthermore, additional numerical fluxes that
+better cope with the discontinuities would also be an option.
+
+One ingredient also necessary for supersonic flows are appropriate boundary
+conditions. As opposed to the subsonic outflow boundaries discussed in the
+introduction and implemented in the program, all characteristics are outgoing
+for supersonic outflow boundaries, so we do not want to prescribe any external
+data,
+@f[
+\mathbf{w}^+ = \mathbf{w}^- = \begin{pmatrix} \rho^-\\
+(\rho \mathbf u)^- \\ E^-\end{pmatrix} \quad
+ \text{(Neumann)}.
+@f]
+
+In the code, we would simply add the additional statement
+@code
+            else if (supersonic_outflow_boundaries.find(boundary_id) !=
+                     supersonic_outflow_boundaries.end())
+              {
+                w_p        = w_m;
+                at_outflow = true;
+              }
+@endcode
+in the `local_apply_boundary_face()` function.
+
+<h4>Extension to the linearized Euler equations</h4>
+
+When the interest with an Euler solution is mostly in the propagation of sound
+waves, it often makes sense to linearize the Euler equations around a
+background state, i.e., a given density, velocity and energy (or pressure)
+field, and only compute the change against these fields. This is the setting
+of the wide field of aeroacoustics. Even though the resolution requirements
+are sometimes considerably reduced, implementation gets somewhat more
+complicated as the linearization gives rise to additional terms. From a code
+perspective, in the operator evaluation we also need to equip the code with
+the state to linearize against. This information can be provided either by
+analytical functions (that are evaluated in terms of the position of the
+quadrature points) or by a vector similar to the solution. Based on that
+vector, we would create an additional FEEvaluation object to read from it and
+provide the values of the field at quadrature points. If the background
+velocity is zero and the density is constant, the linearized Euler equations
+further simplify to the acoustic wave equation.
+
+A challenge in the context of sound propagation is often the definition of
+boundary conditions, as the computational domain needs to be of finite size,
+whereas the actual simulation often spans an infinite (or at least much
+larger) physical domain. Conventional Dirichlet or Neumann boundary conditions
+give rise to reflections of the sound waves that eventually propagate back to
+the region of interest and spoil the solution. Therefore, various variants of
+non-reflecting boundary conditions or sponge layers, often in the form of
+perfectly matched layers -- where the solution is damped without reflection --
+are common.
+
+<h4>Extension to the compressible Navier-Stokes equation</h4>
+
+The solver presented in this tutorial program can also be extended to the
+compressible Navier--Stokes equations by adding viscous terms, as described in
+@cite FehnWallKronbichler2019. To keep as much of the performance obtained
+here despite the additional cost of elliptic terms e.g. via an interior
+penalty method, one can switch the basis from FE_DGQ to FE_DGQHermite like in
+the step-59 tutorial program.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.