Calculate charge by ParticleBinning

[m-rezaei-p]

Dear Experts and Developers of Smilei

Description

I tried to prepare a test for calculating the total charge in AMcylindrical base on this tutorial but with ParticleBinning. https://github.com/SmileiPIC/TP-M2-GI/blob/main/Postprocessing_Scripts/Compute_bunch_parameters.py

Setup: I defined a small cylinder r=10um and h=5 and n0=3E18 cm-3 witch species are:

Species(
    name = 'electron',
    position_initialization = "random",

    momentum_initialization = 'cold',
    particles_per_cell = 4,
    mass = 1.0,
    charge = -1.0,
    charge_density = n0_ncrit_cm,
    mean_velocity = [beta, 0.0, 0.0],
    time_frozen=Lx,

    boundary_conditions =  [
        ["remove", "remove"],
        ["remove", "remove"],
    ],
)

where beta is:

gamma = 100 / 0.511  # 100 MeV electrons
beta = math.sqrt(gamma*gamma - 1.0)/gamma

Also, I defined two diagnostic:

DiagParticleBinning(
    deposited_quantity = "weight",
    every = 10,
    time_average = 1,
    species = ['electron'],
    axes = [["gamma",   190,  200,   10 ] ]
)

DiagTrackParticles(
  species = "electron",
  every = 10,
  attributes = ["x", "y", "z", "px", "py", "pz", "w"]
)

I use the following equation for calculation total charge by DiagTrackParticles Q = total_weight* q * nc * (conversion_factor*1e-6)**3 * 10**(12) # Total charge in pC and I got correct result but when I tried with DiagParticleBinning by following equation: Q_spec =q*count_charge_v * nc * (conversion_factor*1e-6)**3 * 10**(12) where q: electron charge

count_charge_v: sum of DiagParticleBinning for one step 
nc: critical density
conversion_factor: lambda0/2pi (meter)

I got the wrong result. I tried to plot hist with weights of the particle(with panda) with data of TrackParticles. I got the exact number of total_weight of TrackParticles and on another hand, I check the HDF5 file of ParticleBinning and I see the exact value on the table... so that means Happi manipulate data in the wrong way for this purpose...

so can I force Happi, not to covert data? or I must read data directly by h5py? best regards thanks in advance. PS: this is my test file: test.txt

[massim0-1]

Hello, can you provide the commands (including the opening of the diagnostic with happi) that you have used to compute the total charge with Track Particles and with Particle Binning?

What do you see if you make a histogram starting from Track Particles and when you do it plotting the Particle Binning? Can you show the plots with the units (and the scripts to generate them)?

[m-rezaei-p]

Hello, can you provide the commands (including the opening of the diagnostic with happi) that you have used to compute the total charge with Track Particles and with Particle Binning?

What do you see if you make a histogram starting from Track Particles and when you do it plotting the Particle Binning? Can you show the plots with the units (and the scripts to generate them)?

this is my code that I calculate charge and plot histogram: charge_test.txt this is histogram from Track Particles this is HDF5 from DiagParticleBinning

[massim0-1]

Ok, I am reading your script. After the command charge=np.array(S.ParticleBinning(0).getData()), i.e. where you export the data with happi, can you make a print of the array charge? What do you see?

[m-rezaei-p]

Ok, I am reading your script. After the command charge=np.array(S.ParticleBinning(0).getData()), i.e. where you export the data with happi, can you make a print of the array charge? What do you see?

[mccoys]

Happi will give you different units than you think. Plot the results to see the units

[m-rezaei-p]

Happi will give you different units than you think. Plot the results to see the units

when I plot I see this results with this unit:

As I understand mean divided by critical density
Am I right? so your point is I must convert to other unit? but witch unit?

[mccoys]

You should multiply by the bin size and by the box size.

[m-rezaei-p]

You should multiply by the bin size and by the box size.

my bin size in this simulation is 1 but for box size I have a little confuse if h=5um and r=10um what is mean box size: pi h (r^2) or h * r? and in normalize unit ?

I tried this :

L=Q_spec*math.pi*((Radius_plasma)**2)*hr*((1E-4)**3)==> 5.230247758141911e-12
L1=Q_spec*S.namelist.Lr*S.namelist.Lx==>9.985448695563615
L2=math.pi*((S.namelist.Lr)**2)*S.namelist.Lx*Q_spec==>2533.914204616654

and all of them is wrong(correct value is 806.569)

[mccoys]

Ah sorry. I never thought about the AM case. I need some time to figure if happi correctly handles this

[m-rezaei-p]

Ah sorry. I never thought about the AM case. I need some time to figure if happi correctly handles this

Thanks so much

[mccoys]

Ok in the AM case the results are also divided by a volume which is equal to grid_length[0] x grid_length[1]^2

I am adding a line in the happi output which makes this explicit