Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Efficient third order tensor-oriented directional splitting for exponential integrators

Fabio Cassini

TL;DR

The paper tackles stiff ODE systems arising from PDE discretizations with Kronecker-sum structure and develops third-order directional-split approximations for φ-functions to enable efficient tensor-based actions via the Tucker operator. It introduces three splitting strategies (two-term 2D, two-term d-D with complex coefficients, and three-term d-D with real coefficients) and provides practical implementation details for exponential Runge–Kutta time integration, leveraging small φ-function evaluations and hardware-accelerated tensor operations. Numerical experiments on the 2D Schnakenberg and 3D FitzHugh–Nagumo models show that the proposed third-order directional-split methods achieve the expected convergence while offering substantial wall-clock time advantages over non-split approaches, particularly on GPUs. The results demonstrate strong scalability and practical impact for simulating diffusion–advection–reaction systems with Kronecker-sum structure on modern hardware.

Abstract

Suitable discretizations through tensor product formulas of popular multidimensional operators (diffusion or diffusion--advection, for instance) lead to matrices with $d$-dimensional Kronecker sum structure. For evolutionary Partial Differential Equations containing such operators and integrated in time with exponential integrators, it is then of paramount importance to efficiently approximate the actions of $\varphi$-functions of the arising matrices. In this work, we show how to produce directional split approximations of third order with respect to the time step size. They conveniently employ tensor-matrix products (the so-called $μ$-mode product and related Tucker operator, realized in practice with high performance level 3 BLAS), and allow for the effective usage of exponential Runge--Kutta integrators up to order three. The technique can also be efficiently implemented on modern computer hardware such as Graphic Processing Units. The approach has been successfully tested against state-of-the-art techniques on two well-known physical models that lead to Turing patterns, namely the 2D Schnakenberg and the 3D FitzHugh--Nagumo systems, on different hardware and software architectures.

Efficient third order tensor-oriented directional splitting for exponential integrators

TL;DR

The paper tackles stiff ODE systems arising from PDE discretizations with Kronecker-sum structure and develops third-order directional-split approximations for φ-functions to enable efficient tensor-based actions via the Tucker operator. It introduces three splitting strategies (two-term 2D, two-term d-D with complex coefficients, and three-term d-D with real coefficients) and provides practical implementation details for exponential Runge–Kutta time integration, leveraging small φ-function evaluations and hardware-accelerated tensor operations. Numerical experiments on the 2D Schnakenberg and 3D FitzHugh–Nagumo models show that the proposed third-order directional-split methods achieve the expected convergence while offering substantial wall-clock time advantages over non-split approaches, particularly on GPUs. The results demonstrate strong scalability and practical impact for simulating diffusion–advection–reaction systems with Kronecker-sum structure on modern hardware.

Abstract

Suitable discretizations through tensor product formulas of popular multidimensional operators (diffusion or diffusion--advection, for instance) lead to matrices with -dimensional Kronecker sum structure. For evolutionary Partial Differential Equations containing such operators and integrated in time with exponential integrators, it is then of paramount importance to efficiently approximate the actions of -functions of the arising matrices. In this work, we show how to produce directional split approximations of third order with respect to the time step size. They conveniently employ tensor-matrix products (the so-called -mode product and related Tucker operator, realized in practice with high performance level 3 BLAS), and allow for the effective usage of exponential Runge--Kutta integrators up to order three. The technique can also be efficiently implemented on modern computer hardware such as Graphic Processing Units. The approach has been successfully tested against state-of-the-art techniques on two well-known physical models that lead to Turing patterns, namely the 2D Schnakenberg and the 3D FitzHugh--Nagumo systems, on different hardware and software architectures.
Paper Structure (11 sections, 3 theorems, 29 equations, 8 figures, 7 tables, 2 algorithms)

This paper contains 11 sections, 3 theorems, 29 equations, 8 figures, 7 tables, 2 algorithms.

Table of Contents

  1. Introduction
  2. Third order directional splitting
  3. Two-term two-dimensional splitting
  4. Two-term d-dimensional splitting with complex coefficients
  5. Three-term d-dimensional splitting with real coefficients
  6. Implementation details
  7. Numerical examples
  8. Two-dimensional Schnakenberg model
  9. Three-dimensional FitzHugh--Nagumo model
  10. Conclusions
  11. Pseudocodes

Key Result

Theorem 1

The coefficients in Table tab:phi1phi22d are the solutions of system eq:phi1phi22d for $\ell=1,2$ with $\ell_1=1$ and $\ell_2=2$.

Figures (8)

  • Figure 1: Results in MATLAB language (standard laptop) for the simulation of 2D Schnakenberg model \ref{['eq:Schnakenberg_2d']} with $n=150$ spatial discretization points per direction up to final time $T=0.25$. Top plot: error decay, with reference slope lines of order two (dashed) and three (solid). Bottom plot: work-precision diagram.
  • Figure 2: Results in C++ language (standard laptop) for the simulation of 2D Schnakenberg model \ref{['eq:Schnakenberg_2d']} with $n=150$ spatial discretization points per direction up to final time $T=0.25$. Top plot: error decay, with reference slope lines of order two (dashed) and three (solid). Bottom plot: work-precision diagram.
  • Figure 3: Turing pattern ($u$ component) for 2D Schnakenberg model \ref{['eq:Schnakenberg_2d']} obtained at final time $T=2$ with $n=150$ spatial discretization points per direction.
  • Figure 4: Results in C++ and CUDA languages (professional hardware) for the simulation of 2D Schnakenberg model \ref{['eq:Schnakenberg_2d']} with $n=300$ spatial discretization points per direction up to final time $T=0.25$. Top plot: error decay, with reference slope lines of order two (dashed) and three (solid). Bottom plot: work-precision diagram.
  • Figure 5: Results in MATLAB language (standard laptop) for the simulation of 3D FitzHugh--Nagumo model \ref{['eq:FitzHughNagumo_3d']} with $n=64$ spatial discretization points per direction up to final time $T=5$. Top plot: error decay, with reference slope lines of order two (dashed) and three (solid). Bottom plot: work-precision diagram.
  • ...and 3 more figures

Theorems & Definitions (7)

  • Theorem 1
  • proof
  • Theorem 2
  • proof
  • Theorem 3
  • proof
  • Remark 1