using MatrixFreeOperators, LinearAlgebra
base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8))
forest = BlockForest(base; blocksize = (4, 4), maxlevel = 4) # 2×2 root tiling
length(collect(leaves(forest))) # 4 leaf blocks at level 04Adaptive Mesh Refinement
BlockForest adds adaptive resolution while keeping every existing operator working unchanged. This page tells the story of the design — the problem, the approach, and the one rule everything rests on — and shows the current capabilities.
The motivation: stop wasting cells
A uniform grid wastes effort. Most of a simulation is smooth and undemanding, but a few regions — a shock front, a flame, a boundary layer — need fine cells. Adaptive Mesh Refinement (AMR) concentrates resolution only where it is needed. The hard constraint here: do that without rewriting every operator, because the whole package assumes a field is a dense rectangular array you sweep a stencil across.
Uniform grid — fine EVERYWHERE Adaptive forest — fine ONLY where needed
(most cells do nothing useful) (resolution follows the feature ▓)
┌─┬─┬─┬─┬─┬─┬─┬─┐ ┌───────┬───────┐
├─┼─┼─┼─┼─┼─┼─┼─┤ │ │ │
├─┼─┼─┼─┼─┼─┼─┼─┤ │ │ │
├─┼─┼─┼─┼─┼─┼─┼─┤ ├───┬───┼───────┤
├─┼─┼─┼─┼─┼─┼─┼─┤ │▓▓▓│ │ │
├─┼─┼─┼─┼─┼─┼─┼─┤ ├─┬─┼───┤ │
├─┼─┼─┼─┼─┼─┼─┼─┤ │▓│▓│ │ │
└─┴─┴─┴─┴─┴─┴─┴─┘ └─┴─┴───┴───────┘A forest of fixed-size leaf-blocks
Rather than refine individual cells, the domain is tiled by many equal-sized small blocks — e.g. each block is an 8×8 (2-D) or 8×8×8 (3-D) patch of cells, an ordinary CartesianGrid with its own halo. To add resolution somewhere, you replace one block with 2ᴺ child blocks covering the same physical area at half the spacing, and recurse. Each root block and its descendants form a tree; the collection of trees is the forest. A leaf is a block nobody refined further — the leaves tile the domain with no gaps or overlaps, and that is where the solution lives.
Spatial view (2×2 root tiling) Tree view
Level 0 refine BL corner root₀ root₁ root₂ root₃
┌────┬────┐ ┌────┬────┐ │
│ r0 │ r1 │ │ r0 │ r1 │ refine r0
├────┼────┤ ─▶ ├─┬─┬┼────┤ │
│ r2 │ r3 │ ├─┼─┤│ │ ┌───┬──┴─┬───┐
│ │ │ ├─┼─┤│ r3 │ c0 c1 c2 c3 (4 leaves,
└────┴────┘ └─┴─┴┴────┘ half the spacing)Fixed-size blocks keep every leaf a plain array
Every block always holds the same cell count — only the physical spacing shrinks with depth (level ℓ → h₀ / 2^ℓ). Because each leaf is still a dense rectangular array, every existing leaf operator runs unchanged per block — the same stencil sweep, broadcasts, and KernelAbstractions kernels. All of the adaptivity complexity is confined to the block faces (the coarse–fine interfaces).
Level-0 block Level-1 child block
spacing = h spacing = h/2
┌─┬─┬─┬─┐ ┌┬┬┬┐
├─┼─┼─┼─┤ ├┼┼┼┤ same 4×4 cells,
├─┼─┼─┼─┤ ├┼┼┼┤ HALF the physical footprint
├─┼─┼─┼─┤ ├┼┼┼┤ ⇒ twice the resolution
└─┴─┴─┴─┘ └┴┴┴┘
4×4 dense array 4×4 dense array ← still a plain array!This is why the design picks the block-structured school (FLASH / PARAMESH / AMReX) and rejects a cell-based octree (p4est / Trixi style): the octree puts a hanging node at nearly every cell, which would break the dense-array stencil sweeps the whole package is built on. The block-structured choice trades a little redundant memory at block boundaries for keeping the interior of every block a plain array.
The load-bearing rule: 2:1 balance
The invariant everything rests on is 2:1 balance: any two neighboring leaves stay within one refinement level of each other. This reduces every coarse–fine interface to just three cases — same level / one-coarser / one-finer — which is what makes interface interpolation tractable. refine! and coarsen! re-establish balance automatically (you can also call balance! directly).
LEGAL (neighbors ≤ 1 level apart) ILLEGAL (2 levels apart)
┌─────────┬────┬────┐ ┌─────────┬──┬──┐
│ │ L1 │ L1 │ │ ├──┼──┤
│ L0 ├────┼────┤ │ L0 ├──┼──┤ L2 directly
│ │ L1 │ L1 │ │ ├──┼──┤ touching L0
└─────────┴────┴────┘ └─────────┴──┴──┘ → forbiddenWhere the intelligence lives
Operators barely change. The forest does a single halo_update! over the whole forest to fill each block’s interface ghosts from its neighbors, then a single forest-level apply_bc! pass fills the physical domain-boundary ghosts from precomputed face lists; each leaf then runs its ordinary stencil. So the only genuinely new operator code is those two forest-level passes — plus their declared transposes (halo_update_adjoint! and the forest-level fold_bc!), which keep the adjoint identity ⟨Lx,y⟩ = ⟨x,Lᵀy⟩ exact on a forest so Krylov solvers stay correct.
forest
│
▼
┌───────────────────────────────────────────────┐
│ 1. halo_update!(forest) │ ONCE, whole forest —
│ fill each block's INTERFACE ghosts from │ flat descriptor loop,
│ its neighbors │ no topology queries
└───────────────────────────────────────────────┘
│
▼
┌───────────────────────────────────────────────┐
│ 2. apply_bc!(forest) │ ONCE, whole forest —
│ fill each boundary block's PHYSICAL ghosts │ per-(dim, side)
│ from precomputed face lists │ leaf-index lists
└───────────────────────────────────────────────┘
│
├───────────────┬───────────────┐
▼ ▼ ▼
┌───────────┐ ┌───────────┐ ┌───────────┐
│per leaf: │ │per leaf: │ │per leaf: │
│3 stencil │ │3 stencil │ │3 stencil │ ← EXISTING single-grid code
└───────────┘ └───────────┘ └───────────┘Every leaf grid carries the Interface boundary on all faces, so per-leaf BC sweeps are no-ops and every leaf shares one concrete grid type — the physical BCs live on the forest, applied by the face pass above. This is the homogeneous-blocks structure production AMR codes (AMReX, p4est) use: boundary operators are separate passes over face lists, never per-block special cases.
Usage
A BlockForest wraps a base CartesianGrid (the root tiling). Operators and fields are built on it exactly as on a single grid:
Refine a region with a predicate on each block’s center; balance is restored automatically:
refine!(forest, x -> x[1] < 0.3 && x[2] < 0.3) # refine near the origin corner
length(collect(leaves(forest))) # grew, still 2:1 balanced7Parity with a single grid
On a uniformly-refined forest, the operators reproduce a single CartesianGrid of the same resolution exactly. The interior values are identical (just stored in a different, block-by-block order), so any order-independent reduction matches:
g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (8, 8))
bf = BlockForest(g; blocksize = (4, 4), maxlevel = 3) # level 0 ⇒ same 8×8
u = set!(scalar_field(g), x -> sinpi(x[1]) * sinpi(x[2]))
uf = set!(scalar_field(bf), x -> sinpi(x[1]) * sinpi(x[2]))
norm(flatten(laplacian(g) * u)) ≈ norm(flatten(laplacian(bf) * uf))trueIndicator-driven regridding
The center-predicate refine! above is geometric. To adapt to the solution, regrid! evaluates per-block criteria on a live field — each criterion receives the leaf as an ordinary Field, so it is typically a reduction like b -> maximum(abs, interior(b)) > τ — edits the topology in a single pass (refine-marked leaves split, complete fully-marked families coarsen, refine wins conflicts), re-establishes 2:1 balance, and returns freshly allocated fields with the data carried across:
bump = set!(scalar_field(forest), x -> exp(-((x[1] - 0.5)^2 + (x[2] - 0.5)^2) / 0.01))
bump = regrid!(bump; refine = b -> maximum(abs, interior(b)) > 0.5)
length(collect(leaves(forest))) # resolution now follows the bump16The transfer is interior-only and keyed by leaf identity: an unchanged leaf is copied, a refined leaf is filled by per-dimension linear interpolation from its old parent (one-sided at parent-block edges), and a coarsened leaf by the conservative 2⁻ᴺ mean of its old children. Ghosts of the returned fields are scratch, refilled by halo_update!/apply_bc! on the next operator application.
Because a real regrid invalidates every field allocated on the old leaf set (the generation guards throw on stale use), any field needed afterwards must ride the same call — compute indicator fields before regridding and pass them first:
η = magnitude_of(gradient(forest) * u) # indicator, pre-regrid
η, u = regrid!(η, u; refine = b -> maximum(interior(b)) > τ,
coarsen = b -> maximum(interior(b)) < τ / 10)When the marks change nothing, regrid! returns the same field objects and existing prepared operators remain valid — the natural convergence signal for an adaptive loop.
The adaptive solve loop
The loop is deliberately user code — solver, tolerances, indicator, and stopping rule are all problem-specific. The canonical shape (in full in examples/adaptive_poisson.jl):
for cycle in 1:ncycles
P = prepare(laplacian(forest), u) # a regrid invalidates P: rebuild each cycle
sol, stats = Krylov.gmres(P, rhs) # nonsymmetric on an adapted forest — not cg
flat_to_interior!(u, sol)
η = magnitude_of(gradient(forest) * u)
η, u = regrid!(η, u; refine = b -> maximum(interior(b)) > τ)
endTwo rules keep it honest: the indicator is computed before the regrid (its field rides the transfer), and everything tied to the old leaf set — fields, PreparedForests — is rebuilt after; using stale state throws rather than silently reading wrong-shaped storage.
For the same loop driving a time-dependent 3D problem, with the adaptive answer checked against an equal-resolution uniform run and an independent code, see the Niederer Benchmark.
Status
Every AMR phase is built, including the packed-storage performance phase: the forest is adaptive end-to-end (solve → indicator → regrid! → re-prepare), and on GPU backends the packed layout sweeps all leaves in single KernelAbstractions launches with a batched device halo exchange.
Phase 0 ██████████ DONE forest topology + BlockForest grid + BlockField
Phase 1 ██████████ DONE inter-block halo_update! + operator PARITY
Phase 2 ██████████ DONE coarse–fine INTERPOLATION (non-uniform forests)
Phase 3 ██████████ DONE Restriction/Prolongation + V-cycle MULTIGRID
Phase 4 ██████████ DONE indicator-driven REGRID driver (regrid!)
Phase 5 ██████████ DONE packed storage (PackedBlockField) + Laplacian kernel
██████████ DONE per-operator forest kernels + kernelized exchangepack(u) converts a BlockField into a PackedBlockField — every leaf in one contiguous buffer — and prepare(L, pack(u)) selects the packed path by dispatch (unpack converts back). On GPU backends the packed sweep is a single kernel launch over all leaves instead of one per leaf: measured on an RTX 4000 Ada (benchmark/gpu.jl), packed mul! runs 30–95× faster than the per-leaf device path at equal DOFs, host allocations stay flat in the leaf count, and the batched exchange kernels beat the per-descriptor loop on the same device data by 230–1370×. On CPU it runs the same per-leaf fused broadcasts as the reference layout, so packed storage costs nothing extra there. The reference BlockField remains the path for AD, Reactant, and regrid-heavy workflows, since a packed field must be re-packed after each regrid.