Level Set

The LevelSet module implements a fire spread simulator using the level set method. A signed distance function φ(x, y, t) implicitly tracks the fire front:

Region Condition Meaning
Burned φ < 0 Fire has passed
Front φ = 0 Active fire perimeter
Unburned φ > 0 Fuel ahead of fire

The front evolves according to the Hamilton-Jacobi equation:

\[\frac{\partial \phi}{\partial t} + F |\nabla \phi| = 0\]

where \(F(x, y, t) \geq 0\) is the local spread rate (e.g. from the Rothermel model).

LevelSetGrid

grid = LevelSetGrid(100, 100, dx=30.0)
LevelSetGrid{Float64} 100×100 (t=0.0, burned=0/10000)

Boundary Conditions

The bc keyword controls how the solver handles grid edges. Three built-in options are available:

Type Behavior Physical Interpretation
ZeroNeumann() Zero-gradient at edges (default). Missing neighbors contribute zero to finite differences. Open boundaries — the fire front can reach and exit the domain naturally. Most common choice for realistic simulations where the domain is a subset of a larger landscape.
Dirichlet() Fixed-value edges. Edge cells are never updated by advance! or reinitialize!. Firebreak / barrier — edge cells retain their initial values (unburned). Useful for modeling impassable boundaries like rivers, roads, or cleared fire lines.
Periodic() Wrap-around. Edges reference values from the opposite side of the grid. Tiling / infinite landscape — fire exiting one edge enters from the opposite side. Useful for studying fire behavior in a homogeneous landscape without edge effects. 
grid_d = LevelSetGrid(50, 50, dx=30.0, bc=Dirichlet())
grid_p = LevelSetGrid(50, 50, dx=30.0, bc=Periodic())
println("Dirichlet: ", grid_d)
println("Periodic:  ", grid_p)
Dirichlet: LevelSetGrid{Float64}(50×50)
Periodic:  LevelSetGrid{Float64}(50×50)

Custom Boundary Conditions

Define a custom BC by subtyping AbstractBoundaryCondition and implementing the four one-sided finite-difference helpers. Optionally implement _skip_update to exclude certain cells from updates.

using Wildfires.LevelSet

struct MyBC <: AbstractBoundaryCondition end

# One-sided finite differences at boundaries:
#   _Dxm(φ, i, j, dx, bc)        → backward x-difference
#   _Dxp(φ, i, j, nx, dx, bc)    → forward x-difference
#   _Dym(φ, i, j, dy, bc)        → backward y-difference
#   _Dyp(φ, i, j, ny, dy, bc)    → forward y-difference

# Optional: skip updating certain cells (default: false)
#   _skip_update(i, j, ny, nx, bc) → Bool

See the built-in implementations in src/LevelSet.jl for reference.

Coordinates

xs = xcoords(grid)
println("x: $(first(xs)) to $(last(xs)) m  ($(length(xs)) cells)")

ys = ycoords(grid)
println("y: $(first(ys)) to $(last(ys)) m  ($(length(ys)) cells)")
x: 15.0 to 2985.0 m  (100 cells)
y: 15.0 to 2985.0 m  (100 cells)

Ignition

ignite!(grid, 1500.0, 1500.0, 100.0)
grid
LevelSetGrid{Float64} 100×100 (t=0.0, burned=32/10000, ignited=32)

Queries

burn_area(grid)  # m²
28800.0
sum(burned(grid))  # number of burned cells
32

Advancing the Front

The spread rate field F has the same dimensions as the grid. Each cell specifies the local spread rate in m/min.

F = fill(10.0, size(grid))  # uniform 10 m/min everywhere

for _ in 1:20
    advance!(grid, F, 0.5)
end

grid
LevelSetGrid{Float64} 100×100 (t=10.0, burned=264/10000, ignited=264)

Solvers

The advance! function accepts an optional solver argument that controls the numerical scheme for evolving the level set equation. Three solvers are available:

Solver Order Time Stepping Stencil Use Case
Godunov() 1st Forward Euler 1-cell Default — robust, simple, good for most cases
Superbee() 2nd RK2 (Heun) 2-cell Sharper fronts, less numerical diffusion
WENO5() 5th SSP-RK3 3-cell Highest accuracy, sharpest front resolution 
# Explicit solver selection (Godunov is the default)
advance!(grid, F, 0.5, Godunov())

