🔢 mathematics

A Fully Matrix-Free Three-Grid Preconditioner for the Time-Harmonic Maxwell Equations at Extreme Scale

This paper presents a fully matrix-free, factorization-free three-grid preconditioner for the time-harmonic Maxwell equations that achieves extreme scalability, solving systems with over 10 billion unknowns in under 72 seconds on 64 NVIDIA A100 GPUs by strategically confining complex shifts to an auxiliary preconditioner while preserving the physical operator.

Original authors: Shubin Fu

Published 2026-08-25

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Shubin Fu

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine trying to map the invisible ripples of light or radio waves as they travel through a complex world. Whether designing a new antenna, imaging the inside of the human body, or understanding how radar bounces off a stealth aircraft, scientists must solve equations that describe how these electromagnetic waves move. When these waves travel far, the math becomes incredibly difficult. The computer models required to track them accurately are so massive that they often overwhelm even the most powerful supercomputers. The problem is not just the size of the data, but the nature of the waves themselves; they oscillate rapidly, and to capture their path without error, the computer must divide space into tiny, tiny pieces. As the distance the wave travels increases, the number of these pieces grows so fast that the memory required to store the calculations becomes impossible to manage.

For decades, researchers have tried to solve these massive systems using a strategy called multigrid, which works like a team of editors checking a manuscript at different levels of detail. They solve the problem on a coarse, low-resolution grid first to get a rough idea, then refine that answer on finer and finer grids. However, for these specific wave equations, the standard approach hits a wall. The final step usually requires a direct calculation on the coarsest grid that is so memory-intensive it stops the whole process from scaling up. Furthermore, many existing methods try to simplify the physics by altering the equations, which can introduce errors or fail to capture the true behavior of the waves.

In a recent breakthrough, a researcher has developed a new way to solve these problems that avoids these bottlenecks entirely. They created a method that never stores the full mathematical matrix—the giant grid of numbers that usually represents the system—in the computer's memory. Instead, the computer calculates the necessary values on the fly, only when needed. This "matrix-free" approach allows the system to handle problems of a scale previously thought unmanageable. The researcher tested their new solver on a set of four different scenarios, ranging from a simple, uniform space to complex environments with lenses, repeating patterns, and random variations in the material. They ran these tests on a cluster of sixty-four powerful graphics processors, a setup typically used for artificial intelligence, to see if it could handle the extreme demands of electromagnetic simulation.

The results were striking. The researcher successfully solved systems containing approximately 10.89 billion unknowns, which represents the number of points in the grid where the wave's behavior is calculated. In the most complex case, involving a random, three-dimensional material, the entire calculation was completed in just 72 seconds. Even in the simplest case, the solution took only 42 seconds. These times include the setup and the solving process, all while using less than 38 gigabytes of memory per graphics processor. This is a massive reduction compared to traditional methods, which would require terabytes of memory and hours or days of computation for a problem of this size.

The secret to this speed lies in how the researcher structured their solver. Instead of trying to force the computer to solve the hardest part of the equation directly, they split the work into a clever three-step process. They use a primary solver to handle the main physical problem, which remains untouched and accurate. To help this solver along, they employ a secondary, "helper" system that is slightly modified to be easier to solve. This helper system acts as a guide, correcting the main solver's path without ever changing the fundamental physics of the original problem. By keeping the main equations pure and using the helper only for guidance, they avoid the memory crashes that usually occur when trying to solve these specific wave equations at such large scales.

The researcher also verified that their method was not just fast, but accurate. They compared their results against known mathematical solutions for a simple case and found that the simulated waves matched the theoretical predictions with high precision. They tested the method on different types of materials, including a converging lens that bends waves and a random medium that scatters them in all directions. In every case, the solver maintained its speed and stability, proving that the method works regardless of how complex the environment is. The researcher demonstrated that by carefully balancing the computational steps and avoiding the storage of massive data tables, it is possible to simulate the behavior of light and radio waves across vast distances in a matter of seconds.

This work suggests that the limitations of current electromagnetic simulations are not a fundamental barrier, but a solvable engineering challenge. By rethinking how the calculations are organized and leveraging the specific strengths of modern graphics hardware, the researcher has opened the door to simulating much larger and more realistic scenarios. Whether for designing better communication networks, improving medical imaging, or understanding radar signatures, the ability to solve these equations quickly and with modest resources could transform how engineers and scientists approach problems involving wave propagation. The study confirms that with the right algorithmic design, the extreme scale of these problems is no longer an obstacle, but a manageable task.

Technical Summary: A Fully Matrix-Free Three-Grid Preconditioner for Time-Harmonic Maxwell Equations

Problem Statement
Three-dimensional time-harmonic Maxwell simulations generate massive, complex, indefinite linear systems. Maintaining phase accuracy in high-frequency regimes requires fine spatial resolution, causing the number of unknowns to grow cubically with propagation distance. These systems are inherently difficult to solve due to the massive gradient kernel of the curl–curl operator, the indefiniteness of the wave term, and the complex anisotropic coefficients introduced by Perfectly Matched Layer (PML) truncation.

