Dynamical p-adaptivity for LES of compressible flows in a high order DG framework

We investigate the possibility of reducing the computational burden of LES models by employing locally and dynamically adaptive polynomial degrees in the framework of a high order DG method. A degree adaptation technique especially featured to be effective for LES applications, that was previously developed by the authors and tested in the statically adaptive case, is applied here in a dynamically adaptive fashion. Two significant benchmarks are considered, comparing the results of adaptive and non adaptive simulations. The proposed dynamically adaptive approach allows for a significant reduction of the computational cost of representative LES computation, while allowing to maintain the level of accuracy guaranteed by LES carried out with constant, maximum polynomial degree values.


Introduction
During the last decades, Large Eddy Simulation (LES) has been extensively studied and has allowed to achieve results comparable to those of Direct Numerical Simulations (DNS) with a significant reduction in computational cost, enabling the use of high accuracy models to simulate turbulent flows beyond the extremely idealized configurations to which DNS was typically limited. While the increase of the available computing power makes even more fluid flow configurations amenable to LES approaches, two interlinked difficulties prevent LES from being a widespread tool for industrial applications. The first major issue is that LES simulations, due to their space and time resolution requirements, are still very computationally expensive. The second is that selecting the right spatial resolution for an LES of a unsteady turbulent flow in a complex geometry is not a easy task. On one hand, the resolution must be sufficient to solve the equations with adequate accuracy, but limited to contain the computational costs. On the other hand, the implicit LES (ILES) commonly employed uses the discretization as the filter, so that the mesh resolution sets implicitly the filter width and the threshold between the resolved scales and the modelled ones. Moreover, in a turbulent flow in complex geometry, turbulence characteristics change in time and in space throughout the computational domain and cannot be as easily estimated as it was the case in the simple geometries employed for many years in LES studies. Therefore, an adaptive LES, where the resolution and the filter width are automatically adapted to the flow conditions, would be of great benefit for the solution of both problems, in order to have an LES with the correct filter sizing and an efficient spatial discretization. The problem of estimating the right resolution is further complicated by the fact that, as pointed out by the authors in [34], using standard error estimators to adapt the resolution is not the optimal approach for LES, if refinement to DNS resolutions is to be avoided, so that specific strategies must be devised in this context.
The need of an adaptive LES was first stated by [26], but not many attempts have been made in this direction since then. In the numerical solution of PDEs, however, in the last decades considerable efforts have been devoted to the development of refinement indicators. Most of these are based on local estimates of the discretization error, see e.g. [4,5,20,28,29,37,36,35,39]. In the LES context, however, simply increasing the resolution in order to decrease the error leads to a DNS solution [23,29], which is in contrast with the goal of adaptive LES, which consists in adjusting the resolution in order to directly resolve only a prescribed amount of the turbulent scales. An alternative approach, pursued e.g. in [17,19,18], employs indicators aimed at obtaining a good representation of some quantity of interest, rather than of the physical scales of the turbulence. Adjoint-based adaptation methods proposed by other authors [11] have also shown good performance, but entail significant computational costs.
In the present work, we will consider a physically based refinement indicator proposed by the authors in [34,31], which is especially suited for LES. This adaptivity approach is implemented in a numerical framework based on the Discontinuous Galerkin (DG) method, that was already presented and validated in [1], which allows to extend seamlessly LES concept to unstructured meshes. DG methods combine the high order of accuracy and low dissipation/dispersion properties with good parallel performances, since most of the computations are local to the element and only interelement boundary fluxes must be communicated. The adaptation technique we employ is a p−adaptive approach, in which the polynomial degree is varied locally, rather than the mesh size, as in more standard h−adaptive methods. DG methods provide an interesting environment for p−adaptivity, since they do not require to enforce continuity constraints at the interelement boundaries. Furthermore, p−adaptive techniques are appealing since they allow to correct possible shortcomings of the computational mesh, as well as to perform dynamically adaptive simulations without extensive remeshing. Static polynomial adaptivity was applied to LES in [10,34,24,31] to efficiently simulate statistically steady phenomena, while an example of dynamical adaptation is described in [11] in an ILES context.
In this paper, we will show how the adaptive method of [34,31] can be applied locally and dynamically in time, in order to simulate efficiently and accurately transient phenomena. A number of time dependent numerical simulations of two different configurations is considered, namely a temporally evolving mixing layer and the interaction of a vortex impinging on a square cylinder with the cylinder wake. The results obtained show that a significant reduction of the computational cost of LES can be achieved, reducing the number of degrees of freedom employed to values of up to 50% of those required by constant maximum polynomial degree simulations, while maintaining essentially the same level of accuracy.
In Section 2, the model equations and numerical method are reviewed. In Section 3, the dynamical adaptivity approach we propose is outlined. In Sections 4, 5, results of numerical simulations on two significant benchmarks are presented. Finally, in Section 6 some conclusions are drawn and some perspectives of future development of this approach are discussed.