# Second-order Superbee
grid2 = LevelSetGrid(100, 100, dx=30.0)
ignite!(grid2, 1500.0, 1500.0, 100.0)
F2 = fill(10.0, size(grid2))
advance!(grid2, F2, 0.5, Superbee())
grid2
LevelSetGrid{Float64} 100×100 (t=0.5, burned=52/10000, ignited=52)

When no solver is specified, advance!(grid, F, dt) defaults to Godunov() for backward compatibility.

Godunov() — First-Order Upwind

The default solver computes \(|\nabla\phi|\) using Godunov’s upwind scheme, a monotone numerical flux for Hamilton-Jacobi equations. The idea is to use one-sided (upwind) finite differences that respect the direction information is traveling — toward the fire front.

At each cell, four one-sided differences are computed:

\[D_x^- = \frac{\phi_{i,j} - \phi_{i,j-1}}{\Delta x}, \quad D_x^+ = \frac{\phi_{i,j+1} - \phi_{i,j}}{\Delta x}\]

\[D_y^- = \frac{\phi_{i,j} - \phi_{i-1,j}}{\Delta y}, \quad D_y^+ = \frac{\phi_{i+1,j} - \phi_{i,j}}{\Delta y}\]

Since the spread rate \(F \geq 0\), the front moves in the direction of \(\nabla\phi\) (outward from burned regions). Godunov’s scheme selects the appropriate one-sided difference in each coordinate:

\[\hat{D}_x = \max\!\bigl(\max(D_x^-,\, 0),\; -\min(D_x^+,\, 0)\bigr)\]

\[\hat{D}_y = \max\!\bigl(\max(D_y^-,\, 0),\; -\min(D_y^+,\, 0)\bigr)\]

Then \(|\nabla\phi| \approx \sqrt{\hat{D}_x^2 + \hat{D}_y^2}\) and the update is:

\[\phi^{n+1}_{i,j} = \phi^n_{i,j} - \Delta t \; F_{i,j} \; \sqrt{\hat{D}_x^2 + \hat{D}_y^2}\]

This selection rule ensures that only upstream information influences the update, preventing oscillations that would occur with central differences. At domain boundaries, the one-sided differences that would require out-of-bounds data are set to zero (implicit zero-Neumann / reflective boundary condition).

Superbee() — Second-Order TVD

The Superbee solver converts the Hamilton-Jacobi equation into advection form and applies a Superbee flux limiter with RK2 (Heun’s method) time stepping, following the approach used in ELMFIRE. This produces sharper fire fronts with less numerical diffusion than the first-order Godunov scheme.

How it works:

  1. Compute the front normal \(\hat{n} = \nabla\phi / |\nabla\phi|\) via central differences
  2. Set advection velocities \(u_x = F \cdot n_x\), \(u_y = F \cdot n_y\)
  3. Apply the Superbee flux limiter to compute limited gradients \(\partial\phi/\partial x\), \(\partial\phi/\partial y\)
  4. Advance with RK2 (Heun’s method):
    • Stage 1: Forward Euler \(\phi^* = \phi^n - \Delta t \, (u_x \, \partial\phi/\partial x + u_y \, \partial\phi/\partial y)\)
    • Stage 2: Recompute from \(\phi^*\), then average: \(\phi^{n+1} = \tfrac{1}{2}(\phi^n + \phi^{**})\)

The Superbee limiter is total variation diminishing (TVD), preventing spurious oscillations while maintaining second-order accuracy at smooth regions. It uses a 2-cell stencil (accessing φ[i±2, j±2]), wider than Godunov’s 1-cell stencil.

Parameters:

Superbee(phi_clamp=100.0, grad_clamp=1000.0)
Parameter Default Description
phi_clamp 100.0 Clamp φ values to [-phi_clamp, phi_clamp] after each RK stage
grad_clamp 1000.0 Clamp flux-limited gradients to prevent runaway values

These guards prevent numerical blowup when the level set field develops steep gradients or large values far from the front.

WENO5() — Fifth-Order WENO

