Numerical Mathematics Group
Research within the Numerical Mathematics Group at KAUST focuses on the design, analysis, and implementation of numerical methods for ordinary and partial differential equations, as well as application of numerical methods to problems in nonlinear wave propagation. Their work spans the numerical method development pipeline from theoretical numerical analysis, to computational experimentation, to software development and high performance computational science, in collaboration with application domain researchers. Particular emphasis is placed on robustness, wide applicability, and efficiency.
Recent Submissions

Relaxation RungeKutta Methods: Fully Discrete Explicit EntropyStable Schemes for the Compressible Euler and NavierStokes Equations(SIAM Journal on Scientific Computing, Society for Industrial & Applied Mathematics (SIAM), 20200312) [Article]The framework of inner product norm preserving relaxation RungeKutta methods [D. I. Ketcheson, SIAM J. Numer. Anal., 57 (2019), pp. 28502870] is extended to general convex quantities. Conservation, dissipation, or other solution properties with respect to any convex functional are enforced by the addition of a relaxation parameter that multiplies the RungeKutta update at each step. Moreover, other desirable stability (such as strong stability preservation) and efficiency (such as low storage requirements) properties are preserved. The technique can be applied to both explicit and implicit RungeKutta methods and requires only a small modification to existing implementations. The computational cost at each step is the solution of one additional scalar algebraic equation for which a good initial guess is available. The effectiveness of this approach is proved analytically and demonstrated in several numerical examples, including applications to high order entropyconservative and entropystable semidiscretizations on unstructured grids for the compressible Euler and NavierStokes equations.

Relaxation RungeKutta Methods: FullyDiscrete Explicit EntropyStable Schemes for the Euler and NavierStokes Equations(arXiv, 20190522) [Preprint]The framework of inner product norm preserving relaxation RungeKutta methods (David I. Ketcheson, Relaxation RungeKutta Methods: Conservation and Stability for InnerProduct Norms, 2019. arXiv: 1905.09847 [math.NA]) is extended to general convex quantities. Conservation, dissipation, or other solution properties with respect to any convex functional are enforced by the addition of a relaxation parameter that multiplies the RungeKutta update at each step. Moreover, other desirable stability (such as strong stability preservation) and efficiency (such as low storage requirements) properties are preserved. The technique can be applied to both explicit and implicit RungeKutta methods and requires only a small modification to existing implementations. The computational cost at each step is the solution of one additional scalar algebraic equation for which a good initial guess is available. The effectiveness of this approach is proved analytically and demonstrated in several numerical examples, including applications to highorder entropyconservative and entropystable semidiscretizations on unstructured grids for the compressible Euler and NavierStokes equations.

Diffractons: Solitary Waves Created by Diffraction in Periodic Media(Multiscale Modeling & Simulation, Society for Industrial & Applied Mathematics (SIAM), 20150331) [Article]A new class of solitary waves arises in the solution of nonlinear wave equations with constant impedance and no dispersive terms. These solitary waves depend on a balance between nonlinearity and a dispersionlike effect due to spatial variation in the sound speed of the medium. A highorder homogenized model confirms this effective dispersive behavior, and its solutions agree well with those obtained by direct simulation of the variablecoefficient system. These waves are observed to be longtime stable, globally attracting solutions that arise in general as solutions to nonlinear wave problems with periodically varying sound speed. They share some properties with known classes of solitary waves but possess important differences as well.

Twodimensional wave propagation in layered periodic media(SIAM Journal on Applied Mathematics, Society for Industrial & Applied Mathematics (SIAM), 20141203) [Article]We study twodimensional wave propagation in materials whose properties vary periodically in one direction only. High order homogenization is carried out to derive a dispersive effective medium approximation. Onedimensional materials with constant impedance exhibit no effective dispersion. We show that a new kind of effective dispersion may arise in two dimensions, even in materials with constant impedance. This dispersion is a macroscopic effect of microscopic diffraction caused by spatial variation in the sound speed. We analyze this dispersive effect by using highorder homogenization to derive an anisotropic, dispersive effective medium. We generalize to two dimensions a homogenization approach that has been used previously for onedimensional problems. Pseudospectral solutions of the effective medium equations agree to high accuracy with finite volume direct numerical simulations of the variablecoeffi cient equations.