Model equations and numerical method
In this Section, the model equations and the numerical method will be briefly illustrated. For a more detailed description, we refer to [1], where this methodology was first introduced and validated for LES applications. The LES filtered Navier-Stokes equations for compressible flows can be written in dimensionless form as: where U = [ρ, ρ u T , ρ e] T are the prognostic resolved variables and F = F ca − F v + F sgs are the fluxes, composed by the advective fluxes the viscous fluxes and finally the subgrid fluxes Here, · and · represent the grid filter and the Favre filter operators, so that p denotes the filtered pressure, ρ h = ρ e+p the resolved enthalpy and σ and q the momentum and heat diffusive fluxes, respectively. Equations (1) are completed by the dimensionless state equation for a perfect gas where ‹ T is the filtered temperature. To close the system the constitutive equations must also be specified: where the rate of resolved strain tensor is defined as and the dynamic viscosity, according to a power law, is with α = 0.7. In the Local Discontinuous Galerkin approach [8], the equations are rewritten as a first order system where ϕ = [ u T , ‹ T ] T . A tessellation T h composed of non overlapping tetrahedral elements is defined in the domain Ω over which a discontinuous finite element space V h is defined in the following way where P p K (K) denotes the space of polynomial functions of total degree p K . The weak discrete form of equation (9) is derived using test functions belonging to the same space as the numerical solution. A modal DG approach is here applied, by using a hierarchical orthonormal polynomial basis for each element K in the finite dimensional space V h to represent the numerical approximation of the generic variable a 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 p K (K) of degree p K , defined in R 3 as: To compute the advective flux {F} a the Rusanov flux [21] is employed while a centred flux is used for the diffusive fluxes {F} v , {F} sgs and for the gradient variables flux {ϕ}.
In the present framework, as explained in detail in [1], the filtering operators are built in the DG spatial discretization, so that the LES filtering · is equivalent to the projection onto the employed finite dimensional solution subspace while Favre filter operator is defined as In the simulations presented in this paper, two different subgrid models have been employed, the classical Smagorinsky model, in which the filter width is determined as and the anisotropic dynamic model [1,3], in which the filter width is included in the coefficients to be dynamically determined [2] using the Germano identity [15].

