Krylov-based Adaptive-Rank Implicit Time Integrators for Stiff Problems with Application to Nonlinear Fokker-Planck Kinetic Models

TL;DR

The paper addresses solving high-dimensional, stiff time-dependent PDEs with strict conservation using a high-order implicit method that adaptively discovers low-rank representations. It introduces an extended Krylov subspace–based adaptive-rank DIRK framework, augmented with a residual/LTE stopping rule and a LoMaC projection to preserve mass, momentum, and energy, and applies it to the Lenard-Bernstein Fokker-Planck model as well as a prototype heat equation. Key contributions include a detailed complexity analysis showing linear-in-dimension scaling, a Krylov-based reduced Sylvester solver with adaptive rank growth, and a LoMaC-based macroscopic conservation mechanism, achieving third-order temporal accuracy in tests. The results demonstrate accuracy comparable to full-rank solvers with substantial speedups and robust conservation, suggesting a scalable pathway for high-dimensional kinetic problems and other stiff time-dependent PDEs.

Abstract

We propose a high order adaptive-rank implicit integrators for stiff time-dependent PDEs, leveraging extended Krylov subspaces to efficiently and adaptively populate low-rank solution bases. This allows for the accurate representation of solutions with significantly reduced computational costs. We further introduce an efficient mechanism for residual evaluation and an adaptive rank-seeking strategy that optimizes low-rank settings based on a comparison between the residual size and the local truncation errors of the time-stepping discretization. We demonstrate our approach with the challenging Lenard-Bernstein Fokker-Planck (LBFP) nonlinear equation, which describes collisional processes in a fully ionized plasma. The preservation of {the equilibrium state} is achieved through the Chang-Cooper discretization, and strict conservation of mass, momentum and energy via a Locally Macroscopic Conservative (LoMaC) procedure. The development of implicit adaptive-rank integrators, demonstrated here up to third-order temporal accuracy via diagonally implicit Runge-Kutta schemes, showcases superior performance in terms of accuracy, computational efficiency, equilibrium preservation, and conservation of macroscopic moments. This study offers a starting point for developing scalable, efficient, and accurate methods for high-dimensional time-dependent problems.

Figures (8)

• Figure 1: Flow-chart of the extended-Krylov-based implicit adaptive-rank algorithm with the LoMaC projection.
• Figure 2: Simulation of the heat equation under periodic boundary conditions employing the BE, DIRK2, and DIRK3 integrators. In Fig. (a), we illustrate a temporal error convergence study spanning a range of $\Delta t/\Delta x^2$ values [100:900], presenting results for both the adaptive-rank integrator and the classical full-rank integrator counterpart. Fig. (b) presents results for both the BE adaptive-rank integrator with and without LoMaC projection, and the classical full-rank integrator counterpart. Fig. (c) shows the rank evolution as a function of time for the first-order adaptive-rank BE integrator with LoMaC projection.
• Figure 3: Simulation of the Fokker-Planck equation employing three temporal integrators--specifically Backward Euler, DIRK2, and DIRK3. In Fig. (a), we illustrate a temporal error convergence study spanning a range $\lambda=\frac{\Delta t}{{\Delta v}^2}$ values [100:600] for both the low-rank integrator and the classical full-rank counterpart. Fig. (b) shows the rank evolution as a function of time for the Backward Euler simulation.
• Figure 4: Comparison of computational complexity between the adaptive-rank integrator (illustrated in blue) and the full-rank integrator (shown in red) for BE, DIRK2, and DIRK3. Here, $N$ represents the number of grid points for each velocity dimension, denoted as $N_{v_1}$ and $N_{v_2}$, where for simplicity we set $N=N_{v_1}=N_{v_2}$. The simulation time was measured using MATLAB's timeit function.
• Figure 5: Evolution of the extended Krylov iterations per timestep over time for various $\lambda$ numbers for BE. This iteration number is representative for each stage in DIRK.

Theorems & Definitions (2)

• Proposition 3.1
• proof