Simulation

Fire state and simulation functions

This module provides the core fire simulation state and execution functions.

using Elmfire
using Plots

Types

FireState

Complete state of a fire simulation, including the level set field, burned area, and fire behavior outputs.

mutable struct FireState{T<:AbstractFloat}
    # Level set field (with 2-cell padding for stencil operations)
    phi::Matrix{T}
    phi_old::Matrix{T}

    # Output fields (no padding, matches grid dimensions)
    time_of_arrival::Matrix{T}    # -1 if not burned
    burned::BitMatrix             # True if cell has burned
    spread_rate::Matrix{T}        # Final spread rate (ft/min)
    fireline_intensity::Matrix{T} # Fireline intensity (kW/m)
    flame_length::Matrix{T}       # Flame length (ft)

    # Velocity components (with padding)
    ux::Matrix{T}
    uy::Matrix{T}

    # Narrow band tracking
    narrow_band::NarrowBand

    # Grid parameters
    ncols::Int           # Number of columns (without padding)
    nrows::Int           # Number of rows (without padding)
    cellsize::T          # Cell size (ft)
    xllcorner::T         # X coordinate of lower-left corner
    yllcorner::T         # Y coordinate of lower-left corner
    padding::Int         # Padding for stencil operations
end

Constructor:

FireState{T}(ncols, nrows, cellsize;
    xllcorner = zero(T),
    yllcorner = zero(T),
    padding = 2,
    band_thickness = 5
) -> FireState{T}

# Create a 100x100 grid with 30ft cells
state = FireState{Float64}(100, 100, 30.0)

println("Grid size: $(state.ncols) x $(state.nrows)")
println("Cell size: $(state.cellsize) ft")
println("Domain: $(state.ncols * state.cellsize) x $(state.nrows * state.cellsize) ft")

Grid size: 100 x 100
Cell size: 30.0 ft
Domain: 3000.0 x 3000.0 ft

SimulationConfig

Configuration for full simulation with crown fire, spotting, and weather interpolation.

struct SimulationConfig{T<:AbstractFloat}
    enable_crown_fire::Bool        # Enable crown fire modeling
    enable_spotting::Bool          # Enable ember transport
    crown_fire_adj::T              # Crown fire adjustment factor
    critical_canopy_cover::T       # Minimum CC for active crown fire
    foliar_moisture::T             # Foliar moisture content (%)
    spotting_params::Union{Nothing, SpottingParameters{T}}
    use_sardoy::Bool               # Use Sardoy model for spotting
end

Constructor:

SimulationConfig{T}(;
    enable_crown_fire = false,
    enable_spotting = false,
    crown_fire_adj = 1.0,
    critical_canopy_cover = 0.4,
    foliar_moisture = 100.0,
    spotting_params = nothing,
    use_sardoy = false
) -> SimulationConfig{T}

# Basic configuration (surface fire only)
config_basic = SimulationConfig{Float64}()

# Crown fire enabled
config_crown = SimulationConfig{Float64}(
    enable_crown_fire = true,
    foliar_moisture = 100.0
)

# Full configuration with spotting
spotting = SpottingParameters{Float64}(
    mean_distance = 100.0,
    normalized_variance = 0.5,
    ws_exponent = 1.0,
    flin_exponent = 0.5,
    nembers_max = 10,
    surface_spotting_percent = 1.0,
    crown_spotting_percent = 10.0,
    pign = 50.0,
    min_distance = 10.0,
    max_distance = 2000.0
)

config_full = SimulationConfig{Float64}(
    enable_crown_fire = true,
    enable_spotting = true,
    spotting_params = spotting
)

SimulationConfig{Float64}(true, true, 1.0, 0.4, 100.0, SpottingParameters{Float64}(100.0, 0.5, 1.0, 0.5, 10, 1.0, 10.0, 50.0, 10.0, 2000.0), false)

CanopyGrid

Grid of canopy properties for crown fire modeling.

struct CanopyGrid{T<:AbstractFloat}
    cbd::Matrix{T}      # Canopy bulk density (kg/m³)
    cbh::Matrix{T}      # Canopy base height (m)
    cc::Matrix{T}       # Canopy cover (fraction 0-1)
    ch::Matrix{T}       # Canopy height (m)
end

# Create uniform canopy
canopy = CanopyGrid{Float64}(
    fill(0.15, 100, 100),  # CBD
    fill(2.0, 100, 100),   # CBH
    fill(0.6, 100, 100),   # CC
    fill(20.0, 100, 100)   # CH
)

println("Canopy bulk density: $(canopy.cbd[1,1]) kg/m³")
println("Canopy base height: $(canopy.cbh[1,1]) m")
println("Canopy cover: $(canopy.cc[1,1] * 100)%")
println("Canopy height: $(canopy.ch[1,1]) m")

