🤖 AI Summary
This work addresses the high memory and computational complexity typically associated with solving three-dimensional partial differential equations on Cartesian grids. By exploiting tensor-product structure, the proposed method decomposes the 3D operator into one-dimensional banded kernels aligned with coordinate axes, thereby avoiding explicit assembly of the global matrix and enabling a matrix-free solution strategy. Within a unified framework that integrates diverse numerical approaches—including Kronecker product algebra, compact finite differences, isogeometric analysis, and direct diagonalization—the study systematically identifies three key techniques: multi-right-hand-side reshaping, sum factorization, and pencil-style MPI decomposition. These innovations collectively enhance hardware affinity and parallel scalability, reducing algorithmic complexity to O(N) and storage requirements to O(Nₓ + Nᵧ + N_z), thus enabling efficient large-scale 3D PDE simulations.
📝 Abstract
Every Cartesian three-dimensional PDE solver hides a structural secret that production CFD codes have used for half a century and that graduate-level textbooks rarely state plainly. The derivative matrices, the compact Padé line solves, the Galerkin mass inversions, the alternating-direction-implicit substeps, and even the fast Poisson and Helmholtz diagonalization transforms all factor along the coordinate axes and collapse into repeated one-dimensional banded kernels executed along the grid lines. The three-dimensional operator exists only on paper; it is never assembled, factored, or stored. This paper is the manual for that collapse. We derive the Kronecker-product algebra that makes it exact, carry it cleanly through central differences, compact schemes, tensor-product Galerkin, B-spline and isogeometric methods, collocation, ADI time stepping, and direct Poisson and Helmholtz solves, and bring into the open the three production tricks that turn the reduction into hardware-conscious floating-point throughput on real machines: the multi-right-hand-side reshape that exposes a sweep as one batched line kernel (a dense BLAS-3 GEMM when the line factor is dense or element-local, a banded or stencil kernel when it is not), the sum factorization that rescues high-order Galerkin from the $O(p^{2d})$ quadrature trap, and the pencil decomposition that keeps every direction contiguous across an MPI cluster. For fixed stencil width or fixed polynomial degree, the compute cost stays $O(N)$ in the total number of unknowns $N = N_x N_y N_z$; the operator storage drops to $O(N_x + N_y + N_z)$ up to bandwidth constants; direct separable Poisson and Helmholtz solvers add the expected transform cost; the line kernels are embarrassingly parallel. These facts are familiar to practitioners but rarely assembled in one place; this paper collects them and shows how to use them.