The WENO5 solver uses the Jiang-Shu WENO5 (Weighted Essentially Non-Oscillatory) reconstruction for spatial derivatives with SSP-RK3 (Strong Stability Preserving Runge-Kutta) time stepping. This provides 5th-order accuracy in smooth regions while maintaining sharp resolution near discontinuities.

How it works:

  1. Compute left-biased (\(D^-\)) and right-biased (\(D^+\)) WENO5 derivatives in each direction using a 6-point stencil
  2. Each WENO5 derivative blends three candidate 3rd-order stencils weighted by smoothness indicators — smooth regions get near-optimal weights, while oscillatory regions suppress offending stencils
  3. Apply the Godunov numerical Hamiltonian (same as first-order) to combine derivatives
  4. Advance with SSP-RK3 (Shu-Osher form):
    • Stage 1: \(\phi_1 = \phi^n + \Delta t \cdot L(\phi^n)\)
    • Stage 2: \(\phi_2 = \tfrac{3}{4}\phi^n + \tfrac{1}{4}(\phi_1 + \Delta t \cdot L(\phi_1))\)
    • Stage 3: \(\phi^{n+1} = \tfrac{1}{3}\phi^n + \tfrac{2}{3}(\phi_2 + \Delta t \cdot L(\phi_2))\)

The 3-cell stencil (accessing φ[i±3, j±3]) is the widest of the three solvers.

Parameters:

WENO5(phi_clamp=100.0)
Parameter Default Description
phi_clamp 100.0 Clamp φ values to [-phi_clamp, phi_clamp] after each RK stage

Comparing Solvers

using Wildfires.SpreadModel
using Wildfires.Rothermel: SHORT_GRASS, FuelClasses

M = FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)
model = FireSpreadModel(SHORT_GRASS, UniformWind(speed=8.0), UniformMoisture(M), FlatTerrain())

fig = Figure(size=(900, 300))
for (col, (label, solver)) in enumerate([
    ("Godunov (1st order)", Godunov()),
    ("Superbee (2nd order TVD)", Superbee()),
    ("WENO5 (5th order)", WENO5()),
])
    g = LevelSetGrid(200, 200, dx=30.0)
    ignite!(g, 3000.0, 3000.0, 50.0)
    simulate!(g, model, steps=100, dt=0.5, solver=solver)
    ax = Axis(fig[1, col], title=label, aspect=DataAspect())
    fireplot!(ax, g)
    hidedecorations!(ax)
end
fig
TipWhen to use which solver
  • Godunov() is the safe default — it’s simple, robust, and sufficient for most simulations.
  • Superbee() is worth trying when you need sharper fire front resolution, especially for high-resolution grids or when comparing against observed fire perimeters. It requires roughly 2x the computation per time step (two RK stages + wider stencil).
  • WENO5() provides the highest accuracy and sharpest fronts, but at roughly 3x the cost of Godunov (three RK stages + wider stencil). Best for high-resolution research simulations.

CFL Condition

The time step dt must satisfy the Courant-Friedrichs-Lewy (CFL) condition for numerical stability:

\[\text{dt} \leq \text{cfl} \cdot \frac{\min(\text{dx}, \text{dy})}{\max(F)}\]

where cfl is a safety coefficient (default 0.5). If dt is too large, the upwind scheme becomes unstable and φ values will diverge.

The cfl_dt function computes a stable time step:

F = fill(10.0, size(grid))
cfl_dt(grid, F)           # cfl=0.5 by default
1.5
cfl_dt(grid, F; cfl=1.0)  # less conservative
3.0

The simulate! function uses automatic CFL-limited time stepping by default. Pass an explicit dt to override:

simulate!(grid, model, steps=100)                          # auto CFL (recommended)
simulate!(grid, model, steps=100, dt=0.5)                  # fixed dt
simulate!(grid, model, steps=100, solver=Superbee())       # 2nd-order solver
simulate!(grid, model, steps=100, solver=WENO5())          # 5th-order solver

Reinitialization

Over many time steps, φ can drift away from a true signed distance function (gradients become too flat or too steep). Periodic reinitialization restores |∇φ| ≈ 1, which keeps the upwind scheme accurate.

