search

lautenberger / elmfire

Eulerian Level set Model of FIRE spread
https://elmfire.io
Eclipse Public License 2.0
22 stars 11 forks source link

Elliptical dimensions & slope #63

Open vvvvalvalval opened 6 days ago

vvvvalvalval commented 6 days ago

Hi @lautenberger, potential bug report(s) here. I'm looking at how the slope-parallel velocity is computed:

         DXDT_ROTATED = DYDT*C%NORMVECTORX_DMS + DXDT*C%NORMVECTORY_DMS !ft/min, parallel to slope
         C%UX = DXDT_ROTATED * C%UXOUSX * FTPMIN_TO_MPS !m/s, projected

         DYDT_ROTATED = DYDT*C%NORMVECTORY_DMS - DXDT*C%NORMVECTORX_DMS !ft/min, parallel to slope
         C%UY = DYDT_ROTATED * C%UYOUSY * FTPMIN_TO_MPS !m/s, projected

         C%VELOCITY = SQRT(DXDT_ROTATED*DXDT_ROTATED + DYDT_ROTATED*DYDT_ROTATED) ! ft/min, parallel to slope

These calculations (especially the last line) seem to rely on the assumption that the X and Y directions along the slope are perpendicular. That's not true in general - geometry doesn't work like that. If you walk on a slope and take one step North and then one step West, these 2 steps are generally not orthogonal vectors (1). Other parts of ELMFIRE's code might be affected by this (2).

Should this be considered a bug?

(1) Proof: let (nx, ny, nz) and (wx, wy, wz) the vector coordinates of the North and West steps respectively. By definition nx = 0 and wy = 0, so that the dot product between these vectors is nz*wz. In all slopes that are not grid-aligned, this is non-zero, which proves the 2 vectors are not orthogonal.

(2) I'm also concerned with the code that derives DX_DT and DY_DT. Even after reading the technical reference, it's not clear to me whether the elliptical dimensions are defined in the horizontal or slope-parallel plane; until that intention is clarified, I'll have difficulty telling what a correct implementation should do.

lautenberger commented 5 days ago

Val, this is deep in the weeds of how elliptical dimensions are applied. Is it possible to create an verification case similar to https://github.com/lautenberger/elmfire/tree/main/verification/01-elliptical-shape to test this out?

vvvvalvalval commented 5 days ago

Is it possible to create an verification case

@lautenberger happy to, but to do that I need you to answer the following question:

it's not clear to me whether the elliptical dimensions are defined in the horizontal or slope-parallel plane; until that intention is clarified, I'll have difficulty telling what a correct implementation should do.

In other words, when you picture the ellipses in your mind's eye, do you visualize them growing in the horizontal plane or the sloped plane?

vvvvalvalval commented 5 days ago

Empirical confirmation

Hi @lautenberger, I was able to confirm the issue empirically, through an experiment that rotates the landscape. What I did is that I started from the elliptical case study, added slope, and then rotated wind and aspect by 3 angles (0°, 28°, 90°).

This confirms an issue in the spread algorithm, because the 28° ellipse has a different-shape to the grid-aligned ellipses (0° and 90°), when all 3 ellipses should have exactly the same shape since we're only rotating the landscape:

Screenshot 2024-06-28 at 14 31 16

I also confirmed that slope is the cause of the issue, by running the same experiment with 0 slope:

Screenshot 2024-06-28 at 14 32 33

(The ellipses are oriented in reverse because the slope more than compensates for the wind).

See attached ZIP file if you want to run this for yourself.

vw-elmfire-sloppy-slopes.zip

vvvvalvalval commented 4 days ago

@lautenberger if you agree that this confirms the bug, I can work on finding a fix, but I'll need you to answer the above question 1st: https://github.com/lautenberger/elmfire/issues/63#issuecomment-2195126362

lautenberger commented 4 days ago

Thanks, Val - good stuff! Yes, if you could submit a PR send over some code snippets to address this issue that would be fantastic.

I think the answer to your question is I envision the ellipse growing in the sloped plane and then projected to the horizontal plane.

vvvvalvalval commented 4 days ago

@lautenberger thanks, working on it... Another question: is ws.tif meant to be the slope-tangential wind speed or the horizontal wind speed?

vvvvalvalval commented 22 hours ago

@lautenberger update: I did the math to make elliptical dimensions slope-aware. This will probably require heavier surgery on ELMFIRE than I thought, because the math essentially shows that expressing vectors of interest in x and y coordinates is doomed to fail, and it seems the ELMFIRE codebase really wants to do that. I don't have a PR yet, but in the meantime, see the attached document. elliptical-wavelets.pdf