diff --git a/.github/workflows/wasm-pages.yml b/.github/workflows/wasm-pages.yml new file mode 100644 index 0000000..34b6f2e --- /dev/null +++ b/.github/workflows/wasm-pages.yml @@ -0,0 +1,96 @@ +name: 🌐 Build WASM & Deploy Demo + +# Builds the WebAssembly module (issue #171) and publishes the browser demo to +# GitHub Pages. The .wasm is built here rather than committed: it is a megabyte +# of generated output, and a stale binary next to live sources is worse than a +# two minute build. + +# Publishes from dev. The github-pages environment has to allow dev as a +# deployment branch (Settings β†’ Environments β†’ github-pages). +on: + push: + branches: [ "dev" ] + workflow_dispatch: + +permissions: + contents: read + pages: write + id-token: write + +# One deployment at a time; queued runs supersede rather than race. +concurrency: + group: pages + cancel-in-progress: false + +jobs: + build: + runs-on: ubuntu-latest + steps: + - name: Checkout repository + uses: actions/checkout@v6 + # No LFS: the WASM build reads no test data, and emcmake turns the + # integrity check off by default. + + - name: Cache Emscripten SDK + id: emsdk-cache + uses: actions/cache@v4 + with: + path: ~/emsdk + key: emsdk-${{ runner.os }}-6.0.6 + + - name: Install Emscripten + if: steps.emsdk-cache.outputs.cache-hit != 'true' + run: | + # Pinned to the version the demo was developed and tested against. + git clone --depth 1 https://github.com/emscripten-core/emsdk.git ~/emsdk + ~/emsdk/emsdk install 6.0.6 + ~/emsdk/emsdk activate 6.0.6 + + - name: Build the WebAssembly module + run: | + source ~/emsdk/emsdk_env.sh + emcmake cmake -S . -B build-wasm -DCMAKE_BUILD_TYPE=Release + cmake --build build-wasm -j"$(nproc)" + + - name: Report size + run: | + { + echo "### WebAssembly build" + echo "| file | raw | gzip | brotli |" + echo "| --- | --- | --- | --- |" + for f in bin/forefire.wasm bin/forefire.mjs; do + printf '| `%s` | %s | %s | %s |\n' "$(basename "$f")" \ + "$(stat -c%s "$f" | numfmt --to=iec)" \ + "$(gzip -c "$f" | wc -c | numfmt --to=iec)" \ + "$(brotli -c "$f" | wc -c | numfmt --to=iec)" + done + } | tee -a "$GITHUB_STEP_SUMMARY" + + - name: Smoke test under Node + run: node bindings/wasm/smoke.mjs bin/forefire.mjs + + - name: Assemble the site + run: | + mkdir -p _site + cp bin/forefire.mjs bin/forefire.wasm _site/ + cp bindings/wasm/index.html bindings/wasm/worker.mjs \ + bindings/wasm/landscape.mjs bindings/wasm/terrain3d.mjs _site/ + # Pages serves through Jekyll unless told otherwise, which would + # swallow files it does not recognise. + touch _site/.nojekyll + + - name: Upload artifact + uses: actions/upload-pages-artifact@v3 + with: + path: _site + + deploy: + needs: build + runs-on: ubuntu-latest + environment: + name: github-pages + url: ${{ steps.deployment.outputs.page_url }} + steps: + - name: Deploy to GitHub Pages + id: deployment + uses: actions/deploy-pages@v4 diff --git a/CHANGELOG.md b/CHANGELOG.md index f6769e3..7334c65 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -60,6 +60,12 @@ Merged since v2.5.0, not yet released. typed arrays in place of NumPy arrays. `FOREFIRE_BUILD_CLI` and `FOREFIRE_BUILD_WASM` control the two executables. See `bindings/wasm/README.md`. ([#171]) +- A browser demo of the WebAssembly build (`bindings/wasm/index.html`): pick + a square of ground on a map, and the page builds fuel and elevation layers + from satellite and terrain tiles, runs the fire in a worker, and plays it + back in 2D or on 3D terrain. Explore mode lights a fire where you click; + Research mode exposes the model, its tuning factors, the numerics, and the + fuel moisture and classes. ([#171]) ### Fixed diff --git a/bindings/wasm/README.md b/bindings/wasm/README.md index 57eb87c..f5b0308 100644 --- a/bindings/wasm/README.md +++ b/bindings/wasm/README.md @@ -8,9 +8,10 @@ This is an experimental build, not a supported one yet. It answers the three questions the issue left open: does it link, does it run, how big is it. `.github/workflows/wasm.yml` builds the module and runs `smoke.mjs` on every -pull request. The `.wasm` is built in CI rather than committed: it is a -megabyte of generated output, and a stale binary sitting next to live sources -is worse than a two minute build. +pull request, and `.github/workflows/wasm-pages.yml` publishes the browser +demo to GitHub Pages on every push to `dev`. The `.wasm` is built in CI rather +than committed: it is a megabyte of generated output, and a stale binary +sitting next to live sources is worse than a two minute build. ## What made it possible @@ -55,6 +56,143 @@ Headless, under Node: node bindings/wasm/smoke.mjs bin/forefire.mjs ``` +In a browser β€” copy the four page files next to the build outputs and serve +the directory, since a `file://` page can neither fetch the `.wasm` nor start a +module worker: + +```sh +cp bindings/wasm/{index.html,worker.mjs,landscape.mjs,terrain3d.mjs} bin/ +python3 -m http.server -d bin +``` + +The solver runs in a worker (`worker.mjs`), which computes the **whole run up +front** β€” a few hundred milliseconds for a modest fire, a few seconds for one +that burns 800 ha β€” and the page then scrubs the result on a timeline. Off the +main thread is what keeps the page answering while that happens. + +Each run gets a fresh worker. The engine's state is process-wide (#175), so a +second `FireDomain` in the same module keeps the first one's propagation model, +and switching models between runs would silently do nothing. + +### Explore and Research + +The page has two modes over the same map: + +- **Explore**, the default: click the map to light a fire, pick a wind on the + dial or from a few scenarios, and watch it play back. +- **Research** adds the propagation model, `propagationSpeedAdjustmentFactor` + and `windReductionFactor`, the time step, step count, `perimeterResolution` + and `spatialIncrement`, dead and live fuel moisture and bed depths (the `Md`, + `Ml` and `e` columns of the fuel table), and the excess-green cut-offs + between fuel classes. Each control shows the ForeFire name it sets. + +The models on offer are the ones that run on this fuel table with only slope +and wind: Rothermel, BalbiNov2011Curv, BalbiUnsteady, TroisPourcent and Iso. +Balbi2015, Balbi2020 and BalbiNov2011TMdMl also need moisture and temperature +layers and abort without them; BalbiNov2011 places the front far from the +ignition. + +### The landscape comes from map tiles + +`landscape.mjs` fetches two tiled layers for a square of ground and decodes +both into layer arrays. No NetCDF, no GIS stack: tiles in, an `Int32Array` and +a `Float64Array` out, straight into `addIndexLayer` and `addScalarLayer`. + +The page opens on a 4 km square of the Ajaccio hinterland β€” the site of +`tests/runff/run.ff` β€” but the map is live: drag to pan, scroll to resize the +domain between 500 m and 30 km, and click to light a fire. The tile zoom +is chosen from the domain size so the mosaic stays around a thousand pixels and +36 tiles whatever the scale. Panning and zooming keep the run playing: each +run remembers the square it was computed on, and its fronts are drawn wherever +that square now sits on screen. The fetch is debounced, and the canvas previews +the move with the tiles it already has. + +The ignition point is stored as a longitude and latitude rather than as domain +metres, so panning and zooming leave it on the same patch of ground instead of +sliding it across the map. + +**Fuel**, from satellite imagery, classified by excess green `g - (r+b)/2` +into three rows of the repository's own `tests/runff/fuels.csv`: + +| Class | Cut | `fuels.csv` row | Share of the site | +| --- | --- | --- | --- | +| shrub / maquis | β‰₯ 20 | index 5, bed depth 0.6 m | ~44% | +| sparse / grass | β‰₯ 10 | index 4, bed depth 0.19 m | ~32% | +| bare, built, water | below | index 0, no fuel bed | ~24% | + +**Elevation**, from Terrarium-encoded SRTM, `r*256 + g + b/256 - 32768` metres. +The site runs from 11 m to 761 m. + +Both are load-bearing, because the default model is **Rothermel** with the +parameters from `tests/runff/params.ff`, and its rate of spread reads the fuel bed and the +slope. Index 0 has no fuel bed, so the fire stops at clearings, tracks and +rooftops. And on a bare slope with no wind at all, ten minutes of spread grows +with the grade while the downhill side stays put: + +| slope | uphill | downhill | +| --- | --- | --- | +| flat | 47 m | 45 m | +| 30% | 99 m | 45 m | +| 60% | 210 m | 43 m | + +> **Layer order matters.** The `altitude` layer must be registered *before* the +> propagation layer. `DataBroker` derives the slope layer at the moment a layer +> named `altitude` is registered; a model that asks for `slope` first gets a +> constant zero-altitude stand-in, and the terrain silently does nothing. + +If the tile hosts are unreachable the page falls back to uniform fuel on flat +ground and says so. + +### Knowing when the fire stops + +A fire that stops does not leave a stationary perimeter behind: ForeFire +discards the front and `print[]` returns nothing at all, which on a canvas is +indistinguishable from a bug. The worker detects the front list emptying, drops +that frame, and reports why it ended β€” `left-domain` if the last live front was +against the boundary, `burnt-out` otherwise. A slower death, where the +perimeter survives but stops advancing, is caught by a burnt-area plateau over +a twelve-step window (`stalled`), and an ignition on bare ground is reported as +`never-caught` rather than as a fire that stopped. + +### 3D terrain + +`terrain3d.mjs` draws the elevation field as a textured mesh in raw WebGL β€” no +library, since a heightfield with one directional light is not worth a megabyte +of dependency. The texture is whatever the 2D view painted, so the two views +cannot disagree. Drag to orbit, scroll to zoom the ground in or out (which +fetches new terrain, as in 2D), and the vertical exaggeration is adjustable +because 750 m of relief over 4 km is subtle at true scale. The view degrades to +2D if WebGL is unavailable. + +Switching between the two views is a **camera move, not a swap**. `blendCamera` +interpolates between the orbit and a flat pose β€” straight down through a +near-orthographic lens (a 0.12 rad field of view from far enough back), framed +so the mesh's unit square exactly fills a square viewport, north up, with the +vertical exaggeration wound to zero. That pose draws what the plain 2D canvas +draws, to within 0.03% of a half-width, so the GL canvas can take over while +still flat and then tilt up: no cross-fade, no second renderer, and no visible +handover. Going back holds the GL canvas until the camera is flat again. + +The flat pose stops a third of a degree short of vertical on purpose β€” straight +down puts the view direction on the up vector and `lookAt` divides by zero. + +Clicking the terrain lights a fire there, which needs a ray cast rather +than a coordinate transform: `rayFor` rebuilds the world-space ray through a +pixel from the same camera basis `lookAt` uses, and `raycastHeightfield` +marches it until it drops through the surface, then bisects onto the crossing. +It stops at the first ground it meets, so a click cannot land on a slope hidden +behind a ridge. Picking samples the mesh's own vertex heights, not the +full-resolution field, or a click could land somewhere the visible surface is +not. + +`mvpFor`, `blendCamera`, `rayFor`, `raycastHeightfield` and `projectPoint` are +exported so the camera can be checked without a GL context β€” feed the mesh's corners through and confirm they land in front of +the camera and inside the clip volume. That test exists because the first +version of `multiply` indexed its matrices row-major while `perspective`, +`lookAt` and `uniformMatrix4fv` were all column-major. It still returned a +matrix, and the mesh still drew; it just drew as streaks radiating from a +vanishing point. + ## The API `forefire.mjs` is an ES module exporting a `createForeFire()` factory. It @@ -96,6 +234,8 @@ index with, the same shuffle `_pyforefire.cpp` performs. of the worked example in the NumPy documentation. Run against a native `pyforefire` built from the same tree, the two agree exactly: 349 front nodes, x 1980..3905, y 1763..2308, the front stretched downwind as the doc describes. +- **Speed.** The turning-wind demo runs 15 steps of 20 s over a 4 km domain in + about 1.2 s under Node. ## Known limits diff --git a/bindings/wasm/forefire_wasm.cpp b/bindings/wasm/forefire_wasm.cpp index dc403a1..58d9a29 100644 --- a/bindings/wasm/forefire_wasm.cpp +++ b/bindings/wasm/forefire_wasm.cpp @@ -105,11 +105,11 @@ class ForeFireWasm { /*! \brief Attach a model-backed layer, e.g. ("propagation", "WindDriven", * "propagationModel"). Distinct from the addLayer *command*, which * builds constant layers instead. */ - void addLayer(const std::string &type, const std::string &name, + bool addLayer(const std::string &type, const std::string &name, const std::string &key) { FireDomain *domain = executor.getDomain(); - if (domain == 0) return; - domain->addLayer(type, name, key); + if (domain == 0) return false; + return domain->addLayer(type, name, key); } /*! \brief Register a continuous field (altitude, wind, moisture, ...). */ diff --git a/bindings/wasm/index.html b/bindings/wasm/index.html new file mode 100644 index 0000000..d2bdbda --- /dev/null +++ b/bindings/wasm/index.html @@ -0,0 +1,1107 @@ + + + + + + + ForeFire + + + +
+
+ + +
ForeFirewildfire spread, running in this tab
+ +
+ +
click to light a fire Β· drag to move Β· scroll to zoom
+
+
+ + +
β€”elapsed
+
β€”burnt
+
+
+
+ + + + + + diff --git a/bindings/wasm/landscape.mjs b/bindings/wasm/landscape.mjs new file mode 100644 index 0000000..6915c84 --- /dev/null +++ b/bindings/wasm/landscape.mjs @@ -0,0 +1,284 @@ +/** + * Builds a ForeFire landscape out of public map tiles, in the browser. + * + * There is no NetCDF here and no GIS stack. Satellite tiles go onto a canvas + * and their pixels are classified into fuel indices; terrain tiles go onto a + * second canvas and decode to metres. The two arrays that come out are exactly + * what addIndexLayer and addScalarLayer take. + * + * Both are load-bearing. Rothermel's rate of spread reads the fuel bed and the + * slope, so the classified imagery stops the fire at clearings and roads, and + * the terrain makes it run uphill: on a bare slope with no wind at all, spread + * grows from 47 m to 210 m over ten minutes between flat ground and a 60% + * grade, while the downhill side stays put. + */ + +const TILE = 256; +const IMAGERY = + "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile"; +// Terrarium-encoded SRTM/NED, elevation = r*256 + g + b/256 - 32768 metres. +const TERRAIN = "https://s3.amazonaws.com/elevation-tiles-prod/terrarium"; + +/*! The three fuel classes, as rows of tests/runff/fuels.csv β€” the repository's + * own Rothermel table, not invented numbers. Index 0 has a fuel bed depth of + * zero, which is what makes it a barrier. */ +const FUEL_ROWS = { + 5: "5;626.0;600.0;0.1;1.0;4325.0;5844.0;0.6;1.393;0.201;8.3;1.0;300;70000;18802000.0;18802000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3", + 4: "4;613.0;538.0;0.1;1.0;4357.0;6524.0;0.19;1.286;0.085;8.3;1.0;300;70000;18677000.0;18677000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3", + 0: "0;563.0;522.0;0.1;1.0;6099.0;7273.0;0;0.764;0.352;8.3;1.0;300;70000;18169000.0;18167000.0;1800;1000;600;0.3;2.5e-05;4.0;0.3", +}; +const FUEL_HEADER = + "Index;Rhod;Rhol;Md;Ml;sd;sl;e;Sigmad;Sigmal;stoch;RhoA;Ta;Tau0;Deltah;DeltaH;Cp;Cpa;Ti;X0;r00;Blai;me"; + +export const FUELS = [ + { index: 5, label: "shrub / maquis", note: "depth 0.6 m", rgb: [46, 125, 50] }, + { index: 4, label: "sparse / grass", note: "depth 0.19 m", rgb: [190, 180, 70] }, + { index: 0, label: "bare, built, water", note: "no fuel bed", rgb: [110, 100, 95] }, +]; + +export const FUELS_TABLE = [FUEL_HEADER, ...FUELS.map((f) => FUEL_ROWS[f.index])].join("\n"); + +/*! The defaults below are the values in tests/runff/fuels.csv. */ +export const FUEL_DEFAULTS = { deadMoisture: 0.1, liveMoisture: 1.0, depth: { 5: 0.6, 4: 0.19 } }; + +/*! + * The fuel table with dead and live moisture (Md, Ml) and the burnable classes' + * bed depth (e) overridden. Index 0 keeps its zero depth: it is the barrier. + */ +export function fuelsTable({ deadMoisture, liveMoisture, depth }) { + const columns = FUEL_HEADER.split(";"); + const col = (name) => columns.indexOf(name); + const rows = FUELS.map((f) => { + const cells = FUEL_ROWS[f.index].split(";"); + cells[col("Md")] = String(deadMoisture); + cells[col("Ml")] = String(liveMoisture); + if (depth[f.index] !== undefined) cells[col("e")] = String(depth[f.index]); + return cells.join(";"); + }); + return [FUEL_HEADER, ...rows].join("\n"); +} + +/*! Excess-green cut-offs between the classes; see classify(). */ +export const CLASS_DEFAULTS = { shrub: 20, grass: 10 }; + +/*! + * Fuel indices from the per-cell excess green `g - (r+b)/2`. Excess green + * separates vegetation from rock, road and roof about as well as anything this + * cheap can; water lands below both cuts. + */ +export function classify(greenness, { shrub, grass }) { + const fuel = new Int32Array(greenness.length); + const counts = { 5: 0, 4: 0, 0: 0 }; + for (let i = 0; i < greenness.length; i++) { + const g = greenness[i]; + const index = g >= shrub ? 5 : g >= grass ? 4 : 0; + fuel[i] = index; + counts[index]++; + } + return { fuel, counts }; +} + +/*! Web Mercator: fractional pixel coordinates at a given zoom. */ +function project(lon, lat, zoom) { + const n = 2 ** zoom; + const rad = (lat * Math.PI) / 180; + return [ + ((lon + 180) / 360) * n * TILE, + ((1 - Math.log(Math.tan(rad) + 1 / Math.cos(rad)) / Math.PI) / 2) * n * TILE, + ]; +} + +/*! Ground resolution in metres per pixel, which is what ties the tile grid to + * the simulation's Cartesian metres. */ +function resolution(lat, zoom) { + return (156543.03392804097 * Math.cos((lat * Math.PI) / 180)) / 2 ** zoom; +} + +/*! Metres per degree, near enough for panning a few kilometres. */ +export const M_PER_DEG_LAT = 111320; +export const mPerDegLon = (lat) => 111320 * Math.cos((lat * Math.PI) / 180); + +export const SIDE_LIMITS = { min: 500, max: 30000 }; + +/*! + * Pick a tile zoom for a domain of `side` metres. + * + * Chosen so the mosaic stays around a thousand pixels whatever the domain + * size: the tile count per fetch stays roughly constant, and the imagery keeps + * pace with the ground it covers. Capped at 15, past which the terrain tiles + * are upsampled rather than more detailed. + */ +export function zoomFor(lat, side, targetPx = 1024) { + const ideal = Math.log2((156543.03392804097 * Math.cos((lat * Math.PI) / 180) * targetPx) / side); + return Math.max(10, Math.min(15, Math.round(ideal))); +} + +function loadImage(url) { + return new Promise((resolve, reject) => { + const img = new Image(); + // Required to read the pixels back out. Both hosts send + // Access-Control-Allow-Origin, so the canvas stays untainted. + img.crossOrigin = "anonymous"; + img.onload = () => resolve(img); + img.onerror = () => reject(new Error(`tile ${url} failed`)); + img.src = url; + }); +} + +/*! Fetch a tiled layer into one canvas covering the square. `url` differs per + * host: Esri numbers its path z/row/col, the terrain tiles z/x/y. */ +async function mosaic(url, { zoom, left, top, sidePx }) { + const canvas = document.createElement("canvas"); + canvas.width = canvas.height = sidePx; + const ctx = canvas.getContext("2d"); + const jobs = []; + for (let ty = Math.floor(top / TILE); ty <= Math.floor((top + sidePx) / TILE); ty++) + for (let tx = Math.floor(left / TILE); tx <= Math.floor((left + sidePx) / TILE); tx++) + jobs.push( + loadImage(url(zoom, tx, ty)).then((img) => + ctx.drawImage(img, tx * TILE - left, ty * TILE - top) + ) + ); + await Promise.all(jobs); + return canvas; +} + +/*! Resample a canvas to the simulation grid and hand back its pixels. */ +function resample(canvas, grid) { + const small = document.createElement("canvas"); + small.width = small.height = grid; + const ctx = small.getContext("2d"); + // Nearest neighbour for the DEM: averaging terrarium's packed bytes across + // pixels is meaningless, since the red channel carries 256 m per unit. + ctx.imageSmoothingEnabled = false; + ctx.drawImage(canvas, 0, 0, grid, grid); + return ctx.getImageData(0, 0, grid, grid).data; +} + +/*! + * Fetch the imagery and terrain covering a `side`-metre square centred on + * (lon, lat), and derive a fuel map and an elevation field from them. + * + * Returns the display canvas plus the two layer arrays, both in the layers' + * (t,z,y,x) order with row 0 at the domain's south edge. Throws if the tiles + * cannot be fetched, so the caller can fall back. + */ +export async function buildLandscape({ lon, lat, side, grid, zoom = zoomFor(lat, side) }) { + const metresPerPx = resolution(lat, zoom); + const sidePx = Math.round(side / metresPerPx); + const [cx, cy] = project(lon, lat, zoom); + // Whole pixels: a tile drawn at a fractional offset is antialiased along + // its edges, which shows as a seam in the imagery and, in the terrain, + // decodes the half-blended edge pixels to spikes hundreds of metres high. + const box = { zoom, left: Math.round(cx - sidePx / 2), top: Math.round(cy - sidePx / 2), sidePx }; + + const [image, terrain] = await Promise.all([ + mosaic((z, x, y) => `${IMAGERY}/${z}/${y}/${x}`, box), + mosaic((z, x, y) => `${TERRAIN}/${z}/${x}/${y}.png`, box), + ]); + + const rgb = resample(image, grid); + const dem = resample(terrain, grid); + + const greenness = new Float32Array(grid * grid); + const altitude = new Float64Array(grid * grid); + let lo = Infinity; + let hi = -Infinity; + + for (let row = 0; row < grid; row++) { + // Image rows run north to south; the layers' row 0 is the domain's + // south edge, so the vertical order flips here. + const y = grid - 1 - row; + for (let x = 0; x < grid; x++) { + const i = (row * grid + x) * 4; + + greenness[x + grid * y] = rgb[i + 1] - (rgb[i] + rgb[i + 2]) / 2; + + const metres = dem[i] * 256 + dem[i + 1] + dem[i + 2] / 256 - 32768; + altitude[x + grid * y] = metres; + if (metres < lo) lo = metres; + if (metres > hi) hi = metres; + } + } + + // Kept alongside the fuel map so the class cut-offs can move without + // fetching the tiles again. + const { fuel, counts } = classify(greenness, CLASS_DEFAULTS); + + return { + image, + greenness, + fuel, + altitude, + nx: grid, + ny: grid, + counts, + metresPerPx, + sidePx, + zoom, + // The square this landscape actually covers, which the page needs to + // keep straight from the square the user is currently looking at. + site: { lon, lat, side }, + elevation: { lo, hi }, + }; +} + +/*! A translucent raster of the fuel classes, drawn at simulation resolution and + * scaled up by the caller. */ +export function fuelOverlay(fuel, grid) { + const canvas = document.createElement("canvas"); + canvas.width = canvas.height = grid; + const ctx = canvas.getContext("2d"); + const img = ctx.createImageData(grid, grid); + const palette = Object.fromEntries(FUELS.map((f) => [f.index, f.rgb])); + for (let row = 0; row < grid; row++) { + const y = grid - 1 - row; // back to image order + for (let x = 0; x < grid; x++) { + const [r, g, b] = palette[fuel[x + grid * y]]; + const i = (row * grid + x) * 4; + img.data[i] = r; + img.data[i + 1] = g; + img.data[i + 2] = b; + img.data[i + 3] = 255; + } + } + ctx.putImageData(img, 0, 0); + return canvas; +} + +/*! Relief shading from the elevation field, as a multiply layer over the + * imagery. Satellite pictures taken near local noon carry almost no shadow, + * so without this the terrain is invisible in plan view. */ +export function hillshade(altitude, grid, cellMetres, { azimuth = 315, altitudeAngle = 40 } = {}) { + const canvas = document.createElement("canvas"); + canvas.width = canvas.height = grid; + const ctx = canvas.getContext("2d"); + const img = ctx.createImageData(grid, grid); + + const az = ((360 - azimuth + 90) * Math.PI) / 180; + const alt = (altitudeAngle * Math.PI) / 180; + const at = (x, y) => + altitude[Math.min(grid - 1, Math.max(0, x)) + grid * Math.min(grid - 1, Math.max(0, y))]; + + for (let row = 0; row < grid; row++) { + const y = grid - 1 - row; + for (let x = 0; x < grid; x++) { + const dzdx = (at(x + 1, y) - at(x - 1, y)) / (2 * cellMetres); + const dzdy = (at(x, y + 1) - at(x, y - 1)) / (2 * cellMetres); + const slope = Math.atan(Math.hypot(dzdx, dzdy)); + const aspect = Math.atan2(dzdy, -dzdx); + let shade = + Math.cos(alt) * Math.cos(slope) + + Math.sin(alt) * Math.sin(slope) * Math.cos(az - aspect); + shade = Math.max(0, Math.min(1, shade)); + // Rendered as black at varying alpha, to be composited with + // "multiply" β€” mid-grey means no change, so alpha tracks how far + // the shade falls below neutral. + const i = (row * grid + x) * 4; + img.data[i] = img.data[i + 1] = img.data[i + 2] = 0; + img.data[i + 3] = Math.round(255 * Math.max(0, 0.85 - shade) * 0.9); + } + } + ctx.putImageData(img, 0, 0); + return canvas; +} diff --git a/bindings/wasm/smoke.mjs b/bindings/wasm/smoke.mjs index 0f44956..490c61a 100644 --- a/bindings/wasm/smoke.mjs +++ b/bindings/wasm/smoke.mjs @@ -61,10 +61,18 @@ windU.fill(1.0, 0, W * H); const windV = new Float64Array(2 * W * H); windV.fill(1.0, W * H); -ff.addLayer("propagation", "WindDriven", "propagationModel"); -console.log("fuel layer:", ff.addIndexLayer("table", "fuel", 0, 0, 0, L, L, 0, H, W, 1, 1, fuel)); -console.log("windU layer:", ff.addScalarLayer("windScalDir", "windU", 0, 0, 0, L, L, 0, H, W, 2, 1, windU)); -console.log("windV layer:", ff.addScalarLayer("windScalDir", "windV", 0, 0, 0, L, L, 0, H, W, 2, 1, windV)); +const layers = { + propagation: ff.addLayer("propagation", "WindDriven", "propagationModel"), + fuel: ff.addIndexLayer("table", "fuel", 0, 0, 0, L, L, 0, H, W, 1, 1, fuel), + windU: ff.addScalarLayer("windScalDir", "windU", 0, 0, 0, L, L, 0, H, W, 2, 1, windU), + windV: ff.addScalarLayer("windScalDir", "windV", 0, 0, 0, L, L, 0, H, W, 2, 1, windV), +}; +console.log("layers:", layers); +const missing = Object.entries(layers).filter(([, ok]) => !ok); +if (missing.length) { + console.error("layers not registered: " + missing.map(([k]) => k).join(", ")); + process.exit(1); +} ff.execute("startFire[loc=(2000,2000,0.);t=0]"); diff --git a/bindings/wasm/terrain3d.mjs b/bindings/wasm/terrain3d.mjs new file mode 100644 index 0000000..f9d8fae --- /dev/null +++ b/bindings/wasm/terrain3d.mjs @@ -0,0 +1,465 @@ +/** + * A small WebGL terrain view for the browser demo (issue #171). + * + * The elevation field the simulation runs on is drawn as a mesh, textured with + * whatever the 2D view painted β€” imagery, fuel, burnt area, fronts. One + * drawing path, two presentations: the map is rendered once to an offscreen + * canvas and uploaded here as a texture, so the two views can never disagree. + * + * Deliberately no library: the whole thing is a textured heightfield with one + * directional light, which is not worth a megabyte of dependency. + */ + +const VERT = ` +attribute vec3 aPos; +attribute vec3 aNormal; +attribute vec2 aUV; +uniform mat4 uMVP; +varying vec2 vUV; +varying float vShade; +void main() { + vUV = aUV; + vec3 sun = normalize(vec3(-0.6, 0.5, 0.62)); + // Lifted off zero so shadowed faces stay legible rather than going black. + vShade = 0.55 + 0.45 * max(dot(normalize(aNormal), sun), 0.0); + gl_Position = uMVP * vec4(aPos, 1.0); +}`; + +const FRAG = ` +precision mediump float; +uniform sampler2D uMap; +varying vec2 vUV; +varying float vShade; +void main() { + vec4 c = texture2D(uMap, vUV); + gl_FragColor = vec4(c.rgb * vShade, 1.0); +}`; + +function compile(gl, type, src) { + const s = gl.createShader(type); + gl.shaderSource(s, src); + gl.compileShader(s); + if (!gl.getShaderParameter(s, gl.COMPILE_STATUS)) + throw new Error("shader: " + gl.getShaderInfoLog(s)); + return s; +} + +/*! Column-major, like everything else here and like uniformMatrix4fv with + * transpose=false: element (row r, column c) lives at index c*4 + r. Indexing + * this the other way round still type-checks and still produces a matrix; it + * just projects the mesh into streaks radiating from a vanishing point. */ +function multiply(a, b) { + const o = new Float32Array(16); + for (let c = 0; c < 4; c++) + for (let r = 0; r < 4; r++) { + let sum = 0; + for (let k = 0; k < 4; k++) sum += a[k * 4 + r] * b[c * 4 + k]; + o[c * 4 + r] = sum; + } + return o; +} + +function perspective(fovy, aspect, near, far) { + const f = 1 / Math.tan(fovy / 2); + return new Float32Array([ + f / aspect, 0, 0, 0, + 0, f, 0, 0, + 0, 0, (far + near) / (near - far), -1, + 0, 0, (2 * far * near) / (near - far), 0, + ]); +} + +function lookAt(eye, target, up) { + const sub = (a, b) => [a[0] - b[0], a[1] - b[1], a[2] - b[2]]; + const norm = (v) => { + const l = Math.hypot(...v) || 1; + return [v[0] / l, v[1] / l, v[2] / l]; + }; + const cross = (a, b) => [ + a[1] * b[2] - a[2] * b[1], + a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0], + ]; + const dot = (a, b) => a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; + const z = norm(sub(eye, target)); + const x = norm(cross(up, z)); + const y = cross(z, x); + return new Float32Array([ + x[0], y[0], z[0], 0, + x[1], y[1], z[1], 0, + x[2], y[2], z[2], 0, + -dot(x, eye), -dot(y, eye), -dot(z, eye), 1, + ]); +} + +const ORBIT_FOV = Math.PI / 4; +// Not Ο€/2: straight down puts the view direction on the up vector and lookAt +// divides by zero. A third of a degree short of vertical reads as flat and +// stays conditioned. +const FLAT_ELEVATION = Math.PI / 2 - 0.005; +// A long lens from far away is orthographic in all but name, which is what +// makes the flat pose match the plain 2D map rather than merely resemble it. +const FLAT_FOV = 0.12; + +const lerp = (a, b, t) => a + (b - a) * t; + +/*! + * Blend between the flat map pose and the orbit pose. + * + * `t` of 0 looks straight down through a near-orthographic lens, framed so the + * mesh's unit square exactly fills a square viewport β€” pixel for pixel what the + * 2D canvas draws, which is what lets the two views swap unseen. `t` of 1 is + * the orbit the user has dragged to. Everything in between is a camera move, + * so the transition needs no cross-fade and no second renderer. + */ +export function blendCamera({ azimuth, elevation, distance, exaggeration }, t) { + const fov = lerp(FLAT_FOV, ORBIT_FOV, t); + // Half the square subtends half the vertical field of view: the exact fit. + const flatDistance = 0.5 / Math.tan(fov / 2); + return { + // North-up at the flat end, or the map would land rotated. + azimuth: lerp(0, azimuth, t), + elevation: lerp(FLAT_ELEVATION, elevation, t), + distance: lerp(flatDistance, distance, t), + exaggeration: exaggeration * t, + fov, + }; +} + +/*! + * The full model-view-projection for a camera orbiting the origin. + * + * Exported so it can be checked without a GL context: feed it corners of the + * mesh and confirm they land inside the clip volume. That is the test the + * column-major slip above would have failed. + */ +export function mvpFor({ azimuth, elevation, distance, aspect, exaggeration, side, fov = ORBIT_FOV }) { + const eye = [ + distance * Math.cos(elevation) * Math.sin(azimuth), + -distance * Math.cos(elevation) * Math.cos(azimuth), + distance * Math.sin(elevation), + ]; + const view = multiply( + perspective(fov, aspect, 0.02, 40), + lookAt(eye, [0, 0, 0], [0, 0, 1]) + ); + // x and y are a unit square standing for `side` metres, while z is still + // metres; dividing by `side` puts them back in proportion, and the + // exaggeration then lifts the relief to something readable. + const zScale = exaggeration / side; + const model = new Float32Array([1, 0, 0, 0, 0, 1, 0, 0, 0, 0, zScale, 0, 0, 0, 0, 1]); + return multiply(view, model); +} + +/*! + * The world-space ray through a point on the screen. + * + * Rebuilt from the camera basis rather than by inverting the projection: the + * basis is the same one `lookAt` builds, so the two cannot drift apart, and + * there is no 4x4 inverse to get subtly wrong. + */ +export function rayFor({ azimuth, elevation, distance, fov = ORBIT_FOV }, ndcX, ndcY, aspect) { + const eye = [ + distance * Math.cos(elevation) * Math.sin(azimuth), + -distance * Math.cos(elevation) * Math.cos(azimuth), + distance * Math.sin(elevation), + ]; + const norm = (v) => { + const l = Math.hypot(...v) || 1; + return [v[0] / l, v[1] / l, v[2] / l]; + }; + const cross = (a, b) => [ + a[1] * b[2] - a[2] * b[1], + a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0], + ]; + const forward = norm([-eye[0], -eye[1], -eye[2]]); + const right = norm(cross(forward, [0, 0, 1])); + const up = cross(right, forward); + const half = Math.tan(fov / 2); + const dir = norm([ + forward[0] + right[0] * ndcX * half * aspect + up[0] * ndcY * half, + forward[1] + right[1] * ndcX * half * aspect + up[1] * ndcY * half, + forward[2] + right[2] * ndcX * half * aspect + up[2] * ndcY * half, + ]); + return { origin: eye, dir }; +} + +/*! + * March a ray until it drops through a heightfield, and return where. + * + * `heightAt(u, v)` gives the surface in the same units as the ray, over the + * unit square the mesh occupies. Returns normalised `[u, v]` in [0, 1], or + * null if the ray never meets the ground β€” clicking the sky. + * + * Fixed-step march to bracket the crossing, then bisection to land on it. A + * closed-form solve is not available against an interpolated grid, and a coarse + * march alone would let the point slide down steep faces. + */ +export function raycastHeightfield({ origin, dir }, heightAt, { reach = 6, steps = 2048 } = {}) { + const at = (t) => [origin[0] + dir[0] * t, origin[1] + dir[1] * t, origin[2] + dir[2] * t]; + const gap = (t) => { + const p = at(t); + const u = p[0] + 0.5; + const v = p[1] + 0.5; + if (u < 0 || u > 1 || v < 0 || v > 1) return null; // off the square + return p[2] - heightAt(u, v); + }; + + const step = reach / steps; + let prevT = null; + let prevGap = null; + for (let i = 0; i <= steps; i++) { + const t = i * step; + const g = gap(t); + if (g === null) { + prevT = null; + prevGap = null; + continue; + } + if (prevGap !== null && prevGap > 0 && g <= 0) { + let lo = prevT; + let hi = t; + for (let k = 0; k < 24; k++) { + const mid = (lo + hi) / 2; + const gm = gap(mid); + if (gm === null) break; + if (gm > 0) lo = mid; + else hi = mid; + } + const p = at((lo + hi) / 2); + const clamp01 = (x) => Math.max(0, Math.min(1, x)); + return [clamp01(p[0] + 0.5), clamp01(p[1] + 0.5)]; + } + prevT = t; + prevGap = g; + } + return null; +} + +/*! Apply a column-major 4x4 to a point, returning normalised device + * coordinates and the clip-space w (positive means in front of the camera). */ +export function projectPoint(m, [x, y, z]) { + const c = [0, 1, 2, 3].map((r) => m[0 + r] * x + m[4 + r] * y + m[8 + r] * z + m[12 + r]); + return { ndc: [c[0] / c[3], c[1] / c[3], c[2] / c[3]], w: c[3] }; +} + +/*! + * Build the terrain view. + * + * `altitude` is the simulation's elevation field in (t,z,y,x) order β€” row 0 at + * the domain's south edge β€” `grid` its side, `side` the domain in metres. + * Returns { setMap, render, orbit, zoom, setExaggeration, dispose }, or throws + * if WebGL is unavailable so the caller can stay in 2D. + */ +export function createTerrainView(canvas, { altitude, grid, side, resolution = 129 }) { + const gl = canvas.getContext("webgl", { antialias: true, alpha: false }); + if (!gl) throw new Error("WebGL unavailable"); + + const program = gl.createProgram(); + gl.attachShader(program, compile(gl, gl.VERTEX_SHADER, VERT)); + gl.attachShader(program, compile(gl, gl.FRAGMENT_SHADER, FRAG)); + gl.linkProgram(program); + if (!gl.getProgramParameter(program, gl.LINK_STATUS)) + throw new Error("link: " + gl.getProgramInfoLog(program)); + + // Sample the elevation field down to the mesh resolution. + const N = Math.min(resolution, grid); + const at = (ix, iy) => { + const x = Math.min(grid - 1, Math.round((ix / (N - 1)) * (grid - 1))); + const y = Math.min(grid - 1, Math.round((iy / (N - 1)) * (grid - 1))); + return altitude[x + grid * y]; + }; + + let lo = Infinity; + let hi = -Infinity; + for (const v of altitude) { + if (v < lo) lo = v; + if (v > hi) hi = v; + } + const mid = (lo + hi) / 2; + + // Positions are normalised to a unit square in x/y so the camera framing is + // independent of the domain size; heights stay in metres, scaled at draw + // time by the exaggeration. + const positions = new Float32Array(N * N * 3); + const uvs = new Float32Array(N * N * 2); + // The same heights the mesh is drawn from, kept for picking: sampling the + // full-resolution field instead would let a click land somewhere the + // visible surface is not. + const meshZ = new Float32Array(N * N); + for (let iy = 0; iy < N; iy++) + for (let ix = 0; ix < N; ix++) { + const i = iy * N + ix; + meshZ[i] = at(ix, iy) - mid; + positions[i * 3] = ix / (N - 1) - 0.5; + positions[i * 3 + 1] = iy / (N - 1) - 0.5; + positions[i * 3 + 2] = meshZ[i]; + uvs[i * 2] = ix / (N - 1); + // The texture is the map as drawn, whose rows run north to south; + // the mesh's y runs south to north, so v flips. + uvs[i * 2 + 1] = 1 - iy / (N - 1); + } + + const indices = []; + for (let iy = 0; iy < N - 1; iy++) + for (let ix = 0; ix < N - 1; ix++) { + const a = iy * N + ix; + indices.push(a, a + 1, a + N, a + 1, a + N + 1, a + N); + } + const indexArray = N * N > 65535 ? new Uint32Array(indices) : new Uint16Array(indices); + const uintExt = N * N > 65535 ? gl.getExtension("OES_element_index_uint") : true; + if (!uintExt) throw new Error("mesh too large for this WebGL context"); + + // Normals are computed once on the CPU: the terrain never changes, and this + // keeps the vertex shader to a single matrix multiply. + const cell = side / (N - 1); + const normals = new Float32Array(N * N * 3); + for (let iy = 0; iy < N; iy++) + for (let ix = 0; ix < N; ix++) { + const i = iy * N + ix; + const dzdx = (at(Math.min(N - 1, ix + 1), iy) - at(Math.max(0, ix - 1), iy)) / (2 * cell); + const dzdy = (at(ix, Math.min(N - 1, iy + 1)) - at(ix, Math.max(0, iy - 1))) / (2 * cell); + const l = Math.hypot(-dzdx, -dzdy, 1); + normals[i * 3] = -dzdx / l; + normals[i * 3 + 1] = -dzdy / l; + normals[i * 3 + 2] = 1 / l; + } + + const buffer = (data, target = gl.ARRAY_BUFFER) => { + const b = gl.createBuffer(); + gl.bindBuffer(target, b); + gl.bufferData(target, data, gl.STATIC_DRAW); + return b; + }; + const posBuf = buffer(positions); + const normBuf = buffer(normals); + const uvBuf = buffer(uvs); + const idxBuf = buffer(indexArray, gl.ELEMENT_ARRAY_BUFFER); + + const texture = gl.createTexture(); + gl.bindTexture(gl.TEXTURE_2D, texture); + // No mipmaps and clamped edges: the map canvas is not power-of-two. + gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MIN_FILTER, gl.LINEAR); + gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_MAG_FILTER, gl.LINEAR); + gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_S, gl.CLAMP_TO_EDGE); + gl.texParameteri(gl.TEXTURE_2D, gl.TEXTURE_WRAP_T, gl.CLAMP_TO_EDGE); + + const loc = { + pos: gl.getAttribLocation(program, "aPos"), + normal: gl.getAttribLocation(program, "aNormal"), + uv: gl.getAttribLocation(program, "aUV"), + mvp: gl.getUniformLocation(program, "uMVP"), + map: gl.getUniformLocation(program, "uMap"), + }; + + const camera = { azimuth: -0.6, elevation: 0.72, distance: 1.9 }; + let exaggeration = 2.5; + let transition = 1; // 0 = flat map, 1 = orbit + let mapReady = false; + + gl.enable(gl.DEPTH_TEST); + gl.clearColor(0.05, 0.04, 0.03, 1); + + function setMap(source) { + gl.bindTexture(gl.TEXTURE_2D, texture); + gl.pixelStorei(gl.UNPACK_FLIP_Y_WEBGL, false); + gl.texImage2D(gl.TEXTURE_2D, 0, gl.RGBA, gl.RGBA, gl.UNSIGNED_BYTE, source); + mapReady = true; + } + + function render() { + if (!mapReady) return; + const w = canvas.width; + const h = canvas.height; + gl.viewport(0, 0, w, h); + gl.clear(gl.COLOR_BUFFER_BIT | gl.DEPTH_BUFFER_BIT); + + const posed = blendCamera({ ...camera, exaggeration }, transition); + const final = mvpFor({ ...posed, aspect: w / h, side }); + + gl.useProgram(program); + gl.uniformMatrix4fv(loc.mvp, false, final); + gl.activeTexture(gl.TEXTURE0); + gl.bindTexture(gl.TEXTURE_2D, texture); + gl.uniform1i(loc.map, 0); + + const bind = (buf, index, size) => { + gl.bindBuffer(gl.ARRAY_BUFFER, buf); + gl.enableVertexAttribArray(index); + gl.vertexAttribPointer(index, size, gl.FLOAT, false, 0, 0); + }; + bind(posBuf, loc.pos, 3); + bind(normBuf, loc.normal, 3); + bind(uvBuf, loc.uv, 2); + gl.bindBuffer(gl.ELEMENT_ARRAY_BUFFER, idxBuf); + gl.drawElements( + gl.TRIANGLES, + indexArray.length, + indexArray.BYTES_PER_ELEMENT === 4 ? gl.UNSIGNED_INT : gl.UNSIGNED_SHORT, + 0 + ); + } + + return { + setMap, + render, + /*! The orbit pose, so a page rebuilding the mesh for new ground can + * carry the view over instead of resetting it. */ + camera, + orbit(dAz, dEl) { + camera.azimuth += dAz; + // Stop short of straight down, where the up vector degenerates. + camera.elevation = Math.max(0.08, Math.min(1.45, camera.elevation + dEl)); + }, + zoom(factor) { + // The near limit stops short of the point where the mesh's corners + // all leave the frame and the view loses its footing. + camera.distance = Math.max(0.9, Math.min(5, camera.distance * factor)); + }, + setExaggeration(v) { + exaggeration = v; + }, + /*! 0 = the flat map, 1 = the orbit. Driven by the page's animation. */ + setTransition(t) { + transition = Math.max(0, Math.min(1, t)); + }, + /*! + * Where a click on the canvas meets the ground. + * + * `ndcX`/`ndcY` are normalised device coordinates, -1 to 1 with y up. + * Returns the point in the domain's metres, or null for a click that + * missed the terrain entirely. + */ + pick(ndcX, ndcY, aspect) { + const cam = blendCamera({ ...camera, exaggeration }, transition); + const zScale = cam.exaggeration / side; + // Bilinear over the mesh vertices, which is what the triangles + // interpolate between; nearest would quantise a click to a cell. + const heightAt = (u, v) => { + const gx = u * (N - 1); + const gy = v * (N - 1); + const ix = Math.max(0, Math.min(N - 2, Math.floor(gx))); + const iy = Math.max(0, Math.min(N - 2, Math.floor(gy))); + const fx = gx - ix; + const fy = gy - iy; + const z00 = meshZ[iy * N + ix]; + const z10 = meshZ[iy * N + ix + 1]; + const z01 = meshZ[(iy + 1) * N + ix]; + const z11 = meshZ[(iy + 1) * N + ix + 1]; + return ( + ((z00 * (1 - fx) + z10 * fx) * (1 - fy) + + (z01 * (1 - fx) + z11 * fx) * fy) * + zScale + ); + }; + const hit = raycastHeightfield( + rayFor(cam, ndcX, ndcY, aspect), + heightAt, + { reach: cam.distance + 3 } + ); + return hit ? [hit[0] * side, hit[1] * side] : null; + }, + elevationRange: { lo, hi }, + }; +} diff --git a/bindings/wasm/worker.mjs b/bindings/wasm/worker.mjs new file mode 100644 index 0000000..6a59891 --- /dev/null +++ b/bindings/wasm/worker.mjs @@ -0,0 +1,238 @@ +/** + * Simulation worker for the browser demo (issue #171). + * + * A whole run costs a few hundred milliseconds for a modest fire and a few + * seconds for one that burns eight hundred hectares, so the worker computes + * the entire timeline up front and posts it in one go; the page then scrubs + * through frames rather than waiting on the solver. Running it here rather + * than on the main thread is what keeps the page answering during that. + * + * The fuel map and terrain arrive from the main thread, which is where they + * have to be built: decoding map tiles needs a canvas, and workers have no DOM. + * + * Each worker runs exactly one simulation. The engine's state is process-wide + * (#175): a second FireDomain in the same module keeps the first one's + * propagation model, so changing the model between runs would silently do + * nothing. A fresh worker is a fresh module, and costs a few tens of ms. + */ + +import createForeFire from "./forefire.mjs"; + +let Module = null; + +// The parameters from tests/runff/params.ff. The page overrides the ones its +// research panel exposes. +const SOLVER = { + spatialIncrement: 3, + minimalPropagativeFrontDepth: 20, + perimeterResolution: 10, + initialFrontDepth: 0.1, + relax: 0.5, + smoothing: 0, + minSpeed: 0.009, + bmapLayer: 1, + windReductionFactor: 0.4, + propagationSpeedAdjustmentFactor: 0.6, +}; + +function build(land, model, solver) { + const { side, nx, ny } = land; + const ff = new Module.ForeFire(); + + ff.setString("fuelsTable", land.fuelsTable); + ff.setDouble("defaultFuelType", land.defaultFuel); + for (const [key, value] of Object.entries(solver)) ff.setDouble(key, value); + ff.setString("propagationModel", model); + + ff.execute(`FireDomain[sw=(0,0,0);ne=(${side},${side},0);t=0]`); + + const fuel = land.fuel ?? new Int32Array(nx * ny).fill(land.defaultFuel); + // Two-plane wind layers: plane 0 carries the u response, plane 1 the v. + const windU = new Float64Array(2 * nx * ny); + windU.fill(1, 0, nx * ny); + const windV = new Float64Array(2 * nx * ny); + windV.fill(1, nx * ny); + + const S = side; + // Altitude first, and before the propagation layer: registering a layer + // named "altitude" is what makes the DataBroker derive the slope layer, and + // a propagation model that asks for "slope" before one exists gets a + // constant zero-altitude stand-in that never updates. + if (land.altitude) + ff.addScalarLayer("data", "altitude", 0, 0, 0, S, S, 0, nx, ny, 1, 1, land.altitude); + + const ok = [ + ff.addLayer("propagation", model, "propagationModel"), + ff.addIndexLayer("table", "fuel", 0, 0, 0, S, S, 0, nx, ny, 1, 1, fuel), + ff.addScalarLayer("windScalDir", "windU", 0, 0, 0, S, S, 0, nx, ny, 2, 1, windU), + ff.addScalarLayer("windScalDir", "windV", 0, 0, 0, S, S, 0, nx, ny, 2, 1, windV), + ]; + if (ok.some((v) => !v)) throw new Error("a layer failed to register"); + return ff; +} + +/*! Split print[] into one coordinate array per front. + * + * Fronts nest β€” an inner front is printed indented inside its parent β€” so + * every FireFront line starts a fresh polygon. Running them together is what + * draws spurious chords straight across the burn. + */ +function parseFronts(text) { + const fronts = []; + let current = null; + for (const line of text.split("\n")) { + if (line.includes("FireFront[")) { + current = []; + fronts.push(current); + } else if (current && line.includes("FireNode[")) { + const m = /loc=\(([-\d.eE+]+),([-\d.eE+]+),/.exec(line); + if (m) current.push(parseFloat(m[1]), parseFloat(m[2])); + } + } + // A front of one or two nodes is a degenerate sliver, not a perimeter. + return fronts.filter((f) => f.length >= 6).map((f) => Float32Array.from(f)); +} + +/*! Shoelace area, summed over fronts. Inner fronts wind the other way, so + * signed areas cancel and holes are subtracted for free. */ +function burntArea(fronts) { + let total = 0; + for (const f of fronts) { + let a = 0; + for (let i = 0, n = f.length; i < n; i += 2) { + const j = (i + 2) % n; + a += f[i] * f[j + 1] - f[j] * f[i + 1]; + } + total += a / 2; + } + return Math.abs(total); +} + +/*! Did this front run into the edge of the domain? */ +function touchesEdge(fronts, side, margin) { + for (const f of fronts) + for (let i = 0; i < f.length; i += 2) + if ( + f[i] < margin || + f[i] > side - margin || + f[i + 1] < margin || + f[i + 1] > side - margin + ) + return true; + return false; +} + +/*! + * Run the whole timeline and return it, along with how it ended. + * + * A fire that stops does not linger as a stationary perimeter: ForeFire + * discards the front and print[] comes back with nothing at all, which on a + * canvas looks identical to a bug. So the end of the run is detected on the + * front list emptying, the empty frame is dropped, and the cause is reported β€” + * a front that was touching the domain boundary left the map, anything else + * burnt out. A slower death, where the perimeter survives but stops advancing, + * is caught by the area-plateau test below. + */ +function simulate({ land, model, params, steps, dt, wind, ignition }) { + const solver = { ...SOLVER, ...params }; + const ff = build(land, model, solver); + const [ix, iy] = ignition; + ff.execute(`startFire[loc=(${ix},${iy},0.);t=0]`); + + const frames = []; + const started = performance.now(); + let outcome = { kind: "running" }; + // The plateau test is judged over a window, not step to step: a young + // Rothermel fire adds only a few hundred square metres per step, which is + // indistinguishable from noise one step at a time. + const WINDOW = 12; + + const push = (step, ms) => { + const fronts = parseFronts(ff.execute("print[]")); + frames.push({ + step, + t: ff.getTime(), + ms, + fronts, + nodes: fronts.reduce((n, f) => n + f.length / 2, 0), + area: burntArea(fronts), + windDeg: windAt(step, wind).deg, + }); + }; + push(0, 0); + + for (let step = 1; step <= steps; step++) { + const t0 = performance.now(); + const { vx, vy } = windAt(step, wind); + ff.execute(`trigger[wind;loc=(0.,0.,0.);vel=(${vx.toFixed(4)},${vy.toFixed(4)},0)]`); + ff.execute(`step[dt=${dt}]`); + push(step, performance.now() - t0); + + if (frames.at(-1).fronts.length === 0) { + const last = frames.at(-2); + frames.pop(); // nothing to draw, and it would blank the canvas + outcome = { + kind: + last && touchesEdge(last.fronts, land.side, 3 * solver.perimeterResolution) + ? "left-domain" + : "burnt-out", + t: last ? last.t : 0, + }; + break; + } + + if (frames.length > WINDOW) { + const now = frames.at(-1).area; + const then = frames.at(-1 - WINDOW).area; + // Stalled when the burnt area has stopped growing in any meaningful + // sense β€” under half a percent over the window, and under one fuel + // cell's worth of ground. + const cell = (land.side / land.nx) ** 2; + if (now - then < Math.max(cell, 0.005 * then)) { + outcome = { kind: "stalled", t: frames.at(-1 - WINDOW).t }; + break; + } + } + if (step % 25 === 0) + self.postMessage({ type: "progress", done: step, total: steps }); + } + + if (outcome.kind === "running") outcome = { kind: "step-limit", t: frames.at(-1).t }; + // An ignition on bare ground never becomes a fire at all, which is worth + // saying plainly rather than reporting as a fire that stopped. + const peak = Math.max(...frames.map((f) => f.area)); + if (outcome.kind === "stalled" && peak < 10 * (land.side / land.nx) ** 2) + outcome = { kind: "never-caught", t: outcome.t }; + + return { + frames, + outcome, + wallMs: performance.now() - started, + domain: land.side, + peakArea: peak, + }; +} + +function windAt(step, wind) { + // Wrapped, or a turning run reports "864Β°" after an hour of simulated time. + const deg = (wind.turning ? wind.angleDeg + step * wind.turnPerStep : wind.angleDeg) % 360; + const rad = (deg * Math.PI) / 180; + return { vx: wind.speed * Math.cos(rad), vy: wind.speed * Math.sin(rad), deg }; +} + +self.onmessage = ({ data: msg }) => { + if (msg.type === "simulate") { + try { + const result = simulate(msg.run); + self.postMessage( + { type: "timeline", ...result }, + result.frames.flatMap((f) => f.fronts.map((c) => c.buffer)) + ); + } catch (err) { + self.postMessage({ type: "failed", message: String(err && err.message) }); + } + } +}; + +Module = await createForeFire(); +self.postMessage({ type: "ready", version: Module.version() });