An adaptive mesh refinement approach based on optimal sparse sensing

Abstract

We introduce a new approach for adaptive mesh refinement in which adaptivity is driven by low rank decomposition and optimal sensing of the dynamically evolving flow field. This method seeks an ordered set of locations for mesh adaptation from the instantaneous data-driven basis of an online proper orthogonal decomposition of the velocity, which organizes features into sparse optimal orthogonal modes based on an energy norm. The sensing is achieved via a computationally expedient discrete empirical interpolation method using rank-revealing QR factorization (Drmac and Gugercin SIAM J Sci Comput 38(2):A631–A648, 2016). The methodology is applicable to a wide range of numerical discretizations, and is tested on a spatiotemporally evolving incompressible turbulent jet, a complex wind turbine wake, and supersonic flow over a forward-facing step. The proposed approach is demonstrated to predict accurate velocity statistics and yield significantly smaller grids in comparison to gradient-based methods. The algorithm is seen to focus refinement in the vicinity of dynamically significant regions such as those characterized by high turbulence kinetic energy, coherent structures and shock interactions. Moreover, the approach does not require parameters or thresholds, which may be difficult to obtain for complex flows, to be known a priori to facilitate mesh adaptation.

Appendices

Governing equations and setup for spatially evolving turbulent jet

The mean streamwise velocity at the inflow is based on a hyperbolic tangent profile [17] and given by

$$\begin{aligned} u_1(x_1=0, x_2, x_3) = \frac{U_1 + U_2}{2} + \frac{U_1 - U_2}{2} \tanh \left( \frac{h}{4\theta _0} \left( 1 - \frac{2 |x_2|}{h} \right) \right) , \end{aligned}$$

(23)

where $U_1$ is the centerline jet velocity, $U_2$ is the co-flow, $\theta _0$ is momentum thickness, and h is the slot width of the inlet. Both transverse mean velocities $u_2$ and $u_3$ are set to zero at the inlet. The top-hat inflow profile is obtained by mirroring the profile around the $x_2$ plane. A small fluctuating velocity is added to each velocity component of the mean inflow profile to include statistics characteristic of isotropic turbulence. The fluctuations are prescribed to impose an energy spectrum given by

$$\begin{aligned} E(k) = k^s \exp \left( -\frac{s}{2} \left( \frac{k}{k_0} \right) ^2 \right) \end{aligned}$$

(24)

where the wavenumber $k = (k_2^2 + k_3^2)^{1/2}$ and s and peak wavenumber $k_0$ are the selected such that the energy input is dominated at small scales or high $k_0$ and typical of decaying isotropic turbulence $s \le 4$. The fluctuations are only imposed near the shear layers of the inflow.

The modified vorticity transport equations in compact index notation are as follows ($i,j=1,2,3$):

$$\begin{aligned} \frac{\partial \omega _i}{\partial t} + \frac{\partial }{\partial x_j} \left( u_j \omega _i - u_i \omega _j \right) + u_i \frac{\partial \omega _j}{\partial x_j} = \frac{1}{\mathrm{Re}} \frac{\partial ^2 \omega _i}{\partial x_j \partial _j} + \epsilon _{ijk} \frac{\partial F_k}{\partial x_j}, \end{aligned}$$

(25)

where $u_i$ is the velocity, $\omega _i = \epsilon _{ijk} \partial u_k / \partial x_j$ is the vorticity, and $F_k$ is a body force. $\mathrm{Re} = U L/\nu$ is the Reynolds number defined by a characteristic velocity scale U, characteristic length L, and kinematic viscosity $\nu$.

The velocity is obtained from the vorticity through a Poisson equation given by

$$\begin{aligned} \frac{\partial ^2 u_i}{\partial x_j \partial x_j} = - \epsilon _{ijk} \frac{\partial \omega _k}{\partial x_j}. \end{aligned}$$

(26)

The governing equations are discretized with a multi-dimensional upwind finite volume approach as described in Ref. [24]. This approach uses upwind differences corrected by second-order differences around smooth flow regions to produce a solution that is second-order accurate in space and time. The multi-dimensional scheme solved in a series of dimensional sweeps shown as follows (for clarity vectors are shown in bold):

