Quadtree AMR — Letting the Mesh Grow Where It Needs To

Stanford, 1984. Marsha Berger and Joseph Oliger were chasing one shock. The mesh held a million cells. About 1% of them, the ones on the shock, did any real work. The other 99% kept recomputing a uniform free stream at the same precision. The Journal of Computational Physics paper from that year is Berger–Oliger AMR. The same habit is still inside OpenFOAM's dynamicRefineFvMesh, AMReX, Chombo, and p4est.

Drag the white dot and the mesh splits around it; the threshold slider changes the scaled slope $\tau$ that decides which cells cut.

Drop $\tau$ to 0.05 and leaves chase the Gaussian tails. Push it past 0.4 and the mesh stays coarse. add feature drops a second peak, and only those two cores go deep. The upper-left readout is the cell count against a uniform mesh.

Why a Uniform Grid Wastes 90%#

A 3D compressible RANS run makes the trade obvious. Double the spacing and one step gets 8 times cheaper, because $2^3=8$. Shocks, shear layers, and separation smear with that choice. Refine the whole domain for the thin band and the free-stream cells shrink too. Memory and time blow up together.

Refine the band. Leave the rest alone. The ratio is this estimate.

$L$ is the deepest refinement level. $V_\text{feature}$ is the volume that needs the fine spacing. $V_\text{domain}$ is the whole box. In 2D, $4^{L}$ is what one level does to the cell count: four children. At $L=5$, $4^5=1024$. Take the opening 1% as the volume fraction. Against the coarsest grid the cell count is $0.01\times 1024+0.99=11.23$ times larger. A uniform mesh at that same level is 1024 times the coarse grid. Divide the two and AMR sits near 1.1% of the uniform fine mesh. A shock with $V_\text{feature}/V_\text{domain}\sim h$ pushes the share down further. The single-digit percent is that division. In 3D, halving one edge multiplies cells by 8. Spread those 8 across the whole domain and the bill covers everything. Keep them inside the band and the factor of 8 stops there.

h-, p-, r-: Three Flavors of Adaptation#

An adaptive mesh — one that changes with the solution — splits by what actually moves.

• h-adaptivity: shrink the cell size $h$. This is the usual meaning of AMR.
• p-adaptivity: keep the cells, raise the polynomial order $p$. Common in DGM and SEM (spectral element method).
• r-adaptivity (relocation): keep the cell count and the order, slide nodes toward steep regions. Memory stays flat. Distortion does not.

The rest of this note stays on h-AMR, and on the 2D quadtree. In 3D the same tree is an octree.

Quadtree: Cells That Own Their Children#

A uniform grid is an array. An adaptive grid needs a tree. One refine call makes four children, and each child points at its parent. Refine a child and the grandchildren show up.

class QuadCell:
    __slots__ = ("cx", "cy", "size", "level", "children", "parent", "data")
 
    def __init__(self, cx, cy, size, level=0, parent=None):
        self.cx = cx           # cell-center x
        self.cy = cy           # cell-center y
        self.size = size       # edge length
        self.level = level     # root = 0, refine adds 1
        self.children = None   # None if leaf, length-4 list otherwise
        self.parent = parent
        self.data = 0.0        # conserved quantity (e.g. density)
 
    def is_leaf(self):
        return self.children is None
 
    def split(self):
        h = self.size * 0.5
        q = self.size * 0.25
        offsets = ((-q, -q), (+q, -q), (-q, +q), (+q, +q))
        self.children = [
            QuadCell(self.cx + dx, self.cy + dy, h, self.level + 1, self)
            for dx, dy in offsets
        ]
        for ch in self.children:
            ch.data = self.data    # initialize from parent (conservative)

Copying the parent value in split is about conservation. The parent's integral $U \cdot V$ has to equal the sum over the four children. Each child has a quarter of the parent's volume. Copy the value straight across and the four pieces add as $4 \times (U \cdot V/4) = U \cdot V$. That copy is only first-order interpolation. Production codes usually add the parent's slope, MUSCL-style, and call the result second-order prolongation. The slope can go in. The four-way sum still has to match $U \cdot V$. If it does not, the split creates or destroys mass.

Where to Cut — the Refinement Criterion#

A quadtree is an empty box. The criterion decides which cell is cut. Five forms show up constantly.

dynamicRefineFvMesh in OpenFOAM thresholds a field the user picked. AMReX leaves a per-cell Boolean, tagcells(), for the caller to fill. What gets tagged changes the tree more than the container does. Keeping the same answer while cutting the cell count continues in tagging and coarse-fine reflux.

The browser mesh uses only the first row. A central difference of $|\nabla u|$ at the cell center, times the cell size $h$. Above $\tau$, the cell splits. Otherwise it stays a leaf.

2:1 Balance and Hanging Nodes#

The face where a child meets a neighbor is the awkward one. One side is a large cell. The other side is two small cells. That face grows a hanging node. In finite volume a hanging node breaks flux balance.

Two ways through.

1. Enforce 2:1 balance. Neighboring leaves may differ by at most one level. p4est, AMReX, and the Berger codes all do this.
2. Non-conforming interpolation. Build the coarse face value from the two small cells. Usual in DGM-style codes.

