Implicit shock tracking using an optimization-based high-order discontinuous Galerkin method

Highlights

• Implicit shock tracking for high-order accurate resolution of discontinuities.

• DG solution on discontinuity-aligned mesh is the solution of an optimization problem.

• Objective function penalizes DG residual in enriched test space and mesh distortion.

• Novel SQP solver simultaneously converges mesh and DG solution to optimal values.

• Accurately resolves solutions with discontinuities on coarse, high-order meshes.

Abstract

A novel framework for resolving discontinuous solutions of conservation laws, e.g., contact lines, shock waves, and interfaces, using implicit tracking and a high-order discontinuous Galerkin (DG) discretization was introduced in [39]. Central to the framework is an optimization problem whose solution is a discontinuity-aligned mesh and the corresponding high-order approximation to the flow that does not require explicit meshing of the unknown discontinuity surface. The method was shown to deliver highly accurate solutions on coarse, high-order discretizations without nonlinear stabilization and recover optimal convergence rates $\mathcal{O}(h^{p+1})$ even for problems with discontinuous solutions. This work extends the implicit tracking framework such that robustness is improved and convergence accelerated. In particular, we introduce an improved formulation of the central optimization problem and an associated sequential quadratic programming (SQP) solver. The new error-based objective function penalizes violation of the DG residual in an enriched test space and is shown to have excellent tracking properties. The SQP solver simultaneously converges the nodal coordinates of the mesh and DG solution to their optimal values and is equipped with a number of features to ensure robust, fast convergence: Levenberg-Marquardt approximation of the Hessian with weighted elliptic regularization, backtracking line search based on the ${\ell }_1$ merit function, and rigorous convergence criteria. We use the proposed method to solve a range of inviscid conservation laws of varying difficulty. We show the method is able to deliver accurate solutions on coarse, high-order meshes and the SQP solver is robust and usually able to drive the first-order optimality system to tight tolerances.

Introduction

High-order methods, such as the discontinuous Galerkin (DG) method [7], [19], are widely believed to be superior to traditional low-order schemes for simulation of turbulent flow problems. However, in the presence of shocks and other discontinuities such as contact lines or interfaces, the lack of nonlinear stability proves to be a fundamental challenge. These features are ubiquitous in engineering and science applications, particularly in high-speed flow problems. Many solutions have been proposed, but most require excessively refined meshes to resolve all the features, which in practices makes it difficult to accurately predict high Reynolds, high Mach flows that feature shocks, boundary layers, and interactions between them. Therefore, new advances are required to make high-order schemes sufficiently robust and competitive for real-world problems.

Most of the techniques for addressing shocks are based on so-called shock capturing, that is, the numerical discretization somehow incorporates the discontinuities independently of the computational grid. One simple method is to use a sensor that identifies the mesh elements that contain shocks, and reduce their polynomial degrees [3], [6]. For the DG method, this essentially leads to a standard cell-centered finite volume scheme locally, which is well-known to handle shocks robustly. Related, more sophisticated approaches include limiting, such as the weighted essentially non-oscillatory (WENO) schemes [16], [24], [20]. For high-order methods, artificial viscosity has also proven to be highly competitive, since it can smoothly resolve the jumps in the solution without introducing additional discontinuities between the elements [29]. The main problem with all these approaches is that they reduce to first order accuracy in the affected elements, which translates into a globally first order accurate scheme. This can be remedied by using local mesh refinement around the shock (h-adaptivity) [12], although the anisotropic elements that are required for efficiency are difficult to generate and excessively fine elements are needed around the shock.

An alternative approach is shock tracking or shock fitting, where the computational mesh is moved such that their faces are aligned with the discontinuities in the solution [33], [34], [4], [17], [38], [37], [41], [35], [2], [31], [15], [28]. This is very natural in the setting of a DG method since the numerical scheme already incorporates jumps between the elements and the approximate Riemann solvers employed on the element faces handle the discontinuities correctly. However, it is a difficult meshing problem since it essentially requires generating a fitted mesh to the (unknown) shock surface. Also, in the early approaches to shock fitting, it was applied to low-order schemes where the relative advantage over shock capturing is smaller than for high-order methods. For these reasons, shock tracking is largely not used in practical CFD today.

