Results: Shock Detection

With the ANNSI networks for $`\mathcal{N}=5`$ and $`\mathcal{N}=9`$ from Sec. 13.4 in place, we now apply the trained shock indicators to a range of 2D shock problems. We first describe the selected flow cases briefly in Sec. [subsec:testcases]. With these definitions out of the way, we then discuss the results for the ANNSI approach and contrast it against the established indicators from Sec. 14.4. We first focus here on the shock detection aspect part of the proposed indicator, i.e. the marking of grid elements where shocks likely reside. The aspect of shock localization within an element will be discussed in Sec. [sec:NNshock_local]. For all computations presented herein, we rely on the discretization scheme, including the TVD shock capturing method, presented in Sec. 14.2 and 14.3. The polynomial ansatz degree of the scheme is either $`\mathcal{N}=5`$ and $`\mathcal{N}=9`$ as indicated; the scheme is a standard DGSEM collocation formulation on Legendre-Gauss points, with a Roe or HLLE scheme for the numerical flux functions of both DG and FV elements. A minmod limiter ensures the TVD property of the linear FV reconstruction polynomials. The solution is advanced in time by an explicit high order Runge-Kutta scheme ; the typical CFL condition limits the allowable global timestep . The density is chosen as an input quantity into ANNSI.
All simulations with ANNSI or ANNSL shown in this manuscript are run ab initio, i.e. from the given initial conditions at time $`t_{start}=0.0`$ to the indicated $`t_{end}`$, using the proposed indicators only. While we focus on discussing the solution quality at a given time typical for the selected test case, this shows that ANNSI and ANNSL approaches are capable of dealing with highly transient situations without the need for any special treatment.

Test Cases

2D Riemann Problems

In , the definition of the 1D Riemann (RP) problem for the Euler system is extended to two dimensions. The authors identify 16 distinct configurations in which constant states in the four domain quadrants are separated initially by an elementary wave and propose them as canonical test cases. We choose three configurations (4,6,12), which among them include the three typical wave phenomena, shocks, rarefactions and contact discontinuities. The initial conditions and end times of the RPs are given in Tbl. 1; they define the solution in the four quadrants of the domain (numbered from the upper right quadrant counterclockwise). The boundary conditions on all domain faces are set to the exact solutions.

Fig. 1 shows the resulting density contours for the three configurations as a reference, computed on a fine grid with the numerical scheme described above. Note that along the contact discontinuities in configuration 12, a Kelvin-Helmholtz instability triggers the roll-up of vortical structures. For the comparisons of the shock indicators later on, we define a baseline regular Cartesian mesh $`M_{RP}`$ with $`50\times 50`$ elements for all configurations.

Double Mach Reflection Problem

This testcase consists of a strong shock at Mach 10 impinging on a wedge at an oblique angle of $`30^{\circ}`$. Our setup is consistent with the description in , with slip-wall boundary conditions at the top and bottom, a Dirichlet inflow on the left and a supersonic outflow at the right hand side of the domain. Fig. 2 shows the reference solution at $`t_\text{end}=0.2`$ computed on a fine grid. We define a coarse baseline grid $`M_{DMR}`$ with $`49\times 12`$ elements for later use in this study. The numerical flux is approximated by a Roe scheme.

Forward Facing Step Problem

This test setup, introduced in  describes the flow over a forward facing step (FFS) with a Mach number of $`\mathrm{Ma}=3.0`$. The expected flow phenomena are mainly a developing bow shock in front of the step and the reflections of this shock from the upper and lower wall. Inflow and outflow are described with Dirichlet boundary conditions, other boundaries are modeled as walls. In Fig. 3 the density at $`t_\text{end}=4`$ is visualized for a calculation with $`50`$ elements in the $`x`$-direction and $`25`$ and $`20`$ elements in the two blocks in $`y`$-direction using $`\mathcal{N}=9`$ and an HLLE Riemann solver. For more details on the used mesh see also Fig. 14.

2D NACA 0012

As a final testcase, we consider the flow over a NACA 0012 airfoil at $`Ma=2.0, AoA=0^{\circ}`$ on an unstructured grid, depicted in Fig. 4. This setup results in a detached, curved bow shock and a tail shock. The solution is advanced in time from an initial freestream condition until a steady state is achieved. The mesh is coarse, without symmetries and does not take a priori knowledge of the shock front position, strength and shape into account, resulting in a deliberately “bad” mesh for this task. This setup was generated to test the shock indicators in a more realistic setting for high order methods than in the previous test cases, where the shock location might be uncertain or rapidly moving, grid refinement is not available or too costly, and large grid elements are the norm. Fig. 4 shows the computational grid and the flow solution.

