Yes, recently the DG formulation has received a large impulse and many papers in relevant journals appeared. I recently focused on some approaches using DG for DNS/LES/ILES. I don't think that the relevance is mainly on parallelization

I invite Luca Bonaventura and Antonella Abbà to join the discussion and give you details as experts in such approach

Shortly, DG methods can be interpreted as FV methods with an in-built subgrid reconstruction possibility, which allows them to reach high order accuracy without need for large stencils. These ensure parallel efficiency gains not only on GPUs, but also on more conventional architectures. For turbulence studies, DG allows for a very clean and straightforward implementation of Variational Multi Scale (VMS) approaches, as well as with its combination with degree adaption techniques, as we tried to point out in the papers

Article Dynamic models for Large Eddy Simulation of compressible flo...

Article A locally p-adaptive approach for Large Eddy Simulation of c...

DG can handle complex geometries and, in most recent developments, even almost arbitrary polytopal meshes. One of the most often encountered problems is the need for overintegration (e.g., use of quadrature rules and quadrature nodes of higher order than strictly necessary according to theory) when solving some nonlinear equations. This problem and some analysis and cures for it are discussed in works by Munz, Kopriva and Gassner. Also the treatment of diffusive terms requires extra stabilizations, but the underlying theory, due to Arnold, Cockburn, Brezzi and others, is now pretty well established. More recent developments of interest are also the so called hybridizable DG methods.

The interest in DG methods for LES has experienced a dramatic growth over the last 5-10 years. I would argue this was originally due to the fact that the numerical results were surprisingly accurate, but now there is a much better understanding of why this is the case. Roughly speaking, the numerical dissipation properties in wavenumber space of DG methods are very well-suited for transition prediction and for LES of moderately-high-Reynolds number turbulence [1] (*). Not surprisingly, it has been acknowledged that the Riemann solver in DG methods (which is responsible for the numerical dissipation) plays the role of an implicit subgrid-scale model (see e.g. [2,3,4]) and this implicit model somehow behaves like a dynamic model (it vanishes for laminar flows, etc.)

Regarding implementation, theoretically, the potential of DG methods for the incoming computing architectures such as GPUs is huge (**), but one needs to be clever and think out of the box. For example, state-of-the-art algorithms that are great for CPUs (pseudo-time marching, Jacobi-Newton-GMRES, etc.) can be extremely inefficient for GPUs. With the right choice of algorithms (some of them perhaps need to be new), and for some applications, however, the speedup in terms of "accuracy vs. time-to-solution" and "accuracy vs. energy consumption" with DG methods could be massive (as massive as the uncertainty on whether this potential will ever be realized in practice).

(*) The numerical dissipation properties of DG methods are arguably beneficial for LES but not for RANS-type simulations. Since DG methods are more expensive than 2nd-order FV, and there is little or no benefit in the RANS setting, that's probably why the use of DG methods in CFD was so marginal in the past.

(**) This is true not only for DG but also for the majority of high-order methods.

[1] P. Fernandez, R. Moura, G. Mengaldo, J. Peraire, Non-modal analysis of spectral element methods: Towards accurate and robust large-eddy simulations, arXiv preprint arXiv:1804.09712

[2] A.D. Beck, D. Flad, C. Tonhäuser, G. Gassner, C.-D. Munz, On the Influence of Polynomial De-aliasing on Subgrid Scale Models, Flow Turbul. Combust. 97 (2) (2016) 475–511.

[3] P. Fernandez, N.C. Nguyen, J. Peraire, On the ability of discontinuous Galerkin methods to simulate under-resolved turbulent flows, arXiv preprint arXiv:1810.09435

[4] J. Manzanero, E. Ferrer, G. Rubio, E. Valero, On the role of numerical dissipation in stabilising under-resolved turbulent simulations using discontinuous Galerkin methods, arXiv preprint arXiv:1805.10519

One extra advantage of DG for LES is, in my opinion, the fact that Variational Multiscale (VMS) approaches can be very easily included in a modal DG framework. This was first highlighted by Collis several years ago, see e.g.

Collis SS, Chang Y. The DG/VMS method for unified turbulence simulation, AIAApaper 2002; 3124:24–7.

and we have tried to pursue this further e.g. in

Article Dynamic models for Large Eddy Simulation of compressible flo...