$$\begin{aligned} \varvec{\omega }^{n+1}_{i,j,k} = \varvec{\omega }^n_{i,j,k} - \frac{\Delta t}{\Delta _x} \sum ^3_{m=1} \left[ \left( (\varvec{F}^l_m)^n_{i+\frac{1}{2},j,k} - (\varvec{F}^r_m)^n_{i-\frac{1}{2},j,k} \right) - \left( (\varvec{G}_m)^n_{i+\frac{1}{2},j,k} - (\varvec{G}_m)^n_{i-\frac{1}{2},j,k} \right) \right] , \end{aligned}$$

(27)

where n is the time iteration and i, j, k are the indices of the grid cells in three directions. The flux functions $\varvec{F}^l, \varvec{F}^r$ are normal direction flux function are obtained by solving a generalized Riemann problem developed for the VTE [50]. The transverse direction flux functions $\varvec{G}$ are fluxes from the transverse direction and are computed using the flux-based wave propagation approach [24]. This method has been shown to accurately predict a number of vortex-dominated and turbulent flows [23,24,25].

The domain is carefully considered to accurately simulate a spatially evolving jet in the $x_1$ direction with the VTE. The size of the computational domain in each direction $(L_1 \times L_2 \times L_3) = (12h \times 12h \times 3h)$, where inflow and outflow conditions are specified in the $x_1$ direction and periodic in the other two directions. While the inflow conditions for the velocity are prescribed with Eq. (23), the outflow boundary conditions must be considered separately for vorticity and velocity. Due to the hyperbolicity of Eq. (25), the vorticity at the outflow that convects out of the domain due to a positive velocity at the outlet does not have a significant effect on the vorticity inside the domain; however, the velocity is strongly dependent on the vorticity due to the Biot–Savart law as computed through Eq. (26). To ensure that the boundary conditions do not affect the solution, the domain is extended with a buffer region in the $x_1$ direction to 36h. The outflow boundary condition at the boundary of the buffer region for velocity is computed using a fast multi-pole method. The domain with buffer region is shown in Fig. 19.

The domain is initially discretized with several refinement levels of uniform structured grid cells, where the coarsest level uses $(N_1 \times N_2 \times N_3) = (96 \times 32 \times 8)$ across a $(36h \times 12h \times 3h)$ region including a computation and buffer region. This corresponds to a coarsest grid resolution of $\Delta x_i = 0.375h$. A second level of refinement (with a refinement factor of 2 in each direction) is applied in the computation domain $(12h \times 12h \times 3h)$ with $(N_1 \times N_2 \times N_3) = (64 \times 64 \times 16)$, which ensures a uniform grid size of $\Delta x_i = 0.1875h$. Additional levels of refinement are only applied on the domain with a refinement factor of 2.

Fig. 19

The domain outlined in blue and entire computational domain with buffer region outlined red. The buffer region is discretized with the coarsest grid cells (color figure online)

The momentum thickness of the jet at the inlet is set to be $h/\theta _0 = 30$, and a slot width $h = 10$. The centerline velocity is $U_1 = 1.091$ and the co-flow velocity is $U_2 = 0.091$. The Reynolds number $\text {Re}_h = \Delta U h / \nu = 3000$, where $\Delta U = U_1 - U_2$. The magnitude of the inlet velocity fluctuations are given at 10%.

Validation of spatially evolving turbulent jet

The validation of the static grid cases, S1, S2, and S3, for the planar jet using statistics of the flow quantities is detailed. The temporal averaging of the velocity and vorticity flow fields begins after each simulation reaches a statistically steady state, and continued until the averaged flow field reaches statistical convergence. This corresponds to a total elapsed time of temporal averaging over 200 Kelvin–Helmholtz (K-H) vortices to leave the domain if we consider the frequency of K-H vortices to be $St_{sl} = f_{sl} \theta _0/U_v = 0.0333$ where $U_v = 1/2(U_1 + U_2)$.