In [39], we introduced a novel approach to shock tracking that does not require explicitly generating a mesh of the unknown discontinuity surface. Rather, the conservation law is discretized on a mesh without knowledge of the discontinuity surface and an optimization problem is formulated such that its solution is a mesh that aligns with discontinuities in the flow and the corresponding solution of the discrete conservation law. That is, tracking of the discontinuities is implicitly defined through the solution of the optimization problem and will be referred to as implicit shock tracking. While this approach works with any discretization that allows for inter-element discontinuities, we focus on high-order DG methods due to the high degree of accuracy attainable on coarse meshes, proper treatment of discontinuities with approximate Riemann solvers, and the ability to used curved elements to track discontinuities with curvature. The optimization problem is solved by simultaneously converging the mesh and solution to their optimal values, which never requires the fully converged DG solution on non-aligned meshes and does not require nonlinear stabilization. The combination of implicit tracking with a DG discretization is truly high-order accurate, since the solution is smooth within each element, and very accurate solutions can be obtained on coarse meshes.

This paper extends our prior work with a new error-based objective function, a sequential quadratic programming solver for the optimization problem, and a number of practical considerations. The proposed objective function penalizes violation of the DG residual in an enriched test space, which is a surrogate for violation of the weak formulation of the conservation law. Even though a traditional DG solution will oscillate about discontinuities in the solution on a non-aligned mesh, this violates the true conservation law; as a result, the proposed objective function promotes alignment of the mesh with discontinuities. This formulation has the added benefit of r-adaptive behavior; even in smooth regions of the flow, nodes will adjust to improve the approximation of the conservation law. The other main contribution of this work is an SQP solver for the optimization problem that leverages its structure. Due to the minimum-residual structure of the objective function, we employ a Levenberg-Marquardt approximation of the Hessian. We propose to use the stiffness matrix of a linear elliptic partial differential equation (PDE) with the coefficients chosen inversely proportional to the local element size as the regularization matrix. This tends to smooth out the search directions for the mesh coordinates and is particularly important for problems with elements of significantly different size. The SQP method is globalized with a line search based on the ${\ell }_1$ merit function and equipped with termination conditions based on the first-order optimality criteria. For the method to be practical for difficult problems, we initialize the solve with the $p=0$ DG solution and use continuation in the polynomial degree (solution and mesh) for high-order ($p>1$) discretizations. Finally, we identify the critical role of smoothness of the DG numerical flux function with respect to variations in the element normal in the context of implicit shock tracking, mainly with respect to solver convergence, and use smoothed versions of traditional numerical fluxes throughout.

To our knowledge, the only other approach to implicit shock tracking was proposed in [10], [11], [8], [9], [22], where the authors enforce a DG discretization with unconventional numerical fluxes and the Rankine-Hugoniot interface conditions in a minimum-residual sense. Interestingly, enforcement of the interface condition circumvented traditional stability requirements for the DG numerical fluxes, allowing them to solely rely on fluxes interior to an element. Their method was shown to successfully track even complex discontinuity surfaces and provide accurate approximations to the conservation law on traditionally coarse, high-order meshes. The present work incorporates some aspects of their method, in particular topological mesh operations and some aspects of the Hessian approximation and regularization. However, our approach that directly enforces a stable DG discretization (with zero residual) inherits many attractive features of DG methods such as guaranteed conservation and a rigorous framework for high-order convergence. Furthermore, by choosing an error-based objective function rather than a physics-based one, the extension to viscous problems only requires treatment of second-order terms in the DG setting, which has been well-established [1].

The remainder of the paper is organized as follows. Section 2 introduces the governing system of inviscid conservation laws and its discretization using a discontinuous Galerkin method. Section 3 recalls the implicit tracking framework originally proposed in [39], introduces the new error-based objective function, and discusses a parametrization of the mesh deformation that ensures the boundaries of the computational domain remain on the boundaries of the actual domain. Section 4 introduces the proposed SQP solver for the central optimization problem that incorporates a Levenberg-Marquardt approximation of the Hessian with a novel regularization matrix, a line search based on a ${\ell }_1$ merit function, and termination criteria based on the first-order optimality conditions. Section 5 discusses two important details required to make the proposed tracking framework work in practice: initialization of the SQP solver and topological mesh operations to remove small elements. Finally, Section 6 presents a number of numerical experiments that demonstrate the method is able to accurately approximate complex flows using coarse, high-order meshes and the SQP solver is able to quickly converge to a mesh that tracks all discontinuities and exhibits deep convergence to the first-order optimality conditions.