Internal Error Propagation in Explicit RungeKutta Methods(SIAM Journal on Numerical Analysis, Society for Industrial & Applied Mathematics (SIAM), 20140917) [Article]In practical computation with RungeKutta methods, the stage equations are not satisfied exactly, due to roundoff errors, algebraic solver errors, and so forth. We show by example that propagation of such errors within a single step can have catastrophic effects for otherwise practical and wellknown methods. We perform a general analysis of internal error propagation, emphasizing that it depends significantly on how the method is implemented. We show that for a fixed method, essentially any set of internal stability polynomials can be obtained by modifying the implementation details. We provide bounds on the internal error amplification constants for some classes of methods with many stages, including strong stability preserving methods and extrapolation methods. These results are used to prove error bounds in the presence of roundoff or other internal errors.

On the absolute stability regions corresponding to partial sums of the exponential function(IMA Journal of Numerical Analysis, Oxford University Press (OUP), 20140915) [Article]Certain numerical methods for initial value problems have as stability function the nth partial sum of the exponential function. We study the stability region, i.e., the set in the complex plane over which the nth partial sum has at most unit modulus. It is known that the asymptotic shape of the part of the stability region in the left halfplane is a semidisk. We quantify this by providing disks that enclose or are enclosed by the stability region or its left halfplane part. The radius of the smallest disk centered at the origin that contains the stability region (or its portion in the left halfplane) is determined for 1 n 20. Bounds on such radii are proved for n 2; these bounds are shown to be optimal in the limit n ! +1. We prove that the stability region and its complement, restricted to the imaginary axis, consist of alternating intervals of length tending to , as n ! 1. Finally, we prove that a semidisk in the left halfplane with vertical boundary being the imaginary axis and centered at the origin is included in the stability region if and only if n 0 mod 4 or n 3 mod 4. The maximal radii of such semidisks are exactly determined for 1 n 20.

Error Analysis of Explicit Partitioned Runge–Kutta Schemes for Conservation Laws(Journal of Scientific Computing, Springer Nature, 20140827) [Article]An error analysis is presented for explicit partitioned Runge–Kutta methods and multirate methods applied to conservation laws. The interfaces, across which different methods or time steps are used, lead to order reduction of the schemes. Along with cellbased decompositions, also fluxbased decompositions are studied. In the latter case mass conservation is guaranteed, but it will be seen that the accuracy may deteriorate.

A comparison of highorder explicit Runge–Kutta, extrapolation, and deferred correction methods in serial and parallel(Communications in Applied Mathematics and Computational Science, Mathematical Sciences Publishers, 20140613) [Article]We compare the three main types of highorder onestep initial value solvers: extrapolation, spectral deferred correction, and embedded Runge–Kutta pairs. We consider orders four through twelve, including both serial and parallel implementations. We cast extrapolation and deferred correction methods as fixedorder Runge–Kutta methods, providing a natural framework for the comparison. The stability and accuracy properties of the methods are analyzed by theoretical measures, and these are compared with the results of numerical tests. In serial, the eighthorder pair of Prince and Dormand (DOP8) is most efficient. But other highorder methods can be more efficient than DOP8 when implemented in parallel. This is demonstrated by comparing a parallelized version of the wellknown ODEX code with the (serial) DOP853 code. For an Nbody problem with N = 400, the experimental extrapolation code is as fast as the tuned Runge–Kutta pair at loose tolerances, and is up to two times as fast at tight tolerances.

Rational functions with maximal radius of absolute monotonicity(LMS Journal of Computation and Mathematics, Oxford University Press (OUP), 20140519) [Article]We study the radius of absolute monotonicity R of rational functions with numerator and denominator of degree s that approximate the exponential function to order p. Such functions arise in the application of implicit sstage, order p RungeKutta methods for initial value problems and the radius of absolute monotonicity governs the numerical preservation of properties like positivity and maximumnorm contractivity. We construct a function with p=2 and R>2s, disproving a conjecture of van de Griend and Kraaijevanger. We determine the maximum attainable radius for functions in several oneparameter families of rational functions. Moreover, we prove earlier conjectured optimal radii in some families with 2 or 3 parameters via uniqueness arguments for systems of polynomial inequalities. Our results also prove the optimality of some strong stability preserving implicit and singly diagonally implicit RungeKutta methods. Whereas previous results in this area were primarily numerical, we give all constants as exact algebraic numbers.

Propagation of internal errors in explicit Runge–Kutta methods and internal stability of SSP and extrapolation methods(20140411) [Technical Report]In practical computation with RungeKutta methods, the stage equations are not satisfied exactly, due to roundoff errors, algebraic solver errors, and so forth. We show by example that propagation of such errors within a single step can have catastrophic effects for otherwise practical and wellknown methods. We perform a general analysis of internal error propagation, emphasizing that it depends significantly on how the method is implemented. We show that for a fixed method, essentially any set of internal stability polynomials can be obtained by modifying the implementation details. We provide bounds on the internal error amplification constants for some classes of methods with many stages, including strong stability preserving methods and extrapolation methods. These results are used to prove error bounds in the presence of roundoff or other internal errors.