ANNSI Results for $`\mathcal{N}=5`$

We first consider the results for the $`\mathcal{N}=5`$ computations and networks, and discuss selected testcases from Sec. [subsec:testcases] and contrast the results from ANNSI and the troubled cell indicators from Sec. 14.4. It is important to underline that while for the modal and jump indicators $`\mathcal{I}_\text{modal}`$ and $`\mathcal{I}_\text{jump}`$ the respective parameters defining the thresholds need to be set by the user and possibly adapted for each situation, the ANNSI prediction of the neural network established in Sec. 13 is used here ’as is’ throughout this chapter, without the need (or possibility) of user intervention or parameter adjustment.
In a first step, we apply the ANNSI to a range of test problems to establish its suitability, robustness and accuracy for the intended purpose.

Fig. 5 shows the resulting flow fields in the top row and the flagged grid elements in the bottom row for the Riemann problems and the DMR. All simulations using ANNSI are stable and the resulting flow fields are qualitatively as expected. The flagged cells also coincide with the shock locations, the shock fronts are unbroken and well-captured. For the Riemann problem configuration 12, the contact discontinuity lines are not detected, which is likely due to the roll up of the shear layer and the weak gradients.

Based on these initial findings that the ANNSI method delivers a useful and robust shock detection, we now compare its predictions against the $`\mathcal{I}_\text{modal}`$ and $`\mathcal{I}_\text{jump}`$ indicators. As pointed out above, the upper and lower detection thresholds of these indicators are tuneable parameters and no strict quantitative guidelines exist that can help fix them. Instead, the task of the user is to find a setting mainly through trial and error that results in sufficient robustness but do not trigger excessive solution stabilization mechanisms, which can either be computationally expensive or unduly influence smooth regions of the solution - or both. Beyond the academic examples considered here, such parameter variations however are often too costly and cumbersome, and thus often a parameter set based on previous experience is chosen. In order to now compare this approach with the one proposed here, we hand-tuned the $`\mathcal{I}_\text{modal}`$ and $`\mathcal{I}_\text{jump}`$ indicators collectively on the four test cases from Fig. 5 (Riemann problems 4,6,12 and DMR) under the following conditions: For each of the two indicators, the optimal parameter settings must provide a stable simulation of all four test cases and must not result in excessive flagging of cells outside of the expected regions. Our approach resulted in the settings listed in Tbl. 2.

We compare these results with the ANNSI cases in Fig. 6, where we have repeated the relevant results from Fig. 5 as the first column. All three indicators yield comparable results, flagging shock locations but not expansion waves. The ANNSI predictions typically flag two adjacent cells across a shock, while in particular the jump indicator shows a more refined picture. On the other hand, the shock fronts detected by ANNSI are continuous and unbroken in contrast to the $`\mathcal{I}_\text{modal}`$ and $`\mathcal{I}_\text{jump}`$ results, where “holes” in the shock fronts appear. Applying the same indicators to the NACA0012 case leads to a stable simulation in all cases; the results for the flagged cells are shown in Fig. 7. The ANNSI driven method performs well on this suboptimal grid, yielding continuous shock fronts and a more conservative estimate than the other two. We note that as the computational grid is not symmetric w.r.t the $`x-`$axis, so we also do not expect a necessarily symmetric flagging of the shock fronts.

The results presented in this section can thus be summarized as: The developed ANNSI method is robust, suitable for structured and unstructured grids and of comparable accuracy to other well-tuned indicators.

ANNSI Results for $`\mathcal{N}=9`$

In this section we apply the ANNSI network to the same 2D problems as in the previous section. We perform the calculations with a polynomial degree of $`\mathcal{N}=9`$ to show the invariance of the performance with respect to the polynomial degree of the ansatz functions.

