C++17 AMReX MPI / GPU BSD-3-Clause

AMReX Level-set Redistancing

Block-structured adaptive-mesh reinitialization ("redistancing") of an advected level-set scalar $\phi$, built on the AMReX framework. The library restores the signed-distance property $|\nabla\phi| = 1$ without moving the zero contour, and offers two interchangeable algorithms behind a common interface.

Doxygen API reference → Source on GitHub →

Overview

The level-set method tracks a moving interface as the zero contour of a scalar field $\phi$. Advection degrades the signed-distance property, so $\phi$ must be periodically reinitialized so that interface quantities (normals, curvature, band widths) stay accurate. This library performs that step on a full AMReX hierarchy, with correct coarse–fine coupling and dynamic regridding.

Class	Algorithm	Best for
`LevelSetRedistance`	PDE (Sussman): Godunov/WENO5 + TVD-RK, Russo–Smereka, Sussman–Fatemi	Corner-heavy / thin-feature interfaces; contour pinned each step
`GeometricRedistance`	Geometric (Ausas–Dari–Buscaglia): bounded-facet reconstruction + eikonal + band-local mass	When a reconstruction-based distance is preferred over a PDE relaxation

Building

2D and 3D need separate build trees. AMREX_SPACEDIM is a compile-time constant, so configure one tree per dimension.

# 2D, CPU + MPI
cmake -S . -B build   -DREDIST_DIM=2 -DREDIST_MPI=ON
cmake --build build   -j

# 3D
cmake -S . -B build3d -DREDIST_DIM=3
cmake --build build3d -j

# CUDA GPU (3D)
cmake -S . -B buildcuda -DREDIST_DIM=3 -DREDIST_GPU_BACKEND=CUDA
cmake --build buildcuda -j

AMReX is fetched automatically (pinned tag) if not found by find_package. Build the API docs with cmake --build build --target docs.

PDE (Sussman) redistancing — `LevelSetRedistance`

Integrates the Hamilton–Jacobi equation to steady state in a pseudo-time $\tau$:

$$ \frac{\partial \phi}{\partial \tau} = \operatorname{sgn}(\phi_0)\bigl(1 - |\nabla\phi|\bigr), \qquad \operatorname{sgn}(\phi_0) \approx \frac{\phi_0}{\sqrt{\phi_0^2 + \epsilon^2}}, $$

where $\phi_0$ is the field frozen at the start of the step and $\epsilon = \texttt{sign\_eps}\cdot\Delta x$. Fixed points satisfy the eikonal equation $|\nabla\phi| = 1$; because the driving term vanishes where $\phi_0 = 0$, the zero contour is (to discretization error) stationary. See LevelSetRedistance in the API reference.

Godunov Hamiltonian

The gradient magnitude uses a sign-dependent upwind selection. With backward ($a$) and forward ($b$) one-sided differences per axis:

$$ |\nabla\phi|^2 = \begin{cases} \sum_x \max\!\bigl(\max(a,0)^2,\ \min(b,0)^2\bigr) & \operatorname{sgn}(\phi_0)>0,\\[4pt] \sum_x \max\!\bigl(\min(a,0)^2,\ \max(b,0)^2\bigr) & \operatorname{sgn}(\phi_0)<0. \end{cases} $$

Select with Redist.scheme = godunov1 (first order, robust; needs 2 ghost cells).

Fifth-order HJ-WENO5

For higher accuracy in smooth regions the biased derivatives $\phi_x^\pm$ are reconstructed as a convex combination of three third-order stencils,

$$ \phi_x = \sum_{k=1}^{3} w_k p_k,\qquad w_k = \frac{\alpha_k}{\sum_l \alpha_l},\qquad \alpha_k = \frac{d_k}{(\beta_k + \varepsilon)^2}, $$

with ideal weights $d=(0.1,0.6,0.3)$ and smoothness indicators $\beta_k$. Select with Redist.scheme = weno5 (needs $\ge 3$ ghost cells; regularization $\varepsilon$ = Redist.weno_eps).

