Dynamical Degree Adaptivity for DG-LES Models

We present a first dynamically adaptive test of a physically based p-refinement criterion for DG-LES approaches recently proposed by the authors and applied so far to static adaptation only. The results, albeit preliminary, are encouraging and show that, also in the dynamically adaptive case, the proposed criterion allows to reduce significantly the number of degrees of freedom employed in a LES simulation without a major loss of accuracy.

equations and numerical discretization approach. Here, only a short description of the discretization elements necessary to introduce dynamical adaptivity will be reported. On the computational domain ⊂ R 3 a tessellation T h is defined, composed of non overlapping simplicial elements. A discontinuous finite element space V h is defined as where P q K (K) denotes the space of polynomial functions of total degree q K . The degree can vary arbitrarily from element to element, and the definition of a suitable way to assign such polynomial degree will be discussed in the following. The numerical approximation of the generic variable a can be expressed as where φ K l are the basis functions on element K, a (l) are the modal coefficients of the basis functions and n φ (K) + 1 is the number of basis functions required to span the polynomial space P q K (K) of degree q K , defined in R 3 as: It is worth noting that the expression in (2) can be rewritten, thanks to the hierarchical nature of the basis, as where d 0 = {0} and d p = l ∈ 1 . . . n φ (K) | φ l ∈ P p (K)\P p−1 (K) is the set of indices of the basis functions of degree p. Obtaining a more or less accurate approximation can be done through increasing or decreasing the limit q K of the sum over p. It is also worth noticing that the basis normalization implies that the first coefficient of the polynomial expansion a (0) coincides with the mean value of a h | K over K.
In the present DG-LES approach, as discussed extensively in [1], the LES filtering operators are built directly into the DG discretization, in a spirit similar to the VMS approach [4]. Considering V : L 2 ( ) → V the L 2 projector over the subspace V ⊂ L 2 ( ), defined by it is possible to define the LES filtering · as the projection over the finite dimensional solution subspace V h in the following way: The application of the main LES filtering is purely formal, since it coincides with the discretization of the equations. In this way, simply discretizing the equations leads to solving them for the filtered quantities. Another parameter to be defined is the filter characteristic dimension, , employed in the definition of all the eddy-viscosity based subgrid model. The definition of the filter size is constant over each element, since the projection is performed elementwise. While more refined definitions can be employed, see e.g. [2], the simple definition was employed with success. For the time discretization, the five stages, fourth order Strong Stability Preserving Runge-Kutta method proposed in [8] is employed. The numerical implementation of the previously sketched approach is built in the solver dg-comp using the finite elements toolkit FEMilaro [7]. A first attempt to introduce static p-adaptivity in a DG-LES framework has been presented in [10]. In order to overcome the limitations of classical error estimations in LES, a novel indicator based on the classical structure function was proposed. Large values of the structure function calculated inside the element denote a poorly correlated velocity field and the need of higher resolution, while a low structure function value denotes a highly correlated velocity field, which is an indication of a well resolved turbulent region or laminar conditions and of the possibility to employ a lower resolution. However, most of the subgrid models (and in particular the Smagorinsky model) perform adequately in a regime of homogeneous isotropic turbulence, if the filter cut-off length is inside the inertial range. Therefore, in such conditions excessive refinement is not necessary and one can let the subgrid scale model simulate the turbulent dissipation. For this reason, the contribution due to homogeneous isotropic turbulence is removed from the structure function (7). This contribution, as discussed in detail in [10], can be written as where r = r and D LL , D NN are the longitudinal and transverse structure functions, respectively. Once r is known, only D LL and D NN need to be determined. The procedure to compute the error indicator can then be described as follows: 1. choose a pair of points defining x and r in K 2. compute the structure function D ij (K) based on x, r and the simulated velocity field 3. compute D NN and D LL by a least square fit of (8) to the structure function values within the element 4. define the degree adaptation indicator as: The static adaptivity procedure presented in [10] is able to produce accurate results with a significant reduction in computational cost. For the simulation of transient phenomena, however, a dynamic adaptivity approach must be applied. The goal of this work, which summarizes results presented in [9], is to extend the above approach to dynamical adaptation, which was successfully employed in the inviscid case in [11,12].
In those papers, in which special time discretizations approaches were employed that allow the use of very long time steps, the adaptation process was performed at each time step. In the dynamically adaptive simulations presented here, instead, which are carried out with a relatively small time step, the structure function indicator I nd SF (K) is computed every n i (K) time steps and the average of s i (K) subsequent values of this quantity is computed. Then, every n i (K) × s i (K) time steps, based on the resulting indicator value in each element, either the polynomial degree is left unchanged or it is updated along with the solution representation. Since the solution is expressed in terms of a hierarchical basis (4), when lowering the polynomial degree, the contribution bound to the removed modes is simply discarded, while when raising the polynomial degree the contribution of the newly added mode is left to zero, to be populated when the integrals over the element and faces couple the old modes with the newly introduced ones.
Notice that, in the present implementation, no dynamic load balancing has been implemented for parallel runs. This means that, during the parallel execution, the dynamic change of number of degrees of freedom could potentially lead to unbalances between the load of different processors. At the moment the balancing is generally executed using a static polynomial distribution. While avoiding excessive unbalancing, this is definitely not the optimal approach and more effective load balancing techniques will have to be investigated in the future.