Concerning the implicit diffusion introduced by the Riemann solvers, I am less confident than Pablo about the fact that it is generally sufficient to provide an acceptable representation of the unresolved scales. It can be seen that in many situations (although I agree not always) 'no model' simulations which only rely on the numerical flux for subgrid scale modelling are inferior to standard Smagorinsky or dynamical LES closures, this was one of the points of the above paper with Antonella Abbà.

I think that an issue about the ILES idea should be highlighted. Any CFD code used on a grid in under-resolved conditions produces a "no-model LES". It is not uncommon that such procedures gives superior solutions than simply eddy viscosity SGS model. The key is that we should represent the effects of the unresolved scale on the resolved one across the filter cut-off. In LES the filter cut-off lies ideally in the inertial range therefore an excess of eddy viscosity contribution in the explicit SGS model or numerical viscosity implicitly acting in the local truncation error in ILES, can alter this energy mechanism. Often, adding a scale similar model help to have a more correct energy transfer. Of course, one should demonstrate that in ILES such a behaviour, mimiking the scale similar transfer, exists in the local truncation error of some mehtod. I am not aware of papers that showed other than the implicit numerical viscosity, maybe you that are expert in DG approach can give some insight about that.

Luca Bonaventura, I totally agree on both points. Regarding the need/no need for an explicit model, a way I like to think about it is in terms of "spectral numerical dissipation" vs. "actual spectral eddy viscosity". This interpretation has some limitations but it's not a bad starting point. From this perspective, DG methods without an explicit model can be expected to perform well for transition prediction and aeroacoustics, as it seems to be the case in practice, due to their low numerical dissipation at large scales. For LES of low- and moderate-Reynolds-number turbulence, the exact subgrid-scale (SGS) dissipation at large scales is also small and a similar logic applies. On the other hand, the SGS dissipation at large scales is not negligible in the high Reynolds number limit, and this makes the dissipation properties of high-order methods less attractive, and an explicit model arguably needed, for very-high-Reynolds number turbulence.

Filippo Maria Denaro, there are works on modified equation analysis for DG (e.g. [1]), but one issue is that a large number of truncation terms need to be retained for under-resolved computations; which makes the analysis complicated. The non-modal analysis we have recently proposed [2] can be extended to study how energy is transferred between wavenumbers; which can give critical insights on the properties of the energy cascade built-in in the scheme.

[1] R.C. Moura, S. Sherwin, J. Peiró, Modified Equation Analysis for the Discontinuous Galerkin Formulation, In: Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2014 375-383. https://www.springer.com/us/book/9783319197999

[2] P. Fernandez, R. Moura, G. Mengaldo, J. Peraire, Non-modal analysis of spectral element methods: Towards accurate and robust large-eddy simulations, https://arxiv.org/pdf/1804.09712.pdfPreprint Non-modal analysis of spectral element methods: Towards accu...

I cannot add more than what has been discussed in the valuable comments on this question. However, I refer to the following papers [1],[2] as a first attempt to establish some comparisons between DG and FD, CD and DRP schemes. In the first paper[1] we have also proposed a new Fourier analysis formulation (combined-mode analysis) that accounts for all the eigenmodes of a particular DG or multi-degree of freedom scheme in studying their dispersion and dissipation properties. It was found that there is a lack of monotic dissipation behavior with respect to the wavenumber which results in some ranges of less dissipation at high wavenumbers. That said, this is not only the case of DG schemes but to some extent Compact difference schemes experience the same behavior when used with RK time integration.

In addition, even for pure diffusion problems DG methods still have less than exact dissipation for high wavenumbers [3].

I think for more robust behavior high-order than 3 for RKDG methods may need some kind of high wavenumber filtering. That is for robustness rather than accuracy.

[1] M. Alhawwary, Z.J. Wang, "Fourier analysis and evaluation of DG, FD and Compact difference methods for conservation laws", Journal of Computational Physics 373, Article Fourier analysis and evaluation of DG, FD and Compact differ...

[2] M. Alhawwary, Z.J. Wang, "Comparative Fourier Analysis of DG, FD and Compact Difference schemes", in AIAA 2018 Aviation forum, 2018 Fluid Dynamics ConferenceAt: Atlanta, GA, Conference Paper Comparative Fourier Analysis of DG, FD and Compact Difference schemes