The mean velocity flow field, $\langle u_1 \rangle$ (also ensemble-averaged over the spanwise $x_3$ direction), where $\langle \cdot \rangle$ indicates $x_3$ direction averaging, is normalized by the difference in the streamwise centerline velocity $U_c = \langle u_1 \rangle (x_2=0)$ and the local streamwise co-flow velocity $U_\infty = \langle u_1 \rangle (x_2=6h)$. The jet half-width $\delta _{1/2}$, which is calculated using the mean streamwise velocity $\langle u_1 \rangle (x_2=\delta _{1/2}) - U_\infty = \frac{1}{2} \left( U_c - U_\infty \right)$ is used to normalize the streamwise direction $x_1$.

Profiles of the mean streamwise velocity at $x_1/h = 10$ for all three static grid cases are shown in Fig. 20a, and compared with experimental measurements [32, 53]. The location is downstream the onset of the self-similar regime ($x_1/h = 7$) as shown in DNS [17] and experiments [59, 60] of a planar jet with similar inlet fluctuations and Re.

The mean streamwise velocity profile for S1, the coarsest grid resolution, shows slightly different behavior compared to S2 and S3, which have two and four times more grid cells in each direction, respectively. Because of the coarse resolution, the shear layer is not fully resolved, yielding a slightly higher and lower gradient near and far from the centerline, respectively. On the other hand, the S2 and S3 profiles agree well with the experimental measurements. There are slight discrepancies between the present simulation cases S2 and S3, and experimental measurements near $x_2/\delta _{1/2} = 2$ where the velocity approaches zero. However, the velocity magnitude is quite small compared to the centerline velocity, and the two experimental measurements show divergent behavior.

The mean out-of-plane vorticity profiles $\langle \omega _3\rangle$ at $x_1/h = 10$ are shown in Fig. 20b. Similar to the mean velocity profiles, the S1 case does not have a sufficient number of grid cells to resolve the shear layer. The peak vorticity occurs near $x_2/\delta _{1/2} = 0.5$, which is inward compared to the expected vorticity peak around $x_2/\delta _{1/2} = 1$. The S2 and S3 cases have peak vorticity near the expected location. There are some differences between the S2 and S3 cases, where the most resolved S3 case (with sufficient grid cells to be considered a DNS) produces a smooth profile and a vorticity of zero at the origin. Overall, the mean statistics indicate that the S1 case is severely under-resolved, while the grid resolution for the S2 case is able to capture a majority of the mean statistics compared to the fully resolved S3 case.

Fig. 20

Profiles at $x_1/h = 10$ of the a mean streamwise velocity $\langle u_1 \rangle - U_\infty$ normalized by the centerline velocity $U_c - U_\infty$ compared with experimental results from Ref. [32] (circles) and Ref. [53] (squares) and b mean out-of-plane vorticity $\langle \omega _3\rangle$ normalized by the jet half-width $\delta _{1/2}$ and centerline velocity

Figure 21a, b shows further comparisons of the static grid cases with the root-mean-square (RMS) streamwise and transverse velocity profiles, respectively, at $x_1/h = 10$. Similar to the mean velocity, the RMS streamwise velocity profile for the S1 case is significantly different from the other cases and the experimental measurements due to the inadequate resolution. The RMS velocity also indicates further discrepancies between the S2 and S3 cases, where the S3 case aligns with the experimental measurements consistently better. Furthermore, the RMS transverse velocity shows how the grid resolution significantly affects statistics, especially in the transverse direction.

Fig. 21

Profiles at $x_1/h = 10$ of the a RMS streamwise velocity $\langle u^\prime _1 u^\prime _1 \rangle ^{1/2}$ and b RMS transverse velocity $\langle u^\prime _2 u^\prime _2\rangle ^{1/2}$ normalized by the centerline velocity $U_c - U_\infty$ compared with experimental results from Ref. [32] (circles) and Ref. [53] (squares)

Effect of snapshot matrix size

The effect of the size of the snapshot matrix on the sensor location is analyzed in Fig. 22. Three new cases with a snapshot matrix size of $N=50, 200$ and 400 are used with a base grid resolution equal to that of the R1 case and compared with the R1 case, where the snapshot matrix size is $N=100$. Figure 22a shows the eigenvalue spectrum of each case. Similar to the eigenvalue spectrum in Fig. 4a, there is a significant energy drop after the first mode for each case. This is significant because the low-rank representation of the jet can be approximated with a relatively small snapshot matrix. The sensor placement of the grid refinement tagging is shown for the four cases in Fig. 22b. As the number of POD modes increases with the snapshot matrix, the placement of the sensors begins to scatter over more of the domain. However, the placement remains concentrated at the inflow and outflow boundaries and at downstream locations ($x_1/h > 6$). Moreover, the number of total grid cells produced for each case is relatively similar. This is due to the efficiency of the load balancing and grid clustering algorithms employed to refine the grid cells near the sensors. On the other hand, if there is a desire to increase the grid resolution in the larger areas around the sampling points, modifications to the grid clustering and load balancing can be utilized to refine the grid further away from the sampling points.