A dynamically p−adaptive approach
A physically based refinement indicator especially suited for LES has been proposed by the authors in [34,31]. This indicator is based on the classical structure function where · represents the expected value operator. The structure function has been widely used to study turbulence statistics and it is known to be directly related to subgrid stresses [6,7,14]. Larger values of the structure function calculated inside the element denote a poorly correlated velocity field, and the need of a higher resolution, while lower values denote a highly correlated velocity field, which is an indication of a very well resolved turbulent region or laminar conditions, hence the possibility of employing a lower resolution. The degree adaptation indicator is then defined as: where D iso ij is the structure function in isotropic conditions [25]. In [34], this procedure was shown to be effective in a statically p−adaptive framework, producing accurate results with a significant reduction in the computational cost. However, static adaptivity presents some limitations, since the resolution is fixed at the beginning of the simulation and constant in time. In the simulation of a transient phenomenon, for which a time resolved solution is sought rather than a statistical average, employing a dynamic adaptivity framework is necessary. In our implementation of the dynamic degree adaptation, indicator (18) is computed at runtime at time intervals ∆t i and averaged over time. Once the indicator has been averaged on a time interval ∆t a which is sufficiently larger than ∆t i , but still small with respect to the time scale of the motion of the main turbulent structures, a new polynomial distribution is computed based on the indicator values and the flow field approximation is updated, by either reducing or increasing the number of degrees of freedom employed in each element.
In this work, the admissible polynomial degrees ranged from 2 to 4. Two indicator thresholds 1 < 2 were employed. Cells with a value of the indicator (18) lower than the smaller threshold were assigned the lowest allowed polynomial degree 2, those with value higher than the largest threshold were assigned the highest polynomial degree 4 and the intermediate ones were assigned degree 3. Polynomial degrees were only allowed to be increased or decreased by one degree at each adaptation time. Since the solution is represented by a hierarchical basis, when lowering the polynomial degree, the contribution associated to the removed modes is simply discarded, while when raising the polynomial degree the coefficients of the newly added modes are initialized to zero, to be updated by the subsequent time evolution.
It should be remarked that the procedure outlined above can potentially lead to unbalances between the computational load of different processors in a parallel run. The implementation of a full dynamic load balancing is beyond the scope of the present work, whose focus is mostly on the assessment of the accuracy of the above methodology. As a consequence, in all the numerical experiments presented below, the computational load among the processors is approximately balanced at the beginning of the simulation only, which is clearly suboptimal with respect to full parallel efficiency. Indeed, the results will not in general be assessed based of the CPU time required, but rather on the total number of the degrees of freedom employed in the adaptive versus non adaptive simulations. However, it must be noted that,t even with this suboptimal configuration, the adaptive simulation always led to a net reduction in computational effort. In order to assess the effective reduction in the CPU time required, further work on code optimization is needed, along with the development and application of load balancing approaches such as those employed e.g. in [11,40]. 8

