900 *
constexpr LowStorageRungeKuttaScheme lsrk_scheme = stage_5_order_4;
904 * Eventually, we select a detail of the spatial discretization, namely the
905 * numerical flux (Riemann solver) at the faces between cells. For
this
906 * program, we have implemented a modified variant of the Lax--Friedrichs
907 * flux and the Harten--Lax--van Leer (HLL) flux.
910 *
enum EulerNumericalFlux
912 * lax_friedrichs_modified,
913 * harten_lax_vanleer,
915 *
constexpr EulerNumericalFlux numerical_flux_type = lax_friedrichs_modified;
922 * <a name=
"step_67-Equationdata"></a>
923 * <h3>Equation
data</h3>
927 * We now define a
class with the exact solution for the test case 0 and one
928 * with a background flow field
for test
case 1 of the channel. Given that
929 * the Euler equations are a problem with @f$d+2@f$ equations in @f$d@f$ dimensions,
930 * we need to tell the
Function base
class about the correct number of
935 *
class ExactSolution :
public Function<dim>
938 * ExactSolution(
const double time)
943 *
const unsigned int component = 0)
const override;
950 * As far as the actual function implemented is concerned, the analytical
951 * test
case is an isentropic vortex
case (see
e.g. the book by Hesthaven
952 * and Warburton, Example 6.1 in Section 6.6 on page 209) which fulfills the
953 * Euler equations with zero force term on the right hand side. Given that
954 * definition, we
return either the density, the momentum, or the energy
955 * depending on which component is requested. Note that the original
956 * definition of the density involves the @f$\frac{1}{\
gamma -1}@f$-th power of
957 * some expression. Since `
std::pow()` has pretty slow implementations on
958 * some systems, we replace it by logarithm followed by exponentiation (of
959 * base 2), which is mathematically equivalent but usually much better
960 * optimized. This formula might lose accuracy in the last digits
961 *
for very small
numbers compared to `
std::pow()`, but we are happy with
962 * it anyway, since small
numbers map to
data close to 1.
966 * For the channel test
case, we simply select a density of 1, a velocity of
967 * 0.4 in @f$x@f$ direction and zero in the other directions, and an energy that
968 * corresponds to a speed of sound of 1.3 measured against the background
969 * velocity field, computed from the relation @f$E = \frac{c^2}{\gamma (\gamma
970 * -1)} + \frac 12 \rho \|u\|^2@f$.
974 *
double ExactSolution<dim>::value(
const Point<dim> &x,
975 *
const unsigned int component)
const
977 *
const double t = this->
get_time();
983 *
Assert(dim == 2, ExcNotImplemented());
984 *
const double beta = 5;
988 *
const double radius_sqr =
989 * (x - x0).
norm_square() - 2. * (x[0] - x0[0]) * t + t * t;
990 *
const double factor =
992 *
const double density_log = std::log2(
993 *
std::abs(1. - (gamma - 1.) / gamma * 0.25 * factor * factor));
994 *
const double density = std::exp2(density_log * (1. / (gamma - 1.)));
995 *
const double u = 1. - factor * (x[1] - x0[1]);
996 *
const double v = factor * (x[0] - t - x0[0]);
998 *
if (component == 0)
1000 *
else if (component == 1)
1001 *
return density * u;
1002 *
else if (component == 2)
1003 *
return density * v;
1006 *
const double pressure =
1007 * std::exp2(density_log * (gamma / (gamma - 1.)));
1008 *
return pressure / (
gamma - 1.) +
1009 * 0.5 * (density * u * u + density * v * v);
1015 *
if (component == 0)
1017 *
else if (component == 1)
1019 *
else if (component == dim + 1)
1020 *
return 3.097857142857143;
1036 * <a name=
"step_67-LowstorageexplicitRungeKuttatimeintegrators"></a>
1037 * <h3>Low-storage
explicit Runge--Kutta time integrators</h3>
1041 * The next few lines implement a few low-storage variants of Runge--Kutta
1042 * methods. These methods have specific Butcher tableaux with coefficients
1043 * @f$b_i@f$ and @f$a_i@f$ as shown in the introduction. As usual in Runge--Kutta
1044 * method, we can deduce time steps, @f$c_i = \sum_{j=1}^{i-2} b_i + a_{i-1}@f$
1045 * from those coefficients. The main advantage of
this kind of scheme is the
1046 * fact that only two vectors are needed per stage, namely the accumulated
1047 * part of the solution @f$\mathbf{
w}@f$ (that will hold the solution
1048 * @f$\mathbf{
w}^{n+1}@f$ at the
new time @f$t^{n+1}@f$ after the last stage), the
1049 * update vector @f$\mathbf{r}_i@f$ that gets evaluated during the stages, plus
1050 * one vector @f$\mathbf{k}_i@f$ to hold the evaluation of the
operator. Such a
1051 * Runge--Kutta setup reduces the memory storage and memory access. As the
1052 * memory bandwidth is often the performance-limiting factor on modern
1053 * hardware when the evaluation of the differential
operator is
1054 * well-optimized, performance can be improved over standard time
1055 * integrators. This is
true also when taking into account that a
1056 * conventional Runge--Kutta scheme might allow
for slightly larger time
1057 * steps as more free parameters allow
for better stability properties.
1061 * In
this tutorial programs, we concentrate on a few variants of
1062 * low-storage schemes defined in the article by Kennedy, Carpenter, and
1063 * Lewis (2000), as well as one variant described by Tselios and Simos
1064 * (2007). There is a large series of other schemes available, which could
1065 * be addressed by additional sets of coefficients or slightly different
1070 * We define a single
class for the four integrators, distinguished by the
1071 *
enum described above. To each scheme, we then fill the vectors
for the
1072 * @f$b_i@f$ and @f$a_i@f$ to the given variables in the
class.
1075 *
class LowStorageRungeKuttaIntegrator
1078 * LowStorageRungeKuttaIntegrator(const LowStorageRungeKuttaScheme scheme)
1083 * First comes the three-stage scheme of order three by Kennedy et al.
1084 * (2000). While its stability region is significantly smaller than
for
1085 * the other schemes, it only involves three stages, so it is very
1086 * competitive in terms of the work per stage.
1091 *
case stage_3_order_3:
1099 * The next scheme is a five-stage scheme of order four, again
1100 * defined in the paper by Kennedy et al. (2000).
1103 *
case stage_5_order_4:
1111 * The following scheme of seven stages and order four has been
1112 * explicitly derived
for acoustics problems. It is a balance of
1113 * accuracy
for imaginary
eigenvalues among fourth order schemes,
1114 * combined with a large stability region. Since DG schemes are
1115 * dissipative among the highest frequencies,
this does not
1116 * necessarily translate to the highest possible time step per
1117 * stage. In the context of the present tutorial program, the
1118 * numerical flux plays a crucial role in the dissipation and thus
1119 * also the maximal stable time step
size. For the modified
1120 * Lax--Friedrichs flux,
this scheme is similar to the
1121 * `stage_5_order_4` scheme in terms of step
size per stage
if only
1122 * stability is considered, but somewhat less efficient
for the HLL
1126 *
case stage_7_order_4:
1134 * The last scheme included here is the nine-stage scheme of order
1135 * five from Kennedy et al. (2000). It is the most accurate among
1136 * the schemes used here, but the higher order of accuracy
1137 * sacrifices some stability, so the step length normalized per
1138 * stage is less than
for the fourth order schemes.
1141 *
case stage_9_order_5:
1152 * rk_integrator(lsrk);
1153 * rk_integrator.get_coefficients(ai, bi, ci);
1156 *
unsigned int n_stages() const
1163 * The main function of the time integrator is to go through the stages,
1164 * evaluate the
operator, prepare the @f$\mathbf{r}_i@f$ vector
for the next
1165 * evaluation, and update the solution vector @f$\mathbf{
w}@f$. We hand off
1166 * the work to the `pde_operator` involved in order to be able to
merge
1167 * the vector operations of the Runge--Kutta setup with the evaluation of
1168 * the differential
operator for better performance, so all we
do here is
1169 * to delegate the vectors and coefficients.
1173 * We separately call the
operator for the
first stage because we need
1174 * slightly modified arguments there: We evaluate the solution from
1175 * the old solution @f$\mathbf{
w}^n@f$ rather than a @f$\mathbf r_i@f$ vector, so
1176 * the
first argument is `solution`. We here let the stage vector
1177 * @f$\mathbf{r}_i@f$ also hold the temporary result of the evaluation, as it
1178 * is not used otherwise. For all subsequent stages, we use the vector
1179 * `vec_ki` as the
second vector argument to store the result of the
1180 *
operator evaluation. Finally, when we are at the last stage, we must
1181 * skip the computation of the vector @f$\mathbf{r}_{s+1}@f$ as there is no
1182 * coefficient @f$a_s@f$ available (nor will it be used).
1185 *
template <
typename VectorType,
typename Operator>
1186 *
void perform_time_step(
const Operator &pde_operator,
1187 *
const double current_time,
1188 *
const double time_step,
1189 * VectorType &solution,
1190 * VectorType &vec_ri,
1191 * VectorType &vec_ki)
const
1195 * pde_operator.perform_stage(current_time,
1196 * bi[0] * time_step,
1197 * ai[0] * time_step,
1203 *
for (
unsigned int stage = 1; stage < bi.size(); ++stage)
1205 *
const double c_i = ci[stage];
1206 * pde_operator.perform_stage(current_time + c_i * time_step,
1207 * bi[stage] * time_step,
1208 * (stage == bi.size() - 1 ?
1210 * ai[stage] * time_step),
1219 * std::vector<double> bi;
1220 * std::vector<double> ai;
1221 * std::vector<double> ci;
1229 * <a name=
"step_67-ImplementationofpointwiseoperationsoftheEulerequations"></a>
1230 * <h3>Implementation of
point-wise operations of the Euler equations</h3>
1234 * In the following
functions, we implement the various problem-specific
1235 * operators pertaining to the Euler equations. Each function acts on the
1236 * vector of conserved variables @f$[\rho, \rho\mathbf{u},
E]@f$ that we hold in
1237 * the solution vectors, and computes various derived quantities.
1241 * First out is the computation of the velocity, that we derive from the
1242 * momentum variable @f$\rho \mathbf{u}@f$ by division by @f$\rho@f$. One thing to
1243 * note here is that we decorate all those
functions with the keyword
1245 * compiler-specific keyword that tells the compiler to never create a
1246 * function call
for any of those
functions, and instead move the
1248 * href=
"https://en.wikipedia.org/wiki/Inline_function">
inline</a> to where
1249 * they are called. This is critical
for performance because we call into some
1250 * of those
functions millions or billions of times: For example, we both use
1251 * the velocity
for the computation of the flux further down, but also
for the
1252 * computation of the pressure, and both of these places are evaluated at
1253 * every quadrature
point of every cell. Making sure these
functions are
1254 * inlined ensures not only that the processor does not have to execute a jump
1255 * instruction into the function (and the corresponding
return jump), but also
1256 * that the compiler can re-use intermediate information from one function
's
1257 * context in code that comes after the place where the function was called.
1258 * (We note that compilers are generally quite good at figuring out which
1259 * functions to inline by themselves. Here is a place where compilers may or
1260 * may not have figured it out by themselves but where we know for sure that
1261 * inlining is a win.)
1265 * Another trick we apply is a separate variable for the inverse density
1266 * @f$\frac{1}{\rho}@f$. This enables the compiler to only perform a single
1267 * division for the flux, despite the division being used at several
1268 * places. As divisions are around ten to twenty times as expensive as
1269 * multiplications or additions, avoiding redundant divisions is crucial for
1270 * performance. We note that taking the inverse first and later multiplying
1271 * with it is not equivalent to a division in floating point arithmetic due
1272 * to roundoff effects, so the compiler is not allowed to exchange one way by
1273 * the other with standard optimization flags. However, it is also not
1274 * particularly difficult to write the code in the right way.
1278 * To summarize, the chosen strategy of always inlining and careful
1279 * definition of expensive arithmetic operations allows us to write compact
1280 * code without passing all intermediate results around, despite making sure
1281 * that the code maps to excellent machine code.
1284 * template <int dim, typename Number>
1285 * inline DEAL_II_ALWAYS_INLINE
1286 * Tensor<1, dim, Number>
1287 * euler_velocity(const Tensor<1, dim + 2, Number> &conserved_variables)
1289 * const Number inverse_density = Number(1.) / conserved_variables[0];
1291 * Tensor<1, dim, Number> velocity;
1292 * for (unsigned int d = 0; d < dim; ++d)
1293 * velocity[d] = conserved_variables[1 + d] * inverse_density;
1300 * The next function computes the pressure from the vector of conserved
1301 * variables, using the formula @f$p = (\gamma - 1) \left(E - \frac 12 \rho
1302 * \mathbf{u}\cdot \mathbf{u}\right)@f$. As explained above, we use the
1303 * velocity from the `euler_velocity()` function. Note that we need to
1304 * specify the first template argument `dim` here because the compiler is
1305 * not able to deduce it from the arguments of the tensor, whereas the
1306 * second argument (number type) can be automatically deduced.
1309 * template <int dim, typename Number>
1310 * inline DEAL_II_ALWAYS_INLINE
1312 * euler_pressure(const Tensor<1, dim + 2, Number> &conserved_variables)
1314 * const Tensor<1, dim, Number> velocity =
1315 * euler_velocity<dim>(conserved_variables);
1317 * Number rho_u_dot_u = conserved_variables[1] * velocity[0];
1318 * for (unsigned int d = 1; d < dim; ++d)
1319 * rho_u_dot_u += conserved_variables[1 + d] * velocity[d];
1321 * return (gamma - 1.) * (conserved_variables[dim + 1] - 0.5 * rho_u_dot_u);
1326 * Here is the definition of the Euler flux function, i.e., the definition
1327 * of the actual equation. Given the velocity and pressure (that the
1328 * compiler optimization will make sure are done only once), this is
1329 * straight-forward given the equation stated in the introduction.
1332 * template <int dim, typename Number>
1333 * inline DEAL_II_ALWAYS_INLINE
1334 * Tensor<1, dim + 2, Tensor<1, dim, Number>>
1335 * euler_flux(const Tensor<1, dim + 2, Number> &conserved_variables)
1337 * const Tensor<1, dim, Number> velocity =
1338 * euler_velocity<dim>(conserved_variables);
1339 * const Number pressure = euler_pressure<dim>(conserved_variables);
1341 * Tensor<1, dim + 2, Tensor<1, dim, Number>> flux;
1342 * for (unsigned int d = 0; d < dim; ++d)
1344 * flux[0][d] = conserved_variables[1 + d];
1345 * for (unsigned int e = 0; e < dim; ++e)
1346 * flux[e + 1][d] = conserved_variables[e + 1] * velocity[d];
1347 * flux[d + 1][d] += pressure;
1348 * flux[dim + 1][d] =
1349 * velocity[d] * (conserved_variables[dim + 1] + pressure);
1357 * This next function is a helper to simplify the implementation of the
1358 * numerical flux, implementing the action of a tensor of tensors (with
1359 * non-standard outer dimension of size `dim + 2`, so the standard overloads
1360 * provided by deal.II's tensor classes
do not
apply here) with another
1361 * tensor of the same inner dimension, i.e., a matrix-vector product.
1364 * template <int n_components, int dim, typename Number>
1371 *
for (
unsigned int d = 0;
d < n_components; ++
d)
1372 * result[d] = matrix[d] * vector;
1378 * This function implements the numerical flux (Riemann solver). It gets the
1379 * state from the two sides of an
interface and the normal vector, oriented
1380 * from the side of the solution @f$\mathbf{
w}^-@f$ towards the solution
1381 * @f$\mathbf{
w}^+@f$. In finite
volume methods which rely on piece-wise
1382 *
constant data, the numerical flux is the central ingredient as it is the
1383 * only place where the physical information is entered. In DG methods, the
1384 * numerical flux is less central due to the polynomials within the elements
1385 * and the physical flux used there. As a result of higher-degree
1386 * interpolation with consistent
values from both sides in the limit of a
1387 * continuous solution, the numerical flux can be seen as a control of the
1388 * jump of the solution from both sides to weakly impose continuity. It is
1389 * important to realize that a numerical flux alone cannot stabilize a
1390 * high-order DG method in the presence of shocks, and thus any DG method
1391 * must be combined with further shock-capturing techniques to handle those
1392 * cases. In
this tutorial, we focus on wave-like solutions of the Euler
1393 * equations in the subsonic regime without strong discontinuities where our
1394 * basic scheme is sufficient.
1398 * Nonetheless, the numerical flux is decisive in terms of the numerical
1399 * dissipation of the overall scheme and influences the admissible time step
1400 *
size with
explicit Runge--Kutta methods. We consider two choices, a
1401 * modified Lax--Friedrichs scheme and the widely used Harten--Lax--van Leer
1402 * (HLL) flux. For both variants, we
first need to get the velocities and
1403 * pressures from both sides of the interface and evaluate the physical
1408 * For the local Lax--Friedrichs flux, the definition is @f$\hat{\mathbf{F}}
1409 * =\frac{\mathbf{
F}(\mathbf{
w}^-)+\mathbf{
F}(\mathbf{
w}^+)}{2} +
1410 * \frac{\
lambda}{2}\left[\mathbf{
w}^--\mathbf{
w}^+\right]\otimes
1411 * \mathbf{n^-}@f$, where the factor @f$\lambda =
1412 * \max\left(\|\mathbf{u}^-\|+c^-, \|\mathbf{u}^+\|+c^+\right)@f$ gives the
1413 * maximal wave speed and @f$c = \sqrt{\
gamma p / \rho}@f$ is the speed of
1414 * sound. Here, we choose two modifications of that expression
for reasons
1415 * of computational efficiency, given the small impact of the flux on the
1416 * solution. For the above definition of the factor @f$\lambda@f$, we would need
1417 * to take four square roots, two
for the two velocity norms and two
for the
1418 * speed of sound on either side. The
first modification is hence to rather
1419 * use @f$\sqrt{\|\mathbf{u}\|^2+c^2}@f$ as an estimate of the maximal speed
1420 * (which is at most a factor of 2 away from the actual maximum, as shown in
1421 * the introduction). This allows us to pull the square root out of the
1422 * maximum and get away with a single square root computation. The
second
1423 * modification is to further relax on the parameter @f$\lambda@f$---the smaller
1424 * it is, the smaller the dissipation factor (which is multiplied by the
1425 * jump in @f$\mathbf{
w}@f$, which might result in a smaller or bigger
1426 * dissipation in the
end). This allows us to fit the spectrum into the
1427 * stability region of the
explicit Runge--Kutta integrator with bigger time
1428 * steps. However, we cannot make dissipation too small because otherwise
1429 * imaginary
eigenvalues grow larger. Finally, the current conservative
1430 * formulation is not energy-stable in the limit of @f$\lambda\to 0@f$ as it is
1431 * not skew-symmetric, and would need additional measures such as split-form
1432 * DG schemes in that
case.
1436 * For the HLL flux, we follow the formula from literature, introducing an
1437 * additional weighting of the two states from Lax--Friedrichs by a
1438 * parameter @f$s@f$. It is derived from the physical transport directions of
1439 * the Euler equations in terms of the current direction of velocity and
1440 * sound speed. For the velocity, we here choose a simple arithmetic average
1441 * which is sufficient
for DG scenarios and moderate jumps in material
1446 * Since the numerical flux is multiplied by the normal vector in the weak
1447 * form, we multiply by the result by the normal vector
for all terms in the
1448 * equation. In these multiplications, the `
operator*` defined above enables
1449 * a compact notation similar to the mathematical definition.
1453 * In
this and the following functions, we use variable suffixes `_m` and
1454 * `_p` to indicate quantities derived from @f$\mathbf{w}^-@f$ and @f$\mathbf{
w}^+@f$,
1455 * i.e.,
values "here" and
"there" relative to the current cell when looking
1456 * at a neighbor cell.
1459 *
template <
int dim,
typename Number>
1466 *
const auto velocity_m = euler_velocity<dim>(u_m);
1467 *
const auto velocity_p = euler_velocity<dim>(u_p);
1469 *
const auto pressure_m = euler_pressure<dim>(u_m);
1470 *
const auto pressure_p = euler_pressure<dim>(u_p);
1472 *
const auto flux_m = euler_flux<dim>(u_m);
1473 *
const auto flux_p = euler_flux<dim>(u_p);
1475 *
switch (numerical_flux_type)
1477 *
case lax_friedrichs_modified:
1481 * gamma * pressure_p * (1. / u_p[0]),
1482 * velocity_m.norm_square() +
1483 * gamma * pressure_m * (1. / u_m[0])));
1485 *
return 0.5 * (flux_m * normal + flux_p * normal) +
1486 * 0.5 * lambda * (u_m - u_p);
1489 *
case harten_lax_vanleer:
1491 *
const auto avg_velocity_normal =
1492 * 0.5 * ((velocity_m + velocity_p) * normal);
1495 * (pressure_p * (1. / u_p[0]) + pressure_m * (1. / u_m[0]))));
1496 *
const Number s_pos =
1497 *
std::max(Number(), avg_velocity_normal + avg_c);
1498 *
const Number s_neg =
1499 *
std::min(Number(), avg_velocity_normal - avg_c);
1500 *
const Number inverse_s = Number(1.) / (s_pos - s_neg);
1502 *
return inverse_s *
1503 * ((s_pos * (flux_m * normal) - s_neg * (flux_p * normal)) -
1504 * s_pos * s_neg * (u_m - u_p));
1519 * This and the next function are helper
functions to provide compact
1520 * evaluation calls as multiple points get batched together via a
1521 *
VectorizedArray argument (see the @ref step_37
"step-37" tutorial
for details). This
1522 * function is used
for the subsonic outflow boundary conditions where we
1523 * need to set the energy component to a prescribed
value. The next one
1524 * requests the solution on all components and is used
for inflow boundaries
1525 * where all components of the solution are set.
1528 *
template <
int dim,
typename Number>
1532 *
const unsigned int component)
1535 *
for (
unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v)
1538 *
for (
unsigned int d = 0;
d < dim; ++
d)
1539 * p[d] = p_vectorized[d][v];
1540 * result[v] = function.value(p, component);
1546 *
template <
int dim,
typename Number,
int n_components = dim + 2>
1553 *
for (
unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v)
1556 *
for (
unsigned int d = 0;
d < dim; ++
d)
1557 * p[d] = p_vectorized[d][v];
1558 *
for (
unsigned int d = 0;
d < n_components; ++
d)
1559 * result[d][v] = function.value(p, d);
1569 * <a name=
"step_67-TheEulerOperationclass"></a>
1570 * <h3>The EulerOperation
class</h3>
1574 * This
class implements the evaluators for the Euler problem, in analogy to
1575 * the `LaplaceOperator`
class of @ref step_37
"step-37" or @ref step_59
"step-59". Since the present
1576 *
operator is non-linear and does not require a
matrix interface (to be
1577 * handed over to preconditioners), we skip the various `vmult`
functions
1578 * otherwise present in
matrix-free operators and only implement an `
apply`
1579 * function as well as the combination of `
apply` with the required vector
1580 * updates
for the low-storage Runge--Kutta time integrator mentioned above
1581 * (called `perform_stage`). Furthermore, we have added three additional
1582 * functions involving matrix-free routines, namely one to compute an
1583 * estimate of the time step scaling (that is combined with the Courant
1584 * number
for the actual time step
size) based on the velocity and speed of
1585 * sound in the elements, one
for the projection of solutions (specializing
1587 * against a possible analytical solution or norms against some background
1592 * The rest of the
class is similar to other
matrix-free tutorials. As
1593 * discussed in the introduction, we provide a few
functions to allow a user
1594 * to pass in various forms of boundary conditions on different parts of the
1596 * possible body forces.
1599 *
template <
int dim,
int degree,
int n_po
ints_1d>
1600 *
class EulerOperator
1603 *
static constexpr unsigned int n_quadrature_points_1d = n_points_1d;
1613 *
void set_subsonic_outflow_boundary(
1619 *
void set_body_force(std::unique_ptr<
Function<dim>> body_force);
1621 *
void apply(
const double current_time,
1626 * perform_stage(
const Number cur_time,
1627 *
const Number factor_solution,
1628 *
const Number factor_ai,
1637 * std::array<double, 3> compute_errors(
1641 *
double compute_cell_transport_speed(
1652 * std::map<types::boundary_id, std::unique_ptr<Function<dim>>>
1653 * inflow_boundaries;
1654 * std::map<types::boundary_id, std::unique_ptr<Function<dim>>>
1655 * subsonic_outflow_boundaries;
1656 * std::set<types::boundary_id> wall_boundaries;
1657 * std::unique_ptr<Function<dim>> body_force;
1659 *
void local_apply_inverse_mass_matrix(
1663 *
const std::pair<unsigned int, unsigned int> &cell_range)
const;
1665 *
void local_apply_cell(
1669 *
const std::pair<unsigned int, unsigned int> &cell_range)
const;
1671 *
void local_apply_face(
1675 *
const std::pair<unsigned int, unsigned int> &face_range)
const;
1677 *
void local_apply_boundary_face(
1681 *
const std::pair<unsigned int, unsigned int> &face_range)
const;
1686 *
template <
int dim,
int degree,
int n_po
ints_1d>
1687 * EulerOperator<dim, degree, n_points_1d>::EulerOperator(
TimerOutput &timer)
1695 * For the initialization of the Euler
operator, we set up the
MatrixFree
1696 * variable contained in the
class. This can be done given a mapping to
1697 * describe possible curved boundaries as well as a
DoFHandler object
1698 * describing the degrees of freedom. Since we use a discontinuous Galerkin
1699 * discretization in
this tutorial program where no constraints are imposed
1700 * strongly on the solution field, we
do not need to pass in an
1702 * construction. With respect to quadrature, we want to select two different
1703 * ways of computing the underlying integrals: The
first is a flexible one,
1704 * based on a
template parameter `n_points_1d` (that will be assigned the
1705 * `n_q_points_1d`
value specified at the top of
this file). More accurate
1706 * integration is necessary to avoid the aliasing problem due to the
1707 * variable coefficients in the Euler
operator. The
second less accurate
1708 * quadrature formula is a tight one based on `fe_degree+1` and needed
for
1709 * the inverse mass
matrix. While that formula provides an exact inverse
1710 * only on affine element shapes and not on deformed elements, it enables
1711 * the fast inversion of the mass matrix by tensor product techniques,
1712 * necessary to ensure optimal computational efficiency overall.
1715 * template <int dim, int degree, int n_points_1d>
1716 *
void EulerOperator<dim, degree, n_points_1d>::reinit(
1720 *
const std::vector<const DoFHandler<dim> *> dof_handlers = {&dof_handler};
1722 *
const std::vector<const AffineConstraints<double> *> constraints = {&dummy};
1723 *
const std::vector<Quadrature<1>> quadratures = {
QGauss<1>(n_q_points_1d),
1730 * additional_data.mapping_update_flags_inner_faces =
1733 * additional_data.mapping_update_flags_boundary_faces =
1736 * additional_data.tasks_parallel_scheme =
1740 * mapping, dof_handlers, constraints, quadratures, additional_data);
1745 *
template <
int dim,
int degree,
int n_po
ints_1d>
1746 *
void EulerOperator<dim, degree, n_points_1d>::initialize_vector(
1749 *
data.initialize_dof_vector(vector);
1756 * The subsequent four member
functions are the ones that must be called from
1757 *
outside to specify the various
types of boundaries. For an inflow boundary,
1758 * we must specify all components in terms of density @f$\rho@f$, momentum @f$\rho
1759 * \mathbf{u}@f$ and energy @f$E@f$. Given
this information, we then store the
1760 * function alongside the respective boundary
id in a map member variable of
1761 *
this class. Likewise, we proceed
for the subsonic outflow boundaries (where
1762 * we request a function as well, which we use to retrieve the energy) and
for
1763 * wall (no-penetration) boundaries where we impose zero normal velocity (no
1764 * function necessary, so we only request the boundary
id). For the present
1765 * DG code where boundary conditions are solely applied as part of the weak
1766 * form (during time integration), the call to set the boundary conditions
1767 * can appear both before or after the `
reinit()` call to
this class. This
1768 * is different from continuous finite element codes where the boundary
1770 * sent into
MatrixFree for initialization, thus requiring to be set before
1771 * the initialization of the
matrix-free
data structures.
1775 * The checks added in each of the four function are used to
1776 * ensure that boundary conditions are mutually exclusive on the various
1777 * parts of the boundary, i.e., that a user does not accidentally designate a
1778 * boundary as both an inflow and say a subsonic outflow boundary.
1781 *
template <
int dim,
int degree,
int n_po
ints_1d>
1782 *
void EulerOperator<dim, degree, n_points_1d>::set_inflow_boundary(
1786 *
AssertThrow(subsonic_outflow_boundaries.find(boundary_id) ==
1787 * subsonic_outflow_boundaries.end() &&
1788 * wall_boundaries.find(boundary_id) == wall_boundaries.end(),
1789 * ExcMessage(
"You already set the boundary with id " +
1790 * std::to_string(
static_cast<int>(boundary_id)) +
1791 *
" to another type of boundary before now setting " +
1793 *
AssertThrow(inflow_function->n_components == dim + 2,
1794 * ExcMessage(
"Expected function with dim+2 components"));
1796 * inflow_boundaries[
boundary_id] = std::move(inflow_function);
1800 *
template <
int dim,
int degree,
int n_po
ints_1d>
1801 *
void EulerOperator<dim, degree, n_points_1d>::set_subsonic_outflow_boundary(
1805 *
AssertThrow(inflow_boundaries.find(boundary_id) ==
1806 * inflow_boundaries.end() &&
1807 * wall_boundaries.find(boundary_id) == wall_boundaries.end(),
1808 * ExcMessage(
"You already set the boundary with id " +
1809 * std::to_string(
static_cast<int>(boundary_id)) +
1810 *
" to another type of boundary before now setting " +
1811 *
"it as subsonic outflow"));
1812 *
AssertThrow(outflow_function->n_components == dim + 2,
1813 * ExcMessage(
"Expected function with dim+2 components"));
1815 * subsonic_outflow_boundaries[
boundary_id] = std::move(outflow_function);
1819 *
template <
int dim,
int degree,
int n_po
ints_1d>
1820 *
void EulerOperator<dim, degree, n_points_1d>::set_wall_boundary(
1823 *
AssertThrow(inflow_boundaries.find(boundary_id) ==
1824 * inflow_boundaries.end() &&
1825 * subsonic_outflow_boundaries.find(boundary_id) ==
1826 * subsonic_outflow_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 wall boundary"));
1832 * wall_boundaries.insert(boundary_id);
1836 *
template <
int dim,
int degree,
int n_po
ints_1d>
1837 *
void EulerOperator<dim, degree, n_points_1d>::set_body_force(
1842 * this->body_force = std::move(body_force);
1850 * <a name=
"step_67-Localevaluators"></a>
1851 * <h4>Local evaluators</h4>
1855 * Now we proceed to the local evaluators
for the Euler problem. The
1856 * evaluators are relatively simple and follow what has been presented in
1857 * @ref step_37
"step-37", @ref step_48
"step-48", or @ref step_59
"step-59". The
first notable difference is the fact
1858 * that we use an
FEEvaluation with a non-standard number of quadrature
1859 * points. Whereas we previously
always set the number of quadrature points
1860 * to
equal the polynomial degree plus one (ensuring exact integration on
1861 * affine element shapes), we now set the number quadrature points as a
1862 * separate variable (
e.g. the polynomial degree plus two or three halves of
1863 * the polynomial degree) to more accurately handle nonlinear terms. Since
1864 * the evaluator is fed with the appropriate
loop lengths via the
template
1865 * argument and keeps the number of quadrature points in the whole cell in
1867 * the more accurate formula without further changes.
1871 * The
second difference is due to the fact that we are now evaluating a
1872 * multi-component system, as opposed to the
scalar systems considered
1873 * previously. The
matrix-free framework provides several ways to handle the
1874 * multi-component
case. The variant shown here utilizes an
FEEvaluation
1875 *
object with multiple components embedded into it, specified by the fourth
1876 *
template argument `dim + 2`
for the components in the Euler system. As a
1879 * several elements), but a
Tensor of `dim+2` components. The functionality
1880 * is otherwise similar to the scalar case; it is handled by a template
1882 * variant would have been to use several
FEEvaluation objects, a scalar one
1883 * for the density, a vector-valued one with `dim` components for the
1884 * momentum, and another scalar evaluator for the energy. To ensure that
1885 * those components point to the correct part of the solution, the
1886 * constructor of
FEEvaluation takes three optional integer arguments after
1888 * multi-
DoFHandler systems (taking the
first by default), the number of the
1889 * quadrature point in case there are multiple
Quadrature objects (see more
1890 * below), and as a third argument the component within a vector system. As
1891 * we have a single vector for all components, we would go with the third
1892 * argument, and set it to `0` for the density, `1` for the vector-valued
1893 * momentum, and `dim+1` for the energy slot.
FEEvaluation then picks the
1894 * appropriate subrange of the solution vector during
1896 *
FEEvaluation::distributed_local_to_global() or the more compact
1902 * When it comes to the evaluation of the body force vector, we distinguish
1903 * between two cases for efficiency reasons: In case we have a constant
1904 * function (derived from
Functions::ConstantFunction), we can precompute
1905 * the value outside the loop over quadrature points and simply use the
1906 * value everywhere. For a more general function, we instead need to call
1907 * the `evaluate_function()` method we provided above; this path is more
1908 * expensive because we need to access the memory associated with the
1909 * quadrature point
data.
1913 * The rest follows the other tutorial programs. Since we have implemented
1914 * all physics for the Euler equations in the separate `euler_flux()`
1915 * function, all we have to do here is to call this function
1916 * given the current solution evaluated at quadrature points, returned by
1917 * `phi.get_value(q)`, and tell the
FEEvaluation object to queue the flux
1918 * for testing it by the gradients of the shape functions (which is a
Tensor
1919 * of outer `dim+2` components, each holding a tensor of `dim` components
1920 * for the @f$x,y,z@f$ component of the Euler flux). One final thing worth
1921 * mentioning is the order in which we queue the
data for testing by the
1922 * value of the test function, `phi.submit_value()`, in case we are given an
1923 * external function: We must do this after calling `phi.get_value(q)`,
1924 * because `get_value()` (reading the solution) and `submit_value()`
1925 * (queuing the value for multiplication by the test function and summation
1926 * over quadrature points) access the same underlying
data field. Here it
1927 * would be easy to achieve also without temporary variable `w_q` since
1928 * there is no mixing between values and gradients. For more complicated
1929 * setups, one has to
first copy out e.g. both the value and gradient at a
1930 * quadrature point and then queue results again by
1936 * argument of this function, which is a call-back from
MatrixFree::loop().
1937 * The interfaces imposes the present list of arguments, but since we are in
1938 * a member function where the
MatrixFree object is already available as the
1939 * `
data` variable, we stick with that to avoid confusion.
1942 * template <
int dim,
int degree,
int n_points_1d>
1943 *
void EulerOperator<dim, degree, n_points_1d>::local_apply_cell(
1947 * const
std::pair<
unsigned int,
unsigned int> &cell_range) const
1955 *
if (constant_function)
1956 * constant_body_force = evaluate_function<dim, Number, dim>(
1959 *
for (
unsigned int cell = cell_range.first;
cell < cell_range.second; ++
cell)
1967 * phi.submit_gradient(euler_flux<dim>(w_q), q);
1968 *
if (body_force.get() !=
nullptr)
1971 * constant_function ? constant_body_force :
1972 * evaluate_function<dim, Number, dim>(
1973 * *body_force, phi.quadrature_point(q));
1976 *
for (
unsigned int d = 0;
d < dim; ++
d)
1977 * forcing[d + 1] = w_q[0] * force[d];
1978 *
for (
unsigned int d = 0;
d < dim; ++
d)
1979 * forcing[dim + 1] += force[d] * w_q[d + 1];
1981 * phi.submit_value(forcing, q);
1985 * phi.integrate_scatter(((body_force.get() !=
nullptr) ?
1997 * The next function concerns the computation of integrals on interior
1998 * faces, where we need evaluators from both cells adjacent to the face. We
1999 * associate the variable `phi_m` with the solution component @f$\mathbf{
w}^-@f$
2000 * and the variable `phi_p` with the solution component @f$\mathbf{
w}^+@f$. We
2002 *
second argument, with `
true`
for the interior side and `
false`
for the
2003 * exterior side, with interior and exterior denoting the orientation with
2004 * respect to the normal vector.
2010 * and the sum factorization parts. This combined operation not only saves a
2011 * line of code, but also contains an important optimization: Given that we
2012 * use a nodal basis in terms of the Lagrange polynomials in the points of
2013 * the Gauss-Lobatto quadrature formula, only @f$(p+1)^{
d-1}@f$ out of the
2014 * @f$(p+1)^
d@f$ basis
functions evaluate to non-zero on each face. Thus, the
2015 * evaluator only accesses the necessary
data in the vector and skips the
2016 * parts which are multiplied by zero. If we had
first read the vector, we
2017 * would have needed to load all
data from the vector, as the call in
2018 * isolation would not know what
data is required in subsequent
2020 * values and derivatives, indeed all @f$(p+1)^d@f$ vector entries for each
2021 * component are needed, as the normal derivative is nonzero for all basis
2026 * The arguments to the evaluators as well as the procedure is similar to
2027 * the cell evaluation. We again use the more accurate (over-)integration
2028 * scheme due to the nonlinear terms, specified as the third template
2029 * argument in the list. At the quadrature points, we then go to our
2030 * free-standing function for the numerical flux. It receives the solution
2031 * evaluated at quadrature points from both sides (i.e., @f$\mathbf{
w}^-@f$ and
2032 * @f$\mathbf{
w}^+@f$), as well as the normal vector onto the minus side. As
2033 * explained above, the numerical flux is already multiplied by the normal
2034 * vector from the minus side. We need to
switch the sign because the
2035 * boundary term comes with a minus sign in the weak form derived in the
2036 * introduction. The flux is then queued
for testing both on the minus sign
2037 * and on the plus sign, with switched sign as the normal vector from the
2038 * plus side is exactly opposed to the one from the minus side.
2041 * template <int dim, int degree, int n_points_1d>
2042 *
void EulerOperator<dim, degree, n_points_1d>::local_apply_face(
2046 *
const std::pair<unsigned int, unsigned int> &face_range)
const
2053 *
for (
unsigned int face = face_range.first; face < face_range.second; ++face)
2055 * phi_p.reinit(face);
2058 * phi_m.reinit(face);
2061 *
for (
const unsigned int q : phi_m.quadrature_point_indices())
2063 * const auto numerical_flux =
2064 * euler_numerical_flux<dim>(phi_m.get_value(q),
2065 * phi_p.get_value(q),
2066 * phi_m.normal_vector(q));
2067 * phi_m.submit_value(-numerical_flux, q);
2068 * phi_p.submit_value(numerical_flux, q);
2080 * For faces located at the boundary, we need to impose the appropriate
2081 * boundary conditions. In
this tutorial program, we implement four cases as
2082 * mentioned above. (A fifth
case,
for supersonic outflow conditions is
2083 * discussed in the
"Results" section below.) The discontinuous Galerkin
2084 * method imposes boundary conditions not as constraints, but only
2085 * weakly. Thus, the various conditions are imposed by finding an appropriate
2086 * <i>exterior</i> quantity @f$\mathbf{
w}^+@f$ that is then handed to the
2087 * numerical flux function also used
for the interior faces. In essence,
2088 * we
"pretend" a state on the
outside of the domain in such a way that
2089 *
if that were reality, the solution of the PDE would satisfy the boundary
2090 * conditions we want.
2094 * For wall boundaries, we need to impose a no-normal-flux condition on the
2095 * momentum variable, whereas we use a Neumann condition
for the density and
2096 * energy with @f$\rho^+ = \rho^-@f$ and @f$E^+ =
E^-@f$. To achieve the no-normal
2097 * flux condition, we set the exterior
values to the interior
values and
2098 * subtract two times the velocity in wall-normal direction, i.e., in the
2099 * direction of the normal vector.
2103 * For inflow boundaries, we simply set the given Dirichlet
data
2104 * @f$\mathbf{
w}_\mathrm{D}@f$ as a boundary
value. An alternative would have been
2105 * to use @f$\mathbf{
w}^+ = -\mathbf{
w}^- + 2 \mathbf{
w}_\mathrm{D}@f$, the
2106 * so-called mirror principle.
2110 * The imposition of outflow is essentially a Neumann condition, i.e.,
2111 * setting @f$\mathbf{
w}^+ = \mathbf{
w}^-@f$. For the
case of subsonic outflow,
2112 * we still need to impose a
value for the energy, which we derive from the
2113 * respective function. A special step is needed
for the
case of
2114 * <i>backflow</i>, i.e., the
case where there is a momentum flux into the
2115 * domain on the Neumann portion. According to the literature (a fact that can
2116 * be derived by appropriate energy arguments), we must
switch to another
2117 * variant of the flux on inflow parts, see Gravemeier, Comerford,
2118 * Yoshihara, Ismail, Wall,
"A novel formulation for Neumann inflow
2119 * conditions in biomechanics", Int. J. Numer. Meth. Biomed. Eng., vol. 28
2120 * (2012). Here, the momentum term needs to be added once again, which
2121 * corresponds to removing the flux contribution on the momentum
2122 * variables. We
do this in a post-processing step, and only
for the
case
2123 * when we both are at an outflow boundary and the dot product between the
2124 * normal vector and the momentum (or, equivalently, velocity) is
2125 *
negative. As we work on
data of several quadrature points at once
for
2126 * SIMD vectorizations, we here need to explicitly
loop over the array
2127 * entries of the SIMD array.
2131 * In the implementation below, we
check for the various
types
2132 * of boundaries at the
level of quadrature points. Of course, we could also
2133 * have moved the decision out of the quadrature
point loop and treat entire
2134 * faces as of the same kind, which avoids some map/set lookups in the inner
2135 *
loop over quadrature points. However, the loss of efficiency is hardly
2136 * noticeable, so we opt
for the simpler code here. Also note that the
final
2137 * `
else` clause will
catch the
case when some part of the boundary was not
2138 * assigned any boundary condition via `EulerOperator::set_..._boundary(...)`.
2141 *
template <
int dim,
int degree,
int n_po
ints_1d>
2142 *
void EulerOperator<dim, degree, n_points_1d>::local_apply_boundary_face(
2146 *
const std::pair<unsigned int, unsigned int> &face_range)
const
2150 *
for (
unsigned int face = face_range.first; face < face_range.second; ++face)
2155 *
for (
const unsigned int q : phi.quadrature_point_indices())
2157 * const auto w_m = phi.get_value(q);
2158 *
const auto normal = phi.normal_vector(q);
2160 *
auto rho_u_dot_n = w_m[1] * normal[0];
2161 *
for (
unsigned int d = 1;
d < dim; ++
d)
2162 * rho_u_dot_n += w_m[1 + d] * normal[d];
2164 *
bool at_outflow =
false;
2168 *
if (wall_boundaries.find(boundary_id) != wall_boundaries.end())
2171 *
for (
unsigned int d = 0;
d < dim; ++
d)
2172 * w_p[d + 1] = w_m[d + 1] - 2. * rho_u_dot_n * normal[d];
2173 * w_p[dim + 1] = w_m[dim + 1];
2175 *
else if (inflow_boundaries.find(boundary_id) !=
2176 * inflow_boundaries.end())
2178 * evaluate_function(*inflow_boundaries.find(boundary_id)->second,
2179 * phi.quadrature_point(q));
2180 *
else if (subsonic_outflow_boundaries.find(boundary_id) !=
2181 * subsonic_outflow_boundaries.end())
2184 * w_p[dim + 1] = evaluate_function(
2185 * *subsonic_outflow_boundaries.find(boundary_id)->second,
2186 * phi.quadrature_point(q),
2188 * at_outflow =
true;
2192 * ExcMessage(
"Unknown boundary id, did "
2193 *
"you set a boundary condition for "
2194 *
"this part of the domain boundary?"));
2196 *
auto flux = euler_numerical_flux<dim>(w_m, w_p, normal);
2199 *
for (
unsigned int v = 0; v < VectorizedArray<Number>::size(); ++v)
2201 *
if (rho_u_dot_n[v] < -1e-12)
2202 *
for (
unsigned int d = 0;
d < dim; ++
d)
2203 * flux[d + 1][v] = 0.;
2206 * phi.submit_value(-flux, q);
2217 * The next function implements the inverse mass
matrix operation. The
2218 * algorithms and rationale have been discussed extensively in the
2219 * introduction, so we here limit ourselves to the technicalities of the
2221 * operations as the forward evaluation of the mass
matrix, except with a
2222 * different interpolation
matrix, representing the inverse @f$S^{-1}@f$
2223 * factors. These represent a change of basis from the specified basis (in
2224 *
this case, the Lagrange basis in the points of the Gauss--Lobatto
2225 * quadrature formula) to the Lagrange basis in the points of the Gauss
2226 * quadrature formula. In the latter basis, we can
apply the inverse of the
2227 *
point-wise `JxW` factor, i.e., the quadrature weight times the
2228 *
determinant of the Jacobian of the mapping from reference to real
2229 * coordinates. Once
this is done, the basis is changed back to the nodal
2230 * Gauss-Lobatto basis again. All of these operations are done by the
2231 * `
apply()` function below. What we need to provide is the local fields to
2232 * operate on (which we extract from the global vector by an
FEEvaluation
2233 *
object) and write the results back to the destination vector of the mass
2238 * One thing to note is that we added two integer arguments (that are
2241 * only have one) and the
second being 1 to make the quadrature formula
2242 * selection. As we use the quadrature formula 0
for the over-integration of
2243 * nonlinear terms, we use the formula 1 with the
default @f$p+1@f$ (or
2244 * `fe_degree+1` in terms of the variable name) points
for the mass
2245 *
matrix. This leads to square contributions to the mass
matrix and ensures
2246 * exact integration, as explained in the introduction.
2249 *
template <
int dim,
int degree,
int n_po
ints_1d>
2250 *
void EulerOperator<dim, degree, n_points_1d>::local_apply_inverse_mass_matrix(
2254 *
const std::pair<unsigned int, unsigned int> &cell_range)
const
2260 *
for (
unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
2263 * phi.read_dof_values(src);
2265 * inverse.apply(phi.begin_dof_values(), phi.begin_dof_values());
2267 * phi.set_dof_values(dst);
2276 * <a name=
"step_67-Theapplyandrelatedfunctions"></a>
2277 * <h4>The
apply() and related functions</h4>
2281 * We now come to the function which implements the evaluation of the Euler
2282 * operator as a whole, i.e., @f$\mathcal M^{-1} \mathcal L(t, \mathbf{
w})@f$,
2283 * calling into the local evaluators presented above. The steps should be
2284 * clear from the previous code. One thing to note is that we need to adjust
2285 * the time in the functions we have associated with the various parts of
2286 * the boundary, in order to be consistent with the equation in
case the
2288 * perform the cell and face integrals, including the necessary ghost
data
2289 * exchange in the `src` vector. The seventh argument to the function,
2290 * `
true`, specifies that we want to zero the `dst` vector as part of the
2291 *
loop, before we start accumulating integrals into it. This variant is
2292 * preferred over explicitly calling `dst = 0.;` before the
loop as the
2293 * zeroing operation is done on a subrange of the vector in parts that are
2294 * written by the integrals nearby. This enhances
data locality and allows
2295 *
for caching, saving one roundtrip of vector
data to main memory and
2296 * enhancing performance. The last two arguments to the
loop determine which
2298 * one faces, typical of
first-order hyperbolic problems, and since we have
2299 * a nodal basis with nodes at the reference element surface, we only need
2300 * to exchange those parts. This again saves precious memory bandwidth.
2304 * Once the spatial
operator @f$\mathcal L@f$ is applied, we need to make a
2307 * is cheaper than the full loop as access only goes to the degrees of
2308 * freedom associated with the locally owned cells, which is simply the
2309 * locally owned degrees of freedom for DG discretizations. Thus, no ghost
2310 * exchange is needed here.
2314 * Around all these functions, we put timer scopes to record the
2315 * computational time for statistics about the contributions of the various
2319 * template <
int dim,
int degree,
int n_points_1d>
2320 *
void EulerOperator<dim, degree, n_points_1d>::apply(
2321 * const
double current_time,
2328 *
for (
auto &i : inflow_boundaries)
2329 * i.
second->set_time(current_time);
2330 *
for (
auto &i : subsonic_outflow_boundaries)
2331 * i.
second->set_time(current_time);
2333 *
data.loop(&EulerOperator::local_apply_cell,
2334 * &EulerOperator::local_apply_face,
2335 * &EulerOperator::local_apply_boundary_face,
2347 *
data.cell_loop(&EulerOperator::local_apply_inverse_mass_matrix,
2358 * Let us move to the function that does an entire stage of a Runge--Kutta
2359 * update. It calls EulerOperator::apply() followed by some updates
2360 * to the vectors, namely `next_ri = solution + factor_ai * k_i` and
2361 * `solution += factor_solution * k_i`. Rather than performing these
2362 * steps through the vector interfaces, we here present an alternative
2363 * strategy that is faster on cache-based architectures. As the memory
2364 * consumed by the vectors is often much larger than what fits into caches,
2365 * the
data has to effectively come from the slow RAM memory. The situation
2366 * can be improved by loop fusion, i.e., performing both the updates to
2367 * `next_ki` and `solution` within a single sweep. In that case, we would
2368 * read the two vectors `rhs` and `solution` and write into `next_ki` and
2369 * `solution`, compared to at least 4 reads and two writes in the baseline
2370 * case. Here, we go one step further and perform the loop immediately when
2371 * the mass matrix inversion has finished on a part of the
2372 * vector.
MatrixFree::cell_loop() provides a mechanism to attach an
2373 * `
std::function` both before the loop over cells
first touches a vector
2374 * entry (which we do not use here, but is e.g. used for zeroing the vector)
2375 * and a
second `
std::function` to be called after the loop last touches
2376 * an entry. The callback is in form of a range over the given vector (in
2377 * terms of the local index numbering in the
MPI universe) that can be
2378 * addressed by `local_element()` functions.
2382 * For this
second callback, we create a lambda that works on a range and
2383 * write the respective update on this range. Ideally, we would add the
2385 * compiler to SIMD parallelize this loop (which means in practice that we
2386 * ensure that there is no overlap, also called aliasing, between the index
2387 * ranges of the pointers we use inside the loops). It turns out that at the
2388 * time of this writing, GCC 7.2 fails to compile an OpenMP pragma inside a
2389 * lambda function, so we comment this pragma out below. If your compiler is
2390 * newer, you should be able to uncomment these lines again.
2394 * Note that we select a different code path for the last
2395 * Runge--Kutta stage when we do not need to update the `next_ri`
2396 * vector. This strategy gives a considerable speedup. Whereas the inverse
2397 * mass matrix and vector updates take more than 60% of the computational
2398 * time with default vector updates on a 40-core machine, the percentage is
2399 * around 35% with the more optimized variant. In other words, this is a
2400 * speedup of around a third.
2403 * template <
int dim,
int degree,
int n_points_1d>
2404 *
void EulerOperator<dim, degree, n_points_1d>::perform_stage(
2405 * const Number current_time,
2406 * const Number factor_solution,
2407 * const Number factor_ai,
2416 *
for (
auto &i : inflow_boundaries)
2417 * i.
second->set_time(current_time);
2418 *
for (
auto &i : subsonic_outflow_boundaries)
2419 * i.
second->set_time(current_time);
2421 *
data.loop(&EulerOperator::local_apply_cell,
2422 * &EulerOperator::local_apply_face,
2423 * &EulerOperator::local_apply_boundary_face,
2436 * &EulerOperator::local_apply_inverse_mass_matrix,
2440 * std::function<
void(
const unsigned int,
const unsigned int)>(),
2441 * [&](
const unsigned int start_range,
const unsigned int end_range) {
2442 *
const Number ai = factor_ai;
2443 *
const Number bi = factor_solution;
2444 *
if (ai == Number())
2447 *
for (
unsigned int i = start_range; i < end_range; ++i)
2449 *
const Number k_i = next_ri.local_element(i);
2450 *
const Number sol_i = solution.local_element(i);
2451 * solution.local_element(i) = sol_i + bi * k_i;
2457 *
for (
unsigned int i = start_range; i < end_range; ++i)
2459 *
const Number k_i = next_ri.local_element(i);
2460 *
const Number sol_i = solution.local_element(i);
2461 * solution.local_element(i) = sol_i + bi * k_i;
2462 * next_ri.local_element(i) = sol_i + ai * k_i;
2473 * Having discussed the implementation of the
functions that deal with
2474 * advancing the solution by one time step, let us now move to
functions
2475 * that implement other, ancillary operations. Specifically, these are
2476 *
functions that compute projections, evaluate errors, and compute the speed
2477 * of information transport on a cell.
2483 * elements where there is no need to set up and solve a linear system, as
2485 * code here, besides a small speedup of
this non-critical operation, is that
2486 * it shows additional functionality provided by
2491 * The projection operation works as follows: If we denote the
matrix of
2492 * shape
functions evaluated at quadrature points by @f$S@f$, the projection on
2493 * cell @f$K@f$ is an operation of the form @f$\underbrace{S J^
K S^\mathrm
2494 * T}_{\mathcal M^
K} \mathbf{
w}^
K = S J^
K
2495 * \tilde{\mathbf{
w}}(\mathbf{x}_q)_{q=1:n_q}@f$, where @f$J^
K@f$ is the
diagonal
2497 * weight (JxW), @f$\mathcal M^
K@f$ is the cell-wise mass
matrix, and
2498 * @f$\tilde{\mathbf{
w}}(\mathbf{x}_q)_{q=1:n_q}@f$ is the evaluation of the
2499 * field to be projected onto quadrature points. (In reality the
matrix @f$S@f$
2500 * has additional structure through the tensor product, as explained in the
2501 * introduction.) This system can now equivalently be written as
2502 * @f$\mathbf{
w}^
K = \left(S J^K S^\mathrm T\right)^{-1} S J^
K
2503 * \tilde{\mathbf{
w}}(\mathbf{x}_q)_{q=1:n_q} = S^{-\mathrm T}
2504 * \left(J^K\right)^{-1} S^{-1} S J^
K
2505 * \tilde{\mathbf{
w}}(\mathbf{x}_q)_{q=1:n_q}@f$. Now, the term @f$S^{-1} S@f$ and
2506 * then @f$\left(J^K\right)^{-1} J^
K@f$ cancel, resulting in the
final
2507 * expression @f$\mathbf{
w}^
K = S^{-\mathrm T}
2508 * \tilde{\mathbf{
w}}(\mathbf{x}_q)_{q=1:n_q}@f$. This operation is
2511 * The name is derived from the fact that
this projection is simply
2512 * the multiplication by @f$S^{-\mathrm T}@f$, a basis change from the
2513 * nodal basis in the points of the Gaussian quadrature to the given finite
2515 * the result into the vector, overwriting previous content, rather than
2516 * accumulating the results as typical in integration tasks -- we can do
2517 * this because every vector entry has contributions from only a single
2518 * cell for discontinuous Galerkin discretizations.
2521 * template <
int dim,
int degree,
int n_points_1d>
2522 *
void EulerOperator<dim, degree, n_points_1d>::project(
2529 * solution.zero_out_ghost_values();
2537 * inverse.transform_from_q_points_to_basis(dim + 2,
2538 * phi.begin_dof_values(),
2539 * phi.begin_dof_values());
2540 * phi.set_dof_values(solution);
2548 * The next function again repeats functionality also provided by the
2550 * the
explicit code to highlight how the vectorization across several cells
2551 * works and how to accumulate results via that interface: Recall that each
2552 * <i>lane</i> of the
vectorized array holds
data from a different cell. By
2553 * the
loop over all cell batches that are owned by the current
MPI process,
2555 * we would need to further go on and
sum across the entries in the SIMD
2556 * array. However, such a procedure is not stable as the SIMD array could in
2557 * fact not hold
valid data for all its lanes. This happens when the number
2558 * of locally owned cells is not a multiple of the SIMD width. To avoid
2560 * the
data. While one could imagine that we could make it work by simply
2561 * setting the empty lanes to zero (and thus, not contribute to a sum), the
2562 * situation is more complicated than that: What
if we were to compute a
2563 * velocity out of the momentum? Then, we would need to divide by the
2564 * density, which is zero -- the result would consequently be NaN and
2565 * contaminate the result. This trap is avoided by accumulating the results
2566 * from the
valid SIMD range as we
loop through the cell batches,
using the
2569 * most cells, but can be less on the last cell batch if the number of cells
2570 * has a remainder compared to the SIMD width.
2573 * template <
int dim,
int degree,
int n_points_1d>
2574 *
std::array<
double, 3> EulerOperator<dim, degree, n_points_1d>::compute_errors(
2579 *
double errors_squared[3] = {};
2582 *
for (
unsigned int cell = 0; cell <
data.n_cell_batches(); ++cell)
2587 *
for (
const unsigned int q : phi.quadrature_point_indices())
2589 * const auto error =
2590 * evaluate_function(function, phi.quadrature_point(q)) -
2592 *
const auto JxW = phi.JxW(q);
2594 * local_errors_squared[0] += error[0] * error[0] * JxW;
2595 *
for (
unsigned int d = 0;
d < dim; ++
d)
2596 * local_errors_squared[1] += (error[d + 1] * error[d + 1]) * JxW;
2597 * local_errors_squared[2] += (error[dim + 1] * error[dim + 1]) * JxW;
2599 *
for (
unsigned int v = 0; v <
data.n_active_entries_per_cell_batch(cell);
2601 *
for (
unsigned int d = 0;
d < 3; ++
d)
2602 * errors_squared[d] += local_errors_squared[d][v];
2607 * std::array<double, 3> errors;
2608 *
for (
unsigned int d = 0;
d < 3; ++
d)
2609 * errors[d] =
std::sqrt(errors_squared[d]);
2618 * This
final function of the EulerOperator
class is used to estimate the
2619 * transport speed, scaled by the mesh
size, that is relevant
for setting
2620 * the time step
size in the
explicit time integrator. In the Euler
2621 * equations, there are two speeds of transport, namely the convective
2622 * velocity @f$\mathbf{u}@f$ and the propagation of sound waves with sound
2623 * speed @f$c = \
sqrt{\
gamma p/\rho}@f$ relative to the medium moving at
2624 * velocity @f$\mathbf u@f$.
2628 * In the formula
for the time step
size, we are interested not by
2629 * these absolute speeds, but by the amount of time it takes
for
2630 * information to cross a single cell. For information transported along with
2631 * the medium, @f$\mathbf u@f$ is scaled by the mesh
size,
2632 * so an estimate of the maximal velocity can be obtained by computing
2633 * @f$\|J^{-\mathrm T} \mathbf{u}\|_\infty@f$, where @f$J@f$ is the Jacobian of the
2634 * transformation from real to the reference domain. Note that
2636 * Jacobian, representing the metric term from real to reference
2637 * coordinates, so we do not need to
transpose it again. We store this limit
2638 * in the variable `convective_limit` in the code below.
2642 * The sound propagation is isotropic, so we need to take mesh sizes in any
2643 * direction into account. The appropriate mesh
size scaling is then given
2644 * by the minimal singular value of @f$J@f$ or, equivalently, the maximal
2645 * singular value of @f$J^{-1}@f$. Note that one could
approximate this quantity
2646 * by the minimal distance between vertices of a cell when ignoring curved
2647 * cells. To get the maximal singular
value of the Jacobian, the
general
2648 * strategy would be some LAPACK function. Since all we need here is an
2649 * estimate, we can avoid the hassle of decomposing a tensor of
2651 * eigenvalue function without vectorization, and instead use a few
2652 * iterations (five in the code below) of the power method applied to
2653 * @f$J^{-1}J^{-\mathrm T}@f$. The speed of convergence of
this method depends
2654 * on the ratio of the largest to the next largest eigenvalue and the
2655 *
initial guess, which is the vector of all ones. This might suggest that
2656 * we get slow convergence on cells close to a cube shape where all
2657 * lengths are almost the same. However,
this slow convergence means that
2658 * the result will sit between the two largest singular
values, which both
2659 * are close to the maximal
value anyway. In all other cases, convergence
2660 * will be quick. Thus, we can merely hardcode 5 iterations here and be
2661 * confident that the result is good.
2664 *
template <
int dim,
int degree,
int n_po
ints_1d>
2665 *
double EulerOperator<dim, degree, n_points_1d>::compute_cell_transport_speed(
2669 * Number max_transport = 0;
2672 *
for (
unsigned int cell = 0; cell <
data.n_cell_batches(); ++cell)
2677 *
for (
const unsigned int q : phi.quadrature_point_indices())
2679 * const auto solution = phi.get_value(q);
2680 *
const auto velocity = euler_velocity<dim>(solution);
2681 *
const auto pressure = euler_pressure<dim>(solution);
2683 *
const auto inverse_jacobian = phi.inverse_jacobian(q);
2684 *
const auto convective_speed = inverse_jacobian * velocity;
2686 *
for (
unsigned int d = 0;
d < dim; ++
d)
2687 * convective_limit =
2690 *
const auto speed_of_sound =
2691 *
std::sqrt(gamma * pressure * (1. / solution[0]));
2694 *
for (
unsigned int d = 0;
d < dim; ++
d)
2695 * eigenvector[d] = 1.;
2696 *
for (
unsigned int i = 0; i < 5; ++i)
2698 * eigenvector =
transpose(inverse_jacobian) *
2699 * (inverse_jacobian * eigenvector);
2701 *
for (
unsigned int d = 0;
d < dim; ++
d)
2702 * eigenvector_norm =
2704 * eigenvector /= eigenvector_norm;
2706 *
const auto jac_times_ev = inverse_jacobian * eigenvector;
2707 *
const auto max_eigenvalue =
std::sqrt(
2708 * (jac_times_ev * jac_times_ev) / (eigenvector * eigenvector));
2711 * max_eigenvalue * speed_of_sound + convective_limit);
2716 * Similarly to the previous function, we must make sure to accumulate
2717 * speed only on the
valid cells of a cell batch.
2720 *
for (
unsigned int v = 0; v <
data.n_active_entries_per_cell_batch(cell);
2722 * max_transport =
std::max(max_transport, local_max[v]);
2727 *
return max_transport;
2735 * <a name=
"step_67-TheEulerProblemclass"></a>
2736 * <h3>The EulerProblem
class</h3>
2740 * This
class combines the EulerOperator class with the time integrator and
2742 * actually
run the simulations of the Euler problem.
2746 * The member variables are a
triangulation, a finite element, a mapping (to
2747 * create high-order curved surfaces, see
e.g. @ref step_10
"step-10"), and a
DoFHandler to
2748 * describe the degrees of freedom. In addition, we keep an instance of the
2749 * EulerOperator described above around, which will
do all heavy lifting in
2750 * terms of integrals, and some parameters
for time integration like the
2751 * current time or the time step
size.
2755 * Furthermore, we use a PostProcessor instance to write some additional
2756 * information to the output file, in similarity to what was done in
2757 * @ref step_33
"step-33". The
interface of the
DataPostprocessor class is intuitive,
2758 * requiring us to provide information about what needs to be evaluated
2759 * (typically only the
values of the solution, except
for the Schlieren plot
2760 * that we only enable in 2
d where it makes sense), and the names of what
2761 * gets evaluated. Note that it would also be possible to extract most
2762 * information by calculator tools within visualization programs such as
2763 * ParaView, but it is so much more convenient to
do it already when writing
2767 * template <int dim>
2768 *
class EulerProblem
2776 *
void make_grid_and_dofs();
2778 *
void output_results(
const unsigned int result_number);
2784 * #ifdef DEAL_II_WITH_P4EST
2796 * EulerOperator<dim, fe_degree, n_q_points_1d> euler_operator;
2798 *
double time, time_step;
2807 * std::vector<
Vector<double>> &computed_quantities)
const override;
2809 *
virtual std::vector<std::string>
get_names()
const override;
2811 *
virtual std::vector<
2818 *
const bool do_schlieren_plot;
2824 *
template <
int dim>
2825 * EulerProblem<dim>::Postprocessor::Postprocessor()
2826 * : do_schlieren_plot(dim == 2)
2833 * For the main evaluation of the field variables, we
first check that the
2834 * lengths of the arrays
equal the expected
values (the lengths `2*dim+4` or
2835 * `2*dim+5` are derived from the sizes of the names we specify in the
2836 * get_names() function below). Then we
loop over all evaluation points and
2837 * fill the respective information: First we fill the primal solution
2838 * variables of density @f$\rho@f$, momentum @f$\rho \mathbf{u}@f$ and energy @f$E@f$,
2839 * then we compute the derived velocity @f$\mathbf u@f$, the pressure @f$p@f$, the
2840 * speed of sound @f$c=\
sqrt{\
gamma p / \rho}@f$, as well as the Schlieren plot
2841 * showing @f$s = |\nabla \rho|^2@f$ in
case it is enabled. (See @ref step_69
"step-69" for
2842 * another example where we create a Schlieren plot.)
2845 *
template <
int dim>
2846 *
void EulerProblem<dim>::Postprocessor::evaluate_vector_field(
2850 *
const unsigned int n_evaluation_points = inputs.solution_values.size();
2852 *
if (do_schlieren_plot ==
true)
2853 *
Assert(inputs.solution_gradients.size() == n_evaluation_points,
2854 * ExcInternalError());
2856 *
Assert(computed_quantities.size() == n_evaluation_points,
2857 * ExcInternalError());
2858 *
Assert(inputs.solution_values[0].size() == dim + 2, ExcInternalError());
2860 * dim + 2 + (do_schlieren_plot ==
true ? 1 : 0),
2861 * ExcInternalError());
2863 *
for (
unsigned int p = 0; p < n_evaluation_points; ++p)
2866 *
for (
unsigned int d = 0;
d < dim + 2; ++
d)
2867 * solution[d] = inputs.solution_values[p](d);
2869 *
const double density = solution[0];
2871 *
const double pressure = euler_pressure<dim>(solution);
2873 *
for (
unsigned int d = 0;
d < dim; ++
d)
2874 * computed_quantities[p](d) = velocity[
d];
2875 * computed_quantities[p](dim) = pressure;
2876 * computed_quantities[p](dim + 1) =
std::sqrt(gamma * pressure / density);
2878 *
if (do_schlieren_plot ==
true)
2879 * computed_quantities[p](dim + 2) =
2880 * inputs.solution_gradients[p][0] * inputs.solution_gradients[p][0];
2886 *
template <
int dim>
2887 * std::vector<std::string> EulerProblem<dim>::Postprocessor::get_names() const
2889 * std::vector<std::string> names;
2890 *
for (
unsigned int d = 0;
d < dim; ++
d)
2891 * names.emplace_back(
"velocity");
2892 * names.emplace_back(
"pressure");
2893 * names.emplace_back(
"speed_of_sound");
2895 *
if (do_schlieren_plot ==
true)
2896 * names.emplace_back(
"schlieren_plot");
2905 * For the interpretation of quantities, we have
scalar density, energy,
2906 * pressure, speed of sound, and the Schlieren plot, and vectors
for the
2907 * momentum and the velocity.
2910 *
template <
int dim>
2911 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2912 * EulerProblem<dim>::Postprocessor::get_data_component_interpretation() const
2914 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
2916 *
for (
unsigned int d = 0;
d < dim; ++
d)
2917 * interpretation.push_back(
2922 *
if (do_schlieren_plot ==
true)
2923 * interpretation.push_back(
2926 *
return interpretation;
2933 * With respect to the necessary update flags, we only need the
values for
2934 * all quantities but the Schlieren plot, which is based on the density
2938 *
template <
int dim>
2939 *
UpdateFlags EulerProblem<dim>::Postprocessor::get_needed_update_flags() const
2941 *
if (do_schlieren_plot ==
true)
2951 * The constructor
for this class is unsurprising: We set up a
parallel
2952 *
triangulation based on the `MPI_COMM_WORLD` communicator, a vector finite
2953 * element with `dim+2` components for density, momentum, and energy, a
2954 * high-order mapping of the same degree as the underlying finite element,
2955 * and initialize the time and time step to zero.
2958 * template <int dim>
2959 * EulerProblem<dim>::EulerProblem()
2960 * : pcout(std::cout, Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0)
2961 * #ifdef DEAL_II_WITH_P4EST
2964 * , fe(FE_DGQ<dim>(fe_degree) ^ (dim + 2))
2965 * , mapping(fe_degree)
2966 * , dof_handler(triangulation)
2967 * , timer(pcout, TimerOutput::never, TimerOutput::wall_times)
2968 * , euler_operator(timer)
2977 * As a mesh,
this tutorial program implements two options, depending on the
2978 * global variable `testcase`: For the analytical variant (`testcase==0`),
2979 * the domain is @f$(0, 10) \times (-5, 5)@f$, with Dirichlet boundary
2980 * conditions (inflow) all around the domain. For `testcase==1`, we set the
2981 * domain to a
cylinder in a rectangular box, derived from the flow past
2982 *
cylinder testcase
for incompressible viscous flow by Schäfer and
2983 * Turek (1996). Here, we have a larger variety of boundaries. The inflow
2984 * part at the left of the channel is given the inflow type,
for which we
2985 * choose a
constant inflow profile, whereas we set a subsonic outflow at
2986 * the right. For the boundary around the
cylinder (boundary
id equal to 2)
2987 * as well as the channel walls (boundary
id equal to 3) we use the wall
2988 * boundary type, which is no-normal flow. Furthermore,
for the 3
d cylinder
2989 * we also add a gravity force in vertical direction. Having the base mesh
2990 * in place (including the manifolds set by
2992 * specified number of global refinements, create the unknown numbering from
2994 * initialization of the EulerOperator.
2997 *
template <
int dim>
2998 *
void EulerProblem<dim>::make_grid_and_dofs()
3005 *
for (
unsigned int d = 1;
d < dim; ++
d)
3006 * lower_left[d] = -5;
3009 * upper_right[0] = 10;
3010 *
for (
unsigned int d = 1;
d < dim; ++
d)
3011 * upper_right[d] = 5;
3018 * euler_operator.set_inflow_boundary(
3019 * 0, std::make_unique<ExactSolution<dim>>(0));
3029 * euler_operator.set_inflow_boundary(
3030 * 0, std::make_unique<ExactSolution<dim>>(0));
3031 * euler_operator.set_subsonic_outflow_boundary(
3032 * 1, std::make_unique<ExactSolution<dim>>(0));
3034 * euler_operator.set_wall_boundary(2);
3035 * euler_operator.set_wall_boundary(3);
3038 * euler_operator.set_body_force(
3040 * std::vector<double>({0., 0., -0.2})));
3051 * dof_handler.distribute_dofs(fe);
3053 * euler_operator.reinit(mapping, dof_handler);
3054 * euler_operator.initialize_vector(solution);
3058 * In the following, we output some statistics about the problem. Because we
3059 * often
end up with quite large
numbers of cells or degrees of freedom, we
3060 * would like to print them with a comma to separate each set of three
3061 * digits. This can be done via
"locales", although the way
this works is
3062 * not particularly intuitive. @ref step_32
"step-32" explains
this in slightly more
3066 * std::locale s = pcout.get_stream().getloc();
3067 * pcout.get_stream().imbue(std::locale(
""));
3068 * pcout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
3069 * <<
" ( = " << (dim + 2) <<
" [vars] x "
3073 * pcout.get_stream().imbue(s);
3080 * For output, we
first let the Euler
operator compute the errors of the
3081 * numerical results. More precisely, we compute the error against the
3082 * analytical result
for the analytical solution
case, whereas we compute
3083 * the deviation against the background field with
constant density and
3084 * energy and
constant velocity in @f$x@f$ direction
for the
second test
case.
3088 * The next step is to create output. This is similar to what is done in
3089 * @ref step_33
"step-33": We let the postprocessor defined above control most of the
3090 * output, except
for the primal field that we write directly. For the
3091 * analytical solution test
case, we also perform another projection of the
3092 * analytical solution and print the difference between that field and the
3093 * numerical solution. Once we have defined all quantities to be written, we
3094 * build the patches
for output. Similarly to @ref step_65
"step-65", we create a
3095 * high-order VTK output by setting the appropriate flag, which enables us
3096 * to visualize fields of high polynomial degrees. Finally, we call the
3098 * to the given file name. This function uses special
MPI parallel write
3099 * facilities, which are typically more optimized
for parallel file systems
3100 * than the standard library
's `std::ofstream` variants used in most other
3101 * tutorial programs. A particularly nice feature of the
3102 * `write_vtu_in_parallel()` function is the fact that it can combine output
3103 * from all MPI ranks into a single file, making it unnecessary to have a
3104 * central record of all such files (namely, the "pvtu" file).
3108 * For parallel programs, it is often instructive to look at the partitioning
3109 * of cells among processors. To this end, one can pass a vector of numbers
3110 * to DataOut::add_data_vector() that contains as many entries as the
3111 * current processor has active cells; these numbers should then be the
3112 * rank of the processor that owns each of these cells. Such a vector
3113 * could, for example, be obtained from
3114 * GridTools::get_subdomain_association(). On the other hand, on each MPI
3115 * process, DataOut will only read those entries that correspond to locally
3116 * owned cells, and these of course all have the same value: namely, the rank
3117 * of the current process. What is in the remaining entries of the vector
3118 * doesn't actually matter, and so we can just get away with a cheap trick: We
3120 * with the rank of the current
MPI process. The key is that on each process,
3121 * only the entries corresponding to the locally owned cells will be read,
3122 * ignoring the (wrong)
values in other entries. The fact that every process
3123 * submits a vector in which the correct subset of entries is correct is all
3124 * that is necessary.
3128 * @note As of 2023, Visit 3.3.3 can still not deal with higher-order cells.
3129 * Rather, it simply reports that there is no
data to show. To view the
3130 * results of
this program with Visit, you will want to comment out the
3131 * line that sets `flags.write_higher_order_cells =
true;`. On the other
3132 * hand, Paraview is able to understand VTU files with higher order cells
3136 *
template <
int dim>
3137 *
void EulerProblem<dim>::output_results(
const unsigned int result_number)
3139 *
const std::array<double, 3> errors =
3140 * euler_operator.compute_errors(ExactSolution<dim>(time), solution);
3141 *
const std::string quantity_name = testcase == 0 ?
"error" :
"norm";
3143 * pcout <<
"Time:" << std::setw(8) << std::setprecision(3) << time
3144 * <<
", dt: " << std::setw(8) << std::setprecision(2) << time_step
3145 * <<
", " << quantity_name <<
" rho: " << std::setprecision(4)
3146 * << std::setw(10) << errors[0] <<
", rho * u: " << std::setprecision(4)
3147 * << std::setw(10) << errors[1] <<
", energy:" << std::setprecision(4)
3148 * << std::setw(10) << errors[2] << std::endl;
3153 * Postprocessor postprocessor;
3158 * data_out.set_flags(flags);
3160 * data_out.attach_dof_handler(dof_handler);
3162 * std::vector<std::string> names;
3163 * names.emplace_back(
"density");
3164 *
for (
unsigned int d = 0;
d < dim; ++
d)
3165 * names.emplace_back(
"momentum");
3166 * names.emplace_back(
"energy");
3168 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
3170 * interpretation.push_back(
3172 *
for (
unsigned int d = 0;
d < dim; ++
d)
3173 * interpretation.push_back(
3175 * interpretation.push_back(
3178 * data_out.add_data_vector(dof_handler, solution, names, interpretation);
3180 * data_out.add_data_vector(solution, postprocessor);
3183 *
if (testcase == 0 && dim == 2)
3185 * reference.
reinit(solution);
3186 * euler_operator.project(ExactSolution<dim>(time), reference);
3187 * reference.sadd(-1., 1, solution);
3188 * std::vector<std::string> names;
3189 * names.emplace_back(
"error_density");
3190 *
for (
unsigned int d = 0;
d < dim; ++
d)
3191 * names.emplace_back(
"error_momentum");
3192 * names.emplace_back(
"error_energy");
3194 * std::vector<DataComponentInterpretation::DataComponentInterpretation>
3196 * interpretation.push_back(
3198 *
for (
unsigned int d = 0;
d < dim; ++
d)
3199 * interpretation.push_back(
3201 * interpretation.push_back(
3204 * data_out.add_data_vector(dof_handler,
3212 * data_out.add_data_vector(mpi_owner,
"owner");
3214 * data_out.build_patches(mapping,
3218 *
const std::string filename =
3220 * data_out.write_vtu_in_parallel(filename, MPI_COMM_WORLD);
3228 * The EulerProblem::run() function puts all pieces together. It starts off
3229 * by calling the function that creates the mesh and sets up
data structures,
3230 * and then initializing the time integrator and the two temporary vectors of
3231 * the low-storage integrator. We call these vectors `rk_register_1` and
3232 * `rk_register_2`, and use the
first vector to represent the quantity
3233 * @f$\mathbf{r}_i@f$ and the
second one
for @f$\mathbf{k}_i@f$ in the formulas
for
3234 * the Runge--Kutta scheme outlined in the introduction. Before we start the
3235 * time
loop, we compute the time step
size by the
3236 * `EulerOperator::compute_cell_transport_speed()` function. For reasons of
3237 * comparison, we compare the result obtained there with the minimal mesh
3238 *
size and print them to screen. For velocities and speeds of sound close
3239 * to unity as in
this tutorial program, the predicted effective mesh
size
3240 * will be close, but they could vary
if scaling were different.
3243 *
template <
int dim>
3244 *
void EulerProblem<dim>::run()
3248 *
const unsigned int n_vect_bits = 8 *
sizeof(Number) * n_vect_number;
3250 * pcout <<
"Running with "
3252 * <<
" MPI processes" << std::endl;
3253 * pcout <<
"Vectorization over " << n_vect_number <<
' '
3254 * << (std::is_same_v<Number, double> ?
"doubles" :
"floats") <<
" = "
3255 * << n_vect_bits <<
" bits ("
3260 * make_grid_and_dofs();
3262 *
const LowStorageRungeKuttaIntegrator integrator(lsrk_scheme);
3266 * rk_register_1.
reinit(solution);
3267 * rk_register_2.reinit(solution);
3269 * euler_operator.project(ExactSolution<dim>(time), solution);
3271 *
double min_vertex_distance = std::numeric_limits<double>::max();
3272 *
for (
const auto &cell :
triangulation.active_cell_iterators())
3273 * if (cell->is_locally_owned())
3274 * min_vertex_distance =
3275 *
std::
min(min_vertex_distance, cell->minimum_vertex_distance());
3276 * min_vertex_distance =
3279 * time_step = courant_number * integrator.n_stages() /
3280 * euler_operator.compute_cell_transport_speed(solution);
3281 * pcout <<
"Time step size: " << time_step
3282 * <<
", minimal h: " << min_vertex_distance
3283 * <<
", initial transport scaling: "
3284 * << 1. / euler_operator.compute_cell_transport_speed(solution)
3288 * output_results(0);
3292 * Now we are ready to start the time
loop, which we
run until the time
3293 * has reached the desired
end time. Every 5 time steps, we compute a
new
3294 * estimate
for the time step -- since the solution is nonlinear, it is
3295 * most effective to adapt the
value during the course of the
3296 * simulation. In
case the Courant number was chosen too aggressively, the
3297 * simulation will typically blow up with time step NaN, so that is easy
3298 * to detect here. One thing to note is that roundoff errors might
3299 * propagate to the leading digits due to an interaction of slightly
3300 * different time step selections that in turn lead to slightly different
3301 * solutions. To decrease
this sensitivity, it is common practice to round
3302 * or truncate the time step
size to a few digits,
e.g. 3 in
this case. In
3303 *
case the current time is near the prescribed
'tick' value for output
3304 * (
e.g. 0.02), we also write the output. After the end of the time loop,
3305 * we summarize the computation by printing some statistics, which is
3309 *
unsigned int timestep_number = 0;
3311 *
while (time < final_time - 1e-12)
3313 * ++timestep_number;
3314 *
if (timestep_number % 5 == 0)
3316 * courant_number * integrator.n_stages() /
3318 * euler_operator.compute_cell_transport_speed(solution), 3);
3322 * integrator.perform_time_step(euler_operator,
3330 * time += time_step;
3332 *
if (
static_cast<int>(time / output_tick) !=
3333 *
static_cast<int>((time - time_step) / output_tick) ||
3334 * time >= final_time - 1e-12)
3336 *
static_cast<unsigned int>(std::round(time / output_tick)));
3339 * timer.print_wall_time_statistics(MPI_COMM_WORLD);
3340 * pcout << std::endl;
3349 * The main() function is not surprising and follows what was done in all
3350 * previous
MPI programs: As we run an
MPI program, we need to call `MPI_Init()`
3352 *
Utilities::
MPI::MPI_InitFinalize
data structure. Note that we run the program
3353 * only with
MPI, and set the thread count to 1.
3356 *
int main(
int argc,
char **argv)
3358 *
using namespace Euler_DG;
3359 *
using namespace dealii;
3365 * EulerProblem<dimension> euler_problem;
3366 * euler_problem.run();
3368 *
catch (std::exception &exc)
3370 * std::cerr << std::endl
3372 * <<
"----------------------------------------------------"
3374 * std::cerr <<
"Exception on processing: " << std::endl
3375 * << exc.what() << std::endl
3376 * <<
"Aborting!" << std::endl
3377 * <<
"----------------------------------------------------"
3384 * std::cerr << std::endl
3386 * <<
"----------------------------------------------------"
3388 * std::cerr <<
"Unknown exception!" << std::endl
3389 * <<
"Aborting!" << std::endl
3390 * <<
"----------------------------------------------------"
3398<a name=
"step_67-Results"></a><h1>Results</h1>
3401<a name=
"step_67-Programoutput"></a><h3>Program output</h3>
3404Running the program with the
default settings on a machine with 40 processes
3405produces the following output:
3407Running with 40
MPI processes
3408Vectorization over 8 doubles = 512 bits (AVX512)
3409Number of degrees of freedom: 147,456 ( = 4 [vars] x 1,024 [cells] x 36 [dofs/cell/var] )
3410Time step
size: 0.00689325, minimal h: 0.3125,
initial transport scaling: 0.102759
3412Time: 0, dt: 0.0069, error rho: 2.76e-07, rho * u: 1.259e-06, energy: 2.987e-06
3413Time: 1.01, dt: 0.0069, error rho: 1.37e-06, rho * u: 2.252e-06, energy: 4.153e-06
3414Time: 2.01, dt: 0.0069, error rho: 1.561e-06, rho * u: 2.43e-06, energy: 4.493e-06
3415Time: 3.01, dt: 0.0069, error rho: 1.714e-06, rho * u: 2.591e-06, energy: 4.762e-06
3416Time: 4.01, dt: 0.0069, error rho: 1.843e-06, rho * u: 2.625e-06, energy: 4.985e-06
3417Time: 5.01, dt: 0.0069, error rho: 1.496e-06, rho * u: 1.961e-06, energy: 4.142e-06
3418Time: 6, dt: 0.0083, error rho: 1.007e-06, rho * u: 7.119e-07, energy: 2.972e-06
3419Time: 7, dt: 0.0095, error rho: 9.096e-07, rho * u: 3.786e-07, energy: 2.626e-06
3420Time: 8, dt: 0.0096, error rho: 8.439e-07, rho * u: 3.338e-07, energy: 2.43e-06
3421Time: 9, dt: 0.0096, error rho: 7.822e-07, rho * u: 2.984e-07, energy: 2.248e-06
3422Time: 10, dt: 0.0096, error rho: 7.231e-07, rho * u: 2.666e-07, energy: 2.074e-06
3424+-------------------------------------------+------------------+------------+------------------+
3425| Total wallclock time elapsed | 2.249s 30 | 2.249s | 2.249s 8 |
3427| Section | no. calls |
min time rank |
avg time |
max time rank |
3428+-------------------------------------------+------------------+------------+------------------+
3429| compute errors | 11 | 0.008066s 13 | 0.00952s | 0.01041s 20 |
3430| compute transport speed | 258 | 0.01012s 13 | 0.05392s | 0.08574s 25 |
3431| output | 11 | 0.9597s 13 | 0.9613s | 0.9623s 6 |
3432| rk time stepping total | 1283 | 0.9827s 25 | 1.015s | 1.06s 13 |
3433| rk_stage - integrals L_h | 6415 | 0.8803s 26 | 0.9198s | 0.9619s 14 |
3434| rk_stage - inv mass + vec upd | 6415 | 0.05677s 15 | 0.06487s | 0.07597s 13 |
3435+-------------------------------------------+------------------+------------+------------------+
3438The program output shows that all errors are small. This is due to the fact
3439that we use a relatively fine mesh of @f$32^2@f$ cells with polynomials of degree
34405 for a solution that is smooth. An interesting pattern shows for the time
3441step
size: whereas it is 0.0069 up to time 5, it increases to 0.0096 for later
3442times. The step
size increases once the vortex with some motion on top of the
3443speed of sound (and thus faster propagation) leaves the computational domain
3444between times 5 and 6.5. After that
point, the flow is simply uniform
3445in the same direction, and the maximum velocity of the gas is reduced
3446compared to the previous state where the uniform velocity was overlaid
3447by the vortex. Our time step formula recognizes this effect.
3449The final block of output shows detailed information about the timing
3450of individual parts of the programs; it breaks
this down by showing
3451the time taken by the fastest and the slowest processor, and the
3452average time --
this is often useful in very large computations to
3453find whether there are processors that are consistently overheated
3454(and consequently are throttling their clock speed) or consistently
3455slow
for other reasons.
3456The summary shows that 1283 time steps have been performed
3457in 1.02 seconds (looking at the average time among all
MPI processes),
while
3458the output of 11 files has taken additional 0.96 seconds. Broken down per time
3459step and into the five Runge--Kutta stages, the compute time per evaluation is
34600.16 milliseconds. This high performance is typical of
matrix-free evaluators
3461and a reason why
explicit time integration is very competitive against
3462implicit solvers, especially
for large-
scale simulations. The breakdown of
3463computational times at the
end of the program
run shows that the evaluation of
3464integrals in @f$\mathcal L_h@f$ contributes with around 0.92 seconds and the
3465application of the inverse mass
matrix with 0.06 seconds. Furthermore, the
3466estimation of the transport speed
for the time step
size computation
3467contributes with another 0.05 seconds of compute time.
3469If we use three more levels of global refinement and 9.4 million DoFs in total,
3470the
final statistics are as follows (
for the modified Lax--Friedrichs flux,
3471@f$p=5@f$, and the same system of 40 cores of dual-socket Intel Xeon Gold 6230):
3473+-------------------------------------------+------------------+------------+------------------+
3474| Total wallclock time elapsed | 244.9s 12 | 244.9s | 244.9s 34 |
3476| Section | no. calls |
min time rank |
avg time |
max time rank |
3477+-------------------------------------------+------------------+------------+------------------+
3478| compute errors | 11 | 0.4239s 12 | 0.4318s | 0.4408s 9 |
3479| compute transport speed | 2053 | 3.962s 12 | 6.727s | 10.12s 7 |
3480| output | 11 | 30.35s 12 | 30.36s | 30.37s 9 |
3481| rk time stepping total | 10258 | 201.7s 7 | 205.1s | 207.8s 12 |
3482| rk_stage - integrals L_h | 51290 | 121.3s 6 | 126.6s | 136.3s 16 |
3483| rk_stage - inv mass + vec upd | 51290 | 66.19s 16 | 77.52s | 81.84s 10 |
3484+-------------------------------------------+------------------+------------+------------------+
3487Per time step, the solver now takes 0.02 seconds, about 25 times as long as
3488for the small problem with 147k unknowns. Given that the problem involves 64
3489times as many unknowns, the increase in computing time is not
3490surprising. Since we also do 8 times as many time steps, the compute time
3491should in theory increase by a factor of 512. The actual increase is 205 s /
34921.02 s = 202. This is because the small problem
size cannot fully utilize the
349340 cores due to communication overhead. This becomes clear if we look into the
3494details of the operations done per time step. The evaluation of the
3495differential operator @f$\mathcal L_h@f$ with nearest neighbor communication goes
3496from 0.92 seconds to 127 seconds, i.
e., it increases with a factor of 138. On
3497the other hand, the cost for application of the inverse mass
matrix and the
3498vector updates, which do not need to communicate between the
MPI processes at
3499all, has increased by a factor of 1195. The increase is more than the
3500theoretical factor of 512 because the operation is limited by the bandwidth
3501from RAM memory for the larger
size while for the smaller
size, all vectors
3502fit into the caches of the CPU. The
numbers show that the mass
matrix
3503evaluation and vector update part consume almost 40% of the time spent by the
3504Runge--Kutta stages -- despite using a low-storage Runge--Kutta integrator and
3505merging of vector operations! And despite using over-integration for the
3506@f$\mathcal L_h@f$ operator. For simpler differential operators and more expensive
3507time integrators, the proportion spent in the mass
matrix and vector update
3508part can also reach 70%. If we compute a throughput number in terms of DoFs
3509processed per
second and Runge--Kutta stage, we obtain @f[ \text{throughput} =
3510\frac{n_\mathrm{time steps} n_\mathrm{stages}
3511n_\mathrm{dofs}}{t_\mathrm{compute}} = \frac{10258 \cdot 5 \cdot
35129.4\,\text{MDoFs}}{205s} = 2360\, \text{MDoFs/s} @f] This throughput number is
3513very high, given that simply copying one vector to another one runs at
3514only around 10,000 MDoFs/s.
3516If we go to the next-larger
size with 37.7 million DoFs, the overall
3517simulation time is 2196 seconds, with 1978 seconds spent in the time
3518stepping. The increase in
run time is a factor of 9.3
for the L_h
operator
3519(1179 versus 127 seconds) and a factor of 10.3
for the inverse mass matrix and
3520vector updates (797 vs 77.5 seconds). The reason
for this non-optimal increase
3521in
run time can be traced back to cache effects on the given hardware (with 40
3522MB of L2 cache and 55 MB of L3 cache): While not all of the relevant
data fits
3523into caches for 9.4 million DoFs (one vector takes 75 MB and we have three
3524vectors plus some additional
data in
MatrixFree), there is capacity for one and
3525a half vector nonetheless. Given that modern caches are more sophisticated than
3526the naive least-recently-used strategy (where we would have little re-use as
3527the
data is used in a streaming-like fashion), we can assume that a sizeable
3528fraction of
data can indeed be delivered from caches for the 9.4 million DoFs
3529case. For the larger case, even with optimal caching less than 10 percent of
3530data would fit into caches, with an associated loss in performance.
3533<a name=
"step_67-Convergenceratesfortheanalyticaltestcase"></a><h3>Convergence rates for the analytical test case</h3>
3536For the modified Lax--Friedrichs flux and measuring the error in the momentum
3537variable, we obtain the following convergence table (the rates are very
3538similar for the density and energy variables):
3540<table align=
"center" class=
"doxtable">
3543 <th colspan=
"3"><i>p</i>=2</th>
3544 <th colspan=
"3"><i>p</i>=3</th>
3545 <th colspan=
"3"><i>p</i>=5</th>
3560 <td align=
"right">16</td>
3567 <td align=
"right">2,304</td>
3568 <td align=
"center">1.373e-01</td>
3572 <td align=
"right">64</td>
3576 <td align=
"right">4,096</td>
3577 <td align=
"center">9.130e-02</td>
3579 <td align=
"right">9,216</td>
3580 <td align=
"center">8.899e-03</td>
3584 <td align=
"right">256</td>
3585 <td align=
"right">9,216</td>
3586 <td align=
"center">5.577e-02</td>
3588 <td align=
"right">16,384</td>
3589 <td align=
"center">7.381e-03</td>
3591 <td align=
"right">36,864</td>
3592 <td align=
"center">2.082e-04</td>
3596 <td align=
"right">1024</td>
3597 <td align=
"right">36,864</td>
3598 <td align=
"center">4.724e-03</td>
3600 <td align=
"right">65,536</td>
3601 <td align=
"center">3.072e-04</td>
3603 <td align=
"right">147,456</td>
3604 <td align=
"center">2.625e-06</td>
3608 <td align=
"right">4096</td>
3609 <td align=
"right">147,456</td>
3610 <td align=
"center">6.205e-04</td>
3612 <td align=
"right">262,144</td>
3613 <td align=
"center">1.880e-05</td>
3615 <td align=
"right">589,824</td>
3616 <td align=
"center">3.268e-08</td>
3620 <td align=
"right">16,384</td>
3621 <td align=
"right">589,824</td>
3622 <td align=
"center">8.279e-05</td>
3624 <td align=
"right">1,048,576</td>
3625 <td align=
"center">1.224e-06</td>
3627 <td align=
"right">2,359,296</td>
3628 <td align=
"center">9.252e-10</td>
3632 <td align=
"right">65,536</td>
3633 <td align=
"right">2,359,296</td>
3634 <td align=
"center">1.105e-05</td>
3636 <td align=
"right">4,194,304</td>
3637 <td align=
"center">7.871e-08</td>
3639 <td align=
"right">9,437,184</td>
3640 <td align=
"center">1.369e-10</td>
3644 <td align=
"right">262,144</td>
3645 <td align=
"right">9,437,184</td>
3646 <td align=
"center">1.615e-06</td>
3648 <td align=
"right">16,777,216</td>
3649 <td align=
"center">4.961e-09</td>
3651 <td align=
"right">37,748,736</td>
3652 <td align=
"center">7.091e-11</td>
3657If we
switch to the Harten-Lax-van Leer flux, the results are as follows:
3658<table align=
"center" class=
"doxtable">
3661 <th colspan=
"3"><i>p</i>=2</th>
3662 <th colspan=
"3"><i>p</i>=3</th>
3663 <th colspan=
"3"><i>p</i>=5</th>
3678 <td align=
"right">16</td>
3685 <td align=
"right">2,304</td>
3686 <td align=
"center">1.339e-01</td>
3690 <td align=
"right">64</td>
3694 <td align=
"right">4,096</td>
3695 <td align=
"center">9.037e-02</td>
3697 <td align=
"right">9,216</td>
3698 <td align=
"center">8.849e-03</td>
3702 <td align=
"right">256</td>
3703 <td align=
"right">9,216</td>
3704 <td align=
"center">4.204e-02</td>
3706 <td align=
"right">16,384</td>
3707 <td align=
"center">9.143e-03</td>
3709 <td align=
"right">36,864</td>
3710 <td align=
"center">2.501e-04</td>
3714 <td align=
"right">1024</td>
3715 <td align=
"right">36,864</td>
3716 <td align=
"center">4.913e-03</td>
3718 <td align=
"right">65,536</td>
3719 <td align=
"center">3.257e-04</td>
3721 <td align=
"right">147,456</td>
3722 <td align=
"center">3.260e-06</td>
3726 <td align=
"right">4096</td>
3727 <td align=
"right">147,456</td>
3728 <td align=
"center">7.862e-04</td>
3730 <td align=
"right">262,144</td>
3731 <td align=
"center">1.588e-05</td>
3733 <td align=
"right">589,824</td>
3734 <td align=
"center">2.953e-08</td>
3738 <td align=
"right">16,384</td>
3739 <td align=
"right">589,824</td>
3740 <td align=
"center">1.137e-04</td>
3742 <td align=
"right">1,048,576</td>
3743 <td align=
"center">9.400e-07</td>
3745 <td align=
"right">2,359,296</td>
3746 <td align=
"center">4.286e-10</td>
3750 <td align=
"right">65,536</td>
3751 <td align=
"right">2,359,296</td>
3752 <td align=
"center">1.476e-05</td>
3754 <td align=
"right">4,194,304</td>
3755 <td align=
"center">5.799e-08</td>
3757 <td align=
"right">9,437,184</td>
3758 <td align=
"center">2.789e-11</td>
3762 <td align=
"right">262,144</td>
3763 <td align=
"right">9,437,184</td>
3764 <td align=
"center">2.038e-06</td>
3766 <td align=
"right">16,777,216</td>
3767 <td align=
"center">3.609e-09</td>
3769 <td align=
"right">37,748,736</td>
3770 <td align=
"center">5.730e-11</td>
3775The tables show that we get optimal @f$\mathcal O\left(h^{p+1}\right)@f$
3776convergence rates
for both numerical fluxes. The errors are slightly smaller
3777for the Lax--Friedrichs flux
for @f$p=2@f$, but the picture is reversed
for
3778@f$p=3@f$; in any
case, the differences on
this testcase are relatively
3781For @f$p=5@f$, we reach the roundoff accuracy of @f$10^{-11}@f$ with both
3782fluxes on the finest grids. Also note that the errors are absolute with a
3783domain length of @f$10^2@f$, so relative errors are below @f$10^{-12}@f$. The HLL flux
3784is somewhat better
for the highest degree, which is due to a slight inaccuracy
3785of the Lax--Friedrichs flux: The Lax--Friedrichs flux sets a Dirichlet
3786condition on the solution that leaves the domain, which results in a small
3787artificial reflection, which is accentuated
for the Lax--Friedrichs
3788flux. Apart from that, we see that the influence of the numerical flux is
3789minor, as the polynomial part
inside elements is the main driver of the
3790accucary. The limited influence of the flux also has consequences when trying
3791to approach more challenging setups with the higher-order DG setup: Taking
for
3792example the parameters and grid of @ref step_33
"step-33", we get oscillations (which in turn
3793make density negative and make the solution explode) with both fluxes once the
3794high-mass part comes near the boundary, as opposed to the low-order finite
3795volume case (@f$p=0@f$). Thus, any
case that leads to shocks in the solution
3796necessitates some form of limiting or artificial dissipation. For another
3797alternative, see the @ref step_69
"step-69" tutorial program.
3800<a name=
"step_67-Resultsforflowinchannelaroundcylinderin2D"></a><h3>Results
for flow in channel around
cylinder in 2D</h3>
3803For the test
case of the flow around a
cylinder in a channel, we need to
3804change the
first code line to
3806 constexpr unsigned int testcase = 1;
3808This test
case starts with a background field of a
constant velocity
3810to go around an obstacle in the form of a
cylinder. Since we impose a
3811no-penetration condition on the
cylinder walls, the flow that
3812initially impinges head-on onto to
cylinder has to rearrange,
3813which creates a big sound wave. The following pictures show the pressure at
3814times 0.1, 0.25, 0.5, and 1.0 (top left to bottom right)
for the 2D
case with
38155 levels of global refinement,
using 102,400 cells with polynomial degree of
38165 and 14.7 million degrees of freedom over all 4 solution variables.
3817We clearly see the discontinuity that
3818propagates slowly in the upstream direction and more quickly in downstream
3819direction in the
first snapshot at time 0.1. At time 0.25, the sound wave has
3820reached the top and bottom walls and reflected back to the interior. From the
3821different distances of the reflected waves from lower and upper walls we can
3822see the slight asymmetry of the Schäfer-Turek test
case represented by
3824cylinder compared to below. At later times, the picture is more chaotic with
3825many sound waves all over the place.
3827<table align="center" class="doxtable" style="width:85%">
3846The next picture shows an elevation plot of the pressure at time 1.0 looking
3847from the channel inlet towards the outlet at the same resolution -- here,
3848we can see the large number
3849of reflections. In the figure, two
types of waves are visible. The
3850larger-amplitude waves correspond to various reflections that happened as the
3851initial discontinuity hit the walls, whereas the small-amplitude waves of
3852size similar to the elements correspond to numerical artifacts. They have their
3853origin in the finite resolution of the scheme and appear as the discontinuity
3854travels through elements with high-order polynomials. This effect can be cured
3855by increasing resolution. Apart from this effect, the rich wave structure is
3856the result of the transport accuracy of the high-order DG method.
3860If we run the program with degree 2 and 6 levels of global refinements (410k
3861cells, 14.7M unknowns), we get the following evolution of the pressure
3862(elevation plot, colored by the value of the density):
3866 <iframe width="560" height="315" src="https:
3868 allow="accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture"
3869 allowfullscreen></iframe>
3874With 2 levels of global refinement with 1,600 cells, the mesh and its
3875partitioning on 40
MPI processes looks as follows:
3879When we run the code with 4 levels of global refinements on 40 cores, we get
3880the following output:
3882Running with 40
MPI processes
3883Vectorization over 8 doubles = 512 bits (AVX512)
3884Number of degrees of freedom: 3,686,400 ( = 4 [vars] x 25,600 [cells] x 36 [dofs/cell/var] )
3885Time step
size: 7.39876e-05, minimal h: 0.001875, initial transport scaling: 0.00110294
3887Time: 0, dt: 7.4e-05, norm rho: 4.17e-16, rho * u: 1.629e-16, energy: 1.381e-15
3888Time: 0.05, dt: 6.3e-05, norm rho: 0.02075, rho * u: 0.03801, energy: 0.08772
3889Time: 0.1, dt: 5.9e-05, norm rho: 0.02211, rho * u: 0.04515, energy: 0.08953
3890Time: 0.15, dt: 5.7e-05, norm rho: 0.02261, rho * u: 0.04592, energy: 0.08967
3891Time: 0.2, dt: 5.8e-05, norm rho: 0.02058, rho * u: 0.04361, energy: 0.08222
3892Time: 0.25, dt: 5.9e-05, norm rho: 0.01695, rho * u: 0.04203, energy: 0.06873
3893Time: 0.3, dt: 5.9e-05, norm rho: 0.01653, rho * u: 0.0401, energy: 0.06604
3894Time: 0.35, dt: 5.7e-05, norm rho: 0.01774, rho * u: 0.04264, energy: 0.0706
3898Time: 1.95, dt: 5.8e-05, norm rho: 0.01488, rho * u: 0.03923, energy: 0.05185
3899Time: 2, dt: 5.7e-05, norm rho: 0.01432, rho * u: 0.03969, energy: 0.04889
3901+-------------------------------------------+------------------+------------+------------------+
3902| Total wallclock time elapsed | 273.6s 13 | 273.6s | 273.6s 0 |
3904| Section | no. calls | min time rank | avg time | max time rank |
3905+-------------------------------------------+------------------+------------+------------------+
3906| compute errors | 41 | 0.01112s 35 | 0.0672s | 0.1337s 0 |
3907| compute transport speed | 6914 | 5.422s 35 | 15.96s | 29.99s 1 |
3908| output | 41 | 37.24s 35 | 37.3s | 37.37s 0 |
3909| rk time stepping total | 34564 | 205.4s 1 | 219.5s | 230.1s 35 |
3910| rk_stage - integrals L_h | 172820 | 153.6s 1 | 164.9s | 175.6s 27 |
3911| rk_stage - inv mass + vec upd | 172820 | 47.13s 13 | 53.09s | 64.05s 33 |
3912+-------------------------------------------+------------------+------------+------------------+
3915The norms shown here for the various quantities are the deviations
3916@f$\rho'@f$, @f$(\rho u)'@f$, and @f$E'@f$ against the background field (namely, the
3917initial condition). The distribution of run time is overall similar as in the
3918previous test case. The only slight difference is the larger proportion of
3919time spent in @f$\mathcal L_h@f$ as compared to the inverse mass matrix and vector
3920updates. This is because the geometry is deformed and the matrix-free
3921framework needs to load additional arrays for the geometry from memory that
3922are compressed in the affine mesh case.
3924Increasing the number of global refinements to 5, the output becomes:
3926Running with 40
MPI processes
3927Vectorization over 8 doubles = 512 bits (AVX512)
3928Number of degrees of freedom: 14,745,600 ( = 4 [vars] x 102,400 [cells] x 36 [dofs/cell/var] )
3932+-------------------------------------------+------------------+------------+------------------+
3933| Total wallclock time elapsed | 2693s 32 | 2693s | 2693s 23 |
3935| Section | no. calls | min time rank | avg time | max time rank |
3936+-------------------------------------------+------------------+------------+------------------+
3937| compute errors | 41 | 0.04537s 32 | 0.173s | 0.3489s 0 |
3938| compute transport speed | 13858 | 40.75s 32 | 85.99s | 149.8s 0 |
3939| output | 41 | 153.8s 32 | 153.9s | 154.1s 0 |
3940| rk time stepping total | 69284 | 2386s 0 | 2450s | 2496s 32 |
3941| rk_stage - integrals L_h | 346420 | 1365s 32 | 1574s | 1718s 19 |
3942| rk_stage - inv mass + vec upd | 346420 | 722.5s 10 | 870.7s | 1125s 32 |
3943+-------------------------------------------+------------------+------------+------------------+
3946The effect on performance is similar to the analytical test case -- in
3947theory, computation times should increase by a factor of 8, but we actually
3948see an increase by a factor of 11 for the time steps (219.5 seconds versus
39492450 seconds). This can be traced back to caches, with the small case mostly
3950fitting in caches. An interesting effect, typical of programs with a mix of
3951local communication (integrals @f$\mathcal L_h@f$) and global communication (computation of
3952transport speed) with some load imbalance, can be observed by looking at the
3953MPI ranks that encounter the minimal and maximal time of different phases,
3954respectively. Rank 0 reports the fastest throughput for the "rk time stepping
3955total" part. At the same time, it appears to be slowest for the "compute
3956transport speed" part, almost a factor of 2 slower than the
3957average and almost a factor of 4 compared to the faster rank.
3958Since the latter involves global communication, we can attribute the
3959slowness in this part to the fact that the local Runge--Kutta stages have
3960advanced more quickly on this rank and need to wait until the other processors
3961catch up. At this point, one can wonder about the reason for this imbalance:
3962The number of cells is almost the same on all
MPI processes.
3963However, the matrix-free framework is faster on affine and Cartesian
3964cells located towards the outlet of the channel, to which the lower
MPI ranks
3965are assigned. On the other hand, rank 32, which reports the highest run time
3966for the Runga--Kutta stages, owns the curved cells near the cylinder, for
3967which no
data compression is possible. To improve throughput, we could assign
3968different weights to different cell
types when partitioning the
3970few time steps and try to rebalance then.
3972The throughput per Runge--Kutta stage can be computed to 2085 MDoFs/s for the
397314.7 million DoFs test case over the 346,000 Runge--Kutta stages, slightly slower
3974than the Cartesian mesh throughput of 2360 MDoFs/s reported above.
3976Finally, if we add one additional refinement, we record the following output:
3978Running with 40
MPI processes
3979Vectorization over 8 doubles = 512 bits (AVX512)
3980Number of degrees of freedom: 58,982,400 ( = 4 [vars] x 409,600 [cells] x 36 [dofs/cell/var] )
3984Time: 1.95, dt: 1.4e-05, norm rho: 0.01488, rho * u: 0.03923, energy: 0.05183
3985Time: 2, dt: 1.4e-05, norm rho: 0.01431, rho * u: 0.03969, energy: 0.04887
3987+-------------------------------------------+------------------+------------+------------------+
3988| Total wallclock time elapsed | 2.166e+04s 26 | 2.166e+04s | 2.166e+04s 24 |
3990| Section | no. calls | min time rank | avg time | max time rank |
3991+-------------------------------------------+------------------+------------+------------------+
3992| compute errors | 41 | 0.1758s 30 | 0.672s | 1.376s 1 |
3993| compute transport speed | 27748 | 321.3s 34 | 678.8s | 1202s 1 |
3994| output | 41 | 616.3s 32 | 616.4s | 616.4s 34 |
3995| rk time stepping total | 138733 | 1.983e+04s 1 | 2.036e+04s | 2.072e+04s 34 |
3996| rk_stage - integrals L_h | 693665 | 1.052e+04s 32 | 1.248e+04s | 1.387e+04s 19 |
3997| rk_stage - inv mass + vec upd | 693665 | 6404s 10 | 7868s | 1.018e+04s 32 |
3998+-------------------------------------------+------------------+------------+------------------+
4001The "rk time stepping total" part corresponds to a throughput of 2010 MDoFs/s. The
4002overall run time to perform 139k time steps is 20k seconds (5.7 hours) or 7
4003time steps per
second -- not so bad for having nearly 60 million
4004unknowns. More throughput can be achieved by adding more cores to
4008<a name="step_67-Resultsforflowinchannelaroundcylinderin3D"></a><h3>Results for flow in channel around cylinder in 3D</h3>
4011Switching the channel test case to 3D with 3 global refinements, the output is
4013Running with 40
MPI processes
4014Vectorization over 8 doubles = 512 bits (AVX512)
4015Number of degrees of freedom: 221,184,000 ( = 5 [vars] x 204,800 [cells] x 216 [dofs/cell/var] )
4019Time: 1.95, dt: 0.00011, norm rho: 0.01131, rho * u: 0.03056, energy: 0.04091
4020Time: 2, dt: 0.00011, norm rho: 0.0119, rho * u: 0.03142, energy: 0.04425
4022+-------------------------------------------+------------------+------------+------------------+
4023| Total wallclock time elapsed | 1.734e+04s 4 | 1.734e+04s | 1.734e+04s 38 |
4025| Section | no. calls | min time rank | avg time | max time rank |
4026+-------------------------------------------+------------------+------------+------------------+
4027| compute errors | 41 | 0.6551s 34 | 3.216s | 7.281s 0 |
4028| compute transport speed | 3546 | 160s 34 | 393.2s | 776.9s 0 |
4029| output | 41 | 1350s 34 | 1353s | 1357s 0 |
4030| rk time stepping total | 17723 | 1.519e+04s 0 | 1.558e+04s | 1.582e+04s 34 |
4031| rk_stage - integrals L_h | 88615 | 1.005e+04s 32 | 1.126e+04s | 1.23e+04s 11 |
4032| rk_stage - inv mass + vec upd | 88615 | 3056s 11 | 4322s | 5759s 32 |
4033+-------------------------------------------+------------------+------------+------------------+
4036The physics are similar to the 2D case, with a slight motion in the z
4037direction due to the gravitational force. The throughput per Runge--Kutta
4038stage in this case is
4040\text{throughput} = \frac{n_\mathrm{time steps} n_\mathrm{stages}
4041n_\mathrm{dofs}}{t_\mathrm{compute}} =
4042\frac{17723 \cdot 5 \cdot 221.2\,\text{M}}{15580s} = 1258\, \text{MDoFs/s}.
4045The throughput is lower than in 2D because the computation of the @f$\mathcal L_h@f$ term
4046is more expensive. This is due to over-integration with `degree+2` points and
4047the larger fraction of face integrals (worse volume-to-surface ratio) with
4048more expensive flux computations. If we only consider the inverse mass
matrix
4049and vector update part, we record a throughput of 4857 MDoFs/s
for the 2D
case
4050of the isentropic vortex with 37.7 million unknowns, whereas the 3D
case
4051runs with 4535 MDoFs/s. The performance is similar because both cases are in
4052fact limited by the memory bandwidth.
4054If we go to four levels of global refinement, we need to increase the number
4055of processes to fit everything in memory -- the computation needs around 350
4056GB of RAM memory in
this case. Also, the time it takes to complete 35k time
4057steps becomes more tolerable by adding additional resources. We therefore use
40586 nodes with 40 cores each, resulting in a computation with 240
MPI processes:
4060Running with 240
MPI processes
4061Vectorization over 8 doubles = 512 bits (AVX512)
4062Number of degrees of freedom: 1,769,472,000 ( = 5 [vars] x 1,638,400 [cells] x 216 [dofs/cell/var] )
4066Time: 1.95, dt: 5.6e-05,
norm rho: 0.01129, rho * u: 0.0306, energy: 0.04086
4067Time: 2, dt: 5.6e-05,
norm rho: 0.01189, rho * u: 0.03145, energy: 0.04417
4069+-------------------------------------------+------------------+------------+------------------+
4070| Total wallclock time elapsed | 5.396e+04s 151 | 5.396e+04s | 5.396e+04s 0 |
4072| Section | no. calls |
min time rank |
avg time |
max time rank |
4073+-------------------------------------------+------------------+------------+------------------+
4074| compute errors | 41 | 2.632s 178 | 7.221s | 16.56s 0 |
4075| compute transport speed | 7072 | 714s 193 | 1553s | 3351s 0 |
4076| output | 41 | 8065s 176 | 8070s | 8079s 0 |
4077| rk time stepping total | 35350 | 4.25e+04s 0 | 4.43e+04s | 4.515e+04s 193 |
4078| rk_stage - integrals L_h | 176750 | 2.936e+04s 134 | 3.222e+04s | 3.67e+04s 99 |
4079| rk_stage - inv mass + vec upd | 176750 | 7004s 99 | 1.207e+04s | 1.55e+04s 132 |
4080+-------------------------------------------+------------------+------------+------------------+
4082This simulation had nearly 2 billion unknowns -- quite a large
4083computation indeed, and still only needed around 1.5 seconds per time
4087<a name=
"step_67-Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
4090The code presented here straight-forwardly extends to adaptive meshes, given
4091appropriate indicators for setting the refinement flags. Large-
scale
4092adaptivity of a similar solver in the context of the acoustic wave equation
4093has been achieved by the <a href=
"https://github.com/kronbichler/exwave">exwave
4094project</a>. However, in the present context, the benefits of adaptivity are often
4095limited to early times and effects close to the origin of sound waves, as the
4096waves eventually reflect and diffract. This leads to steep
gradients all over
4097the place, similar to turbulent flow, and a more or less globally
4100Another topic that we did not discuss in the results section is a comparison
4101of different time integration schemes. The program provides four variants of
4102low-storage Runga--Kutta integrators that each have slightly different
4103accuracy and stability behavior. Among the schemes implemented here, the
4104higher-order ones provide additional accuracy but come with slightly lower
4105efficiency in terms of step
size per stage before they violate the CFL
4106condition. An interesting extension would be to compare the low-storage
4107variants proposed here with standard Runge--Kutta integrators or to use vector
4108operations that are
run separate from the mass
matrix operation and compare
4112<a name=
"step_67-Moreadvancednumericalfluxfunctionsandskewsymmetricformulations"></a><h4>More advanced numerical flux
functions and skew-
symmetric formulations</h4>
4115As mentioned in the introduction, the modified Lax--Friedrichs flux and the
4116HLL flux employed in this program are only two variants of a large body of
4117numerical fluxes available in the literature on the Euler equations. One
4118example is the HLLC flux (Harten-Lax-van Leer-Contact) flux which adds the
4119effect of rarefaction waves missing in the HLL flux, or the Roe flux. As
4120mentioned in the introduction, the effect of numerical fluxes on high-order DG
4121schemes is debatable (unlike for the case of low-order discretizations).
4123A related improvement to increase the stability of the solver is to also
4124consider the spatial integral terms. A shortcoming in the rather naive
4125implementation used above is the fact that the energy conservation of the
4126original Euler equations (in the absence of shocks) only holds up to a
4127discretization error. If the solution is under-resolved, the discretization
4128error can give rise to an increase in the numerical energy and eventually
4129render the discretization unstable. This is because of the inexact numerical
4130integration of the terms in the Euler equations, which both contain rational
4131nonlinearities and higher-degree content from curved cells. A way out of this
4132dilemma are so-called skew-
symmetric formulations, see @cite Gassner2013 for a
4133simple variant. Skew symmetry means that switching the role of the solution
4134@f$\mathbf{
w}@f$ and test
functions @f$\mathbf{v}@f$ in the weak form produces the
4135exact
negative of the original quantity, apart from some boundary terms. In
4136the discrete setting, the challenge is to keep this skew symmetry also when
4137the integrals are only computed approximately (in the continuous case,
4138skew-symmetry is a consequence of integration by parts). Skew-
symmetric
4139numerical schemes balance spatial derivatives in the conservative form
4140@f$(\nabla \mathbf v, \mathbf{
F}(\mathbf
w))_{
K}@f$ with contributions in the
4141convective form @f$(\mathbf v, \tilde{\mathbf{
F}}(\mathbf
w)\nabla
4142\mathbf{
w})_{
K}@f$ for some @f$\tilde{\mathbf{
F}}@f$. The precise terms depend on
4143the equation and the integration formula, and can in some cases by understood
4144by special skew-
symmetric finite difference schemes.
4146To get started, interested readers could take a look at
4148skew-
symmetric DG formulation is implemented with deal.II for a simple advection
4151<a name=
"step_67-Equippingthecodeforsupersoniccalculations"></a><h4>Equipping the code for supersonic calculations</h4>
4154As mentioned in the introduction, the solution to the Euler equations develops
4155shocks as the Mach number increases, which require additional mechanisms to
4156stabilize the scheme,
e.g. in the form of limiters. The main challenge besides
4157actually implementing the limiter or artificial viscosity approach would be to
4158load-balance the computations, as the additional computations involved for
4159limiting the oscillations in troubled cells would make them more expensive than the
4160plain DG cells without limiting. Furthermore, additional numerical fluxes that
4161better cope with the discontinuities would also be an option.
4163One ingredient also necessary for supersonic flows are appropriate boundary
4164conditions. As opposed to the subsonic outflow boundaries discussed in the
4165introduction and implemented in the program, all characteristics are outgoing
4166for supersonic outflow boundaries, so we do not want to prescribe any external
4169\mathbf{
w}^+ = \mathbf{
w}^- = \
begin{pmatrix} \rho^-\\
4170(\rho \mathbf u)^- \\
E^-\
end{pmatrix} \quad
4174In the code, we would simply add the additional statement
4176 else if (supersonic_outflow_boundaries.find(
boundary_id) !=
4177 supersonic_outflow_boundaries.
end())
4183in the `local_apply_boundary_face()` function.
4185<a name=
"step_67-ExtensiontothelinearizedEulerequations"></a><h4>Extension to the linearized Euler equations</h4>
4188When the interest with an Euler solution is mostly in the propagation of sound
4189waves, it often makes sense to linearize the Euler equations around a
4190background state, i.e., a given density, velocity and energy (or pressure)
4191field, and only compute the change against these fields. This is the setting
4192of the wide field of aeroacoustics. Even though the resolution requirements
4193are sometimes considerably reduced, implementation gets somewhat more
4194complicated as the linearization gives rise to additional terms. From a code
4195perspective, in the
operator evaluation we also need to equip the code with
4196the state to linearize against. This information can be provided either by
4197analytical
functions (that are evaluated in terms of the position of the
4198quadrature points) or by a vector similar to the solution. Based on that
4199vector, we would create an additional
FEEvaluation object to read from it and
4200provide the
values of the field at quadrature points. If the background
4201velocity is zero and the density is
constant, the linearized Euler equations
4202further simplify and can equivalently be written in the form of the
4203acoustic wave equation.
4205A challenge in the context of sound propagation is often the definition of
4206boundary conditions, as the computational domain needs to be of finite
size,
4207whereas the actual simulation often spans an infinite (or at least much
4208larger) physical domain. Conventional Dirichlet or Neumann boundary conditions
4209give rise to reflections of the sound waves that eventually propagate back to
4210the region of interest and spoil the solution. Therefore, various variants of
4211non-reflecting boundary conditions or sponge layers, often in the form of
4213href=
"https://en.wikipedia.org/wiki/Perfectly_matched_layer">perfectly
4214matched layers</a> -- where the solution is damped without reflection
4218<a name=
"step_67-ExtensiontothecompressibleNavierStokesequations"></a><h4>Extension to the compressible Navier-Stokes equations</h4>
4221The solver presented in
this tutorial program can also be extended to the
4222compressible Navier--Stokes equations by adding viscous terms, as described in
4223@cite FehnWallKronbichler2019. To keep as much of the performance obtained
4224here despite the additional cost of elliptic terms,
e.g. via an interior
4226the @ref step_59
"step-59" tutorial program.
4229<a name=
"step_67-Usingcellcentricloopsandsharedmemory"></a><h4>Using cell-centric loops and shared memory</h4>
4232In
this tutorial, we used face-centric loops. Here, cell and face integrals
4233are treated in separate loops, resulting in multiple writing accesses into the
4234result vector, which is relatively expensive on modern hardware since writing
4235operations generally result also in an implicit read operation. Element-centric
4236loops, on the other hand, are processing a cell and in direct succession
4237processing all its 2
d faces. Although
this kind of
loop implies that fluxes have
4238to be computed twice (
for each side of an interior face), the fact that the
4239result vector has to accessed only once might - and the fact that the resulting
4240algorithm is free of race-conditions and as such perfectly suitable
for
4241shared memory - already give a performance
boost. If you are interested in these
4242advanced topics, you can take a look at @ref step_76
"step-76" where we take the present
4243tutorial and modify it so that we can use these features.
4246<a name=
"step_67-PlainProg"></a>
4247<h1> The plain program</h1>
4248@include
"step-67.cc"
void write_vtu_in_parallel(const std::string &filename, const MPI_Comm comm) const
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={})
virtual UpdateFlags get_needed_update_flags() const =0
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
virtual std::vector< std::string > get_names() const =0
virtual std::vector< DataComponentInterpretation::DataComponentInterpretation > get_data_component_interpretation() const
void submit_dof_value(const value_type val_in, const unsigned int dof)
void set_dof_values(VectorType &dst, const unsigned int first_index=0, const std::bitset< n_lanes > &mask=std::bitset< n_lanes >().flip()) const
value_type get_value(const unsigned int q_point) const
const ShapeInfoType * data
Tensor< 2, dim, VectorizedArrayType > inverse_jacobian(const unsigned int q_point) const
std_cxx20::ranges::iota_view< unsigned int, unsigned int > quadrature_point_indices() const
Point< dim, VectorizedArrayType > quadrature_point(const unsigned int q) const
const unsigned int n_q_points
void gather_evaluate(const VectorType &input_vector, const EvaluationFlags::EvaluationFlags evaluation_flag)
void evaluate(const EvaluationFlags::EvaluationFlags evaluation_flag)
virtual RangeNumberType value(const Point< dim > &p, const unsigned int component=0) const
void reinit(const size_type size, const bool omit_zeroing_entries=false)
Abstract base class for mapping classes.
void transform_from_q_points_to_basis(const unsigned int n_actual_components, const VectorizedArrayType *in_array, VectorizedArrayType *out_array) const
unsigned int n_active_entries_per_cell_batch(const unsigned int cell_batch_index) const
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 > &)> &inner_face_operation, const std::function< void(const MatrixFree< dim, Number, VectorizedArrayType > &, OutVector &, const InVector &, const std::pair< unsigned int, unsigned int > &)> &boundary_face_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
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
void reinit(const MappingType &mapping, const DoFHandler< dim > &dof_handler, const AffineConstraints< number2 > &constraint, const QuadratureType &quad, const AdditionalData &additional_data=AdditionalData())
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
void print_wall_time_statistics(const MPI_Comm mpi_comm, const double print_quantile=0.) const
static constexpr std::size_t size()
#define DEAL_II_ALWAYS_INLINE
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertThrow(cond, exc)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_NOT_IMPLEMENTED()
std::vector< index_type > data
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
DataComponentInterpretation
@ component_is_part_of_vector
void approximate(const SynchronousIterators< std::tuple< typename DoFHandler< dim, spacedim >::active_cell_iterator, Vector< float >::iterator > > &cell, const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof_handler, const InputVector &solution, const unsigned int component)
The namespace for the EvaluationFlags enum.
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
void cylinder(Triangulation< dim > &tria, const double radius=1., const double half_length=1.)
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)
@ valid
Iterator points to a valid object.
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
@ general
No special properties.
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > E(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
@ LOW_STORAGE_RK_STAGE9_ORDER5
@ LOW_STORAGE_RK_STAGE3_ORDER3
@ LOW_STORAGE_RK_STAGE7_ORDER4
@ LOW_STORAGE_RK_STAGE5_ORDER4
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm mpi_communicator)
T max(const T &t, const MPI_Comm mpi_communicator)
T min(const T &t, const MPI_Comm mpi_communicator)
unsigned int this_mpi_process(const MPI_Comm mpi_communicator)
std::string get_current_vectorization_level()
Number truncate_to_n_digits(const Number number, const unsigned int n_digits)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
constexpr T pow(const T base, const int iexp)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
long double gamma(const unsigned int n)
DEAL_II_HOST constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
unsigned int n_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
int(&) functions(const void *v1, const void *v2)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
static constexpr double PI
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
****code * * MPI_Finalize()
bool write_higher_order_cells
UpdateFlags mapping_update_flags
DEAL_II_HOST constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
1.9.8