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

  • Two-dimensional wave propagation in layered periodic media

    Quezada de Luna, Manuel; Ketcheson, David I. (Society for Industrial & Applied Mathematics (SIAM), 2014-09-16)
    We study two-dimensional 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. One-dimensional 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 one-dimensional problems. Pseudospectral solutions of the effective medium equations agree to high accuracy with finite volume direct numerical simulations of the variable-coeffi cient equations.
  • Internal Error Propagation in Explicit Runge--Kutta Methods

    Ketcheson, David I.; Loczi, Lajos; Parsani, Matteo (Society for Industrial & Applied Mathematics (SIAM), 2014-09-11)
    In practical computation with Runge--Kutta 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 well-known 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.
  • Error Analysis of Explicit Partitioned Runge–Kutta Schemes for Conservation Laws

    Hundsdorfer, Willem; Ketcheson, David I.; Savostianov, Igor (Springer Nature, 2014-08-27)
    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 cell-based decompositions, also flux-based decompositions are studied. In the latter case mass conservation is guaranteed, but it will be seen that the accuracy may deteriorate.
  • A comparison of high-order explicit Runge–Kutta, extrapolation, and deferred correction methods in serial and parallel

    Ketcheson, David I.; Waheed, Umair bin (Mathematical Sciences Publishers, 2014-06-13)
    We compare the three main types of high-order one-step 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 fixed-order 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 eighth-order pair of Prince and Dormand (DOP8) is most efficient. But other high-order 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 N-body 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.
  • Propagation of internal errors in explicit Runge–Kutta methods and internal stability of SSP and extrapolation methods

    Ketcheson, David I.; Loczi, Lajos; Parsani, Matteo (2014-04-11)
    In practical computation with Runge--Kutta 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 well-known 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

    Ketcheson, David I.; Kocsis, Tihamer; Loczi, Lajos (Oxford University Press (OUP), 2013-12-03)
    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 half-plane is a semi-disk. We quantify this by providing disks that enclose or are enclosed by the stability region or its left half-plane part. The radius of the smallest disk centered at the origin that contains the stability region (or its portion in the left half-plane) 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 semi-disk in the left half-plane 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 semi-disks are exactly determined for 1 n 20.
  • Spatially Partitioned Embedded Runge--Kutta Methods

    Ketcheson, David I.; MacDonald, Colin B.; Ruuth, Steven J. (Society for Industrial & Applied Mathematics (SIAM), 2013-10-30)
    We study spatially partitioned embedded Runge--Kutta (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 fifth-order weighted nonoscillatory spatial discretizations. Numerical experiments are provided to support the theory.
  • Optimal Strong-Stability-Preserving Runge–Kutta Time Discretizations for Discontinuous Galerkin Methods

    Kubatko, Ethan J.; Yeager, Benjamin A.; Ketcheson, David I. (Springer Science + Business Media, 2013-10-29)
    Discontinuous Galerkin (DG) spatial discretizations are often used in a method-of-lines approach with explicit strong-stability-preserving (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 “DG-optimized” 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 fourth-order methods with up to eight stages are presented, and their stability properties are verified through application to numerical test cases.
  • Strong Stability Preserving Explicit Runge--Kutta Methods of Maximal Effective Order

    Hadjimichael, Yiannis; Macdonald, Colin B.; Ketcheson, David I.; Verner, James H. (Society for Industrial & Applied Mathematics (SIAM), 2013-07-23)
    We apply the concept of effective order to strong stability preserving (SSP) explicit Runge--Kutta methods. Relative to classical Runge--Kutta 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 four-stage SSP methods with effective order four (such methods cannot have classical order four). However, we also prove that effective order five methods---like classical order five methods---require the use of nonpositive weights and so cannot be SSP. By numerical optimization, we construct explicit SSP Runge--Kutta 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

    Quezada de Luna, Manuel; Ketcheson, David I. (Springer Nature, 2013-07-14)
    We study the behavior of nonlinear waves in a two-dimensional medium with density and stress relation that vary periodically in space. Efficient approximate Riemann solvers are developed for the corresponding variable-coefficient first-order 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 Runge--Kutta Schemes for the Spectral Difference Method Applied to Wave Propagation Problems

    Parsani, Matteo; Ketcheson, David I.; Deconinck, W. (Society for Industrial & Applied Mathematics (SIAM), 2013-04-10)
    Explicit Runge--Kutta schemes with large stable step sizes are developed for integration of high-order 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 Runge--Kutta schemes available in the literature. Furthermore, they have a small principal error norm and admit a low-storage implementation. The advantages of the new schemes are demonstrated through application to the Euler equations and the linearized Euler equations.
  • High-Order Wave Propagation Algorithms for Hyperbolic Systems

    Ketcheson, David I.; Parsani, Matteo; LeVeque, Randall J. (Society for Industrial & Applied Mathematics (SIAM), 2013-01-22)
    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 well-known 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 Runge--Kutta integration in time. The method can be extended to arbitrarily high order of accuracy and allows a well-balanced implementation for capturing solutions of balance laws near steady state. This well-balancing is achieved through the $f$-wave Riemann solver and a novel wave-slope 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 near-equilibrium solutions of balance laws.
  • Optimal stability polynomials for numerical integration of initial value problems

    Ketcheson, David I.; Ahmadia, Aron (Mathematical Sciences Publishers, 2013-01-08)
    We consider the problem of finding optimally stable polynomial approximations to the exponential for application to one-step 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

    Ketcheson, David I.; Mandli, Kyle; Ahmadia, Aron; Alghamdi, Amal; de Luna, Manuel Quezada; Parsani, Matteo; Knepley, Matthew G.; Emmett, Matthew (Society for Industrial & Applied Mathematics (SIAM), 2012-08-15)
    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 hand-coded 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 Python-based 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 high-order 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.
  • Shock dynamics in layered periodic media

    Ketcheson, David I.; Leveque, Randall J. (International Press of Boston, 2012)
    Solutions of constant-coeffcient 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 one-dimensional periodic layered medium by a computational study of time-reversibility 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 long-term 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.
  • Optimized low-order explicit Runge-Kutta schemes for high- order spectral difference method

    Parsani, Matteo; Ketcheson, David I.; Deconinck, Willem (University of Oulu, 2012)
    Optimal explicit Runge-Kutta (ERK) schemes with large stable step sizes are developed for method-of-lines 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 low-storage 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 5th-order spatial discretization.
  • Strong Stability Preserving Two-step Runge–Kutta Methods

    Ketcheson, David I.; Gottlieb, Sigal; Macdonald, Colin B. (Society for Industrial & Applied Mathematics (SIAM), 2011-12-22)
    We investigate the strong stability preserving (SSP) property of two-step Runge–Kutta (TSRK) methods. We prove that all SSP TSRK methods belong to a particularly simple subclass of TSRK methods, in which stages from the previous step are not used. We derive simple order conditions for this subclass. Whereas explicit SSP Runge–Kutta methods have order at most four, we prove that explicit SSP TSRK methods have order at most eight. We present explicit TSRK methods of up to eighth order that were found by numerical search. These methods have larger SSP coefficients than any known methods of the same order of accuracy and may be implemented in a form with relatively modest storage requirements. The usefulness of the TSRK methods is demonstrated through numerical examples, including integration of very high order weighted essentially non-oscillatory discretizations.
  • Step Sizes for Strong Stability Preservation with Downwind-Biased Operators

    Ketcheson, David I. (Society for Industrial & Applied Mathematics (SIAM), 2011-08-04)
    Strong stability preserving (SSP) integrators for initial value ODEs preserve temporal monotonicity solution properties in arbitrary norms. All existing SSP methods, including implicit methods, either require small step sizes or achieve only first order accuracy. It is possible to achieve more relaxed step size restrictions in the discretization of hyperbolic PDEs through the use of both upwind- and downwind-biased semidiscretizations. We investigate bounds on the maximum SSP step size for methods that include negative coefficients and downwind-biased semi-discretizations. We prove that the downwind SSP coefficient for linear multistep methods of order greater than one is at most equal to two, while the downwind SSP coefficient for explicit Runge–Kutta methods is at most equal to the number of stages of the method. In contrast, the maximal downwind SSP coefficient for second order Runge–Kutta methods is shown to be unbounded. We present a class of such methods with arbitrarily large SSP coefficient and demonstrate that they achieve second order accuracy for large CFL number.
  • PetClaw: A scalable parallel nonlinear wave propagation solver for Python

    Alghamdi, Amal; Ahmadia, Aron; Ketcheson, David I.; Knepley, Matthew; Mandli, Kyle; Dalcin, Lisandro (19th High Performance Computing Symposium (HPC 2011), 2011)
    We present PetClaw, a scalable distributed-memory solver for time-dependent nonlinear wave propagation. PetClaw unifies two well-known scientific computing packages, Clawpack and PETSc, using Python interfaces into both. We rely on Clawpack to provide the infrastructure and kernels for time-dependent nonlinear wave propagation. Similarly, we rely on PETSc to manage distributed data arrays and the communication between them.We describe both the implementation and performance of PetClaw as well as our challenges and accomplishments in scaling a Python-based code to tens of thousands of cores on the BlueGene/P architecture. The capabilities of PetClaw are demonstrated through application to a novel problem involving elastic waves in a heterogeneous medium. Very finely resolved simulations are used to demonstrate the suppression of shock formation in this system.
  • Using Python to Construct a Scalable Parallel Nonlinear Wave Solver

    Mandli, Kyle; Alghamdi, Amal; Ahmadia, Aron; Ketcheson, David I.; Scullin, William (THE 10th PYTHON IN SCIENCE CONF. (SCIPY 2011), 2011)
    Computational scientists seek to provide efficient, easy-to-use tools and frameworks that enable application scientists within a specific discipline to build and/or apply numerical models with up-to-date computing technologies that can be executed on all available computing systems. Although many tools could be useful for groups beyond a specific application, it is often difficult and time consuming to combine existing software, or to adapt it for a more general purpose. Python enables a high-level approach where a general framework can be supplemented with tools written for different fields and in different languages. This is particularly important when a large number of tools are necessary, as is the case for high performance scientific codes. This motivated our development of PetClaw, a scalable distributed-memory solver for time-dependent nonlinear wave propagation, as a case-study for how Python can be used as a highlevel framework leveraging a multitude of codes, efficient both in the reuse of code and programmer productivity. We present scaling results for computations on up to four racks of Shaheen, an IBM BlueGene/P supercomputer at King Abdullah University of Science and Technology. One particularly important issue that PetClaw has faced is the overhead associated with dynamic loading leading to catastrophic scaling. We use the walla library to solve the issue which does so by supplanting high-cost filesystem calls with MPI operations at a low enough level that developers may avoid any changes to their codes.

View more