FireDynamics / fdsreader

Python reader for FDS data
GNU General Public License v3.0
44 stars 18 forks source link

Problem with particle quantities #71

Closed r-broek closed 5 months ago

r-broek commented 5 months ago

Dear developers,

First of all, thank you for the useful data extraction tool! I am currently investigating the behaviour of individual particle trajectories and quantities. I have written a function that finds the (x y z) position of a particle from class small roundwood with its tag particle_id:

sim = fds.Simulation("..")

def getPositions(particle_id):
    quant = sim.particles['small roundwood'].positions
    partquant = np.zeros((len(quant),3))
    for i in range(len(quant)):
        index = np.where(sim.particles['small roundwood'].tags[i]==particle_id)
        if len(index[0])==1:
            partquant[i] = quant[i][index[0][0]]
        else:
            print('particle left simulation at i=%i'%i)
            return partquant,i 
            break
    return partquant, -1

This seems to work properly, as the trajectories are continuous and similar to Smokeview. However, I am having trouble with the particle quantities, which I obtain through a similar function:

def getQuantity(particle_id, quantity):
    quant = sim.particles['small roundwood'].data[quantity]
    partquant = np.zeros(len(quant))
    for i in range(len(quant)):
        index = np.where(sim.particles['small roundwood'].tags[i]==particle_id)
        if len(index[0])==1:
            partquant[i] = quant[i][index[0][0]]
        else:
            print('particle left simulation at i=%i'%i)
            break
    return partquant

The returned data (e.g. 'PARTICLE TEMPERATURE') displays discontinuous behaviour and differs from what is observed with Smokeview. Furthermore, other quantities such as 'PARTICLE DIAMETER' seem to load temperature data, as the values are in the same range as those from the temperature, while they should be orders of magnitude different. It should be noted that Smokeview does properly load the other particle quantities. I have tested this for FDS version 6.7.7 and 6.8.0. Also, I have not altered the storage format, so they should be 4-byte and therefore the default FDSReader settings should work.

I am not sure whether I have simply made a mistake in my code, or if there is a bug in the FDSReader module. The example particle file did not give me any new insights, and although the Fire Simulation lecture notes were helpful, I was still not able to fix it. Hopefully, you could help me resolve the issue.

Any help is much appreciated!

JanVogelsang commented 5 months ago

Thanks for bringing this up! In order to test particles in more detail, I replaced the particle example by the water droplets verification example from the FDS repository and unfortunately was not able to reproduce your issue (FDS 6.8.0). Could you share your fds-case and python-script for me to reproduce the issue?
Also be aware of the filter_by_tag function that we provide. I have to admit I never used that function yet, so I never tested it on any example other than a small example case which I used for debugging purposes when writing the particle code. So while I would expect the function to work, please properly test it (verify the correction of results) before extensively using it.

JanVogelsang commented 5 months ago

I tested the filter_by_tag function on the water droplets verification case and it seems to do the correct thing. Please test this function in your case as well and in case it still doesn't work just proceed by sending your scripts. As everything seems to work fine in a typical example case, one should be able to easily find the differences between the two cases and see why that causes issues in the fdsreader.

r-broek commented 5 months ago

Thank you for looking into the issue. Due to your comment, I realized that the filter_by_tag function is exactly what I needed and fully replaced the functions I made myself.

My Python script is now as follows:

import numpy as np
import matplotlib.pyplot as plt
import fdsreader as fds

rng = np.random.default_rng(seed=None)
part_id = 'small roundwood'

sim = fds.Simulation("..")
times = np.array(sim.particles[part_id].times)

tags = sim.particles[part_id].tags[0]
tag = tags[rng.integers(0,len(tags))]
print(tag)

PART1 = sim.particles[part_id].filter_by_tag(tag)

TEMP = PART1.data['PARTICLE TEMPERATURE']
MASS = PART1.data['PARTICLE MASS']

plt.plot(times,TEMP)
plt.plot(times,MASS)
plt.show()