Canopy bulk density: 0.15 kg/m³
Canopy base height: 2.0 m
Canopy cover: 60.0%
Canopy height: 20.0 m

Simulation Functions

simulate!

Run fire simulation with spatially-varying conditions.

simulate!(
    state::FireState{T},
    fuel_ids::AbstractMatrix{Int},
    fuel_table::FuelModelTable{T},
    weather::ConstantWeather{T},
    slope::AbstractMatrix{T},
    aspect::AbstractMatrix{T},
    t_start::T,
    t_stop::T;
    dt_initial = 1.0,
    target_cfl = 0.45,
    dt_max = 10.0,
    spread_rate_adj = 1.0,
    accel_time_constant = 0.0,
    diurnal = nothing,
    dampening = SpreadRateDampeningConfig{T}(),
    callback = nothing
)

Argument	Description
state	Fire state (modified in place)
fuel_ids	Matrix of fuel model IDs
fuel_table	Fuel model lookup table
weather	Constant weather conditions
slope	Slope in degrees for each cell
aspect	Aspect direction in degrees
t_start	Start time (minutes)
t_stop	Stop time (minutes)
accel_time_constant	Fire acceleration time constant (minutes). 0 = no acceleration ramp.
diurnal	DiurnalConfig for day/night spread rate modulation, or nothing
dampening	SpreadRateDampeningConfig for spread rate dampening

fuel_table = create_standard_fuel_table(Float64)
weather = ConstantWeather{Float64}(
    wind_speed_mph = 10.0,
    wind_direction = 270.0,
    M1 = 0.06, M10 = 0.08, M100 = 0.10,
    MLH = 0.60, MLW = 0.90
)

state = FireState{Float64}(100, 100, 30.0)
fuel_ids = fill(1, 100, 100)
slope = zeros(Float64, 100, 100)
aspect = zeros(Float64, 100, 100)

ignite!(state, 50, 50, 0.0)
simulate!(state, fuel_ids, fuel_table, weather, slope, aspect, 0.0, 30.0)

println("Burned area: $(round(get_burned_area_acres(state), digits=2)) acres")

