PATH OBSCURATION beam is truncated across meshes when XB endpoints are descending

[ManuelOsburg]

PATH OBSCURATION beam is truncated across meshes when XB endpoints are descending

FDS version: FDS-6.11.0-0-g369a20b-release

Summary

A PATH OBSCURATION (also TRANSMISSION) beam device is integrated across all meshes it crosses only if its XB endpoints are given in ascending order (x1≤x2, y1≤y2, z1≤z2). If the same physical beam is written descending in an axis it spans (first endpoint greater than second), it is silently truncated to the single mesh containing the beam midpoint — every other mesh along the path is dropped. The obscuration is a line integral and should be invariant to endpoint order; in a multi-mesh domain it is not. No warning is printed.

Minimal reproduction

Uniform, quiescent smoke field (no fire), so the optical depth is analytic:
b = K_m·rho_s·L, Obscuration = 100·(1−exp(−b)). Two meshes split at x=1.0.
K_m=8700, Y=3e-5 → K = K_m·rho_s ≈ 0.3119 m⁻¹.

issue_pathobsc.fds.txt

&HEAD CHID='issue_pathobsc'/
&MESH ID='A', IJK=10,5,5, XB=0.0,1.0,0.0,1.0,0.0,1.0/
&MESH ID='B', IJK=10,5,5, XB=1.0,2.0,0.0,1.0,0.0,1.0/
&TIME T_END=2.0/
&DUMP DT_DEVC=0.5/
&MISC TMPA=20.0/
&SPEC ID='MY SMOKE', MASS_EXTINCTION_COEFFICIENT=8700/
&INIT SPEC_ID(1)='MY SMOKE', MASS_FRACTION(1)=3.0E-5/
&PROP ID='smoke source', SPEC_ID='MY SMOKE'/
&DEVC XYZ=0.5,0.5,0.5, QUANTITY='EXTINCTION COEFFICIENT', SPEC_ID='MY SMOKE', ID='K_field'/
&DEVC XB=0.05,0.95,0.5,0.5,0.5,0.5, QUANTITY='PATH OBSCURATION', ID='beam_in_A',      PROP_ID='smoke source'/
&DEVC XB=0.05,1.95,0.5,0.5,0.5,0.5, QUANTITY='PATH OBSCURATION', ID='beam_ascending',  PROP_ID='smoke source'/
&DEVC XB=1.95,0.05,0.5,0.5,0.5,0.5, QUANTITY='PATH OBSCURATION', ID='beam_descending', PROP_ID='smoke source'/
&TAIL/

beam_ascending and beam_descending are the same physical line (endpoints swapped).

Expected vs. observed

device	expected	observed
beam_in_A (L=0.9 m, single mesh)	24.48 %	24.48 % ✓
beam_ascending (L=1.9 m, crosses A→B)	44.72 %	44.72 % ✓
beam_descending (same line, XB swapped)	44.72 %	25.65 % ✗ (only the 0.95 m in mesh B)

Single-mesh domain (un-split): both crossing beams give 44.72 %. Behaviour is the same in serial and under MPI (mpiexec -n 2), and is axis-general (a y- or z-split shows the same with descending y/z).

