Extracting data on multiple 2d plane from PL3D data

[Sagar-Singhal]

I have PL3D data of fire simulation (which includes several number of meshes around 300). I wish to export data on different 2d planes for further processing and visualisation. Previously, I have used fdsreader for exporting data from 2D slice planes. I tried using this Python script :

import os
import fdsreader
import matplotlib.pyplot as plt

# Load FDS response
path_to_sim= os.path.join("path", "to", "simulation", "results")
sim_data = fdsreader.Simulation(path_to_sim)

# Check which Plot3D quantities are available
quantities_avail = sim_data.data_3d.quantities
print(quantities_avail )

# Pick a desired quantity
quantity = sim_data.data_3d.quantities[0]

# Read Plot3D data by desired quantity
plot_3d = sim_data.data_3d.get_by_quantity(quantity)

# Get coordinates and merge data from multiple meshes
data, coordinates = plot_3d.to_global(masked=True, return_coordinates=True, fill=np.nan)

# Create convenience function
def get_2d_slice(array, axis, slice_index, timestep):
    if axis in 'xyz':
        axis_dict = {'x': 0, 'y': 1, 'z': 2}
        axis = axis_dict[axis]
    return np.take(array[timestep], indices=slice_index, axis=axis)

# Get last time step, as an example
t = -1

# Define slice through Plot3D data along the x axis, i.e. PBX
# Use 'slice_index' to select the slice, i.e. n-th cell along the axis
slice_2d = get_2d_slice(data, axis='x', slice_index=5, timestep=t)

# Plot result
plt.imshow(slice_2d.T,
           origin='lower',
           cmap='viridis')

---

I am running this simulation.

---

&HEAD CHID='wind_PF_100mm_FP_6_meshes_TI_0_1', TITLE='Wind Blown 100 mm n-heptane pool fire' / 
! Multiple meshes
! Parallel processing
!&MESH IJK=325,75,25, XB=-0.3,1,-0.15,0.15,0.2,0.3 /
&MESH IJK=600,80,8, XB=-0.3,1.2,-0.1,0.1,0.0,0.02  /
&MESH IJK=600,80,8, XB=-0.3,1.2,-0.1,0.1,0.02,0.04  /
&MESH IJK=600,80,8, XB=-0.3,1.2,-0.1,0.1,0.04,0.06  /
&MESH IJK=150,20,6, XB=-0.3,1.2,-0.1,0.1,0.06,0.12  /
&MESH IJK=150,20,6, XB=-0.3,1.2,-0.1,0.1,0.12,0.18  /
&MESH IJK=150,20,6, XB=-0.3,1.2,-0.1,0.1,0.18,0.24  /
!&MESH IJK=650,50,25, XB=-0.3,1,-0.05,0.05,0.05,0.1 /
!&MESH IJK=650,50,50, XB=-0.3,1,-0.05,0.05,0.1,0.2 /
!&MESH IJK=325,25,50, XB=-0.3,1,-0.15,0.05,0.0,0.2 /
!&MESH IJK=325,25,50, XB=-0.3,1,0.05,0.15,0.0,0.2 /
&TIME T_END=5/
&DUMP DT_RESTART=5, DT_SLCF=0.01, DT_PL3D = 0.01, PLOT3D_QUANTITY(1:5)='TEMPERATURE','U-VELOCITY','V-VELOCITY','W-VELOCITY','HRRPUV', WRITE_XYZ=.TRUE. /
&MISC RESTART=F, SIMULATION_MODE='LES',IBLANK_SMV=.FALSE./
&REAC FUEL='N-HEPTANE',	SOOT_YIELD=0.015/
&SURF ID='FUEL'
      COLOR='RED'
      HRRPUA=1275.0/
&SURF ID='WIND'
      COLOR='BLUE'
      VEL=-2.5/