[3] M. Alhawwary, Z.J. Wang, "On the Accuracy and Stability of Various DG Formulations for Diffusion", Submitted to journal, preprint arXiv:1810.03095Preprint On the Accuracy and Stability of Various DG Formulations for Diffusion

These papers are part of the project " Fourier Analysis of High-Order Methods in the Context of LES", https://www.researchgate.net/project/Fourier-Analysis-of-High-Order-Methods-in-the-Context-of-LES, and preprint copies are available online.

Additional interesting paper on the effect of SGS models is:

[4] Y. Li, Z.J. Wang, "A Priori and a Posteriori Evaluations of Sub-grid Scale Models for the Burgers’ Equation", Computers & Fluids Volume 139, https://www.sciencedirect.com/science/article/abs/pii/S0045793016301190

Abolfazl Karchani

1) Yes this is indeed a challenge for high-order methods, yet there are some solutions right now. Gmsh can generate high-order meshes and other software are in the way. There is also a handy tool developed in our group that can curve the wall elements for you if you have an original linear mesh. This tool is called MeshCurve: https://sites.google.com/site/meshcurvesoftware/meshcurve

5) and part of 6), In my no. [3] paper in the previous comment I tried to provide such a detailed flux formulation for implementation in 1D. This is very close to 2D quadrilateral and 3D hexahedral grids. For other element types it is indeed a little bit complicated. I intend to provide more details in the revised version of that paper but still in 1D. In addition, I have done some work on 2D and I may publish that some time with more details.

Regarding, the LDG and CDG, I think in 1D they are exactly the same and in 2D and 3D CDG provides a natural extension to what was intended in the proposal of the original LDG method. That is keeping the compactness of the scheme. It is yet a bit tricky to do this implementation in higher dimensions because of the switches used.

8) I kindly refer you to the following implementation of RKDG/LDG/HDG in a MATLAB/OCTAVE framework that can be useful for implementations.

https://github.com/FESTUNG/project

I agree with your point Abolfazl Karchani

At the risk of being moderately off-topic...

Luca Bonaventura, I don't get why people say that high-order FVM doesn't parallelize well. This isn't a gripe at you personally, it's a general gripe: lots of people say this.

The communication costs for reconstruction and flux evaluation are negligible / easy to hide; more data needs to be passed, sure, but not so much that bandwidth starts to matter. Communication costs for matrix assembly are harder to hide, but still relatively small. Unless there's evidence I don't know of about the parallel efficiency of the linear algebra for implicit solvers being different for FVM versus FEM / DG because of the differences in matrix structure, I don't see any basis for saying FVM doesn't parallelize well.

Our experience says that weak scaling is pretty close to perfect from two cores out as far as we're tested. IMO, weak scaling is a better model for the way CFD users use big iron: give me more resources, and I'll solve a bigger problem in the same time, not the same problem faster.

Having said that, I really need to do the rigorous numerical experiments to support what I just said (and what our experience shows to be roughly true), and write the paper. :-)

My opinion is that because the high-order mesh is not that easy to obtain for a general complex shape. The DG solver needs a curved mesh as an input, so a CFD commercial software company probably lack of motivation to add the DG solver as a part of the package since there is no general way to generate an input mesh. Assume a researcher figures out how to operate an open-source DG solver. He still need to create a high-order mesh but how? He may need to figure out how to operate an open-source high-order mesh generation software which may only handle a limit number of shapes. And these shapes may not be what this researcher want. So as a researcher focus on some CFD applications, he may just chose FVM to save a few months, even DG methods may provide a better result for his application.

Dear Carl, I understand your point and I believe myself that high order FV methods and DG are conceptually very similar. Indeed, I did not say that FVM do not parallelize well. Only, in DG methods the arrangement of the degrees of freedom is such that less communication is required, especially for hyperbolic problems, for which only interelement fluxes need to be exchanged.

If parabolic terms are included, the amount of communication required is larger, but on the other hand achieving an accuracy of second order and above for diffusion terms is more straightforward for DG than for FV methods. Comparisons between DG and FVM methods have been carried out by several

groups, one example is reported here

https://link.springer.com/chapter/10.1007/978-3-319-40528-5_17

I agree that DG is mathematically rich but like any other high-order methods DG is expensive and requires a significantly fast time discretization for practical CFD problems. In this regards, we have shown some obvious advantages in applying our exponential time schemes to DG both for steady and unsteady problems.

