Sylvester-Preconditioned Adaptive-Rank Implicit Time Integrators for Advection-Diffusion Equations with Variable Coefficients

TL;DR

The paper addresses efficient time integration of multi-dimensional, variable-coefficient advection-diffusion PDEs by reformulating the implicit discretization as a generalized Sylvester equation (GSE) and solving it with a projection-based, adaptive-rank Krylov method. It introduces an averaged-coefficient Sylvester (ACS) preconditioner and constructs dimension-wise xKrylov subspaces to capture low-rank solution structure, using SVD/HOSVD-based truncation to control growth. The main theoretical contributions are explicit complexity scalings: $\mathcal{O}(N r^2) + \mathcal{O}(r^{d+1})$ operations and $\mathcal{O}(N r) + \mathcal{O}(r^d)$ memory (with $d=2,3$), validated through extensive 2D/3D numerical experiments showing rank largely independent of 1D grid size $N$ and mesh-independent GMRES convergence under ACS. The practical impact is the ability to perform large-scale 3D advection-diffusion simulations on standard hardware with near-constant rank, mitigating the curse of dimensionality for low-rank solutions.

Abstract

We consider the adaptive-rank integration of {2D and 3D} time-dependent advection-diffusion partial differential equations (PDEs) with variable coefficients. We employ a standard finite-difference method for spatial discretization coupled with diagonally implicit Runge-Kutta temporal schemes. The discrete equation is a generalized Sylvester equation (GSE), which we solve with an adaptive-rank algorithm structured around three key strategies: {(i) constructing dimension-wise subspaces based on an extended Krylov strategy, (ii) developing an effective preconditioner for the reduced coefficient matrix, and (iii) efficiently computing the residual of the equation without explicitly reverting to the full-rank form. {The low-rank decomposition is performed in 2D using SVD, and with high-order SVD (HOSVD) in 3D to represent the tensor in a compressed Tucker format.} The computational complexity of the proposed approach {is demonstrated numerically to} be comparable to the constant-coefficient case [El Kahza et al, J. Comput. Phys., 518 (2024)], scaling as $\mathcal{O}(N {r^2} + {r^{d+1}})$ for $d$-dimensional problems (here, $d = 2$ or $3$), with $N$ the resolution in one dimension and $r$ the maximal rank during the Krylov iteration (which we find to be largely independent of $N$). We present numerical examples that illustrate the computational efficacy and complexity of our algorithm.}

Figures (8)

• Figure 1: Example 4.1. (a) Temporal convergence study of BE, DIRK2, and DIRK3 schemes. (b) Rank evolution of the solution in the $x$, $y$, and $z$ directions over time for $\lambda_D = 140$.
• Figure 2: Example 4.1. Convergence study using the xKrylov subspace method with DIRK3 adaptive-rank integrator for several outer basis tolerances, $\epsilon_{tol}$.
• Figure 3: Example 4.1. Log-log plots of GMRES residuals versus iterations, with and without preconditioning. Simulations are conducted for a single time step with $\Delta t = 0.01$ at the third outer-basis augmentation, using grid resolutions $N_x = N_y = N_z = 200$ in (a), $1000$ in (b), and $2500$ in (c).
• Figure 4: Example 4.1. (a) Average condition number for Galerkin-projected advection and diffusion operators. (b)-(d) Log-log plots of GMRES residuals versus iterations, with and without preconditioning for various subspace sizes $r_\kappa = 37$, $r_\kappa = 73$, and $r_\kappa = 217$, respectively. These cases correspond to the highlighted blue data points in (a). Simulations are conducted for a single time step with $\Delta t = 0.01$ and a grid of $N_x = N_y = N_z = 1000$.
• Figure 5: Example 4.1. (a) Profiling study depicting the computational cost of various components of the algorithm, demonstrating overall${\cal O}(N)$ computational complexity of DIRK3 adaptive-rank integrator in 3D. (b) Demonstration of ${\cal O}(r^{4})$ computational complexity of the GMRES-ACS solver in 3D. For this test only, we allow the rank to grow indefinitely by setting the outer residual tolerance, $\epsilon_{\texttt{tol}}$, to machine precision. (c) Demonstration of ${\cal O}(N)$ memory storage complexity of DIRK3 adaptive-rank integrator in 3D. The execution times were recorded using MATLAB's timeit function. The memory usage was recorded using MATLAB’s whos command at the end of each time step, and the values reported correspond to the average memory usage across all time steps in megabytes.