Temporally evolving mixing layer
In a first assessment of the adaptive LES approach outlined in the previous sections, an isothermal time developing mixing layer was simulated. The temporal mixing layer represents an interesting test case for LES, due to its simplicity and advantages from the computational viewpoint, as well as the complexity of the physics related to the mixing, see e.g. the discussions in [13,30,38,41]. The commonly employed configuration to represent a mixing layer flow was employed, imposing periodic boundary conditions on the mean flow and in the span wise directions. Thanks to the periodicity in the stream wise direction, uncertainties related to the imposition of inlet and outlet boundary conditions are avoided. In the simulations presented here, the anisotropic dynamic subgrid model [1,3] has been used, in which the filter width is included in the coefficients to be dynamically determined. A further application of this approach is presented in [2].
A sketch of the flow configuration is shown in Figure 1. The mixing layer is characterized by two parallel flows with different velocities U 1 and U 2 . The convection velocity of the isothermal mixing layer, defined as U c = (U 1 + U 2 )/2, is the velocity that transports the cortical structures at the centreline. In the configuration we considered, the convection velocity is assumed to be zero, so that the eddies do not travel inside the domain by means of the flow, but only move by mutual interaction. As a consequence, U 1 = −U 2 and |U 1 | = |U 2 | = (U 1 − U 2 )/2. The reference frame is chosen so that the x axis is aligned with the main direction of the flow centred at the middle of the domain, while the y and z axis denote transverse and lateral directions, respectively. The initial vorticity thickness, defined as is chosen as reference length for the non-dimensionalization. The initial velocity jump ∆U r = (U 1 − U 2 ) is used as reference velocity, the initial uniform density ρ r and temperature T r are used as reference density and temperature, respectively. The fundamental dimensionless groups are then The temporally evolving mixing layer develops from a specified initial condition. In these computations, we use a hyperbolic tangent where U 1 = 0.5 and U 2 = −0.5 so that U c = 0. Similarly to [16,12], a 3D incompressible disturbance is added to the base velocity profile to initiate the transition process. It consists of harmonic disturbances expressed as where A = 0.025(U 1 −U 2 ) and the decay in y-direction is governed by σ = 0.05. Moreover, a 3D white noise perturbation is added to the initial velocity field, generated for each mesh node by a logistic map rescaled to take values on [−A, A] and multiplied by the attenuation factor e −σy 2 . Uniform pressure p = 1 and density ρ = 1 are used to initialize the simulation. The longitudinal dimension of the computational domain must be large enough to allow at least for the merging of two principal vortical structure, so it must be taken as a multiple of the wavelength λ a characterizing the most unstable perturbation, which according to the analysis in [30] is given in this context by λ a = 7.66. In the present simulation, the size in the streamwise direction has been chosen L x = 8×λ a = 60.8, so that up to three vortices are allowed to merge. In the normal direction, the computational domain extends for −30 ≤ y ≤ 30, while in spanwise direction L z = 30.4. Periodic boundary conditions are applied in the statistically homogeneous x and z directions, while a sponge layer is employed at the top and bottom boundaries, so that the effective size of the computational domain in the y direction is equal to 80.
The computational mesh is built first subdividing the domain in N hx × N hy × N hz = 20 × 36 × 10 hexahedra. Each hexahedron is then split into 6 tetrahedra, yielding a total of 43200 tetrahedral elements. The actual resolution depends on the local polynomial order and can be computed as For p = 4, we get ∆ x = 5.11, ∆ y = 0.37, ∆ z = 0.5.
A preliminary study of the mixing layer in the present configuration was presented in [27], where the sensitivity of the results to the choice of spatial resolution and subgrid scale model was checked. The objective of the present work is to test the p−adaptivity technique in a time evolving flow. For this purpose, the results obtained with polynomial degree dynamically adapting in space and time are compared with the results obtained with constant fourth and third degree polynomials.
To set the values of the adaptation thresholds 1 , 2 , the structure function indicator was evaluated on the velocity field obtained with constant p = 4 at the final dimensionless time t = 150. Thresholds were then chosen in order to obtain, for that field, an adapted number of degrees of freedom (dofs) lower than 80% of the corresponding number of dofs in the constant p = 3 case. The resulting threshold values are then 1 = 10 −6 and 2 = 2 × 10 −3 . The dynamical p−adaptation was carried out by computing the indicator value at intervals ∆t i = 0.0375 and averaging them over intervals ∆t a = 0.075, while the value ∆t = 0.0075 was used for time integration.
The polynomial degree distribution, the density and the streamwise component of the velocity, at different instants, in the plane z = 0, are represented in Figure 2. In the initial condition, the elements with highest degree p = 4 are concentrated along the middle of the domain, where high resolution is required. Then, the region with higher polynomial degree extends in space during time, adapting to the merging of the vortices, to the diffusion of turbulent structures and to the growth of the layer thickness. An important quantity for  the characterization of the mixing layer is the momentum thickness, defined as The time evolution of this quantity is represented in Figure 3. Notice that slope variations in the time evolution are associated to the merging of the vortices. A difference in the behaviour of the p = 3 simulation with respect to the p = 4 one is evident after the second merging, while the p−adaptive simulations well reproduces the p = 4 results. No relevant differences between the results obtained using adaptive or constant polynomial degree distribution are observed in the mean density and velocity profiles, shown in Figure 4, nor in the mean normal stresses profiles, shown in Figures 5(a)-5(c). Instead, looking at the shear stresses in Figure 5(d), although the peak is slightly overestimated, the p−adaptive solution is in good agreement with that of the p = 4 case. From the same picture, the more dissipative character of the p = 3 solution is also apparent.
Looking at the turbulent kinetic energy profile in Figure 6(a), we can affirm that the turbulent energy amount of the flow is reproduced by the p−adaptive simulation with the same accuracy as in the p = 4 case. In Figure 6(b), the subgrid scale (SGS) kinetic energy in the p = 3 case is clearly larger than in the other cases, due to the lower resolution that requires more intense contribution from the model. The SGS energy for the p−adaptive case is closer to that of the p = 4 case, demonstrating that the two simulations have an equivalent effective resolution.
All the simulations were carried out on the Marconi cluster at     Core hours p = 4 12600 p = 3 9000 p−adaptive 9300 Table 1: Core hours required for the mixing layer simulations using uniform polynomial degree distribution p = 4, p = 3 and p−adaptive, respectively.
CINECA, using 272 KNL processors. In Figure 7, the total number of dofs for the adaptive simulation is shown as a function of time and compared with the values corresponding to uniform polynomial distributions, that are about 1.51 × 10 6 and 8.64 × 10 5 for p = 4 and p = 3, respectively. The gain using p−adaptivity is evident, since it requires less than 80% of the dofs for p = 3 and less then half of the dofs for the p = 4 case. The core hours necessary to perform the simulations are reported in Table 1. In spite of the absence of a specific dynamic load balancing procedure, the p−adaptive technique also yields a reduction of the required CPU time, which has a value only marginally larger than that of the p = 3 simulation. On the other hand, the sub-optimal nature of the present implementation is highlighted by the fact that the p−adaptive simulation required more CPU time than the p = 3 case, even though it involved less dofs.
Concluding this analysis, we can affirm that the p−adaptive solutions display almost the same accuracy as those obtained with constant degree four, but require less then the half the number of dofs and about 29% less CPU time for this test case, in spite of a sub-optimal implementation. Moreover, the structure function indicator which was proven successful in static adaptation in [34] is shown to be suitable also for simulation of transitional flows using a dynamic adaptation in time.

