Space-time adaptive methods for parabolic evolution equations

TL;DR

The paper tackles efficient numerical solution of parabolic evolution equations (linear heat, reaction-diffusion, unsteady Stokes, Navier-Stokes) by developing a space-time adaptive, integral equation–based framework on a 2D periodic domain. It leverages heat kernel potentials, continuous fast Gauss transforms, and Helmholtz decompositions to form explicit or semi-implicit marching schemes, with automatic refinement/coarsening on level-restricted quadtrees. Key contributions include a unified, high-order, linear-scaling approach that avoids large implicit systems for linear problems and reduces nonlinear solves to pointwise scalar problems in the semilinear case, demonstrated across multiple canonical problems with strong adaptivity. The work lays a foundation for efficient, boundary-aware extensions and open-source implementations, enabling accurate simulations of complex multiscale parabolic phenomena on nonuniform meshes.

Abstract

We present a family of integral equation-based solvers for the heat equation, reaction-diffusion systems, the unsteady Stokes equation and the incompressible Navier-Stokes equations in two space dimensions. Our emphasis is on the development of methods that can efficiently follow complex solution features in space-time by refinement and coarsening at each time step on an adaptive quadtree. For simplicity, we focus on problems posed in a square domain with periodic boundary conditions. The performance and robustness of the methods are illustrated with several numerical examples.

Figures (5)

• Figure 1: [Adapted from fgt2024] The new FGT computes volume integrals discretized on a level-restricted quadtree with arbitrary order accuracy, determined by the order of the polynomial approximation on each leaf node. (a) For any precision $\epsilon$, the FGT defines the cutoff level $l$ in the tree hierarchy to be that where the box dimension $D_{cut}= 2^{-l}$ satisfies $e^{-D_{cut}^2/(4\delta)} < \epsilon$. In the figure, $B_2$ is a leaf node at the cutoff level, $B_1$ is a leaf node at a finer level and $B_3$ is a leaf node at a coarser level. [For readers familiar with hierarchical fast transforms: plane wave expansions are merged in an upward pass (from fine to coarse) until the cutoff level is reached. They are then translated to near neighbors at the cutoff level and communicated to finer levels in a downward pass. For coarse leaf nodes such as $B_3$, the continuous Gauss transform is evaluated in the near field directly, using precomputed tables for optimal performance.] An interesting aspect of the Gauss transform on non-uniform data structures is that the grid resolving the output may need to be finer than the grid resolving the input! In (b) we consider $u_0({\bf x})$ as the sum of two sharply peaked Gaussians. In (c), we plot the initial potential (the Gauss transform of $u_0$) for $\delta = 4 \times 10^{-3}$. In (d), we show the quadtree that resolves the input data to $12$ digits of accuracy. In (e), we show the quadtree that resolves the output data to $12$ digits of accuracy after convolving with a sharp Gaussian. The reason for this surprising behavior is that a sharp feature can diffused a short distance and reach regions where it is too sharp to be well resolved by the original adaptive grid, indicating that careful monitoring of resolution is essential. The total number of boxes in (e) is actually slightly less than in (d), after refinement where more resolution is needed and coarsening where the resolution is determined to be excessive. See section \ref{['coarseningsec']} for further details. \newlabelfig:adaptree0
• Figure 1: Evolution of the solution to the inhomogeneous, periodic heat equation with forcing function defined in \ref{['heat_forcing']}. Each subfigure displays both the numerical solution and the corresponding level-restricted quadtree at selected time steps. The adaptive trees are generated using our automatic refinement and coarsening strategy, and functions on each leaf node are discretized using a scaled $8 \times 8$ Chebyshev grid.
• Figure 2: Evolution of the $v$-component of the Gray-Scott reaction-diffusion system at selected time steps. Each panel shows the numerical solution overlaid with the corresponding adaptive quadtree mesh. The adaptive meshes are generated using our automatic refinement and coarsening strategy with an error tolerance $\epsilon = 10^{-9}$. The spatial discretization employs an eighth-order tensor product Chebyshev approximation on each leaf node.
• Figure 3: Evolution of the vorticity for the unsteady Stokes equations\ref{['eq-us2']} on $D = [-0.5, 0.5]^2$ at selected time steps up to the final time $T=0.1$. Each subfigure shows the numerical solution overlaid with the corresponding adaptive quadtree mesh.
• Figure 4: Evolution of the vorticity field for the double shear layer problem\ref{['double sheer layer_eq']} at selected time steps with $\Delta t=0.0008s$. Each subfigure shows the numerical solution overlaid with the corresponding adaptive quadtree.