&VENT MB ='XMIN',SURF_ID='WIND',VEL_RMS=.2/
!&VENT PBX=-.30, SURF_ID='WIND'/
&VENT PBX= 1.2, SURF_ID='OPEN'/
&VENT PBY=-.1, SURF_ID='OPEN'/
&VENT PBY= .1, SURF_ID='OPEN'/
&VENT PBZ= .24, SURF_ID='OPEN'/
&MULT ID='cube array', DX=0.005,DY=0.005,DZ=0.005,
      I_LOWER=-50,I_UPPER=50,
      J_LOWER=-50,J_UPPER=50,
      K_UPPER=50/
! add synthetic turbulence to break symmetry
&VENT XB=-.05,.05,-.05,.05,0,0, XYZ=0,0,0.0, RADIUS=0.05, SPREAD_RATE=10000,IOR=3, SURF_ID='FUEL'/
!&VENT PBZ=0, SURF_ID='INERT'/
&SLCF PBZ=0, QUANTITY='VELOCITY', VECTOR=.TRUE./
&SLCF PBY=0.001, QUANTITY='VELOCITY', VECTOR=.TRUE./
&SLCF PBY=0.001, QUANTITY='MASS FRACTION', SPEC_ID='N-HEPTANE', CELL_CENTERED=.TRUE./
&SLCF PBY=0.001, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE./
&SLCF PBY=0.001, QUANTITY='HRRPUV', CELL_CENTERED=.TRUE./
&SLCF PBY=0.001, QUANTITY='DIVERGENCE', CELL_CENTERED=.TRUE./
&SLCF PBY=0.001, QUANTITY='MIXTURE FRACTION', CELL_CENTERED=.TRUE./
&SLCF PBY=0.001, QUANTITY='HRRPUV', CELL_CENTERED=.TRUE./
&SLCF PBZ=0.1, QUANTITY='HRRPUV', CELL_CENTERED=.TRUE./
&SLCF PBZ=0.1, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE./
&SLCF PBZ=0.1, QUANTITY='VELOCITY', VECTOR=.TRUE./
&SLCF PBX=0.000, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE./
&SLCF PBX=0.000,  QUANTITY='VELOCITY', VECTOR=.TRUE./
&SLCF PBX=-0.1, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE./
&SLCF PBX=-0.1,  QUANTITY='VELOCITY', VECTOR=.TRUE./
&SLCF PBX=0.1, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE./
&SLCF PBX=0.1,  QUANTITY='VELOCITY', VECTOR=.TRUE./
&SLCF PBX=0.2, QUANTITY='TEMPERATURE', CELL_CENTERED=.TRUE./
&SLCF PBX=0.2,  QUANTITY='VELOCITY', VECTOR=.TRUE./
&TAIL /

---

But there seems to be some issue with merging data from multiple meshes.

---

ValueError                                Traceback (most recent call last)
Cell In[9], line 22
     19 plot_3d = sim_data.data_3d.get_by_quantity(quantity)
     21 # Get coordinates and merge data from multiple meshes
---> 22 data, coordinates = plot_3d.to_global(masked=True, return_coordinates=True, fill=np.nan)
     25 # Create convenience function
     26 def get_2d_slice(array, axis, slice_index, timestep):

File c:\Users\Admin\anaconda3\Lib\site-packages\fdsreader\pl3d\pl3d.py:234, in Plot3D.to_global(self, masked, fill, return_coordinates)
    231     # If the slice should be masked, we set all cells at which an obstruction is in the
    232     # simulation space to the fill value set by the user
    233     if masked:
--> 234         subplot_data = np.where(mask, subplot_data, fill)
    236     grid[:, start_idx['x']: end_idx['x'], start_idx['y']: end_idx['y'],
    237     start_idx['z']: end_idx['z']] = subplot_data.reshape(
    238         (self.n_t, end_idx['x'] - start_idx['x'], end_idx['y'] - start_idx['y'],
    239          end_idx['z'] - start_idx['z']))
    241 if return_coordinates:

ValueError: operands could not be broadcast together with shapes (9,601,81,8) (11,601,81,8) () 