TVD / SSP Runge–Kutta

Pseudo-time uses strong-stability-preserving RK. TVD-RK3 (Shu–Osher), with $\Delta\tau = \texttt{cfl}\cdot\Delta x$:

$$ \begin{aligned} \phi^{(1)} &= \phi^n + \Delta\tau\,L(\phi^n),\\ \phi^{(2)} &= \tfrac34\phi^n + \tfrac14\phi^{(1)} + \tfrac14\Delta\tau\,L(\phi^{(1)}),\\ \phi^{n+1} &= \tfrac13\phi^n + \tfrac23\phi^{(2)} + \tfrac23\Delta\tau\,L(\phi^{(2)}). \end{aligned} $$

Select with Redist.integrator = rk2 | rk3.

Russo–Smereka subcell constraint

Pins the zero contour by advancing interface-adjacent cells with a neighbour-independent distance,

$$ L_{ij} = -\frac{1}{\Delta x} \Bigl(\operatorname{sgn}(\phi^0_{ij})\,|\phi_{ij}| - D_{ij}\Bigr),\qquad D_{ij} = \frac{\phi^0_{ij}}{|\nabla\phi^0|_{ij}}. $$

Because $L_{ij}$ never reads a neighbour, integrating it cannot move the interface into an adjacent cell. Enable with Redist.subcell_fix = 1.

Sussman–Fatemi volume correction

Adds a per-cell Lagrange multiplier enforcing $\frac{d}{d\tau}\!\int H(\phi)\,dV = 0$ in the band. With $f = \delta_\epsilon(\phi_0)|\nabla\phi|$,

$$ \lambda_{ij} = \frac{-\sum_{\Omega_{ij}} \delta_\epsilon(\phi_0)\,L} {\sum_{\Omega_{ij}} \delta_\epsilon(\phi_0)^2\,|\nabla\phi|}, \qquad \phi_\tau = L(\phi) + \lambda_{ij} f_{ij}. $$

Enable with Redist.volume_fix = 1.

Geometric redistancing — `GeometricRedistance`

Computes distances directly from a reconstructed interface, in four stages. See GeometricRedistance.

Bounded-facet reconstruction. Each cut cell yields a segment (2D) / disc (3D) of half-extent $L$ through the crossing centroid, normal $\mathbf n = \nabla\phi/|\nabla\phi|$. For a point at signed perpendicular distance $p$ and in-plane distance $s$, $$ d_{\text{facet}} = \begin{cases} |p| & s\le L,\\ \sqrt{p^2+(s-L)^2} & s>L.\end{cases} $$ Bounding the facet stops an unbounded plane from chamfering convex corners.
First-band distances. Minimum bounded-facet distance over the $3^d$ neighbourhood, frozen as Dirichlet data.
Propagation. Parallel Rouy–Tourin (Jacobi) eikonal sweeps solve $|\nabla d| = 1$; the band is clamped and $\operatorname{sgn}(\phi_0)$ reapplied.
Band-local mass correction. A per-neighbourhood shift $$ \delta_{ij} = \frac{\sum_{\Omega_{ij}}(H_\epsilon(\phi_0)-H_\epsilon(\phi))} {\sum_{\Omega_{ij}} \delta_\epsilon(\phi)},\qquad \phi \mathrel{-}= \delta_{ij}, $$ restores lost area at corners without translating intact thin walls (so it does not close the Zalesak slot).

Advection & AMR

Coarse–fine consistency uses FillPatchTwoLevels and average_down. The companion AdvectionSolver adds Berger–Colella time subcycling, $\Delta t_\ell = \Delta t_{\ell-1}/r$, with time-interpolated coarse–fine boundary data and a volume-conserving coarse-band sync, exercised on the ZalesakProblem rotating slotted-disk / sphere benchmark.

`Redist.*` — PDE driver keys

