TrackParticles related issue

[Sonuphd]

Respected Sir,

I am facing post-processing related issue for my simulation. I am sending you input file and commands. import math

l0 = 2.0math.pi # laser wavelength [in code units] t0 = l0 # optical cycle Lsimx = 300.l0 # length of the simulation box in x (= 240 mum) Lsimy = 80.l0 # length of the simulation box in y (= 64 mum) Tsim = 1500.t0 # duartion of the simulation (= 4 ps) resx = 16. resy = 2. # nb of cells in on laser wavelength rest = resx/0.8 # nb of timesteps in one optical cycle dt = t0/rest nx = 300*16 # nb of cells along x-axis npatchx = 32 # nb of patches along x-axis

-----------

n0max = 5.734e-3 # electron density (code units =>1=plasma critical density)

initial density profile of electrons

def n0_electron(x,y):

return 1./(1.+math.exp((abs(x-x0)-widthx)/Lx))*n0max

------------

DEFINING SMILEI's VARIABLES

All in "bloacks"

Main( geometry = "2Dcartesian", interpolation_order = 2, timestep = dt, simulation_time = Tsim, cell_length = [l0/resx,l0/resy], grid_length = [Lsimx,Lsimy], number_of_patches = [32,16], clrw = nx/npatchx, reference_angular_frequency_SI = 2.0math.pi3e8/0.8e-6, EM_boundary_conditions = [["silver-muller","silver-muller"],["silver-muller","silver-muller"],], solve_poisson = False, print_every = 100, random_seed = smilei_mpi_rank )

MovingWindow( time_start = 0.9*Main.grid_length[0], velocity_x = 0.9997 )

LoadBalancing( initial_balance = False, every = 20, cell_load = 1.,

)

