'cannot convert float NaN to integer' error when opening Track Particles Diagnostic

[Nkehoe-QUB]

Description

When trying to open a Track Particles diagnostic using happi I get the following error:
Error in the ordering of the tracked particles
cannot convert float NaN to integer

>>> import happi
>>> import os
>>> Sim=happi.Open(os.getcwd())
WARNING: 'show=False' was set earlier. Restart python if you want figures to appear.
Loaded simulation '<Path_To_simulation>'
Scanning for Scalar diagnostics
Scanning for Field diagnostics
Scanning for Probe diagnostics
Scanning for ParticleBinning diagnostics
Scanning for RadiationSpectrum diagnostics
Scanning for Performance diagnostics
Scanning for Screen diagnostics
Scanning for Tracked particle diagnostics
Scanning for new particle diagnostics
>>>Sim.TrackParticles('carbon')
Ordering particles ... (this could take a while)
Number of particles: 7500
    Ordering @ timestep = 0
Error in the ordering of the tracked particles
cannot convert float NaN to integer
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/users/40230511/software/Smilei/happi/_Factories.py", line 570, in __call__
    return TrackParticles(self._simulation, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/users/40230511/software/Smilei/happi/_Diagnostics/Diagnostic.py", line 74, in __init__
    remaining_kwargs = self._init(*args, **kwargs)
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/users/40230511/software/Smilei/happi/_Diagnostics/TrackParticles.py", line 115, in _init
    self._orderFiles(orderedfile, chunksize, sort)
  File "/users/40230511/software/Smilei/happi/_Diagnostics/TrackParticles.py", line 457, in _orderFiles
    else        : ordered.fill(self._np.nan)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: cannot convert float NaN to integer

Steps to reproduce the problem

If relevant, provide a step-by-step guide

Parameters

• Smilei version 5.0-65-g5b206cc4e-develop
• Input file ?
• Commands or script used to run the faulty simulation ?
• Computing resources
  • type of computer: HPC
  • OS: Linux
  • Python version: 3.12.0

[Nkehoe-QUB]

DiagTrackParticles(
	species = "carbon",
    every =every,
    flush_every = every * 5,
    attributes = ("x", "y", "px", "py", "w", "q")
)

[charlesprouveur]

The information provided is insufficient to reproduce the ouput above. Please provide the entire input file.

[Nkehoe-QUB]

Full input file:

# ----------------------------------------------------------------------------------------#
#                     SIMULATION PARAMETERS FOR THE PIC-CODE SMILEI                       #
# ----------------------------------------------------------------------------------------#
# All parameters are calculated in SI first, then converted to the normalised units as used by Smilei

import math
import numpy as np

########### Constants ##################################
c = 299792458. 
me = 9.11e-31
epsilon0 = 8.854187e-12
e = 1.602176e-19
amu = 1.673776e-27
Na = 6.022e23
massNeutron = 1838. # in units of electron mass
massProton = 1836.

micro = 1e-6
nano = 1e-9
pico = 1e-12
femto = 1e-15
########################################################

#----------------------------------------------------------------------------------------#
#                           Simulation parameters in SI                                  #
#----------------------------------------------------------------------------------------#

#################### Lasers ############################
lambda_las = 0.8 * micro # wavelength of laser 
omega_las = 2.*math.pi*c / lambda_las # laser frequency - we will be using this as our reference frequency

I_0 = 1.0e21 # Peak intensity of laser in W/cm2

ped_a0 = 0.85 * math.sqrt(((1e-6*I_0)/1e18)*((lambda_las*1e6)**2))
ramp_up = 5 * femto
ped_dur = 2 * pico

Tau_I = 25 * femto # FWHM duration in intensity space
fwhm_I = 3.5 * micro # FWHM spot size in intensity space
AoI_deg = 10 # angle of incidence of laser in degrees
l_ellipticity = 1 # polarisation ellipticity of pulse (0 = LP, +/-1 = CP )
polarization_ang = 0 # angle of the polarization ellipse major axis relative to the X-Y plane, in radians
laser_a0 = ( (l_ellipticity+1)**-0.5) * 0.85 * math.sqrt((I_0/1e18)*((lambda_las*1e6)**2))

#fwhm_E = fwhm_I * math.sqrt(2) # FWHM of spot in field space
w_E = (math.sqrt(2) * fwhm_I ) / (2 * math.sqrt(math.log(2)) ) # beam waist radius
Tau_E = math.sqrt(2) * Tau_I
########################################################


############# Normalisation constants ##################
L_r = c / omega_las # length
T_r = 1 / omega_las # time
den_crit = (me * epsilon0 * omega_las**2) / e**2 # critical density of laser in m^-3 - also density unit in Smilei
E_r = me * (c**2)
P_r = me * c
########################################################


########### General simulation parameters ##############

box_x = 60. * micro # length of simulation box in x
box_y = 20. * micro

res_x = 5. * nano # cell size/resolution in x
res_y = 20. * nano

CFL_factor = 0.99
time_step = (res_x / L_r)*CFL_factor
timestep_SI = time_step * T_r #(CFL_factor * res_x) / c # timestep

# rough number of cells in each direction (will be automatically changed to closest multiple of numpatches, as required by Smilei)
numcell_x = round(box_x / res_x)
numcell_y = round(box_y / res_y)

# SimTime = (box_x/c)+(50 * femto)  # Total simulation time
SimTime = (2. * pico )+ (box_x/c) + (5*femto)
########################################################


##################### Target ###########################
target_thickness = 10. * nano 
contaminant_thickness = 10. * nano
totden = 2000. # mass density in kg/m3
Znum = 6
Anum = 12
x_vac = 10. * micro # position of front edge of target from xmin boundary
y_vac = 2.5 * micro
y_len = box_y - (2 * y_vac)

temp_ev = 1 # temperature in eV (if temp is not set to 'cold')

hydrogen_num = 1 # number of H atoms in molecule
carbon_num = 1 # number of C atoms in molecule

electron_density_SI = (((totden/1000) * Na * Znum)/Anum)/micro
electron_density = electron_density_SI / den_crit

carbon_density_SI = (carbon_num * electron_density_SI) / Znum
carbon_density = carbon_density_SI / den_crit

# proton_density_SI = (hydrogen_num * electron_density_SI) / Znum
# proton_density = proton_density_SI / den_crit
proton_density = 10

# particles per cell for each species
npart_carbon = 5
npart_proton = 5
########################################################


################# Computation parameters ###############
# A function that will find the closest multiple of a number, X, to your chosen number
# You can either add the remainder (+1) or subtract the remainder (-1)
def closest_multiple(num, X, add_or_sub=-1):
    return num + (add_or_sub * (num % X))

# Define number of patches in each direction (must be a power of 2!)
patches_x = 512
patches_y = 128

# Total patches must divide evenly into the number of cells in each direction
true_numcell_x = closest_multiple(numcell_x, patches_x)
true_numcell_y = closest_multiple(numcell_y, patches_y)

box_ycentre = 0.5*true_numcell_y*res_y # in m
########################################################


#----------------------------------------------------------------------------------------#
#                        Simulation blocks - remember to convert units                   #
#----------------------------------------------------------------------------------------#


########### Simulation parameters ######################
Main(
    geometry = "2Dcartesian",

    interpolation_order = 2 ,
    timestep_over_CFL = CFL_factor,#timestep = time_step, # in 2D: dT < dx/sqrt(2)
    simulation_time = SimTime / T_r,

    cell_length = [res_x / L_r, res_y / L_r],     # wavelength/(no. of cells in one wavelength)
    number_of_cells = [true_numcell_x, true_numcell_y],
    #grid_length = [1600.*lambda_las],

    number_of_patches = [ patches_x, patches_y ],       # distribution of processors in the decided geometry

    EM_boundary_conditions = [
        ['silver-muller'],              # open boundary condtion for lasers
        ['silver-muller'],
    ],

	reference_angular_frequency_SI = omega_las,

)

LoadBalancing(
    initial_balance = True,
    every = 150,
)

CurrentFilter(
    model = "binomial",
    passes = [1],
)

Collisions(
    species1 = ["electron"],
    species2 = ["protonF","protonR", "carbon"],
)
######################################################


############## laser definition #####################
LaserGaussian2D(
    box_side        = "xmin",                 # laser originates at left surface
    a0              = ped_a0, # normalized amplitude
    omega           = 1.,  # already in code units, so frequency is just 1
    focus           = [x_vac/L_r, box_ycentre/L_r], # coordinates of laser focus (middle of sim box)
    waist           = w_E / L_r,      # beam waist
    incidence_angle = math.radians(AoI_deg),
    polarization_phi = polarization_ang,      # polarized in the z direction
    ellipticity      = l_ellipticity,
    time_envelope   = ttrapezoidal(start=0., plateau=None, slope1=ramp_up/T_r, slope2=0.),
)

# LaserGaussian2D(
#     box_side        = "xmin",                 # laser originates at left surface
#     a0              = laser_a0, # normalized amplitude
#     omega           = 1.,  # already in code units, so frequency is just 1
#     focus           = [x_vac/L_r, box_ycentre/L_r], # coordinates of laser focus (middle of sim box)
#     waist           = w_E / L_r,      # beam waist
#     incidence_angle = math.radians(AoI_deg),
#     polarization_phi = polarization_ang,      # polarized in the z direction
#     ellipticity      = l_ellipticity,
#     time_envelope   = tgaussian(start=(ped_dur+ramp_up)/T_r, duration=2.4*Tau_I/T_r, fwhm=Tau_E/T_r),
# )
######################################################


################## Target ############################
Species(
    name = 'protonF',
    position_initialization = 'random',       
	momentum_initialization = 'maxwell-juettner',      # becomes maxwell-boltzmann for low temperatures
    particles_per_cell = npart_proton,
    mass = massProton,
    atomic_number = 1,
	ionization_model = 'tunnel',
    ionization_electrons = 'electron',
    charge = 0,
    number_density = trapezoidal(proton_density,
                                 xvacuum=x_vac/L_r, xplateau=contaminant_thickness/L_r,
                                 yvacuum=y_vac/L_r, yplateau=y_len/L_r),
    temperature = [temp_ev/(0.511 * 1e6)],    
	boundary_conditions = [
		["remove", "remove"],           # particles are removed when they cross the boundaries
	],
)

Species(
	name = 'carbon',
	position_initialization = 'random',
	momentum_initialization = 'maxwell-juettner',
	particles_per_cell = npart_carbon,
	mass = (6*massNeutron)+(6*massProton),
	atomic_number = 6,
	ionization_model = 'tunnel',
    ionization_electrons = 'electron',
	charge = 0,
	number_density = trapezoidal(carbon_density,
                                 xvacuum=(x_vac + contaminant_thickness)/L_r, xplateau=target_thickness/L_r,
                                 yvacuum=y_vac/L_r, yplateau=y_len/L_r),
	temperature = [temp_ev/(0.511 * 1e6)],
	boundary_conditions = [
		["remove", "remove"],
	],
)

Species(
    name = 'protonR',
    position_initialization = 'random',       
	momentum_initialization = 'maxwell-juettner',      # becomes maxwell-boltzmann for low temperatures
    particles_per_cell = npart_proton,
    mass = massProton,
    atomic_number = 1,
	ionization_model = 'tunnel',
    ionization_electrons = 'electron',
    charge = 0,
    number_density = trapezoidal(proton_density,
                                 xvacuum=(x_vac + contaminant_thickness + target_thickness)/L_r, xplateau=contaminant_thickness/L_r,
                                 yvacuum=y_vac/L_r, yplateau=y_len/L_r),
    temperature = [temp_ev/(0.511 * 1e6)],    
	boundary_conditions = [
		["remove", "remove"],           # particles are removed when they cross the boundaries
	],
)

Species(
	name = 'electron',
	position_initialization = 'random',
	momentum_initialization = 'maxwell-juettner',
	particles_per_cell = 0,
	mass = 1.,
	charge = -1.0,
	number_density = 0,
	temperature = [temp_ev/(0.511 * 1e6)],
	boundary_conditions = [
		["remove", "remove"],
	],
)
#####################################################


############### Diagnostics ##########################
every_SI = 10 * femto # How often do you want to output the data
every = int(every_SI / timestep_SI) # convert to number of timesteps
every_hi_res = every/2

# Max energy limit detected by diagnostics, numbers in MeV/n and converted to SI
MeV_to_J = 1.6e-13

maxEnP = 400 * MeV_to_J
maxEnC = 100 * 12 * MeV_to_J
maxEnE = 100 * MeV_to_J

# momentum limits
maxP_P = math.sqrt(2 * amu * maxEnP)
maxP_C = math.sqrt(2 * 12 * amu * maxEnC)
maxP_E = math.sqrt(2 * me * maxEnE)

xspacebins = 1500
yspacebins = 500
phasebins = 500
specbins = 500

from numpy import s_ # This enables us to select a subgrid for our fields, reducing disk usage for large grids


DiagTrackParticles(
	species = "protonF",
    every =every,
    flush_every = every * 5,
    attributes = ("x", "y", "px", "py", "w", "q")
)

DiagTrackParticles(
	species = "protonR",
    every =every,
    flush_every = every * 5,
    attributes = ("x", "y", "px", "py", "w", "q")
)

DiagTrackParticles(
	species = "carbon",
    every =every,
    flush_every = every * 5,
    attributes = ("x", "y", "px", "py", "w", "q")
)

DiagTrackParticles(
	species = "electron",
    every =every,
    flush_every = every * 5,
    attributes = ("x", "y", "px", "py", "w", "q")
)

[charlesprouveur]

I don't think this is related to track particles. If i try to run the above input file on one core i see:

Stack trace (most recent call last):
#5    Object "[0xffffffffffffffff]", at 0xffffffffffffffff, in 
#4    Object "./smilei", at 0x559b942a697d, in _start
#3    Object "/lib/x86_64-linux-gnu/libc.so.6", at 0x7f468dd2d082, in __libc_start_main
#2    Object "./smilei", at 0x559b942a4f42, in main
#1    Object "./smilei", at 0x559b9469f335, in VectorPatch::applyBinaryProcesses(Params&, int, Timers&)
#0    Object "./smilei", at 0x559b942c0487, in BinaryProcesses::apply(Params&, Patch*, int, std::vector<Diagnostic*, std::allocator<Diagnostic*> >&)
Floating point exception (Integer divide by zero [0x559b942c0487])

Which is what probably created the NaN you are seing. Now why there is a division by zero is the question.
After running on 4 cores: computation run until a Floating point exception is encountered on rank 0.

[mccoys]

@Nkehoe-QUB You might have run into issue #712. Please update and see if it works now.

[charlesprouveur]

indeed running with the latest version on my side seems to solve the issue

[Nkehoe-QUB]

Will update and run it again

[Nkehoe-QUB]

Upgraded to version: 5.0-206-g456444957-master
I still get the error:

Ordering particles ... (this could take a while)
Number of particles: 7500
   Ordering @ timestep = 0
Error in the ordering of the tracked particles
cannot convert float NaN to integer 

When running:

Sim = happi.Open(Folder)
Sim.TrackParticles('carbon')

[charlesprouveur]

For clarification:

There are no error message during the execution.
I was able to reproduce the error with this downscaled version of the input file:
input.pdf
( rename the .pdf in a .py, this is not an actual pdf)
Among the diagnostics, carbon and both proton have the nan error; for electron this raise an exemption due to no particles.
Looking at the datasets: the "w" property is the issue (with h5ls -d TrackParticles_carbon.h5/w we see a bunch of NaN.
Adding a field diagnostics provides no more information, although i am surprised to only see zeros for iteration 1 and 2

[mccoys]

Maybe the particles really have nans due to overflows?

[Nkehoe-QUB]

For clarification:

There are no error message during the execution. I was able to reproduce the error with this downscaled version of the input file: input.pdf ( rename the .pdf in a .py, this is not an actual pdf) Among the diagnostics, carbon and both proton have the nan error; for electron this raise an exemption due to no particles. Looking at the datasets: the "w" property is the issue (with h5ls -d TrackParticles_carbon.h5/w we see a bunch of NaN. Adding a field diagnostics provides no more information, although i am surprised to only see zeros for iteration 1 and 2

The electron species come from ionisation so iteration 1 and 2 will have no electrons present yet. I'm not sure as to why the "w" property has NaN values. I used one of the standard number density profiles.

[mccoys]

Ok it turns out it is an issue about numpy. Somewhere between version 1.20 and 1.24, it has lost the capability to handle nan integers. I will look into it

[mccoys]

Should be solved in branch develop. Please reopen if necessary