Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
step-67.h
Go to the documentation of this file.
1 
898  * stage_5_order_4, /* Kennedy, Carpenter, Lewis, 2000 */
899  * stage_7_order_4, /* Tselios, Simos, 2007 */
900  * stage_9_order_5, /* Kennedy, Carpenter, Lewis, 2000 */
901  * };
902  * constexpr LowStorageRungeKuttaScheme lsrk_scheme = stage_5_order_4;
903  *
904  * @endcode
905  *
906  * Eventually, we select a detail of the spatial discretization, namely the
907  * numerical flux (Riemann solver) at the faces between cells. For this
908  * program, we have implemented a modified variant of the Lax--Friedrichs
909  * flux and the Harten--Lax--van Leer (HLL) flux.
910  *
911  * @code
912  * enum EulerNumericalFlux
913  * {
914  * lax_friedrichs_modified,
915  * harten_lax_vanleer,
916  * };
917  * constexpr EulerNumericalFlux numerical_flux_type = lax_friedrichs_modified;
918  *
919  *
920  *
921  * @endcode
922  *
923  *
924  * <a name="Equationdata"></a>
925  * <h3>Equation data</h3>
926  *
927 
928  *
929  * We now define a class with the exact solution for the test case 0 and one
930  * with a background flow field for test case 1 of the channel. Given that
931  * the Euler equations are a problem with @f$d+2@f$ equations in @f$d@f$ dimensions,
932  * we need to tell the Function base class about the correct number of
933  * components.
934  *
935  * @code
936  * template <int dim>
937  * class ExactSolution : public Function<dim>
938  * {
939  * public:
940  * ExactSolution(const double time)
941  * : Function<dim>(dim + 2, time)
942  * {}
943  *
944  * virtual double value(const Point<dim> & p,
945  * const unsigned int component = 0) const override;
946  * };
947  *
948  *
949  *
950  * @endcode
951  *
952  * As far as the actual function implemented is concerned, the analytical
953  * test case is an isentropic vortex case (see e.g. the book by Hesthaven
954  * and Warburton, Example 6.1 in Section 6.6 on page 209) which fulfills the
955  * Euler equations with zero force term on the right hand side. Given that
956  * definition, we return either the density, the momentum, or the energy
957  * depending on which component is requested. Note that the original
958  * definition of the density involves the @f$\frac{1}{\gamma -1}@f$-th power of
959  * some expression. Since `std::pow()` has pretty slow implementations on
960  * some systems, we replace it by logarithm followed by exponentiation (of
961  * base 2), which is mathematically equivalent but usually much better
962  * optimized. This formula might lose accuracy in the last digits
963  * for very small numbers compared to `std::pow()`, but we are happy with
964  * it anyway, since small numbers map to data close to 1.
965  *
966 
967  *
968  * For the channel test case, we simply select a density of 1, a velocity of
969  * 0.4 in @f$x@f$ direction and zero in the other directions, and an energy that
970  * corresponds to a speed of sound of 1.3 measured against the background
971  * velocity field, computed from the relation @f$E = \frac{c^2}{\gamma (\gamma
972  * -1)} + \frac 12 \rho \|u\|^2@f$.
973  *
974  * @code
975  * template <int dim>
976  * double ExactSolution<dim>::value(const Point<dim> & x,
977  * const unsigned int component) const
978  * {
979  * const double t = this->get_time();
980  *
981  * switch (testcase)
982  * {
983  * case 0:
984  * {
985  * Assert(dim == 2, ExcNotImplemented());
986  * const double beta = 5;
987  *
988  * Point<dim> x0;
989  * x0[0] = 5.;
990  * const double radius_sqr =
991  * (x - x0).norm_square() - 2. * (x[0] - x0[0]) * t + t * t;
992  * const double factor =
993  * beta / (numbers::PI * 2) * std::exp(1. - radius_sqr);
994  * const double density_log = std::log2(
995  * std::abs(1. - (gamma - 1.) / gamma * 0.25 * factor * factor));
996  * const double density = std::exp2(density_log * (1. / (gamma - 1.)));
997  * const double u = 1. - factor * (x[1] - x0[1]);
998  * const double v = factor * (x[0] - t - x0[0]);
999  *
1000  * if (component == 0)
1001  * return density;
1002  * else if (component == 1)
1003  * return density * u;
1004  * else if (component == 2)
1005  * return density * v;
1006  * else
1007  * {
1008  * const double pressure =
1009  * std::exp2(density_log * (gamma / (gamma - 1.)));
1010  * return pressure / (gamma - 1.) +
1011  * 0.5 * (density * u * u + density * v * v);
1012  * }
1013  * }
1014  *
1015  * case 1:
1016  * {
1017  * if (component == 0)
1018  * return 1.;
1019  * else if (component == 1)
1020  * return 0.4;
1021  * else if (component == dim + 1)
1022  * return 3.097857142857143;
1023  * else
1024  * return 0.;
1025  * }
1026  *
1027  * default:
1028  * Assert(false, ExcNotImplemented());
1029  * return 0.;
1030  * }
1031  * }
1032  *
1033  *
1034  *
1035  * @endcode
1036  *
1037  *
1038  * <a name="LowstorageexplicitRungeKuttatimeintegrators"></a>
1039  * <h3>Low-storage explicit Runge--Kutta time integrators</h3>
1040  *
1041 
1042  *
1043  * The next few lines implement a few low-storage variants of Runge--Kutta
1044  * methods. These methods have specific Butcher tableaux with coefficients
1045  * @f$b_i@f$ and @f$a_i@f$ as shown in the introduction. As usual in Runge--Kutta
1046  * method, we can deduce time steps, @f$c_i = \sum_{j=1}^{i-2} b_i + a_{i-1}@f$
1047  * from those coefficients. The main advantage of this kind of scheme is the
1048  * fact that only two vectors are needed per stage, namely the accumulated
1049  * part of the solution @f$\mathbf{w}@f$ (that will hold the solution
1050  * @f$\mathbf{w}^{n+1}@f$ at the new time @f$t^{n+1}@f$ after the last stage), the
1051  * update vector @f$\mathbf{r}_i@f$ that gets evaluated during the stages, plus
1052  * one vector @f$\mathbf{k}_i@f$ to hold the evaluation of the operator. Such a
1053  * Runge--Kutta setup reduces the memory storage and memory access. As the
1054  * memory bandwidth is often the performance-limiting factor on modern
1055  * hardware when the evaluation of the differential operator is
1056  * well-optimized, performance can be improved over standard time
1057  * integrators. This is true also when taking into account that a
1058  * conventional Runge--Kutta scheme might allow for slightly larger time
1059  * steps as more free parameters allow for better stability properties.
1060  *
1061 
1062  *
1063  * In this tutorial programs, we concentrate on a few variants of
1064  * low-storage schemes defined in the article by Kennedy, Carpenter, and
1065  * Lewis (2000), as well as one variant described by Tselios and Simos
1066  * (2007). There is a large series of other schemes available, which could
1067  * be addressed by additional sets of coefficients or slightly different
1068  * update formulas.
1069  *
1070 
1071  *
1072  * We define a single class for the four integrators, distinguished by the
1073  * enum described above. To each scheme, we then fill the vectors for the
1074  * @f$b_i@f$ and @f$a_i@f$ to the given variables in the class.
1075  *
1076  * @code
1077  * class LowStorageRungeKuttaIntegrator
1078  * {
1079  * public:
1080  * LowStorageRungeKuttaIntegrator(const LowStorageRungeKuttaScheme scheme)
1081  * {
1082  * @endcode
1083  *
1084  * First comes the three-stage scheme of order three by Kennedy et al.
1085  * (2000). While its stability region is significantly smaller than for
1086  * the other schemes, it only involves three stages, so it is very
1087  * competitive in terms of the work per stage.
1088  *
1089  * @code
1090  * switch (scheme)
1091  * {
1092  * case stage_3_order_3:
1093  * {
1094  * bi = {{0.245170287303492, 0.184896052186740, 0.569933660509768}};
1095  * ai = {{0.755726351946097, 0.386954477304099}};
1096  *
1097  * break;
1098  * }
1099  *
1100  * @endcode
1101  *
1102  * The next scheme is a five-stage scheme of order four, again
1103  * defined in the paper by Kennedy et al. (2000).
1104  *
1105  * @code
1106  * case stage_5_order_4:
1107  * {
1108  * bi = {{1153189308089. / 22510343858157.,
1109  * 1772645290293. / 4653164025191.,
1110  * -1672844663538. / 4480602732383.,
1111  * 2114624349019. / 3568978502595.,
1112  * 5198255086312. / 14908931495163.}};
1113  * ai = {{970286171893. / 4311952581923.,
1114  * 6584761158862. / 12103376702013.,
1115  * 2251764453980. / 15575788980749.,
1116  * 26877169314380. / 34165994151039.}};
1117  *
1118  * break;
1119  * }
1120  *
1121  * @endcode
1122  *
1123  * The following scheme of seven stages and order four has been
1124  * explicitly derived for acoustics problems. It is a balance of
1125  * accuracy for imaginary eigenvalues among fourth order schemes,
1126  * combined with a large stability region. Since DG schemes are
1127  * dissipative among the highest frequencies, this does not
1128  * necessarily translate to the highest possible time step per
1129  * stage. In the context of the present tutorial program, the
1130  * numerical flux plays a crucial role in the dissipation and thus
1131  * also the maximal stable time step size. For the modified
1132  * Lax--Friedrichs flux, this scheme is similar to the
1133  * `stage_5_order_4` scheme in terms of step size per stage if only
1134  * stability is considered, but somewhat less efficient for the HLL
1135  * flux.
1136  *
1137  * @code
1138  * case stage_7_order_4:
1139  * {
1140  * bi = {{0.0941840925477795334,
1141  * 0.149683694803496998,
1142  * 0.285204742060440058,
1143  * -0.122201846148053668,
1144  * 0.0605151571191401122,
1145  * 0.345986987898399296,
1146  * 0.186627171718797670}};
1147  * ai = {{0.241566650129646868 + bi[0],
1148  * 0.0423866513027719953 + bi[1],
1149  * 0.215602732678803776 + bi[2],
1150  * 0.232328007537583987 + bi[3],
1151  * 0.256223412574146438 + bi[4],
1152  * 0.0978694102142697230 + bi[5]}};
1153  *
1154  * break;
1155  * }
1156  *
1157  * @endcode
1158  *
1159  * The last scheme included here is the nine-stage scheme of order
1160  * five from Kennedy et al. (2000). It is the most accurate among
1161  * the schemes used here, but the higher order of accuracy
1162  * sacrifices some stability, so the step length normalized per
1163  * stage is less than for the fourth order schemes.
1164  *
1165  * @code
1166  * case stage_9_order_5:
1167  * {
1168  * bi = {{2274579626619. / 23610510767302.,
1169  * 693987741272. / 12394497460941.,
1170  * -347131529483. / 15096185902911.,
1171  * 1144057200723. / 32081666971178.,
1172  * 1562491064753. / 11797114684756.,
1173  * 13113619727965. / 44346030145118.,
1174  * 393957816125. / 7825732611452.,
1175  * 720647959663. / 6565743875477.,
1176  * 3559252274877. / 14424734981077.}};
1177  * ai = {{1107026461565. / 5417078080134.,
1178  * 38141181049399. / 41724347789894.,
1179  * 493273079041. / 11940823631197.,
1180  * 1851571280403. / 6147804934346.,
1181  * 11782306865191. / 62590030070788.,
1182  * 9452544825720. / 13648368537481.,
1183  * 4435885630781. / 26285702406235.,
1184  * 2357909744247. / 11371140753790.}};
1185  *
1186  * break;
1187  * }
1188  *
1189  * default:
1190  * AssertThrow(false, ExcNotImplemented());
1191  * }
1192  * }
1193  *
1194  * unsigned int n_stages() const
1195  * {
1196  * return bi.size();
1197  * }
1198  *
1199  * @endcode
1200  *
1201  * The main function of the time integrator is to go through the stages,
1202  * evaluate the operator, prepare the @f$\mathbf{r}_i@f$ vector for the next
1203  * evaluation, and update the solution vector @f$\mathbf{w}@f$. We hand off
1204  * the work to the `pde_operator` involved in order to be able to merge
1205  * the vector operations of the Runge--Kutta setup with the evaluation of
1206  * the differential operator for better performance, so all we do here is
1207  * to delegate the vectors and coefficients.
1208  *
1209 
1210  *
1211  * We separately call the operator for the first stage because we need
1212  * slightly modified arguments there: We evaluate the solution from
1213  * the old solution @f$\mathbf{w}^n@f$ rather than a @f$\mathbf r_i@f$ vector, so
1214  * the first argument is `solution`. We here let the stage vector
1215  * @f$\mathbf{r}_i@f$ also hold the temporary result of the evaluation, as it
1216  * is not used otherwise. For all subsequent stages, we use the vector
1217  * `vec_ki` as the second vector argument to store the result of the
1218  * operator evaluation. Finally, when we are at the last stage, we must
1219  * skip the computation of the vector @f$\mathbf{r}_{s+1}@f$ as there is no
1220  * coefficient @f$a_s@f$ available (nor will it be used).
1221  *
1222  * @code
1223  * template <typename VectorType, typename Operator>
1224  * void perform_time_step(const Operator &pde_operator,
1225  * const double current_time,
1226  * const double time_step,
1227  * VectorType & solution,
1228  * VectorType & vec_ri,
1229  * VectorType & vec_ki) const
1230  * {
1231  * AssertDimension(ai.size() + 1, bi.size());
1232  *
1233  * pde_operator.perform_stage(current_time,
1234  * bi[0] * time_step,
1235  * ai[0] * time_step,
1236  * solution,
1237  * vec_ri,
1238  * solution,
1239  * vec_ri);
1240  * double sum_previous_bi = 0;
1241  * for (unsigned int stage = 1; stage < bi.size(); ++stage)
1242  * {
1243  * const double c_i = sum_previous_bi + ai[stage - 1];
1244  * pde_operator.perform_stage(current_time + c_i * time_step,
1245  * bi[stage] * time_step,
1246  * (stage == bi.size() - 1 ?
1247  * 0 :
1248  * ai[stage] * time_step),
1249  * vec_ri,
1250  * vec_ki,
1251  * solution,
1252  * vec_ri);
1253  * sum_previous_bi += bi[stage - 1];
1254  * }
1255  * }
1256  *
1257  * private:
1258  * std::vector<double> bi;
1259  * std::vector<double> ai;
1260  * };
1261  *
1262  *
1263  *
1264  * @endcode
1265  *
1266  *
1267  * <a name="ImplementationofpointwiseoperationsoftheEulerequations"></a>
1268  * <h3>Implementation of point-wise operations of the Euler equations</h3>
1269  *
1270 
1271  *
1272  * In the following functions, we implement the various problem-specific
1273  * operators pertaining to the Euler equations. Each function acts on the
1274  * vector of conserved variables @f$[\rho, \rho\mathbf{u}, E]@f$ that we hold in
1275  * the solution vectors, and computes various derived quantities.
1276  *
1277 
1278  *
1279  * First out is the computation of the velocity, that we derive from the
1280  * momentum variable @f$\rho \mathbf{u}@f$ by division by @f$\rho@f$. One thing to
1281  * note here is that we decorate all those functions with the keyword
1282  * `DEAL_II_ALWAYS_INLINE`. This is a special macro that maps to a
1283  * compiler-specific keyword that tells the compiler to never create a
1284  * function call for any of those functions, and instead move the
1285  * implementation <a
1286  * href="https://en.wikipedia.org/wiki/Inline_function">inline</a> to where
1287  * they are called. This is critical for performance because we call into some
1288  * of those functions millions or billions of times: For example, we both use
1289  * the velocity for the computation of the flux further down, but also for the
1290  * computation of the pressure, and both of these places are evaluated at
1291  * every quadrature point of every cell. Making sure these functions are
1292  * inlined ensures not only that the processor does not have to execute a jump
1293  * instruction into the function (and the corresponding return jump), but also
1294  * that the compiler can re-use intermediate information from one function's
1295  * context in code that comes after the place where the function was called.
1296  * (We note that compilers are generally quite good at figuring out which
1297  * functions to inline by themselves. Here is a place where compilers may or
1298  * may not have figured it out by themselves but where we know for sure that
1299  * inlining is a win.)
1300  *
1301 
1302  *
1303  * Another trick we apply is a separate variable for the inverse density
1304  * @f$\frac{1}{\rho}@f$. This enables the compiler to only perform a single
1305  * division for the flux, despite the division being used at several
1306  * places. As divisions are around ten to twenty times as expensive as
1307  * multiplications or additions, avoiding redundant divisions is crucial for
1308  * performance. We note that taking the inverse first and later multiplying
1309  * with it is not equivalent to a division in floating point arithmetic due
1310  * to roundoff effects, so the compiler is not allowed to exchange one way by
1311  * the other with standard optimization flags. However, it is also not
1312  * particularly difficult to write the code in the right way.
1313  *
1314 
1315  *
1316  * To summarize, the chosen strategy of always inlining and careful
1317  * definition of expensive arithmetic operations allows us to write compact
1318  * code without passing all intermediate results around, despite making sure
1319  * that the code maps to excellent machine code.
1320  *
1321  * @code
1322  * template <int dim, typename Number>
1323  * inline DEAL_II_ALWAYS_INLINE
1324  * Tensor<1, dim, Number>
1325  * euler_velocity(const Tensor<1, dim + 2, Number> &conserved_variables)
1326  * {
1327  * const Number inverse_density = Number(1.) / conserved_variables[0];
1328  *
1329  * Tensor<1, dim, Number> velocity;
1330  * for (unsigned int d = 0; d < dim; ++d)
1331  * velocity[d] = conserved_variables[1 + d] * inverse_density;
1332  *
1333  * return velocity;
1334  * }
1335  *
1336  * @endcode
1337  *
1338  * The next function computes the pressure from the vector of conserved
1339  * variables, using the formula @f$p = (\gamma - 1) \left(E - \frac 12 \rho
1340  * \mathbf{u}\cdot \mathbf{u}\right)@f$. As explained above, we use the
1341  * velocity from the `euler_velocity()` function. Note that we need to
1342  * specify the first template argument `dim` here because the compiler is
1343  * not able to deduce it from the arguments of the tensor, whereas the
1344  * second argument (number type) can be automatically deduced.
1345  *
1346  * @code
1347  * template <int dim, typename Number>
1348  * inline DEAL_II_ALWAYS_INLINE
1349  * Number
1350  * euler_pressure(const Tensor<1, dim + 2, Number> &conserved_variables)
1351  * {
1352  * const Tensor<1, dim, Number> velocity =
1353  * euler_velocity<dim>(conserved_variables);
1354  *
1355  * Number rho_u_dot_u = conserved_variables[1] * velocity[0];
1356  * for (unsigned int d = 1; d < dim; ++d)
1357  * rho_u_dot_u += conserved_variables[1 + d] * velocity[d];
1358  *
1359  * return (gamma - 1.) * (conserved_variables[dim + 1] - 0.5 * rho_u_dot_u);
1360  * }
1361  *
1362  * @endcode
1363  *
1364  * Here is the definition of the Euler flux function, i.e., the definition
1365  * of the actual equation. Given the velocity and pressure (that the
1366  * compiler optimization will make sure are done only once), this is
1367  * straight-forward given the equation stated in the introduction.
1368  *
1369  * @code
1370  * template <int dim, typename Number>
1371  * inline DEAL_II_ALWAYS_INLINE
1372  * Tensor<1, dim + 2, Tensor<1, dim, Number>>
1373  * euler_flux(const Tensor<1, dim + 2, Number> &conserved_variables)
1374  * {
1375  * const Tensor<1, dim, Number> velocity =
1376  * euler_velocity<dim>(conserved_variables);
1377  * const Number pressure = euler_pressure<dim>(conserved_variables);
1378  *
1379  * Tensor<1, dim + 2, Tensor<1, dim, Number>> flux;
1380  * for (unsigned int d = 0; d < dim; ++d)
1381  * {
1382  * flux[0][d] = conserved_variables[1 + d];
1383  * for (unsigned int e = 0; e < dim; ++e)
1384  * flux[e + 1][d] = conserved_variables[e + 1] * velocity[d];
1385  * flux[d + 1][d] += pressure;
1386  * flux[dim + 1][d] =
1387  * velocity[d] * (conserved_variables[dim + 1] + pressure);
1388  * }
1389  *
1390  * return flux;
1391  * }
1392  *
1393  * @endcode
1394  *
1395  * This next function is a helper to simplify the implementation of the
1396  * numerical flux, implementing the action of a tensor of tensors (with
1397  * non-standard outer dimension of size `dim + 2`, so the standard overloads
1398  * provided by deal.II's tensor classes do not apply here) with another
1399  * tensor of the same inner dimension, i.e., a matrix-vector product.
1400  *
1401  * @code
1402  * template <int n_components, int dim, typename Number>
1403  * inline DEAL_II_ALWAYS_INLINE
1406  * const Tensor<1, dim, Number> & vector)
1407  * {
1409  * for (unsigned int d = 0; d < n_components; ++d)
1410  * result[d] = matrix[d] * vector;
1411  * return result;
1412  * }
1413  *
1414  * @endcode
1415  *
1416  * This function implements the numerical flux (Riemann solver). It gets the
1417  * state from the two sides of an interface and the normal vector, oriented
1418  * from the side of the solution @f$\mathbf{w}^-@f$ towards the solution
1419  * @f$\mathbf{w}^+@f$. In finite volume methods which rely on piece-wise
1420  * constant data, the numerical flux is the central ingredient as it is the
1421  * only place where the physical information is entered. In DG methods, the
1422  * numerical flux is less central due to the polynomials within the elements
1423  * and the physical flux used there. As a result of higher-degree
1424  * interpolation with consistent values from both sides in the limit of a
1425  * continuous solution, the numerical flux can be seen as a control of the
1426  * jump of the solution from both sides to weakly impose continuity. It is
1427  * important to realize that a numerical flux alone cannot stabilize a
1428  * high-order DG method in the presence of shocks, and thus any DG method
1429  * must be combined with further shock-capturing techniques to handle those
1430  * cases. In this tutorial, we focus on wave-like solutions of the Euler
1431  * equations in the subsonic regime without strong discontinuities where our
1432  * basic scheme is sufficient.
1433  *
1434 
1435  *
1436  * Nonetheless, the numerical flux is decisive in terms of the numerical
1437  * dissipation of the overall scheme and influences the admissible time step
1438  * size with explicit Runge--Kutta methods. We consider two choices, a
1439  * modified Lax--Friedrichs scheme and the widely used Harten--Lax--van Leer
1440  * (HLL) flux. For both variants, we first need to get the velocities and
1441  * pressures from both sides of the interface and evaluate the physical
1442  * Euler flux.
1443  *
1444 
1445  *
1446  * For the local Lax--Friedrichs flux, the definition is @f$\hat{\mathbf{F}}
1447  * =\frac{\mathbf{F}(\mathbf{w}^-)+\mathbf{F}(\mathbf{w}^+)}{2} +
1448  * \frac{\lambda}{2}\left[\mathbf{w}^--\mathbf{w}^+\right]\otimes
1449  * \mathbf{n^-}@f$, where the factor @f$\lambda =
1450  * \max\left(\|\mathbf{u}^-\|+c^-, \|\mathbf{u}^+\|+c^+\right)@f$ gives the
1451  * maximal wave speed and @f$c = \sqrt{\lambda p / \rho}@f$ is the speed of
1452  * sound. Here, we choose two modifications of that expression for reasons
1453  * of computational efficiency, given the small impact of the flux on the
1454  * solution. For the above definition of the factor @f$\lambda@f$, we would need
1455  * to take four square roots, two for the two velocity norms and two for the
1456  * speed of sound on either side. The first modification is hence to rather
1457  * use @f$\sqrt{\|\mathbf{u}\|^2+c^2}@f$ as an estimate of the maximal speed
1458  * (which is at most a factor of 2 away from the actual maximum, as shown in
1459  * the introduction). This allows us to pull the square root out of the
1460  * maximum and get away with a single square root computation. The second
1461  * modification is to further relax on the parameter @f$\lambda@f$---the smaller
1462  * it is, the smaller the dissipation factor (which is multiplied by the
1463  * jump in @f$\mathbf{w}@f$, which might result in a smaller or bigger
1464  * dissipation in the end). This allows us to fit the spectrum into the
1465  * stability region of the explicit Runge--Kutta integrator with bigger time
1466  * steps. However, we cannot make dissipation too small because otherwise
1467  * imaginary eigenvalues grow larger. Finally, the current conservative
1468  * formulation is not energy-stable in the limit of @f$\lambda\to 0@f$ as it is
1469  * not skew-symmetric, and would need additional measures such as split-form
1470  * DG schemes in that case.
1471  *
1472 
1473  *
1474  * For the HLL flux, we follow the formula from literature, introducing an
1475  * additional weighting of the two states from Lax--Friedrichs by a
1476  * parameter @f$s@f$. It is derived from the physical transport directions of
1477  * the Euler equations in terms of the current direction of velocity and
1478  * sound speed. For the velocity, we here choose a simple arithmetic average
1479  * which is sufficient for DG scenarios and moderate jumps in material
1480  * parameters.
1481  *
1482 
1483  *
1484  * Since the numerical flux is multiplied by the normal vector in the weak
1485  * form, we multiply by the result by the normal vector for all terms in the
1486  * equation. In these multiplications, the `operator*` defined above enables
1487  * a compact notation similar to the mathematical definition.
1488  *
1489 
1490  *
1491  * In this and the following functions, we use variable suffixes `_m` and
1492  * `_p` to indicate quantities derived from @f$\mathbf{w}^-@f$ and @f$\mathbf{w}^+@f$,
1493  * i.e., values "here" and "there" relative to the current cell when looking
1494  * at a neighbor cell.
1495  *
1496  * @code
1497  * template <int dim, typename Number>
1498  * inline DEAL_II_ALWAYS_INLINE
1500  * euler_numerical_flux(const Tensor<1, dim + 2, Number> &u_m,
1501  * const Tensor<1, dim + 2, Number> &u_p,
1502  * const Tensor<1, dim, Number> & normal)
1503  * {
1504  * const auto velocity_m = euler_velocity<dim>(u_m);
1505  * const auto velocity_p = euler_velocity<dim>(u_p);
1506  *
1507  * const auto pressure_m = euler_pressure<dim>(u_m);
1508  * const auto pressure_p = euler_pressure<dim>(u_p);
1509  *
1510  * const auto flux_m = euler_flux<dim>(u_m);
1511  * const auto flux_p = euler_flux<dim>(u_p);
1512  *
1513  * switch (numerical_flux_type)
1514  * {
1515  * case lax_friedrichs_modified:
1516  * {
1517  * const auto lambda =
1518  * 0.5 * std::sqrt(std::max(velocity_p.norm_square() +
1519  * gamma * pressure_p * (1. / u_p[0]),
1520  * velocity_m.norm_square() +
1521  * gamma * pressure_m * (1. / u_m[0])));
1522  *
1523  * return 0.5 * (flux_m * normal + flux_p * normal) +
1524  * 0.5 * lambda * (u_m - u_p);
1525  * }
1526  *
1527  * case harten_lax_vanleer:
1528  * {
1529  * const auto avg_velocity_normal =
1530  * 0.5 * ((velocity_m + velocity_p) * normal);
1531  * const auto avg_c = std::sqrt(std::abs(
1532  * 0.5 * gamma *
1533  * (pressure_p * (1. / u_p[0]) + pressure_m * (1. / u_m[0]))));
1534  * const Number s_pos =
1535  * std::max(Number(), avg_velocity_normal + avg_c);
1536  * const Number s_neg =
1537  * std::min(Number(), avg_velocity_normal - avg_c);
1538  * const Number inverse_s = Number(1.) / (s_pos - s_neg);
1539  *
1540  * return inverse_s *
1541  * ((s_pos * (flux_m * normal) - s_neg * (flux_p * normal)) -
1542  * s_pos * s_neg * (u_m - u_p));
1543  * }
1544  *
1545  * default:
1546  * {
1547  * Assert(false, ExcNotImplemented());
1548  * return {};
1549  * }
1550  * }
1551  * }
1552  *
1553  *
1554  *
1555  * @endcode
1556  *
1557  * This and the next function are helper functions to provide compact
1558  * evaluation calls as multiple points get batched together via a
1559  * VectorizedArray argument (see the @ref step_37 "step-37" tutorial for details). This
1560  * function is used for the subsonic outflow boundary conditions where we
1561  * need to set the energy component to a prescribed value. The next one
1562  * requests the solution on all components and is used for inflow boundaries
1563  * where all components of the solution are set.
1564  *
1565  * @code
1566  * template <int dim, typename Number>
1568  * evaluate_function(const Function<dim> & function,
1569  * const Point<dim, VectorizedArray<Number>> &p_vectorized,
1570  * const unsigned int component)
1571  * {
1572  * VectorizedArray<Number> result;
1573  * for (unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v)
1574  * {
1575  * Point<dim> p;
1576  * for (unsigned int d = 0; d < dim; ++d)
1577  * p[d] = p_vectorized[d][v];
1578  * result[v] = function.value(p, component);
1579  * }
1580  * return result;
1581  * }
1582  *
1583  *
1584  * template <int dim, typename Number, int n_components = dim + 2>
1586  * evaluate_function(const Function<dim> & function,
1587  * const Point<dim, VectorizedArray<Number>> &p_vectorized)
1588  * {
1591  * for (unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v)
1592  * {
1593  * Point<dim> p;
1594  * for (unsigned int d = 0; d < dim; ++d)
1595  * p[d] = p_vectorized[d][v];
1596  * for (unsigned int d = 0; d < n_components; ++d)
1597  * result[d][v] = function.value(p, d);
1598  * }
1599  * return result;
1600  * }
1601  *
1602  *
1603  *
1604  * @endcode
1605  *
1606  *
1607  * <a name="TheEulerOperationclass"></a>
1608  * <h3>The EulerOperation class</h3>
1609  *
1610 
1611  *
1612  * This class implements the evaluators for the Euler problem, in analogy to
1613  * the `LaplaceOperator` class of @ref step_37 "step-37" or @ref step_59 "step-59". Since the present
1614  * operator is non-linear and does not require a matrix interface (to be
1615  * handed over to preconditioners), we skip the various `vmult` functions
1616  * otherwise present in matrix-free operators and only implement an `apply`
1617  * function as well as the combination of `apply` with the required vector
1618  * updates for the low-storage Runge--Kutta time integrator mentioned above
1619  * (called `perform_stage`). Furthermore, we have added three additional
1620  * functions involving matrix-free routines, namely one to compute an
1621  * estimate of the time step scaling (that is combined with the Courant
1622  * number for the actual time step size) based on the velocity and speed of
1623  * sound in the elements, one for the projection of solutions (specializing
1624  * VectorTools::project() for the DG case), and one to compute the errors
1625  * against a possible analytical solution or norms against some background
1626  * state.
1627  *
1628 
1629  *
1630  * The rest of the class is similar to other matrix-free tutorials. As
1631  * discussed in the introduction, we provide a few functions to allow a user
1632  * to pass in various forms of boundary conditions on different parts of the
1633  * domain boundary marked by types::boundary_id variables, as well as
1634  * possible body forces.
1635  *
1636  * @code
1637  * template <int dim, int degree, int n_points_1d>
1638  * class EulerOperator
1639  * {
1640  * public:
1641  * static constexpr unsigned int n_quadrature_points_1d = n_points_1d;
1642  *
1643  * EulerOperator(TimerOutput &timer_output);
1644  *
1645  * void reinit(const Mapping<dim> & mapping,
1646  * const DoFHandler<dim> &dof_handler);
1647  *
1648  * void set_inflow_boundary(const types::boundary_id boundary_id,
1649  * std::unique_ptr<Function<dim>> inflow_function);
1650  *
1651  * void set_subsonic_outflow_boundary(
1653  * std::unique_ptr<Function<dim>> outflow_energy);
1654  *
1655  * void set_wall_boundary(const types::boundary_id boundary_id);
1656  *
1657  * void set_body_force(std::unique_ptr<Function<dim>> body_force);
1658  *
1659  * void apply(const double current_time,
1662  *
1663  * void
1664  * perform_stage(const Number cur_time,
1665  * const Number factor_solution,
1666  * const Number factor_ai,
1667  * const LinearAlgebra::distributed::Vector<Number> &current_ri,
1671  *
1672  * void project(const Function<dim> & function,
1673  * LinearAlgebra::distributed::Vector<Number> &solution) const;
1674  *
1675  * std::array<double, 3> compute_errors(
1676  * const Function<dim> & function,
1677  * const LinearAlgebra::distributed::Vector<Number> &solution) const;
1678  *
1679  * double compute_cell_transport_speed(
1680  * const LinearAlgebra::distributed::Vector<Number> &solution) const;
1681  *
1682  * void
1683  * initialize_vector(LinearAlgebra::distributed::Vector<Number> &vector) const;
1684  *
1685  * private:
1686  * MatrixFree<dim, Number> data;
1687  *
1688  * TimerOutput &timer;
1689  *
1690  * std::map<types::boundary_id, std::unique_ptr<Function<dim>>>
1691  * inflow_boundaries;
1692  * std::map<types::boundary_id, std::unique_ptr<Function<dim>>>
1693  * subsonic_outflow_boundaries;
1694  * std::set<types::boundary_id> wall_boundaries;
1695  * std::unique_ptr<Function<dim>> body_force;
1696  *
1697  * void local_apply_inverse_mass_matrix(
1698  * const MatrixFree<dim, Number> & data,
1701  * const std::pair<unsigned int, unsigned int> & cell_range) const;
1702  *
1703  * void local_apply_cell(
1704  * const MatrixFree<dim, Number> & data,
1707  * const std::pair<unsigned int, unsigned int> & cell_range) const;
1708  *
1709  * void local_apply_face(
1710  * const MatrixFree<dim, Number> & data,
1713  * const std::pair<unsigned int, unsigned int> & cell_range) const;
1714  *
1715  * void local_apply_boundary_face(
1716  * const MatrixFree<dim, Number> & data,
1719  * const std::pair<unsigned int, unsigned int> & cell_range) const;
1720  * };
1721  *
1722  *
1723  *
1724  * template <int dim, int degree, int n_points_1d>
1725  * EulerOperator<dim, degree, n_points_1d>::EulerOperator(TimerOutput &timer)
1726  * : timer(timer)
1727  * {}
1728  *
1729  *
1730  *
1731  * @endcode
1732  *
1733  * For the initialization of the Euler operator, we set up the MatrixFree
1734  * variable contained in the class. This can be done given a mapping to
1735  * describe possible curved boundaries as well as a DoFHandler object
1736  * describing the degrees of freedom. Since we use a discontinuous Galerkin
1737  * discretization in this tutorial program where no constraints are imposed
1738  * strongly on the solution field, we do not need to pass in an
1739  * AffineConstraints object and rather use a dummy for the
1740  * construction. With respect to quadrature, we want to select two different
1741  * ways of computing the underlying integrals: The first is a flexible one,
1742  * based on a template parameter `n_points_1d` (that will be assigned the
1743  * `n_q_points_1d` value specified at the top of this file). More accurate
1744  * integration is necessary to avoid the aliasing problem due to the
1745  * variable coefficients in the Euler operator. The second less accurate
1746  * quadrature formula is a tight one based on `fe_degree+1` and needed for
1747  * the inverse mass matrix. While that formula provides an exact inverse
1748  * only on affine element shapes and not on deformed elements, it enables
1749  * the fast inversion of the mass matrix by tensor product techniques,
1750  * necessary to ensure optimal computational efficiency overall.
1751  *
1752  * @code
1753  * template <int dim, int degree, int n_points_1d>
1755  * const Mapping<dim> & mapping,
1756  * const DoFHandler<dim> &dof_handler)
1757  * {
1758  * const std::vector<const DoFHandler<dim> *> dof_handlers = {&dof_handler};
1760  * const std::vector<const AffineConstraints<double> *> constraints = {&dummy};
1761  * const std::vector<Quadrature<1>> quadratures = {QGauss<1>(n_q_points_1d),
1762  * QGauss<1>(fe_degree + 1)};
1763  *
1764  * typename MatrixFree<dim, Number>::AdditionalData additional_data;
1765  * additional_data.mapping_update_flags =
1767  * update_values);
1768  * additional_data.mapping_update_flags_inner_faces =
1770  * update_values);
1771  * additional_data.mapping_update_flags_boundary_faces =
1773  * update_values);
1774  * additional_data.tasks_parallel_scheme =
1776  *
1777  * data.reinit(
1778  * mapping, dof_handlers, constraints, quadratures, additional_data);
1779  * }
1780  *
1781  *
1782  *
1783  * template <int dim, int degree, int n_points_1d>
1784  * void EulerOperator<dim, degree, n_points_1d>::initialize_vector(
1786  * {
1787  * data.initialize_dof_vector(vector);
1788  * }
1789  *
1790  *
1791  *
1792  * @endcode
1793  *
1794  * The subsequent four member functions are the ones that must be called from
1795  * outside to specify the various types of boundaries. For an inflow boundary,
1796  * we must specify all components in terms of density @f$\rho@f$, momentum @f$\rho
1797  * \mathbf{u}@f$ and energy @f$E@f$. Given this information, we then store the
1798  * function alongside the respective boundary id in a map member variable of
1799  * this class. Likewise, we proceed for the subsonic outflow boundaries (where
1800  * we request a function as well, which we use to retrieve the energy) and for
1801  * wall (no-penetration) boundaries where we impose zero normal velocity (no
1802  * function necessary, so we only request the boundary id). For the present
1803  * DG code where boundary conditions are solely applied as part of the weak
1804  * form (during time integration), the call to set the boundary conditions
1805  * can appear both before or after the `reinit()` call to this class. This
1806  * is different from continuous finite element codes where the boundary
1807  * conditions determine the content of the AffineConstraints object that is
1808  * sent into MatrixFree for initialization, thus requiring to be set before
1809  * the initialization of the matrix-free data structures.
1810  *
1811 
1812  *
1813  * The checks added in each of the four function are used to
1814  * ensure that boundary conditions are mutually exclusive on the various
1815  * parts of the boundary, i.e., that a user does not accidentally designate a
1816  * boundary as both an inflow and say a subsonic outflow boundary.
1817  *
1818  * @code
1819  * template <int dim, int degree, int n_points_1d>
1820  * void EulerOperator<dim, degree, n_points_1d>::set_inflow_boundary(
1822  * std::unique_ptr<Function<dim>> inflow_function)
1823  * {
1824  * AssertThrow(subsonic_outflow_boundaries.find(boundary_id) ==
1825  * subsonic_outflow_boundaries.end() &&
1826  * wall_boundaries.find(boundary_id) == wall_boundaries.end(),
1827  * ExcMessage("You already set the boundary with id " +
1828  * std::to_string(static_cast<int>(boundary_id)) +
1829  * " to another type of boundary before now setting " +
1830  * "it as inflow"));
1831  * AssertThrow(inflow_function->n_components == dim + 2,
1832  * ExcMessage("Expected function with dim+2 components"));
1833  *
1834  * inflow_boundaries[boundary_id] = std::move(inflow_function);
1835  * }
1836  *
1837  *
1838  * template <int dim, int degree, int n_points_1d>
1839  * void EulerOperator<dim, degree, n_points_1d>::set_subsonic_outflow_boundary(
1841  * std::unique_ptr<Function<dim>> outflow_function)
1842  * {
1843  * AssertThrow(inflow_boundaries.find(boundary_id) ==
1844  * inflow_boundaries.end() &&
1845  * wall_boundaries.find(boundary_id) == wall_boundaries.end(),
1846  * ExcMessage("You already set the boundary with id " +
1847  * std::to_string(static_cast<int>(boundary_id)) +
1848  * " to another type of boundary before now setting " +
1849  * "it as subsonic outflow"));
1850  * AssertThrow(outflow_function->n_components == dim + 2,
1851  * ExcMessage("Expected function with dim+2 components"));
1852  *
1853  * subsonic_outflow_boundaries[boundary_id] = std::move(outflow_function);
1854  * }
1855  *
1856  *
1857  * template <int dim, int degree, int n_points_1d>
1858  * void EulerOperator<dim, degree, n_points_1d>::set_wall_boundary(
1860  * {
1861  * AssertThrow(inflow_boundaries.find(boundary_id) ==
1862  * inflow_boundaries.end() &&
1863  * subsonic_outflow_boundaries.find(boundary_id) ==
1864  * subsonic_outflow_boundaries.end(),
1865  * ExcMessage("You already set the boundary with id " +
1866  * std::to_string(static_cast<int>(boundary_id)) +
1867  * " to another type of boundary before now setting " +
1868  * "it as wall boundary"));
1869  *
1870  * wall_boundaries.insert(boundary_id);
1871  * }
1872  *
1873  *
1874  * template <int dim, int degree, int n_points_1d>
1875  * void EulerOperator<dim, degree, n_points_1d>::set_body_force(
1876  * std::unique_ptr<Function<dim>> body_force)
1877  * {
1878  * AssertDimension(body_force->n_components, dim);
1879  *
1880  * this->body_force = std::move(body_force);
1881  * }
1882  *
1883  *
1884  *
1885  * @endcode
1886  *
1887  *
1888  * <a name="Localevaluators"></a>
1889  * <h4>Local evaluators</h4>
1890  *
1891 
1892  *
1893  * Now we proceed to the local evaluators for the Euler problem. The
1894  * evaluators are relatively simple and follow what has been presented in
1895  * @ref step_37 "step-37", @ref step_48 "step-48", or @ref step_59 "step-59". The first notable difference is the fact
1896  * that we use an FEEvaluation with a non-standard number of quadrature
1897  * points. Whereas we previously always set the number of quadrature points
1898  * to equal the polynomial degree plus one (ensuring exact integration on
1899  * affine element shapes), we now set the number quadrature points as a
1900  * separate variable (e.g. the polynomial degree plus two or three halves of
1901  * the polynomial degree) to more accurately handle nonlinear terms. Since
1902  * the evaluator is fed with the appropriate loop lengths via the template
1903  * argument and keeps the number of quadrature points in the whole cell in
1904  * the variable FEEvaluation::n_q_points, we now automatically operate on
1905  * the more accurate formula without further changes.
1906  *
1907 
1908  *
1909  * The second difference is due to the fact that we are now evaluating a
1910  * multi-component system, as opposed to the scalar systems considered
1911  * previously. The matrix-free framework provides several ways to handle the
1912  * multi-component case. The variant shown here utilizes an FEEvaluation
1913  * object with multiple components embedded into it, specified by the fourth
1914  * template argument `dim + 2` for the components in the Euler system. As a
1915  * consequence, the return type of FEEvaluation::get_value() is not a scalar
1916  * any more (that would return a VectorizedArray type, collecting data from
1917  * several elements), but a Tensor of `dim+2` components. The functionality
1918  * is otherwise similar to the scalar case; it is handled by a template
1919  * specialization of a base class, called FEEvaluationAccess. An alternative
1920  * variant would have been to use several FEEvaluation objects, a scalar one
1921  * for the density, a vector-valued one with `dim` components for the
1922  * momentum, and another scalar evaluator for the energy. To ensure that
1923  * those components point to the correct part of the solution, the
1924  * constructor of FEEvaluation takes three optional integer arguments after
1925  * the required MatrixFree field, namely the number of the DoFHandler for
1926  * multi-DoFHandler systems (taking the first by default), the number of the
1927  * quadrature point in case there are multiple Quadrature objects (see more
1928  * below), and as a third argument the component within a vector system. As
1929  * we have a single vector for all components, we would go with the third
1930  * argument, and set it to `0` for the density, `1` for the vector-valued
1931  * momentum, and `dim+1` for the energy slot. FEEvaluation then picks the
1932  * appropriate subrange of the solution vector during
1933  * FEEvaluationBase::read_dof_values() and
1934  * FEEvaluation::distributed_local_to_global() or the more compact
1935  * FEEvaluation::gather_evaluate() and FEEvaluation::integrate_scatter()
1936  * calls.
1937  *
1938 
1939  *
1940  * When it comes to the evaluation of the body force vector, we distinguish
1941  * between two cases for efficiency reasons: In case we have a constant
1942  * function (derived from Functions::ConstantFunction), we can precompute
1943  * the value outside the loop over quadrature points and simply use the
1944  * value everywhere. For a more general function, we instead need to call
1945  * the `evaluate_function()` method we provided above; this path is more
1946  * expensive because we need to access the memory associated with the
1947  * quadrature point data.
1948  *
1949 
1950  *
1951  * The rest follows the other tutorial programs. Since we have implemented
1952  * all physics for the Euler equations in the separate `euler_flux()`
1953  * function, all we have to do here is to call this function
1954  * given the current solution evaluated at quadrature points, returned by
1955  * `phi.get_value(q)`, and tell the FEEvaluation object to queue the flux
1956  * for testing it by the gradients of the shape functions (which is a Tensor
1957  * of outer `dim+2` components, each holding a tensor of `dim` components
1958  * for the @f$x,y,z@f$ component of the Euler flux). One final thing worth
1959  * mentioning is the order in which we queue the data for testing by the
1960  * value of the test function, `phi.submit_value()`, in case we are given an
1961  * external function: We must do this after calling `phi.get_value(q)`,
1962  * because `get_value()` (reading the solution) and `submit_value()`
1963  * (queuing the value for multiplication by the test function and summation
1964  * over quadrature points) access the same underlying data field. Here it
1965  * would be easy to achieve also without temporary variable `w_q` since
1966  * there is no mixing between values and gradients. For more complicated
1967  * setups, one has to first copy out e.g. both the value and gradient at a
1968  * quadrature point and then queue results again by
1969  * FEEvaluationBase::submit_value() and FEEvaluationBase::submit_gradient().
1970  *
1971 
1972  *
1973  * As a final note, we mention that we do not use the first MatrixFree
1974  * argument of this function, which is a call-back from MatrixFree::loop().
1975  * The interfaces imposes the present list of arguments, but since we are in
1976  * a member function where the MatrixFree object is already available as the
1977  * `data` variable, we stick with that to avoid confusion.
1978  *
1979  * @code
1980  * template <int dim, int degree, int n_points_1d>
1981  * void EulerOperator<dim, degree, n_points_1d>::local_apply_cell(
1982  * const MatrixFree<dim, Number> &,
1983  * LinearAlgebra::distributed::Vector<Number> & dst,
1984  * const LinearAlgebra::distributed::Vector<Number> &src,
1985  * const std::pair<unsigned int, unsigned int> & cell_range) const
1986  * {
1988  *
1989  * Tensor<1, dim, VectorizedArray<Number>> constant_body_force;
1990  * const Functions::ConstantFunction<dim> *constant_function =
1991  * dynamic_cast<Functions::ConstantFunction<dim> *>(body_force.get());
1992  *
1993  * if (constant_function)
1994  * constant_body_force = evaluate_function<dim, Number, dim>(
1995  * *constant_function, Point<dim, VectorizedArray<Number>>());
1996  *
1997  * for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
1998  * {
1999  * phi.reinit(cell);
2000  * phi.gather_evaluate(src, true, false);
2001  *
2002  * for (unsigned int q = 0; q < phi.n_q_points; ++q)
2003  * {
2004  * const auto w_q = phi.get_value(q);
2005  * phi.submit_gradient(euler_flux<dim>(w_q), q);
2006  * if (body_force.get() != nullptr)
2007  * {
2009  * constant_function ? constant_body_force :
2010  * evaluate_function<dim, Number, dim>(
2011  * *body_force, phi.quadrature_point(q));
2012  *
2014  * for (unsigned int d = 0; d < dim; ++d)
2015  * forcing[d + 1] = w_q[0] * force[d];
2016  * for (unsigned int d = 0; d < dim; ++d)
2017  * forcing[dim + 1] += force[d] * w_q[d + 1];
2018  *
2019  * phi.submit_value(forcing, q);
2020  * }
2021  * }
2022  *
2023  * phi.integrate_scatter(body_force.get() != nullptr, true, dst);
2024  * }
2025  * }
2026  *
2027  *
2028  *
2029  * @endcode
2030  *
2031  * The next function concerns the computation of integrals on interior
2032  * faces, where we need evaluators from both cells adjacent to the face. We
2033  * associate the variable `phi_m` with the solution component @f$\mathbf{w}^-@f$
2034  * and the variable `phi_p` with the solution component @f$\mathbf{w}^+@f$. We
2035  * distinguish the two sides in the constructor of FEFaceEvaluation by the
2036  * second argument, with `true` for the interior side and `false` for the
2037  * exterior side, with interior and exterior denoting the orientation with
2038  * respect to the normal vector.
2039  *
2040 
2041  *
2042  * Note that the calls FEFaceEvaluation::gather_evaluate() and
2043  * FEFaceEvaluation::integrate_scatter() combine the access to the vectors
2044  * and the sum factorization parts. This combined operation not only saves a
2045  * line of code, but also contains an important optimization: Given that we
2046  * use a nodal basis in terms of the Lagrange polynomials in the points of
2047  * the Gauss-Lobatto quadrature formula, only @f$(p+1)^{d-1}@f$ out of the
2048  * @f$(p+1)^d@f$ basis functions evaluate to non-zero on each face. Thus, the
2049  * evaluator only accesses the necessary data in the vector and skips the
2050  * parts which are multiplied by zero. If we had first read the vector, we
2051  * would have needed to load all data from the vector, as the call in
2052  * isolation would not know what data is required in subsequent
2053  * operations. If the subsequent FEFaceEvaluation::evaluate() call requests
2054  * values and derivatives, indeed all @f$(p+1)^d@f$ vector entries for each
2055  * component are needed, as the normal derivative is nonzero for all basis
2056  * functions.
2057  *
2058 
2059  *
2060  * The arguments to the evaluators as well as the procedure is similar to
2061  * the cell evaluation. We again use the more accurate (over-)integration
2062  * scheme due to the nonlinear terms, specified as the third template
2063  * argument in the list. At the quadrature points, we then go to our
2064  * free-standing function for the numerical flux. It receives the solution
2065  * evaluated at quadrature points from both sides (i.e., @f$\mathbf{w}^-@f$ and
2066  * @f$\mathbf{w}^+@f$), as well as the normal vector onto the minus side. As
2067  * explained above, the numerical flux is already multiplied by the normal
2068  * vector from the minus side. We need to switch the sign because the
2069  * boundary term comes with a minus sign in the weak form derived in the
2070  * introduction. The flux is then queued for testing both on the minus sign
2071  * and on the plus sign, with switched sign as the normal vector from the
2072  * plus side is exactly opposed to the one from the minus side.
2073  *
2074  * @code
2075  * template <int dim, int degree, int n_points_1d>
2076  * void EulerOperator<dim, degree, n_points_1d>::local_apply_face(
2077  * const MatrixFree<dim, Number> &,
2080  * const std::pair<unsigned int, unsigned int> & face_range) const
2081  * {
2083  * true);
2085  * false);
2086  *
2087  * for (unsigned int face = face_range.first; face < face_range.second; ++face)
2088  * {
2089  * phi_p.reinit(face);
2090  * phi_p.gather_evaluate(src, true, false);
2091  *
2092  * phi_m.reinit(face);
2093  * phi_m.gather_evaluate(src, true, false);
2094  *
2095  * for (unsigned int q = 0; q < phi_m.n_q_points; ++q)
2096  * {
2097  * const auto numerical_flux =
2098  * euler_numerical_flux<dim>(phi_m.get_value(q),
2099  * phi_p.get_value(q),
2100  * phi_m.get_normal_vector(q));
2101  * phi_m.submit_value(-numerical_flux, q);
2102  * phi_p.submit_value(numerical_flux, q);
2103  * }
2104  *
2105  * phi_p.integrate_scatter(true, false, dst);
2106  * phi_m.integrate_scatter(true, false, dst);
2107  * }
2108  * }
2109  *
2110  *
2111  *
2112  * @endcode
2113  *
2114  * For faces located at the boundary, we need to impose the appropriate
2115  * boundary conditions. In this tutorial program, we implement four cases as
2116  * mentioned above. (A fifth case, for supersonic outflow conditions is
2117  * discussed in the "Results" section below.) The discontinuous Galerkin
2118  * method imposes boundary conditions not as constraints, but only
2119  * weakly. Thus, the various conditions are imposed by finding an appropriate
2120  * <i>exterior</i> quantity @f$\mathbf{w}^+@f$ that is then handed to the
2121  * numerical flux function also used for the interior faces. In essence,
2122  * we "pretend" a state on the outside of the domain in such a way that
2123  * if that were reality, the solution of the PDE would satisfy the boundary
2124  * conditions we want.
2125  *
2126 
2127  *
2128  * For wall boundaries, we need to impose a no-normal-flux condition on the
2129  * momentum variable, whereas we use a Neumann condition for the density and
2130  * energy with @f$\rho^+ = \rho^-@f$ and @f$E^+ = E^-@f$. To achieve the no-normal
2131  * flux condition, we set the exterior values to the interior values and
2132  * subtract two times the velocity in wall-normal direction, i.e., in the
2133  * direction of the normal vector.
2134  *
2135 
2136  *
2137  * For inflow boundaries, we simply set the given Dirichlet data
2138  * @f$\mathbf{w}_\mathrm{D}@f$ as a boundary value. An alternative would have been
2139  * to use @f$\mathbf{w}^+ = -\mathbf{w}^- + 2 \mathbf{w}_\mathrm{D}@f$, the
2140  * so-called mirror principle.
2141  *
2142 
2143  *
2144  * The imposition of outflow is essentially a Neumann condition, i.e.,
2145  * setting @f$\mathbf{w}^+ = \mathbf{w}^-@f$. For the case of subsonic outflow,
2146  * we still need to impose a value for the energy, which we derive from the
2147  * respective function. A special step is needed for the case of
2148  * <i>backflow</i>, i.e., the case where there is a momentum flux into the
2149  * domain on the Neumann portion. According to the literature (a fact that can
2150  * be derived by appropriate energy arguments), we must switch to another
2151  * variant of the flux on inflow parts, see Gravemeier, Comerford,
2152  * Yoshihara, Ismail, Wall, "A novel formulation for Neumann inflow
2153  * conditions in biomechanics", Int. J. Numer. Meth. Biomed. Eng., vol. 28
2154  * (2012). Here, the momentum term needs to be added once again, which
2155  * corresponds to removing the flux contribution on the momentum
2156  * variables. We do this in a post-processing step, and only for the case
2157  * when we both are at an outflow boundary and the dot product between the
2158  * normal vector and the momentum (or, equivalently, velocity) is
2159  * negative. As we work on data of several quadrature points at once for
2160  * SIMD vectorizations, we here need to explicitly loop over the array
2161  * entries of the SIMD array.
2162  *
2163 
2164  *
2165  * In the implementation below, we check for the various types
2166  * of boundaries at the level of quadrature points. Of course, we could also
2167  * have moved the decision out of the quadrature point loop and treat entire
2168  * faces as of the same kind, which avoids some map/set lookups in the inner
2169  * loop over quadrature points. However, the loss of efficiency is hardly
2170  * noticeable, so we opt for the simpler code here. Also note that the final
2171  * `else` clause will catch the case when some part of the boundary was not
2172  * assigned any boundary condition via `EulerOperator::set_..._boundary(...)`.
2173  *
2174  * @code
2175  * template <int dim, int degree, int n_points_1d>
2176  * void EulerOperator<dim, degree, n_points_1d>::local_apply_boundary_face(
2177  * const MatrixFree<dim, Number> &,
2180  * const std::pair<unsigned int, unsigned int> & face_range) const
2181  * {
2183  *
2184  * for (unsigned int face = face_range.first; face < face_range.second; ++face)
2185  * {
2186  * phi.reinit(face);
2187  * phi.gather_evaluate(src, true, false);
2188  *
2189  * for (unsigned int q = 0; q < phi.n_q_points; ++q)
2190  * {
2191  * const auto w_m = phi.get_value(q);
2192  * const auto normal = phi.get_normal_vector(q);
2193  *
2194  * auto rho_u_dot_n = w_m[1] * normal[0];
2195  * for (unsigned int d = 1; d < dim; ++d)
2196  * rho_u_dot_n += w_m[1 + d] * normal[d];
2197  *
2198  * bool at_outflow = false;
2199  *
2201  * const auto boundary_id = data.get_boundary_id(face);
2202  * if (wall_boundaries.find(boundary_id) != wall_boundaries.end())
2203  * {
2204  * w_p[0] = w_m[0];
2205  * for (unsigned int d = 0; d < dim; ++d)
2206  * w_p[d + 1] = w_m[d + 1] - 2. * rho_u_dot_n * normal[d];
2207  * w_p[dim + 1] = w_m[dim + 1];
2208  * }
2209  * else if (inflow_boundaries.find(boundary_id) !=
2210  * inflow_boundaries.end())
2211  * w_p =
2212  * evaluate_function(*inflow_boundaries.find(boundary_id)->second,
2213  * phi.quadrature_point(q));
2214  * else if (subsonic_outflow_boundaries.find(boundary_id) !=
2215  * subsonic_outflow_boundaries.end())
2216  * {
2217  * w_p = w_m;
2218  * w_p[dim + 1] = evaluate_function(
2219  * *subsonic_outflow_boundaries.find(boundary_id)->second,
2220  * phi.quadrature_point(q),
2221  * dim + 1);
2222  * at_outflow = true;
2223  * }
2224  * else
2225  * AssertThrow(false,
2226  * ExcMessage("Unknown boundary id, did "
2227  * "you set a boundary condition for "
2228  * "this part of the domain boundary?"));
2229  *
2230  * auto flux = euler_numerical_flux<dim>(w_m, w_p, normal);
2231  *
2232  * if (at_outflow)
2233  * for (unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v)
2234  * {
2235  * if (rho_u_dot_n[v] < -1e-12)
2236  * for (unsigned int d = 0; d < dim; ++d)
2237  * flux[d + 1][v] = 0.;
2238  * }
2239  *
2240  * phi.submit_value(-flux, q);
2241  * }
2242  *
2243  * phi.integrate_scatter(true, false, dst);
2244  * }
2245  * }
2246  *
2247  *
2248  *
2249  * @endcode
2250  *
2251  * The next function implements the inverse mass matrix operation. The
2252  * algorithms and rationale have been discussed extensively in the
2253  * introduction, so we here limit ourselves to the technicalities of the
2254  * MatrixFreeOperators::CellwiseInverseMassMatrix class. It does similar
2255  * operations as the forward evaluation of the mass matrix, except with a
2256  * different interpolation matrix, representing the inverse @f$S^{-1}@f$
2257  * factors. These represent a change of basis from the specified basis (in
2258  * this case, the Lagrange basis in the points of the Gauss--Lobatto
2259  * quadrature formula) to the Lagrange basis in the points of the Gauss
2260  * quadrature formula. In the latter basis, we can apply the inverse of the
2261  * point-wise `JxW` factor, i.e., the quadrature weight times the
2262  * determinant of the Jacobian of the mapping from reference to real
2263  * coordinates. Once this is done, the basis is changed back to the nodal
2264  * Gauss-Lobatto basis again. All of these operations are done by the
2265  * `apply()` function below. What we need to provide is the local fields to
2266  * operate on (which we extract from the global vector by an FEEvaluation
2267  * object) and write the results back to the destination vector of the mass
2268  * matrix operation.
2269  *
2270 
2271  *
2272  * One thing to note is that we added two integer arguments (that are
2273  * optional) to the constructor of FEEvaluation, the first being 0
2274  * (selecting among the DoFHandler in multi-DoFHandler systems; here, we
2275  * only have one) and the second being 1 to make the quadrature formula
2276  * selection. As we use the quadrature formula 0 for the over-integration of
2277  * nonlinear terms, we use the formula 1 with the default @f$p+1@f$ (or
2278  * `fe_degree+1` in terms of the variable name) points for the mass
2279  * matrix. This leads to square contributions to the mass matrix and ensures
2280  * exact integration, as explained in the introduction.
2281  *
2282  * @code
2283  * template <int dim, int degree, int n_points_1d>
2284  * void EulerOperator<dim, degree, n_points_1d>::local_apply_inverse_mass_matrix(
2285  * const MatrixFree<dim, Number> &,
2288  * const std::pair<unsigned int, unsigned int> & cell_range) const
2289  * {
2292  * inverse(phi);
2293  *
2294  * for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2295  * {
2296  * phi.reinit(cell);
2297  * phi.read_dof_values(src);
2298  *
2299  * inverse.apply(phi.begin_dof_values(), phi.begin_dof_values());
2300  *
2301  * phi.set_dof_values(dst);
2302  * }
2303  * }
2304  *
2305  *
2306  *
2307  * @endcode
2308  *
2309  *
2310  * <a name="Theapplyandrelatedfunctions"></a>
2311  * <h4>The apply() and related functions</h4>
2312  *
2313 
2314  *
2315  * We now come to the function which implements the evaluation of the Euler
2316  * operator as a whole, i.e., @f$\mathcal M^{-1} \mathcal L(t, \mathbf{w})@f$,
2317  * calling into the local evaluators presented above. The steps should be
2318  * clear from the previous code. One thing to note is that we need to adjust
2319  * the time in the functions we have associated with the various parts of
2320  * the boundary, in order to be consistent with the equation in case the
2321  * boundary data is time-dependent. Then, we call MatrixFree::loop() to
2322  * perform the cell and face integrals, including the necessary ghost data
2323  * exchange in the `src` vector. The seventh argument to the function,
2324  * `true`, specifies that we want to zero the `dst` vector as part of the
2325  * loop, before we start accumulating integrals into it. This variant is
2326  * preferred over explicitly calling `dst = 0.;` before the loop as the
2327  * zeroing operation is done on a subrange of the vector in parts that are
2328  * written by the integrals nearby. This enhances data locality and allows
2329  * for caching, saving one roundtrip of vector data to main memory and
2330  * enhancing performance. The last two arguments to the loop determine which
2331  * data is exchanged: Since we only access the values of the shape functions
2332  * one faces, typical of first-order hyperbolic problems, and since we have
2333  * a nodal basis with nodes at the reference element surface, we only need
2334  * to exchange those parts. This again saves precious memory bandwidth.
2335  *
2336 
2337  *
2338  * Once the spatial operator @f$\mathcal L@f$ is applied, we need to make a
2339  * second round and apply the inverse mass matrix. Here, we call
2340  * MatrixFree::cell_loop() since only cell integrals appear. The cell loop
2341  * is cheaper than the full loop as access only goes to the degrees of
2342  * freedom associated with the locally owned cells, which is simply the
2343  * locally owned degrees of freedom for DG discretizations. Thus, no ghost
2344  * exchange is needed here.
2345  *
2346 
2347  *
2348  * Around all these functions, we put timer scopes to record the
2349  * computational time for statistics about the contributions of the various
2350  * parts.
2351  *
2352  * @code
2353  * template <int dim, int degree, int n_points_1d>
2354  * void EulerOperator<dim, degree, n_points_1d>::apply(
2355  * const double current_time,
2356  * const LinearAlgebra::distributed::Vector<Number> &src,
2357  * LinearAlgebra::distributed::Vector<Number> & dst) const
2358  * {
2359  * {
2360  * TimerOutput::Scope t(timer, "apply - integrals");
2361  *
2362  * for (auto &i : inflow_boundaries)
2363  * i.second->set_time(current_time);
2364  * for (auto &i : subsonic_outflow_boundaries)
2365  * i.second->set_time(current_time);
2366  *
2367  * data.loop(&EulerOperator::local_apply_cell,
2368  * &EulerOperator::local_apply_face,
2369  * &EulerOperator::local_apply_boundary_face,
2370  * this,
2371  * dst,
2372  * src,
2373  * true,
2376  * }
2377  *
2378  * {
2379  * TimerOutput::Scope t(timer, "apply - inverse mass");
2380  *
2381  * data.cell_loop(&EulerOperator::local_apply_inverse_mass_matrix,
2382  * this,
2383  * dst,
2384  * dst);
2385  * }
2386  * }
2387  *
2388  *
2389  *
2390  * @endcode
2391  *
2392  * Let us move to the function that does an entire stage of a Runge--Kutta
2393  * update. It calls EulerOperator::apply() followed by some updates
2394  * to the vectors, namely `next_ri = solution + factor_ai * k_i` and
2395  * `solution += factor_solution * k_i`. Rather than performing these
2396  * steps through the vector interfaces, we here present an alternative
2397  * strategy that is faster on cache-based architectures. As the memory
2398  * consumed by the vectors is often much larger than what fits into caches,
2399  * the data has to effectively come from the slow RAM memory. The situation
2400  * can be improved by loop fusion, i.e., performing both the updates to
2401  * `next_ki` and `solution` within a single sweep. In that case, we would
2402  * read the two vectors `rhs` and `solution` and write into `next_ki` and
2403  * `solution`, compared to at least 4 reads and two writes in the baseline
2404  * case. Here, we go one step further and perform the loop immediately when
2405  * the mass matrix inversion has finished on a part of the
2406  * vector. MatrixFree::cell_loop() provides a mechanism to attach an
2407  * `std::function` both before the loop over cells first touches a vector
2408  * entry (which we do not use here, but is e.g. used for zeroing the vector)
2409  * and a second `std::function` to be called after the loop last touches
2410  * an entry. The callback is in form of a range over the given vector (in
2411  * terms of the local index numbering in the MPI universe) that can be
2412  * addressed by `local_element()` functions. For this second callback, we
2413  * create a lambda that works on a range and write the respective update on
2414  * this range. We add the `DEAL_II_OPENMP_SIMD_PRAGMA` before the local loop
2415  * to suggest to the compiler to SIMD parallelize this loop (which means in
2416  * practice that we ensure that there is no overlap, also called
2417  * aliasing, between the index ranges of the pointers we use inside the
2418  * loops). Note that we select a different code path for the last
2419  * Runge--Kutta stage when we do not need to update the `next_ri`
2420  * vector. This strategy gives a considerable speedup. Whereas the inverse
2421  * mass matrix and vector updates take more than 60% of the computational
2422  * time with default vector updates on a 40-core machine, the percentage is
2423  * around 35% with the more optimized variant. In other words, this is a
2424  * speedup of around a third.
2425  *
2426  * @code
2427  * template <int dim, int degree, int n_points_1d>
2428  * void EulerOperator<dim, degree, n_points_1d>::perform_stage(
2429  * const Number current_time,
2430  * const Number factor_solution,
2431  * const Number factor_ai,
2432  * const LinearAlgebra::distributed::Vector<Number> &current_ri,
2433  * LinearAlgebra::distributed::Vector<Number> & vec_ki,
2434  * LinearAlgebra::distributed::Vector<Number> & solution,
2435  * LinearAlgebra::distributed::Vector<Number> & next_ri) const
2436  * {
2437  * {
2438  * TimerOutput::Scope t(timer, "rk_stage - integrals L_h");
2439  *
2440  * for (auto &i : inflow_boundaries)
2441  * i.second->set_time(current_time);
2442  * for (auto &i : subsonic_outflow_boundaries)
2443  * i.second->set_time(current_time);
2444  *
2445  * data.loop(&EulerOperator::local_apply_cell,
2446  * &EulerOperator::local_apply_face,
2447  * &EulerOperator::local_apply_boundary_face,
2448  * this,
2449  * vec_ki,
2450  * current_ri,
2451  * true,
2454  * }
2455  *
2456  *
2457  * {
2458  * TimerOutput::Scope t(timer, "rk_stage - inv mass + vec upd");
2459  * data.cell_loop(
2460  * &EulerOperator::local_apply_inverse_mass_matrix,
2461  * this,
2462  * next_ri,
2463  * vec_ki,
2464  * std::function<void(const unsigned int, const unsigned int)>(),
2465  * [&](const unsigned int start_range, const unsigned int end_range) {
2466  * const Number ai = factor_ai;
2467  * const Number bi = factor_solution;
2468  * if (ai == Number())
2469  * {
2471  * for (unsigned int i = start_range; i < end_range; ++i)
2472  * {
2473  * const Number k_i = next_ri.local_element(i);
2474  * const Number sol_i = solution.local_element(i);
2475  * solution.local_element(i) = sol_i + bi * k_i;
2476  * }
2477  * }
2478  * else
2479  * {
2481  * for (unsigned int i = start_range; i < end_range; ++i)
2482  * {
2483  * const Number k_i = next_ri.local_element(i);
2484  * const Number sol_i = solution.local_element(i);
2485  * solution.local_element(i) = sol_i + bi * k_i;
2486  * next_ri.local_element(i) = sol_i + ai * k_i;
2487  * }
2488  * }
2489  * });
2490  * }
2491  * }
2492  *
2493  *
2494  *
2495  * @endcode
2496  *
2497  * Having discussed the implementation of the functions that deal with
2498  * advancing the solution by one time step, let us now move to functions
2499  * that implement other, ancillary operations. Specifically, these are
2500  * functions that compute projections, evaluate errors, and compute the speed
2501  * of information transport on a cell.
2502  *
2503 
2504  *
2505  * The first of these functions is essentially equivalent to
2506  * VectorTools::project(), just much faster because it is specialized for DG
2507  * elements where there is no need to set up and solve a linear system, as
2508  * each element has independent basis functions. The reason why we show the
2509  * code here, besides a small speedup of this non-critical operation, is that
2510  * it shows additional functionality provided by
2512  *
2513 
2514  *
2515  * The projection operation works as follows: If we denote the matrix of
2516  * shape functions evaluated at quadrature points by @f$S@f$, the projection on
2517  * cell @f$K@f$ is an operation of the form @f$\underbrace{S J^K S^\mathrm
2518  * T}_{\mathcal M^K} \mathbf{w}^K = S J^K
2519  * \tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q}@f$, where @f$J^K@f$ is the diagonal
2520  * matrix containing the determinant of the Jacobian times the quadrature
2521  * weight (JxW), @f$\mathcal M^K@f$ is the cell-wise mass matrix, and
2522  * @f$\tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q}@f$ is the evaluation of the
2523  * field to be projected onto quadrature points. (In reality the matrix @f$S@f$
2524  * has additional structure through the tensor product, as explained in the
2525  * introduction.) This system can now equivalently be written as
2526  * @f$\mathbf{w}^K = \left(S J^K S^\mathrm T\right)^{-1} S J^K
2527  * \tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q} = S^{-\mathrm T}
2528  * \left(J^K\right)^{-1} S^{-1} S J^K
2529  * \tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q}@f$. Now, the term @f$S^{-1} S@f$ and
2530  * then @f$\left(J^K\right)^{-1} J^K@f$ cancel, resulting in the final
2531  * expression @f$\mathbf{w}^K = S^{-\mathrm T}
2532  * \tilde{\mathbf{w}}(\mathbf{x}_q)_{q=1:n_q}@f$. This operation is
2533  * implemented by
2535  * The name is derived from the fact that this projection is simply
2536  * the multiplication by @f$S^{-\mathrm T}@f$, a basis change from the
2537  * nodal basis in the points of the Gaussian quadrature to the given finite
2538  * element basis. Note that we call FEEvaluation::set_dof_values() to write
2539  * the result into the vector, overwriting previous content, rather than
2540  * accumulating the results as typical in integration tasks -- we can do
2541  * this because every vector entry has contributions from only a single
2542  * cell for discontinuous Galerkin discretizations.
2543  *
2544  * @code
2545  * template <int dim, int degree, int n_points_1d>
2546  * void EulerOperator<dim, degree, n_points_1d>::project(
2547  * const Function<dim> & function,
2548  * LinearAlgebra::distributed::Vector<Number> &solution) const
2549  * {
2552  * inverse(phi);
2553  * solution.zero_out_ghosts();
2554  * for (unsigned int cell = 0; cell < data.n_macro_cells(); ++cell)
2555  * {
2556  * phi.reinit(cell);
2557  * for (unsigned int q = 0; q < phi.n_q_points; ++q)
2558  * phi.submit_dof_value(evaluate_function(function,
2559  * phi.quadrature_point(q)),
2560  * q);
2561  * inverse.transform_from_q_points_to_basis(dim + 2,
2562  * phi.begin_dof_values(),
2563  * phi.begin_dof_values());
2564  * phi.set_dof_values(solution);
2565  * }
2566  * }
2567  *
2568  *
2569  *
2570  * @endcode
2571  *
2572  * The next function again repeats functionality also provided by the
2573  * deal.II library, namely VectorTools::integrate_difference(). We here show
2574  * the explicit code to highlight how the vectorization across several cells
2575  * works and how to accumulate results via that interface: Recall that each
2576  * <i>lane</i> of the vectorized array holds data from a different cell. By
2577  * the loop over all cell batches that are owned by the current MPI process,
2578  * we could then fill a VectorizedArray of results; to obtain a global sum,
2579  * we would need to further go on and sum across the entries in the SIMD
2580  * array. However, such a procedure is not stable as the SIMD array could in
2581  * fact not hold valid data for all its lanes. This happens when the number
2582  * of locally owned cells is not a multiple of the SIMD width. To avoid
2583  * invalid data, we must explicitly skip those invalid lanes when accessing
2584  * the data. While one could imagine that we could make it work by simply
2585  * setting the empty lanes to zero (and thus, not contribute to a sum), the
2586  * situation is more complicated than that: What if we were to compute a
2587  * velocity out of the momentum? Then, we would need to divide by the
2588  * density, which is zero -- the result would consequently be NaN and
2589  * contaminate the result. This trap is avoided by accumulating the results
2590  * from the valid SIMD range as we loop through the cell batches, using the
2591  * function MatrixFree::n_active_entries_per_cell_batch() to give us the
2592  * number of lanes with valid data. It equals VectorizedArray::size() on
2593  * most cells, but can be less on the last cell batch if the number of cells
2594  * has a remainder compared to the SIMD width.
2595  *
2596  * @code
2597  * template <int dim, int degree, int n_points_1d>
2598  * std::array<double, 3> EulerOperator<dim, degree, n_points_1d>::compute_errors(
2599  * const Function<dim> & function,
2600  * const LinearAlgebra::distributed::Vector<Number> &solution) const
2601  * {
2602  * TimerOutput::Scope t(timer, "compute errors");
2603  * double errors_squared[3] = {};
2605  *
2606  * for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
2607  * {
2608  * phi.reinit(cell);
2609  * phi.gather_evaluate(solution, true, false);
2610  * VectorizedArray<Number> local_errors_squared[3] = {};
2611  * for (unsigned int q = 0; q < phi.n_q_points; ++q)
2612  * {
2613  * const auto error =
2614  * evaluate_function(function, phi.quadrature_point(q)) -
2615  * phi.get_value(q);
2616  * const auto JxW = phi.JxW(q);
2617  *
2618  * local_errors_squared[0] += error[0] * error[0] * JxW;
2619  * for (unsigned int d = 0; d < dim; ++d)
2620  * local_errors_squared[1] += (error[d + 1] * error[d + 1]) * JxW;
2621  * local_errors_squared[2] += (error[dim + 1] * error[dim + 1]) * JxW;
2622  * }
2623  * for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell);
2624  * ++v)
2625  * for (unsigned int d = 0; d < 3; ++d)
2626  * errors_squared[d] += local_errors_squared[d][v];
2627  * }
2628  *
2629  * Utilities::MPI::sum(errors_squared, MPI_COMM_WORLD, errors_squared);
2630  *
2631  * std::array<double, 3> errors;
2632  * for (unsigned int d = 0; d < 3; ++d)
2633  * errors[d] = std::sqrt(errors_squared[d]);
2634  *
2635  * return errors;
2636  * }
2637  *
2638  *
2639  *
2640  * @endcode
2641  *
2642  * This final function of the EulerOperator class is used to estimate the
2643  * transport speed, scaled by the mesh size, that is relevant for setting
2644  * the time step size in the explicit time integrator. In the Euler
2645  * equations, there are two speeds of transport, namely the convective
2646  * velocity @f$\mathbf{u}@f$ and the propagation of sound waves with sound
2647  * speed @f$c = \sqrt{\gamma p/\rho}@f$ relative to the medium moving at
2648  * velocity @f$\mathbf u@f$.
2649  *
2650 
2651  *
2652  * In the formula for the time step size, we are interested not by
2653  * these absolute speeds, but by the amount of time it takes for
2654  * information to cross a single cell. For information transported along with
2655  * the medium, @f$\mathbf u@f$ is scaled by the mesh size,
2656  * so an estimate of the maximal velocity can be obtained by computing
2657  * @f$\|J^{-\mathrm T} \mathbf{u}\|_\infty@f$, where @f$J@f$ is the Jacobian of the
2658  * transformation from real to the reference domain. Note that
2659  * FEEvaluationBase::inverse_jacobian() returns the inverse and transpose
2660  * Jacobian, representing the metric term from real to reference
2661  * coordinates, so we do not need to transpose it again. We store this limit
2662  * in the variable `convective_limit` in the code below.
2663  *
2664 
2665  *
2666  * The sound propagation is isotropic, so we need to take mesh sizes in any
2667  * direction into account. The appropriate mesh size scaling is then given
2668  * by the minimal singular value of @f$J@f$ or, equivalently, the maximal
2669  * singular value of @f$J^{-1}@f$. Note that one could approximate this quantity
2670  * by the minimal distance between vertices of a cell when ignoring curved
2671  * cells. To get the maximal singular value of the Jacobian, the general
2672  * strategy would be some LAPACK function. Since all we need here is an
2673  * estimate, we can avoid the hassle of decomposing a tensor of
2674  * VectorizedArray numbers into several matrices and go into an (expensive)
2675  * eigenvalue function without vectorization, and instead use a few
2676  * iterations (five in the code below) of the power method applied to
2677  * @f$J^{-1}J^{-\mathrm T}@f$. The speed of convergence of this method depends
2678  * on the ratio of the largest to the next largest eigenvalue and the
2679  * initial guess, which is the vector of all ones. This might suggest that
2680  * we get slow convergence on cells close to a cube shape where all
2681  * lengths are almost the same. However, this slow convergence means that
2682  * the result will sit between the two largest singular values, which both
2683  * are close to the maximal value anyway. In all other cases, convergence
2684  * will be quick. Thus, we can merely hardcode 5 iterations here and be
2685  * confident that the result is good.
2686  *
2687  * @code
2688  * template <int dim, int degree, int n_points_1d>
2689  * double EulerOperator<dim, degree, n_points_1d>::compute_cell_transport_speed(
2690  * const LinearAlgebra::distributed::Vector<Number> &solution) const
2691  * {
2692  * TimerOutput::Scope t(timer, "compute transport speed");
2693  * Number max_transport = 0;
2695  *
2696  * for (unsigned int cell = 0; cell < data.n_cell_batches(); ++cell)
2697  * {
2698  * phi.reinit(cell);
2699  * phi.gather_evaluate(solution, true, false);
2700  * VectorizedArray<Number> local_max = 0.;
2701  * for (unsigned int q = 0; q < phi.n_q_points; ++q)
2702  * {
2703  * const auto solution = phi.get_value(q);
2704  * const auto velocity = euler_velocity<dim>(solution);
2705  * const auto pressure = euler_pressure<dim>(solution);
2706  *
2707  * const auto inverse_jacobian = phi.inverse_jacobian(q);
2708  * const auto convective_speed = inverse_jacobian * velocity;
2709  * VectorizedArray<Number> convective_limit = 0.;
2710  * for (unsigned int d = 0; d < dim; ++d)
2711  * convective_limit =
2712  * std::max(convective_limit, std::abs(convective_speed[d]));
2713  *
2714  * const auto speed_of_sound =
2715  * std::sqrt(gamma * pressure * (1. / solution[0]));
2716  *
2718  * for (unsigned int d = 0; d < dim; ++d)
2719  * eigenvector[d] = 1.;
2720  * for (unsigned int i = 0; i < 5; ++i)
2721  * {
2722  * eigenvector = transpose(inverse_jacobian) *
2723  * (inverse_jacobian * eigenvector);
2724  * VectorizedArray<Number> eigenvector_norm = 0.;
2725  * for (unsigned int d = 0; d < dim; ++d)
2726  * eigenvector_norm =
2727  * std::max(eigenvector_norm, std::abs(eigenvector[d]));
2728  * eigenvector /= eigenvector_norm;
2729  * }
2730  * const auto jac_times_ev = inverse_jacobian * eigenvector;
2731  * const auto max_eigenvalue = std::sqrt(
2732  * (jac_times_ev * jac_times_ev) / (eigenvector * eigenvector));
2733  * local_max =
2734  * std::max(local_max,
2735  * max_eigenvalue * speed_of_sound + convective_limit);
2736  * }
2737  *
2738  * @endcode
2739  *
2740  * Similarly to the previous function, we must make sure to accumulate
2741  * speed only on the valid cells of a cell batch.
2742  *
2743  * @code
2744  * for (unsigned int v = 0; v < data.n_active_entries_per_cell_batch(cell);
2745  * ++v)
2746  * for (unsigned int d = 0; d < 3; ++d)
2747  * max_transport = std::max(max_transport, local_max[v]);
2748  * }
2749  *
2750  * max_transport = Utilities::MPI::max(max_transport, MPI_COMM_WORLD);
2751  *
2752  * return max_transport;
2753  * }
2754  *
2755  *
2756  *
2757  * @endcode
2758  *
2759  *
2760  * <a name="TheEulerProblemclass"></a>
2761  * <h3>The EulerProblem class</h3>
2762  *
2763 
2764  *
2765  * This class combines the EulerOperator class with the time integrator and
2766  * the usual global data structures such as FiniteElement and DoFHandler, to
2767  * actually run the simulations of the Euler problem.
2768  *
2769 
2770  *
2771  * The member variables are a triangulation, a finite element, a mapping (to
2772  * create high-order curved surfaces, see e.g. @ref step_10 "step-10"), and a DoFHandler to
2773  * describe the degrees of freedom. In addition, we keep an instance of the
2774  * EulerOperator described above around, which will do all heavy lifting in
2775  * terms of integrals, and some parameters for time integration like the
2776  * current time or the time step size.
2777  *
2778 
2779  *
2780  * Furthermore, we use a PostProcessor instance to write some additional
2781  * information to the output file, in similarity to what was done in
2782  * @ref step_33 "step-33". The interface of the DataPostprocessor class is intuitive,
2783  * requiring us to provide information about what needs to be evaluated
2784  * (typically only the values of the solution, except for the Schlieren plot
2785  * that we only enable in 2D where it makes sense), and the names of what
2786  * gets evaluated. Note that it would also be possible to extract most
2787  * information by calculator tools within visualization programs such as
2788  * ParaView, but it is so much more convenient to do it already when writing
2789  * the output.
2790  *
2791  * @code
2792  * template <int dim>
2793  * class EulerProblem
2794  * {
2795  * public:
2796  * EulerProblem();
2797  *
2798  * void run();
2799  *
2800  * private:
2801  * void make_grid_and_dofs();
2802  *
2803  * void output_results(const unsigned int result_number);
2804  *
2806  *
2807  * ConditionalOStream pcout;
2808  *
2809  * #ifdef DEAL_II_WITH_P4EST
2811  * #else
2813  * #endif
2814  *
2815  * FESystem<dim> fe;
2816  * MappingQGeneric<dim> mapping;
2817  * DoFHandler<dim> dof_handler;
2818  *
2819  * TimerOutput timer;
2820  *
2821  * EulerOperator<dim, fe_degree, n_q_points_1d> euler_operator;
2822  *
2823  * double time, time_step;
2824  *
2825  * class Postprocessor : public DataPostprocessor<dim>
2826  * {
2827  * public:
2828  * Postprocessor();
2829  *
2830  * virtual void evaluate_vector_field(
2831  * const DataPostprocessorInputs::Vector<dim> &inputs,
2832  * std::vector<Vector<double>> &computed_quantities) const override;
2833  *
2834  * virtual std::vector<std::string> get_names() const override;
2835  *
2836  * virtual std::vector<
2838  * get_data_component_interpretation() const override;
2839  *
2840  * virtual UpdateFlags get_needed_update_flags() const override;
2841  *
2842  * private:
2843  * const bool do_schlieren_plot;
2844  * };
2845  * };
2846  *
2847  *
2848  *
2849  * template <int dim>
2850  * EulerProblem<dim>::Postprocessor::Postprocessor()
2851  * : do_schlieren_plot(dim == 2)
2852  * {}
2853  *
2854  *
2855  *
2856  * @endcode
2857  *
2858  * For the main evaluation of the field variables, we first check that the
2859  * lengths of the arrays equal the expected values (the lengths `2*dim+4` or
2860  * `2*dim+5` are derived from the sizes of the names we specify in the
2861  * get_names() function below). Then we loop over all evaluation points and
2862  * fill the respective information: First we fill the primal solution
2863  * variables of density @f$\rho@f$, momentum @f$\rho \mathbf{u}@f$ and energy @f$E@f$,
2864  * then we compute the derived velocity @f$\mathbf u@f$, the pressure @f$p@f$, the
2865  * speed of sound @f$c=\sqrt{\gamma p / \rho}@f$, as well as the Schlieren plot
2866  * showing @f$s = |\nabla \rho|^2@f$ in case it is enabled. (See @ref step_69 "step-69" for
2867  * another example where we create a Schlieren plot.)
2868  *
2869  * @code
2870  * template <int dim>
2871  * void EulerProblem<dim>::Postprocessor::evaluate_vector_field(
2872  * const DataPostprocessorInputs::Vector<dim> &inputs,
2873  * std::vector<Vector<double>> & computed_quantities) const
2874  * {
2875  * const unsigned int n_evaluation_points = inputs.solution_values.size();
2876  *
2877  * if (do_schlieren_plot == true)
2878  * Assert(inputs.solution_gradients.size() == n_evaluation_points,
2879  * ExcInternalError());
2880  *
2881  * Assert(computed_quantities.size() == n_evaluation_points,
2882  * ExcInternalError());
2883  * Assert(inputs.solution_values[0].size() == dim + 2, ExcInternalError());
2884  * Assert(computed_quantities[0].size() ==
2885  * dim + 2 + (do_schlieren_plot == true ? 1 : 0),
2886  * ExcInternalError());
2887  *
2888  * for (unsigned int q = 0; q < n_evaluation_points; ++q)
2889  * {
2890  * Tensor<1, dim + 2> solution;
2891  * for (unsigned int d = 0; d < dim + 2; ++d)
2892  * solution[d] = inputs.solution_values[q](d);
2893  *
2894  * const double density = solution[0];
2895  * const Tensor<1, dim> velocity = euler_velocity<dim>(solution);
2896  * const double pressure = euler_pressure<dim>(solution);
2897  *
2898  * for (unsigned int d = 0; d < dim; ++d)
2899  * computed_quantities[q](d) = velocity[d];
2900  * computed_quantities[q](dim) = pressure;
2901  * computed_quantities[q](dim + 1) = std::sqrt(gamma * pressure / density);
2902  *
2903  * if (do_schlieren_plot == true)
2904  * computed_quantities[q](dim + 2) =
2905  * inputs.solution_gradients[q][0] * inputs.solution_gradients[q][0];
2906  * }
2907  * }
2908  *
2909  *
2910  *
2911  * template <int dim>
2912  * std::vector<std::string> EulerProblem<dim>::Postprocessor::get_names() const
2913  * {
2914  * std::vector<std::string> names;
2915  * for (unsigned int d = 0; d < dim; ++d)
2916  * names.emplace_back("velocity");
2917  * names.emplace_back("pressure");
2918  * names.emplace_back("speed_of_sound");
2919  *
2920  * if (do_schlieren_plot == true)
2921  * names.emplace_back("schlieren_plot");
2922  *
2923  * return names;
2924  * }
2925  *
2926  *
2927  *
2928  * @endcode
2929  *
2930  * For the interpretation of quantities, we have scalar density, energy,
2931  * pressure, speed of sound, and the Schlieren plot, and vectors for the
2932  * momentum and the velocity.
2933  *
2934  * @code
2935  * template <int dim>
2936  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2937  * EulerProblem<dim>::Postprocessor::get_data_component_interpretation() const
2938  * {
2939  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2940  * interpretation;
2941  * for (unsigned int d = 0; d < dim; ++d)
2942  * interpretation.push_back(
2944  * interpretation.push_back(DataComponentInterpretation::component_is_scalar);
2945  * interpretation.push_back(DataComponentInterpretation::component_is_scalar);
2946  *
2947  * if (do_schlieren_plot == true)
2948  * interpretation.push_back(
2950  *
2951  * return interpretation;
2952  * }
2953  *
2954  *
2955  *
2956  * @endcode
2957  *
2958  * With respect to the necessary update flags, we only need the values for
2959  * all quantities but the Schlieren plot, which is based on the density
2960  * gradient.
2961  *
2962  * @code
2963  * template <int dim>
2964  * UpdateFlags EulerProblem<dim>::Postprocessor::get_needed_update_flags() const
2965  * {
2966  * if (do_schlieren_plot == true)
2967  * return update_values | update_gradients;
2968  * else
2969  * return update_values;
2970  * }
2971  *
2972  *
2973  *
2974  * @endcode
2975  *
2976  * The constructor for this class is unsurprising: We set up a parallel
2977  * triangulation based on the `MPI_COMM_WORLD` communicator, a vector finite
2978  * element with `dim+2` components for density, momentum, and energy, a
2979  * high-order mapping of the same degree as the underlying finite element,
2980  * and initialize the time and time step to zero.
2981  *
2982  * @code
2983  * template <int dim>
2984  * EulerProblem<dim>::EulerProblem()
2985  * : pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
2986  * #ifdef DEAL_II_WITH_P4EST
2987  * , triangulation(MPI_COMM_WORLD)
2988  * #endif
2989  * , fe(FE_DGQ<dim>(fe_degree), dim + 2)
2990  * , mapping(fe_degree)
2991  * , dof_handler(triangulation)
2992  * , timer(pcout, TimerOutput::never, TimerOutput::wall_times)
2993  * , euler_operator(timer)
2994  * , time(0)
2995  * , time_step(0)
2996  * {}
2997  *
2998  *
2999  *
3000  * @endcode
3001  *
3002  * As a mesh, this tutorial program implements two options, depending on the
3003  * global variable `testcase`: For the analytical variant (`testcase==0`),
3004  * the domain is @f$(0, 10) \times (-5, 5)@f$, with Dirichlet boundary
3005  * conditions (inflow) all around the domain. For `testcase==1`, we set the
3006  * domain to a cylinder in a rectangular box, derived from the flow past
3007  * cylinder testcase for incompressible viscous flow by Sch&auml;fer and
3008  * Turek (1996). Here, we have a larger variety of boundaries. The inflow
3009  * part at the left of the channel is given the inflow type, for which we
3010  * choose a constant inflow profile, whereas we set a subsonic outflow at
3011  * the right. For the boundary around the cylinder (boundary id equal to 2)
3012  * as well as the channel walls (boundary id equal to 3) we use the wall
3013  * boundary type, which is no-normal flow. Furthermore, for the 3D cylinder
3014  * we also add a gravity force in vertical direction. Having the base mesh
3015  * in place (including the manifolds set by
3016  * GridGenerator::channel_with_cylinder()), we can then perform the
3017  * specified number of global refinements, create the unknown numbering from
3018  * the DoFHandler, and hand the DoFHandler and Mapping objects to the
3019  * initialization of the EulerOperator.
3020  *
3021  * @code
3022  * template <int dim>
3023  * void EulerProblem<dim>::make_grid_and_dofs()
3024  * {
3025  * switch (testcase)
3026  * {
3027  * case 0:
3028  * {
3029  * Point<dim> lower_left;
3030  * for (unsigned int d = 1; d < dim; ++d)
3031  * lower_left[d] = -5;
3032  *
3033  * Point<dim> upper_right;
3034  * upper_right[0] = 10;
3035  * for (unsigned int d = 1; d < dim; ++d)
3036  * upper_right[d] = 5;
3037  *
3039  * lower_left,
3040  * upper_right);
3041  * triangulation.refine_global(2);
3042  *
3043  * euler_operator.set_inflow_boundary(
3044  * 0, std_cxx14::make_unique<ExactSolution<dim>>(0));
3045  *
3046  * break;
3047  * }
3048  *
3049  * case 1:
3050  * {
3052  * triangulation, 0.03, 1, 0, true);
3053  *
3054  * euler_operator.set_inflow_boundary(
3055  * 0, std_cxx14::make_unique<ExactSolution<dim>>(0));
3056  * euler_operator.set_subsonic_outflow_boundary(
3057  * 1, std_cxx14::make_unique<ExactSolution<dim>>(0));
3058  *
3059  * euler_operator.set_wall_boundary(2);
3060  * euler_operator.set_wall_boundary(3);
3061  *
3062  * if (dim == 3)
3063  * euler_operator.set_body_force(
3064  * std_cxx14::make_unique<Functions::ConstantFunction<dim>>(
3065  * std::vector<double>({0., 0., -0.2})));
3066  *
3067  * break;
3068  * }
3069  *
3070  * default:
3071  * Assert(false, ExcNotImplemented());
3072  * }
3073  *
3074  * triangulation.refine_global(n_global_refinements);
3075  *
3076  * dof_handler.distribute_dofs(fe);
3077  *
3078  * euler_operator.reinit(mapping, dof_handler);
3079  * euler_operator.initialize_vector(solution);
3080  *
3081  * @endcode
3082  *
3083  * In the following, we output some statistics about the problem. Because we
3084  * often end up with quite large numbers of cells or degrees of freedom, we
3085  * would like to print them with a comma to separate each set of three
3086  * digits. This can be done via "locales", although the way this works is
3087  * not particularly intuitive. @ref step_32 "step-32" explains this in slightly more
3088  * detail.
3089  *
3090  * @code
3091  * std::locale s = pcout.get_stream().getloc();
3092  * pcout.get_stream().imbue(std::locale(""));
3093  * pcout << "Number of degrees of freedom: " << dof_handler.n_dofs()
3094  * << " ( = " << (dim + 2) << " [vars] x "
3095  * << triangulation.n_global_active_cells() << " [cells] x "
3096  * << Utilities::pow(fe_degree + 1, dim) << " [dofs/cell/var] )"
3097  * << std::endl;
3098  * pcout.get_stream().imbue(s);
3099  * }
3100  *
3101  *
3102  *
3103  * @endcode
3104  *
3105  * For output, we first let the Euler operator compute the errors of the
3106  * numerical results. More precisely, we compute the error against the
3107  * analytical result for the analytical solution case, whereas we compute
3108  * the deviation against the background field with constant density and
3109  * energy and constant velocity in @f$x@f$ direction for the second test case.
3110  *
3111 
3112  *
3113  * The next step is to create output. This is similar to what is done in
3114  * @ref step_33 "step-33": We let the postprocessor defined above control most of the
3115  * output, except for the primal field that we write directly. For the
3116  * analytical solution test case, we also perform another projection of the
3117  * analytical solution and print the difference between that field and the
3118  * numerical solution. Once we have defined all quantities to be written, we
3119  * build the patches for output. Similarly to @ref step_65 "step-65", we create a
3120  * high-order VTK output by setting the appropriate flag, which enables us
3121  * to visualize fields of high polynomial degrees. Finally, we call the
3122  * `DataOutInterface::write_vtu_in_parallel()` function to write the result
3123  * to the given file name. This function uses special MPI parallel write
3124  * facilities, which are typically more optimized for parallel file systems
3125  * than the standard library's `std::ofstream` variants used in most other
3126  * tutorial programs. A particularly nice feature of the
3127  * `write_vtu_in_parallel()` function is the fact that it can combine output
3128  * from all MPI ranks into a single file, obviating a VTU master file (the
3129  * "pvtu" file).
3130  *
3131 
3132  *
3133  * For parallel programs, it is often instructive to look at the partitioning
3134  * of cells among processors. To this end, one can pass a vector of numbers
3135  * to DataOut::add_data_vector() that contains as many entries as the
3136  * current processor has active cells; these numbers should then be the
3137  * rank of the processor that owns each of these cells. Such a vector
3138  * could, for example, be obtained from
3139  * GridTools::get_subdomain_association(). On the other hand, on each MPI
3140  * process, DataOut will only read those entries that correspond to locally
3141  * owned cells, and these of course all have the same value: namely, the rank
3142  * of the current process. What is in the remaining entries of the vector
3143  * doesn't actually matter, and so we can just get away with a cheap trick: We
3144  * just fill *all* values of the vector we give to DataOut::add_data_vector()
3145  * with the rank of the current MPI process. The key is that on each process,
3146  * only the entries corresponding to the locally owned cells will be read,
3147  * ignoring the (wrong) values in other entries. The fact that every process
3148  * submits a vector in which the correct subset of entries is correct is all
3149  * that is necessary.
3150  *
3151  * @code
3152  * template <int dim>
3153  * void EulerProblem<dim>::output_results(const unsigned int result_number)
3154  * {
3155  * const std::array<double, 3> errors =
3156  * euler_operator.compute_errors(ExactSolution<dim>(time), solution);
3157  * const std::string quantity_name = testcase == 0 ? "error" : "norm";
3158  *
3159  * pcout << "Time:" << std::setw(8) << std::setprecision(3) << time
3160  * << ", dt: " << std::setw(8) << std::setprecision(2) << time_step
3161  * << ", " << quantity_name << " rho: " << std::setprecision(4)
3162  * << std::setw(10) << errors[0] << ", rho * u: " << std::setprecision(4)
3163  * << std::setw(10) << errors[1] << ", energy:" << std::setprecision(4)
3164  * << std::setw(10) << errors[2] << std::endl;
3165  *
3166  * {
3167  * TimerOutput::Scope t(timer, "output");
3168  *
3169  * Postprocessor postprocessor;
3170  * DataOut<dim> data_out;
3171  *
3172  * DataOutBase::VtkFlags flags;
3173  * flags.write_higher_order_cells = true;
3174  * data_out.set_flags(flags);
3175  *
3176  * data_out.attach_dof_handler(dof_handler);
3177  * {
3178  * std::vector<std::string> names;
3179  * names.emplace_back("density");
3180  * for (unsigned int d = 0; d < dim; ++d)
3181  * names.emplace_back("momentum");
3182  * names.emplace_back("energy");
3183  *
3184  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
3185  * interpretation;
3186  * interpretation.push_back(
3188  * for (unsigned int d = 0; d < dim; ++d)
3189  * interpretation.push_back(
3191  * interpretation.push_back(
3193  *
3194  * data_out.add_data_vector(dof_handler, solution, names, interpretation);
3195  * }
3196  * data_out.add_data_vector(solution, postprocessor);
3197  *
3199  * if (testcase == 0 && dim == 2)
3200  * {
3201  * reference.reinit(solution);
3202  * euler_operator.project(ExactSolution<dim>(time), reference);
3203  * reference.sadd(-1., 1, solution);
3204  * std::vector<std::string> names;
3205  * names.emplace_back("error_density");
3206  * for (unsigned int d = 0; d < dim; ++d)
3207  * names.emplace_back("error_momentum");
3208  * names.emplace_back("error_energy");
3209  *
3210  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
3211  * interpretation;
3212  * interpretation.push_back(
3214  * for (unsigned int d = 0; d < dim; ++d)
3215  * interpretation.push_back(
3217  * interpretation.push_back(
3219  *
3220  * data_out.add_data_vector(dof_handler,
3221  * reference,
3222  * names,
3223  * interpretation);
3224  * }
3225  *
3226  * Vector<double> mpi_owner(triangulation.n_active_cells());
3227  * mpi_owner = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
3228  * data_out.add_data_vector(mpi_owner, "owner");
3229  *
3230  * data_out.build_patches(mapping,
3231  * fe.degree,
3233  *
3234  * const std::string filename =
3235  * "solution_" + Utilities::int_to_string(result_number, 3) + ".vtu";
3236  * data_out.write_vtu_in_parallel(filename, MPI_COMM_WORLD);
3237  * }
3238  * }
3239  *
3240  *
3241  *
3242  * @endcode
3243  *
3244  * The EulerProblem::run() function puts all pieces together. It starts off
3245  * by calling the function that creates the mesh and sets up data structures,
3246  * and then initializing the time integrator and the two temporary vectors of
3247  * the low-storage integrator. We call these vectors `rk_register_1` and
3248  * `rk_register_2`, and use the first vector to represent the quantity
3249  * @f$\mathbf{r}_i@f$ and the second one for @f$\mathbf{k}_i@f$ in the formulas for
3250  * the Runge--Kutta scheme outlined in the introduction. Before we start the
3251  * time loop, we compute the time step size by the
3252  * `EulerOperator::compute_cell_transport_speed()` function. For reasons of
3253  * comparison, we compare the result obtained there with the minimal mesh
3254  * size and print them to screen. For velocities and speeds of sound close
3255  * to unity as in this tutorial program, the predicted effective mesh size
3256  * will be close, but they could vary if scaling were different.
3257  *
3258  * @code
3259  * template <int dim>
3260  * void EulerProblem<dim>::run()
3261  * {
3262  * {
3263  * const unsigned int n_vect_number = VectorizedArray<Number>::size();
3264  * const unsigned int n_vect_bits = 8 * sizeof(Number) * n_vect_number;
3265  *
3266  * pcout << "Running with "
3267  * << Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD)
3268  * << " MPI processes" << std::endl;
3269  * pcout << "Vectorization over " << n_vect_number << " "
3270  * << (std::is_same<Number, double>::value ? "doubles" : "floats")
3271  * << " = " << n_vect_bits << " bits ("
3273  * << std::endl;
3274  * }
3275  *
3276  * make_grid_and_dofs();
3277  *
3278  * const LowStorageRungeKuttaIntegrator integrator(lsrk_scheme);
3279  *
3282  * rk_register_1.reinit(solution);
3283  * rk_register_2.reinit(solution);
3284  *
3285  * euler_operator.project(ExactSolution<dim>(time), solution);
3286  *
3287  * double min_vertex_distance = std::numeric_limits<double>::max();
3288  * for (const auto &cell : triangulation.active_cell_iterators())
3289  * if (cell->is_locally_owned())
3290  * min_vertex_distance =
3291  * std::min(min_vertex_distance, cell->minimum_vertex_distance());
3292  * min_vertex_distance =
3293  * Utilities::MPI::min(min_vertex_distance, MPI_COMM_WORLD);
3294  *
3295  * time_step = courant_number * integrator.n_stages() /
3296  * euler_operator.compute_cell_transport_speed(solution);
3297  * pcout << "Time step size: " << time_step
3298  * << ", minimal h: " << min_vertex_distance
3299  * << ", initial transport scaling: "
3300  * << 1. / euler_operator.compute_cell_transport_speed(solution)
3301  * << std::endl
3302  * << std::endl;
3303  *
3304  * output_results(0);
3305  *
3306  * @endcode
3307  *
3308  * Now we are ready to start the time loop, which we run until the time
3309  * has reached the desired end time. Every 5 time steps, we compute a new
3310  * estimate for the time step -- since the solution is nonlinear, it is
3311  * most effective to adapt the value during the course of the
3312  * simulation. In case the Courant number was chosen too aggressively, the
3313  * simulation will typically blow up with time step NaN, so that is easy
3314  * to detect here. One thing to note is that roundoff errors might
3315  * propagate to the leading digits due to an interaction of slightly
3316  * different time step selections that in turn lead to slightly different
3317  * solutions. To decrease this sensitivity, it is common practice to round
3318  * or truncate the time step size to a few digits, e.g. 3 in this case. In
3319  * case the current time is near the prescribed 'tick' value for output
3320  * (e.g. 0.02), we also write the output. After the end of the time loop,
3321  * we summarize the computation by printing some statistics, which is
3322  * mostly done by the TimerOutput::print_wall_time_statistics() function.
3323  *
3324  * @code
3325  * unsigned int timestep_number = 0;
3326  *
3327  * while (time < final_time - 1e-12)
3328  * {
3329  * ++timestep_number;
3330  * if (timestep_number % 5 == 0)
3331  * time_step =
3332  * courant_number * integrator.n_stages() /
3334  * euler_operator.compute_cell_transport_speed(solution), 3);
3335  *
3336  * {
3337  * TimerOutput::Scope t(timer, "rk time stepping total");
3338  * integrator.perform_time_step(euler_operator,
3339  * time,
3340  * time_step,
3341  * solution,
3342  * rk_register_1,
3343  * rk_register_2);
3344  * }
3345  *
3346  * time += time_step;
3347  *
3348  * if (static_cast<int>(time / output_tick) !=
3349  * static_cast<int>((time - time_step) / output_tick) ||
3350  * time >= final_time - 1e-12)
3351  * output_results(
3352  * static_cast<unsigned int>(std::round(time / output_tick)));
3353  * }
3354  *
3355  * timer.print_wall_time_statistics(MPI_COMM_WORLD);
3356  * pcout << std::endl;
3357  * }
3358  *
3359  * } // namespace Euler_DG
3360  *
3361  *
3362  *
3363  * @endcode
3364  *
3365  * The main() function is not surprising and follows what was done in all
3366  * previous MPI programs: As we run an MPI program, we need to call `MPI_Init()`
3367  * and `MPI_Finalize()`, which we do through the
3368  * Utilities::MPI::MPI_InitFinalize data structure. Note that we run the program
3369  * only with MPI, and set the thread count to 1.
3370  *
3371  * @code
3372  * int main(int argc, char **argv)
3373  * {
3374  * using namespace Euler_DG;
3375  * using namespace dealii;
3376  *
3377  * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
3378  *
3379  * try
3380  * {
3381  * deallog.depth_console(0);
3382  *
3383  * EulerProblem<dimension> euler_problem;
3384  * euler_problem.run();
3385  * }
3386  * catch (std::exception &exc)
3387  * {
3388  * std::cerr << std::endl
3389  * << std::endl
3390  * << "----------------------------------------------------"
3391  * << std::endl;
3392  * std::cerr << "Exception on processing: " << std::endl
3393  * << exc.what() << std::endl
3394  * << "Aborting!" << std::endl
3395  * << "----------------------------------------------------"
3396  * << std::endl;
3397  *
3398  * return 1;
3399  * }
3400  * catch (...)
3401  * {
3402  * std::cerr << std::endl
3403  * << std::endl
3404  * << "----------------------------------------------------"
3405  * << std::endl;
3406  * std::cerr << "Unknown exception!" << std::endl
3407  * << "Aborting!" << std::endl
3408  * << "----------------------------------------------------"
3409  * << std::endl;
3410  * return 1;
3411  * }
3412  *
3413  * return 0;
3414  * }
3415  * @endcode
3416 <a name="Results"></a><h1>Results</h1>
3417 
3418 
3419 <a name="Programoutput"></a><h3>Program output</h3>
3420 
3421 
3422 Running the program with the default settings on a machine with 40 processes
3423 produces the following output:
3424 @code
3425 Running with 40 MPI processes
3426 Vectorization over 8 doubles = 512 bits (AVX512)
3427 Number of degrees of freedom: 147,456 ( = 4 [vars] x 1,024 [cells] x 36 [dofs/cell/var] )
3428 Time step size: 0.00689325, minimal h: 0.3125, initial transport scaling: 0.102759
3429 
3430 Time: 0, dt: 0.0069, error rho: 2.76e-07, rho * u: 1.259e-06, energy: 2.987e-06
3431 Time: 1.01, dt: 0.0069, error rho: 1.37e-06, rho * u: 2.252e-06, energy: 4.153e-06
3432 Time: 2.01, dt: 0.0069, error rho: 1.561e-06, rho * u: 2.43e-06, energy: 4.493e-06
3433 Time: 3.01, dt: 0.0069, error rho: 1.714e-06, rho * u: 2.591e-06, energy: 4.762e-06
3434 Time: 4.01, dt: 0.0069, error rho: 1.843e-06, rho * u: 2.625e-06, energy: 4.985e-06
3435 Time: 5.01, dt: 0.0069, error rho: 1.496e-06, rho * u: 1.961e-06, energy: 4.142e-06
3436 Time: 6, dt: 0.0083, error rho: 1.007e-06, rho * u: 7.119e-07, energy: 2.972e-06
3437 Time: 7, dt: 0.0095, error rho: 9.096e-07, rho * u: 3.786e-07, energy: 2.626e-06
3438 Time: 8, dt: 0.0096, error rho: 8.439e-07, rho * u: 3.338e-07, energy: 2.43e-06
3439 Time: 9, dt: 0.0096, error rho: 7.822e-07, rho * u: 2.984e-07, energy: 2.248e-06
3440 Time: 10, dt: 0.0096, error rho: 7.231e-07, rho * u: 2.666e-07, energy: 2.074e-06
3441 
3442 +-------------------------------------------+------------------+------------+------------------+
3443 | Total wallclock time elapsed | 2.249s 30 | 2.249s | 2.249s 8 |
3444 | | | |
3445 | Section | no. calls | min time rank | avg time | max time rank |
3446 +-------------------------------------------+------------------+------------+------------------+
3447 | compute errors | 11 | 0.008066s 13 | 0.00952s | 0.01041s 20 |
3448 | compute transport speed | 258 | 0.01012s 13 | 0.05392s | 0.08574s 25 |
3449 | output | 11 | 0.9597s 13 | 0.9613s | 0.9623s 6 |
3450 | rk time stepping total | 1283 | 0.9827s 25 | 1.015s | 1.06s 13 |
3451 | rk_stage - integrals L_h | 6415 | 0.8803s 26 | 0.9198s | 0.9619s 14 |
3452 | rk_stage - inv mass + vec upd | 6415 | 0.05677s 15 | 0.06487s | 0.07597s 13 |
3453 +-------------------------------------------+------------------+------------+------------------+
3454 @endcode
3455 
3456 The program output shows that all errors are small. This is due to the fact
3457 that we use a relatively fine mesh of @f$32^2@f$ cells with polynomials of degree
3458 5 for a solution that is smooth. An interesting pattern shows for the time
3459 step size: whereas it is 0.0069 up to time 5, it increases to 0.0096 for later
3460 times. The step size increases once the vortex with some motion on top of the
3461 speed of sound (and thus faster propagation) leaves the computational domain
3462 between times 5 and 6.5. After that point, the flow is simply uniform
3463 in the same direction, and the maximum velocity of the gas is reduced
3464 compared to the previous state where the uniform velocity was overlaid
3465 by the vortex. Our time step formula recognizes this effect.
3466 
3467 The final block of output shows detailed information about the timing
3468 of individual parts of the programs; it breaks this down by showing
3469 the time taken by the fastest and the slowest processor, and the
3470 average time -- this is often useful in very large computations to
3471 find whether there are processors that are consistently overheated
3472 (and consequently are throttling their clock speed) or consistently
3473 slow for other reasons.
3474 The summary shows that 1283 time steps have been performed
3475 in 1.02 seconds (looking at the average time among all MPI processes), while
3476 the output of 11 files has taken additional 0.96 seconds. Broken down per time
3477 step and into the five Runge--Kutta stages, the compute time per evaluation is
3478 0.16 milliseconds. This high performance is typical of matrix-free evaluators
3479 and a reason why explicit time integration is very competitive against
3480 implicit solvers, especially for large-scale simulations. The breakdown of
3481 computational times at the end of the program run shows that the evaluation of
3482 integrals in @f$\mathcal L_h@f$ contributes with around 0.92 seconds and the
3483 application of the inverse mass matrix with 0.06 seconds. Furthermore, the
3484 estimation of the transport speed for the time step size computation
3485 contributes with another 0.05 seconds of compute time.
3486 
3487 If we use three more levels of global refinement and 9.4 million DoFs in total,
3488 the final statistics are as follows (for the modified Lax--Friedrichs flux,
3489 @f$p=5@f$, and the same system of 40 cores of dual-socket Intel Xeon Gold 6230):
3490 @code
3491 +-------------------------------------------+------------------+------------+------------------+
3492 | Total wallclock time elapsed | 244.9s 12 | 244.9s | 244.9s 34 |
3493 | | | |
3494 | Section | no. calls | min time rank | avg time | max time rank |
3495 +-------------------------------------------+------------------+------------+------------------+
3496 | compute errors | 11 | 0.4239s 12 | 0.4318s | 0.4408s 9 |
3497 | compute transport speed | 2053 | 3.962s 12 | 6.727s | 10.12s 7 |
3498 | output | 11 | 30.35s 12 | 30.36s | 30.37s 9 |
3499 | rk time stepping total | 10258 | 201.7s 7 | 205.1s | 207.8s 12 |
3500 | rk_stage - integrals L_h | 51290 | 121.3s 6 | 126.6s | 136.3s 16 |
3501 | rk_stage - inv mass + vec upd | 51290 | 66.19s 16 | 77.52s | 81.84s 10 |
3502 +-------------------------------------------+------------------+------------+------------------+
3503 @endcode
3504 
3505 Per time step, the solver now takes 0.02 seconds, about 25 times as long as
3506 for the small problem with 147k unknowns. Given that the problem involves 64
3507 times as many unknowns, the increase in computing time is not
3508 surprising. Since we also do 8 times as many time steps, the compute time
3509 should in theory increase by a factor of 512. The actual increase is 205 s /
3510 1.02 s = 202. This is because the small problem size cannot fully utilize the
3511 40 cores due to communication overhead. This becomes clear if we look into the
3512 details of the operations done per time step. The evaluation of the
3513 differential operator @f$\mathcal L_h@f$ with nearest neighbor communication goes
3514 from 0.92 seconds to 127 seconds, i.e., it increases with a factor of 138. On
3515 the other hand, the cost for application of the inverse mass matrix and the
3516 vector updates, which do not need to communicate between the MPI processes at
3517 all, has increased by a factor of 1195. The increase is more than the
3518 theoretical factor of 512 because the operation is limited by the bandwidth
3519 from RAM memory for the larger size while for the smaller size, all vectors
3520 fit into the caches of the CPU. The numbers show that the mass matrix
3521 evaluation and vector update part consume almost 40% of the time spent by the
3522 Runge--Kutta stages -- despite using a low-storage Runge--Kutta integrator and
3523 merging of vector operations! And despite using over-integration for the
3524 @f$\mathcal L_h@f$ operator. For simpler differential operators and more expensive
3525 time integrators, the proportion spent in the mass matrix and vector update
3526 part can also reach 70%. If we compute a throughput number in terms of DoFs
3527 processed per second and Runge--Kutta stage, we obtain @f[ \text{throughput} =
3528 \frac{n_\mathrm{time steps} n_\mathrm{stages}
3529 n_\mathrm{dofs}}{t_\mathrm{compute}} = \frac{10258 \cdot 5 \cdot
3530 9.4\,\text{MDoFs}}{205s} = 2360\, \text{MDoFs/s} @f] This throughput number is
3531 very high, given that simply copying one vector to another one runs at
3532 only around 10,000 MDoFs/s.
3533 
3534 If we go to the next-larger size with 37.7 million DoFs, the overall
3535 simulation time is 2196 seconds, with 1978 seconds spent in the time
3536 stepping. The increase in run time is a factor of 9.3 for the L_h operator
3537 (1179 versus 127 seconds) and a factor of 10.3 for the inverse mass matrix and
3538 vector updates (797 vs 77.5 seconds). The reason for this non-optimal increase
3539 in run time can be traced back to cache effects on the given hardware (with 40
3540 MB of L2 cache and 55 MB of L3 cache): While not all of the relevant data fits
3541 into caches for 9.4 million DoFs (one vector takes 75 MB and we have three
3542 vectors plus some additional data in MatrixFree), there is capacity for one and
3543 a half vector nonetheless. Given that modern caches are more sophisticated than
3544 the naive least-recently-used strategy (where we would have little re-use as
3545 the data is used in a streaming-like fashion), we can assume that a sizeable
3546 fraction of data can indeed be delivered from caches for the 9.4 million DoFs
3547 case. For the larger case, even with optimal caching less than 10 percent of
3548 data would fit into caches, with an associated loss in performance.
3549 
3550 
3551 <a name="Convergenceratesfortheanalyticaltestcase"></a><h3>Convergence rates for the analytical test case</h3>
3552 
3553 
3554 For the modified Lax--Friedrichs flux and measuring the error in the momentum
3555 variable, we obtain the following convergence table (the rates are very
3556 similar for the density and energy variables):
3557 
3558 <table align="center" class="doxtable">
3559  <tr>
3560  <th>&nbsp;</th>
3561  <th colspan="3"><i>p</i>=2</th>
3562  <th colspan="3"><i>p</i>=3</th>
3563  <th colspan="3"><i>p</i>=5</th>
3564  </tr>
3565  <tr>
3566  <th>n_cells</th>
3567  <th>n_dofs</th>
3568  <th>error mom</th>
3569  <th>rate</th>
3570  <th>n_dofs</th>
3571  <th>error mom</th>
3572  <th>rate</th>
3573  <th>n_dofs</th>
3574  <th>error mom</th>
3575  <th>rate</th>
3576  </tr>
3577  <tr>
3578  <td align="right">16</td>
3579  <td>&nbsp;</td>
3580  <td>&nbsp;</td>
3581  <td>&nbsp;</td>
3582  <td>&nbsp;</td>
3583  <td>&nbsp;</td>
3584  <td>&nbsp;</td>
3585  <td align="right">2,304</td>
3586  <td align="center">1.373e-01</td>
3587  <td>&nbsp;</td>
3588  </tr>
3589  <tr>
3590  <td align="right">64</td>
3591  <td>&nbsp;</td>
3592  <td>&nbsp;</td>
3593  <td>&nbsp;</td>
3594  <td align="right">4,096</td>
3595  <td align="center">9.130e-02</td>
3596  <td>&nbsp;</td>
3597  <td align="right">9,216</td>
3598  <td align="center">8.899e-03</td>
3599  <td>3.94</td>
3600  </tr>
3601  <tr>
3602  <td align="right">256</td>
3603  <td align="right">9,216</td>
3604  <td align="center">5.577e-02</td>
3605  <td>&nbsp;</td>
3606  <td align="right">16,384</td>
3607  <td align="center">7.381e-03</td>
3608  <td>3.64</td>
3609  <td align="right">36,864</td>
3610  <td align="center">2.082e-04</td>
3611  <td>5.42</td>
3612  </tr>
3613  <tr>
3614  <td align="right">1024</td>
3615  <td align="right">36,864</td>
3616  <td align="center">4.724e-03</td>
3617  <td>3.56</td>
3618  <td align="right">65,536</td>
3619  <td align="center">3.072e-04</td>
3620  <td>4.59</td>
3621  <td align="right">147,456</td>
3622  <td align="center">2.625e-06</td>
3623  <td>6.31</td>
3624  </tr>
3625  <tr>
3626  <td align="right">4096</td>
3627  <td align="right">147,456</td>
3628  <td align="center">6.205e-04</td>
3629  <td>2.92</td>
3630  <td align="right">262,144</td>
3631  <td align="center">1.880e-05</td>
3632  <td>4.03</td>
3633  <td align="right">589,824</td>
3634  <td align="center">3.268e-08</td>
3635  <td>6.33</td>
3636  </tr>
3637  <tr>
3638  <td align="right">16,384</td>
3639  <td align="right">589,824</td>
3640  <td align="center">8.279e-05</td>
3641  <td>2.91</td>
3642  <td align="right">1,048,576</td>
3643  <td align="center">1.224e-06</td>
3644  <td>3.94</td>
3645  <td align="right">2,359,296</td>
3646  <td align="center">9.252e-10</td>
3647  <td>5.14</td>
3648  </tr>
3649  <tr>
3650  <td align="right">65,536</td>
3651  <td align="right">2,359,296</td>
3652  <td align="center">1.105e-05</td>
3653  <td>2.91</td>
3654  <td align="right">4,194,304</td>
3655  <td align="center">7.871e-08</td>
3656  <td>3.96</td>
3657  <td align="right">9,437,184</td>
3658  <td align="center">1.369e-10</td>
3659  <td>2.77</td>
3660  </tr>
3661  <tr>
3662  <td align="right">262,144</td>
3663  <td align="right">9,437,184</td>
3664  <td align="center">1.615e-06</td>
3665  <td>2.77</td>
3666  <td align="right">16,777,216</td>
3667  <td align="center">4.961e-09</td>
3668  <td>3.99</td>
3669  <td align="right">37,748,736</td>
3670  <td align="center">7.091e-11</td>
3671  <td>0.95</td>
3672  </tr>
3673 </table>
3674 
3675 If we switch to the Harten-Lax-van Leer flux, the results are as follows:
3676 <table align="center" class="doxtable">
3677  <tr>
3678  <th>&nbsp;</th>
3679  <th colspan="3"><i>p</i>=2</th>
3680  <th colspan="3"><i>p</i>=3</th>
3681  <th colspan="3"><i>p</i>=5</th>
3682  </tr>
3683  <tr>
3684  <th>n_cells</th>
3685  <th>n_dofs</th>
3686  <th>error mom</th>
3687  <th>rate</th>
3688  <th>n_dofs</th>
3689  <th>error mom</th>
3690  <th>rate</th>
3691  <th>n_dofs</th>
3692  <th>error mom</th>
3693  <th>rate</th>
3694  </tr>
3695  <tr>
3696  <td align="right">16</td>
3697  <td>&nbsp;</td>
3698  <td>&nbsp;</td>
3699  <td>&nbsp;</td>
3700  <td>&nbsp;</td>
3701  <td>&nbsp;</td>
3702  <td>&nbsp;</td>
3703  <td align="right">2,304</td>
3704  <td align="center">1.339e-01</td>
3705  <td>&nbsp;</td>
3706  </tr>
3707  <tr>
3708  <td align="right">64</td>
3709  <td>&nbsp;</td>
3710  <td>&nbsp;</td>
3711  <td>&nbsp;</td>
3712  <td align="right">4,096</td>
3713  <td align="center">9.037e-02</td>
3714  <td>&nbsp;</td>
3715  <td align="right">9,216</td>
3716  <td align="center">8.849e-03</td>
3717  <td>3.92</td>
3718  </tr>
3719  <tr>
3720  <td align="right">256</td>
3721  <td align="right">9,216</td>
3722  <td align="center">4.204e-02</td>
3723  <td>&nbsp;</td>
3724  <td align="right">16,384</td>
3725  <td align="center">9.143e-03</td>
3726  <td>3.31</td>
3727  <td align="right">36,864</td>
3728  <td align="center">2.501e-04</td>
3729  <td>5.14</td>
3730  </tr>
3731  <tr>
3732  <td align="right">1024</td>
3733  <td align="right">36,864</td>
3734  <td align="center">4.913e-03</td>
3735  <td>3.09</td>
3736  <td align="right">65,536</td>
3737  <td align="center">3.257e-04</td>
3738  <td>4.81</td>
3739  <td align="right">147,456</td>
3740  <td align="center">3.260e-06</td>
3741  <td>6.26</td>
3742  </tr>
3743  <tr>
3744  <td align="right">4096</td>
3745  <td align="right">147,456</td>
3746  <td align="center">7.862e-04</td>
3747  <td>2.64</td>
3748  <td align="right">262,144</td>
3749  <td align="center">1.588e-05</td>
3750  <td>4.36</td>
3751  <td align="right">589,824</td>
3752  <td align="center">2.953e-08</td>
3753  <td>6.79</td>
3754  </tr>
3755  <tr>
3756  <td align="right">16,384</td>
3757  <td align="right">589,824</td>
3758  <td align="center">1.137e-04</td>
3759  <td>2.79</td>
3760  <td align="right">1,048,576</td>
3761  <td align="center">9.400e-07</td>
3762  <td>4.08</td>
3763  <td align="right">2,359,296</td>
3764  <td align="center">4.286e-10</td>
3765  <td>6.11</td>
3766  </tr>
3767  <tr>
3768  <td align="right">65,536</td>
3769  <td align="right">2,359,296</td>
3770  <td align="center">1.476e-05</td>
3771  <td>2.95</td>
3772  <td align="right">4,194,304</td>
3773  <td align="center">5.799e-08</td>
3774  <td>4.02</td>
3775  <td align="right">9,437,184</td>
3776  <td align="center">2.789e-11</td>
3777  <td>3.94</td>
3778  </tr>
3779  <tr>
3780  <td align="right">262,144</td>
3781  <td align="right">9,437,184</td>
3782  <td align="center">2.038e-06</td>
3783  <td>2.86</td>
3784  <td align="right">16,777,216</td>
3785  <td align="center">3.609e-09</td>
3786  <td>4.01</td>
3787  <td align="right">37,748,736</td>
3788  <td align="center">5.730e-11</td>
3789  <td>-1.04</td>
3790  </tr>
3791 </table>
3792 
3793 The tables show that we get optimal @f$\mathcal O\left(h^{p+1}\right)@f$
3794 convergence rates for both numerical fluxes. The errors are slightly smaller
3795 for the Lax--Friedrichs flux for @f$p=2@f$, but the picture is reversed for
3796 @f$p=3@f$; in any case, the differences on this testcase are relatively
3797 small.
3798 
3799 For @f$p=5@f$, we reach the roundoff accuracy of @f$10^{-11}@f$ with both
3800 fluxes on the finest grids. Also note that the errors are absolute with a
3801 domain length of @f$10^2@f$, so relative errors are below @f$10^{-12}@f$. The HLL flux
3802 is somewhat better for the highest degree, which is due to a slight inaccuracy
3803 of the Lax--Friedrichs flux: The Lax--Friedrichs flux sets a Dirichlet
3804 condition on the solution that leaves the domain, which results in a small
3805 artificial reflection, which is accentuated for the Lax--Friedrichs
3806 flux. Apart from that, we see that the influence of the numerical flux is
3807 minor, as the polynomial part inside elements is the main driver of the
3808 accucary. The limited influence of the flux also has consequences when trying
3809 to approach more challenging setups with the higher-order DG setup: Taking for
3810 example the parameters and grid of @ref step_33 "step-33", we get oscillations (which in turn
3811 make density negative and make the solution explode) with both fluxes once the
3812 high-mass part comes near the boundary, as opposed to the low-order finite
3813 volume case (@f$p=0@f$). Thus, any case that leads to shocks in the solution
3814 necessitates some form of limiting or artificial dissipation. For another
3815 alternative, see the @ref step_69 "step-69" tutorial program.
3816 
3817 
3818 <a name="Resultsforflowinchannelaroundcylinderin2D"></a><h3>Results for flow in channel around cylinder in 2D</h3>
3819 
3820 
3821 For the test case of the flow around a cylinder in a channel, we need to
3822 change the first code line to
3823 @code
3824  constexpr unsigned int testcase = 1;
3825 @endcode
3826 This test case starts with a background field of a constant velocity
3827 of Mach number 0.31 and a constant initial density; the flow will have
3828 to go around an obstacle in the form of a cylinder. Since we impose a
3829 no-penetration condition on the cylinder walls, the flow that
3830 initially impinges head-on onto to cylinder has to rearrange,
3831 which creates a big sound wave. The following pictures show the pressure at
3832 times 0.1, 0.25, 0.5, and 1.0 (top left to bottom right) for the 2D case with
3833 5 levels of global refinement, using 102,400 cells with polynomial degree of
3834 5 and 14.7 million degrees of freedom over all 4 solution variables.
3835 We clearly see the discontinuity that
3836 propagates slowly in the upstream direction and more quickly in downstream
3837 direction in the first snapshot at time 0.1. At time 0.25, the sound wave has
3838 reached the top and bottom walls and reflected back to the interior. From the
3839 different distances of the reflected waves from lower and upper walls we can
3840 see the slight asymmetry of the Sch&auml;fer-Turek test case represented by
3841 GridGenerator::channel_with_cylinder() with somewhat more space above the
3842 cylinder compared to below. At later times, the picture is more chaotic with
3843 many sound waves all over the place.
3844 
3845 <table align="center" class="doxtable" style="width:85%">
3846  <tr>
3847  <td>
3848  <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_010.png" alt="" width="100%">
3849  </td>
3850  <td>
3851  <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_025.png" alt="" width="100%">
3852  </td>
3853  </tr>
3854  <tr>
3855  <td>
3856  <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_050.png" alt="" width="100%">
3857  </td>
3858  <td>
3859  <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_100.png" alt="" width="100%">
3860  </td>
3861  </tr>
3862 </table>
3863 
3864 The next picture shows an elevation plot of the pressure at time 1.0 looking
3865 from the channel inlet towards the outlet at the same resolution -- here,
3866 we can see the large number
3867 of reflections. In the figure, two types of waves are visible. The
3868 larger-amplitude waves correspond to various reflections that happened as the
3869 initial discontinuity hit the walls, whereas the small-amplitude waves of
3870 size similar to the elements correspond to numerical artifacts. They have their
3871 origin in the finite resolution of the scheme and appear as the discontinuity
3872 travels through elements with high-order polynomials. This effect can be cured
3873 by increasing resolution. Apart from this effect, the rich wave structure is
3874 the result of the transport accuracy of the high-order DG method.
3875 
3876 <img src="https://www.dealii.org/images/steps/developer/step-67.pressure_elevated.jpg" alt="" width="40%">
3877 
3878 With 2 levels of global refinement with 1,600 cells, the mesh and its
3879 partitioning on 40 MPI processes looks as follows:
3880 
3881 <img src="https://www.dealii.org/images/steps/developer/step-67.grid-owner.png" alt="" width="70%">
3882 
3883 When we run the code with 4 levels of global refinements on 40 cores, we get
3884 the following output:
3885 @code
3886 Running with 40 MPI processes
3887 Vectorization over 8 doubles = 512 bits (AVX512)
3888 Number of degrees of freedom: 3,686,400 ( = 4 [vars] x 25,600 [cells] x 36 [dofs/cell/var] )
3889 Time step size: 7.39876e-05, minimal h: 0.001875, initial transport scaling: 0.00110294
3890 
3891 Time: 0, dt: 7.4e-05, norm rho: 4.17e-16, rho * u: 1.629e-16, energy: 1.381e-15
3892 Time: 0.05, dt: 6.3e-05, norm rho: 0.02075, rho * u: 0.03801, energy: 0.08772
3893 Time: 0.1, dt: 5.9e-05, norm rho: 0.02211, rho * u: 0.04515, energy: 0.08953
3894 Time: 0.15, dt: 5.7e-05, norm rho: 0.02261, rho * u: 0.04592, energy: 0.08967
3895 Time: 0.2, dt: 5.8e-05, norm rho: 0.02058, rho * u: 0.04361, energy: 0.08222
3896 Time: 0.25, dt: 5.9e-05, norm rho: 0.01695, rho * u: 0.04203, energy: 0.06873
3897 Time: 0.3, dt: 5.9e-05, norm rho: 0.01653, rho * u: 0.0401, energy: 0.06604
3898 Time: 0.35, dt: 5.7e-05, norm rho: 0.01774, rho * u: 0.04264, energy: 0.0706
3899 
3900 ...
3901 
3902 Time: 1.95, dt: 5.8e-05, norm rho: 0.01488, rho * u: 0.03923, energy: 0.05185
3903 Time: 2, dt: 5.7e-05, norm rho: 0.01432, rho * u: 0.03969, energy: 0.04889
3904 
3905 +-------------------------------------------+------------------+------------+------------------+
3906 | Total wallclock time elapsed | 273.6s 13 | 273.6s | 273.6s 0 |
3907 | | | |
3908 | Section | no. calls | min time rank | avg time | max time rank |
3909 +-------------------------------------------+------------------+------------+------------------+
3910 | compute errors | 41 | 0.01112s 35 | 0.0672s | 0.1337s 0 |
3911 | compute transport speed | 6914 | 5.422s 35 | 15.96s | 29.99s 1 |
3912 | output | 41 | 37.24s 35 | 37.3s | 37.37s 0 |
3913 | rk time stepping total | 34564 | 205.4s 1 | 219.5s | 230.1s 35 |
3914 | rk_stage - integrals L_h | 172820 | 153.6s 1 | 164.9s | 175.6s 27 |
3915 | rk_stage - inv mass + vec upd | 172820 | 47.13s 13 | 53.09s | 64.05s 33 |
3916 +-------------------------------------------+------------------+------------+------------------+
3917 @endcode
3918 
3919 The norms shown here for the various quantities are the deviations
3920 @f$\rho'@f$, @f$(\rho u)'@f$, and @f$E'@f$ against the background field (namely, the
3921 initial condition). The distribution of run time is overall similar as in the
3922 previous test case. The only slight difference is the larger proportion of
3923 time spent in @f$\mathcal L_h@f$ as compared to the inverse mass matrix and vector
3924 updates. This is because the geometry is deformed and the matrix-free
3925 framework needs to load additional arrays for the geometry from memory that
3926 are compressed in the affine mesh case.
3927 
3928 Increasing the number of global refinements to 5, the output becomes:
3929 @code
3930 Running with 40 MPI processes
3931 Vectorization over 8 doubles = 512 bits (AVX512)
3932 Number of degrees of freedom: 14,745,600 ( = 4 [vars] x 102,400 [cells] x 36 [dofs/cell/var] )
3933 
3934 ...
3935 
3936 +-------------------------------------------+------------------+------------+------------------+
3937 | Total wallclock time elapsed | 2693s 32 | 2693s | 2693s 23 |
3938 | | | |
3939 | Section | no. calls | min time rank | avg time | max time rank |
3940 +-------------------------------------------+------------------+------------+------------------+
3941 | compute errors | 41 | 0.04537s 32 | 0.173s | 0.3489s 0 |
3942 | compute transport speed | 13858 | 40.75s 32 | 85.99s | 149.8s 0 |
3943 | output | 41 | 153.8s 32 | 153.9s | 154.1s 0 |
3944 | rk time stepping total | 69284 | 2386s 0 | 2450s | 2496s 32 |
3945 | rk_stage - integrals L_h | 346420 | 1365s 32 | 1574s | 1718s 19 |
3946 | rk_stage - inv mass + vec upd | 346420 | 722.5s 10 | 870.7s | 1125s 32 |
3947 +-------------------------------------------+------------------+------------+------------------+
3948 @endcode
3949 
3950 The effect on performance is similar to the analytical test case -- in
3951 theory, computation times should increase by a factor of 8, but we actually
3952 see an increase by a factor of 11 for the time steps (219.5 seconds versus
3953 2450 seconds). This can be traced back to caches, with the small case mostly
3954 fitting in caches. An interesting effect, typical of programs with a mix of
3955 local communication (integrals @f$\mathcal L_h@f$) and global communication (computation of
3956 transport speed) with some load imbalance, can be observed by looking at the
3957 MPI ranks that encounter the minimal and maximal time of different phases,
3958 respectively. Rank 0 reports the fastest throughput for the "rk time stepping
3959 total" part. At the same time, it appears to be slowest for the "compute
3960 transport speed" part, almost a factor of 2 slower than the
3961 average and almost a factor of 4 compared to the faster rank.
3962 Since the latter involves global communication, we can attribute the
3963 slowness in this part to the fact that the local Runge--Kutta stages have
3964 advanced more quickly on this rank and need to wait until the other processors
3965 catch up. At this point, one can wonder about the reason for this imbalance:
3966 The number of cells is almost the same on all MPI processes.
3967 However, the matrix-free framework is faster on affine and Cartesian
3968 cells located towards the outlet of the channel, to which the lower MPI ranks
3969 are assigned. On the other hand, rank 32, which reports the highest run time
3970 for the Runga--Kutta stages, owns the curved cells near the cylinder, for
3971 which no data compression is possible. To improve throughput, we could assign
3972 different weights to different cell types when partitioning the
3973 parallel::distributed::Triangulation object, or even measure the run time for a
3974 few time steps and try to rebalance then.
3975 
3976 The throughput per Runge--Kutta stage can be computed to 2085 MDoFs/s for the
3977 14.7 million DoFs test case over the 346,000 Runge--Kutta stages, slightly slower
3978 than the Cartesian mesh throughput of 2360 MDoFs/s reported above.
3979 
3980 Finally, if we add one additional refinement, we record the following output:
3981 @code
3982 Running with 40 MPI processes
3983 Vectorization over 8 doubles = 512 bits (AVX512)
3984 Number of degrees of freedom: 58,982,400 ( = 4 [vars] x 409,600 [cells] x 36 [dofs/cell/var] )
3985 
3986 ...
3987 
3988 Time: 1.95, dt: 1.4e-05, norm rho: 0.01488, rho * u: 0.03923, energy: 0.05183
3989 Time: 2, dt: 1.4e-05, norm rho: 0.01431, rho * u: 0.03969, energy: 0.04887
3990 
3991 +-------------------------------------------+------------------+------------+------------------+
3992 | Total wallclock time elapsed | 2.166e+04s 26 | 2.166e+04s | 2.166e+04s 24 |
3993 | | | |
3994 | Section | no. calls | min time rank | avg time | max time rank |
3995 +-------------------------------------------+------------------+------------+------------------+
3996 | compute errors | 41 | 0.1758s 30 | 0.672s | 1.376s 1 |
3997 | compute transport speed | 27748 | 321.3s 34 | 678.8s | 1202s 1 |
3998 | output | 41 | 616.3s 32 | 616.4s | 616.4s 34 |
3999 | rk time stepping total | 138733 | 1.983e+04s 1 | 2.036e+04s | 2.072e+04s 34 |
4000 | rk_stage - integrals L_h | 693665 | 1.052e+04s 32 | 1.248e+04s | 1.387e+04s 19 |
4001 | rk_stage - inv mass + vec upd | 693665 | 6404s 10 | 7868s | 1.018e+04s 32 |
4002 +-------------------------------------------+------------------+------------+------------------+
4003 @endcode
4004 
4005 The "rk time stepping total" part corresponds to a throughput of 2010 MDoFs/s. The
4006 overall run time to perform 139k time steps is 20k seconds (5.7 hours) or 7
4007 time steps per second -- not so bad for having nearly 60 million
4008 unknowns. More throughput can be achieved by adding more cores to
4009 the computation.
4010 
4011 
4012 <a name="Resultsforflowinchannelaroundcylinderin3D"></a><h3>Results for flow in channel around cylinder in 3D</h3>
4013 
4014 
4015 Switching the channel test case to 3D with 3 global refinements, the output is
4016 @code
4017 Running with 40 MPI processes
4018 Vectorization over 8 doubles = 512 bits (AVX512)
4019 Number of degrees of freedom: 221,184,000 ( = 5 [vars] x 204,800 [cells] x 216 [dofs/cell/var] )
4020 
4021 ...
4022 
4023 Time: 1.95, dt: 0.00011, norm rho: 0.01131, rho * u: 0.03056, energy: 0.04091
4024 Time: 2, dt: 0.00011, norm rho: 0.0119, rho * u: 0.03142, energy: 0.04425
4025 
4026 +-------------------------------------------+------------------+------------+------------------+
4027 | Total wallclock time elapsed | 1.734e+04s 4 | 1.734e+04s | 1.734e+04s 38 |
4028 | | | |
4029 | Section | no. calls | min time rank | avg time | max time rank |
4030 +-------------------------------------------+------------------+------------+------------------+
4031 | compute errors | 41 | 0.6551s 34 | 3.216s | 7.281s 0 |
4032 | compute transport speed | 3546 | 160s 34 | 393.2s | 776.9s 0 |
4033 | output | 41 | 1350s 34 | 1353s | 1357s 0 |
4034 | rk time stepping total | 17723 | 1.519e+04s 0 | 1.558e+04s | 1.582e+04s 34 |
4035 | rk_stage - integrals L_h | 88615 | 1.005e+04s 32 | 1.126e+04s | 1.23e+04s 11 |
4036 | rk_stage - inv mass + vec upd | 88615 | 3056s 11 | 4322s | 5759s 32 |
4037 +-------------------------------------------+------------------+------------+------------------+
4038 @endcode
4039 
4040 The physics are similar to the 2D case, with a slight motion in the z
4041 direction due to the gravitational force. The throughput per Runge--Kutta
4042 stage in this case is
4043 @f[
4044 \text{throughput} = \frac{n_\mathrm{time steps} n_\mathrm{stages}
4045 n_\mathrm{dofs}}{t_\mathrm{compute}} =
4046 \frac{17723 \cdot 5 \cdot 221.2\,\text{M}}{15580s} = 1258\, \text{MDoFs/s}.
4047 @f]
4048 
4049 The throughput is lower than in 2D because the computation of the @f$\mathcal L_h@f$ term
4050 is more expensive. This is due to over-integration with `degree+2` points and
4051 the larger fraction of face integrals (worse volume-to-surface ratio) with
4052 more expensive flux computations. If we only consider the inverse mass matrix
4053 and vector update part, we record a throughput of 4857 MDoFs/s for the 2D case
4054 of the isentropic vortex with 37.7 million unknowns, whereas the 3D case
4055 runs with 4535 MDoFs/s. The performance is similar because both cases are in
4056 fact limited by the memory bandwidth.
4057 
4058 If we go to four levels of global refinement, we need to increase the number
4059 of processes to fit everything in memory -- the computation needs around 350
4060 GB of RAM memory in this case. Also, the time it takes to complete 35k time
4061 steps becomes more tolerable by adding additional resources. We therefore use
4062 6 nodes with 40 cores each, resulting in a computation with 240 MPI processes:
4063 @code
4064 Running with 240 MPI processes
4065 Vectorization over 8 doubles = 512 bits (AVX512)
4066 Number of degrees of freedom: 1,769,472,000 ( = 5 [vars] x 1,638,400 [cells] x 216 [dofs/cell/var] )
4067 
4068 ...
4069 
4070 Time: 1.95, dt: 5.6e-05, norm rho: 0.01129, rho * u: 0.0306, energy: 0.04086
4071 Time: 2, dt: 5.6e-05, norm rho: 0.01189, rho * u: 0.03145, energy: 0.04417
4072 
4073 +-------------------------------------------+------------------+------------+------------------+
4074 | Total wallclock time elapsed | 5.396e+04s 151 | 5.396e+04s | 5.396e+04s 0 |
4075 | | | |
4076 | Section | no. calls | min time rank | avg time | max time rank |
4077 +-------------------------------------------+------------------+------------+------------------+
4078 | compute errors | 41 | 2.632s 178 | 7.221s | 16.56s 0 |
4079 | compute transport speed | 7072 | 714s 193 | 1553s | 3351s 0 |
4080 | output | 41 | 8065s 176 | 8070s | 8079s 0 |
4081 | rk time stepping total | 35350 | 4.25e+04s 0 | 4.43e+04s | 4.515e+04s 193 |
4082 | rk_stage - integrals L_h | 176750 | 2.936e+04s 134 | 3.222e+04s | 3.67e+04s 99 |
4083 | rk_stage - inv mass + vec upd | 176750 | 7004s 99 | 1.207e+04s | 1.55e+04s 132 |
4084 +-------------------------------------------+------------------+------------+------------------+
4085 @endcode
4086 This simulation had nearly 2 billion unknowns -- quite a large
4087 computation indeed, and still only needed around 1.5 seconds per time
4088 step.
4089 
4090 
4091 <a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
4092 
4093 
4094 The code presented here straight-forwardly extends to adaptive meshes, given
4095 appropriate indicators for setting the refinement flags. Large-scale
4096 adaptivity of a similar solver in the context of the acoustic wave equation
4097 has been achieved by the <a href="https://github.com/kronbichler/exwave">exwave
4098 project</a>. However, in the present context, the benefits of adaptivity are often
4099 limited to early times and effects close to the origin of sound waves, as the
4100 waves eventually reflect and diffract. This leads to steep gradients all over
4101 the place, similar to turbulent flow, and a more or less globally
4102 refined mesh.
4103 
4104 Another topic that we did not discuss in the results section is a comparison
4105 of different time integration schemes. The program provides four variants of
4106 low-storage Runga--Kutta integrators that each have slightly different
4107 accuracy and stability behavior. Among the schemes implemented here, the
4108 higher-order ones provide additional accuracy but come with slightly lower
4109 efficiency in terms of step size per stage before they violate the CFL
4110 condition. An interesting extension would be to compare the low-storage
4111 variants proposed here with standard Runge--Kutta integrators or to use vector
4112 operations that are run separate from the mass matrix operation and compare
4113 performance.
4114 
4115 
4116 <a name="Moreadvancednumericalfluxfunctionsandskewsymmetricformulations"></a><h4>More advanced numerical flux functions and skew-symmetric formulations</h4>
4117 
4118 
4119 As mentioned in the introduction, the modified Lax--Friedrichs flux and the
4120 HLL flux employed in this program are only two variants of a large body of
4121 numerical fluxes available in the literature on the Euler equations. One
4122 example is the HLLC flux (Harten-Lax-van Leer-Contact) flux which adds the
4123 effect of rarefaction waves missing in the HLL flux, or the Roe flux. As
4124 mentioned in the introduction, the effect of numerical fluxes on high-order DG
4125 schemes is debatable (unlike for the case of low-order discretizations).
4126 
4127 A related improvement to increase the stability of the solver is to also
4128 consider the spatial integral terms. A shortcoming in the rather naive
4129 implementation used above is the fact that the energy conservation of the
4130 original Euler equations (in the absence of shocks) only holds up to a
4131 discretization error. If the solution is under-resolved, the discretization
4132 error can give rise to an increase in the numerical energy and eventually
4133 render the discretization unstable. This is because of the inexact numerical
4134 integration of the terms in the Euler equations, which both contain rational
4135 nonlinearities and higher-degree content from curved cells. A way out of this
4136 dilemma are so-called skew-symmetric formulations, see @cite Gassner2013 for a
4137 simple variant. Skew symmetry means that switching the role of the solution
4138 @f$\mathbf{w}@f$ and test functions @f$\mathbf{v}@f$ in the weak form produces the
4139 exact negative of the original quantity, apart from some boundary terms. In
4140 the discrete setting, the challenge is to keep this skew symmetry also when
4141 the integrals are only computed approximately (in the continuous case,
4142 skew-symmetry is a consequence of integration by parts). Skew-symmetric
4143 numerical schemes balance spatial derivatives in the conservative form
4144 @f$(\nabla \mathbf v, \mathbf{F}(\mathbf w))_{K}@f$ with contributions in the
4145 convective form @f$(\mathbf v, \tilde{\mathbf{F}}(\mathbf w)\nabla
4146 \mathbf{w})_{K}@f$ for some @f$\tilde{\mathbf{F}}@f$. The precise terms depend on
4147 the equation and the integration formula, and can in some cases by understood
4148 by special skew-symmetric finite difference schemes.
4149 
4150 
4151 <a name="Equippingthecodeforsupersoniccalculations"></a><h4>Equipping the code for supersonic calculations</h4>
4152 
4153 
4154 As mentioned in the introduction, the solution to the Euler equations develops
4155 shocks as the Mach number increases, which require additional mechanisms to
4156 stabilize the scheme, e.g. in the form of limiters. The main challenge besides
4157 actually implementing the limiter or artificial viscosity approach would be to
4158 load-balance the computations, as the additional computations involved for
4159 limiting the oscillations in troubled cells would make them more expensive than the
4160 plain DG cells without limiting. Furthermore, additional numerical fluxes that
4161 better cope with the discontinuities would also be an option.
4162 
4163 One ingredient also necessary for supersonic flows are appropriate boundary
4164 conditions. As opposed to the subsonic outflow boundaries discussed in the
4165 introduction and implemented in the program, all characteristics are outgoing
4166 for supersonic outflow boundaries, so we do not want to prescribe any external
4167 data,
4168 @f[
4169 \mathbf{w}^+ = \mathbf{w}^- = \begin{pmatrix} \rho^-\\
4170 (\rho \mathbf u)^- \\ E^-\end{pmatrix} \quad
4171  \text{(Neumann)}.
4172 @f]
4173 
4174 In the code, we would simply add the additional statement
4175 @code
4176  else if (supersonic_outflow_boundaries.find(boundary_id) !=
4177  supersonic_outflow_boundaries.end())
4178  {
4179  w_p = w_m;
4180  at_outflow = true;
4181  }
4182 @endcode
4183 in the `local_apply_boundary_face()` function.
4184 
4185 <a name="ExtensiontothelinearizedEulerequations"></a><h4>Extension to the linearized Euler equations</h4>
4186 
4187 
4188 When the interest with an Euler solution is mostly in the propagation of sound
4189 waves, it often makes sense to linearize the Euler equations around a
4190 background state, i.e., a given density, velocity and energy (or pressure)
4191 field, and only compute the change against these fields. This is the setting
4192 of the wide field of aeroacoustics. Even though the resolution requirements
4193 are sometimes considerably reduced, implementation gets somewhat more
4194 complicated as the linearization gives rise to additional terms. From a code
4195 perspective, in the operator evaluation we also need to equip the code with
4196 the state to linearize against. This information can be provided either by
4197 analytical functions (that are evaluated in terms of the position of the
4198 quadrature points) or by a vector similar to the solution. Based on that
4199 vector, we would create an additional FEEvaluation object to read from it and
4200 provide the values of the field at quadrature points. If the background
4201 velocity is zero and the density is constant, the linearized Euler equations
4202 further simplify and can equivalently be written in the form of the
4203 acoustic wave equation.
4204 
4205 A challenge in the context of sound propagation is often the definition of
4206 boundary conditions, as the computational domain needs to be of finite size,
4207 whereas the actual simulation often spans an infinite (or at least much
4208 larger) physical domain. Conventional Dirichlet or Neumann boundary conditions
4209 give rise to reflections of the sound waves that eventually propagate back to
4210 the region of interest and spoil the solution. Therefore, various variants of
4211 non-reflecting boundary conditions or sponge layers, often in the form of
4212 <a
4213 href="https://en.wikipedia.org/wiki/Perfectly_matched_layer">perfectly
4214 matched layers</a> -- where the solution is damped without reflection
4215 -- are common.
4216 
4217 
4218 <a name="ExtensiontothecompressibleNavierStokesequation"></a><h4>Extension to the compressible Navier-Stokes equation</h4>
4219 
4220 
4221 The solver presented in this tutorial program can also be extended to the
4222 compressible Navier--Stokes equations by adding viscous terms, as described in
4223 @cite FehnWallKronbichler2019. To keep as much of the performance obtained
4224 here despite the additional cost of elliptic terms, e.g. via an interior
4225 penalty method, one can switch the basis from FE_DGQ to FE_DGQHermite like in
4226 the @ref step_59 "step-59" tutorial program.
4227  *
4228  *
4229 <a name="PlainProg"></a>
4230 <h1> The plain program</h1>
4231 @include "step-67.cc"
4232 */
Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
Utilities::pow
constexpr T pow(const T base, const int iexp)
Definition: utilities.h:476
internal::MatrixFreeFunctions::affine
@ affine
Definition: mapping_info.h:62
LAPACKSupport::diagonal
@ diagonal
Matrix is diagonal.
Definition: lapack_support.h:121
internal::QGaussLobatto::gamma
long double gamma(const unsigned int n)
Definition: quadrature_lib.cc:96
DataPostprocessor
Definition: data_postprocessor.h:502
DoFTools::nonzero
@ nonzero
Definition: dof_tools.h:244
Functions
Definition: flow_function.h:28
Utilities::System::get_time
std::string get_time()
Definition: utilities.cc:1020
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
FEEvaluationAccess
Definition: fe_evaluation.h:1187
TimerOutput::Scope
Definition: timer.h:554
Utilities::int_to_string
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:474
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:133
dealii
Definition: namespace_dealii.h:25
internal::SymmetricTensorAccessors::merge
constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
Definition: symmetric_tensor.h:146
LinearAlgebra
Definition: communication_pattern_base.h:30
LinearAlgebra::distributed::Vector< Number >
DataComponentInterpretation::component_is_scalar
@ component_is_scalar
Definition: data_component_interpretation.h:55
Differentiation::SD::OptimizerType::lambda
@ lambda
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Triangulation< dim >
VectorizedArray::max
VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
Definition: vectorization.h:5466
MappingQGeneric
Definition: mapping_q_generic.h:135
DoFRenumbering::downstream
void downstream(DoFHandlerType &dof_handler, const Tensor< 1, DoFHandlerType::space_dimension > &direction, const bool dof_wise_renumbering=false)
Definition: dof_renumbering.cc:1735
VectorizedArray::min
VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &x, const ::VectorizedArray< Number, width > &y)
Definition: vectorization.h:5483
FEEvaluationBase::inverse_jacobian
Tensor< 2, dim, VectorizedArrayType > inverse_jacobian(const unsigned int q_index) const
Functions::ConstantFunction
Definition: function.h:412
VectorType
Patterns::Tools::to_string
std::string to_string(const T &t)
Definition: patterns.h:2360
LAPACKSupport::L
static const char L
Definition: lapack_support.h:171
VectorizedArray
Definition: vectorization.h:395
MatrixFree::AdditionalData::mapping_update_flags
UpdateFlags mapping_update_flags
Definition: matrix_free.h:360
VectorTools::integrate_difference
void integrate_difference(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const InVector &fe_function, const Function< spacedim, typename InVector::value_type > &exact_solution, OutVector &difference, const Quadrature< dim > &q, const NormType &norm, const Function< spacedim, double > *weight=nullptr, const double exponent=2.)
GridGenerator::channel_with_cylinder
void channel_with_cylinder(Triangulation< dim > &tria, const double shell_region_width=0.03, const unsigned int n_shells=2, const double skewness=2.0, const bool colorize=false)
GridGenerator::hyper_rectangle
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
IteratorState::valid
@ valid
Iterator points to a valid object.
Definition: tria_iterator_base.h:38
DoFHandler::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorTools::project
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >(0)), const bool project_to_boundary_first=false)
MatrixFree
Definition: matrix_free.h:117
types::boundary_id
unsigned int boundary_id
Definition: types.h:129
DoFTools::always
@ always
Definition: dof_tools.h:236
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
second
Point< 2 > second
Definition: grid_out.cc:4353
DEAL_II_ALWAYS_INLINE
#define DEAL_II_ALWAYS_INLINE
Definition: config.h:99
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1071
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:140
internal::p4est::functions
int(&) functions(const void *v1, const void *v2)
Definition: p4est_wrappers.cc:339
DoFHandler< dim >
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
deallog
LogStream deallog
Definition: logstream.cc:37
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
DoFHandler::distribute_dofs
virtual void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
Definition: dof_handler.cc:1247
FEEvaluation::n_q_points
const unsigned int n_q_points
Definition: fe_evaluation.h:2624
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
Utilities::CUDA::free
void free(T *&pointer)
Definition: cuda.h:96
MatrixFree::AdditionalData::mapping_update_flags_boundary_faces
UpdateFlags mapping_update_flags_boundary_faces
Definition: matrix_free.h:381
level
unsigned int level
Definition: grid_out.cc:4355
TimerOutput
Definition: timer.h:546
TensorAccessors::extract
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
Definition: tensor_accessors.h:226
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
DataPostprocessor::evaluate_vector_field
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double >> &computed_quantities) const
Definition: data_postprocessor.cc:37
GridGenerator::cylinder
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
LAPACKSupport::T
static const char T
Definition: lapack_support.h:163
Physics::Elasticity::Kinematics::w
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Mapping< dim >
Point::operator*
Point< dim, typename ProductType< Number, typename EnableIfScalar< OtherNumber >::type >::type > operator*(const OtherNumber) const
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition: lapack_support.h:115
MatrixFree::AdditionalData
Definition: matrix_free.h:182
Algorithms::Events::initial
const Event initial
Definition: event.cc:65
TimeStepping::invalid
@ invalid
Definition: time_stepping.h:71
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
internal::reinit
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition: matrix_block.h:621
Tensor
Definition: tensor.h:450
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
GridTools::scale
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:837
update_gradients
@ update_gradients
Shape function gradients.
Definition: fe_update_flags.h:84
MatrixFree::n_active_entries_per_cell_batch
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_number) const
SymmetricTensor::sum
SymmetricTensor< rank, dim, Number > sum(const SymmetricTensor< rank, dim, Number > &local, const MPI_Comm &mpi_communicator)
VectorizedArrayBase< VectorizedArray< Number, width >, 1 >::size
static constexpr std::size_t size()
Definition: vectorization.h:261
FEFaceEvaluation::evaluate
void evaluate(const bool evaluate_values, const bool evaluate_gradients)
MatrixFree::AdditionalData::mapping_update_flags_inner_faces
UpdateFlags mapping_update_flags_inner_faces
Definition: matrix_free.h:402
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
MatrixFreeOperators::CellwiseInverseMassMatrix
Definition: operators.h:624
parallel::distributed::Triangulation< dim >
numbers::E
static constexpr double E
Definition: numbers.h:212
std_cxx17::apply
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std_cxx14::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:40
Threads::internal::call
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
Definition: thread_management.h:607
LAPACKSupport::general
@ general
No special properties.
Definition: lapack_support.h:113
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
Utilities::MPI::sum
T sum(const T &t, const MPI_Comm &mpi_communicator)
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1579
LinearAlgebra::distributed::Vector::sadd
virtual void sadd(const Number s, const Number a, const VectorSpaceVector< Number > &V) override
Function::value
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
numbers
Definition: numbers.h:207
types
Definition: types.h:31
FiniteElement
Definition: fe.h:648
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:66
DataOutBase::VtkFlags
Definition: data_out_base.h:1095
QGauss
Definition: quadrature_lib.h:40
SIMDComparison::equal
@ equal
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
TimerOutput::print_wall_time_statistics
void print_wall_time_statistics(const MPI_Comm mpi_comm, const double print_quantile=0.) const
Definition: timer.cc:844
LocalIntegrators::L2::L2
void L2(Vector< number > &result, const FEValuesBase< dim > &fe, const std::vector< double > &input, const double factor=1.)
Definition: l2.h:171
internal::dummy
const types::global_dof_index * dummy()
Definition: dof_handler.cc:59
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandlerType &)
LogStream::depth_console
unsigned int depth_console(const unsigned int n)
Definition: logstream.cc:349
unsigned int
value
static const bool value
Definition: dof_tools_constraints.cc:433
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:174
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
AffineConstraints
Definition: affine_constraints.h:180
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
update_normal_vectors
@ update_normal_vectors
Normal vectors.
Definition: fe_update_flags.h:136
AdaptationStrategies::Refinement::split
std::vector< value_type > split(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
std::sqrt
inline ::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5412
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
DataOutInterface::write_vtu_in_parallel
void write_vtu_in_parallel(const std::string &filename, MPI_Comm comm) const
Definition: data_out_base.cc:6886
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
SymmetricTensor::determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2645
VectorizedArray::sqrt
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Definition: vectorization.h:5412
DataOutInterface::set_flags
void set_flags(const FlagType &flags)
Definition: data_out_base.cc:7837
FEFaceEvaluation::gather_evaluate
void gather_evaluate(const VectorType &input_vector, const bool evaluate_values, const bool evaluate_gradients)
Utilities::truncate_to_n_digits
Number truncate_to_n_digits(const Number number, const unsigned int n_digits)
Definition: utilities.cc:582
MatrixFree::AdditionalData::tasks_parallel_scheme
TasksParallelScheme tasks_parallel_scheme
Definition: matrix_free.h:336
std::pow
inline ::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &x, const Number p)
Definition: vectorization.h:5428
LinearAlgebra::distributed::Vector::reinit
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
ConditionalOStream
Definition: conditional_ostream.h:81
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
MatrixFreeOperators::CellwiseInverseMassMatrix::transform_from_q_points_to_basis
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
Definition: operators.h:1079
Utilities::MPI::min
T min(const T &t, const MPI_Comm &mpi_communicator)
DataOutBase::VtkFlags::write_higher_order_cells
bool write_higher_order_cells
Definition: data_out_base.h:1178
DerivativeApproximation::internal::approximate
void approximate(SynchronousIterators< std::tuple< TriaActiveIterator< ::DoFCellAccessor< DoFHandlerType< dim, spacedim >, false >>, Vector< float >::iterator >> const &cell, const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
Definition: derivative_approximation.cc:924
Point< dim >
SymmetricTensor::eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
Utilities::System::get_current_vectorization_level
const std::string get_current_vectorization_level()
Definition: utilities.cc:945
Differentiation::SD::sign
Expression sign(const Expression &x)
Definition: symengine_math.cc:280
FEEvaluationBase
Definition: fe_evaluation.h:97
DataPostprocessor::get_needed_update_flags
virtual UpdateFlags get_needed_update_flags() const =0
FEEvaluation
Definition: fe_evaluation.h:57
Function
Definition: function.h:151
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
internal::TriangulationImplementation::n_cells
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12618
DataComponentInterpretation::DataComponentInterpretation
DataComponentInterpretation
Definition: data_component_interpretation.h:49
FEFaceEvaluation
Definition: fe_evaluation.h:2681
MatrixFree::loop
void loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false, const DataAccessOnFaces dst_vector_face_access=DataAccessOnFaces::unspecified, const DataAccessOnFaces src_vector_face_access=DataAccessOnFaces::unspecified) const
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
DataPostprocessor::get_names
virtual std::vector< std::string > get_names() const =0
Quadrature
Definition: quadrature.h:85
first
Point< 2 > first
Definition: grid_out.cc:4352
DataOut< dim >
numbers::PI
static constexpr double PI
Definition: numbers.h:237
Vector
Definition: mapping_q1_eulerian.h:32
FESystem
Definition: fe.h:44
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
parallel
Definition: distributed.h:416
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:470
DataPostprocessorInputs::Vector
Definition: data_postprocessor.h:318
Utilities
Definition: cuda.h:31
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
DataPostprocessor::get_data_component_interpretation
virtual std::vector< DataComponentInterpretation::DataComponentInterpretation > get_data_component_interpretation() const
Definition: data_postprocessor.cc:48
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:67
MatrixFree::cell_loop
void cell_loop(const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &cell_operation, OutVector &dst, const InVector &src, const bool zero_dst_vector=false) const
DataComponentInterpretation::component_is_part_of_vector
@ component_is_part_of_vector
Definition: data_component_interpretation.h:61
DoFHandler::n_dofs
types::global_dof_index n_dofs() const
Utilities::MPI::MPI_InitFinalize
Definition: mpi.h:828
MatrixFree::reinit
void reinit(const Mapping< dim > &mapping, const DoFHandlerType &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
DataOut_DoFData< DoFHandler< dim >, DoFHandler< dim > ::dimension, DoFHandler< dim > ::space_dimension >::add_data_vector
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation=std::vector< DataComponentInterpretation::DataComponentInterpretation >())
Definition: data_out_dof_data.h:1090