Fig. 22

The a singular value spectra and b sensor placement for the grid cell for refinement with the instantaneous out-of-plane vorticity based on the size of the snapshot matrix ($N = 50, 100, 200,$ and 400). Overall grid resolution is based on the R1 case

Comparison to leverage score sampling

In this section, we introduce leverage score sampling [38, 39, 49], an alternative sampling approach to selecting columns of a matrix. The leverage scores $l_i^{(K)} \in {\mathcal {R}}^{N}$ of a rank-K matrix $\varPsi _K \in {\mathcal {R}}^{N \times K}$, which contains the top right singular vectors of $U \in {\mathcal {R}}^{N \times M}$, is given by

$$\begin{aligned} l_i^{(K)} = \left\| [\varPsi _K]_{i,:} \right\| _2^2, \end{aligned}$$

(28)

where i denotes the ith row of $\varPsi _K$. The sampling points are determined by sorting the leverage scores $l_0^{(K)} \le l_1^{(K)} \le \dots \le l_{N-1}^{(K)}$ and finding the index $c \in \{0,1,\dots ,N-1\}$ such that

$$\begin{aligned} c = {\mathop {{{\,\mathrm{arg\,min}\,}}}\limits _x} \left( \sum _{i=0}^c l_i^{(K)} > \theta \right) , \end{aligned}$$

(29)

where $\theta = K - \epsilon$, for $\epsilon \in (0,1)$ is the stopping criteria [49]. The stopping criteria is such that $c = \min (c,K)$. The sampling points are the locations of the top c columns of $\varPsi$. The points obtained via leverage score sampling are determined in a similar manner as Algorithm 1, the basis $\varPsi$ is obtained using the POD of a dense snapshot matrix. The leverage scores are (i) computed with Eq. (28), (ii) sorted via a quicksort algorithm, and (iii) selected as the top c scores found with Eq. (29). The total complexity is as follows: step (i) is ${\mathcal {O}}(NK)$, step (ii) is on average ${\mathcal {O}}(N \log N)$ or worst case ${\mathcal {O}}(N^2)$, and step (iii) is ${\mathcal {O}}(c)$. In general, for the types of problems presented herein, $K = M$ and $N \gg M$, such that the complexity of computing the leverage score sampling is similar or more operations than the sparse sampling algorithm.

The leverage score sampling is tested on the supersonic flow over a forward facing step presented in Sect. 4.3. Figure 23 shows a comparison of the location of the sampling points obtained by sparse sampling and leverage score sampling. The locations are relatively similar, but the sparse sampling selects a few locations downstream.

Fig. 23

Pressure $P/P_1$ contours with grid lines for base grid and locations sampling points using a the sparse sampling method and b leveraged score method ($|\nabla P| > 100$)

The error with respect to grid size is shown in Fig. 24a. At the lowest grid resolutions, the sparse sampling outperforms leverage score sampling. The errors converge at higher grid resolutions. At the highest grid resolution, the sparse sampling method slightly outperforms against the leverage score approach. The total number of cells after clustering and cell refinement produced by each algorithm is shown in Fig. 24b. The leverage score sampling produces significantly more cells at low grid resolution, but as grid resolution increases, the leverage score sampling produces slightly fewer compared to the sparse sampling method. Overall, at lower resolutions, the sparse sampling method performed well, but at higher resolutions, both algorithms perform similarly.

Fig. 24

a L2 error with respect to the ratio of the truth case grid size to case grid size $\Delta x_s/\Delta x_i$ and b total number of grid cells of the steady-state condition of the forward facing step simulations using the sparse sampling method (solid) and leveraged scores method (dashed). The red line shows $(\Delta x_s/\Delta x_i)^2$ (color figure online)