ComputationalRadiationPhysics / picongpu

Performance-Portable Particle-in-Cell Simulations for the Exascale Era :sparkles:
https://picongpu.readthedocs.io
Other
689 stars 216 forks source link

Energy Density : e_all_energyDensity #3789

Open cbontoiu opened 2 years ago

cbontoiu commented 2 years ago

Hello,

I am trying to make use of the field 'e_all_energyDensity'. If I understand correctly, the manual says that within each mesh cell, this is the sum of the kinetic energies divided by the cell volume. That is all particles that happen to be in a cell will have their kinetic energy added up and the result will be divided by the cell volume.

In the figure below, plotting 'e_all_energyDensity' shows a tiny region where this field has a maximum of about 0.7e+17 (I assume in Joule/m^3)

image

and because my cell volume is CELL_VOLUME = 1.33504e-30 m^3 I can multiply 'e_all_energyDensity' by CELL_VOLUME to get the energy to about 1.0e-13 (I assume in Joule)

image

This can then be divided by the electron charge e to get the energy in eV so I can finally say that the tiny electron bunch has an average kinetic energy of ~0.60 MeV.

image

Is my reasoning right? Does the choice of particle per cell (6 PPC here) influence the result somehow, I mean do I have to use 6 in my conversions ? Does it make a difference if I dump all particles or only probes, as it is the case here with using ProbeEveryTwentiethCell = EveryNthCellImp l< mCT::UInt32 < 20, 20, 20 > >; ?

Many thanks for attention. Cristian

sbastrakov commented 2 years ago

Hello @cbontoiu . All data we output is openPMD-annotated. So each field should have a unitSI attribute to convert them to a value in SI. Out of my head I do not know if your interpretation happen to match it in this case. However when working with PIConGPU output in general, it's best to not assume any normalization (as it depends on the internally used units) and use this metainformation to interpret the data,

Number of particles per cell should not matter for it in the SI system, as it is only the "sampling resolution" of PIC method, and not a physical quantity like the distribution density.

sbastrakov commented 2 years ago

The second part of the question I do not quite understand. You are now only outputting the energy density for electrons (via prefix e_). So unless ProbeEveryTwentiethCell is somehow used for electrons, it should have no influence on it at all? Sorry if i'm missing something

cbontoiu commented 2 years ago

About the second part of the question I realize now that the probes or particles dumped through the fileOutput.param are not related to the 'e_all_energyDensity'. This is my file

image

cbontoiu commented 2 years ago

About the metadata, I think this refers to the second field called info_Epe in the line below

image

but I don't know how to use it, so I dangerously convert like this

image

where CA[16] is 1e-06 to convert the energy into MeV.

sbastrakov commented 2 years ago

Ah, originally I didn't realize you were using openPMD viewer and not reading "raw" data from our output. The viewer seems to internally use that unitSI before displaying the data to you. And so then your conversion from J/m^3 to MeV should be fine.

However, keep in mind I am no physicist. Does it look good to you @PrometheusPi ?

PrometheusPi commented 2 years ago

@cbontoiu Please excuse the late reply. Regarding:

I am trying to make use of the field 'e_all_energyDensity'. If I understand correctly, the manual says that within each mesh cell, this is the sum of the kinetic energies divided by the cell volume. That is all particles that happen to be in a cell will have their kinetic energy added up and the result will be divided by the cell volume.