---

Edits

• Tristan: fixed Python and FORTRAN code formatting
• Tristan: fixed error message formatting

[TristanHehnen]

We'll have a look at it.

[Haukiy]

I had a similar problem. It was due to different cell sizes. If I'm not mistaken, we've got two cell sizes here: 0.0025 mm and 0.01 mm. So maybe the issue is related to that.

[TristanHehnen]

That could be.
@Sagar-Singhal could you please load the p3d data of a single mesh and see if this works?

[TristanHehnen]

But then I wonder, if you load the SLCF which are crossing different resolution do you see the same error?

[Haukiy]

@TristanHehnen For SLCF, I get a tuple (related to issue #49). For p3d, I get array shape problems similar to that above. With uniform cell size for all meshes, no tuple and no "shaping problems".

[TristanHehnen]

@Haukiy Okay, in the tuple there are then essentially the globals sorted by cell size?

Edit: yes that seem to be the case

[Sagar-Singhal]

@TristanHehnen
I agree with @Haukiy. I also do see tuple while working with SLCF and similar shape problems when working with p3d. However, when working with SLCF, boolean operation is sufficient to get combined data.

I tried a case with uniform cell size, however number of cells were around 80 millions, so I get an MemoryError. Unable to allocate 524. GiB for an array with shape (910, 401, 401, 481) and data type float64. Is it possible to change data type to float32?

[TristanHehnen]

Okay, this tuple-thing for SLCF makes sense to me, because there may be slice configurations where the sections across different resolutions do not align. Then the user has to combine them on their own, because we cannot know what is the desired outcome in such a case.

It may make sense to have such tuples also with the p3d.

@Sagar-Singhal Just to be clear, you see the memory error raised from Python and not from FDS?

[JanVogelsang]

Please check out the to_global function here: https://firedynamics.github.io/fdsreader/slcf.html
The way FDS handles slices is to produce one slice array for each mesh and output those to individual files. The fdsreader collects all these sub-slices and stores them in a single Slice object. As you have pointed out already, the meshes might have different cell sizes, it is therefore not trivial to decide how to further process these sub-slices. The to_global method creates a single array out of all sub-slices, where sub-slices with higher cell sizes are just repeated based on the smallest cell size of any of the meshes the slice cuts through. This is great for visualization, but for any sorts of analysis I would suggest to just iterate over all sub-slices (e.g. using the subslices property) and handle them individually.

PL3D objects are handled in the same way, just with one additional dimension. They are basically handled equivalently to 3D slices in FDS (at least from what is output, not sure about the internals).

Edit: Nevermind, I didn't see you were actually using the to_global helper function already. As there are so many edge cases in FDS, it seems like you actually found a new one which wasn't handled by our reader yet. If SLCF works and PL3D doesn't, it might be the case that this edge case was already fixed for SLCF but wasn't added to PL3D files. So the best way to debug this issue would be to diff the SLCF and PL3D to_global implementations and find any differences. They should be nearly identical.

[JanVogelsang]

@TristanHehnen I agree with @Haukiy. I also do see tuple while working with SLCF and similar shape problems when working with p3d. However, when working with SLCF, boolean operation is sufficient to get combined data.

I tried a case with uniform cell size, however number of cells were around 80 millions, so I get an MemoryError. Unable to allocate 524. GiB for an array with shape (910, 401, 401, 481) and data type float64. Is it possible to change data type to float32?

It's currently not possible to tell the fdsreader which datatype to use internally. This limitation is because the fdsreader just uses the same format as FDS uses as this allows to just directly read the data arrays from disk to memory. However, FDS uses 32-bit floating points as far as I am aware, so I am quite surprised to see a np.float64 there.

[TristanHehnen]

@JanVogelsang Thank you very much for your input!

We will likely not be able to look at it this week, but rather towards the end of the next one.
In the mean time @Sagar-Singhal you could loop over the individual meshes that have a common cell size, extract the 2d array you are interested in and then merge these together.