Spatially Partitioned Embedded RungeKutta Methods(SIAM Journal on Numerical Analysis, Society for Industrial & Applied Mathematics (SIAM), 20131030) [Article]We study spatially partitioned embedded RungeKutta (SPERK) schemes for partial differential equations (PDEs), in which each of the component schemes is applied over a different part of the spatial domain. Such methods may be convenient for problems in which the smoothness of the solution or the magnitudes of the PDE coefficients vary strongly in space. We focus on embedded partitioned methods as they offer greater efficiency and avoid the order reduction that may occur in nonembedded schemes. We demonstrate that the lack of conservation in partitioned schemes can lead to nonphysical effects and propose conservative additive schemes based on partitioning the fluxes rather than the ordinary differential equations. A variety of SPERK schemes are presented, including an embedded pair suitable for the time evolution of fifthorder weighted nonoscillatory spatial discretizations. Numerical experiments are provided to support the theory.

Optimal StrongStabilityPreserving Runge–Kutta Time Discretizations for Discontinuous Galerkin Methods(Journal of Scientific Computing, Springer Nature, 20131029) [Article]Discontinuous Galerkin (DG) spatial discretizations are often used in a methodoflines approach with explicit strongstabilitypreserving (SSP) Runge–Kutta (RK) time steppers for the numerical solution of hyperbolic conservation laws. The time steps that are employed in this type of approach must satisfy Courant–Friedrichs–Lewy stability constraints that are dependent on both the region of absolute stability and the SSP coefficient of the RK method. While existing SSPRK methods have been optimized with respect to the latter, it is in fact the former that gives rise to stricter constraints on the time step in the case of RKDG stability. Therefore, in this work, we present the development of new “DGoptimized” SSPRK methods with stability regions that have been specifically designed to maximize the stable time step size for RKDG methods of a given order in one space dimension. These new methods represent the best available RKDG methods in terms of computational efficiency, with significant improvements over methods using existing SSPRK time steppers that have been optimized with respect to SSP coefficients. Second, third, and fourthorder methods with up to eight stages are presented, and their stability properties are verified through application to numerical test cases.

Strong Stability Preserving Explicit RungeKutta Methods of Maximal Effective Order(SIAM Journal on Numerical Analysis, Society for Industrial & Applied Mathematics (SIAM), 20130723) [Article]We apply the concept of effective order to strong stability preserving (SSP) explicit RungeKutta methods. Relative to classical RungeKutta methods, methods with an effective order of accuracy are designed to satisfy a relaxed set of order conditions but yield higher order accuracy when composed with special starting and stopping methods. We show that this allows the construction of fourstage SSP methods with effective order four (such methods cannot have classical order four). However, we also prove that effective order five methodslike classical order five methodsrequire the use of nonpositive weights and so cannot be SSP. By numerical optimization, we construct explicit SSP RungeKutta methods up to effective order four and establish the optimality of many of them. Numerical experiments demonstrate the validity of these methods in practice.

Numerical Simulation of Cylindrical Solitary Waves in Periodic Media(Journal of Scientific Computing, Springer Nature, 20130714) [Article]We study the behavior of nonlinear waves in a twodimensional medium with density and stress relation that vary periodically in space. Efficient approximate Riemann solvers are developed for the corresponding variablecoefficient firstorder hyperbolic system. We present direct numerical simulations of this multiscale problem, focused on the propagation of a single localized perturbation in media with strongly varying impedance. For the conditions studied, we find little evidence of shock formation. Instead, solutions consist primarily of solitary waves. These solitary waves are observed to be stable over long times and to interact in a manner approximately like solitons. The system considered has no dispersive terms; these solitary waves arise due to the material heterogeneity, which leads to strong reflections and effective dispersion.

Optimized Explicit RungeKutta Schemes for the Spectral Difference Method Applied to Wave Propagation Problems(SIAM Journal on Scientific Computing, Society for Industrial & Applied Mathematics (SIAM), 20130410) [Article]Explicit RungeKutta schemes with large stable step sizes are developed for integration of highorder spectral difference spatial discretizations on quadrilateral grids. The new schemes permit an effective time step that is substantially larger than the maximum admissible time step of standard explicit RungeKutta schemes available in the literature. Furthermore, they have a small principal error norm and admit a lowstorage implementation. The advantages of the new schemes are demonstrated through application to the Euler equations and the linearized Euler equations.

