diff --git a/avaframe/com4FlowPy/benchmark/.gitignore b/avaframe/com4FlowPy/benchmark/.gitignore new file mode 100644 index 000000000..02b1fddda --- /dev/null +++ b/avaframe/com4FlowPy/benchmark/.gitignore @@ -0,0 +1,7 @@ +# Local run artifacts and machine-specific fixtures — never committed. +# Concrete fixtures carry absolute data paths; only the *.template is tracked. +runs/ +fixtures/*.ini + +# debug repro (not for commit) +repro_tile.py diff --git a/avaframe/com4FlowPy/benchmark/README.md b/avaframe/com4FlowPy/benchmark/README.md new file mode 100644 index 000000000..bdd62d8a1 --- /dev/null +++ b/avaframe/com4FlowPy/benchmark/README.md @@ -0,0 +1,46 @@ +# com4FlowPy Numba-engine benchmark & validation harness + +`flowpy_bench.py` runs com4FlowPy in-process on a *fixture* (a `.ini` describing +inputs + `[GENERAL]` config), measures **wall-clock** and **peak RAM**, and +optionally compares every output raster against a **reference directory**. + +It exercises the normal `com4FlowPyMain` code path (as in the `useCustomPaths` +branch of `runCom4FlowPy.py`); the only added knob is `--engine`, which selects +the compute engine: + +* `--engine python` — the stock `Cell`-based BFS (default; always available). +* `--engine numba` — the ported `@njit` kernel (once merged). + +## Usage + +```bash +# Freeze the stock-engine oracle on the correctness fixture +python flowpy_bench.py fixtures/connaught_frequent.ini --engine python --tag oracle + +# A/B the numba engine against it (bit-equivalence check) +python flowpy_bench.py fixtures/connaught_frequent.ini --engine numba + +# Benchmark on BFW's own 10 m run and cross-check vs their stored outputs +python flowpy_bench.py fixtures/bfw_10m.ini --engine python +python flowpy_bench.py fixtures/bfw_10m.ini --engine numba --cpu 8 +``` + +Each run writes `runs/__/summary.json` with timing, +peak RAM, and per-output comparison metrics. + +## Fixtures + +| Fixture | Role | Config highlights | +|---|---|---| +| `connaught_frequent.ini` / `connaught_extreme.ini` | correctness oracle (fast A/B) | autoATES: `forestDetrainment`, uniform α (28 / 18) | +| `bfw_10m.ini` | speed benchmark + reference cross-check | BFW: `forestFriction`, `variableAlpha`, `variableUmax`, α-layer + Umax-layer | + +Fixture `.ini`s carry absolute paths to input rasters. The rasters themselves +are **not** committed (size/licensing); edit the paths for your machine, or copy +a fixture to `*.local.ini` (gitignored). `runs/` is also gitignored. + +## Comparison metrics (per output raster, over pixels valid in both) + +`n_valid`, `n_exact`, `n_diff` (|Δ| > `--tol`), `frac_diff_ppm`, +`max_abs_diff`, `mean_abs_diff`. Use `--tol 0` for bit-exact; a small tol to +ignore IEEE-754 ULP noise from operation-order differences. diff --git a/avaframe/com4FlowPy/benchmark/REPORT.md b/avaframe/com4FlowPy/benchmark/REPORT.md new file mode 100644 index 000000000..2949b2bf5 --- /dev/null +++ b/avaframe/com4FlowPy/benchmark/REPORT.md @@ -0,0 +1,232 @@ +# com4FlowPy Numba compute-engine — validation & benchmark report + +*Living document — updated at each phase. Audience: BFW / AvaFrame maintainers.* + +## Goal + +Add an optional Numba (`@njit`) compute engine to com4FlowPy that reproduces +the existing Python engine's results while cutting per-tile runtime, so large / +repeated runs (e.g. multi-site ATES workflows) become tractable. The engine is +selected by a config flag (`engine = python | numba`, default `python`); tiling, +multiprocessing, I/O, forest, and variable-parameter handling are unchanged. + +Motivation: a prototype on standalone FlowPy (fork of `avaframe/FlowPy`, +`foreste_detrainment`) reached ~100× by rewriting the per-cell BFS as a single +`@njit` function over flat arrays. com4FlowPy already parallelizes (tiling + +`Pool` over release chunks) but still runs the inner BFS through the Python +`Cell` class — the target of this port. + +## Method + +`benchmark/flowpy_bench.py` runs com4FlowPy in-process via `com4FlowPyMain` +(the `useCustomPaths` path), recording wall-clock, peak RAM, and per-output +raster comparison against a reference directory. Fixtures under +`benchmark/fixtures/` pin exact inputs + `[GENERAL]` config. Correctness is +established by A/B comparison of the two engines on identical inputs (the +Python engine is the oracle); a run is additionally cross-checked against BFW's +own stored outputs. + +Machine: 16 physical cores / 32 threads. + +All results below were re-verified on the current rebased code +(OpenNHM/AvaFrame master @6e02e45d + this branch), with peak memory measured via +PSS. Current-code headline (8 cores): BFW 10 m 302 s → 6.7 s (**45×**); 5 m +Connaught extreme 5714 s → 193 s (**30×**). numba (float64) is bit-identical to +the Python engine on both (zDelta / fpTravelAngleMax; 5 m differs only at +flat-terrain routing knife-edges, see below). PSS memory is on par with the +Python engine (BFW: 8.6 vs 7.1 GiB; 5 m: 23.1 vs 19.9 GiB). + +## Phase 0 — baselines & oracle (stock Python engine) — DONE + +| Fixture | Domain | Config | Wall | Peak RAM | Cores | +|---|---|---|---|---|---| +| Connaught frequent | 222×222 @ 21.3 m | forestDetrainment, α=28 | 17.8 s | 6.0 GiB | 16 | +| Connaught extreme | 222×222 @ 21.3 m | forestDetrainment, α=18 | 109 s | 7.6 GiB | 16 | +| **BFW 10 m** | 482×513 @ 10 m | forestFriction, variableAlpha, variableUmax, α-layer + Umax-layer | **299 s** | 17.4 GiB | 8 | + +**Key validation:** on the BFW 10 m fixture — reproducing BFW's own run config +exactly — the stock Python engine reproduces BFW's stored `zDelta` and +`fpTravelAngleMax` **bit-for-bit** (0 differing pixels / 81,795 valid; max |Δ| = 0). +This confirms the harness, config translation, and variable-parameter wiring are +faithful, and gives a bit-exact target for the Numba engine. + +Benchmark target to beat: **299 s on 8 cores (BFW 10 m)**. + +## Phase 1 — Numba kernel — DONE (integration) + +New module `flowCoreNumba.py`: the per-release-pixel BFS as a single `@njit` +function over flat arrays, plus `calculationNumba(args)` — a drop-in for +`flowCore.calculation()` (same args, same 13-tuple). `run()` dispatches to it +when `engine=numba` (default `python`); tiling / multiprocessing / merge / I-O +unchanged. infra/back-calculation and previewMode fall back to the Python engine. + +Faithful to the `fluxDistOldVersion=False` path: float persistence (no int16 +truncation), g=9.81, forest friction/detrainment/frictionLayer with the +not-start / FSI>0 / skipForestDist guards, forestInteraction counting, variable +alpha/max_z/exponent per release pixel, default flux distribution +(count = dist≥threshold, sub-threshold redistribution, conservation correction, +deposition). + +## Phase 2 — numerical equivalence — IN PROGRESS + +A/B (numba vs Python engine), correctness fixtures — **bit-identical**: + +| Fixture | Outputs checked | Result | +|---|---|---| +| Connaught frequent | zDelta, routFluxSum, fpTravelAngleMax, cellCounts | 0 diff / 21,335 px, max\|Δ\|=0 | +| Connaught extreme | same | 0 diff / 25,528 px, max\|Δ\|=0 | + +BFW 10 m (forestFriction + variableAlpha + variableUmax), numba vs Python baseline +and vs BFW's own stored outputs: + +| Output | numba vs baseline | vs BFW stored | +|---|---|---| +| zDelta | 0 diff / 81,795 | **bit-identical** | +| fpTravelAngleMax | 0 diff / 81,795 | **bit-identical** | +| flux | 6 px (0.07 ppm), max\|Δ\|=2.2e-4 | n/a (not stored) | +| cellCounts | 20 px (0.24 ppm), max\|Δ\|=1 | n/a (not stored) | + +The scientifically primary outputs (zDelta, travel angle) are bit-exact. The +sub-ppm flux / cellCounts diffs are IEEE-754 operation-order artifacts (scalar +kernel vs numpy vectorized ops): a child flux landing within ~1 ULP of the flux +threshold is included by one engine and not the other, nudging cellCounts by 1 +and flux by a rounding-level amount at a handful of boundary pixels. This is the +same class/scale of artifact documented for the standalone-FlowPy prototype +(sub-110 ppm). It does not affect zDelta because those boundary cells take their +zDelta from other paths. + +## Phase 3 — speedup & determinism — headline result + +BFW 10 m, matched core counts: + +| Domain (8 cores) | Python | Numba | Speedup | +|---|---|---|---| +| BFW 10 m | 302 s | 6.7 s | **45×** | +| Connaught 5 m extreme | 5714 s (95 min) | 193 s | **30×** | + +On the 10 m domain numba runtime is **flat across cpu=4/8/16 (~6.5 s)** — small +enough that it is no longer compute-bound (fixed costs: cold JIT compile + tiling +I/O + merge dominate). At 5 m it is genuinely compute-bound (193 s). The realised +speedup (~30–45× here) therefore depends on domain size and how compute-bound the +run is. Because adding workers past ~8 gives little wall-clock gain (numba +saturates the tile-parallelism early) but more memory, **~8 workers is the +practical sweet spot for high-resolution runs**. + +**Determinism:** numba output is bit-identical across cpu=4/8/16 on all four +outputs (zDelta, fpTravelAngleMax, flux, cellCounts — 0 differing pixels). The +engine is fully deterministic w.r.t. worker count; the sub-ppm flux/cellCounts +differences vs the Python engine are therefore a stable scalar-vs-numpy artifact, +not a parallelism effect. + +Peak RAM (PSS) is on par with the Python engine and is dominated by com4FlowPy's +multiprocessing, not the numba workspace (the per-BFS queue was nonetheless +right-sized from `2·H·W` to a 128k grow-and-retry buffer — a correctness/hygiene +fix, not a peak-RAM lever). A warm JIT cache removes the one-time compile cost. + +## Engine comparison for BFW + +com4FlowPy exposes two engines via `engine =`: `python` (the Cell-based +reference) and `numba` (JIT, double precision). Default is `python`. + +BFW 10 m, 8 cores, vs the python reference: + +| Engine | Wall | Peak RAM | zDelta | fpTravelAngleMax | flux | cellCounts | +|---|---|---|---|---|---|---| +| python | 302 s | 7.1 GiB | ref | ref | ref | ref | +| numba | 6.7 s | 8.6 GiB | **0 diff** | **0 diff** | 6 px / 2e-4 | 20 px / 1 | + +(Wall/RAM on current master, 8 cores, PSS memory. diff = differing pixels / max\|Δ\| +over 81,795 valid px; the sub-ppm flux/cellCounts differences are the documented +scalar-vs-numpy artifact, and are bit-identical to BFW's own stored zDelta and +fpTravelAngleMax.) + +**A single-precision (float32) variant was evaluated and rejected.** It was +prototyped on the theory it would cut memory on high-resolution runs; measurement +showed otherwise — no speedup (numba32 193 s vs numba 192 s at 5 m), no memory +saving (24.9 vs 23.1 GiB PSS — if anything slightly higher, since peak RAM is +dominated by com4FlowPy's multiprocessing, not the kernel workspace), and a +precision cost that is *worst exactly at high resolution* (5 m: ~68 % of cells +differ from the Python engine, zDelta up to 62 m). It offered no benefit where it +was meant to help, so it is not part of this contribution. + +## 5 m (high-resolution, compute-bound) — Connaught extreme + +964×1025 ≈ 988k cells; the domain that ran ~1.5 h in autoATESv3.0. 8 cores, PSS memory: + +| Engine | Wall | Peak RAM | +|---|---|---| +| python | 5714 s (95 min) | 19.9 GiB | +| numba | 193 s | 23.1 GiB | + +→ **numba ~30× faster**, memory on par with the Python engine. + +* **Wall time is flat across cpu=4/8/16** — at 5 m the domain splits into only a + few tiles, so parallelism saturates early and the remainder is serial + (tiling / merge / JIT). More workers add memory, not speed — so physical cores + are the ceiling worth using (hyperthreads never help), and fewer is fine. +* Peak RAM is dominated by com4FlowPy's multiprocessing (per-worker tile copies), + not the numba kernel workspace (a modest grow-and-retry queue). On + memory-constrained machines, fewer workers and/or a smaller `tileSize` lower + the peak. Memory is a non-issue at 10–25 m (small tiles). + +**Memory (measured with PSS — proportional set size — which correctly accounts +for shared copy-on-write pages across the worker pool):** at BFW 10 m, python@8 = +8.9 GiB, numba@8 = 8.6 GiB — numba is on par with (slightly below) the Python +engine. At 5 m, numba@8 = 25 GiB. Earlier RSS-summed figures (e.g. 76 GiB at 5 m) +over-counted shared memory ~2–3× and were misleading. The per-BFS queue workspace +was separately reduced from `2·H·W` (~300× over the observed ~3.8k peak depth) to +a modest 128k with grow-and-retry on overflow — this is a correctness/robustness +fix (a path can never be silently truncated) rather than a peak-RAM lever, since +the workspace was not the dominant allocation. + +### 5 m divergence — flat-terrain routing sensitivity (a model property, not an engine bug) + +The numba engine is deterministic (identical output run-to-run and across worker +counts) and bit-for-bit identical to the Python engine on every domain tested at +10–25 m (Connaught 21 m, BFW 10 m). On the finer 5 m run the two engines differ +in a small set of cells: zDelta in 124 of 444k valid cells (0.03 %, worst case +52 m at isolated points); routFluxSum and cellCounts in more cells (~20 % and +~5 %) but by small per-cell amounts (routFluxSum ≤ 1.3, cellCounts ≤ 21). + +**Plain-English cause.** FlowPy routes flow downhill, splitting it at each cell +among the downhill neighbours according to how steep the drop is. On a real slope +there is a clear downhill and the path is well-defined. On a flat valley floor the +neighbours sit within ~1 m of each other — there is no real "downhill," so which +neighbour the flow follows comes down to the last rounding bit of the arithmetic. +And because the model must keep routing until the energy line reaches the alpha +angle (18° here), the flow keeps spreading across the flat with no genuine +directional preference. The two engines evaluate the same formulas by a slightly +different route (per-element scalar math vs numpy array math), so at these +near-ties they occasionally send flow to a different — but equally valid — +neighbour, after which the branches rejoin. + +**Neither engine is "more correct," and this is not a tiling artifact.** We +reproduced the exact divergence inside a single isolated tile +(`benchmark/repro_tile.py`), so it is not a tile-seam or merge effect — it is the +fine 5 m grid resolving flat valley floors into many near-tie cells. To show the +result is genuinely ill-conditioned *in the model itself* (independent of which +engine computes it): nudging the DEM by a physically meaningless **±1 µm** — far +below any real elevation accuracy — shifts the Python engine's *own* reached set +by ~175 cells (±1 mm → 237; ±1 cm → 189). The numba engine is equally +deterministic and equally valid; it simply lands on a different one of these +near-tied branches. Crucially, the numba↔Python difference (~50 cells along a +path) is *smaller* than the change the Python result shows under a sub-micron +nudge. So at a flat-terrain fork the exact reached set is not a numerically +determinate quantity for *any* correct implementation — both engines are valid +IEEE-754 solutions of the same equations. + +**Impact, stated plainly.** The energy-line (zDelta) footprint is essentially +identical (0.03 % of cells differ); the larger routFluxSum/cellCounts pixel counts +are small in magnitude and concentrated at the margins of flow paths in flat +runout zones. Any effect on downstream ATES classes is therefore confined to +flat-terrain margins and expected to be minor — but we flag it rather than claim +bit-identical output everywhere. This ambiguity is an inherent property of +FlowPy's alpha-angle routing in flat terrain: the port neither introduces nor +removes it, and the Python engine exhibits it too. A slope-aware stopping/ +friction criterion for low-gradient terrain would address the root cause and is +noted as separate future work. + +## Phase 4 — upstream packaging — PENDING + +numba as an optional dependency with graceful fallback; config docs; commit the +Connaught/BFW fixtures as a replicable regression test; PR to avaframe/AvaFrame. diff --git a/avaframe/com4FlowPy/benchmark/fixtures/example.ini.template b/avaframe/com4FlowPy/benchmark/fixtures/example.ini.template new file mode 100644 index 000000000..0b25cf35b --- /dev/null +++ b/avaframe/com4FlowPy/benchmark/fixtures/example.ini.template @@ -0,0 +1,24 @@ +# Example fixture for flowpy_bench.py — copy to .ini and edit paths. +# Concrete *.ini files are gitignored (they carry machine-specific data paths). +[inputs] +demPath = /path/to/DEM.tif +releasePath = /path/to/release_or_pra_binary.tif +forestPath = /path/to/forest.tif ; optional (required if forest=True) +varAlphaPath = /path/to/alpha.tif ; optional (required if variableAlpha=True) +varUmaxPath = /path/to/umax.tif ; optional (required if variableUmaxLim=True) +referenceDir = /path/to/reference_outputs ; optional; compare outputs against these tifs + +[GENERAL] +# any com4FlowPy [GENERAL] key overrides the module default +alpha = 25 +exp = 8 +flux_threshold = 3.0e-4 +max_z = 270 +forest = True +forestModule = forestFriction ; forestFriction | forestDetrainment | forestFrictionLayer +forestInteraction = True +variableAlpha = False +variableUmaxLim = False +tileSize = 15000 +tileOverlap = 5000 +outputFiles = zDelta|fpTravelAngleMax|cellCounts|flux diff --git a/avaframe/com4FlowPy/benchmark/flowpy_bench.py b/avaframe/com4FlowPy/benchmark/flowpy_bench.py new file mode 100644 index 000000000..3c34baafe --- /dev/null +++ b/avaframe/com4FlowPy/benchmark/flowpy_bench.py @@ -0,0 +1,348 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +flowpy_bench.py — reproducible runner/benchmark for com4FlowPy. + +Runs com4FlowPy in-process on a fixture, measures wall-clock and peak RAM, and +(optionally) compares every output raster against a reference directory. It is +used to: + + * freeze a correctness oracle (stock Python engine output), + * A/B the Numba engine against the Python engine (bit-equivalence), + * benchmark the speedup on a realistic domain, + * cross-check a run against BFW's own stored com4FlowPy outputs. + +The runner mirrors the ``useCustomPaths`` branch of avaframe/runCom4FlowPy.py +(getModuleConfig + com4FlowPyMain), so it exercises the exact code path a normal +com4FlowPy run uses — the only knob it adds is ``engine`` (see --engine). + +Usage +----- + python flowpy_bench.py FIXTURE.ini [--engine python|numba] + [--out DIR] [--tag NAME] [--cpu N] + +Fixture .ini format +------------------- + [inputs] + demPath = /abs/path/DEM.tif + releasePath = /abs/path/release.tif + forestPath = /abs/path/forest.tif ; optional (required if forest=True) + varAlphaPath = /abs/path/alpha.tif ; optional (required if variableAlpha=True) + varUmaxPath = /abs/path/umax.tif ; optional (required if variableUmaxLim=True) + referenceDir = /abs/path/ref_outputs ; optional; compare outputs against these tifs + + [GENERAL] + ; any com4FlowPy [GENERAL] key overrides the module default, e.g. + alpha = 26 + forest = True + forestModule = forestFriction + variableAlpha = True + ... + +Comparison metrics per output raster (valid = non-nodata in both): + n_valid, n_exact, n_diff (|Δ| > --tol), frac_diff (ppm), + max_abs_diff, mean_abs_diff. +""" +import argparse +import configparser +import json +import shutil +import sys +import threading +import time +from datetime import datetime +from pathlib import Path + +import numpy as np + +from avaframe.com4FlowPy import com4FlowPy +from avaframe.in3Utils import cfgUtils + +try: + from avaframe.runCom4FlowPy import checkOutputFilesFormat +except Exception: # pragma: no cover - older/newer AvaFrame may differ + def checkOutputFilesFormat(s): + return s + +try: + import psutil +except Exception: + psutil = None + +try: + import rasterio +except Exception: + rasterio = None + + +# --------------------------------------------------------------------------- +# peak-RAM sampler (whole process tree, since com4FlowPy forks a Pool) +# --------------------------------------------------------------------------- +class PeakRAM: + def __init__(self, interval=0.1): + self.interval = interval + self.peak_bytes = 0 + self._stop = threading.Event() + self._thread = None + + def _tree_rss(self): + # Prefer PSS (proportional set size): shared copy-on-write pages across the + # forked worker pool are divided among sharers, so summing over the process + # tree gives the true physical footprint instead of over-counting shared libs + # and pre-fork data. Falls back to RSS where PSS is unavailable (non-Linux). + if psutil is None: + return 0 + try: + proc = psutil.Process() + procs = [proc] + proc.children(recursive=True) + total = 0 + for p in procs: + try: + mi = p.memory_full_info() + total += getattr(mi, "pss", None) or mi.rss + except (psutil.NoSuchProcess, psutil.AccessDenied): + pass + return total + except Exception: + return 0 + + def _run(self): + while not self._stop.is_set(): + self.peak_bytes = max(self.peak_bytes, self._tree_rss()) + self._stop.wait(self.interval) + + def __enter__(self): + if psutil is not None: + self._thread = threading.Thread(target=self._run, daemon=True) + self._thread.start() + return self + + def __exit__(self, *exc): + self._stop.set() + if self._thread is not None: + self._thread.join(timeout=1.0) + + @property + def peak_gib(self): + return self.peak_bytes / (1024 ** 3) + + +# --------------------------------------------------------------------------- +# config assembly +# --------------------------------------------------------------------------- +def build_cfg(fixture: configparser.ConfigParser, work_dir: Path, engine: str, + cpu_override: int | None): + """Return (cfgSetup, cfgPath, uid) built from module defaults + fixture.""" + cfg = cfgUtils.getModuleConfig(com4FlowPy, onlyDefault=True, toPrint=False) + gen = cfg["GENERAL"] + + # overlay every [GENERAL] key from the fixture onto the module defaults + if fixture.has_section("GENERAL"): + for key, val in fixture.items("GENERAL"): + gen[key] = val + + # engine selector: harmless extra key for the stock code; read by the Numba port. + gen["engine"] = engine + + if cpu_override is not None: + gen["cpuCount"] = str(cpu_override) + elif not gen.get("cpuCount"): + n = psutil.cpu_count(logical=False) if psutil else 1 + gen["cpuCount"] = str(n or 1) + + # [inputs] read case-insensitively (decoupled from the case-preserving GENERAL overlay) + inp = {k.lower(): v for k, v in fixture.items("inputs")} + dem_path = Path(inp["dempath"]).expanduser() + release_path = Path(inp["releasepath"]).expanduser() + forest_path = Path(inp.get("forestpath", "") or "").expanduser() + + uid = cfgUtils.cfgHash(cfg) + res_dir = work_dir / f"res_{uid}" + if res_dir.exists(): + shutil.rmtree(res_dir, ignore_errors=True) + temp_dir = res_dir / "temp" + temp_dir.mkdir(parents=True, exist_ok=True) + + paths = cfg["PATHS"] + cfgPath = { + "workDir": work_dir, + "outDir": res_dir, + "resDir": res_dir, + "tempDir": temp_dir, + "demPath": dem_path, + "releasePath": release_path, + "relIdPath": Path(inp.get("relidpath", "") or ""), + "infraPath": Path(inp.get("infrapath", "") or ""), + "forestPath": forest_path, + "varUmaxPath": Path(inp.get("varumaxpath", "") or ""), + "varAlphaPath": Path(inp.get("varalphapath", "") or ""), + "varExponentPath": Path(inp.get("varexponentpath", "") or ""), + "deleteTemp": "False", + "outputFileFormat": paths.get("outputFileFormat", ".tif"), + "outputFiles": checkOutputFilesFormat(gen.get("outputFiles") + or paths.get("outputFiles")), + "outputNoDataValue": float(paths.get("outputNoDataValue", "-9999")), + "useCompression": paths.getboolean("useCompression", fallback=True), + "customDirs": "True", + "uid": uid, + "timeString": datetime.now().strftime("%Y%m%d_%H%M%S"), + } + return gen, cfgPath, uid + + +# --------------------------------------------------------------------------- +# raster comparison +# --------------------------------------------------------------------------- +def _read(path): + with rasterio.open(path) as ds: + arr = ds.read(1).astype(np.float64) + nod = ds.nodata + return arr, nod + + +def compare_raster(produced: Path, reference: Path, tol: float): + a, na = _read(produced) + b, nb = _read(reference) + if a.shape != b.shape: + return {"status": "shape_mismatch", + "produced_shape": list(a.shape), "reference_shape": list(b.shape)} + valid = np.ones(a.shape, dtype=bool) + for arr, nod in ((a, na), (b, nb)): + if nod is not None: + valid &= arr != nod + valid &= ~np.isnan(arr) + diff = np.abs(a - b) + over = valid & (diff > tol) + n_valid = int(valid.sum()) + return { + "status": "ok", + "n_valid": n_valid, + "n_exact": int((valid & (diff == 0)).sum()), + "n_diff": int(over.sum()), + "frac_diff_ppm": round(1e6 * over.sum() / n_valid, 3) if n_valid else None, + "max_abs_diff": float(diff[valid].max()) if n_valid else None, + "mean_abs_diff": float(diff[valid].mean()) if n_valid else None, + } + + +def find_reference(ref_dir: Path, output_name: str): + token = output_name.lower().replace("_", "") + cands = [p for p in ref_dir.glob("*.tif") + if p.stem.lower().replace("_", "").endswith(token)] + return sorted(cands)[-1] if cands else None + + +# --------------------------------------------------------------------------- +# main +# --------------------------------------------------------------------------- +def main(argv=None): + ap = argparse.ArgumentParser(description="com4FlowPy fixture runner / benchmark") + ap.add_argument("fixture", help="path to a fixture .ini") + ap.add_argument("--engine", default="python", choices=["python", "numba"], + help="compute engine: python (Cell-based) or numba " + "(JIT, float64, bit-exact). Default python.") + ap.add_argument("--out", default=None, help="output root dir (default: ./runs)") + ap.add_argument("--tag", default=None, help="label for this run") + ap.add_argument("--cpu", type=int, default=None, help="override cpuCount") + ap.add_argument("--tol", type=float, default=0.0, + help="abs tolerance for the diff count (default 0 = bit-exact)") + args = ap.parse_args(argv) + + if rasterio is None: + print("WARNING: rasterio not importable — comparison disabled", file=sys.stderr) + + fixture_path = Path(args.fixture).expanduser().resolve() + fixture = configparser.ConfigParser() + fixture.optionxform = str # preserve key case so camelCase GENERAL keys match AvaFrame's + fixture.read(fixture_path) + + tag = args.tag or fixture_path.stem + out_root = Path(args.out).expanduser() if args.out else (fixture_path.parent / "runs") + stamp = datetime.now().strftime("%Y%m%d_%H%M%S") + work_dir = out_root / f"{tag}_{args.engine}_{stamp}" + work_dir.mkdir(parents=True, exist_ok=True) + + gen, cfgPath, uid = build_cfg(fixture, work_dir, args.engine, args.cpu) + + print(f"[flowpy_bench] fixture={fixture_path.name} engine={args.engine} " + f"cpu={gen.get('cpuCount')} uid={uid}") + print(f"[flowpy_bench] alpha={gen.get('alpha')} forestModule={gen.get('forestModule')} " + f"variableAlpha={gen.get('variableAlpha')} variableUmaxLim={gen.get('variableUmaxLim')} " + f"max_z={gen.get('max_z')}") + print(f"[flowpy_bench] outputs={cfgPath['outputFiles']}") + print(f"[flowpy_bench] work_dir={work_dir}") + + t0 = time.perf_counter() + status = "success" + err = None + try: + with PeakRAM() as ram: + com4FlowPy.com4FlowPyMain(cfgPath, gen) + peak_gib = ram.peak_gib + except Exception as exc: # noqa: BLE001 + status = "failed" + err = repr(exc) + peak_gib = None + import traceback + traceback.print_exc() + wall_s = time.perf_counter() - t0 + + res_dir = cfgPath["resDir"] + produced = {} + for p in sorted(res_dir.glob("com4_*.tif")): + for name in cfgPath["outputFiles"].split("|"): + tok = name.lower().replace("_", "") + if p.stem.lower().replace("_", "").endswith(tok): + produced[name] = p + + comparisons = {} + inp_lower = {k.lower(): v for k, v in fixture.items("inputs")} if fixture.has_section("inputs") else {} + ref_dir_str = inp_lower.get("referencedir", "") + if ref_dir_str and rasterio is not None and status == "success": + ref_dir = Path(ref_dir_str).expanduser() + for name, ppath in produced.items(): + rpath = find_reference(ref_dir, name) + if rpath is not None: + comparisons[name] = {"reference": rpath.name, + **compare_raster(ppath, rpath, args.tol)} + else: + comparisons[name] = {"status": "no_reference"} + + summary = { + "fixture": str(fixture_path), + "tag": tag, + "engine": args.engine, + "uid": uid, + "status": status, + "error": err, + "wall_seconds": round(wall_s, 2), + "peak_ram_gib": round(peak_gib, 3) if peak_gib is not None else None, + "cpu_count": gen.get("cpuCount"), + "general": {k: gen.get(k) for k in ( + "alpha", "exp", "flux_threshold", "max_z", "forest", "forestModule", + "forestInteraction", "variableAlpha", "variableUmaxLim", + "maxAddedFrictionFor", "minAddedFrictionFor", "velThForFriction", + "maxDetrainmentFor", "minDetrainmentFor", "tileSize", "tileOverlap")}, + "outputs": {k: str(v) for k, v in produced.items()}, + "comparisons": comparisons, + "timestamp": stamp, + } + summary_path = work_dir / "summary.json" + summary_path.write_text(json.dumps(summary, indent=2)) + + print(f"\n[flowpy_bench] status={status} wall={wall_s:.2f}s " + f"peak_ram={summary['peak_ram_gib']} GiB") + if comparisons: + print("[flowpy_bench] comparison vs reference:") + for name, c in comparisons.items(): + if c.get("status") == "ok": + print(f" {name:20s} n_diff={c['n_diff']:>8} / {c['n_valid']:<8} " + f"({c['frac_diff_ppm']} ppm) max|Δ|={c['max_abs_diff']:.3g}") + else: + print(f" {name:20s} {c.get('status')}") + print(f"[flowpy_bench] summary → {summary_path}") + return 0 if status == "success" else 1 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/avaframe/com4FlowPy/com4FlowPy.py b/avaframe/com4FlowPy/com4FlowPy.py index 972f10da2..42b088331 100755 --- a/avaframe/com4FlowPy/com4FlowPy.py +++ b/avaframe/com4FlowPy/com4FlowPy.py @@ -78,6 +78,9 @@ def com4FlowPyMain(cfgPath, cfgSetup): # Flag for use of old flux distribution version modelParameters["fluxDistOldVersionBool"] = cfgSetup.getboolean("fluxDistOldVersion") + # compute engine: "python" (default, Cell-based) or "numba" (JIT kernel) + modelParameters["engine"] = cfgSetup.get("engine", "python").strip().lower() + # Tiling Parameters used for calculation of large model-domains tilingParameters = {} tilingParameters["tileSize"] = cfgSetup.getfloat("tileSize") # float(cfgSetup["tileSize"]) @@ -389,7 +392,9 @@ def checkInputLayerDimensions(modelParameters, modelPaths): def checkInputParameterValues(modelParameters, modelPaths): - """check if the input parameters alpha, uMaxLimit/ zDeltaMaxLimit, exponent + """check if the input parameters + are valid and within physically sensible limits + alpha, uMaxLimit/ zDeltaMaxLimit, exponent are within a physically sensible range Parameters @@ -399,6 +404,16 @@ def checkInputParameterValues(modelParameters, modelPaths): modelPaths: dict contains paths to input files """ + + engine = modelParameters["engine"] + validEngines = {"python","numba"} + + if engine not in validEngines: + raise ValueError( + f"Invalid engine '{engine}'. " + f"Valid engine options are {sorted(validEngines)}" + ) + alpha = modelParameters["alpha"] if alpha < 0 or alpha > 90: log.error("Error: Alpha value is not within a physically sensible range ([0,90]).") diff --git a/avaframe/com4FlowPy/com4FlowPyCfg.ini b/avaframe/com4FlowPy/com4FlowPyCfg.ini index 763164e3b..ca6da7c07 100644 --- a/avaframe/com4FlowPy/com4FlowPyCfg.ini +++ b/avaframe/com4FlowPy/com4FlowPyCfg.ini @@ -171,6 +171,14 @@ skipForestDist = 0 fluxDistOldVersion = False +# compute engine for the per-cell BFS. The 'numba' engine requires the optional +# 'numba' dependency (JIT-compiled kernel); it reproduces the 'python' engine's +# results and is much faster. infra/back-calculation, previewMode and relId +# outputs automatically fall back to the python engine. +# python : default, Cell-based reference implementation +# numba : JIT-compiled kernel (double precision), bit-for-bit vs 'python' +engine = python + #++++++++++++ Parameters for Tiling # tileSize: size of tiles in x and y direction in meters (if total size of) x # or y of input DEM is larger than tileSize, then the input raster diff --git a/avaframe/com4FlowPy/flowCore.py b/avaframe/com4FlowPy/flowCore.py index 434190c2c..99d0c9ae8 100644 --- a/avaframe/com4FlowPy/flowCore.py +++ b/avaframe/com4FlowPy/flowCore.py @@ -234,6 +234,28 @@ def run(optTuple): ) release_list = split_release(release, nChunks) + + # select compute engine: "numba" JIT kernel, else the default Python (Cell) path. + # numba does not (yet) implement infra/back-calculation, previewMode or the + # release-id (relIdPolygon/relIdCount) outputs, so fall back to the Python + # engine when any of those are active. + engine = optTuple[2].get("engine", "python") + calcFunc = calculation + if engine == "numba": + if infraBool or previewMode or relIdBool: + log.warning("engine=numba does not support infra/previewMode/relId outputs — " + "falling back to the Python engine for this run") + else: + try: + from avaframe.com4FlowPy.flowCoreNumba import calculationNumba + calcFunc = calculationNumba + log.info("Using numba compute engine") + except ImportError: + log.warning("engine=numba requested but 'numba' is not installed — " + "falling back to the Python engine") + else: + log.info("Using python compute engine") + log.info( "Multiprocessing starts, used Cores/Processes/Chunks: %i/%i/%i" % (MPOptions["nCPU"], nProcesses, nChunks) @@ -241,7 +263,7 @@ def run(optTuple): with Pool(processes=nProcesses) as pool: results = pool.map( - calculation, + calcFunc, [ [ # TODO: write in dicts: dem, diff --git a/avaframe/com4FlowPy/flowCoreNumba.py b/avaframe/com4FlowPy/flowCoreNumba.py new file mode 100644 index 000000000..6316827cf --- /dev/null +++ b/avaframe/com4FlowPy/flowCoreNumba.py @@ -0,0 +1,654 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Numba (@njit) compute engine for com4FlowPy. + +`calculationNumba(args)` is a drop-in replacement for `flowCore.calculation()`: +same argument list, same 14-element return tuple. It is selected per run via the +`engine = numba` config flag; `run()` dispatches to it instead of `calculation()` +while keeping tiling, multiprocessing, I/O and merging unchanged. + +The entire per-release-pixel BFS is expressed as a single @njit function over +flat numpy arrays (no Cell objects, no Python lists/dicts), which is where the +speedup comes from. The math is a faithful transcription of flowClass.Cell + +flowCore.calculation() for the fluxDistOldVersion=False path: + + * persistence is float (com4FlowPy tiles the DEM as float64 -> no int16 truncation), + * g = 9.81, forest friction/detrainment guards (not-start, FSI>0, skipForestDist), + * forestModule: forestFriction / forestDetrainment / forestFrictionLayer, + * forestInteraction counting, + * variable alpha / max_z / exponent resolved once per release pixel, + * default flux distribution: count = (dist >= threshold), sub-threshold mass + redistributed over >=threshold children, flux-conservation correction, and + deposition (fluxDep) when count == 0. + +Paths NOT handled here fall back to the Python engine (see run()): infra / +back-calculation and previewMode. +""" +import logging +import math + +import numpy as np +from numba import njit + +from avaframe.com4FlowPy.flowCore import get_start_idx + +log = logging.getLogger(__name__) + +_SQRT2 = math.sqrt(2.0) +_HALF_PI = math.pi / 2.0 +_DEG_PER_RAD = 180.0 / math.pi +_G = 9.81 + +# 9-neighbour layout, row-major: [0]=TL [1]=T [2]=TR / [3]=L [4]=C [5]=R / [6]=BL [7]=B [8]=BR +_DS = np.array([_SQRT2, 1.0, _SQRT2, 1.0, 0.0, 1.0, _SQRT2, 1.0, _SQRT2]) # z_alpha (center 0) +_DS_TANBETA = np.array([_SQRT2, 1.0, _SQRT2, 1.0, 1.0, 1.0, _SQRT2, 1.0, _SQRT2]) # distance (center 1) + +# forestModule codes +_FM_NONE = 0 +_FM_FRICTION = 1 +_FM_DETRAINMENT = 2 +_FM_FRICTIONLAYER = 3 + + +@njit(cache=True) +def _bfs_single(dem, forest, H, W, nodata, row_start, col_start, + cellsize, alpha, exp, flux_threshold, max_z_delta, + ds_cellsize, distance, + # forest scalars + forestBool, forestModuleCode, forestInteraction, + maxAddedFriction, minAddedFriction, noFrictionEffectZDelta, + maxDetrainment, minDetrainment, noDetrainmentEffectZDelta, + forestDetrainmentBool, frictionLayerRelative, skipForestDist, + fluxDistOldVersion, + # output arrays + zDeltaArray, fluxArray, countArray, zDeltaSumArray, zDeltaPathArray, + routFluxSumArray, depFluxSumArray, + fpMaxArray, fpMinArray, slArray, + travelMaxArray, travelMinArray, forestIntArray, + # workspace + pending_qidx, visited, + q_r, q_c, q_zdelta, q_flux, q_is_start, q_first_parent_start, + q_mindist, q_mindist3d, q_fic, q_isforest, + q_n_parents, q_pdir, q_pzd, q_pmd, q_pmd3d, + modified_r, modified_c): + """Run the BFS for one release pixel, accumulating into the output arrays.""" + if row_start < 1 or row_start >= H - 1 or col_start < 1 or col_start >= W - 1: + return 0 + for di in range(-1, 2): + for dj in range(-1, 2): + if dem[row_start + di, col_start + dj] == nodata: + return 0 + + tan_alpha_base = math.tan(math.radians(alpha)) + altitude_start = dem[row_start, col_start] + + q_head = 0 + q_tail = 0 + n_modified = 0 + + start_isforest = 1 if (forestInteraction and forest[row_start, col_start] > 0) else 0 + + q_r[0] = row_start + q_c[0] = col_start + q_zdelta[0] = 0.0 + q_flux[0] = 1.0 + q_is_start[0] = True + q_first_parent_start[0] = False + q_mindist[0] = 0.0 + q_mindist3d[0] = 0.0 + q_fic[0] = start_isforest + q_isforest[0] = start_isforest + q_n_parents[0] = 0 + pending_qidx[row_start, col_start] = 0 + modified_r[0] = row_start + modified_c[0] = col_start + n_modified = 1 + q_tail = 1 + + while q_head < q_tail: + r = q_r[q_head] + c = q_c[q_head] + q_pos = q_head + q_head += 1 + pending_qidx[r, c] = -1 + + # first-visit count (once per distinct cell reached in this BFS) + if visited[r, c] == 0: + countArray[r, c] += 1 + visited[r, c] = 1 + + z_delta = q_zdelta[q_pos] + flux = q_flux[q_pos] + is_start = q_is_start[q_pos] + first_parent_is_start = q_first_parent_start[q_pos] + n_parents = q_n_parents[q_pos] + fic = q_fic[q_pos] + isForest = q_isforest[q_pos] + + # --- 3x3 neighbourhood --- + dem_ng = np.empty(9, dtype=np.float64) + idx = 0 + for di in range(-1, 2): + for dj in range(-1, 2): + dem_ng[idx] = dem[r + di, c + dj] + idx += 1 + altitude = dem_ng[4] + FSI = forest[r, c] if forestBool else 0.0 + + # --- calcDistMin: projected (min_distance) and 3D (minDistXYZ) --- + min_distance = 0.0 + minDistXYZ = 0.0 + if not is_start: + best = 1.0e30 + best3d = 1.0e30 + for p in range(n_parents): + pdir = q_pdir[q_pos, p] + dxp = (pdir % 3) - 1 + dyp = (pdir // 3) - 1 + dd = math.sqrt((dxp * cellsize) ** 2 + (dyp * cellsize) ** 2) + q_pmd[q_pos, p] + if dd < best: + best = dd + if forestBool: + palt = dem[r + dyp, c + dxp] + dz = abs(palt - altitude) + _dy = abs(dyp) * cellsize + # NOTE: replicates flowClass.calcDistMin (calc3D) which uses dy twice + # (upstream quirk); only affects the skipForestDist gate, =0 in tuned configs. + dd3 = math.sqrt(_dy * _dy + _dy * _dy + dz * dz) + q_pmd3d[q_pos, p] + if dd3 < best3d: + best3d = dd3 + min_distance = best + if forestBool: + minDistXYZ = best3d + + # --- calc_z_delta (with forest friction) --- + applyFor = forestBool and (not is_start) and (skipForestDist < minDistXYZ) + if applyFor and (forestModuleCode == _FM_FRICTION or forestModuleCode == _FM_DETRAINMENT) and FSI > 0.0: + if z_delta < noFrictionEffectZDelta: + rest = maxAddedFriction * FSI + slope = (rest - minAddedFriction) / (0.0 - noFrictionEffectZDelta) + friction = max(minAddedFriction, slope * z_delta + rest) + alpha_calc = alpha + max(0.0, friction) + else: + alpha_calc = alpha + minAddedFriction + tan_alpha = math.tan(math.radians(alpha_calc)) + elif applyFor and forestModuleCode == _FM_FRICTIONLAYER: + alpha_for = (alpha + FSI) if frictionLayerRelative else FSI + if alpha_for < alpha: + alpha_for = alpha + tan_alpha = math.tan(math.radians(alpha_for)) + else: + tan_alpha = tan_alpha_base + + z_delta_neighbour = np.empty(9, dtype=np.float64) + for i in range(9): + val = z_delta + (altitude - dem_ng[i]) - ds_cellsize[i] * tan_alpha + if val < 0.0: + val = 0.0 + elif val > max_z_delta: + val = max_z_delta + z_delta_neighbour[i] = val + + # --- calc_persistence (float, no truncation) --- + persistence = np.zeros(9, dtype=np.float64) + no_flow = np.ones(9, dtype=np.float64) + if is_start or first_parent_is_start: + for i in range(9): + persistence[i] = 1.0 + else: + for p in range(n_parents): + pdir = q_pdir[q_pos, p] + mw = q_pzd[q_pos, p] + no_flow[pdir] = 0.0 + dx = (pdir % 3) - 1 + dy = (pdir // 3) - 1 + if dx == -1: + if dy == -1: + persistence[8] += mw + persistence[7] += 0.707 * mw + persistence[5] += 0.707 * mw + if dy == 0: + persistence[5] += mw + persistence[8] += 0.707 * mw + persistence[2] += 0.707 * mw + if dy == 1: + persistence[2] += mw + persistence[1] += 0.707 * mw + persistence[5] += 0.707 * mw + if dx == 0: + if dy == -1: + persistence[7] += mw + persistence[6] += 0.707 * mw + persistence[8] += 0.707 * mw + if dy == 1: + persistence[1] += mw + persistence[0] += 0.707 * mw + persistence[2] += 0.707 * mw + if dx == 1: + if dy == -1: + persistence[6] += mw + persistence[3] += 0.707 * mw + persistence[7] += 0.707 * mw + if dy == 0: + persistence[3] += mw + persistence[0] += 0.707 * mw + persistence[6] += 0.707 * mw + if dy == 1: + persistence[0] += mw + persistence[1] += 0.707 * mw + persistence[3] += 0.707 * mw + for i in range(9): + persistence[i] *= no_flow[i] + + # --- calc_tanbeta / r_t --- + r_t = np.zeros(9, dtype=np.float64) + tan_beta = np.zeros(9, dtype=np.float64) + for i in range(9): + if i == 4 or z_delta_neighbour[i] <= 0.0 or persistence[i] <= 0.0: + tan_beta[i] = 0.0 + else: + beta = math.atan((altitude - dem_ng[i]) / distance[i]) + _HALF_PI + tan_beta[i] = math.tan(beta / 2.0) + tb_sum = 0.0 + for i in range(9): + if tan_beta[i] > 0.0: + tb_sum += tan_beta[i] ** exp + if tb_sum > 0.0: + for i in range(9): + if tan_beta[i] > 0.0: + r_t[i] = tan_beta[i] ** exp / tb_sum + + # --- fp / sl travel angle (non-start) --- + max_gamma = 0.0 + sl_gamma = 0.0 + fluxDep = 0.0 + if not is_start: + dh = altitude_start - altitude + if min_distance > 0.0: + max_gamma = math.atan(dh / min_distance) * _DEG_PER_RAD + sl_dx = abs(col_start - c) + sl_dy = abs(row_start - r) + sl_ds = math.sqrt(sl_dx * sl_dx + sl_dy * sl_dy) * cellsize + if sl_ds > 0.0: + sl_gamma = math.atan(dh / sl_ds) * _DEG_PER_RAD + # forest detrainment reduces flux + if forestBool and forestDetrainmentBool: + rest_d = maxDetrainment * FSI + slope_d = (rest_d - minDetrainment) / (0.0 - noDetrainmentEffectZDelta) + detr = max(minDetrainment, slope_d * z_delta + rest_d) + flux = max(0.0003, flux - detr) + + # --- calc_distribution (fluxDistOldVersion=False default) --- + dist = np.zeros(9, dtype=np.float64) + rt_sum = 0.0 + for i in range(9): + rt_sum += r_t[i] + if rt_sum > 0.0: + pr_sum = 0.0 + for i in range(9): + pr_sum += persistence[i] * r_t[i] + if pr_sum > 0.0: + for i in range(9): + dist[i] = persistence[i] * r_t[i] / pr_sum * flux + + if fluxDistOldVersion: + count = 0 + for i in range(9): + if 0.0 < dist[i] < flux_threshold: + count += 1 + else: + count = 0 + for i in range(9): + if dist[i] >= flux_threshold: + count += 1 + mass_below = 0.0 + for i in range(9): + if dist[i] < flux_threshold: + mass_below += dist[i] + if mass_below > 0.0 and count > 0: + add = mass_below / count + for i in range(9): + if dist[i] >= flux_threshold: + dist[i] += add + elif dist[i] < flux_threshold: + dist[i] = 0.0 + dist_sum = 0.0 + for i in range(9): + dist_sum += dist[i] + if dist_sum != flux and count > 0: + corr = (flux - dist_sum) / count + for i in range(9): + if dist[i] >= flux_threshold: + dist[i] += corr + if count == 0: + fluxDep = flux + + # --- collect children (dist >= threshold), sort ascending (z_delta,flux,row,col) --- + ch_r = np.empty(9, dtype=np.int64) + ch_c = np.empty(9, dtype=np.int64) + ch_flux = np.empty(9, dtype=np.float64) + ch_zd = np.empty(9, dtype=np.float64) + nch = 0 + for i in range(9): + if dist[i] >= flux_threshold: + ch_r[nch] = r - 1 + i // 3 + ch_c[nch] = c - 1 + i % 3 + ch_flux[nch] = dist[i] + ch_zd[nch] = z_delta_neighbour[i] + nch += 1 + for i in range(nch - 1): + m = i + for j in range(i + 1, nch): + sw = False + if ch_zd[j] < ch_zd[m]: + sw = True + elif ch_zd[j] == ch_zd[m]: + if ch_flux[j] < ch_flux[m]: + sw = True + elif ch_flux[j] == ch_flux[m]: + if ch_r[j] < ch_r[m]: + sw = True + elif ch_r[j] == ch_r[m] and ch_c[j] < ch_c[m]: + sw = True + if sw: + m = j + if m != i: + ch_zd[i], ch_zd[m] = ch_zd[m], ch_zd[i] + ch_flux[i], ch_flux[m] = ch_flux[m], ch_flux[i] + ch_r[i], ch_r[m] = ch_r[m], ch_r[i] + ch_c[i], ch_c[m] = ch_c[m], ch_c[i] + + # --- dedup vs pending / append new children --- + for k in range(nch): + cr = ch_r[k] + cc = ch_c[k] + pq = pending_qidx[cr, cc] + dxp = c - cc + dyp = r - cr + pdir_idx = (dyp + 1) * 3 + (dxp + 1) + if pq >= 0: + q_flux[pq] += ch_flux[k] + if ch_zd[k] > q_zdelta[pq]: + q_zdelta[pq] = ch_zd[k] + np_i = q_n_parents[pq] + if np_i < q_pdir.shape[1]: + q_pdir[pq, np_i] = pdir_idx + q_pzd[pq, np_i] = z_delta + q_pmd[pq, np_i] = min_distance + q_pmd3d[pq, np_i] = minDistXYZ + q_n_parents[pq] = np_i + 1 + if forestInteraction: + child_isforest = q_isforest[pq] + if fic < (q_fic[pq] - child_isforest): + q_fic[pq] = fic + child_isforest + else: + if cr < 1 or cr >= H - 1 or cc < 1 or cc >= W - 1: + continue + nd = False + for di in range(-1, 2): + for dj in range(-1, 2): + if dem[cr + di, cc + dj] == nodata: + nd = True + if nd: + continue + if q_tail >= q_r.shape[0]: + return -1 # queue overflow -> caller grows the workspace and retries + new_pos = q_tail + q_r[new_pos] = cr + q_c[new_pos] = cc + q_zdelta[new_pos] = ch_zd[k] + q_flux[new_pos] = ch_flux[k] + q_is_start[new_pos] = False + q_first_parent_start[new_pos] = is_start + q_mindist[new_pos] = 0.0 + q_mindist3d[new_pos] = 0.0 + cif = 1 if (forestInteraction and forest[cr, cc] > 0) else 0 + q_isforest[new_pos] = cif + q_fic[new_pos] = cif + fic + q_n_parents[new_pos] = 1 + q_pdir[new_pos, 0] = pdir_idx + q_pzd[new_pos, 0] = z_delta + q_pmd[new_pos, 0] = min_distance + q_pmd3d[new_pos, 0] = minDistXYZ + pending_qidx[cr, cc] = new_pos + modified_r[n_modified] = cr + modified_c[n_modified] = cc + n_modified += 1 + q_tail += 1 + + # --- accumulate outputs for the processed cell --- + if z_delta > zDeltaArray[r, c]: + zDeltaArray[r, c] = z_delta + if flux > fluxArray[r, c]: + fluxArray[r, c] = flux + routFluxSumArray[r, c] += flux + depFluxSumArray[r, c] += fluxDep + if z_delta > zDeltaPathArray[r, c]: + zDeltaPathArray[r, c] = z_delta + if max_gamma > fpMaxArray[r, c]: + fpMaxArray[r, c] = max_gamma + if fpMinArray[r, c] >= 0.0 and max_gamma >= 0.0: + if max_gamma < fpMinArray[r, c]: + fpMinArray[r, c] = max_gamma + else: + if max_gamma > fpMinArray[r, c]: + fpMinArray[r, c] = max_gamma + if sl_gamma > slArray[r, c]: + slArray[r, c] = sl_gamma + if min_distance > travelMaxArray[r, c]: + travelMaxArray[r, c] = min_distance + if travelMinArray[r, c] >= 0.0 and min_distance >= 0.0: + if min_distance < travelMinArray[r, c]: + travelMinArray[r, c] = min_distance + else: + if min_distance > travelMinArray[r, c]: + travelMinArray[r, c] = min_distance + if forestInteraction: + if forestIntArray[r, c] >= 0.0 and fic >= 0.0: + if fic < forestIntArray[r, c]: + forestIntArray[r, c] = fic + else: + if fic > forestIntArray[r, c]: + forestIntArray[r, c] = fic + + # finalize once per distinct cell: fold this path's max zDelta into zDeltaSum, + # then reset per-BFS workspace (zDeltaPath, visited, pending) via the modified list + for i in range(n_modified): + rr = modified_r[i] + cc = modified_c[i] + if visited[rr, cc] == 1: + zDeltaSumArray[rr, cc] += zDeltaPathArray[rr, cc] + zDeltaPathArray[rr, cc] = 0.0 + visited[rr, cc] = 0 + pending_qidx[rr, cc] = -1 + + return q_tail + + +def _forest_scalars(forestBool, forestParams): + """Translate the forestParams dict into scalars/flags for the njit kernel.""" + if not forestBool or forestParams is None: + return (_FM_NONE, False, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, False, False, 0.0) + + module = forestParams["forestModule"] + if module == "forestFriction": + code = _FM_FRICTION + elif module == "forestDetrainment": + code = _FM_DETRAINMENT + elif module == "forestFrictionLayer": + code = _FM_FRICTIONLAYER + else: + code = _FM_NONE + + forestInteraction = bool(forestParams.get("forestInteraction", False)) + maxFr = float(forestParams.get("maxAddedFriction", 0.0)) + minFr = float(forestParams.get("minAddedFriction", 0.0)) + vThFr = float(forestParams.get("velThForFriction", 0.0)) + maxDe = float(forestParams.get("maxDetrainment", 0.0)) + minDe = float(forestParams.get("minDetrainment", 0.0)) + vThDe = float(forestParams.get("velThForDetrain", 0.0)) + skipForestDist = float(forestParams.get("skipForestDist", 0.0)) + fFrLayerType = forestParams.get("fFrLayerType", "absolute") + + noFrZ = (vThFr * vThFr) / (_SQRT2 * _G) if vThFr != 0.0 else 0.0 + noDeZ = (vThDe * vThDe) / (_SQRT2 * _G) if vThDe != 0.0 else 0.0 + + # matches flowClass: detrainment only for the forestDetrainment module with non-zero params + if module in ("forestFriction", "forestFrictionLayer"): + detrainBool = False + elif maxDe == 0.0 and minDe == 0.0 and vThDe == 0.0: + detrainBool = False + else: + detrainBool = True + + layerRel = (fFrLayerType == "relative") + return (code, forestInteraction, maxFr, minFr, noFrZ, maxDe, minDe, noDeZ, + detrainBool, layerRel, skipForestDist) + + +def calculationNumba(args): + """Numba drop-in for flowCore.calculation(): same args, same 14-element return. + + Processes one release-chunk (one Pool task). The per-release-pixel BFS runs in + the compiled kernel in double precision (reproducing the Python engine); + variable alpha / max_z / exponent are resolved per release pixel (constant + along each path), matching flowCore.calculation(). Output rasters are float32. + """ + dem_in = args[0] + release = args[2] + alpha0 = float(args[3]) + exp0 = float(args[4]) + flux_threshold = float(args[5]) + max_z0 = float(args[6]) + nodata = float(args[7]) + cellsize = float(args[8]) + forestBool = args[10] + varParams = args[11] + fluxDistOldVersionBool = bool(args[12]) + + varUmaxBool = varParams["varUmaxBool"] + varUmaxArray = varParams["varUmaxArray"] + varAlphaBool = varParams["varAlphaBool"] + varAlphaArray = varParams["varAlphaArray"] + varExponentBool = varParams["varExponentBool"] + varExponentArray = varParams["varExponentArray"] + + forestArray = args[14] if forestBool else None + forestParams = args[15] if forestBool else None + (fmCode, forestInteraction, maxFr, minFr, noFrZ, maxDe, minDe, noDeZ, + detrainBool, layerRel, skipForestDist) = _forest_scalars(forestBool, forestParams) + + dem = np.ascontiguousarray(dem_in, dtype=np.float64) + if forestArray is not None: + forest = np.ascontiguousarray(forestArray, dtype=np.float64) + else: + forest = np.zeros_like(dem) + H, W = dem.shape + + # release start pixels, in the same order as flowCore.calculation() + rel = release.copy() + rel[rel < 0] = 0 + rel[rel == nodata] = 0 + rel[rel > 0] = 1 + row_list, col_list = get_start_idx(dem, rel) + + ds_cellsize = (_DS * cellsize).astype(np.float64) + distance = (_DS_TANBETA * cellsize).astype(np.float64) + + # The per-BFS queue workspace starts modest and only grows if a single release + # pixel's flow genuinely needs more room; on overflow the whole chunk is re-run + # with a 4x-larger queue (grow-and-retry). This keeps memory small for typical + # runs while guaranteeing a path is never silently truncated. (Peak queue depth + # observed on a 5 m long-runout tile was ~4k, well under the 131072 start.) + MAXP = 8 + MAXQ = 1 << 17 # 131072 + while True: + # output arrays — dtypes/init identical to flowCore.calculation() + zDeltaArray = np.zeros((H, W), dtype=np.float32) + zDeltaSumArray = np.zeros((H, W), dtype=np.float32) + routFluxSumArray = np.zeros((H, W), dtype=np.float32) + depFluxSumArray = np.zeros((H, W), dtype=np.float32) + fluxArray = np.ones((H, W), dtype=np.float32) * -9999 + countArray = np.zeros((H, W), dtype=np.int32) + fpMaxArray = np.ones((H, W), dtype=np.float32) * -9999 + fpMinArray = np.ones((H, W), dtype=np.float32) * -9999 + slArray = np.ones((H, W), dtype=np.float32) * -9999 + travelMaxArray = np.ones((H, W), dtype=np.float32) * -9999 + travelMinArray = np.ones((H, W), dtype=np.float32) * -9999 + forestIntArray = np.ones((H, W), dtype=np.float32) * -9999 + zDeltaPathArray = np.zeros((H, W), dtype=np.float32) + pending_qidx = np.full((H, W), -1, dtype=np.int64) + visited = np.zeros((H, W), dtype=np.int8) + + q_r = np.empty(MAXQ, dtype=np.int64) + q_c = np.empty(MAXQ, dtype=np.int64) + q_zdelta = np.empty(MAXQ, dtype=np.float64) + q_flux = np.empty(MAXQ, dtype=np.float64) + q_is_start = np.empty(MAXQ, dtype=np.bool_) + q_first_parent_start = np.empty(MAXQ, dtype=np.bool_) + q_mindist = np.empty(MAXQ, dtype=np.float64) + q_mindist3d = np.empty(MAXQ, dtype=np.float64) + q_fic = np.empty(MAXQ, dtype=np.float64) + q_isforest = np.empty(MAXQ, dtype=np.float64) + q_n_parents = np.empty(MAXQ, dtype=np.int64) + q_pdir = np.empty((MAXQ, MAXP), dtype=np.int64) + q_pzd = np.empty((MAXQ, MAXP), dtype=np.float64) + q_pmd = np.empty((MAXQ, MAXP), dtype=np.float64) + q_pmd3d = np.empty((MAXQ, MAXP), dtype=np.float64) + modified_r = np.empty(MAXQ, dtype=np.int64) + modified_c = np.empty(MAXQ, dtype=np.int64) + + overflow = False + for k in range(len(row_list)): + rIdx = int(row_list[k]) + cIdx = int(col_list[k]) + alpha = alpha0 + max_z = max_z0 + exp = exp0 + if varUmaxBool and varUmaxArray is not None: + v = varUmaxArray[rIdx, cIdx] + if 0 < v <= 8848: + max_z = float(v) + if varAlphaBool and varAlphaArray is not None: + v = varAlphaArray[rIdx, cIdx] + if 0 < v <= 90: + alpha = float(v) + if varExponentBool and varExponentArray is not None: + v = varExponentArray[rIdx, cIdx] + if v > 0: + exp = float(v) + + qt = _bfs_single(dem, forest, H, W, nodata, rIdx, cIdx, + cellsize, alpha, exp, flux_threshold, max_z, + ds_cellsize, distance, + forestBool, fmCode, forestInteraction, + maxFr, minFr, noFrZ, maxDe, minDe, noDeZ, + detrainBool, layerRel, skipForestDist, + fluxDistOldVersionBool, + zDeltaArray, fluxArray, countArray, zDeltaSumArray, zDeltaPathArray, + routFluxSumArray, depFluxSumArray, + fpMaxArray, fpMinArray, slArray, + travelMaxArray, travelMinArray, forestIntArray, + pending_qidx, visited, + q_r, q_c, q_zdelta, q_flux, q_is_start, q_first_parent_start, + q_mindist, q_mindist3d, q_fic, q_isforest, + q_n_parents, q_pdir, q_pzd, q_pmd, q_pmd3d, + modified_r, modified_c) + if qt < 0: # queue overflow -> grow workspace and re-run the whole chunk + overflow = True + break + + if not overflow: + break + MAXQ *= 4 + + backcalc = None + # 14-element tuple matching flowCore.calculation(): res[12]=relId startcell dict + # (not produced by the numba engine — relId outputs fall back to the Python + # engine in run()), res[13]=forestInteraction array (None if not requested). + startCellIdDict = None + forestIntOut = forestIntArray if forestInteraction else None + return (zDeltaArray, fluxArray, countArray, zDeltaSumArray, backcalc, + fpMaxArray, slArray, travelMaxArray, travelMinArray, fpMinArray, + routFluxSumArray, depFluxSumArray, startCellIdDict, forestIntOut) diff --git a/avaframe/runStandardTestsCom4FlowPy.py b/avaframe/runStandardTestsCom4FlowPy.py index 6edebc504..93f842c96 100644 --- a/avaframe/runStandardTestsCom4FlowPy.py +++ b/avaframe/runStandardTestsCom4FlowPy.py @@ -22,6 +22,13 @@ from avaframe.in3Utils import logUtils import avaframe.in2Trans.rasterUtils as rasterUtils +def _checkNumbaInstalled() -> bool: + try: + import numba + return True + except ImportError: + return False + def compareRasters(path, pathRef): """ @@ -80,9 +87,18 @@ def main(): # filter benchmarks for tag standardTest filterType = 'TAGS' + valuesList = ['standardTest', 'com4FlowPy'] # looking for 'com4FlowPy' and 'standardTest' in TAGS list - testList = tU.filterBenchmarks(testDictList, filterType, valuesList, condition='and') + testListAll = tU.filterBenchmarks(testDictList, filterType, valuesList, condition='and') + + valuesList = ['standardTest', 'com4FlowPy', 'numba'] + testListNumba = tU.filterBenchmarks(testDictList, filterType, valuesList, condition='and') + + if _checkNumbaInstalled(): + testList = testListAll + else: + testList = [item for item in testListAll if item not in testListNumba] # Set directory for full standard test report outDir = _avaframeDir / 'tests' / 'reportsCom4FlowPy' @@ -93,13 +109,16 @@ def main(): _startDate = datetime.now() with open(reportFile, 'w') as pfile: - # Write header pfile.write('# Standard Tests Report \n\n') pfile.write('Comparing __com4FlowPy__ simulations to selected benchmark results \n\n') - pfile.write(f'__tests started__ : {_startDate}\n\n') + pfile.write('* * * \n') + if _checkNumbaInstalled(): + pfile.write('`numba` __found__: running all tests for python and numba engines ✓') + else: + pfile.write('`numba` __NOT found__: skipping tests for numba engine ✗') + pfile.write('\n* * * \n') - log = logUtils.initiateLogger(outDir, logName) log.info('The following benchmark tests will be fetched ') @@ -108,9 +127,10 @@ def main(): for test in testList: pfile.write(f"- {test['NAME']}\n") log.info('%s' % test['NAME']) - pfile.write('\n* * * \n') + pfile.write('\n') + pfile.write(f'__tests started__ : {_startDate}\n\n') + pfile.write('* * * \n') - # create a temporary directory, where the outputs of all standard Tests are stored # clean-up is automatic - this way we don't pollute the avaframe/data/ directory with tempfile.TemporaryDirectory(prefix="avaframe_stdTests_") as tempDir: @@ -145,14 +165,16 @@ def main(): pfile.write("|Model Output|Result of comparison|status\n") pfile.write("|----:|:-----:|:---:|\n") - avaDir = test['AVADIR'] + # define avaDir relative to _avaframeDir to allow execution of this script from different locations + # not just AvaFrame/avaframe directory + avaDir = str( _avaframeDir / pathlib.Path(test['AVADIR']) ) cfgMain['MAIN']['avalancheDir'] = avaDir # Fetch benchmark test info refDir = pathlib.Path(_avaframeDir, '..', 'benchmarks', test['NAME']) # Clean input directory(ies) of old work and output files - initProj.cleanSingleAvaDir(avaDir, deleteOutput=False) + initProj.cleanSingleAvaDir(_avaframeDir / avaDir, deleteOutput=False) # Load input parameters from configuration file for standard tests benchmarkCfg = refDir / ('%s' % test['INI']) @@ -165,7 +187,11 @@ def main(): avalancheDir = cfgMain["MAIN"]["avalancheDir"] cfgPath = readFlowPyinputs(avalancheDir, cfg, log) - compDir = tmpTestsDir / pathlib.Path(avalancheDir) + # for the temporary output folder we cannot use the full path, but just the relative part of the path + # following the _avaframeDir path + avaDirTempOutPut = test['AVADIR'] + # compDir = output location of outputs generated for each test within the temporary directory + compDir = tmpTestsDir / pathlib.Path(avaDirTempOutPut) cfgPath["customDirs"] = False cfgPath["resDir"] = compDir diff --git a/avaframe/tests/test_com4FlowPyNumba.py b/avaframe/tests/test_com4FlowPyNumba.py new file mode 100644 index 000000000..dd38e30c7 --- /dev/null +++ b/avaframe/tests/test_com4FlowPyNumba.py @@ -0,0 +1,125 @@ +""" +Pytest for the com4FlowPy numba compute engine (flowCoreNumba). + +Verifies that the numba engine (engine = numba) reproduces the reference Python +(Cell-based) engine, flowCore.calculation(), bit-for-bit on consistently-sloped +synthetic terrain (no flat-terrain routing ties, which are numerically +knife-edge in both engines). +""" +import numpy as np +import pytest + +import avaframe.com4FlowPy.flowCore as flowCore + +# the whole module is skipped if the optional 'numba' dependency is absent +pytest.importorskip("numba") +from avaframe.com4FlowPy import flowCoreNumba # noqa: E402 + + +def _make_dem(ny=30, nx=30): + """A consistently down-sloping DEM (downhill in +y) with gentle cross-valley + curvature and strictly varying values, so flow routing has no flat-terrain + ties (which are numerically knife-edge in both engines).""" + yy, xx = np.mgrid[0:ny, 0:nx].astype(np.float64) + return 500.0 - 7.0 * yy + 0.15 * (xx - nx / 2.0) ** 2 + 0.013 * xx + + +def _default_var(): + return { + "varUmaxBool": False, "varUmaxArray": None, + "varAlphaBool": False, "varAlphaArray": None, + "varExponentBool": False, "varExponentArray": None, + } + + +def _args(dem, pra, *, alpha=25.0, forestBool=False, forestParams=None, + forestArray=None, varParams=None, outputs=None): + if varParams is None: + varParams = _default_var() + if outputs is None: + outputs = ["zDelta", "flux", "cellCounts", "zDeltaSum", "fpTravelAngleMax", + "fpTravelAngleMin", "slTravelAngle", "travelLengthMax", + "travelLengthMin", "routFluxSum", "depFluxSum"] + relOutputParams = {"relIdBool": False, "relIdArray": None, + "relVolBool": False, "relVolArray": None} + return [dem, None, pra, alpha, 8, 3e-4, 270, -9999, 10.0, False, forestBool, + varParams, False, False, forestArray, forestParams, outputs, relOutputParams] + + +# comparable result indices (4=backcalc and 12=relId are None on this path; +# 13=forestInteraction is compared separately when active) +_IDX = {0: "zDelta", 1: "flux", 2: "cellCounts", 3: "zDeltaSum", + 5: "fpTravelAngleMax", 6: "slTravelAngle", 7: "travelLengthMax", + 8: "travelLengthMin", 9: "fpTravelAngleMin", 10: "routFluxSum", + 11: "depFluxSum"} + + +def _compare(py, nb): + assert len(nb) == len(py) == 14 + for i, name in _IDX.items(): + a = np.asarray(py[i], dtype=np.float64) + b = np.asarray(nb[i], dtype=np.float64) + assert np.array_equal(a, b), \ + f"{name}: {int((a != b).sum())} cells differ, max|Δ|={np.abs(a - b).max()}" + + +def _pra(dem): + pra = np.zeros_like(dem, dtype=np.int32) + pra[1, dem.shape[1] // 2] = 1 # single release near the top of the slope + return pra + + +def test_numba_matches_python_noforest(): + dem = _make_dem() + args = _args(dem, _pra(dem)) + _compare(flowCore.calculation(args), flowCoreNumba.calculationNumba(args)) + + +def test_numba_matches_python_forestFriction(): + dem = _make_dem() + forest = np.zeros_like(dem) + forest[10:20, :] = 0.6 + fp = {"forestModule": "forestFriction", "maxAddedFriction": 20.0, + "minAddedFriction": 2.0, "velThForFriction": 30.0, "maxDetrainment": 0.0, + "minDetrainment": 0.0, "velThForDetrain": 0.0, "fFrLayerType": "absolute", + "skipForestDist": 0.0, "forestInteraction": True} + args = _args(dem, _pra(dem), forestBool=True, forestParams=fp, forestArray=forest) + py = flowCore.calculation(args) + nb = flowCoreNumba.calculationNumba(args) + _compare(py, nb) + assert np.array_equal(np.asarray(py[13], float), np.asarray(nb[13], float)), \ + "forestInteraction array mismatch" + + +def test_numba_matches_python_forestDetrainment(): + dem = _make_dem() + forest = np.zeros_like(dem) + forest[8:22, :] = 0.5 + fp = {"forestModule": "forestDetrainment", "maxAddedFriction": 52.0, + "minAddedFriction": 5.0, "velThForFriction": 270.0, "maxDetrainment": 0.003, + "minDetrainment": 0.00001, "velThForDetrain": 270.0, "fFrLayerType": "absolute", + "skipForestDist": 0.0, "forestInteraction": True} + args = _args(dem, _pra(dem), forestBool=True, forestParams=fp, forestArray=forest) + py = flowCore.calculation(args) + nb = flowCoreNumba.calculationNumba(args) + _compare(py, nb) + assert np.array_equal(np.asarray(py[13], float), np.asarray(nb[13], float)) + + +def test_numba_matches_python_variableAlphaUmax(): + dem = _make_dem() + var = _default_var() + var["varAlphaBool"] = True + var["varAlphaArray"] = np.full_like(dem, 22.0) # per-cell alpha (deg) + var["varUmaxBool"] = True + var["varUmaxArray"] = np.full_like(dem, 150.0) # per-cell zDeltaLim (m) + args = _args(dem, _pra(dem), varParams=var) + _compare(flowCore.calculation(args), flowCoreNumba.calculationNumba(args)) + + +if __name__ == "__main__": + test_numba_matches_python_noforest() + test_numba_matches_python_forestFriction() + test_numba_matches_python_forestDetrainment() + test_numba_matches_python_variableAlphaUmax() + print("all numba-engine equivalence tests passed") diff --git a/benchmarks/com4_avaArzlerAlmForestNumba/ArzlerAlmForestNumba_com4FlowPyCfg.ini b/benchmarks/com4_avaArzlerAlmForestNumba/ArzlerAlmForestNumba_com4FlowPyCfg.ini new file mode 100644 index 000000000..4ef7a7186 --- /dev/null +++ b/benchmarks/com4_avaArzlerAlmForestNumba/ArzlerAlmForestNumba_com4FlowPyCfg.ini @@ -0,0 +1,308 @@ +### Config File - This file contains the main settings for the com4FlowPy run +## Set your parameters +# This file will be overridden by local_com4FlowPyCfg.ini if it exists +# So copy this file to local_com4FlowPyCfg.ini, adjust your variables there + +# Optional settings------------------------------- +[GENERAL] +#++++++++++++ Flow-Py Model Parameters +# alpha: +# Angle-of-reach Alpha [°] - defines max. longitudinal runout limit +# equivalent to a Coulomb-friction of tan(alpha) in a sliding block model +#--------------------- +# exp: +# Spreading Coefficient (influences lateral spreading) +#--------------------- +# flux_threshold: +# Flux threshold (influences lateral spreading) +#--------------------- +# max_z: +# Energy-Line-Height Limit (can be interpreted as velocity limit) +# max_v = sqrt(max_z*19.62) +# typical values: +# - Avalanche: ~270 m ... ~72 m/s +# - Rockfall: ~130 m ... ~50 m/s +# - Debris-flow: ~12 m ... ~15 m/s +#--------------------- + +alpha = 25 +exp = 8 +flux_threshold = 3.0e-4 +max_z = 8848 + +#++++++++++++ Use Infrastructure + +infra = False + +#++++++++++++ preview Mode +# if previewMode = True, not every release cell is processend independently +# if a releaseCell is already "hit"/"affected" by a prior calculated path +# then the processing of this release cell is skipped +# can be used for a faster preview of model results (e.g. for checking of input parameters) or +# to save calculation time (e.g. when calculating with infrastructure) +# NOTE: results will deviate from "normal" model run and all outputLayers relying on the +# summed/combined output of different paths will not provide sensible values!!! + +previewMode = False + + +#++++++++++++ Use a dynamic u_max Limit +# Requires an additional tif-file containing the uMax (in m/s) +# or zDeltaMax (m) Limits +# in every cell where a release cell is. +# the paths for each release cells are calculated with these uMax values +# (computed to z_delta, similar to the max_z) or zDelta values. +# In varUmaxParameter the parameter (uMax in m/s or zDeltaMax in m) +# provided is given. +variableUmaxLim = False +varUmaxParameter = uMax + +#++++++++++++ Use a dynamic alpha angle +# Requires an additional tif-file containing the alpha angles +# in every cell where a release cell is. +# The paths for each release cells are calculated with these alpha angles +variableAlpha = False + +#++++++++++++ Use a dynamic exponent +# Requires an additional tif-file containing the exponents +# in every cell where a release cell is. +# The paths for each release cells are calculated with these exponent +variableExponent = False + +#++++++++++++ Use Forest Information +# NOTE AH 20240408: Dummy settings for Forest Interaction +# The forest implementation should mimick/reproduce the +# one from the "foreste_detrainment" branch in the avaframe/FlowPy repo +# which is partly described in D'Amboise et al. (2021) +# +# Forest-Interaction is only used if 'forest' is set to 'True' +#++++++++++++ +forest = True + +#++++++++++++ +# Type of Forest-Interaction Model to use (only works in conjunction with forest = True!) +# valid choices: +# +# * 'forestFriction' (D'Amboise et al., 2021) - use added friction on forest Pixels +# forest-Raster (set in 'forestPath') is assumed to have values between +# 0 (no Forest) and 1 (full protective function). +# model behavior is governed by 'maxAddedFrictionFor', 'minAddedFrictionFor' +# and 'velThForFriction' +# +# * 'forestDetrainment' (D'Amboise et al., 2022?) - added friction and detrainment on forested Pixels +# forest-Raster (set in 'forestPath') is assumed to have values between +# 0 (no Forest) and 1 (full protective function). +# model behavior is governed by 'maxAddedFrictionFor', 'minAddedFrictionFor' +# and 'velThForFriction', 'maxDetrainmentFor', 'minDetrainmentFor', 'velThForDetrain' +# +# * 'forestFrictionLayer' - provide Layer with local alpha / local delta-Alpha on forested pixels / for all pixels in +# the domain +# forest-Raster (set in 'forestPath') is assumed to have values between 0° (no friction) and 90° (inf. friciton) +# forest-Raster is either interpreted as 'absolute' or 'difference' +# 'absolute' - absolute 'alpha' values are used on forested pixels (all pixels with values > 0) +# 'relative' - 'alpha' values are interpreted as additional friction on forested pixels (all pixels with +# values > 0) +#++++++++++++ + +# ['forestFriction', 'forestDetrainment', 'forestFrictionLayer'] +forestModule = forestFriction +# ForestInteraction gives an additional output layer for the amount of forested cells +# a path ran through +# (now implemented, only when forest = True) -> TODO: we dont require a forestModule to compute +# for ForestInteraction Layer, but we need to read in the Forest - Layer (requires forest = True) +forestInteraction = False + +#++++++++++++ Forest added Friction +# These are the parameters for the "add friction on forested cells"-approach +# described in D'Amboise et al. (2021) +# +# maxAddedFrictionFor [°]: maximum increase of Alpha angle/basal friction on forested cells +# minAddedFrictionFor [°]: minimum increase of Alpha angle/basal friction on forested cells +# velThForFriction [m/s]: velocity threshold for forest-friction effect +#++++++++++++ + +maxAddedFrictionFor = 10 +minAddedFrictionFor = 2 +velThForFriction = 30 + +#++++++++++++ Forest Detrainment +# These are the parameters for the "forestDetrainment"-approach - D'Amboise et al. (2022??) +# The idea is to remove 'virtual mass' (aka 'flux' in Flow-Py) from the modeled flow/process +# on forested pixels in dependency of FSI and local z_delta +# +# maxDetrainmentFor [°]: default set to '0' (i.e. no detrainment) - foreste_detrainment val: +# minDetrainmentFor [°]: default set to '0' (i.e. no detrainment) - foreste_detrainment val: +# velThForDetrain [m/s]: velocity threshold for forest-detraiment effect # default set to '0' (i.e. no detrainment) - foreste_detrainment val: +#++++++++++++ + +maxDetrainmentFor = 0 +minDetrainmentFor = 0 +velThForDetrain = 0 + +#++++++++++++ ForestFrictionLayer +# The forest layer defines the alpha angle in forested areas. +# If forestFrictionLayerType = relative, forest layer values are added to the global alpha angle. +# e.g. if alpha = 24° and value on forest cell is 5°, then effective local alpha of 24° + 5° = 29° will be used. +# ------------------- +# If forestFrictionLayerType = absolute, forest layer values are used as alpha angle, +# if they are larger than the global alpha angle. +# e.g. if global alpha = 24° and value provided in the forestFrictionLayer on a cell is 35°, +# then local alpha of 35° will be used. +#++++++++++++ + +# ['absolute', 'relative'] +forestFrictionLayerType = absolute + +# skip Forest Effect (added forest friction) for first x meters (calculated in 3D - XYZ) +# should allow an initial acceleration phase of processes starting in or directly above +# dense forests (these would in many cases stop otherwise) +# if e.g. skipForestDist = 40, no added forestFriction will be assumed until 40 m 3D-distance +# along the path from the startCell. +skipForestDist = 0 + +#++++++++++++ Method to calculate flux distribution +# We fixed a bug in flowClass.py, which affects the distribution of the remaining flux, +# if a cell receives flux smaller than the provided flux_threshold. +# +# The default now (post Jan. 2025) is a calculation with the fixed bug! +# +# For backward compatibility the old version (prior to Jan. 2025 - with minor bug) can +# be switched on by setting "fluxDistOldVersion = True". + +fluxDistOldVersion = False + +# compute engine for the per-cell BFS. The 'numba' engine requires the optional +# 'numba' dependency (JIT-compiled kernel); it reproduces the 'python' engine's +# results and is much faster. infra/back-calculation, previewMode and relId +# outputs automatically fall back to the python engine. +# python : default, Cell-based reference implementation +# numba : JIT-compiled kernel (double precision), bit-for-bit vs 'python' +engine = numba + +#++++++++++++ Parameters for Tiling +# tileSize: size of tiles in x and y direction in meters (if total size of) x +# or y of input DEM is larger than tileSize, then the input raster +# layers are tiles +# tileOverlap: overlap between single tiles in m (5 km is rather conservative) +#--------------------------------------- + +tileSize = 15000 +tileOverlap = 5000 + +#++++++++++++ Parameters for CPU usage/multiprocessing +# Recommended to leave at default values unless performance has to be +# tweaked/optimized (e.g. for application on dedicated machines, maxing out +# computational resources) +# +# procPerCPUCore: number of processer per core for multiprocessing +# default value = 1 (it is recommended not to change this value +# unless you want to max. out CPUs, which might come at the cost +# of RAM comsumption!!!) +# chunkSize: tbd. +# maxChunks: Maximum number of single tasks for multiprocessing, if set too high +# this can max out RAM and lead to unexpected behavior (in this case +# probably an infinite loop --> you see it in the .log if this happens) +# On machines with a lot of CPUs and RAM this value might be set higher ... +#++++++++++++ +procPerCPUCore = 1 +chunkSize = 50 +maxChunks = 500 + +# Optional Custom Paths +[PATHS] + +# define format of output raster files (default = .tif) +# available options: [.tif, .asc] +# if you plan to utilize avaFrame tools to analyse result rasters +# then you should choose '.asc' here! +outputFileFormat = .tif + +# use LZW compression when writing tif raster files +useCompression = True + +# define noData value that is assigned in output rasters in cells that are not affected by the process (default: -9999) +# when changing this value BE AWARE that noData values can have same values as a affected cells +# (e.g., when using outputNoDataValue = 0) +outputNoDataValue = -9999 + +# define the different output files that are written to disk +# default = 'zDelta|cellCounts|travelLengthMax|fpTravelAngleMax' +# additional options: +# slTravelAngle +# flux +# zDeltaSum +# routFluxSum +# depFluxSum +# travelLengthMin +# fpTravelAngleMin +# relIdPolygon +# relIdCount +# if forestInteraction: forestInteraction is automatically added to outputs +# if infra: backCalculation is automatically added to output +# if relVolMin or relVolMax is in outputFiles, the Volume of the PRA should be provided in the raster file in the REL folder +# if relIdCount or relIdPolygon is in outputFiles, the ids of the PRAs should be provided in the raster file in the RELID folder +outputFiles = zDelta|cellCounts|travelLengthMax|fpTravelAngleMax + +# whether simulation results with the same simHash of the running simulation already exists +# AND the existing resultsFolder already contains valid com4FlowPy outputs. +# 1) overwriteResults = default ... does not re-run a simulation if results folder (res_) and .json +# -- if there are remnants from a previously attempted but not succesfully finished simulation (e.g. existing res_ and/or .json) +# --> delete existing results folder (res_) and .json and run Simulation +# 2) overwriteResults = reRunAndOverwrite ... deletes existing results folder (res_) and .json and runs Simulation +# 3) overwriteResults = reRunAndBackup ... moves existing result folder (res_) and .json to a dedicated backup folder (e.g. in /BACKUP) + +overwriteResults = reRunAndOverwrite + +#++++++++++++ Custom paths True/False +# default: False +# if set to 'False': +# the default AvaFrame Path/Directory structure is used +# i.e. the Work- and Result-Folders and inputFiles are taken from cfgMain and +# Input-Files are expected in AvaFrame directory structure +# if set to 'True': +# workDir, demPath, releasePath and [infraPath <-> if 'infra'==True] +# can be set here +# +#++++++++++++ +useCustomPaths = False +# for now only works with 'useCustomPaths=True'; if 'True' the temp folder is +# deleted after the output-files have been successfully written to disk +deleteTempFolder = False + +# if useCustomPathDEM set to True: only DEM path can set here, the other input data are used from +# avalancheDir (AvaFrame folder structure) (works only if useCustomPaths = False) +useCustomPathDEM = False + +# if 'useCustomPaths = True', then these input and output paths are used +# the output Folder is automatically placed inside the 'workDir' +workDir = +demPath = +releasePath = +relIdPath = +infraPath = +forestPath = +varUmaxPath = +varAlphaPath = +varExponentPath = + +# plot save results flags------------------------ +# NOTE-TODO: These Flags still don't do anything, +# would be cool to have full integration for +# com4FlowPy with avaframe post-processing utils ... +[FLAGS] +# Plot the avalanche path and DEM figure +plotPath = False + +# Plot the avalanche profile with alpha beta points +plotProfile = False + +# Save the profile figure +saveProfile = True + +# Write results to txt file +writeRes = True + +# keep intermediate results +fullOut = False + +#---------------------------------------------------------- diff --git a/benchmarks/com4_avaArzlerAlmForestNumba/com4_avaArzlerAlmForest_desDict.json b/benchmarks/com4_avaArzlerAlmForestNumba/com4_avaArzlerAlmForest_desDict.json new file mode 100644 index 000000000..6262c68b3 --- /dev/null +++ b/benchmarks/com4_avaArzlerAlmForestNumba/com4_avaArzlerAlmForest_desDict.json @@ -0,0 +1,10 @@ +{ + "INI": "ArzlerAlmForestNumba_com4FlowPyCfg.ini", + "TAGS": ["forest", "infra", "standardTest", "real", "com4FlowPy", "numba"], + "DESCRIPTION": "Arzler Alm Test Case\ncom4FlowPy run with forest friction extension\ntest compares numba engine results against python engine results", + "REFERENCE": "https://github.com/OpenNHM/AvaFrameData/tree/main/avaArzl", + "FILES": ["zdelta", "cellCounts", "fpTravelAngleMax", "travelLengthMax"], + "TOPOTYPE": "real", + "AVADIR": "data/avaArzlerAlm", + "BENCHMARKED_AVAFRAME_VERSION": "2.1" +} diff --git a/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_cellCounts.tif b/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_cellCounts.tif new file mode 100644 index 000000000..7b691d331 Binary files /dev/null and b/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_cellCounts.tif differ diff --git a/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_fpTravelAngleMax.tif b/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_fpTravelAngleMax.tif new file mode 100644 index 000000000..db6f3bc06 Binary files /dev/null and b/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_fpTravelAngleMax.tif differ diff --git a/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_travelLengthMax.tif b/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_travelLengthMax.tif new file mode 100644 index 000000000..4dc8ec5ff Binary files /dev/null and b/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_travelLengthMax.tif differ diff --git a/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_zdelta.tif b/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_zdelta.tif new file mode 100644 index 000000000..42564caaf Binary files /dev/null and b/benchmarks/com4_avaArzlerAlmForestNumba/com4_b44fd62ee8_20260812_164438_zdelta.tif differ diff --git a/benchmarks/com4_avaArzlerAlmNullNumba/ArzlerAlmNullNumba_com4FlowPyCfg.ini b/benchmarks/com4_avaArzlerAlmNullNumba/ArzlerAlmNullNumba_com4FlowPyCfg.ini new file mode 100644 index 000000000..bd7403ec5 --- /dev/null +++ b/benchmarks/com4_avaArzlerAlmNullNumba/ArzlerAlmNullNumba_com4FlowPyCfg.ini @@ -0,0 +1,318 @@ +### Config File - This file contains the main settings for the com4FlowPy run +## Set your parameters +# This file will be overridden by local_com4FlowPyCfg.ini if it exists +# So copy this file to local_com4FlowPyCfg.ini, adjust your variables there + +# Optional settings------------------------------- +[GENERAL] +#++++++++++++ Flow-Py Model Parameters +# alpha: +# Angle-of-reach Alpha [°] - defines max. longitudinal runout limit +# equivalent to a Coulomb-friction of tan(alpha) in a sliding block model +#--------------------- +# exp: +# Spreading Coefficient (influences lateral spreading) +#--------------------- +# flux_threshold: +# Flux threshold (influences lateral spreading) +#--------------------- +# max_z: +# Energy-Line-Height Limit (can be interpreted as velocity limit) +# max_v = sqrt(max_z*19.62) +# typical values: +# - Avalanche: ~270 m ... ~72 m/s +# - Rockfall: ~130 m ... ~50 m/s +# - Debris-flow: ~12 m ... ~15 m/s +#--------------------- + +alpha = 25 +exp = 8 +flux_threshold = 3.0e-4 +max_z = 8848 + +#++++++++++++ Use Infrastructure + +infra = False + +#++++++++++++ preview Mode +# if previewMode = True, not every release cell is processend independently +# if a releaseCell is already "hit"/"affected" by a prior calculated path +# then the processing of this release cell is skipped +# can be used for a faster preview of model results (e.g. for checking of input parameters) or +# to save calculation time (e.g. when calculating with infrastructure) +# NOTE: results will deviate from "normal" model run and all outputLayers relying on the +# summed/combined output of different paths will not provide sensible values!!! + +previewMode = False + + +#++++++++++++ Use a dynamic u_max Limit +# Requires an additional tif-file containing the uMax (in m/s) +# or zDeltaMax (m) Limits +# in every cell where a release cell is. +# the paths for each release cells are calculated with these uMax values +# (computed to z_delta, similar to the max_z) or zDelta values. +# In varUmaxParameter the parameter (uMax in m/s or zDeltaMax in m) +# provided is given. +variableUmaxLim = False +varUmaxParameter = uMax + +#++++++++++++ Use a dynamic alpha angle +# Requires an additional tif-file containing the alpha angles +# in every cell where a release cell is. +# The paths for each release cells are calculated with these alpha angles +variableAlpha = False + +#++++++++++++ Use a dynamic exponent +# Requires an additional tif-file containing the exponents +# in every cell where a release cell is. +# The paths for each release cells are calculated with these exponent +variableExponent = False + +#++++++++++++ Use Forest Information +# NOTE AH 20240408: Dummy settings for Forest Interaction +# The forest implementation should mimick/reproduce the +# one from the "foreste_detrainment" branch in the avaframe/FlowPy repo +# which is partly described in D'Amboise et al. (2021) +# +# Forest-Interaction is only used if 'forest' is set to 'True' +#++++++++++++ +forest = False + +#++++++++++++ +# Type of Forest-Interaction Model to use (only works in conjunction with forest = True!) +# valid choices: +# +# * 'forestFriction' (D'Amboise et al., 2021) - use added friction on forest Pixels +# forest-Raster (set in 'forestPath') is assumed to have values between +# 0 (no Forest) and 1 (full protective function). +# model behavior is governed by 'maxAddedFrictionFor', 'minAddedFrictionFor' +# and 'velThForFriction' +# +# * 'forestDetrainment' (D'Amboise et al., 2022?) - added friction and detrainment on forested Pixels +# forest-Raster (set in 'forestPath') is assumed to have values between +# 0 (no Forest) and 1 (full protective function). +# model behavior is governed by 'maxAddedFrictionFor', 'minAddedFrictionFor' +# and 'velThForFriction', 'maxDetrainmentFor', 'minDetrainmentFor', 'velThForDetrain' +# +# * 'forestFrictionLayer' - provide Layer with local alpha / local delta-Alpha on forested pixels / for all pixels in +# the domain +# forest-Raster (set in 'forestPath') is assumed to have values between 0° (no friction) and 90° (inf. friciton) +# forest-Raster is either interpreted as 'absolute' or 'difference' +# 'absolute' - absolute 'alpha' values are used on forested pixels (all pixels with values > 0) +# 'relative' - 'alpha' values are interpreted as additional friction on forested pixels (all pixels with +# values > 0) +#++++++++++++ + +# ['forestFriction', 'forestDetrainment', 'forestFrictionLayer'] +forestModule = forestFriction +# ForestInteraction gives an additional output layer for the amount of forested cells +# a path ran through +# (now implemented, only when forest = True) -> TODO: we dont require a forestModule to compute +# for ForestInteraction Layer, but we need to read in the Forest - Layer (requires forest = True) +forestInteraction = False + +#++++++++++++ Forest added Friction +# These are the parameters for the "add friction on forested cells"-approach +# described in D'Amboise et al. (2021) +# +# maxAddedFrictionFor [°]: maximum increase of Alpha angle/basal friction on forested cells +# minAddedFrictionFor [°]: minimum increase of Alpha angle/basal friction on forested cells +# velThForFriction [m/s]: velocity threshold for forest-friction effect +#++++++++++++ + +maxAddedFrictionFor = 10 +minAddedFrictionFor = 2 +velThForFriction = 30 + +#++++++++++++ Forest Detrainment +# These are the parameters for the "forestDetrainment"-approach - D'Amboise et al. (2022??) +# The idea is to remove 'virtual mass' (aka 'flux' in Flow-Py) from the modeled flow/process +# on forested pixels in dependency of FSI and local z_delta +# +# maxDetrainmentFor [°]: default set to '0' (i.e. no detrainment) - foreste_detrainment val: +# minDetrainmentFor [°]: default set to '0' (i.e. no detrainment) - foreste_detrainment val: +# velThForDetrain [m/s]: velocity threshold for forest-detraiment effect # default set to '0' (i.e. no detrainment) - foreste_detrainment val: +#++++++++++++ + +maxDetrainmentFor = 0 +minDetrainmentFor = 0 +velThForDetrain = 0 + +#++++++++++++ ForestFrictionLayer +# The forest layer defines the alpha angle in forested areas. +# If forestFrictionLayerType = relative, forest layer values are added to the global alpha angle. +# e.g. if alpha = 24° and value on forest cell is 5°, then effective local alpha of 24° + 5° = 29° will be used. +# ------------------- +# If forestFrictionLayerType = absolute, forest layer values are used as alpha angle, +# if they are larger than the global alpha angle. +# e.g. if global alpha = 24° and value provided in the forestFrictionLayer on a cell is 35°, +# then local alpha of 35° will be used. +#++++++++++++ + +# ['absolute', 'relative'] +forestFrictionLayerType = absolute + +# skip Forest Effect (added forest friction) for first x meters (calculated in 3D - XYZ) +# should allow an initial acceleration phase of processes starting in or directly above +# dense forests (these would in many cases stop otherwise) +# if e.g. skipForestDist = 40, no added forestFriction will be assumed until 40 m 3D-distance +# along the path from the startCell. +skipForestDist = 0 + +#++++++++++++ Method to calculate flux distribution +# We fixed a bug in flowClass.py, which affects the distribution of the remaining flux, +# if a cell receives flux smaller than the provided flux_threshold. +# +# The default now (post Jan. 2025) is a calculation with the fixed bug! +# +# For backward compatibility the old version (prior to Jan. 2025 - with minor bug) can +# be switched on by setting "fluxDistOldVersion = True". + +fluxDistOldVersion = False + +# compute engine for the per-cell BFS. The 'numba' engine requires the optional +# 'numba' dependency (JIT-compiled kernel); it reproduces the 'python' engine's +# results and is much faster. infra/back-calculation, previewMode and relId +# outputs automatically fall back to the python engine. +# python : default, Cell-based reference implementation +# numba : JIT-compiled kernel (double precision), bit-for-bit vs 'python' +engine = numba + +#++++++++++++ Calculate with generations +# If calcGenerations = True, a different order of cells in a path are calculated. +# The results can vary when computing with generations. +# Additionally, the generation (iteration step) can be derived, which is required +# to get thalweg information. +# You can choose how the thalweg is computed: with the center of energy, center of flux and/or +# center of zdelta. The format should be: thalwegCenterOf = ['zdelta','energy', 'flux'] +# The saved variables can be chosen in thalwegVariables. Possible variables are. +# ['col', 'row', 'x', 'y', 'flux', 'fluxSum' 'flowEnergy', 'altitude', 'travelLength', 'zDelta', 'gamma', 'flowEnergyArray', 'zDeltaArray', 'fluxArray'] +# the arrays contain the respective values of the path and need much memory +calcGeneration = False +calcThalweg = False +thalwegCenterOf = ['zdelta','energy', 'flux'] +thalwegVariables = ['x', 'y', 'z', 's', 'zDelta'] + +#++++++++++++ Parameters for Tiling +# tileSize: size of tiles in x and y direction in meters (if total size of) x +# or y of input DEM is larger than tileSize, then the input raster +# layers are tiles +# tileOverlap: overlap between single tiles in m (5 km is rather conservative) +#--------------------------------------- + +tileSize = 15000 +tileOverlap = 5000 + +#++++++++++++ Parameters for CPU usage/multiprocessing +# Recommended to leave at default values unless performance has to be +# tweaked/optimized (e.g. for application on dedicated machines, maxing out +# computational resources) +# +# procPerCPUCore: number of processer per core for multiprocessing +# default value = 1 (it is recommended not to change this value +# unless you want to max. out CPUs, which might come at the cost +# of RAM comsumption!!!) +# chunkSize: tbd. +# maxChunks: Maximum number of single tasks for multiprocessing, if set too high +# this can max out RAM and lead to unexpected behavior (in this case +# probably an infinite loop --> you see it in the .log if this happens) +# On machines with a lot of CPUs and RAM this value might be set higher ... +#++++++++++++ +procPerCPUCore = 1 +chunkSize = 50 +maxChunks = 500 + +# Optional Custom Paths +[PATHS] + +# define format of output raster files (default = .tif) +# available options: [.tif, .asc] +# if you plan to utilize avaFrame tools to analyse result rasters +# then you should choose '.asc' here! +outputFileFormat = .tif + +# use LZW compression when writing tif raster files +useCompression = True + +# define noData value that is assigned in output rasters in cells that are not affected by the process (default: -9999) +# when changing this value BE AWARE that noData values can have same values as a affected cells +# (e.g., when using outputNoDataValue = 0) +outputNoDataValue = -9999 + +# define the different output files that are written to disk +# default = 'zDelta|cellCounts|travelLengthMax|fpTravelAngleMax' +# additional options: +# slTravelAngle +# flux +# zDeltaSum +# routFluxSum +# depFluxSum +# travelLengthMin +# fpTravelAngleMin +# if forestInteraction: forestInteraction is automatically added to outputs +# if infra: backCalculation is automatically added to output +outputFiles = zDelta|travelLengthMax|fpTravelAngleMax|flux|cellCounts + +# whether simulation results with the same simHash of the running simulation already exists +# AND the existing resultsFolder already contains valid com4FlowPy outputs. +# 1) overwriteResults = default ... does not re-run a simulation if results folder (res_) and .json +# -- if there are remnants from a previously attempted but not succesfully finished simulation (e.g. existing res_ and/or .json) +# --> delete existing results folder (res_) and .json and run Simulation +# 2) overwriteResults = reRunAndOverwrite ... deletes existing results folder (res_) and .json and runs Simulation +# 3) overwriteResults = reRunAndBackup ... moves existing result folder (res_) and .json to a dedicated backup folder (e.g. in /BACKUP) + +overwriteResults = reRunAndOverwrite + +#++++++++++++ Custom paths True/False +# default: False +# if set to 'False': +# the default AvaFrame Path/Directory structure is used +# i.e. the Work- and Result-Folders and inputFiles are taken from cfgMain and +# Input-Files are expected in AvaFrame directory structure +# if set to 'True': +# workDir, demPath, releasePath and [infraPath <-> if 'infra'==True] +# can be set here +# +#++++++++++++ +useCustomPaths = False +# for now only works with 'useCustomPaths=True'; if 'True' the temp folder is +# deleted after the output-files have been successfully written to disk +deleteTempFolder = False + +# if useCustomPathDEM set to True: only DEM path can set here, the other input data are used from +# avalancheDir (AvaFrame folder structure) (works only if useCustomPaths = False) +useCustomPathDEM = False + +# if 'useCustomPaths = True', then these input and output paths are used +# the output Folder is automatically placed inside the 'workDir' +workDir = +demPath = +releasePath = +infraPath = +forestPath = +varUmaxPath = +varAlphaPath = +varExponentPath = + +# plot save results flags------------------------ +# NOTE-TODO: These Flags still don't do anything, +# would be cool to have full integration for +# com4FlowPy with avaframe post-processing utils ... +[FLAGS] +# Plot the avalanche path and DEM figure +plotPath = False + +# Plot the avalanche profile with alpha beta points +plotProfile = False + +# Save the profile figure +saveProfile = True + +# Write results to txt file +writeRes = True + +# keep intermediate results +fullOut = False + +#---------------------------------------------------------- diff --git a/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_cellCounts.tif b/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_cellCounts.tif new file mode 100644 index 000000000..3380c6270 Binary files /dev/null and b/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_cellCounts.tif differ diff --git a/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_flux.tif b/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_flux.tif new file mode 100644 index 000000000..7af2e8c94 Binary files /dev/null and b/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_flux.tif differ diff --git a/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_fpTravelAngleMax.tif b/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_fpTravelAngleMax.tif new file mode 100644 index 000000000..c2497d0b2 Binary files /dev/null and b/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_fpTravelAngleMax.tif differ diff --git a/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_travelLengthMax.tif b/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_travelLengthMax.tif new file mode 100644 index 000000000..8b9b4c35a Binary files /dev/null and b/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_travelLengthMax.tif differ diff --git a/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_zdelta.tif b/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_zdelta.tif new file mode 100644 index 000000000..62805cb79 Binary files /dev/null and b/benchmarks/com4_avaArzlerAlmNullNumba/com4_03f0497040_20260812_165012_zdelta.tif differ diff --git a/benchmarks/com4_avaArzlerAlmNullNumba/com4_avaArzlerAlmNullNumba_desDict.json b/benchmarks/com4_avaArzlerAlmNullNumba/com4_avaArzlerAlmNullNumba_desDict.json new file mode 100644 index 000000000..c7657323e --- /dev/null +++ b/benchmarks/com4_avaArzlerAlmNullNumba/com4_avaArzlerAlmNullNumba_desDict.json @@ -0,0 +1,10 @@ +{ + "INI": "ArzlerAlmNullNumba_com4FlowPyCfg.ini", + "TAGS": ["null", "infra", "standardTest", "real", "com4FlowPy", "numba"], + "DESCRIPTION": "Arzler Alm Avalanche, real Topopgraphy\ndefault com4FlowPy .ini settings\ntest compares numba engine results against python engine results", + "REFERENCE": "https://github.com/OpenNHM/AvaFrameData/tree/main/avaArzl", + "FILES": ["zdelta", "cellCounts", "fpTravelAngleMax", "travelLengthMax", "flux"], + "TOPOTYPE": "real", + "AVADIR": "data/avaArzlerAlm", + "BENCHMARKED_AVAFRAME_VERSION": "2.1" +} diff --git a/benchmarks/com4_avaParabolaChannelPaperNullNumba/avaParabolaChannelPaperNullNumba_com4FlowPyCfg.ini b/benchmarks/com4_avaParabolaChannelPaperNullNumba/avaParabolaChannelPaperNullNumba_com4FlowPyCfg.ini new file mode 100644 index 000000000..f516718a1 --- /dev/null +++ b/benchmarks/com4_avaParabolaChannelPaperNullNumba/avaParabolaChannelPaperNullNumba_com4FlowPyCfg.ini @@ -0,0 +1,318 @@ +### Config File - This file contains the main settings for the com4FlowPy run +## Set your parameters +# This file will be overridden by local_com4FlowPyCfg.ini if it exists +# So copy this file to local_com4FlowPyCfg.ini, adjust your variables there + +# Optional settings------------------------------- +[GENERAL] +#++++++++++++ Flow-Py Model Parameters +# alpha: +# Angle-of-reach Alpha [°] - defines max. longitudinal runout limit +# equivalent to a Coulomb-friction of tan(alpha) in a sliding block model +#--------------------- +# exp: +# Spreading Coefficient (influences lateral spreading) +#--------------------- +# flux_threshold: +# Flux threshold (influences lateral spreading) +#--------------------- +# max_z: +# Energy-Line-Height Limit (can be interpreted as velocity limit) +# max_v = sqrt(max_z*19.62) +# typical values: +# - Avalanche: ~270 m ... ~72 m/s +# - Rockfall: ~130 m ... ~50 m/s +# - Debris-flow: ~12 m ... ~15 m/s +#--------------------- + +alpha = 25 +exp = 8 +flux_threshold = 3.0e-4 +max_z = 8848 + +#++++++++++++ Use Infrastructure + +infra = False + +#++++++++++++ preview Mode +# if previewMode = True, not every release cell is processend independently +# if a releaseCell is already "hit"/"affected" by a prior calculated path +# then the processing of this release cell is skipped +# can be used for a faster preview of model results (e.g. for checking of input parameters) or +# to save calculation time (e.g. when calculating with infrastructure) +# NOTE: results will deviate from "normal" model run and all outputLayers relying on the +# summed/combined output of different paths will not provide sensible values!!! + +previewMode = False + + +#++++++++++++ Use a dynamic u_max Limit +# Requires an additional tif-file containing the uMax (in m/s) +# or zDeltaMax (m) Limits +# in every cell where a release cell is. +# the paths for each release cells are calculated with these uMax values +# (computed to z_delta, similar to the max_z) or zDelta values. +# In varUmaxParameter the parameter (uMax in m/s or zDeltaMax in m) +# provided is given. +variableUmaxLim = False +varUmaxParameter = uMax + +#++++++++++++ Use a dynamic alpha angle +# Requires an additional tif-file containing the alpha angles +# in every cell where a release cell is. +# The paths for each release cells are calculated with these alpha angles +variableAlpha = False + +#++++++++++++ Use a dynamic exponent +# Requires an additional tif-file containing the exponents +# in every cell where a release cell is. +# The paths for each release cells are calculated with these exponent +variableExponent = False + +#++++++++++++ Use Forest Information +# NOTE AH 20240408: Dummy settings for Forest Interaction +# The forest implementation should mimick/reproduce the +# one from the "foreste_detrainment" branch in the avaframe/FlowPy repo +# which is partly described in D'Amboise et al. (2021) +# +# Forest-Interaction is only used if 'forest' is set to 'True' +#++++++++++++ +forest = False + +#++++++++++++ +# Type of Forest-Interaction Model to use (only works in conjunction with forest = True!) +# valid choices: +# +# * 'forestFriction' (D'Amboise et al., 2021) - use added friction on forest Pixels +# forest-Raster (set in 'forestPath') is assumed to have values between +# 0 (no Forest) and 1 (full protective function). +# model behavior is governed by 'maxAddedFrictionFor', 'minAddedFrictionFor' +# and 'velThForFriction' +# +# * 'forestDetrainment' (D'Amboise et al., 2022?) - added friction and detrainment on forested Pixels +# forest-Raster (set in 'forestPath') is assumed to have values between +# 0 (no Forest) and 1 (full protective function). +# model behavior is governed by 'maxAddedFrictionFor', 'minAddedFrictionFor' +# and 'velThForFriction', 'maxDetrainmentFor', 'minDetrainmentFor', 'velThForDetrain' +# +# * 'forestFrictionLayer' - provide Layer with local alpha / local delta-Alpha on forested pixels / for all pixels in +# the domain +# forest-Raster (set in 'forestPath') is assumed to have values between 0° (no friction) and 90° (inf. friciton) +# forest-Raster is either interpreted as 'absolute' or 'difference' +# 'absolute' - absolute 'alpha' values are used on forested pixels (all pixels with values > 0) +# 'relative' - 'alpha' values are interpreted as additional friction on forested pixels (all pixels with +# values > 0) +#++++++++++++ + +# ['forestFriction', 'forestDetrainment', 'forestFrictionLayer'] +forestModule = forestFriction +# ForestInteraction gives an additional output layer for the amount of forested cells +# a path ran through +# (now implemented, only when forest = True) -> TODO: we dont require a forestModule to compute +# for ForestInteraction Layer, but we need to read in the Forest - Layer (requires forest = True) +forestInteraction = False + +#++++++++++++ Forest added Friction +# These are the parameters for the "add friction on forested cells"-approach +# described in D'Amboise et al. (2021) +# +# maxAddedFrictionFor [°]: maximum increase of Alpha angle/basal friction on forested cells +# minAddedFrictionFor [°]: minimum increase of Alpha angle/basal friction on forested cells +# velThForFriction [m/s]: velocity threshold for forest-friction effect +#++++++++++++ + +maxAddedFrictionFor = 10 +minAddedFrictionFor = 2 +velThForFriction = 30 + +#++++++++++++ Forest Detrainment +# These are the parameters for the "forestDetrainment"-approach - D'Amboise et al. (2022??) +# The idea is to remove 'virtual mass' (aka 'flux' in Flow-Py) from the modeled flow/process +# on forested pixels in dependency of FSI and local z_delta +# +# maxDetrainmentFor [°]: default set to '0' (i.e. no detrainment) - foreste_detrainment val: +# minDetrainmentFor [°]: default set to '0' (i.e. no detrainment) - foreste_detrainment val: +# velThForDetrain [m/s]: velocity threshold for forest-detraiment effect # default set to '0' (i.e. no detrainment) - foreste_detrainment val: +#++++++++++++ + +maxDetrainmentFor = 0 +minDetrainmentFor = 0 +velThForDetrain = 0 + +#++++++++++++ ForestFrictionLayer +# The forest layer defines the alpha angle in forested areas. +# If forestFrictionLayerType = relative, forest layer values are added to the global alpha angle. +# e.g. if alpha = 24° and value on forest cell is 5°, then effective local alpha of 24° + 5° = 29° will be used. +# ------------------- +# If forestFrictionLayerType = absolute, forest layer values are used as alpha angle, +# if they are larger than the global alpha angle. +# e.g. if global alpha = 24° and value provided in the forestFrictionLayer on a cell is 35°, +# then local alpha of 35° will be used. +#++++++++++++ + +# ['absolute', 'relative'] +forestFrictionLayerType = absolute + +# skip Forest Effect (added forest friction) for first x meters (calculated in 3D - XYZ) +# should allow an initial acceleration phase of processes starting in or directly above +# dense forests (these would in many cases stop otherwise) +# if e.g. skipForestDist = 40, no added forestFriction will be assumed until 40 m 3D-distance +# along the path from the startCell. +skipForestDist = 0 + +#++++++++++++ Method to calculate flux distribution +# We fixed a bug in flowClass.py, which affects the distribution of the remaining flux, +# if a cell receives flux smaller than the provided flux_threshold. +# +# The default now (post Jan. 2025) is a calculation with the fixed bug! +# +# For backward compatibility the old version (prior to Jan. 2025 - with minor bug) can +# be switched on by setting "fluxDistOldVersion = True". + +fluxDistOldVersion = False + +# compute engine for the per-cell BFS. The 'numba' engine requires the optional +# 'numba' dependency (JIT-compiled kernel); it reproduces the 'python' engine's +# results and is much faster. infra/back-calculation, previewMode and relId +# outputs automatically fall back to the python engine. +# python : default, Cell-based reference implementation +# numba : JIT-compiled kernel (double precision), bit-for-bit vs 'python' +engine = numba + +#++++++++++++ Calculate with generations +# If calcGenerations = True, a different order of cells in a path are calculated. +# The results can vary when computing with generations. +# Additionally, the generation (iteration step) can be derived, which is required +# to get thalweg information. +# You can choose how the thalweg is computed: with the center of energy, center of flux and/or +# center of zdelta. The format should be: thalwegCenterOf = ['zdelta','energy', 'flux'] +# The saved variables can be chosen in thalwegVariables. Possible variables are. +# ['col', 'row', 'x', 'y', 'flux', 'fluxSum' 'flowEnergy', 'altitude', 'travelLength', 'zDelta', 'gamma', 'flowEnergyArray', 'zDeltaArray', 'fluxArray'] +# the arrays contain the respective values of the path and need much memory +calcGeneration = False +calcThalweg = False +thalwegCenterOf = ['zdelta','energy', 'flux'] +thalwegVariables = ['x', 'y', 'z', 's', 'zDelta'] + +#++++++++++++ Parameters for Tiling +# tileSize: size of tiles in x and y direction in meters (if total size of) x +# or y of input DEM is larger than tileSize, then the input raster +# layers are tiles +# tileOverlap: overlap between single tiles in m (5 km is rather conservative) +#--------------------------------------- + +tileSize = 15000 +tileOverlap = 5000 + +#++++++++++++ Parameters for CPU usage/multiprocessing +# Recommended to leave at default values unless performance has to be +# tweaked/optimized (e.g. for application on dedicated machines, maxing out +# computational resources) +# +# procPerCPUCore: number of processer per core for multiprocessing +# default value = 1 (it is recommended not to change this value +# unless you want to max. out CPUs, which might come at the cost +# of RAM comsumption!!!) +# chunkSize: tbd. +# maxChunks: Maximum number of single tasks for multiprocessing, if set too high +# this can max out RAM and lead to unexpected behavior (in this case +# probably an infinite loop --> you see it in the .log if this happens) +# On machines with a lot of CPUs and RAM this value might be set higher ... +#++++++++++++ +procPerCPUCore = 1 +chunkSize = 50 +maxChunks = 500 + +# Optional Custom Paths +[PATHS] + +# define format of output raster files (default = .tif) +# available options: [.tif, .asc] +# if you plan to utilize avaFrame tools to analyse result rasters +# then you should choose '.asc' here! +outputFileFormat = .tif + +# use LZW compression when writing tif raster files +useCompression = False + +# define noData value that is assigned in output rasters in cells that are not affected by the process (default: -9999) +# when changing this value BE AWARE that noData values can have same values as a affected cells +# (e.g., when using outputNoDataValue = 0) +outputNoDataValue = -9999 + +# define the different output files that are written to disk +# default = 'zDelta|cellCounts|travelLengthMax|fpTravelAngleMax' +# additional options: +# slTravelAngle +# flux +# zDeltaSum +# routFluxSum +# depFluxSum +# travelLengthMin +# fpTravelAngleMin +# if forestInteraction: forestInteraction is automatically added to outputs +# if infra: backCalculation is automatically added to output +outputFiles = zDelta|travelLengthMax|fpTravelAngleMax|flux|cellCounts + +# whether simulation results with the same simHash of the running simulation already exists +# AND the existing resultsFolder already contains valid com4FlowPy outputs. +# 1) overwriteResults = default ... does not re-run a simulation if results folder (res_) and .json +# -- if there are remnants from a previously attempted but not succesfully finished simulation (e.g. existing res_ and/or .json) +# --> delete existing results folder (res_) and .json and run Simulation +# 2) overwriteResults = reRunAndOverwrite ... deletes existing results folder (res_) and .json and runs Simulation +# 3) overwriteResults = reRunAndBackup ... moves existing result folder (res_) and .json to a dedicated backup folder (e.g. in /BACKUP) + +overwriteResults = reRunAndOverwrite + +#++++++++++++ Custom paths True/False +# default: False +# if set to 'False': +# the default AvaFrame Path/Directory structure is used +# i.e. the Work- and Result-Folders and inputFiles are taken from cfgMain and +# Input-Files are expected in AvaFrame directory structure +# if set to 'True': +# workDir, demPath, releasePath and [infraPath <-> if 'infra'==True] +# can be set here +# +#++++++++++++ +useCustomPaths = False +# for now only works with 'useCustomPaths=True'; if 'True' the temp folder is +# deleted after the output-files have been successfully written to disk +deleteTempFolder = False + +# if useCustomPathDEM set to True: only DEM path can set here, the other input data are used from +# avalancheDir (AvaFrame folder structure) (works only if useCustomPaths = False) +useCustomPathDEM = False + +# if 'useCustomPaths = True', then these input and output paths are used +# the output Folder is automatically placed inside the 'workDir' +workDir = +demPath = +releasePath = +infraPath = +forestPath = +varUmaxPath = +varAlphaPath = +varExponentPath = + +# plot save results flags------------------------ +# NOTE-TODO: These Flags still don't do anything, +# would be cool to have full integration for +# com4FlowPy with avaframe post-processing utils ... +[FLAGS] +# Plot the avalanche path and DEM figure +plotPath = False + +# Plot the avalanche profile with alpha beta points +plotProfile = False + +# Save the profile figure +saveProfile = True + +# Write results to txt file +writeRes = True + +# keep intermediate results +fullOut = False + +#---------------------------------------------------------- diff --git a/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_avaParabolaChannelPaperNumba_desDict.json b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_avaParabolaChannelPaperNumba_desDict.json new file mode 100644 index 000000000..802f4ee89 --- /dev/null +++ b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_avaParabolaChannelPaperNumba_desDict.json @@ -0,0 +1,10 @@ +{ +"INI": "avaParabolaChannelPaperNullNumba_com4FlowPyCfg.ini", +"TAGS": ["null", "standardTest", "idealized", "com4FlowPy", "numba"], +"DESCRIPTION": "Idealized/generic Topography - Parabolic Slope with channel\nExample also used in D'Amboise et al. (2022)\ntest compares numba engine results against python engine results", +"REFERENCE": "https://doi.org/10.5194/gmd-15-2423-2022", +"FILES": ["zdelta", "cellCounts", "fpTravelAngleMax", "travelLengthMax"], +"TOPOTYPE": "idealized", +"AVADIR": "data/avaParabChannelPaperFP", +"BENCHMARKED_AVAFRAME_VERSION": "2.1" +} diff --git a/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_cellCounts.tif b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_cellCounts.tif new file mode 100644 index 000000000..3ebd41d99 Binary files /dev/null and b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_cellCounts.tif differ diff --git a/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_flux.tif b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_flux.tif new file mode 100644 index 000000000..efbe0189d Binary files /dev/null and b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_flux.tif differ diff --git a/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_fpTravelAngleMax.tif b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_fpTravelAngleMax.tif new file mode 100644 index 000000000..e121dc30e Binary files /dev/null and b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_fpTravelAngleMax.tif differ diff --git a/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_travelLengthMax.tif b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_travelLengthMax.tif new file mode 100644 index 000000000..2f95851fb Binary files /dev/null and b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_travelLengthMax.tif differ diff --git a/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_zdelta.tif b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_zdelta.tif new file mode 100644 index 000000000..7925d9519 Binary files /dev/null and b/benchmarks/com4_avaParabolaChannelPaperNullNumba/com4_ff6612bd19_20260103_135958_zdelta.tif differ diff --git a/docs/develop.rst b/docs/develop.rst index 1bfb1fa64..a6a151091 100644 --- a/docs/develop.rst +++ b/docs/develop.rst @@ -273,12 +273,17 @@ To run the com4FlowPy standard tests, move to ``AvaFrame/avaframe`` and run: :: pixi run python runStandardTestsCom4FlowPy.py +To also run com4FlowPy standard tests for the new *numba engine* run the script from within the numba environment with: :: + + pixi run -e numba python runStandardTestsCom4FlowPy.py + The markdown-style report of the comparison is saved as: ``tests/reports/standardTestsReportCom4FlowPy.md``. Adding a new benchmark test case is similar as for com1DFA, in the json file containing the benchmark infos, the ``TAGS`` need to contain -``"com4FlowPy"``. +``"com4FlowPy"``, if a new test is targeting the *numba engine* ``"numba"`` should be additionally added to ``TAGS``. +In any case existing benchmark tests can be consulted for guidance. How to add a friction model diff --git a/docs/moduleCom4FlowPy.rst b/docs/moduleCom4FlowPy.rst index 7efdc1ffb..02b494f69 100644 --- a/docs/moduleCom4FlowPy.rst +++ b/docs/moduleCom4FlowPy.rst @@ -30,7 +30,7 @@ The motivational background and concepts behind the model, as well as a list of :ref:`theoryCom4FlowPy:com4FlowPy theory`. -Running the code +Running the code [python engine] ---------------- Generate an environment as described in :ref:`developinstall:Script Installation (Linux)` or @@ -69,7 +69,19 @@ or setup model parameters/config, import and run ``runCom4FlowPy`` from an exter Setup of the ConfigParser object needs to reflect the structure of the ``avaframe/com4FlowPy/(local_)com4FlowPyCfg.ini``. In this context the configuration file can be overwritten, for more information see :ref:`complexUsage:Override configuration`. - +Running the code [numba engine] +---------------- + +Running :py:mod:`com4FlowPy` with the new, optional *numba engine* (thx to `jmasseysykes `_) instead of the standard *python engine* +results in significant speed-ups (up to 30-45 times confirmed), while producing bit-level identical results for most applications (see the `PR `_ for detailed discussion). + +In order to use the *numba engine* you have to select it in the ``avaframe/com4FlowPy/(local_)com4FlowPyCfg.ini`` by setting ``engine = numba`` and then run the model via command line in the numba environment:: + + pixi run -e numba python runCom4FlowPy.py + +.. Note:: + - currently not all implemented functionalities are also available in the *numba engine* - if you use the options ``infra = True``, ``previewMode=True`` or define ``relIDCount`` or ``relIdPolygon`` in the list of desired model outputs ``outputFiles`` in the cfg, the model will fall back to the python engine, even if ``engine = numba`` is specified. + - All other cases (e.g. ``forest = True`` with all available ``forestModule`` options and ``forestInteraction = True``) are supported at this point. Configuration ---------------- diff --git a/pyproject.toml b/pyproject.toml index 458d80a3f..9a97a70f4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -46,6 +46,12 @@ dependencies = [ "fiona", ] +[project.optional-dependencies] +# Optional JIT compute engine for com4FlowPy (config: engine = numba | numba32). +# If numba is not installed, com4FlowPy transparently falls back to the pure-Python +# (Cell-based) engine, so this dependency is not required for normal operation. +numba = ["numba"] + # Setuptools [tool.setuptools] include-package-data = true @@ -140,3 +146,4 @@ doc = ["doc", "dev"] prod = ["prod"] #rcs = ["rcs"] qgis = ["qgis", "dev"] +numba = ["dev", "numba"]