Species( name = "electron", position_initialization = "regular", momentum_initialization = "cold", particles_per_cell = 4, mass = 1.0, charge = -1.0, number_density = 0.005734,# electron density (code units =>1=plasma critical density) mean_velocity = [0.0,0.0,0.0], pusher = "boris", time_frozen = 0.0, boundary_conditions = [["remove","remove"],["remove","remove"],], )

FWHMinI = 12.5l0 #(= 10 mum) waistinI = FWHMinI/(2.0math.sqrt(math.log(2.0))) waistinE = FWHMinI/math.sqrt(2.0math.log(2.0)) FWHMtinI = 8.0t0 #(= 21.3 fs) FWHMtinE = FWHMtinI*math.sqrt(2.0) #(= 30 fs)

diagEvery = int(37.5*t0/dt) # frequency of outputs(= 100 fs)

LaserGaussian2D( box_side = "xmin", a0 = 5.2, # intensity 5.8 e+19 omega = 1., focus = [0.,Main.grid_length[1]/2.], waist = waistinE,

time_envelope = tgaussian(start=0.,duration=2.*FWHMtinE,fwhm=FWHMtinE,center=FWHMtinE)

) list_fields = ['Ex','Ey','Bz','Rho','Jx']

DiagFields( every = diagEvery, fields = list_fields )

DiagProbe( every = diagEvery, origin = [0.,Main.grid_length[1]/2.], corners = [[Main.grid_length[0],Main.grid_length[1]/2.], ], number = [nx], fields = list_fields )

DiagScalar( every = int(diagEvery/10), vars=['Uelm','Ukin_electron', 'ExMax','ExMaxCell','EyMax','EyMaxCell','RhoMin', 'RhoMinCell','Ukin_bnd','Uelm_bnd','Ukin_out_mvw','Ukin_inj_mvw','Uelm_out_mvw','Uelm_inj_mvw','Utot'] )

DiagParticleBinning( deposited_quantity = "weight", every = diagEvery, species = ["electron"], axes = [ ["moving_x",0., Lsimx, 300], ["ekin",1.,500.,200] ] ) DiagParticleBinning( deposited_quantity = "weight", every = diagEvery, species = ["electron"], axes = [ ["moving_x",0., Lsimx, 300], ["px",-1.,1.,100] ] ) DiagParticleBinning( deposited_quantity = "weight", every = diagEvery, species = ["electron"], axes = [ ["y",0., Lsimy, 300], ["py",-1.,1.,100] ] )

DiagParticleBinning( deposited_quantity = "weight", every = diagEvery, species = ["electron"], axes = [ ["moving_x",0., Lsimx, 300], ["px",-1.,1.,100] ] ) DiagParticleBinning( deposited_quantity = "weight", every = diagEvery, species = ["electron"], axes = [ ["y",0., Lsimy, 300], ["py",-1.,1.,100]]

    )

DiagTrackParticles( species = "electron", every = diagEvery,

flush_every = 100,

filter = my_filter,

attributes = ["x", "px", "py", "Ex", "Ey", "Bz"]

) DiagPerformances( every = diagEvery

Post-Processing Command: In [8]: diag=S.TrackParticles("electron",axes=["x","y"])
/home/sonu/Smilei/happi/_Diagnostics/TrackParticles.py:45: H5pyDeprecationWarning: The default file mode will change to 'r' (read-only) in h5py 3.0. To suppress this warning, pass the mode you need to h5py.File(), or set the global default h5.get_config().default_file_mode, or set the environment variable H5PY_DEFAULT_READONLY=1. Available modes are: 'r', 'r+', 'w', 'w-'/'x', 'a'. See the docs for details. f = self._h5py.File(orderedfile) Selecting particles ... (this may take a while) Removing dead particles ... Kept 15744000 particles

In [9]: diag.plot()
Loading data ... Killed

[mccoys]

You don't have enough memory for plotting 16 billion particles !!!

[Sonuphd]

Respected Sir,

How much memory need for track particles diagnostic.

[mccoys]

I don't know how much memory is needed. The problem occurs in plotting which is difficult to guess. But you are trying to plot 16 billion lines which is completely impossible in the (relatively) basic utility matplotlib that smilei uses. Usually it is impossible to plot more than 1000 lines.

Do you REALLY need to plot those lines ? What do you expect to see ? Even 20 lines is something very hard to understand on one figure. So what is your intent with 16 billion ????

[Z10Frank]

@Sonuphd the tracked particles are too many because you probably need to add a filter in your TrackParticle diagnostic, otherwise it will track all the particles of Species "electron" in the simulation domain, which is completely useless and memory consuming if you simulate LWFA:

DiagTrackParticles(
species = "electron",
every = diagEvery,
flush_every = 100,
filter = my_filter,
attributes = ["x", "px", "py", "Ex", "Ey", "Bz","weight"]
)

Typically in LWFA you want to track only particles above a certain energy (the accelerated bunch for example), so you can add a threshold on px for the particles to be tracked. These lines must be added before the TrackParticles block using my_filter:

def my_filter(particles):
    return (particles.px>10)

With this example you will track only particles with a longitudinal momentum px> 10 m_e*c. This way you will ignore most of the electrons in the plasma background, and keep only the most energetic ones.

If they are still too many try to increase the px threshold or just sample a subset of the particles, e.g. every 5 particles, to see only representative trajectories and not all of them, which is again useless and time consuming.

Please read carefully the doc on TrackParticles and relative filtering https://smileipic.github.io/Smilei/namelist.html#diagtrackparticles

[Sonuphd]

Thank you sir, problem solved.

[tanjhiop]

@Z10Frank Hi Francesco, so can we track tunneling ionized high energy electron ? Say ionization injection LWFA, filtering to obtain the IDs and restart to track them from the 'birth' to end ? Thanks in advance !

[Z10Frank]

@tanjhiop yes, the electrons generated from ionisation are a separated Species (in the benchmark folder see e.g. the namelist tstAM_14_envelope_ionization_linear_polarization.py, where they are called electronfromion), so for LWFA I normally put a Track diagnostic only on that Species.

Again be careful because without defining an appropriate filter the diagnostic will store all the electrons created from ionization, often saturating the memory. As filter you may want to define a certain threshold on energy, and/or a certain maximum distance from the laser propagation axis. If you don't put a filter you will already see them from birth to end.

[tanjhiop]

@Z10Frank Thanks for you timely reply! I did set a proper filter to obtain a memory acceptable amout of electrons, however, i fail to ensure exactly same IDs in two individual runs for electrons pass the same filter. To my understanding of Smilei particle identification, paritcles will only get their IDs when they pass the filter. So, how can i manage to ensure same electronsfromion identification with two different filters (one for energy threshold, second for IDs) ?

[mccoys]

@tanjhiop You can filter all particles on the first iteration so that they all have an ID, then make your filter on the subsequent iterations.

[tanjhiop]

@mccoys Thanks for your reply !

I am still confused and have two more questions. 1. Electrons to be tracked are from tunneling ionization, which are created continuouly along the laser propagation. So it is hard for me to find the time when all the electrons interested in have been created thus they can be filtered and obtain IDs. 2. And the IDs filtered in the first iteration will still be effective in subsequent iterations with different filter ?

[mccoys]

1. I don't know what is your criterion on which electrons you are interested in so I don't know how to answer. If you want all ionised electrons, simply track only the ionised electron species (which can be made separate)

2. Don't do a different filter. Only 1 filter can be made anyways. You must tell the filter to select ALL particle initially (so that they all have an ID that they will keep), then continue the filtering on the particles you want. Example:

   def myfilter(particles):
   if Main.iteration == 0:
       return particles.x == particles.x # always true
   else:
       return particles.px > 0.5 # Only particles with px > 0.5 are selected

[tanjhiop]

@mccoys

1. Indeed i am interested in electron with energy above a certain threshold, which can be easily selected.
2. I think the underlying basis of the filter you suggested is ALL electrons involved have been created at the start. However, in terms of tunnling ionized electrons, when Main.iteration == 0, most electrons with energy above a threshold at some later time Main.iteration == XXXX are not created at all when Main.iteration ==0, thus they can not pass the filter and acquire an ID. So i think this prevent me getting ALL electrons identified at the start of the simulation. I doubt whether i am right, correct me please !

[Z10Frank]

@tanjhiop of course you will need to adapt @mccoys 's suggestion to your needs, and again it depends on what you want to do with those electrons.

If you want to know the properties like energy spread, emittance, etc of the electrons that are trapped and accelerated in the plasma wave, just put a filter as a px threshold and you'll get all the high energy electrons of the bunch, which will be equivalent to all the bunch particles normally. You won't see them from their creation, but when you study those properties it is useless to compute them from the start, so it should give you the relevant information.

If instead you want to do something else which requires to store them since their creation, like plotting the trajectories of the electrons from their creation, it is useless to keep all of them. You must have a general idea of when injection should start or has already started (no need to be precise), let's say iteration Y. Then, use that iteration in @mccoys 's suggestion instead of iteration 0, adapting the px threshold to your needs for iterations after Y.

[tanjhiop]

@Z10Frank Yes, figuring out a proper iteration Y is what i am trying. Thank you and @mccoys so much !

One small piece of possible idea to facilitate the particle identification and track process, maybe we can construct a filter return two boolean arrays at the same time, one for particle ID generation (all that are interested.), another for a subset that are really tracked.

[mccoys]

Ok you are right concerning new electrons being created. You can also detect those because they have a default ID == 0. I agree that it makes the process a bit complicated.