Interaction of Vortex and Square Section Cylinder
A body-vortex interaction flow represents another interesting time dependent problem to test dynamic p−adaptivity. Here, a vortex interacting with a square section cylinder at Re = 22000 and M a = 0.3 is considered. This configuration was chosen because the analogous case without a vortex has already been studied with statically adaptive techniques in [34]. The classical Smagorinsky SGS model was applied for this test. The results obtained with polynomial adaptivity are compared with the solution obtained with uniform in space and constant in time degree distribution. The constant degree four solution is taken as reference.
The computational domain employed is box shaped and a 2D sketch of the geometry is represented in Figure 8. Denoting the cylinder side with H, which is used as reference length, the inflow length has been taken equal to L f = 10H, while the outflow length is equal to L r = 20H. The cylinder is vertically centred with a distance of L s = 10H from the upper and lower boundaries. The domain is extruded in the spanwise direction by a length equal to L z = 4H. This configuration was employed e.g. in [9] and has been used to assess statically adaptive LES simulations of the standard flow around a square section cylinder in [34,31,32]. The mesh employed consists of 23816 tetrahedra arranged in an outer unstructured area and a structured O-grid mesh around the cylinder, which is then extruded in the spanwise direction. Considering a polynomial basis of degree 4, the resolution in space around the cylinder is ∆ n = 0.008H in the wall normal direction, ∆ t = 0.024H along the cylinder sides and ∆ z = 0.135 along the spanwise direction. As in [34], the thresholds for the structure function (SF) indicator are taken to be 1 = 5×10 −4 , 2 = 1×10 −2 , so as to obtain a number of dofs just lower than the one obtained with uniform degree 3.
To simulate the cylinder-vortex interaction, a vortex is superimposed to an initial condition corresponding to a developed flow around the cylinder, obtained by a previous LES simulation. The introduced vortex is parallel to the z direction and defined by the velocity components where C v denotes the vortex strength, R v is the vortex radius, r 2 = x 2 + y 2 is the distance in the x − y plane from the vortex centre and U ∞ is the uniform flow over which the vortex is superimposed [22]. The maximum radial velocity is obtained at r = R v and is given by v θmax = Cv R √ e . Regarding the other variables, it is known from [9] that for a viscous compressible vortex the radial pressure distribution given by the solution of Assuming a constant temperature distribution T = T ∞ , one obtains from (24) the profile Finally the density distribution follows the equation of state ρ = P/(R gas T ).
In [31,33] simple advection of this vortex was tested. The maximum radial velocity chosen is v θmax = 0.5 and the vortex radius is R v = 0.41. A reference simulation without vortex was used to calibrate the introduction time in order for the vortex to reach the cylinder in the instant of maximum lift. The dynamical p−adaptation strategy consisted in this case of adapting the polynomial degrees approximately three times during the time needed by a fluid particle to pass through the smallest elements on the cylinder walls at the free stream velocity U ∞ . This choice led to computing the indicator value at intervals ∆t i = 0.0016 and averaging them over intervals ∆t a = 0.05, while the value ∆t = 0.0004 was used for time integration.
The dynamically adaptive procedure was able also in this case to effectively represent the structures of the flow, both in the advected vortex region and in the wake of the obstacle, as can be seen in Figures 10-13. It is possible to note that a higher polynomial degree is employed in the advected vortex region and in the shear layers around the cylinder. The dynamic adaptation procedure is able to effectively increase the polynomial degree around the vortex and to follow it as it is advected downstream, leaving all the elements with no vortex activity at the lowest resolution. Notice that the centre of the vortex is represented with degree 4 polynomials in the first instants, but only degree 3 polynomials are employed later, since the vortex enters an area with a more refined mesh in which a lower polynomial degree is sufficient.
Among the different effects caused by the vortex-cylinder interaction, the most interesting are those on the forces acting on the cylinder. The history of the load coefficients before, during and after the interaction is reported in Figure 14.
In the developed flow around the square cylinder, the periodic oscillation of the force coefficients are related to the unsteadiness of the recirculation bubble on the upper and lower sides and to the vortex shedding. Reaching the upper corner of the cylinder at t = 5, the vortex changes the pressure distribution on the front of the cylinder and strongly modifies the recirculating bubble on the upper side. In this condition, it is very important to use an adequate resolution on the upper side and in the separation bubble. In Figure  14 it can be observed that the peaks of the lift and drag coefficients at t = 5 with p−adaptivity equal those in the reference solution p = 4, demonstrating the capability of the proposed dynamically adaptive approach to well represent this unsteady phenomenon. The values of the differences of the peaks for the lift and drag coefficients with respect to the vortex free solution are reported in Table 2. Also in this case the values obtained with p−adaptivity are in very good agreement with the reference solution.      configuration Table 2: Differences, in percentage, of the force coefficients with respect to the simulation without vortex, during the vortex interaction time, t = 4 − 8.5. For the lift, only the change in the subsequent peak is presented, for the drag both the difference in maximum and minimum.
The simulations were carried out on the Marconi cluster at CI-NECA, using 272 Brodwell processors. The adaptive simulation allows a reduction of about 60% of the required CPU time with respect to the constant p = 4 degree case. In Figure 15, the time history of the dofs number is shown. It is worth noting that, even though the number of dofs is affected by the flow regime inside the domain, its variations are anyway fairly limited and the total number remains less than the one employed for the uniform p = 3 simulation. The core hours necessary to perform the simulations are reported in Table 3. Also in this case,the p−adaptive technique yields a reduction of the required CPU time, which has a value only marginally larger than that of the p = 3 simulation. Again, the sub-optimal nature of the present implementation is highlighted by the fact that the p−adaptive simulation required more CPU time than the p = 3 case, even though it involved less dofs.
Core hours p = 4 756 p = 3 388 p−adaptive 462 Table 3: Core hours required for the vortex-cylinder simulations using uniform polynomial degree distribution p = 4, p = 3 and p−adaptive, respectively.

