ParticleBinning diagnostic for pressure

[Tissot11]

I have found inconsistency while calculating Pxx for calculating the temperatures using the Particle Binning diagnostics. I have initialised electrons and ions using the same temperature and density profiles. While I get the correct results for electrons (Txx, Tyy, Tzz) but for ions, I get incorrect results only for Pxx or Txx. Tyy and Tzz for ions are correct too. I have already adjusted for ion mass factor. Both the values and the profile (x<0) for Txx component of ions is incorrect at the initialisation. Txx is supposed to have the same profile like for other components and other species.

I am a bit perplexed by this inconsistency. How to understand this? I do have a drift in x-direction for both electron and ions. Below are the snippets for electrons and ions diagnostics

DiagParticleBinning(
     name = "densEon",
     deposited_quantity = "weight",
     every = everyDiag,
     time_average = 1,
     species = ["eon1", "eon2"],
     axes = [
            ["x", 0.0, Lx, 4096],
            ],
)
DiagParticleBinning(
     name = "TxxEon",
     deposited_quantity = "weight_vx_px",
     every = everyDiag,
     time_average = 1,
     species = ["eon1", "eon2"],
     axes = [
        ["x", 0.0, Lx, 4096],
           ],
)
DiagParticleBinning(
     name = "TyyEon",
     deposited_quantity = "weight_vy_py",
     every = everyDiag,
     time_average = 1,
     species = ["eon1", "eon2"],
     axes = [
          ["x", 0.0, Lx, 4096],
            ],
)
DiagParticleBinning(
     name = "TzzEon",
     deposited_quantity = "weight_vz_pz",
     every = everyDiag,
     time_average = 1,
     species = ["eon1", "eon2"],
     axes = [
          ["x", 0.0, Lx, 4096],
            ],
)
DiagParticleBinning(
     name = "densIon",
     deposited_quantity = "weight",
     every = everyDiag,
     time_average = 1,
     species = ["ion1", "ion2"],
     axes = [
          ["x", 0.0, Lx, 4096],
            ],
)
DiagParticleBinning(
     name = "TxxIon",
     deposited_quantity = "weight_vx_px",
     every = everyDiag,
     time_average = 1,
     species = ["ion1", "ion2"],
     axes = [
          ["x", 0.0, Lx, 4096],
            ],
)
DiagParticleBinning(
     name = "TyyIon",
     deposited_quantity = "weight_vy_py",
     every = everyDiag,
     time_average = 1,
     species = ["ion1", "ion2"],
     axes = [
          ["x", 0.0, Lx, 4096],
           ],
)
DiagParticleBinning(
     name = "TzzIon",
     deposited_quantity = "weight_vz_pz",
     every = everyDiag,
     time_average = 1,
     species = ["ion1", "ion2"],
     axes = [
          ["x", 0.0, Lx, 4096],
            ],
)

[mccoys]

Hi @Tissot11

I moved this to the discussions as I don't think this is a bug. I believe your result is related to the fact that the definition you are using in your Pxx diagnostic (weight_vx_px) does not provide exactly what you seek.

weight_vx_px divided by weight essentially provides the average of $v_{x} p_{x}$ whereas the pressure is correctly defined as the average of $(v_{x}-\langle v_x\rangle) (p_{x}-\langle p_x\rangle)$. As this latter average is equal to $\langle v_x p_x\rangle - \langle v_x\rangle\langle p_x\rangle$, one solution to your problem would be to make separate diagnostics for $\langle v_x\rangle$ and $\langle p_x\rangle$. You can perform the operation during post-processing.

Note of course that $\langle v_x\rangle$ and $\langle p_x\rangle$ are equal to a constant factor when non-relativistic.

[Tissot11]

Thanks. Since the default categories are either Bug or Feature Request for submitting a new issue, I choose Bug.

Both electrons and ions have the same non-relativistic drifts (mean velocities) in x direction, so it is a bit perplexing that electrons diagnostic gives correct results while ion diagnostics does not. Is it because, the ion momentum is scaled in terms of electrons mass scale, making it relativistic in terms of electron momentum? Of course, I can create separate diagnostics and do the stuff in post-processing.

[mccoys]

Yes the difference with electrons might be related to relativistic velocity. I am not sure though.

[Tissot11]

Do you think, this would work?

def stuff(particles):
    return particles.weight * (particle.vx*particles.px)

and subtract this stuff from weight_vx_px in post-processing. Or this stuff is actually weight_vx_px diagnostic?

