Per-Olof Persson

Matthew J. Zahr, Per-Olof Persson

This work introduces a novel discontinuity-tracking framework for resolving discontinuous solutions of conservation laws with high-order numerical discretizations that support inter-element solution discontinuities, such as discontinuous Galerkin methods. The proposed method aims to align inter-element boundaries with discontinuities in the solution by deforming the computational mesh. A discontinuity-aligned mesh ensures the discontinuity is represented through inter-element jumps while smooth basis functions interior to elements are only used to approximate smooth regions of the solution, thereby avoiding Gibbs' phenomena that create well-known stability issues. Therefore, very coarse high-order discretizations accurately resolve the piecewise smooth solution throughout the domain, provided the discontinuity is tracked. Central to the proposed discontinuity-tracking framework is a discrete PDE-constrained optimization formulation that simultaneously aligns the computational mesh with discontinuities in the solution and solves the discretized conservation law on this mesh. The optimization objective is taken as a combination of the the deviation of the finite-dimensional solution from its element-wise average and a mesh distortion metric to simultaneously penalize Gibbs' phenomena and distorted meshes. We advocate a gradient-based, full space solver where the mesh and conservation law solution converge to their optimal values simultaneously and therefore never require the solution of the discrete conservation law on a non-aligned mesh. The merit of the proposed method is demonstrated on a number of one- and two-dimensional model problems including 2D supersonic flow around a bluff body. We demonstrate optimal $\mathcal{O}(h^{p+1})$ convergence rates in the $L^1$ norm for up to polynomial order $p=6$ and show that accurate solutions can be obtained on extremely coarse meshes.

Will Pazner, Per-Olof Persson

In this paper, we develop new techniques for solving the large, coupled linear systems that arise from fully implicit Runge-Kutta methods. This method makes use of the iterative preconditioned GMRES algorithm for solving the linear systems, which has seen success for fluid flow problems and discontinuous Galerkin discretizations. By transforming the resulting linear system of equations, one can obtain a method which is much less computationally expensive than the untransformed formulation, and which compares competitively with other time-integration schemes, such as diagonally implicit Runge-Kutta (DIRK) methods. We develop and test several ILU-based preconditioners effective for these large systems. We additionally employ a parallel-in-time strategy to compute the Runge-Kutta stages simultaneously. Numerical experiments are performed on the Navier-Stokes equations using Euler vortex and 2D and 3D NACA airfoil test cases in serial and in parallel settings. The fully implicit Radau IIA Runge-Kutta methods compare favorably with equal-order DIRK methods in terms of accuracy, number of GMRES iterations, number of matrix-vector multiplications, and wall-clock time, for a wide range of time steps.

Will Pazner, Per-Olof Persson

In this paper, we develop a new tensor-product based preconditioner for discontinuous Galerkin methods with polynomial degrees higher than those typically employed. This preconditioner uses an automatic, purely algebraic method to approximate the exact block Jacobi preconditioner by Kronecker products of several small, one-dimensional matrices. Traditional matrix-based preconditioners require $\mathcal{O}(p^{2d})$ storage and $\mathcal{O}(p^{3d})$ computational work, where $p$ is the degree of basis polynomials used, and $d$ is the spatial dimension. Our SVD-based tensor-product preconditioner requires $\mathcal{O}(p^{d+1})$ storage, $\mathcal{O}(p^{d+1})$ work in two spatial dimensions, and $\mathcal{O}(p^{d+2})$ work in three spatial dimensions. Combined with a matrix-free Newton-Krylov solver, these preconditioners allow for the solution of DG systems in linear time in $p$ per degree of freedom in 2D, and reduce the computational complexity from $\mathcal{O}(p^9)$ to $\mathcal{O}(p^5)$ in 3D. Numerical results are shown in 2D and 3D for the advection and Euler equations, using polynomials of degree up to $p=15$. For many test cases, the preconditioner results in similar iteration counts when compared with the exact block Jacobi preconditioner, and performance is significantly improved for high polynomial degrees $p$.