In Fig. 8 we repeat the calculations from Fig. 5 with $`\mathcal{N}=9`$ instead of $`\mathcal{N}=5`$. A very similar behavior can be observed: again all simulations using the ANNSI are stable and the resulting flow fields are qualitatively as expected. The shock localization is correctly recognized and captured with the finite volume sub-cell scheme. Differences can be seen at the contact discontinuities in RP12 which are partially detected by the indicator. We also compared the ANNSI results against those for the selected problems in a repetition of Fig. 6. We found the results from $`\mathcal{N}=5`$ transferable to the $`\mathcal{N}=9`$ case.

Influence of Grid Resolution

Having established the properties of the ANNSI approach on the baseline grids, we now investigate the influence of grid resolution changes by either halving or doubling the grid cells for the $`\mathcal{N}=5`$ case per direction. Since the training data for ANNSI was constructed in reference space, we expect its results to be robust against grid resolution changes. For $`\mathcal{I}_\text{jump}`$ and $`\mathcal{I}_\text{modal}`$, we again use the found parameter settings from Sec. [subsec:n5].

Fig. 9 and 10 repeat the investigations for the coarse and fine grids. Missing plots indicate that the simulation was not stable. The coarse grid results confirm our assumption that the ANNSI approach is rather robust against a resolution drop. Comparing Fig. 9 and Fig. 6, the shock position and extension are predicted in a consistent manner. On the fine grid in Fig. 10, the ANNSI results converge and become sharper as desirable.

Summarizing the findings in this section, the developed ANNSI approach results in robust solutions for the problems considered. It compares well against other solution-based indicators in terms of accuracy and does not excessively flag elements outside of shock regions and can deal robustly with the initial solution transients and shock movement. It has shown to be insensitive to grid resolution and consistently flags the same regions / features, thus working as intended as a discretization independent shock detector, and not as a troubled cell indicator. While results of comparable or possibly better quality for specific cases can be achieved by parameter studies and tuning with other a priori indicators, once trained, ANNSI is parameter-free. It can be argued that the time spent training the neural network should be weighted against the time one might invest in parameter tuning, but with large scale application with huge computational costs in mind, it is certainly more economic to do invest the costs “offline” before the actual computation. This aspect will become even more important in the next section, where we discuss the results of the ANNSL extension of ANNSI, where a localization of the shock front within an element is proposed.

Results: Shock Localization

In this section, we illustrate how the information about the shock location inside an element generated by the ANNSL network can be exploited. First, the test cases from the previous section are used to illustrate the capabilities of the shock localization procedure. Next, the ANNSL indicator is applied to a new testcase and the possible exploitation of the shock localization in an adaptive mesh refinement (AMR) framework is illustrated. All ANNSL computations are again started from ab initio, with ANNSL being active in every timestep.

Shock Localization

In the previous section, we have shown that the ANNSI shock detector performs well in terms of robustness, accuracy and insensitivity to spatial resolution changes for a range of classical test problems. Here, we use those applications again to illustrate the capabilities of the ANNSL network of identifying the shock position inside an element, i.e. of creating an inner-element bounding box around the shock front and flagging the solution points adjacent to the front. Since this is particularly useful for high order discretizations such as $`\mathcal{N}=9`$, we will focus our investigation on this case, using the network trained in Sec. 13.4 on the data from Sec. 13.3.2.

Fig. 11 and Fig. 12 apply the ANNSL network to the established test cases. The binary edge map is shown alongside the flow solution, with dark colors indicating that a pixel / solution point is considered to be in the direct vicinity of or on the shock front. One can observe that with the proposed indicator, a very sharp identification of discontinuities is possible. The identification is generally continuous within a cell (no holes or bumps in the shock fronts) and consistent across element boundaries - neighboring element “agree” in the identification of the fronts. All of these properties are desirable in edge detection, however, many more simple edge detection algorithms struggle with one or more of them . For the intentionally “bad” grid in Fig. 12, the ANNSL prediction also works very well on coarse cells. The infrequently occurring ’two-stripes’ patterns (see e.g. lower region of the bow shock) are due to the discrete resolution of a shock: on the discrete level, it consists of two regions with kinks at the foot and head of the shock and an almost linear part in between. As the linear part is not distinguishable from other smooth solutions, it is not recognized as a shock, and only the shock edges are recognized as parts of a shock.
In Fig. 13, we analyze the convergence behaviour of the shock map with decreasing element size. We observe that the shock localization is stable with respect to the grid size, and the edge map predictions become successively more refined. This feature is an important prerequisite for successful h/p/r-adaption strategies to improve accuracy at the shock. In the next section, we discuss a possibility for taking advantage of this knowledge during the simulation.