TRANSMISSION for the same beams reports 73.20 %/m (= exp(-K)·100) for both ascending and descending beams, i.e. it looks correct. That is because TRANSMISSION normalizes the (truncated) Σρ_s·Δx by the (equally truncated) integrated path length L, so for a uniform field the per-metre value is independent of how much length was integrated — it masks the truncation rather than avoiding it. (This matches the original report's observation, see below.)

Apparent cause in the source

(read.f90 at commit 369a20b.) A beam is integrated only over meshes that have a subdevice. Two places flag meshes:

• the midpoint mesh loop MESH_LOOP (read.f90:13977+) flags the mesh containing the beam midpoint (XYZ, set to the XB midpoint at 13873);
• CHECK_MESH_LOOP (read.f90:13900–13935) flags the other crossed meshes using an AABB overlap test, e.g. for x (13906):

  IF (XB(1)/=XB(2) .AND. (MIN(XB(2),M%XF)-MAX(XB(1),M%XS)<TWENTY_EPSILON_EB)) OVERLAPPING_X = .FALSE.

  which assumes XB(1) ≤ XB(2).

For beam quantities, CALL CHECK_XB(XB) (which sorts XB ascending) is skipped (read.f90:13898), so a user-supplied descending XB is preserved. The overlap test then returns "no overlap" for every mesh, so CHECK_MESH_LOOP adds no subdevices and only the midpoint mesh is integrated → truncation.

Suggested fix

Make the overlap test independent of endpoint order (the beam direction used by the path walk is stored separately in DV%X1..X2), e.g. use MIN(XB(1),XB(2)) / MAX(XB(1),XB(2)) in the three overlap tests in CHECK_MESH_LOOP — equivalently, restore the order-independent form that the 2022 fix introduced (see below).

Relationship to #10225 (appears to be a regression)

This is the same problem reported in #10225 (Jan 2022, FDS 6.7.7): PATH OBSCURATION across multiple meshes wrong "when the x or y position of the endpoint of a line DEVC is lower than the start point", while TRANSMISSION for the same coordinates looked correct. #10225 was fixed and closed in 2022 by commit 152e1592 ("Issue #10225. Rework previous fix"), which kept CHECK_XB skipped for beams but made the overlap test order-independent:

! 152e1592 (2022 fix): symmetric in XB(1)/XB(2) -> descending beams handled
IF (XB(1)/=XB(2) .AND. ((XB(1)>=M%XF.AND.XB(2)>=M%XF) .OR. (XB(1)<=M%XS.AND.XB(2)<=M%XS))) OVERLAPPING_X = .FALSE.

That order-independent form was later replaced (commit 4f4aca5c, 2024-08-28, "FDS Source: Fix < vs <= problem or SPATIAL_STATs"; first released in FDS-6.10.0) by a form that again assumes XB(1) ≤ XB(2):

! since 4f4aca5c (6.10.0+): assumes XB(1) <= XB(2) -> descending beams overlap no mesh
IF (XB(1)/=XB(2) .AND. (MIN(XB(2),M%XF)-MAX(XB(1),M%XS)<TWENTY_EPSILON_EB)) OVERLAPPING_X = .FALSE.

so the #10225 behaviour has returned. (Confirmed by tag bisection: the order-independent form is present in 6.7.8–6.9.1; the MIN/MAX form appears in 6.10.0, 6.10.1 and 6.11.0.) The SPATIAL_STATS fix in 4f4aca5c appears to have inadvertently dropped the order-independence added for #10225.

Note

The verification case Verification/Detectors/beam_detector.fds exercises only ascending beams, so this is not currently caught by the test suite.

[drjfloyd]

Is the solution to change to:

IF (XB(1)/=XB(2) .AND. (MIN(MAX(XB(1),XB(2)),M%XF)-MAX(MIN(XB(1),XB(2)),M%XS)<TWENTY_EPSILON_EB)) OVERLAPPING_X = .FALSE.

[rmcdermo]

I replaced this (starting at current 13897) with your line:

      IF (XB(1)>-1.E5_EB) THEN
         IF (QUANTITY/='PATH OBSCURATION' .AND. QUANTITY/='TRANSMISSION' .AND. QUANTITY/='FIRE DEPTH') CALL CHECK_XB(XB)
         BAD = .TRUE.
         CHECK_MESH_LOOP: DO NM=1,NMESHES
            M=>MESHES(NM)
            OVERLAPPING_X = .TRUE.
            OVERLAPPING_Y = .TRUE.
            OVERLAPPING_Z = .TRUE.
            IF (XB(1)/=XB(2) .AND. (MIN(MAX(XB(1),XB(2)),M%XF)-MAX(MIN(XB(1),XB(2)),M%XS)<TWENTY_EPSILON_EB)) &
               OVERLAPPING_X = .FALSE.
            ! IF (XB(1)==XB(2) .AND. (XB(1)> M%XF .OR. XB(2)< M%XS)) OVERLAPPING_X = .FALSE.

but still get the same result.

[ManuelOsburg]

Thanks @drjfloyd and @rmcdermo for the quick look.

Up front: the source-level analysis and build verification below were done with AI assistance — I'm not a Fortran expert, so please sanity-check the code specifics. The numbers, though, are from actual FDS runs against an analytic field.

I built FDS from source (gfortran 13.3 + OpenMPI, commit ge29dd78, i.e. the 6.11.0 line; HYPRE/SUNDIALS disabled) to test this end-to-end, and I can confirm:

@drjfloyd's order-independent overlap test fully fixes it. The unpatched build reproduces the bug (descending beam = 25.65 %); with the fix, every descending beam matches its ascending twin and the analytic value, on both axes and for multi-mesh paths:

case	unpatched	with fix	analytic
2-mesh, X descending	25.65 %	44.72 %	44.72 %
4-mesh, X descending	26.80 %	70.38 %	70.38 %
4-mesh, X descending starting in mesh 3	–	59.53 %	59.53 %
Y descending	(truncated)	70.38 %	70.38 %

Ascending/control beams are unchanged. The cross-mesh combine already does the right thing — EXCHANGE_DEVICE in main.f90 sums the per-mesh subdevice VALUE_1 (MPI_SUM) — so once a subdevice is created for each crossed mesh, fixing the overlap test in CHECK_MESH_LOOP is sufficient.

@rmcdermo — re "still get the same result": I think the rebuild may not have picked up read.f90. The MIN(MAX(...)) line is longer than 132 columns, and with the GNU flags (-std=f2018) that triggers -Werror=line-truncation, so read.o fails to compile and the previous binary in the build dir is silently reused. I hit exactly that — as a single long line the build errors out; splitting it with a continuation compiles and the fix works:

IF (XB(1)/=XB(2) .AND. &
    (MIN(MAX(XB(1),XB(2)),M%XF)-MAX(MIN(XB(1),XB(2)),M%XS)<TWENTY_EPSILON_EB)) OVERLAPPING_X = .FALSE.
IF (XB(3)/=XB(4) .AND. &
    (MIN(MAX(XB(3),XB(4)),M%YF)-MAX(MIN(XB(3),XB(4)),M%YS)<TWENTY_EPSILON_EB)) OVERLAPPING_Y = .FALSE.
IF (XB(5)/=XB(6) .AND. &
    (MIN(MAX(XB(5),XB(6)),M%ZF)-MAX(MIN(XB(5),XB(6)),M%ZS)<TWENTY_EPSILON_EB)) OVERLAPPING_Z = .FALSE.

A clean rebuild (or rm read.o fds_*) and a check that read.o actually recompiles should do it.

[drjfloyd]

should be fixed now.

[rmcdermo]

Yes, works for me. I had been compiling the linux version but running the macOS version (doh!).