Per-Olof Persson

We present a new line-based discontinuous Galerkin (DG) discretization scheme for first- and second-order systems of partial differential equations. The scheme is based on fully unstructured meshes of quadrilateral or hexahedral elements, and it is closely related to the standard nodal DG scheme as well as several of its variants such as the collocation-based DG spectral element method (DGSEM) or the spectral difference (SD) method. However, our motivation is to maximize the sparsity of the Jacobian matrices, since this directly translates into higher performance in particular for implicit solvers, while maintaining many of the good properties of the DG scheme. To achieve this, our scheme is based on applying one-dimensional DG solvers along each coordinate direction in a reference element. This reduces the number of connectivities drastically, since the scheme only connects each node to a line of nodes along each direction, as opposed to the standard DG method which connects all nodes inside the element and many nodes in the neighboring ones. The resulting scheme is similar to a collocation scheme, but it uses fully consistent integration along each 1-D coordinate direction which results in different properties for nonlinear problems and curved elements. Also, the scheme uses solution points along each element face, which further reduces the number of connections with the neighboring elements. Second-order terms are handled by an LDG-type approach, with an upwind/downwind flux function based on a switch function at each element face. We demonstrate the accuracy of the method and compare it to the standard nodal DG method for problems including Poisson's equation, Euler's equations of gas dynamics, and both the steady-state and the transient compressible Navier-Stokes equations.

Daniel Z. Huang, Per-Olof Persson, Matthew J. Zahr

This work introduces a general framework for constructing high-order, linearly stable, partitioned solvers for multiphysics problems from a monolithic implicit-explicit Runge-Kutta (IMEX-RK) discretization of the semi-discrete equations. The generic multiphysics problem is modeled as a system of n systems of partial differential equations where the ith subsystem is coupled to the other subsystems through a coupling term that can depend on the state of all the other subsystems. This coupled system of partial differential equations reduces to a coupled system of ordinary differential equations via the method of lines where an appropriate spatial discretization is applied to each subsystem. The coupled system of ordinary differential equations is taken as a monolithic system and discretized using an IMEX-RK discretization with a specific implicit-explicit decomposition that introduces the concept of a predictor for the coupling term. We propose four coupling predictors that enable the monolithic system to be solved in a partitioned manner and preserve the IMEX-RK structure and therefore the design order of accuracy of the monolithic scheme. The four partitioned solvers that result from these predictors are high-order accurate, allow for maximum re-use of existing single-physics software, and two of the four solvers allow the subsystems to be solved in parallel at a given stage and time step. We also analyze the stability of a coupled, linear model problem and show that one of the partitioned solvers achieves unconditional linear stability, while the others are unconditionally stable only for certain values of the coupling strength. We demonstrate the performance of the proposed partitioned solvers on several classes of multiphysics problems including a simple linear system of ODEs, advection-diffusion-reaction systems, FSI problems, and particle-laden flows.

Will Pazner, Per-Olof Persson

We study the convergence of iterative linear solvers for discontinuous Galerkin discretizations of systems of hyperbolic conservation laws with polygonal mesh elements compared with that of traditional triangular elements. We solve the semi-discrete system of equations by means of an implicit time discretization method, using iterative solvers such as the block Jacobi method and GMRES. We perform a von Neumann analysis to analytically study the convergence of the block Jacobi method for the two-dimensional advection equation on four classes of regular meshes: hexagonal, square, equilateral-triangular, and right-triangular. We find that hexagonal and square meshes give rise to smaller eigenvalues, and thus result in faster convergence of Jacobi's method. We perform numerical experiments with variable velocity fields, irregular, unstructured meshes, and the Euler equations of gas dynamics to confirm and extend these results. We additionally study the effect of polygonal meshes on the performance of block ILU(0) and Jacobi preconditioners for the GMRES method.

Daniel Z. Huang, Per-Olof Persson, Matthew J. Zahr

