U vector not front-normal when it should?

[vvvvalvalval]

Summary: Potential bug report. AFAICT the calculation of the U vector is incorrect: it's not orthogonal to the fire front, and that causes systematic bias in computing fireline intensity (FLI). This is easily verified in the flanking-fire case.

---

Because we calculate FLI based on it, I understand that the U vector is supposed to be normal to the fire front:

Yet the U vector is computed as follows:

In particular, if ω=90° (flanking fire), then we find U*y = c. Except that for a flanking fire, we should have U*y = 0.

You might object: "it doesn't matter, because ultimately we compute the dot-product with ∇φ, so only the front-normal component of U has to be correct". That's true as far as the level-set PDE is concerned, but it very much matters for FLI calculations - with the above formula, the spread rate and consequently the FLI of a flanking fire are systematically overestimated (by a factor of (1 + (c/b)^2)^0.5).

AFAICT the code is faithful to the docs, and therefore makes this error too:

            DYDT     = (RDENOM * AACOSANG ) + 0.5 * (C%VELOCITY_DMS - C%VBACK)

[vvvvalvalval]

Btw, when I do the math, the formulas I end up with are:

C1 = 0.5 * (C%VELOCITY_DMS - C%VBACK)

DY_DT = C1*COSANG*COSANG + DENOM*COSANG
DX_DT = (C1*COSANG + DENOM)*SINANG

C%VELOCITY = DENOM + C1*COSANG

It is readily verified that this yields the correct values for VELOCITY in the heading, flanking and backing cases.

What's more, it's easy enough to see that another valid formula for the front-normal spread rate (VELOCITY) is DY_DT*COSANG + DX_DT*SINANG (Hint: take the dot product of the U vector to the front-normal unit vector, which has coordinates nx=sin(ω) and ny=cos(ω)).
Applying this formula to my U vector and yours yields the same result for VELOCITY (the one I wrote above). This shows that my proposed U vector and yours agree on the front-normal component, but yours also as a non-zero front-tangential component. This corroborates my hypothesis that the current U computation comes from a mathematical derivation which is suitable for solving the level-set PDE, but not suitable for computing the spread rate.

[lautenberger]

Val, same question as #63 - can we create a verification case with an exact solution to check this?

[vvvvalvalval]

@lautenberger I'll see what I can do - hopefully I'll have some empirical evidence for you in a few days.

In the meantime, as further theoretical evidence, I'll add that Richards (1995) agrees with the calculation of my previous comment:

(Note that I derived my calculation independently - it's probably not a coincidence that Richards and I agree!)

I now refine my hypothesis: could it be that ELMFIRE's current implementation arose from misinterpreting the meaning of the X and Y quantities defined equations (16) and (17) of Richards (1995)? They're not the same thing as a front-normal spread rate vector.

[vvvvalvalval]

@lautenberger here you go, here's an empirical case study.

Summary

1. I actually just reused the elliptical verification case, which contains enough information to reveal the issue.
2. I compared the reported spread rate (from outputs/vs_0000001_0025170.tif) to the empirical spread rate (derived from the Time-of-Arrival map, at outputs/time_of_arrival_0000001_0025170.tif)
3. As my theory predicted, the reported spread rate is correct for heading and backing fires, but overestimated for flanking fires. Same for Fire-Line Intensity.

Detailed repro

I ran the case study with:

root@job-big-nrepl-0-dft56:/elmfire/elmfire/verification/01-elliptical-shape# export PATH=/elmfire/elmfire/build/linux/bin:$PATH && ./01-elliptical-shape.sh

Then I downloaded the outputs to my laptop, and loaded them into NumPy arrays:

import matplotlib.pyplot as plt
import numpy
import rasterio

def fetch_raster_from_geotiff__rio(path_to_geotiff):
    with rasterio.open(path_to_geotiff) as src:
        ret = src.read(1).astype(float)
        ret[ret == src.nodata] = numpy.nan
        return ret

ell0_flin = fetch_raster_from_geotiff__rio('/Users/val/projects/SIG/elm-ellipse-0/outputs/flin_0000001_0025170.tif')
ell0_toa = fetch_raster_from_geotiff__rio('/Users/val/projects/SIG/elm-ellipse-0/outputs/time_of_arrival_0000001_0025170.tif')
ell0_vs = fetch_raster_from_geotiff__rio('/Users/val/projects/SIG/elm-ellipse-0/outputs/vs_0000001_0025170.tif')

For the record, here are some visualizations of the outputs:

I then computed the empirical spread rate using the ToA gradient:

def empirical_spread_rate(toa_arr2d, pixel_length=1.0):
    """
    Computes the front-normal spread rate from a Time-of-Arrival map.
    Accepts a 2D-array of ToA values (e.g. in min),
    a pixel length (e.g. in m),
    and returns a 2D-array of spread rates (e.g. in m/min).
    """
    gy, gx = numpy.gradient(toa_arr2d, pixel_length)
    # The empirical spread rate vector is |∇ToA|^-2 * ∇ToA;
    # the front-normal spread rate is its magnitude, i.e. |∇ToA|^-1:
    return (gx**2 + gy**2)**-0.5
    
ft_to_m = 0.3048
s_to_min = 60.0**-1
ell0_ros = empirical_spread_rate(ell0_toa, pixel_length=5) / (ft_to_m * s_to_min)

From there I displayed the empirical/reported ratio:

The above map should be 1 everywhere. We see flanking spread rate getting overestimated by a factor of 2.

And as expected, the reported Fire-Line Intensity has the same error patterns:

If the outputs were correct, the above map would be constant, since FLI is supposed to be proportional to RoS.