[mccoys]

Your function will not work as vx does not exist. Note also that the px is a little tricky, as it does not include the mass. In fact, this px is only the lorentz factor multiplied by the velocity.

Anyways, what you have here is the same as weight_vx_px. Instead, you need another diagnostic that simply uses weight_px.

[Tissot11]

Yeah, I realised it after trying this. I dump a 2D diagnostic xPx which have both px and x and I save weight. Doing the same for xVx so x and vx and saving weight should work as the second part in your expression above?

[mccoys]

I don't understand. You want x and px as axes? This is not what you should do. You must have all diagnostics with the same axes (x in your case). Then one diagnostic has deposited_quantity=weight_vx_px and another one has deposited_quantity=weight_px

[Tissot11]

This will help only getting <px>, but how to get <vx> then? I need then another diagnostic for gamma vs x? To be precise, could you please tell how to calculate the second term in the expression you wrote above and whether I should use it a a general procedure for all species?

[mccoys]

If this is non relativistic, px and vx are the same. If you are considering a species with fixed charge, you can actually use weight_charge_vx to get vx. Otherwise you need a custom function like lambda p: p.weight * p.px / (1+p.px**2+p.py**2+p.pz**2)**0.5

[mccoys]

Note that using a custom function prevents changing units in happi

[Tissot11]

I designed two more diagnostics (numDiags=256)

DiagParticleBinning(
     name = "densIon",
     deposited_quantity = "weight",
     every = everyDiag,
     time_average = 1,
     species = ["ion"],
     axes = [
          ["x", 0.0, Lx, numDiags],
            ],
)
DiagParticleBinning(
     name = "TxxIon",
     deposited_quantity = "weight_vx_px",
     every = everyDiag,
     time_average = 1,
     species = ["ion"],
     axes = [
          ["x", 0.0, Lx, numDiags],
            ],
)

DiagParticleBinning(
     name = "pxIon",
     deposited_quantity = "weight_px",
     every = everyDiag,
     time_average = 1,
     species = ["ion"],
     axes = [
          ["x", 0.0, Lx, numDiags],
            ],
)
DiagParticleBinning(
     name = "vxIon",
     deposited_quantity = "weight_charge_vx",
     every = everyDiag,
     time_average = 1,
     species = ["ion"],
     axes = [
          ["x", 0.0, Lx, numDiags],
            ],
)```


and use this 


```pressure = S.ParticleBinning("TxxIon", timesteps=timeSteps[timesteps])
        pressureData = np.array(pressure.getData()[0]).T
        
        pxDiag = S.ParticleBinning("pxIon", timesteps=timeSteps[timesteps])
        pxData = np.array(pxDiag.getData()[0]).T / massratio
        
        vxDiag = S.ParticleBinning("vxIon", timesteps=timeSteps[timesteps])
        vxData = np.array(vxDiag.getData()[0]).T
        
        temperature = (pressureData - pxData* vxData) / (massRatio*densData)

It again produces

it is the same as before.

if I use

        pressureData = np.array(pressure.getData()[0]).T
        
        pxDiag = S.ParticleBinning("pxIon", timesteps=timeSteps[timesteps])
        pxData = np.array(pxDiag.getData()[0]).T 
        
        vxDiag = S.ParticleBinning("vxIon", timesteps=timeSteps[timesteps])
        vxData = np.array(vxDiag.getData()[0]).T
        
        temperature = (pressureData - pxData* vxData) / (massRatio*densData)

I get even a funny result

So something is really strange here for ion species with a non-relativistic drift in x direction. I request you to kindly look into it carefully and I hope it is not a bug..

Please note that I made slight changes in the Namelist, so the fluctuations on the right side are a bit smaller than before.

[mccoys]

The mass must divide the diags with px. You forgot to divide the px and vx dishes individually by the weight diag. And you also forgot to divide one by the charge

[Tissot11]

I am sorry but I am not following you here. Charge is +1, so there is not need to worry about. But what do you mean divide by weight?

Could you please take my expression and write it clearly...

[mccoys]

temperature = (pressureData/densdata - pxData/densdata* vxData/densdata) / massRatio

[Tissot11]

Since I initialize temperature in the units of electron rest mass energy, so I do not need to divide it by mass ratio for checking the results. After realising weight_charge_vx is indeed current density in x, I rescaled it (as you suggested) as temperature = ( pressureData - pxData* (vxData/densData)) / (densData). This indeed gives the correct result.