Neural Network informed H-Adaptation

We use the flow over a forward facing step problem on a regular Cartesian mesh defined in Sec. [subsec:ffs] to illustrate how the information about the inner-cell shock position can be exploited during the calculation to improve the accuracy of the shock capturing approach. We restrict the discussion to a three-level h-refinement strategy in which a donor cell on the original grid can be halved in any coordinate direction separately or in both at the same time, leading to either 2 or 4 new child cells replacing the original element (anisotropic or isotropic split). This process can be repeated twice, so three levels of grid elements can be created. The original regular grid is shown in the top left column in Fig. 14. We conduct all simulations with a polynomial degree of $`\mathcal{N}=9`$ and the ANNSL network from the previous sections and visualize the results at $`t_\text{end}=4`$.

The left column of Fig. 14 shows the solution on the baseline mesh, which is direction-wise equidistant. The corresponding shock localization indicator and the resulting density field are shown below, again, the ANNSL approach leads to a sharp and consistent prediction of the front locations, confirming the findings from Sec. [subsec:annsl1]. From the binary edge map, two pieces of information can now be extracted: First, the shock width is characterized by the amount of flagged solution points / pixels within one element. Secondly, the orientation of the shock front is available. Several methods to take advantage of this information to inform an h-refinement are possible; we choose a very simple strategy for demonstration purposes here which already produces remarkable results: We define an indicator for anisotropic mesh refinement which exploits both information from the binary edge map

\begin{align}
\label{eq:MeshRefIndicator}
    \mathcal{I}_\text{meshref}^{\text{dir}}=\sum_{i=0}^{\mathcal{N}}\sum_{s=1}^{r^{\text{dir}}_i-1}\left(1.7^{s-1}+1.7^s\right),
\end{align}

independently for both the $`x`$ and $`y`$-directions, with $`r^{\text{dir}}_i`$ being the number of solution points with class $`1`$ of layer $`i`$ in the corresponding $`x`$ or $`y`$-direction, i.e. we count the number of flagged pixels line-wise. After evaluating the indicator we apply the mesh refinement in each direction independently. Last, a constant interpolation on the new grid is done. We apply the indicator again on the refined solution, choosing the threshold values $`\mathcal{I}^\text{threshold}_\text{meshref}=33`$ and $`\mathcal{I}^\text{threshold}_\text{meshref}=172`$ for the first and the second refinement. These three levels of grid elements provide a sharp resolution of the shocks. Note that during the refinement it is taken care that only non-conforming interfaces with a 2:1 ratio occur. Details on DGSEM and the finite volume sub-cell method on meshes with non-conforming interfaces can be found in . Fig. 14 (right column) shows the resulting mesh, the results on the refined mesh and the binary edge map, illustrating the sharper resolution of the shocks due to the refined mesh.

Fig. 15 (top row) provides a zoomed-in view on the upper crossing point of the shocks, showing both the density field and the associated mesh. It is easy to recognize how the refinement informed by ANNSL has improved the sharpness of the shock while avoiding unnecessary refinement of cells or along directions. Also, the refinement criterion proposed in Eq. [eq:MeshRefIndicator] has guided the selection between anisotropic and isotropic splits, which additionally helps to keep the number of newly introduced cells and thus computational cost low. In the bottom row of Fig. 15, the computational grid and the binary edge map in this region are depicted. Cells of the computational grid are denoted by black lines; in elements with white backgrounds, no shocks have been detected and the DG representation of the solution is used. In cells with a light blue background, shocks are present, and the FV representation of the solution on equispaced grid cells is active. In FV cells / pixels of the binary edge map marked with a dark color, the ANNSL method has predicted the occurrence of a shock. As discussed before, the network has been designed to mark the two adjacent pixels to the shock front per direction, i.e. to provide a bounding box of the discontinuity on the scale of smallest resolved length scale within a DG element. The bottom plots in Fig. 15 shows that ANNSL method successfully provides such a bounding box along the shock front, given a consistent and sharp inner-cell shock localization on both the original (left) and refined grid (right).
Concluding, this example illustrates how the information obtained by the novel indicator for shock localization can be used in an anisotropic adaptive mesh refinement framework to improve the resolution of the shock. Due to the separation of shock detection and shock capturing in our approach, the indicator can be transferred to other discretization schemes and can then help to inform other shock capturing methods.