Key	Default	Meaning
`scheme`	`godunov1`	`godunov1` \| `weno5`
`integrator`	`rk3`	`rk2` \| `rk3`
`subcell_fix`	`1`	Russo–Smereka interface constraint
`volume_fix`	`0`	Sussman–Fatemi volume conservation
`cfl`	`0.5`	$\Delta\tau = \texttt{cfl}\cdot\Delta x$
`sign_eps`	`1.0`	Sign-smoothing width (cells)
`weno_eps`	`1e-6`	WENO5 regularization
`vc_eps`, `vc_radius`	`1.5`, `1`	Volume-fix delta width / neighbourhood radius
`n_iter`	`5`	Default pseudo-time iterations

`Geom.*` — geometric driver keys

Key	Default	Meaning
`mass_preserving`	`1`	Band-local volume-conserving correction
`band`	`5.0` / `6.0` for Zalesak	Narrow-band half-width (cells)
`prop_iters`	`0`	Eikonal iterations (0 → auto)
`mass_eps`	`1.5`	Smoothed-Heaviside width (cells)
`vc_radius`	`2`	Local mass-correction neighbourhood (cells)

redist.method = levelset | geometric selects the driver in the standalone executables.

`inputs/inputs.regression`

Feeds the combined regression driver (all namespaces in one file):

# Shared grid (used by unit / volume / geometric / zalesak / sync).
geometry.prob_lo = 0.0 0.0 0.0
geometry.prob_hi = 1.0 1.0 1.0
amr.n_cell       = 64 64 64
amr.max_level    = 1
amr.ref_ratio    = 2 2

# PDE reinitialization driver.
Redist.scheme      = weno5
Redist.integrator  = rk3
Redist.subcell_fix = 1
Redist.volume_fix  = 0
Redist.cfl         = 0.5

# Convergence test (value-error metric; converges ~2).
conv.band_lo     = 0.0
conv.band_factor = 4.0
conv.min_order   = 1.5
conv.resolutions = 32 64 128

# Volume conservation (smoothed-Heaviside A/B).
vol.n_iter        = 40
vol.heaviside_eps = 1.5
vol.min_reduction = 1.5

# Geometric driver + geometric_test.
Geom.mass_preserving = 1
Geom.band            = 5.0
geom.grad_tol        = 8.0e-2
geom.vol_tol         = 1.0e-2

# Conservative-sync A/B.
sync.min_reduction = 1.3
sync.min_nosync    = 1.0e-7

`inputs/inputs.zdisk2d`

Standalone 2D Zalesak disk (namespace zd2d; redist.method chooses the algorithm):

amr.n_cell       = 128 128
amr.max_level    = 2
amr.ref_ratio    = 2 2

redist.method    = levelset       # levelset | geometric

# PDE and geometric driver blocks (whichever method is chosen reads its own).
Redist.scheme    = weno5
Redist.integrator= rk3
Geom.mass_preserving = 1
Geom.band            = 6.0

# Advection / rotation / coupling.
zd2d.adv_scheme   = weno5         # godunov | weno5
zd2d.revs         = 1.0
zd2d.omega        = 6.283185307179586    # 2*pi
zd2d.cfl          = 0.5
zd2d.reinit_every = 5
zd2d.regrid_every = 10

# Slotted-disk geometry.
zd2d.disk_cx = 0.5
zd2d.disk_cy = 0.75
zd2d.radius  = 0.15
zd2d.slot_width  = 0.05
zd2d.slot_bottom = 0.60
zd2d.slot_top    = 0.85

# Output.
zd2d.error_tol   = 0.15
zd2d.area_tol    = 0.05
zd2d.write_plot  = 1
zd2d.plot_int    = 0              # >0 -> per-step animation frames
zd2d.volume_file = volume.txt     # per-step enclosed-area log

`inputs/inputs.zsphere3d`

Standalone 3D Zalesak sphere (namespace zs3d):

amr.n_cell    = 64 64 64
amr.max_level = 1
amr.ref_ratio = 2

redist.method = geometric         # levelset | geometric

zs3d.adv_scheme = weno5
zs3d.revs       = 1.0
zs3d.cfl        = 0.5
zs3d.reinit_every = 5