heatmap(state.burned', title = "Burned Area", color = :YlOrRd, aspect_ratio = 1)

Burned area: 5.19 acres

simulate_uniform!

Run simulation with uniform fuel, slope, and aspect across the entire domain.

simulate_uniform!(
    state::FireState{T},
    fuel_id::Int,
    fuel_table::FuelModelTable{T},
    weather::ConstantWeather{T},
    slope_deg::T,
    aspect_deg::T,
    t_start::T,
    t_stop::T;
    kwargs...
)

state = FireState{Float64}(100, 100, 30.0)
ignite!(state, 50, 50, 0.0)

# Uniform fuel model 1 (short grass), 10° slope facing east
simulate_uniform!(state, 1, fuel_table, weather, 10.0, 90.0, 0.0, 30.0)

println("Burned area: $(round(get_burned_area_acres(state), digits=2)) acres")

Burned area: 7.09 acres

simulate_full!

Run full simulation with crown fire, spotting, and weather interpolation support.

simulate_full!(
    state::FireState{T},
    fuel_ids::AbstractMatrix{Int},
    fuel_table::FuelModelTable{T},
    weather_interp::WeatherInterpolator{T},
    slope::AbstractMatrix{T},
    aspect::AbstractMatrix{T},
    t_start::T,
    t_stop::T;
    canopy = nothing,
    config = SimulationConfig{T}(),
    dt_initial = 1.0,
    target_cfl = 0.45,
    dt_max = 10.0,
    spread_rate_adj = 1.0,
    accel_time_constant = 0.0,
    diurnal = nothing,
    dampening = SpreadRateDampeningConfig{T}(),
    callback = nothing,
    rng = Random.default_rng()
) -> SpotFireTracker

When a CanopyGrid is provided, simulate_full! uses the sheltered wind adjustment factor (3-argument wind_adjustment_factor) to account for canopy cover and height, reducing mid-flame wind speed under forest canopy.

using Random

state = FireState{Float64}(100, 100, 30.0)
fuel_ids = fill(10, 100, 100)  # Timber understory
slope = zeros(Float64, 100, 100)
aspect = zeros(Float64, 100, 100)

# Create canopy
canopy = CanopyGrid{Float64}(
    fill(0.15, 100, 100),
    fill(2.0, 100, 100),
    fill(0.6, 100, 100),
    fill(20.0, 100, 100)
)

# Configuration with crown fire
config = SimulationConfig{Float64}(
    enable_crown_fire = true,
    foliar_moisture = 100.0
)

weather_interp = create_constant_interpolator(weather, 100, 100, 30.0)

ignite!(state, 50, 50, 0.0)
simulate_full!(state, fuel_ids, fuel_table, weather_interp, slope, aspect,
    0.0, 30.0; canopy = canopy, config = config)

println("Burned area: $(round(get_burned_area_acres(state), digits=2)) acres")
println("Max fireline intensity: $(round(maximum(state.fireline_intensity), digits=0)) kW/m")

Burned area: 0.17 acres
Max fireline intensity: 204.0 kW/m

Ignition Functions

ignite!

Ignite a cell at grid coordinates.

ignite!(state::FireState{T}, ix::Int, iy::Int, t::T)

Sets up the level set as a signed distance function near the ignition point.

state = FireState{Float64}(50, 50, 30.0)
ignite!(state, 25, 25, 0.0)

println("Cell (25,25) burned: $(state.burned[25, 25])")
println("Time of arrival: $(state.time_of_arrival[25, 25]) min")

Cell (25,25) burned: true
Time of arrival: 0.0 min

ignite_point!

Ignite a cell at world coordinates.

ignite_point!(state::FireState{T}, x::T, y::T, t::T)

Converts world coordinates to grid indices and ignites the corresponding cell.

ignite_circle!

Ignite all cells within a circular area.

ignite_circle!(state::FireState{T}, center_x::Int, center_y::Int, radius_cells::T, t::T)

state = FireState{Float64}(50, 50, 30.0)
ignite_circle!(state, 25, 25, 3.0, 0.0)  # 3-cell radius

println("Burned cells: $(count(state.burned))")

Burned cells: 29

State Management

reset!

Reset a FireState to initial conditions for reuse in ensemble simulations.

reset!(state::FireState{T})

state = FireState{Float64}(50, 50, 30.0)
ignite!(state, 25, 25, 0.0)
simulate_uniform!(state, 1, fuel_table, weather, 0.0, 0.0, 0.0, 10.0)

println("Before reset: $(count(state.burned)) burned cells")
reset!(state)
println("After reset: $(count(state.burned)) burned cells")

Before reset: 52 burned cells
After reset: 0 burned cells

copy

Create a deep copy of a FireState for thread-safe parallel execution.

Base.copy(state::FireState{T}) -> FireState{T}

Output Functions

get_fire_perimeter

Get the grid coordinates of cells on the fire perimeter.

get_fire_perimeter(state::FireState) -> Vector{Tuple{Int, Int}}

Returns burned cells adjacent to unburned cells.

state = FireState{Float64}(50, 50, 30.0)
ignite!(state, 25, 25, 0.0)
simulate_uniform!(state, 1, fuel_table, weather, 0.0, 0.0, 0.0, 15.0)

perimeter = get_fire_perimeter(state)
println("Perimeter cells: $(length(perimeter))")

Perimeter cells: 43

get_burned_area

Get the total burned area in square feet.

get_burned_area(state::FireState{T}) -> T

get_burned_area_acres

Get the total burned area in acres.

get_burned_area_acres(state::FireState{T}) -> T

state = FireState{Float64}(100, 100, 30.0)
ignite!(state, 50, 50, 0.0)
simulate_uniform!(state, 1, fuel_table, weather, 0.0, 0.0, 0.0, 30.0)

println("Burned area: $(round(get_burned_area(state), digits=0)) ft²")
println("Burned area: $(round(get_burned_area_acres(state), digits=2)) acres")

Burned area: 225900.0 ft²
Burned area: 5.19 acres

Helper Functions

wind_adjustment_factor

Calculate wind adjustment factor (WAF) from 20-ft wind to mid-flame wind.

Unsheltered (no canopy):

wind_adjustment_factor(fuel_bed_depth::T) -> T

Sheltered (under canopy):

wind_adjustment_factor(fuel_bed_depth::T, canopy_cover::T, canopy_height::T) -> T

Uses the Albini & Baughman (1979) sheltered formula when canopy cover and height are non-trivial. The 3-argument form is used automatically by simulate_full! when a CanopyGrid is provided.

depths = [0.5, 1.0, 2.0, 3.0, 4.0]
wafs = [wind_adjustment_factor(d) for d in depths]
wafs_sheltered = [wind_adjustment_factor(d, 0.6, 20.0) for d in depths]

for (d, w, ws) in zip(depths, wafs, wafs_sheltered)
    println("Depth $d ft: WAF = $(round(w, digits=3)) (unsheltered), $(round(ws, digits=3)) (sheltered)")
end

Depth 0.5 ft: WAF = 0.319 (unsheltered), 0.163 (sheltered)
Depth 1.0 ft: WAF = 0.362 (unsheltered), 0.163 (sheltered)
Depth 2.0 ft: WAF = 0.418 (unsheltered), 0.163 (sheltered)
Depth 3.0 ft: WAF = 0.459 (unsheltered), 0.163 (sheltered)
Depth 4.0 ft: WAF = 0.492 (unsheltered), 0.163 (sheltered)

calculate_tanslp2

Calculate tan²(slope) from slope in degrees.

calculate_tanslp2(slope_degrees::T) -> T

slopes = [0, 10, 20, 30, 45]
for s in slopes
    println("Slope $(s)°: tan²(slope) = $(round(calculate_tanslp2(Float64(s)), digits=4))")
end

Slope 0°: tan²(slope) = 0.0
Slope 10°: tan²(slope) = 0.0311
Slope 20°: tan²(slope) = 0.1325
Slope 30°: tan²(slope) = 0.3333
Slope 45°: tan²(slope) = 1.0

Coordinate Conversion

grid_to_padded(state::FireState, ix::Int, iy::Int) -> Tuple{Int, Int}
padded_to_grid(state::FireState, px::Int, py::Int) -> Tuple{Int, Int}

Convert between grid coordinates (1:ncols, 1:nrows) and padded array coordinates.

Fire Acceleration

acceleration_factor

Compute the fire acceleration factor 1 - exp(-t/τ), modeling the buildup from ignition to steady-state spread rate. When τ ≤ 0.5 or t/τ > 7, returns 1.0 (full spread).

acceleration_factor(t::T, tau::T) -> T

Pass accel_time_constant to simulate! or simulate_full! to apply this automatically.

tau = 10.0  # 10-minute time constant
times = 0.0:0.5:60.0
factors = [acceleration_factor(t, tau) for t in times]

plot(times, factors,
    xlabel = "Time since ignition (min)",
    ylabel = "Acceleration Factor",
    title = "Fire Acceleration (τ = $tau min)",
    linewidth = 2,
    legend = false,
    ylims = (0, 1.1)
)

# Effect on burned area
state1 = FireState{Float64}(100, 100, 30.0)
ignite!(state1, 50, 50, 0.0)
simulate_uniform!(state1, 1, fuel_table, weather, 0.0, 0.0, 0.0, 30.0)

state2 = FireState{Float64}(100, 100, 30.0)
ignite!(state2, 50, 50, 0.0)
simulate_uniform!(state2, 1, fuel_table, weather, 0.0, 0.0, 0.0, 30.0;
    accel_time_constant = 10.0)

println("No acceleration: $(round(get_burned_area_acres(state1), digits=2)) acres")
println("With τ=10 min:   $(round(get_burned_area_acres(state2), digits=2)) acres")

No acceleration: 5.19 acres
With τ=10 min:   2.36 acres

Diurnal Adjustment

DiurnalConfig

Configuration for day/night spread rate modulation. During the burn period (daytime), spread rate is unmodified. Outside the burn period, it is multiplied by overnight_factor (default 0.1).

struct DiurnalConfig{T<:AbstractFloat}
    forecast_start_hour::T          # Hour of day when simulation starts (default: 12)
    sunrise_hour::T                 # Hour of sunrise (default: 6)
    sunset_hour::T                  # Hour of sunset (default: 20)
    burn_period_center_frac::T      # Fraction of daylight for burn center (default: 0.667)
    burn_period_length::T           # Length of burn period in hours (default: 10)
    overnight_factor::T             # Spread rate multiplier at night (default: 0.1)
end

diurnal_adjustment

Returns the spread rate multiplier at simulation time t (minutes).

diurnal_adjustment(config::DiurnalConfig{T}, t::T) -> T

config = DiurnalConfig{Float64}()

times = 0.0:5.0:1440.0  # 24 hours in minutes
factors = [diurnal_adjustment(config, t) for t in times]

plot(times ./ 60, factors,
    xlabel = "Hours from simulation start",
    ylabel = "Diurnal Factor",
    title = "Diurnal Spread Rate Adjustment",
    linewidth = 2,
    legend = false,
    ylims = (0, 1.2)
)

Callbacks

The simulation functions accept an optional callback that is called after each timestep:

callback(state, t, dt, iteration)

state = FireState{Float64}(50, 50, 30.0)
ignite!(state, 25, 25, 0.0)

areas = Float64[]
times = Float64[]

function track_area(state, t, dt, iter)
    push!(times, t)
    push!(areas, get_burned_area_acres(state))
end

simulate_uniform!(state, 1, fuel_table, weather, 0.0, 0.0, 0.0, 20.0;
    callback = track_area)

plot(times, areas,
    xlabel = "Time (min)",
    ylabel = "Burned Area (acres)",
    title = "Fire Growth Over Time",
    linewidth = 2,
    legend = false
)