Dynamical Adaptivity Experiments
The proposed dynamic adaptation criterion has been tested in the simulation of a isolated vortex superimposed on a uniform horizontal flow [5]. This simple test has been chosen for the preliminary study reported here, in anticipation the more complex tests already discussed in [9], in which the same isolated vortex impinges on an obstacle. The DG-LES approach described in [1] was applied, as in [10], with a standard Smagorinsky model for the subgrid stresses. A coarser and a finer mesh have been employed, both based on fully unstructured tetrahedra of constant characteristic length equal to l h = 1 and l h = 0.5, respectively. The indicator (9) is computed every n i (K) = 2 time steps and s i (K) = 10 subsequent values are averaged, in order to adapt the resolution every 20 time steps. The sensitivity analysis of the results with respect to these parameters has not yet been carried out and will be the focus of future study. As in [10], two threshold values 1 , 2 are used to determine p-refinement and p-derefinement. More specifically. the cells with indicator values smaller than 1 are assigned polynomial degree 2, those with indicator values larger than 2 are assigned polynomial degree 4, while the others are polynomial degree 3. The threshold values employed are given by 1 = 1 × 10 −4 , 2 = 1 × 10 −2 . Following [10], these values were chosen so as to achieve on average a total number of degrees of freedom slightly smaller than that required by a uniform degree simulation with p = 3. The dynamic adaptation procedure is able to effectively increase the polynomial degree around the vortex and follow it as it is advected downstream, leaving all the elements with no vortex activity at the lowest resolution. A map of the polynomial degrees in the domain during the advection of the vortex is shown in Fig. 1.
The profiles of velocity magnitude recorded during time, along the path of the vortex, at different distances from the vortex starting point, employing the coarsest mesh, are presented in Fig. 2. The simulations obtained at different uniform polynomial orders are compared with the adaptive results. It can be observed that, even at the highest uniform resolution of degree 4 the velocity profile is distorted during the advection, due to the very limited grid resolution. However, the vortex does not diffuse and dissipate excessively, as opposed to the low resolution uniform The comparison with the uniform high degree simulations can be more easily observed in Fig. 3, which show the difference of the velocity magnitude profiles with respect to the uniform degree 4 results, still for the coarse mesh case. In the locations nearer to the starting position of the vortex the adaptive simulation appears close to the degree 4 solution when the first part of the vortex is passing, while a slight difference appears in the second part of the vortex, which is however always within the error of the uniform degree 3 simulation. In the locations farther from the initial starting point of the vortex, which sense the vortex passage after a longer advection time, the adapted simulation is always very close to the uniform degree 4 solution. It has to be noted that the average number of degrees of freedom of the adaptive simulation is 41,488, which remain almost constant throughout the simulation. This is 10.8% more than the 37,430 degrees of freedom needed for the uniform degree 2 solution, 44.6% less than the 74,860 degrees of freedom of the uniform degree 3 resolution, which is always outperformed by the adaptive one, and 68.3% less than the uniform degree 4 simulation.
To correctly assess the effects of adaptivity in the case of the refined mesh, we study the difference of the various results with respect to the uniform degree 4 one, presented in Fig. 4. The differences are generally very small, even for the lowest resolution, however it is possible to note how the adaptive results are always comparable to the uniform degree 3 results, and in many points better. Nonetheless, the improvement created by the adaptivity is more limited than in the coarse case, mainly due to the fact that the mesh by itself sufficient to resolve the vortex. In this case the average number of degrees of freedom of the adaptive case is 170,470, which is 5.7% more than the 161,320 degrees of freedom of the uniform degree 2 case, 47.2% less than the uniform degree 3 case and 70.0% less than the uniform degree 4 case. Also the difference in vorticity profiles between the simulation at uniform degree 4 and the lower resolution simulations are presented in Fig. 5 for the coarse resolution and in Fig. 6 for the finer resolution. By comparing the results at the two different resolution is possible to note also for the vorticity that, at the finer resolution, the large scale phenomenon is correctly represented by almost all polynomial degrees, with a minimal vorticity dissipation, while at the coarser resolution only the higher polynomial degree, as well as the adaptive simulation, avoid an excessive dissipation of vorticity. At the coarser resolution, the difference of the adaptive simulations with respect to the uniform degree 4 ones is smaller than the differences between the other uniform degree simulations (Fig. 5), showing that with the adaptation is also possible to obtain a better resolution of the vorticity profiles. The same is true also at the finer resolution (Fig. 5). In the dynamically adaptive simulations spurious acoustic waves seem to be produced by the dynamical adaptation process, see    Fig. 7. These spurious disturbances were not observed in the dynamically adaptive tests presented in [11,12], which employed an implicit time discretization, thus strongly damping these high frequency solution components. However, as it can be seen inspecting the time series of the pressure values (not reported here due to the limited space available), these disturbances decrease rapidly in amplitude on the finer mesh and do not seem to propagate through the domain but rather follow the advected vortex. This spurious feature warrants further investigation of the dynamical adaptation approach if a correct approximation of acoustic waves is desired.

Conclusions
The novel degree adaptation criterion for LES simulations in adaptive DG frameworks proposed in [10] and tested so far only in statically adaptive simulations has been also employed in dynamically adaptive simulations. Numerical results in the benchmark case of the advection of an isolated vortex have been presented. These results are meant to be a preliminary for the study of more complex configurations in which the same isolated vortex impinges on an obstacle. The presented results show that the proposed criterion is also effective in the dynamical case. With a coarse basic mesh resolution the effects of p-adaptivity are significant, leading to results close to the ones obtained with the maximum resolution allowed to the polynomial base, while when the mesh resolution is already suitable to represent the vortex even with the lowest polynomial degrees the adaptivity leads anyway to accurate results, but with an even higher reduction of the number of degrees of freedom with respect to the non-adaptive solutions. In a subsequent work, the results obtained in [9] for the case of the isolated vortex impinging on an obstacle will be presented, along with other application to fully three-dimensional turbulent flows. The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.