# Slotted-sphere geometry (adds z-centre and slot z-extent).
zs3d.cx = 0.5
zs3d.cy = 0.75
zs3d.cz = 0.5
zs3d.radius      = 0.15
zs3d.slot_width  = 0.05
zs3d.slot_bottom = 0.60
zs3d.slot_top    = 0.85
zs3d.slot_zwidth = 0.30

zs3d.error_tol = 0.20
zs3d.area_tol  = 0.05
zs3d.write_plot= 1

C++ API

Minimal usage of each driver:

#include "LevelSetRedistance.H"    // PDE driver
#include "GeometricRedistance.H"   // geometric driver

amrex::Vector<amrex::MultiFab*> phi = /* per-level level-set pointers */;

// PDE (Sussman) redistancing
LevelSetRedistance ls(&amrcore, /*nghost=*/3);
ls.setScheme(RedistScheme::WENO5);
ls.setIntegrator(RedistIntegrator::RK3);
ls.setInterfaceConstraint(true);      // Russo–Smereka
ls.readParameters();                  // "Redist" inputs override
ls.redistance(phi, /*n_iter=*/5);

// Geometric (Ausas–Dari–Buscaglia) redistancing
GeometricRedistance geo(&amrcore, /*nghost=*/3);
geo.setMassPreserving(true);
geo.readParameters();                 // "Geom" inputs override
geo.redistance(phi, /*n_iter=*/0);    // 0 -> auto propagation iterations

Full class and method documentation, including every parameter and the kernel mathematics, is in the Doxygen API reference (LevelSetRedistance, GeometricRedistance, AdvectionSolver).

Tests & metrics

cd build/Test
./regression inputs.regression --output regression_results.xml   # full suite
./regression inputs.regression --only geometric_test             # one case
ctest --output-on-failure                                        # via CTest

Test	Verifies	Metric
`unit_test`	$\|\nabla\phi\|\to1$; interface preserved	perturbed circle/sphere
`convergence_test`	reinitialization order	value error vs exact distance (order ≈ 2)
`volume_test`	Sussman–Fatemi reduces drift	smoothed-Heaviside volume A/B
`geometric_test`	geometric driver accuracy + volume	band gradient + volume
`zalesak_test`	shape/area after one revolution	finest-data masked sym. difference
`zalesak_sync_test`	conservative sync reduces mismatch	smoothed coarse–fine A/B

Metrics deliberately use a smoothed Heaviside, not a sharp cell count, and a value-error convergence measure, so they resolve sub-cell interface motion and are invariant across 1–8 MPI ranks in 2D and 3D.

Python / Jupyter harness

zalesak_harness.py, a Jupytext-compatible script runs the standalone executables and provides static visualisation, time-series animation, level-set-vs-geometric comparison, parameter sweeps, and a mass-conservation study: volume, relative drift, and instantaneous drift rate.

pip install numpy matplotlib yt scikit-image jupytext   # + ffmpeg for MP4
jupytext --to notebook zalesak_harness.py

References

Sussman, Smereka & Osher (1994). J. Comput. Phys. 114(1), 146–159. doi:10.1006/jcph.1994.1155
Russo & Smereka (2000). J. Comput. Phys. 163(1), 51–67. doi:10.1006/jcph.2000.6553
Sussman & Fatemi (1999). SIAM J. Sci. Comput. 20(4), 1165–1191. doi:10.1137/S1064827596298245
Jiang & Peng (2000). SIAM J. Sci. Comput. 21(6), 2126–2143. doi:10.1137/S106482759732455X
Ausas, Dari & Buscaglia (2011). Int. J. Numer. Methods Fluids 65(8), 989–1010. doi:10.1002/fld.2227
Zalesak (1979). J. Comput. Phys. 31(3), 335–362. doi:10.1016/0021-9991(79)90051-2
Zhang et al. (2019). J. Open Source Softw. 4(37), 1370. doi:10.21105/joss.01370

The full bibliography and BibTeX are in the project README and the API reference.