The enforcement loop is short. Walk the refine candidates. If a neighbor leaf is two or more levels coarser, push it onto the queue. Repeat until the queue is empty.

def enforce_2to1(roots, refine_queue):
    """ensure every cell in refine_queue is 2:1 with its neighbors"""
    pending = list(refine_queue)
    queued = set(id(c) for c in pending)
    while pending:
        cell = pending.pop()
        for nb in face_neighbors(cell, roots):
            if nb.is_leaf() and cell.level - nb.level >= 1:
                if id(nb) not in queued:
                    pending.append(nb)
                    queued.add(id(nb))
    return [c for c in (refine_queue + pending) if c.is_leaf()]

face_neighbors returns leaves that share a face. The walk is $O(\log N)$.

A 1D Adaptive Grid in Python#

A full 2D quadtree listing gets long. The same idea fits in a 1D binary tree for one pass. Plant a Gaussian pulse. Bisect only the cells that clear the threshold.

import numpy as np
 
def gaussian_pulse(x, x0=0.6, sigma=0.04):
    return np.exp(-((x - x0) ** 2) / (2 * sigma ** 2))
 
def refine_indicator(cell_lo, cell_hi, fn, h_scale=1.0):
    """scaled first derivative = absolute difference across cell"""
    u_lo = fn(cell_lo)
    u_hi = fn(cell_hi)
    return abs(u_hi - u_lo) * h_scale
 
def adaptive_bisect_1d(fn, tau=0.05, max_level=8):
    """recursively bisect [0, 1] based on the threshold"""
    leaves = []
    def recurse(lo, hi, level):
        ind = refine_indicator(lo, hi, fn)
        if level >= max_level or ind < tau:
            leaves.append((lo, hi, level))
            return
        mid = 0.5 * (lo + hi)
        recurse(lo, mid, level + 1)
        recurse(mid, hi, level + 1)
    recurse(0.0, 1.0, 0)
    return leaves
 
cells = adaptive_bisect_1d(gaussian_pulse, tau=0.05, max_level=7)
print(f"leaf cells: {len(cells)}")
print(f"vs uniform level 7: {len(cells) / 2**7 * 100:.1f}%")
# example output: leaf cells: 41, vs uniform level 7: 32.0%

One Gaussian occupies a thin slice, so the run stops at 41 leaves: 32.0% of a uniform level-7 line. Two pulses send only the two cores deep. The gap stays coarse.

Between 41 Leaves and 128 Cells#

Divide the printed 41 by $2^7=128$ and the quotient is 0.3203. The comment's 32.0% is that number, cut to one decimal place. Subtract and 87 cells never get built at the finest spacing, about 68% of the uniform level-7 line. In the code, $\sigma=0.04$ is 4% of the unit interval. Only a neighborhood of that pulse clears $\tau=0.05$. Flat free stream has almost no difference between the two ends, so bisection stops early. One 1D cut makes two children. One 2D cut makes four, so the same thin band drops faster once $4^{L}$ from the estimate is in play. Raise $\tau$ from 0.05 to 0.4 and the tails that sat between those two numbers stop splitting.

A Pitfall in Parallel — Load Balancing and SFCs#

Distributed memory hits a different wall. One rank holding every deep branch sets the wall-clock time for everyone. The cell count is no longer balanced.

Line the tree up in 1D and cut the line. A space-filling curve (SFC) — Morton (Z-order) or Hilbert — keeps spatially neighboring cells close on that line. p4est's Morton index interleaves the bits of the coordinates.

def morton_2d(ix, iy, level):
    """interleave integer coordinates (ix, iy) into a Morton code"""
    code = 0
    for b in range(level):
        code |= ((ix >> b) & 1) << (2 * b)
        code |= ((iy >> b) & 1) << (2 * b + 1)
    return code
 
# Z-order over a 4x4 grid
for iy in range(4):
    for ix in range(4):
        print(morton_2d(ix, iy, 2), end=" ")
    print()
# 0  1  4  5
# 2  3  6  7
# 8  9 12 13
# 10 11 14 15

Sort by Morton index, then slice the list across ranks. Neighbors land together and the messages shrink. Square domains still jump at the corners. A Hilbert curve avoids the jump and costs more to index, so Morton is the default in practice.

Tag the Wrong Field and 32% Is Still Expensive#

The tree shape is mostly settled. Four children, then a Morton cut so neighbors share a rank. On a 4×4 block the index starts 0, 1, 4, 5 because the coordinate bits are interleaved. The waste shows up after that. Tag density while the error sits in a shear layer, and those 41 leaves at level 7 pile up on the wrong band. The 32% belongs to a criterion aimed at $x_0=0.6$ with $\sigma=0.04$. Change the tagged field and $V_\text{feature}$ in the estimate moves with it, including the factor of 11.23 against the coarse grid.

Skipping 2:1 is a real choice. Cap the level gap at 1 and a hanging node sees at most two small cells on a face. Drop the cap and every face needs its own non-conforming interpolation. The face flux then drifts away from the parent integral $U\cdot V$. The tree can be settled and the time-step question is still open: which line should call refine. How far one misplaced line moves the error is the subject of the recursive time cycle in LBM AMR.