For 3D complicated industrial CFD problems, I don't think DG will do better than FVM method.

I would again highlight that FVM can be seen in the framework of the same weak formulation used in FEM. Both formulations share the integral theoretical framework

I would like to point out that high-order FDM does not have a consistent functional space interpretation. Sliding neighbour stencils are used in FDM, whereas DG has a constant functional space support. This may lead to the instability of high-order FDM.

not getting results

There have been some excellent answers in this thread. I did not know about DG naturally giving a dynamic subgrid model.

This is a very broad question specifically about the limitations of DG in CFD, so I am going to try and present a view which is as general as possible. DG is a way to improve the local accuracy of the discretized space by increasing the polynomial order within the element, and FVM then becomes a special case with these polynomials being piecewise constants for FVM.

The expense of the methods would be highly dependent on which PDEs within the realm of CFD you are trying to solve. If you are going for a compressible PDEs which are much closer to purely hyperbolic PDEs (Euler equations are hyperbolic), one could use fully explicit timeschemes and one only has to invert the mass matrix of the temporal terms which can be lumped and locally inverted cheaply, so everything is local, and expense is not that high. This is an close to ideal scenario as the heavy lifting of higher order polynomials is completely local to the element and possibly be offloaded to GPUs. However, we have to keep in mind that the CFL restriction controls the timestep and gets more and more stringent as we increase the order of the polynomials in the DG discretization.

Now, if we are interested in solving incompressible Navier-Stokes with finite diffusion. There are two main challenges in my view with a DG formulation. Because of the pressure-velocity coupling you would either have to use stabilized methods (e.g VMS like Luca Bonaventura mentioned can be easily formulated) or Pressure Poisson solve. Both of these approaches would require you to assemble the matrix to use iterative techniques to solve it. This is a real challenge, because as you increase the order of your DG formulation the number of degrees of freedom explodes, and the cost of matrix assembly becomes huge. You can partially circumvent this problem by using Partial assembly methods, but that has its own limitations.

The second challenge is handling the higher order diffusion terms (because it is difficult to handle derivatives of functions on the boundaries of the elements), like Luca Bonaventura mentioned there has been a lot of work (ideas like Symmetric Interior Penalty Method) on this over the last thirty years by many influential names which will take a long list to mention (Di Pietro and Ern 's book on DG talks about this in some detail). Now, if you write a time-scheme where you treat diffusion implicitly and all the other terms explicitly (these are quite common IMEX schemes) you would need to invert a matrix which goes back to the challenge of pressure Poisson. In general if your timescheme requires you to invert a matrix once or multiple times, it is going to become expensive with DG because of the sheer number of degrees of freedom as we increase the order of the polynomials. If you decide to treat diffusion terms explicitly then you have an even more restrictive timestep condition than CFL which will get worse as you increase order of your polynomials.

Another way to overcome this challenge (of high matrix assembly cost) is design function spaces which are tailored such that the global matrix is small and most of the work is local (e.g. Recent work from Beatrice Riviere's group, there are many more). There is lot of ongoing work on this, with one of the main challenges being implementation of these discretizations. Efficient parallel implementation is not as straightforward (this does not mean it is not doable). This opens up a whole new discussion though and I am going to stick to standard DG spaces.

By the way increasing spatial order also creates other problems (local oscillations) if you are solving for variables which have positivity/realizability requirements (concentration, temperature), needing you to use efficient limiters (also a very rich field). However that is also a discussion for another thread.

Therefore, coming up with an efficient highly parallel implementation of a standard DG scheme is possible for CFD problems but it has its own challenges which vary depending on the PDE you are trying to solve. This is a very rich field with a lot of challenges in the mathematical aspects and the computer science aspects.

Sorry for the long reply.

-Makrand

I would like to add the HDiv mixed element approximations as an alternative. They are considerably more precise than DG methods and in their condensed form have high parallel content. The precision of DG per dof is very disapointing.

My view: dG is of course a very good method for transport and related problems. Like all methods it has pluses and minuses. For the minuses: it needs inter-element jumps. This adds significantly to the data structure and complexity (unless one is doing adaptive mesh refinement). In may cases it's not significant but in cases where one is computing near machine limits, it can be. It is hard for me to believe in it as eliminating the need for a SGS model. But, it is an effective discretization for ''getting linear transport right'' and that goes a long way to getting good results. -Bill