Fire: opt-in WRF-style Rothermel wind coupling for the directional projection formula

[RaymndH]

Summary

The current directional projection formula evaluates Rothermel's wind
factor with the wind projected onto the front normal before raising it
to the B power, R(theta) ~ (U cos theta)^B. This is non-convex once
phi_w(B-1) > 1, so a finite ignition line's endpoints grow a Wulff wedge
whose tip falls well below the head rate over time - verified against a
Munoz-Esparza-style finite-line test case, where the head rate degraded
from ~100% to ~40% of the theoretical Rf over 2100s.

erf.fire.directional_wind_coupling = "wrf" (default: "projection",
unchanged) instead exponentiates Rothermel's wind/slope factor from the
raw, unprojected wind speed and scales the result by the cosine to the
spread direction afterward, matching WRF-Fire's fire_ros
(module_fr_fire_phys.F). This keeps R(theta) linear in cos(theta) --
the convex support function of a circle -- so the wedge does not form;
the same test case now tracks the theoretical head rate to within 0.5%
for the full run.

Default behavior is verified bit-for-bit identical to the prior code
(max|delta phi| = 0 over the whole domain on the same test case).

Also includes two small pieces of test infrastructure used to validate
this:

• erf.fire.prescribed_wind[_x/_y]: bypass atmospheric wind
  interpolation with an exact constant wind, for controlled
  idealized-wind test cases.
• erf.fire.use_wind_limit: allow disabling the MEWS wind-speed cap
  (default on, unchanged) -- needed when a test's wind speed exceeds the
  cap.

Test plan

• Rebuilt clean against current ERF-Hazard
• Munoz-Esparza-style finite-line test: directional_wind_coupling =
  "wrf" tracks theoretical Rf to within 0.5% for the full 2100s run
  (vs. degrading to ~40% of Rf with the default "projection" mode)
• Default ("projection") mode verified bit-for-bit identical to
  pre-change code

[hgopalan]

Thanks for this. I audited it in detail: test-merged onto current ERF-Hazard (clean, no conflicts), built Release with unit tests on (zero warnings, all 20 gtests in the DirectionalShape / WindLimit / LevelSetAdvection suites pass), and ran the existing Exec/RegTests/FireDirectionalShape decks with base and merged binaries.

Verdict: the WRF coupling itself is correct and does what the description says. But the PR also bundles a second, undocumented change to the front-normal stencil that alters every default directional run, so the "bit-for-bit identical" claim does not hold for the default scheme. Requesting changes to split that out.

Blocking

1. Hidden change to the front normal

fill_directional_ros now computes the normal as 0.5*(dm + dp) of one_sided_derivatives with the run's gradient scheme, instead of central differences. The default scheme is weno5z_front, so inside the front band the normal comes from two HJ-WENO5-Z reconstructions: a different stencil. This is not mentioned in the title, description or commit message.

Existing FireDirectionalShape projection deck, base vs merged binary, 1500 s:

deck	base head [m/s]	merged head [m/s]	arrival times
projection, default weno5z_front	0.18487 (40 % of the way from the Wulff tip to Rothermel)	0.14949 (7 %)	differ everywhere
projection, erf.fire.levelset.gradient = upwind	0.15279	0.15279	identical (phi differs by 1e-9, rounding of the averaged stencil)

So the bit-for-bit claim holds only on an upwind deck; the test that produced "max|delta phi| = 0" must have used levelset.gradient = upwind. The sharper normal pushes the projection closer to its exact Wulff tip, which is defensible in isolation, but it changes the README table of that regtest and every default directional run. It also leaves the code inconsistent: the hybrid path (ERF_HybridRos.H) calls fill_directional_ros positionally and so keeps a first-order normal, and Balbi keeps central differences.

Please revert the normal change in this PR and open it separately, with the re-measured FireDirectionalShape README table and the three places that still say "central differences" (file comment in ERF_DirectionalRos.H, Docs/sphinx_doc/theory/fire_propagation.rst, the checker docstring) updated.

2. Stale description

The description and commit message list erf.fire.use_wind_limit as new test infrastructure. That flag is already in ERF-Hazard from #381 and this diff contains no change to it.

Should fix

3. Validate the option combination at parse time

directional_wind_coupling = "wrf" is silently ignored with directional_shape = "ellipse" and with every model other than Rothermel. The surrounding parse block aborts on contradictory combinations (Anderson flanks without the ellipse), so this should abort too. BEHAVE (behave_ros_cell) uses the identical power-law form, so extending the coupling to it is a few lines and removes one of the silent cases.

4. No regression test

The finite-line case that motivated the change is not in the PR, and the convention here is opt-in default plus a deck that compares old and new. Cheapest route: add inputs_wrf (the projection deck plus erf.fire.directional_wind_coupling = "wrf") to Exec/RegTests/FireDirectionalShape and a wrf kind to check_firedirectionalshape.py whose exact normal speed is R0 (1 + phi_w(U) max(cos theta, 0)). I did that in a scratch copy with the merged binary:

wrf:
  head (+x)          PASS  0.24846 m/s from 181 cells vs Hopf 0.24945 m/s (-0.39 %)
  head vs Rothermel  PASS  0.24846 vs 0.24943 m/s (-0.39 %)
  back (-x)          PASS  0.02396 m/s vs Hopf 0.02405 m/s
  flank (+y)/(-y)    PASS  0.02403 m/s vs Hopf 0.02403 m/s

That is a one-deck proof of the claim on a point fire. The finite-line deck would be a welcome second test, since FireLineFire is periodic in y and so an infinite line.

5. Documentation and decks

No rows for the four new keys in Docs/sphinx_doc/Inputs.rst, nothing next to the directional_shape paragraph in the theory doc, and Exec/CanonicalTests/Fire/inputs_fire_master_reference does not list them. rothermel_ros_cell also needs @param entries for the two new cosine arguments.

6. Position the option against the existing fix

directional_shape = "ellipse" already removes the wedge and gives -0.67 % on the same deck. The file comment presents the ellipse as the only remedy; it should say what "wrf" adds (the model's own oblique rates, linear in cos theta, rather than an imposed ellipse, and the form WRF-Fire users will compare against).

Minor

• The comments call R0 (1 + phi_w cos theta) "the support function of a circle". With the clamp at zero for backing it is the support function of a stadium (a disc of radius R0 swept along the wind vector). Still convex, so the conclusion stands, but the wording is off.
• prescribed_wind bypasses the interpolation cleanly and fire_wind_ref has no ghost cells, so the setVal is fine. WAF and the FARSITE terrain corrections still apply downstream, so "exact constant wind" is only true with use_waf = false and use_terrain_wind = false; worth saying in the parameter comment.
• Checked against WRF-SFIRE fire_ros (module_fr_fire_phys.F): it caps the raw speed at 30 m/s, exponentiates, then multiplies by max(0, (v.n)/|v|) and max(0, (dz.n)/tanphi). The PR matches this apart from the 30 m/s cap, which ERF covers with the MEWS limit.

🤖 Generated with Claude Code