The filter_by_tag function seems to work as intended; it can always find the tag for the specified particle class. It will not take a particle from a different class. I have tested this by plotting the PARTICLE TEMPERATURE for large roundwood, which FDS does not store, and the script would therefore throw an error. Furthermore, for small roundwood the function returns the expected data: Particle(name=small roundwood, quantities=[Quantity('PARTICLE TEMPERATURE'), Quantity('PARTICLE MASS'), Quantity('PARTICLE DIAMETER'), Quantity('PARTICLE BULK DENSITY'), Quantity('PARTICLE HEAT TRANSFER COEFFICIENT'), Quantity('PARTICLE CONVECTIVE HEAT FLUX'), Quantity('PARTICLE RADIATIVE HEAT FLUX')])

I have also run the Verification FDS simulation water_evaporation_2 and plotting the particle quantities worked without any of the issues. This is what you observed as well. However, plotting quantities of my simulation resulted in similar problems as before. The data seems to be incorrect. For example, in the graph below the supposed particle temperature and mass are shown. The behaviour (and order of magnitude) of the particle mass is completely wrong. We would expect it to decrease. This figure was obtained using the Python code on the results from the simulation that is included as an attachment. The particles are set up in the _veg.txt file (the extension should .txt). The FDS file extension was altered in order to upload it to GitHub. The simulation has 4 meshes and should ideally be run on 4 nodes.

I thought that maybe the EMBER_PARTICLE=T tag given to the particles would change their particle tag or that if particles escaped the simulation domain, the tags would shift. However, neither is the case (I asked this on the FDS Q&A forum). I do not know what else could be an issue, besides having more than one particle class in the simulation.

Hopefully, you have some additional insights. I really appreciate your help.

image

briannak7 commented 5 months ago

Hey all,

I am seeing a similar problem with the particle quantities as r-broek.

This simulation contains two particles "TREATMENT PARTICLE" and "CONTROL PARTICLE". A fire is ignited after 30 seconds and the time starts at -30 to account for this.

The "CONTROL PARTICLE" is a rectangle filling out most of the boundary ground (height of 0.6 m) with an empty circle in its center.

"TREATMENT PARTICLE" (height of 0.1 m) is particles located in the empty circle of "CONTROL PARTICLE".

The mass and temperature of both particles have very similar data, and is not what is expected/seen within Smokeview.

import matplotlib.pyplot as plt

treat_part = sim.particles[0]
control_part = sim.particles[1]

treat_data = treat_part.data
control_data = control_part.data

treat_temp_data = treat_data["PARTICLE TEMPERATURE"] 
treat_mass_data = treat_data["PARTICLE MASS"] 

max_ = 0
min_ = 0

mass = []
temp = []

for i in range(len(treat_mass_data)):
    mass_array = treat_mass_data[i]
    temp_array = treat_temp_data[i]

    mass.append(mass_array.sum())
    temp.append(temp_array.sum())

plt.plot(mass, label="mass")
plt.plot(temp, label="temp", alpha=0.5)
plt.legend()
plt.show()

particle_mass_v_temp

I went ahead and plotted the particle temperature data for a time step 0 and similarly found that the data does not match what is expected.

import matplotlib.pyplot as plt
%matplotlib inline

i = 0

treat_part = sim.particles[0]
control_part = sim.particles[1]

treat_data = treat_part.data
control_data = control_part.data

treat_temp_data = treat_data["PARTICLE TEMPERATURE"] 
control_temp_data = control_data["PARTICLE TEMPERATURE"]

x_control, y_control = control_positions[i][:, 0], control_positions[i][:, 1]
x_treat, y_treat = treat_positions[i][:, 0], treat_positions[i][:, 1]

# treatment data plot
plt.scatter(x_treat, y_treat, c=treat_temp_data[i], cmap='viridis', s=100, marker='o')
plt.colorbar(label='Temperature')
plt.xlabel('X Position')
plt.ylabel('Y Position')
plt.title('Temperature Distribution treatment')
plt.show()

# control data plot
plt.scatter(x_control, y_control, c=control_temp_data[i], cmap='viridis', s=100, marker='o')
plt.colorbar(label='Temperature')
plt.xlabel('X Position')
plt.ylabel('Y Position')
plt.title('Temperature Distribution control')
plt.show()

treatment_temp_timestep_0

control_temp_timestep_0

Since a fire ignition does not start for 30 seconds after the simulation begins, we do not expect to see difference in temperature of any of the particles at timestep 0. Smokeview also shows no temperature difference until the fire is ignited.

Smokeview image at time step 0:

Screenshot 2024-03-14 at 16 02 46