Generation of 1D Training Data

In the following, we detail the generation of training samples for the different analytical functions with $`\eps=0.1`$, chosen from 7 families:

• No. 1: Linear functions, initialized on three meshes $`n_e=1,10,20`$. The probability of choosing each of the meshes is $`p_e=0.5,0.3,0.2`$, respectively.

• No. 2: Superposition of oscillations with different amplitudes and frequencies, initialized on meshes with $`n_e=1,10,20`$ and probabilities of $`p_e=0.3,0.4,0.3`$. Note that $`f_\text{Nyquist}=\mathcal{N}/2`$ approximates the maximum possible frequency which can be resolved by the Lagrange polynomials.

• No. 3: Exponential functions in both directions, initialized on $`n_e=10,20`$ with probabilities $`p_e=0.6,0.4`$. No function evaluation with $`n_e=1`$ is done in order not to indicate too steep gradients with class $`0`$.

• No. 4: Four constant values separated with straight and curved lines, initialized on $`n_e=1`$. The assignment to class $`1`$ is done if the maximum jump height fulfills the condition $`|\text{max}(|u_i|)-\text{min}(u_i)|>\eps~\text{max}(|u_i|)`$ and $`\text{max}(|u_i|)>0.01`$ holds. Degrees of freedom neighboring or being directly at the discontinuity are labeled as $`1`$. The parameter $`d`$ is chosen to give a ratio of straight to curved separation lines of 7:3.

• No. 5: Magnitude function with linear gradients connected by kinks, initialized on $`n_e=1`$. The assignment to class $`1`$ is done if the condition $`a>\eps~\text{max}(|u_i|)`$ for the gradient $`a`$ and $`\text{max}(|u_i|)>0.01`$ for the maximum value holds and is set for degrees of freedom neighboring or being directly at the discontinuity.

• No. 6: Linear or constant states connected by kinks, initialized on $`n_e=1`$. The assignment to class $`1`$ is done if the maximum gradient fulfills $`\text{max}(b_1,b_2)>\eps~\text{max}(|u_i|)`$ and $`\text{max}(|u_i|)>0.01`$ holds and is set for degrees of freedom neighboring or being directly at the discontinuity.

• No. 7: Decaying high frequency oscillations in analogy to Gibb’s instability, initialized on $`n_e=1`$. The assignment to class $`1`$ is done if the maximum difference in the element fulfills the condition $`|\text{max}(|u_i|)-\text{min}(u_i)|>0.2~\text{max}(|u_i|)`$ and $`\text{max}(|u_i|)>0.01`$ holds. Then degrees of freedom having the highest absolute value in the element are labeled as $`1`$.

Tbl. 3 summarizes the function families and their parameters. After the generation of the training samples, the inputs are normalized to enhance the convergence of the training process. We first shift the data

\begin{align*}
u_i = \begin{cases}
u_i - \text{min}_x u_i & \text{if} \ \text{min}_x u_i < 0, \\
u_i                    & \text{otherwise}, \end{cases}
\end{align*}

and then scale them to the interval $`\mathcal{I}=[0,1]`$

\begin{align*}
u_i = \begin{cases}
\frac{u_i}{\text{max}_x \vert{u_i\vert}} & \text{if} \ \text{max}_x \vert{u_i\vert} > 1, \\
u_i                    & \text{otherwise}. \end{cases}
\end{align*}

This enables to train with data in the positive interval $`\mathcal{I}=[0,1]`$ while retaining the information about relative differences in the solution. To obtain training samples from the analytical functions, they are initialized on Cartesian meshes covering the interval $`[-1,1]^2`$ with $`n_e^2`$ elements. To account for under resolved phenomena, we use a polynomial degree in each grid cell of $`2\mathcal{N}`$ for the evaluation of the analytical functions. Next, the analytical functions given in the Lagrange basis representation of degree $`2\mathcal{N}`$ are projected onto the Lagrange basis representation of degree $`\mathcal{N}`$, which is used for the simulation with DGSEM. Summing up, the input image $`\vec{X}`$ is build up with the normalized nodal values of the Lagrange polynomials of degree $`\mathcal{N}`$ resulting in $`(\mathcal{N}+1)^2`$ values. The true output $`\vec{Y}`$ is the binary edge map.