Conclusion
The need for adaptive LES approaches was first stated in [26], but adaptive LES simulations are not very common yet. DG methods provide a appropriate environment for adaptive approaches, since they do not require to enforce continuity constraints at the interelement boundaries. Furthermore, p−adaptive techniques are appealing, since they allow to correct possible shortcomings of the computational mesh as well as to perform dynamically adaptive simulations without extensive remeshing. Static polynomial adaptivity has been applied to LES in a DG context in [10,34,24,31]. An essential tool for such simulations is an adaptation criterion that, rather than simply increasing the resolution in order to decrease the error, which is known to lead to a DNS solution [23,29], tries instead to adjusting the resolution in order to directly resolve only a prescribed amount of the turbulent scales.
In this work, the physically based refinement indicator proposed by the authors in [34,31] was shown to be applicable locally and dynamically in time, in order to simulate efficiently and accurately transient phenomena. Numerical simulations of a temporally evolving mixing layer and of a vortex impinging on a square cylinder were performed. The results obtained show that a significant reduction of the computational cost of LES is achieved, reducing the number of degrees of freedom employed to values down to 50% of those required by constant maximum polynomial degree simulations, while maintaining essentially the same level of accuracy.
A main limitation of the present results is that, in this preliminary implementation of a dynamically adaptive approach, no dynamic load balancing has been performed. In future work on code optimization, load balancing approaches such as those employed in [11,40] will be considered and tested in order to achieve maximum parallel efficiency.