Heat equation — explicit Euler

//! 2-D heat equation — explicit forward-Euler stencil.
//!
//!   ∂u/∂t = α ∇² u
//!
//! Discretised on a w × h grid with periodic boundary. One step per call.
//! Stable when `α · dt · (1/dx² + 1/dy²) ≤ 0.5`. V-class L0-closed.

use serde::{Deserialize, Serialize};

#[derive(Debug, Clone, Copy, Serialize, Deserialize)]
pub struct HeatParams {
    pub alpha: f32,
    pub dt: f32,
    pub dx: f32,
}

impl Default for HeatParams {
    fn default() -> Self {
        Self { alpha: 0.25, dt: 0.2, dx: 1.0 }
    }
}

#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct HeatKernelReport {
    pub grid_w: u32,
    pub grid_h: u32,
    pub bit_ops: u64,
}

/// One explicit forward-Euler step on a periodic w × h grid.
/// Reads from `u`, writes to `u_next`. Returns the kernel report.
pub fn heat_step(u: &[f32], u_next: &mut [f32], w: u32, h: u32, p: &HeatParams) -> HeatKernelReport {
    let w_us = w as usize;
    let h_us = h as usize;
    debug_assert_eq!(u.len(), w_us * h_us);
    debug_assert_eq!(u_next.len(), w_us * h_us);
    let k = p.alpha * p.dt / (p.dx * p.dx);
    for y in 0..h_us {
        let yp = if y == 0 { h_us - 1 } else { y - 1 };
        let yn = if y + 1 == h_us { 0 } else { y + 1 };
        for x in 0..w_us {
            let xp = if x == 0 { w_us - 1 } else { x - 1 };
            let xn = if x + 1 == w_us { 0 } else { x + 1 };
            let c = u[y * w_us + x];
            let lap = u[y * w_us + xp]
                + u[y * w_us + xn]
                + u[yp * w_us + x]
                + u[yn * w_us + x]
                - 4.0 * c;
            u_next[y * w_us + x] = c + k * lap;
        }
    }
    let cells = (w_us * h_us) as u64;
    HeatKernelReport {
        grid_w: w,
        grid_h: h,
        bit_ops: cells * 8 * 32,
    }
}

/// Total energy (∫ u dx ≈ sum) — should be conserved exactly under
/// periodic BC with the linear stencil. Used for the convergence test.
pub fn integral(u: &[f32]) -> f32 {
    u.iter().sum()
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn periodic_bc_conserves_integral() {
        let w = 32;
        let h = 32;
        let mut u = vec![0.0f32; (w * h) as usize];
        // Localised hot spot.
        u[(h / 2 * w + w / 2) as usize] = 100.0;
        let mut u_next = vec![0.0f32; (w * h) as usize];
        let p = HeatParams::default();
        let s0 = integral(&u);
        for _ in 0..1000 {
            heat_step(&u, &mut u_next, w, h, &p);
            std::mem::swap(&mut u, &mut u_next);
        }
        let s1 = integral(&u);
        let rel = (s1 - s0).abs() / s0.abs();
        assert!(rel < 1e-3, "heat: integral drift {rel}");
    }

    #[test]
    fn diffusion_spreads_hot_spot() {
        let w = 32;
        let h = 32;
        let mut u = vec![0.0f32; (w * h) as usize];
        u[(h / 2 * w + w / 2) as usize] = 100.0;
        let mut u_next = vec![0.0f32; (w * h) as usize];
        let p = HeatParams::default();
        for _ in 0..50 {
            heat_step(&u, &mut u_next, w, h, &p);
            std::mem::swap(&mut u, &mut u_next);
        }
        let centre = u[(h / 2 * w + w / 2) as usize];
        let neighbour = u[(h / 2 * w + w / 2 + 1) as usize];
        assert!(centre < 100.0 && neighbour > 0.0, "heat should diffuse");
    }
}

This is the exact Rust file compiled into the WASM module the page just loaded. Every line that runs on your device is here. The receipt is a function of this code, not a bespoke benchmark.