Existing solvers face a scalability bottleneck. Sparse direct solvers suffer from superlinear memory fill-in and communication overhead in 3D. While domain-decomposition and sweeping preconditioners offer low iteration counts, their directional dependencies and factorization storage requirements are poorly aligned with the massive parallel throughput of modern GPUs. Standard geometric multigrid methods, while naturally suited for GPUs, are ultimately stalled by the memory and communication costs of exact coarse-grid factorizations, which become prohibitive at extreme scales.

Methodology
The author presents a fully matrix-free, factorization-free three-grid preconditioner designed specifically for curl-conforming Nédelec finite element discretizations with optimally blended quadrature. The core innovation is a strategic decoupling of the physical problem from the stabilization mechanism:

Three-Grid Hierarchy: The solver operates on three nested hexahedral meshes ( $h$ , $2h$ , $4h$ ).
- Fine Grid ( $h$ ): Solves the unshifted physical Maxwell equation ( $A_h$ ) using an outer FGMRES solver.
- Intermediate Grid ( $2h$ ): Computes an inexact correction via an inner FGMRES process. Crucially, the Krylov operator here remains the unshifted physical operator ( $A_{2h}$ ), but it is preconditioned by a shifted auxiliary cycle.
- Coarse Grid ( $4h$ ): The bottom level of the auxiliary cycle uses a shifted operator ( $B_{4h}$ ) to dampen under-resolved components.
Shifted Auxiliary Cycle: To stabilize the indefinite wave equation without altering the physical spectrum, a complex shift ( $\sigma$ ) is applied only to the auxiliary hierarchy ( $B_\ell = C_\ell - (1+i\sigma)k_0^2 M_\ell$ ). This shift is strictly confined to the preconditioner for the intermediate-grid correction. The primary physical operator ( $A_h$ ) remains untouched, preserving the physical Maxwell operator for the outer solver.
Matrix-Free Implementation: The method avoids assembling global matrices or performing exact coarse-grid factorizations.
- Operators: All levels use rediscretized operators via optimally blended quadrature (suppressing dispersion error to fourth order) rather than Galerkin products.
- Transfers: Edge transfers (prolongation and restriction) are performed geometrically using coordinate-based interpolation without forming global transfer matrices.
- Smoothers: Fixed-work Krylov smoothers (GMRES with damped Jacobi) are used at all levels, eliminating the need for inner convergence checks and keeping memory footprint proportional strictly to the number of edge unknowns.
Parameter Selection: Local Fourier Analysis (LFA) is employed to derive blended Maxwell branches and compatible edge transfers. This analysis identifies robust parameters for the complex shift ( $\sigma = 0.4$ ) and Jacobi damping factors, ensuring the shifted auxiliary inverse remains bounded and free of real-frequency singularities while confining the preconditioned spectrum to a favorable region.

Key Contributions

Architectural Design: A novel three-grid preconditioner that isolates the complex shift to an auxiliary inner cycle, allowing the outer solver to operate on the unshifted physical operator. This resolves the tension between stabilizing indefinite wave equations and preserving physical accuracy.
Extreme Scalability: The algorithm achieves extreme scalability by eliminating factorization storage and global assembly. It relies solely on geometric edge transfers and fixed-work Krylov smoothers.
Optimized Discretization: The use of optimally blended quadrature with lowest-order Nédelec elements provides superior phase fidelity without increasing the system size or introducing auxiliary unknowns.
GPU Efficiency: The implementation is tailored for GPU throughput, utilizing element-based matrix-free kernels that minimize global memory bandwidth usage compared to stored-stencil approaches.

Results
The method was validated against analytical Maxwell Green tensors and tested on four distinct media: homogeneous, converging lens, periodic, and random isotropic.

Scale: The solver successfully handled systems with approximately 10.89 billion complex edge unknowns.
Performance: Using a single, uniform solver configuration (no retuning for specific media), the largest systems were solved in 42.0 to 72.0 seconds on just 64 NVIDIA A100 GPUs.
Efficiency:
- Weak Scaling: Efficiency ranged from 69.6% to 91.0% as the problem size and GPU count increased (up to 64 GPUs).
- Strong Scaling: Parallel efficiency remained robust (72.0%–93.8%) when doubling GPU counts for fixed problem sizes.
- Memory: Peak memory usage remained below 37.7 GiB per GPU.
- Setup: Matrix-free setup times were under 0.8 seconds for all configurations.
Precision: Complex single precision was found sufficient, preserving numerical results with negligible deviation from double precision while offering significant performance benefits.

Significance
The paper claims that these calculations represent the largest reported 3D time-harmonic Maxwell simulations utilizing curl-conforming finite elements. The significance lies in demonstrating that extreme-scale electromagnetic simulations can be achieved without the traditional bottlenecks of coarse-grid factorization or massive memory footprints. By confining the complex shift to an auxiliary preconditioner and leveraging a fully matrix-free, factorization-free architecture, the method achieves unprecedented speed and scalability on modern GPU hardware, solving problems that were previously intractable at this scale with this level of discretization fidelity.

Technical Summary: A Fully Matrix-Free Three-Grid Preconditioner for Time-Harmonic Maxwell Equations

More like this