This it totally correct. I would just add, that since we are using macro-particles with a spatial shape, neighboring particles will also influence this sum over all particles. (see https://picongpu.readthedocs.io/en/latest/models/shapes.html for more details)

Regarding:

(1) In the figure below, plotting 'e_all_energyDensity' shows a tiny region where this field has a maximum of about 0.7e+17 (I assume in Joule/m^3) (2) and because my cell volume is CELL_VOLUME = 1.33504e-30 m^3 I can multiply 'e_all_energyDensity' by CELL_VOLUME to get the energy to about 1.0e-13 (I assume in Joule) (3) This can then be divided by the electron charge e to get the energy in eV so I can finally say that the tiny electron bunch has an average kinetic energy of ~0.60 MeV. (4) Is my reasoning right?

To (1): as far as I understand your code in your later comment you did not convert to SI units - thus your energy density is still in PIConGPU units. (Assuming you are using the openPMD-api) Perhaps this is defined in CA[4]- I do not know what this variable means.

To (2): If your previous result would have been in SI units, that step would be correct.

To (3): (Assuming you where is SI units before.) No, the value is not an average value. The number you calculated is the total kinetic energy within the volume of a cell. To convert this to an average energy of the electrons you need to dived that total energy per cell by the number of electrons per cell. This is however a dynamic quantity. I would recommend to use the charge density here. This would also allow you to avoid the multiplication with the cell volume and the division by the elementary charge (to convert from J to eV) .

To (4): your general reasoning is correct, but there were some minor errors.

Regarding:

Does the choice of particle per cell (6 PPC here) influence the result somehow, I mean do I have to use 6 in my conversions ? Does it make a difference if I dump all particles or only probes, as it is the case here with using ProbeEveryTwentiethCell = EveryNthCellImp l< mCT::UInt32 < 20, 20, 20 > >; ?

The numbers per cell that you initialize should not matter if the analysis is set up correctly. However, the number of electrons per cell (dynamic quantity) is essential for computing the average (please see "To (3)".)

sbastrakov commented 2 years ago

@PrometheusPi i thought the openPMD viewer always gives SI, as opposed to manually opening our output files via openPMD API?

cbontoiu commented 2 years ago

@PrometheusPi CA[4] is only a numeric conversion and equals 1e15 for time, that is I think in time = 2 meaning femtosecond but for openPMD viewer time needs to be in seconds as time = 2e-15. It shouldn't affect the results at all.

Thank you for your kind reply. I will read it carefully. Kind reagrds

PrometheusPi commented 2 years ago

@sbastrakov You are right :+1: - I thought @cbontoiu is using the openPMD-api and not the openPMD-viewer. (Is this correct @cbontoiu?)

If you are using the openPMD-viewer interface, my comment on "To (1):" can be ignored and your only error occurs at computing the average (see "To (3):").

cbontoiu commented 2 years ago

Yes, I always use openPMD viewer.

PrometheusPi commented 2 years ago

Perfect :+1: - than the only issue I see is how you define "average" - but perhaps I misunderstand your goal here?

cbontoiu commented 2 years ago

As I understood the average is given by 'e_all_energyDensity' itself with respect to a 3D mesh cell and I am happy with the concept. My problem was two fold:

Thank you

sbastrakov commented 2 years ago

As I mentioned above, the number of macroparticles per cell and shape are merely "sampling resolution" of the PIC method. So they should not significantly affect the physical picture you see. If they do, it hints at either a mistake at interpreting the data, or that the method with chosen parameters is not well-applicable to the problem.

PrometheusPi commented 2 years ago

Okay, I thought you were interested in the average energy per real electron resolved in space (as a scalar field).

PrometheusPi commented 2 years ago

Thus the increase in energy density you are seeing between y=550nm and 650nm comes from the energy gain in the laser field. The energy density between y=420nm and 500nm might have two reasons: either there are just more particles as the back of the plasma response to the laser and thus the sum over particles produces a higher total energy, or you have energetic particles (injection?) or both.

cbontoiu commented 2 years ago

Yes, so now I can safely state the energy gained by injected charge density cloud, and from the integration of this in volume I can infer the number of electrons captured and accelerated.

2021-09-17_10-30

PrometheusPi commented 2 years ago

An easier way to determine the "the number of electrons (captured)" would be to use the e_chargeDensity, do the same steps as above (1)-(3) (but now dividing by charge converts charge density to number density) and than you now exactly how many electrons are in each PIC cell.

PrometheusPi commented 2 years ago

@cbontoiu What did you change between your second-to-last and last plot? (the color scale changed by 7 orders of magnitude).

cbontoiu commented 2 years ago

@PrometheusPi this last image can be compared to the last from the batch of three shown previously. The latter is in MeV and the former in eV, hence the orders of magnitude difference. However, they are from different simulations.

cbontoiu commented 2 years ago

An easier way to determine the "the number of electrons (captured)" would be to use the e_chargeDensity, do the same steps as above (1)-(3) (but now dividing by charge converts charge density to number density) and than you now exactly how many electrons are in each PIC cell.

Yes, This is the usual way I do it. Not shown here.

cbontoiu commented 2 years ago

Hello I need to reopen this issue because I obtained some inconsistencies. For the two cases below I expected to obtain similar values for the average energy in mesh cell but the differences are large. The images below show the raw data as represented using the openPMD viewer and then after being multiplied by CELL_VOLUME/e

CASE 1 (low mesh) CELL_VOLUME = 1.004281e-30 [m^-3] image image image

CASE 2 (high mesh) CELL_VOLUME = 2.163462e-31 [m^-3] image image image

I can see that the histograms obtained using built in plug-in (EnergyHistogramData) image give similar results but I am confused how to obtain the colour density plot correctly. Any help will be greatly useful and appreciated. Thank you.

PrometheusPi commented 2 years ago

@cbontoiuCould you please explain, what the issue with the

the colour density plot

is in detail?

Why did you plot your energy density twice and with a 90° rotation?

cbontoiu commented 2 years ago

Hello @PrometheusPi and thank you for looking to this issue. The situation is like this. I have two 2D simulations for exactly the same setup: laser, target, ionization, particle per cell etc except for the mesh. More details are shown together with the colour plots provided. My problem is that I get different energy densities.

The two colour plots with the y-axis vertical are directly from the openPMD viewer just to check the raw e_all_energyDensity. Already, the growth of the peak value from ~7e16 [J/m^3] to ~ 8.5e16 [J/m^3] is of some concern.

The other two colour plots are my conversion to keV average energy after multiplying by CELL_VOLUME/e. However in both case the energy histogram looks similar, which makes me think that my colour plots should be smoothen. I just wanted to check with you how can mesh impact the energy density.

PrometheusPi commented 2 years ago

@cbontoiu If you say, they are exactly the same, I assume you mean from the point of few of physical parameter setup (including laser parameters, density, etc.) not physically-equivalent that they should behave the same due to symmetry reasons, is this correct?

Assuming the previous assumption is correct: How do you initialize your particle distribution? Do you use random our quite start in your initialization pipeline? One reason for different plasma dynamics despite the same global physics parameters might be that the initial random particle position when using random start positions results in different global physics effects. Another could simply come from the non-commutability of float additions which becomes a source for randomness even more with parallel processing.
Both cases would pinpoint to the fact that the system you are trying to simulate is intrinsically dynamic (at least at your resolution scale) and thus small changes might have a big effect.

Could you provide a difference plot to estimate the effect. Ideally by something like abs( sim_1_value - sim_2_value ) / abs(max(sim_1_value, sim_2_value) + epsilon) with epsilon being a minimal value of certainty you think should be achieved.

If your system is intrinsically dynamic, only statical predictions are possible. Bt before you assume that, you could test whether increasing the ppc and resolution results in more stable predictions.

cbontoiu commented 2 years ago

Hell @PrometheusPi and thank you so much for your reply. Your assumption is correct, I meant the setup. I don't use any symmetry, not even for x-direction though the laser is a plane wave and it is recommended to apply the x-symmetry. However I don't want to get the charge density cloud swapping sides; I mean what leave from max(x) shouldn't enter back from x = 0;

Here is how I initialize the Carbon atoms in the speciesInitialization.param, three times ionized and their released electrons

image

with fragments from the particle.param

image image image

I will do the comparison with a larger PPC. For now I can say that my injected charge moves non relativistically as shown below. I am not sure if this complicates the stability of the results.

image

cbontoiu commented 2 years ago

I have some results on this aspect. With the same mesh changing the PPC as 6,9,12 there is no difference on the histograms. image There is however a difference on the colour density plot. image I can understand that with larger PPC the figures become more accurate and converge, but I am not sure how to relate the colour plots to the histograms. Looking at the colours I can see that the injected bunch shows an energy spread in the range 0 - 100 keV at most, while there are some pixels going higher in energy. At the same time, the lower lines in the histograms (inside the bubble) show a sharp drop in the range ~75 - 100 keV and this is consistent.

The problem is that there are Counts for > 200 keV, that is energy values which don't appear in the colour plots. Their colour bar refers to the whole simulation domain and it is not restricted to the bubble. The histograms are binned in 1000 bins for a maximum energy of 1000 keV.

image

I guess I need to reduce the number of bins and thus get rid of the high energy tail.

PrometheusPi commented 2 years ago

Thanks for providing all the details. From a first glance your setup looks good. And the fact that you reproduce the(global and bubbly only) energy histogram independent of the ppc is great. That points to a stable simulation. (The same is true for the energy plot - they are within the statistical fluctuation very similar.

I guess that the difference between your energy histogram and your local energy plot originates from the fact that the code of the local energy uses the energy density times the cell volume (times the charge) and thus sums over all particles within a small volume (to be precises a few cells - depending on the particle shape used). Thus each pixel/voxel is a sum of all particle energies.

There might be now two reasons why your local energy is below the max value of the energy histogram (>1000MeV). 1.) Your weighting is below 1 thus your kinetic energy of a macro particle is the kinetic energy of a single electron times the weighting. Assuming you have a 1000MeV particle but its weighting is <<0.5, the kinetic energy is distributes in the energy density data is only 500MeV while the energy histogram takes it as a 1000MeV electron with weighting 0.5 (Here, I assumed there is only single particle that contributes to the histogram, where counts is at 0.5 < 1 electron - more likely there a multiple macro particles contributing to this energy bin and thus the weighting is much lower. ) . 2.) Depending on your particle shape, the maximum contribution you can get for a single macro particle is below E_kin * w since it is distributed over multiple cells. Assuming your particle is perfectly centered on a energy-density grid point, and your simulation is 2D, using a TSC shape: your maximum contribution would be 0.75^2=0.5625, for PCS it would be 0.22.

I hope this helped understand the discrepancy between both plots. I think both are valid, there interpretation with macro particles representing less than one real electron is however challenging/difficult.