HighOrder Wave Propagation Algorithms for Hyperbolic Systems(SIAM Journal on Scientific Computing, Society for Industrial & Applied Mathematics (SIAM), 20130122) [Article]We present a finite volume method that is applicable to hyperbolic PDEs including spatially varying and semilinear nonconservative systems. The spatial discretization, like that of the wellknown Clawpack software, is based on solving Riemann problems and calculating fluctuations (not fluxes). The implementation employs weighted essentially nonoscillatory reconstruction in space and strong stability preserving RungeKutta integration in time. The method can be extended to arbitrarily high order of accuracy and allows a wellbalanced implementation for capturing solutions of balance laws near steady state. This wellbalancing is achieved through the $f$wave Riemann solver and a novel waveslope WENO reconstruction procedure. The wide applicability and advantageous properties of the method are demonstrated through numerical examples, including problems in nonconservative form, problems with spatially varying fluxes, and problems involving nearequilibrium solutions of balance laws.

Optimal stability polynomials for numerical integration of initial value problems(Communications in Applied Mathematics and Computational Science, Mathematical Sciences Publishers, 20121231) [Article]We consider the problem of finding optimally stable polynomial approximations to the exponential for application to onestep integration of initial value ordinary and partial differential equations. The objective is to find the largest stable step size and corresponding method for a given problem when the spectrum of the initial value problem is known. The problem is expressed in terms of a general least deviation feasibility problem. Its solution is obtained by a new fast, accurate, and robust algorithm based on convex optimization techniques. Global convergence of the algorithm is proven in the case that the order of approximation is one and in the case that the spectrum encloses a starlike region. Examples demonstrate the effectiveness of the proposed algorithm even when these conditions are not satisfied.

PyClaw: Accessible, Extensible, Scalable Tools for Wave Propagation Problems(SIAM Journal on Scientific Computing, Society for Industrial & Applied Mathematics (SIAM), 20120815) [Article]Development of scientific software involves tradeoffs between ease of use, generality, and performance. We describe the design of a general hyperbolic PDE solver that can be operated with the convenience of MATLAB yet achieves efficiency near that of handcoded Fortran and scales to the largest supercomputers. This is achieved by using Python for most of the code while employing automatically wrapped Fortran kernels for computationally intensive routines, and using Python bindings to interface with a parallel computing library and other numerical packages. The software described here is PyClaw, a Pythonbased structured grid solver for general systems of hyperbolic PDEs [K. T. Mandli et al., PyClaw Software, Version 1.0, http://numerics.kaust.edu.sa/pyclaw/ (2011)]. PyClaw provides a powerful and intuitive interface to the algorithms of the existing Fortran codes Clawpack and SharpClaw, simplifying code development and use while providing massive parallelism and scalable solvers via the PETSc library. The package is further augmented by use of PyWENO for generation of efficient highorder weighted essentially nonoscillatory reconstruction code. The simplicity, capability, and performance of this approach are demonstrated through application to example problems in shallow water flow, compressible flow, and elasticity.

Optimized loworder explicit RungeKutta schemes for high order spectral difference method(Proceedings of the 11th Finnish Mechanics Days, University of Oulu, 2012) [Conference Paper]Optimal explicit RungeKutta (ERK) schemes with large stable step sizes are developed for methodoflines discretizations based on the spectral difference (SD) spatial discretization on quadrilateral grids. These methods involve many stages and provide the optimal linearly stable time step for a prescribed SD spectrum and the minimum leading truncation error coefficient, while admitting a lowstorage implementation. Using a large number of stages, the new ERK schemes lead to efficiency improvements larger than 60% over standard ERK schemes for 4th and 5thorder spatial discretization.

Shock dynamics in layered periodic media(Communications in Mathematical Sciences, International Press of Boston, 2012) [Article]Solutions of constantcoeffcient nonlinear hyperbolic PDEs generically develop shocks, even if the initial data is smooth. Solutions of hyperbolic PDEs with variable coeffcients can behave very differently. We investigate formation and stability of shock waves in a onedimensional periodic layered medium by a computational study of timereversibility and entropy evolution. We find that periodic layered media tend to inhibit shock formation. For small initial conditions and large impedance variation, no shock formation is detected even after times much greater than the time of shock formation in a homogeneous medium. Furthermore, weak shocks are observed to be dynamically unstable in the sense that they do not lead to significant longterm entropy decay. We propose a characteristic condition for admissibility of shocks in heterogeneous media that generalizes the classical Lax entropy condition and accurately predicts the formation or absence of shocks in these media.