I saw the recent comments about filter_by_tag and quickly tried this. I still found inconsistent data when plotting the mass vs temperature sums. This may be error on my part and I can dive more into this and provide feedback/images if it would be helpful!

I have a suspicion that this error may be related to an update with FDS and found an issue on firemodels/fds that may point to the specific update: https://github.com/firemodels/fds/issues/12529

FDS version used:

Current Date : March 12, 2024 09:58:39 Revision : FDS-6.8.0-1630-gbc9da78-master Revision Date : Tue Feb 27 12:28:41 2024 -0500 Compiler : GCC version 10.2.0 Compilation Date : Feb 28, 2024 08:59:25

JanVogelsang commented 5 months ago

@r-broek Thanks for sending me your simulation script! Unfortunately, if I run your python-script the temperature and mass seem to be quite high, but fall steadily, which I think is what you expect. Could you tell me the seed for a particle that produces results you wouldn't expect? Unfortunately, if I try to load the particles in SmokeView, it crashes with a BadAccess error message and I couldn't find any solution for that crash.
However, when looking at the raw data that the fdsreader loads, it's easy to spot a mistake. I tried to find differences to the FDS verification cases and while there are several differences, I think the number of quantities might be the one with the highest potential of causing issues. I removed all other quantities apart from temperature and the data looks much more reasonable. Could you verify that this is indeed the case? As SmokeView still crashes with a single quantity, it's difficult for me to verify the results. As soon as I also add mass to the output quantities, the results don't look right again. In fact, it was trivial to show that different data for temperature is output depending on the number of other quantities that are also output. This clearly shows that adding quantities is indeed the issue (as the number of quantities also influences how the binary data FDS outputs is interpreted). To now continue debugging, it would be great to know if the temperature data is correct if no other quantities are output. If that is the case, that will help me to scan the binary output of FDS for the correct values to interpret where the data extraction method currently implemented is buggy.
If possible, could either you @r-broek or you @briannak7 create a minimal example case which also produces incorrect output? Ideally, with just a single particle, that would be perfect, however simply having lower numbers of particles already helps a lot.

r-broek commented 5 months ago

@JanVogelsang I have quickly made two test cases, one contains a single particle of class 'small roundwood' and the other contains two. In both cases, these particles are above a fire. The expected behaviour is, therefore, loss of particle mass and increase of particle temperature, which are both tracked. This is observed in Smokeview. In the 1 particle case, this is also obtained via FDSReader. In the 2 particle case, the issue arises.

Regarding the Smokeview issue, I have not encountered this issue myself. Are you using SMV version 6.8.0?

(the simulations have now only 1 mesh, so run it using a single node) 2m_M14_veg.txt 2part.fds.txt 1part.fds.txt

JanVogelsang commented 5 months ago

Awesome, thanks for preparing these test cases, they will certainly help in identifying the bug!

JanVogelsang commented 5 months ago

@briannak7 @r-broek I fixed the issue, please update to the 1.10.3 and try out if it works for you! Thanks for finding the bug and helping me fix it, that was much appreciated!

r-broek commented 5 months ago

Thank you for looking into it! Unfortunately, I could not update FDSReader to the newer version via pip. After using pip install --upgrade fdsreader, some of the .py files would be empty.

I have tried installing FDSReader manually by downloading the fdsreader directory from GitHub, however, testing the two case simulations still had the same results as before.

My Python script is now as follows:

import numpy as np
import matplotlib.pyplot as plt
import fdsreader as fds

rng = np.random.default_rng(seed=None)
part_id = 'small roundwood'

sim = fds.Simulation("..")
times = np.array(sim.particles[part_id].times)

tags = sim.particles[part_id].tags[0]
tag = tags[rng.integers(0,len(tags))]
print(tag)

PART1 = sim.particles[part_id].filter_by_tag(tag)

TEMP = PART1.data['PARTICLE TEMPERATURE']
MASS = PART1.data['PARTICLE MASS']

plt.plot(times,TEMP)
plt.plot(times,MASS)
plt.show()
JanVogelsang commented 5 months ago

Installing from source was expected to not work due to an oversight from my side. But the pypi package should have worked, no idea why all files in the wheel are empty.

JanVogelsang commented 5 months ago

Found the issue, please try 1.10.4!

r-broek commented 5 months ago

It works as intended now, thanks so much for fixing this quickly!