Floating point error leading to extra particle insertion

Attached is an example where many INITs are used to insert particles, with the aim of each INIT covering exactly one cell. This is relevant when building complex arrangements of particles (such as with vegetation using a BULK_DENSITY_FILE) and we want to have INIT regions align exactly with grid lines to minimize cost.

In the example, the region of particles should cross each mesh, giving exactly 32 particles per mesh. However, we can see that the third mesh has extra particles:

 Run Time Diagnostics

       Time Step        1   March 22, 2024  15:16:04
       Step Size:    0.100E-01 s, Total Time:      0.010 s
       Pressure Iterations: 1
       Maximum Velocity Error:  0.24E-04 on Mesh 2 at (32,3,13)
       Maximum Pressure Error:  0.31E-06 on Mesh 3 at (3,26,27)
       ---------------------------------------------------------------
       Mesh    1
       Max CFL number:  0.52E-02 at (16,6,10)
       Max divergence:  0.19E-05 at (19,3,2)
       Min divergence: -0.20E-05 at (21,26,31)
       Max VN number:   0.60E-02 at (9,1,23)
       No. of Lagrangian Particles:  32
       Mesh    2
       Max CFL number:  0.50E-02 at (17,8,3)
       Max divergence:  0.19E-05 at (25,21,2)
       Min divergence: -0.19E-05 at (9,6,31)
       Max VN number:   0.55E-02 at (30,32,8)
       No. of Lagrangian Particles:  32
       Mesh    3
       Max CFL number:  0.51E-02 at (5,28,30)
       Max divergence:  0.20E-05 at (13,9,2)
       Min divergence: -0.21E-05 at (3,11,31)
       Max VN number:   0.50E-02 at (11,32,17)
       No. of Lagrangian Particles:  34

Marcos and I tracked this down to a floating point error in the M%CELLSI array, which is used in GET_IJK to determine which cells to include in each INIT based on XB:

https://github.com/firemodels/fds/blob/dd59a937836555b9eef5dc75351ccc8e0c6f7ae1/Source/part.f90#L1227-L1228

In the example case, I1 and I2 above should always be the same (one cell per INIT), but theres a couple of cases where it catches two cells and so adds two particles. The logic for mass conservation works out, so it's not catastrophic, but in large vegetation cases with potentially millions of INIT regions implicitly defined in a BULK_DENSITY_FILE, these extra particles can build up unwanted cost.

A potential solution is something like this in read.f90, where we ensure there is no floating point error in the CELLSI/J/K arrays for the locations that are exactly at grid lines. In other words, we know the value should be an integer so we enforce it:

   DO I=M%CELLSI_LO,M%CELLSI_HI
      M%CELLSI(I) = GINV(REAL(I,EB)/M%RDXINT,1,NM)*M%RDXI
      IF(MOD(I,500)==0) M%CELLSI(I) = REAL(NINT(M%CELLSI(I)),EB)
      M%CELLSI(I) = MAX(M%CELLSI(I),-0.9_EB)
      M%CELLSI(I) = MIN(M%CELLSI(I),REAL(M%IBAR)+0.9_EB)
   ENDDO

Unless there are better ideas, I'll run firebot to see if this tweak causes any issues voxel_offset_error.txt

firemodels / fds

Floating point error leading to extra particle insertion #12673