reinitialize!(grid)
grid
LevelSetGrid{Float64} 100×100 (t=10.5, burned=280/10000, ignited=280)

Two reinitialization methods are available:

Method Accuracy Cost Use Case
IterativeReinit(iterations=5) Approximate Low Default — fast, good enough for most simulations
NewtonReinit(upsample=8) Machine-precision Higher Research — exact signed distance restoration

NewtonReinit — Newton Closest-Point

The Newton method produces machine-precision signed distance fields by:

  1. Building a cubic interpolant of φ using Interpolations.jl
  2. Locating zero-contour sample points via Newton projection
  3. Building a KD-tree for efficient nearest-neighbor queries
  4. Computing exact signed distance from each grid point to the nearest interface point
reinitialize!(grid, NewtonReinit())
reinitialize!(grid, NewtonReinit(upsample=16, maxiters=30))

In simulate!, pass the reinit keyword:

simulate!(grid, model, steps=100, reinit=NewtonReinit())

Curvature Regularization

The curvature keyword in simulate! adds curvature-dependent speed regularization:

\[F_{\text{eff}} = \max(F - b \cdot \kappa, 0)\]

where \(b\) is the curvature coefficient and \(\kappa\) is the mean curvature of the level set. This smooths the fire front by reducing speed at convex bulges and increasing it at concave indentations.

simulate!(grid, model, steps=100, curvature=5.0)

The mean curvature is computed via central differences:

\[\kappa = \frac{\phi_{xx}\phi_y^2 - 2\phi_x\phi_y\phi_{xy} + \phi_{yy}\phi_x^2}{(\phi_x^2 + \phi_y^2)^{3/2}}\]

When curvature > 0, the CFL time step also enforces the parabolic constraint \(\Delta t \leq \Delta x^2 / (2b)\). For typical fire simulation parameters, this is never the bottleneck.

Plotting

Standard Makie Plots

The LevelSetGrid integrates with Makie’s heatmap, contour, contourf, and surface via convert_arguments:

fig = Figure(size=(700, 300))
ax1 = Axis(fig[1, 1], title="heatmap", aspect=DataAspect())
heatmap!(ax1, grid, colormap=:RdYlGn)

ax2 = Axis(fig[1, 2], title="contourf", aspect=DataAspect())
contourf!(ax2, grid, colormap=:RdYlGn)
fig

fireplot

The fireplot recipe combines a heatmap of φ with a contour line at the fire front:

fig = Figure()
ax = Axis(fig[1, 1], title="t = $(grid.t) min", aspect=DataAspect(),
    xlabel="x (m)", ylabel="y (m)")
fireplot!(ax, grid)
fig

Full Example

Simulate fire spread with a spatially varying spread rate and visualize snapshots:

# Fresh grid with ignition
g = LevelSetGrid(200, 200, dx=15.0)
ignite!(g, 1500.0, 1500.0, 50.0)

# Spread rate: faster on the right side
F = [5.0 + 10.0 * (j / 200)
    for i in 1:200, j in 1:200]

fig = Figure(size=(700, 250))
times = [0, 50, 100]
for (col, target_t) in enumerate(times)
    while g.t < target_t
        advance!(g, F, 0.25)
    end
    col > 1 && reinitialize!(g)
    t = round(g.t, digits=1)
    ax = Axis(fig[1, col],
        title="t = $t min",
        aspect=DataAspect())
    fireplot!(ax, g)
    hidedecorations!(ax)
end
fig

References

  • Osher, S. & Sethian, J.A. (1988). Fronts propagating with curvature-dependent speed. J. Computational Physics, 79(1), 12–49.
  • Mallet, V., Keyes, D.E., & Fendell, F.E. (2009). Modeling wildland fire propagation with level set methods. Computers & Mathematics with Applications.
  • Lautenberger, C. (2013). Wildland fire modeling with an Eulerian level set method and automated calibration. Fire Safety Journal, 62, 289–298. (ELMFIRE — Superbee flux limiter approach)
  • Jiang, G.-S. & Peng, D. (2000). Weighted ENO schemes for Hamilton-Jacobi equations. SIAM J. Scientific Computing, 21(6), 2126–2143. (WENO5 for level set methods)