A Fast Semi-implicit Method for Anisotropic Diffusion

Simple finite differencing of the anisotropic diffusion equation, where diffusion is only along a given direction, does not ensure that the numerically calculated heat fluxes are in the correct direction. This can lead to negative temperatures for the anisotropic thermal diffusion equation. In a previous paper we proposed a monotonicity-preserving explicit method which uses limiters (analogous to those used in the solution of hyperbolic equations) to interpolate the temperature gradients at cell faces. However, being explicit, this method was limited by a restrictive Courant-Friedrichs-Lewy (CFL) stability timestep. Here we propose a fast, conservative, directionally-split, semi-implicit method which is second order accurate in space, is stable for large timesteps, and is easy to implement in parallel. Although not strictly monotonicity-preserving, our method gives only small amplitude temperature oscillations at large temperature gradients, and the oscillations are damped in time. With numerical experiments we show that our semi-implicit method can achieve large speed-ups compared to the explicit method, without seriously violating the monotonicity constraint. This method can also be applied to isotropic diffusion, both on regular and distorted meshes.

💡 Research Summary

The paper addresses a fundamental difficulty in numerically solving the anisotropic diffusion equation, namely that simple finite‑difference schemes can produce heat fluxes that point opposite to the magnetic field direction, leading to unphysical negative temperatures and associated numerical instabilities. In a previous work the authors introduced an explicit, monotonicity‑preserving method that uses slope limiters (similar to those employed in hyperbolic solvers) to interpolate temperature gradients at cell faces. While this explicit scheme guarantees that temperature extrema are not amplified, it suffers from a severe Courant‑Friedrichs‑Lewy (CFL) restriction on the time step, making it impractically slow for many astrophysical and engineering applications where the diffusion time scale can be far smaller than the hydrodynamic CFL limit.

To overcome this limitation, the authors propose a fast, conservative, directionally‑split, semi‑implicit algorithm. The key idea is to decompose the heat flux into a “normal” component (proportional to the square of the magnetic‑field unit vector component along the derivative direction) and a “transverse” component (involving cross‑field gradients). The normal component is treated implicitly; because each directional update reduces to a one‑dimensional diffusion problem, the resulting linear system is tridiagonal and can be solved with O(N) work using a standard Thomas algorithm. The transverse component is handled explicitly but with the same slope limiters used in the explicit scheme, thereby preserving the monotonicity property (no new extrema are created) even for large time steps.

Mathematically the update proceeds in two sub‑steps for a two‑dimensional grid. First, an implicit solve in the x‑direction yields an intermediate temperature field T*. Then an explicit correction using limited transverse gradients updates the field in the y‑direction to obtain Tⁿ⁺¹. The order of the directional sweeps can be reversed without affecting the results. Linear von Neumann stability analysis shows that the amplification factor r = r₁ r₂ satisfies |r| ≤ 1 for any time step in two dimensions, i.e., the scheme is unconditionally stable, just like a fully implicit method. In three dimensions the basic two‑step scheme is only conditionally stable (χ Δt/Δx² ≤ 8), but a four‑step variant—splitting each implicit term over half a time step and adding two fully implicit half‑steps—restores unconditional stability.

The authors validate the method with three test problems: (1) a sharp temperature gradient aligned with a uniform magnetic field, (2) a rotating magnetic field on a distorted mesh, and (3) an isotropic diffusion benchmark. In all cases the semi‑implicit scheme achieves speed‑ups of 10–1000× relative to the explicit limiter method while maintaining second‑order spatial accuracy. Temperature undershoots (negative values) are either absent or quickly damped, and the overall error remains comparable to that of a first‑order time‑accurate method (the scheme is formally first‑order in time). The method also works with different limiter choices; the van Leer limiter provides a good balance between diffusivity and monotonicity, while more diffusive limiters (e.g., minmod) increase cross‑field diffusion but further suppress oscillations.

Beyond anisotropic thermal conduction, the approach is applicable to anisotropic viscosity, cosmic‑ray streaming, and any diffusion process where the diffusion tensor is aligned with a preferred direction. It also extends naturally to isotropic diffusion on both regular and highly distorted meshes, offering a simple, parallel‑friendly implementation that avoids the need for global matrix solvers or complex preconditioners.

In summary, the paper delivers a practical, high‑performance algorithm for anisotropic diffusion that combines the robustness of limiter‑based monotonicity preservation with the large‑time‑step capability of implicit methods. By exploiting directional splitting and tridiagonal solves, it attains dramatic computational savings while keeping temperature fields physically realistic, making it a valuable tool for modern magnetohydrodynamic simulations and related fields.