Author: Martin Vymazal, Imperial College London
High-order methods on unstructured grids are now increasingly being used to improve the accuracy of flow simulations since they simultaneously provide geometric flexibility and high fidelity. We are particularly interested in efficient algorithms for incompressible Navier-Stokes equations that employ high-order space discretization and a time splitting scheme. The cost of one step in time is largely determined by the amount of work needed to obtain the pressure field, which is defined as a solution to a scalar elliptic problem. Several Galerkin-type methods are available for this task, each of them have specific advantages and drawbacks.
High-order continuous Galerkin (CG) method is the oldest. Compared to its discontinuous counterparts, it involves a smaller number of unknowns (figure 1), especially in a low-order setting. The CG solution can be accelerated by means of static condensation, which produces a globally coupled system involving only those degrees of freedom on the mesh skeleton. The element interior unknowns are subsequently obtained from the mesh skeleton data by solving independent local problems that do not require any parallel communication.
Figure 1: Distribution of unknowns for continuous and discontinuous Galerkin methods.
The amount of information interchanged while constructing and solving the statically condensed system, however, is determined by the topology of the underlying grid. Unstructured mesh generators often produce meshes with high vertex valency (number of elements incident to given vertex) and CG therefore has rather complex communication patterns in parallel runs, which has a negative impact on scaling [3].
Discontinuous Galerkin (DG) methods [1], on the other hand, duplicate discrete variables on element boundaries, thus decoupling mesh elements and requiring at most pairwise communication between them. This is at the expense of larger linear system and more time spent in the linear solver. Discontinuous discretization is therefore expected to scale better on parallel computers, but the improved scaling is not necessarily reflected in significantly smaller CPU times when compared to a CG solver.
Hybrid discontinuous Galerkin (HDG) methods [2] address this problem by introducing an additional (hybrid) variable on the mesh skeleton. The hybrid degrees of freedom determine the rank of the global system matrix and HDG therefore produces a statically condensed system that is similar in size to the CG case. In contrast with CG, the static condensation in HDG takes place by construction rather than being an optional iterative technique. Similarly to the classical DG method, HDG scales favourably in comparison with CG, but the work-to-communication ratio is once again improved due to increased amount of intra-node work rather than due to better overall efficiency.
To maximize the potential of each Galerkin variant in a unified setting, we study a finite element discretization that combines the continuous and discontinuous approach by considering a hybrid discontinuous Galerkin method applied to connected groups of elements supporting a globally continuous polynomial basis. This settings leads naturally to a formulation of weak Dirichlet boundary conditions for the CG method.
Combining continuous and discontinuous discretization
Consider an unstructured grid divided into multiple subdomains. Each partition can be regarded as a ’macro-element’ with a piecewise-continuous basis, while the patches are coupled together weakly as in HDG and the scalar flux (hybrid variable) is only defined on inter-partition boundaries (figure 2). The mixed variational form of this combined CG-DG method applied
Figure 2: Combination of CG and HDG - full domain and a single-partition problem.
to groups of agglomerated elements remains formally identical to HDG. We thus obtained a Galerkin method with weakly imposed Dirichlet boundary conditions in each mesh subdomain, and limited (at most pairwise) coupling between discontinuous subdomains. When a single partition is considered, this can be further reduced to a CG method with weakly imposed Dirichlet data (Figure 2b).
Weak vs. strong boundary conditions - a comparison
We evaluated the accuracy of weakly imposed Dirichlet boundary conditions on a scalar Helmholtz problem solved in a square domain (−1,1)2 with zero Dirichlet boundary conditions. To verify the convergence rates, the problem was first solved on a series of meshes with increasing number of unknowns and P3 elements. We then repeated the test on a fixed triangular and quadrilateral grid with polynomial degrees increasing from 2 to 14.
Figure (3) shows that convergence rate in strong and weak settings is virtually identical. The solution of linear system with weak boundary conditions included is slightly slower because the Dirichlet degrees of freedom are not eliminated a-priori as in the strong case.
Figure 3: Convergence to the exact solution in L2 norm.
Our current effort is concentrated on building a theoretical performance model of the combined CG-DG scheme that would assess its feasibility in massively parallel setting.
If you want to find out more about the contents of this article, please contact Martin Vymazal, David Moxey, Chris Cantwell, Mike Kirby and Spencer Sherwin of the Imperial College London.
References
[1] Douglas N Arnold, Franco Brezzi, Bernardo Cockburn, and L Donatella Marini. Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM journal on numerical analysis, 39(5):1749–1779, 2002.
[2] Bernardo Cockburn, Jayadeep Gopalakrishnan, and Raytcho Lazarov. Unified Hybridization of Discontinuous Galerkin, Mixed, and Continuous Galerkin Methods for Second Order Elliptic Problems. SIAM Journal on Numerical Analysis, 47(2):1319–1365, 2009.
[3] Sergey Yakovlev, David Moxey, Robert M. Kirby, and Spencer J. Sherwin. To CG or to HDG: A comparative study in 3D. Journal of Scientific Computing, 67(1):192–220, 2016.