A high-order accurate adjoint-based optimization framework is presented for unsteady multiphysics problems. The fully discrete adjoint solver relies on the high-order, linearly stable, partitioned solver introduced in [1], where different subsystems are modeled and discretized separately. The coupled system of semi-discretized ordinary differential equations is taken as a monolithic system and partitioned using an implicit-explicit Runge-Kutta (IMEX-RK) discretization [2]. Quantities of interest (QoI) that take the form of space-time integrals are discretized in a solver-consistent manner. The corresponding adjoint equations are derived to compute exact gradients of QoI, which can be solved in a partitioned manner, i.e. subsystem-by-subsystem and substage-by-substage, thanks to the partitioned primal solver. These quantities of interest and their gradients are then used in the context of gradient-based PDE-constrained optimization. The present optimization framework is applied to two fluid-structure interaction problems: 1D piston problem with a three-field formulation and a 2D energy harvesting problem with a two-field formulation.

Arjun Narayanan, Yulong Pan, Per-Olof Persson

We present a learning based framework for mesh quality improvement on unstructured triangular and quadrilateral meshes. Our model learns to improve mesh quality according to a prescribed objective function purely via self-play reinforcement learning with no prior heuristics. The actions performed on the mesh are standard local and global element operations. The goal is to minimize the deviation of the node degrees from their ideal values, which in the case of interior vertices leads to a minimization of irregular nodes.

Will Pazner, Per-Olof Persson

We develop a discretely entropy-stable line-based discontinuous Galerkin method for hyperbolic conservation laws based on a flux differencing technique. By using standard entropy-stable and entropy-conservative numerical flux functions, this method guarantees that the discrete integral of the entropy is non-increasing. This nonlinear entropy stability property is important for the robustness of the method, in particular when applied to problems with discontinuous solutions or when the mesh is under-resolved. This line-based method is significantly less computationally expensive than a standard DG method. Numerical results are shown demonstrating the effectiveness of the method on a variety of test cases, including Burgers' equation and the Euler equations, in one, two, and three spatial dimensions.

Jaume Peraire, Per-Olof Persson

We present a compact discontinuous Galerkin (CDG) method for an elliptic model problem. The problem is first cast as a system of first order equations by introducing the gradient of the primal unknown, or flux, as an additional variable. A standard discontinuous Galerkin (DG) method is then applied to the resulting system of equations. The numerical interelement fluxes are such that the equations for the additional variable can be eliminated at the element level, thus resulting in a global system that involves only the original unknown variable. The proposed method is closely related to the local discontinuous Galerkin (LDG) method [B. Cockburn and C.-W. Shu, SIAM J. Numer. Anal., 35 (1998), pp. 2440-2463], but, unlike the LDG method, the sparsity pattern of the CDG method involves only nearest neighbors. Also, unlike the LDG method, the CDG method works without stabilization for an arbitrary orientation of the element interfaces. The computation of the numerical interface fluxes for the CDG method is slightly more involved than for the LDG method, but this additional complication is clearly offset by increased compactness and flexibility. Compared to the BR2 [F. Bassi and S. Rebay, J. Comput. Phys., 131 (1997), pp. 267-279] and IP [J. Douglas, Jr., and T. Dupont, in Computing Methods in Applied Sciences (Second Internat. Sympos., Versailles, 1975), Lecture Notes in Phys. 58, Springer, Berlin, 1976, pp. 207-216] methods, which are known to be compact, the present method produces fewer nonzero elements in the matrix and is computationally more efficient.

Alan Edelman, Per-Olof Persson

We present efficient numerical techniques for calculation of eigenvalue distributions of random matrices in the beta-ensembles. We compute histograms using direct simulations on very large matrices, by using tridiagonal matrices with appropriate simplifications. The distributions are also obtained by numerical solution of the Painleve II and V equations with high accuracy. For the spacings we show a technique based on the Prolate matrix and Richardson extrapolation, and we compare the distributions with the zeros of the Riemann zeta function.