Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Averaged over an axis in a 2D simulation #5460

Closed
Tissot11 opened this issue Nov 14, 2024 · 10 comments
Closed

Averaged over an axis in a 2D simulation #5460

Tissot11 opened this issue Nov 14, 2024 · 10 comments
Labels
enhancement New feature or request

Comments

@Tissot11
Copy link

Is it possible to write in the diagnostics data averaged over one direction e.g. y for a 2D simulation?

Also I wanted to save particle momentum p, so I used (#5131 input deck)

#diagCombined.particle_fields_to_plot = u
#diagCombined.particle_fields.u(x,y,z,ux,uy,uz) = sqrt( ux*ux + uy*uy + uz*uz)

But doing this gives u vs x and y. Is it possible to save p vs x for a 2D simulation with p being averaged over y in the diagnostic?

@Tissot11 Tissot11 added the enhancement New feature or request label Nov 14, 2024
@n01r
Copy link
Member

n01r commented Nov 14, 2024

Hi @Tissot11,

For FullDiagnostics field output, we do not have a direct option yet to average over a defined region, I believe.
However, there are ways to hack this.

You could use a combination of the options:

  • <diag_name>.coarsening_ratio = <Nx Ny Nz> which will average over every Nx, Ny, Nz cells
  • <diag_name>.diag_lo = <lo_x lo_y lo_z> which sets a lower limit to the box for field output in simulation coordinates (in units of meters)
  • <diag_name>.diag_hi = <hi_x hi_y hi_z> which sets the upper limits accordingly

You could play around and set the number of cells to average over to exactly the number of cells in the field output you are getting.

However, what you are describing with the momenta sounds to me like you would want to output phase spaces of particles. That is easily possible with the ParticleHistogram2D reduced diagnostics.
Check out this example for reference: https://warpx.readthedocs.io/en/latest/usage/examples/laser_ion/README.html#run
You need to switch over to the Executable: Input File tab because our PICMI input does not yet implement this feature, as far as I know. But the input deck you shared uses the same format.

PhaseSpaceElectrons.type                                 = ParticleHistogram2D
PhaseSpaceElectrons.intervals                            = 100
PhaseSpaceElectrons.species                              = electrons
PhaseSpaceElectrons.bin_number_abs                       = 1000
PhaseSpaceElectrons.bin_number_ord                       = 1000
PhaseSpaceElectrons.bin_min_abs                          = -5.e-6
PhaseSpaceElectrons.bin_max_abs                          = 25.e-6
PhaseSpaceElectrons.bin_min_ord                          = -197
PhaseSpaceElectrons.bin_max_ord                          = 197
PhaseSpaceElectrons.histogram_function_abs(t,x,y,z,ux,uy,uz,w) = "z"
PhaseSpaceElectrons.histogram_function_ord(t,x,y,z,ux,uy,uz,w) = "uz"
PhaseSpaceElectrons.value_function(t,x,y,z,ux,uy,uz,w) = "w"
PhaseSpaceElectrons.filter_function(t,x,y,z,ux,uy,uz,w) = "sqrt(x*x+y*y) < 1e-6"

You can choose axes for the 2D histograms freely (e.g., ux and x or total momentum u=sqrt( ux*ux + uy*uy + uz*uz) and x) how many bins and with which limits (in units of the quantity you chose, which for u means p = gamma * beta * m * c, and u = beta * gamma)

Does that answer your question?

@Tissot11
Copy link
Author

Hi @n01r

Thanks for the fast reply! I am already saving phase space using ParticleHistogram2D. This p vs x, I need for the spectra of particles. I need it with x, because I need to focus on a given region. I can also do it in the post-processing but it is always faster in the code itself.

I had overlooked diag_lo and diag_hi options. I had instead used the following

diagSpectra.particle_fields_to_plot = u
diagSpectra.particle_fields.u(x,y,z,ux,uy,uz) = sqrt( ux*ux + uy*uy + uz*uz)
diagSpectra.particle_fields.u.filter(x,y,z,ux,uy,uz) = (x>=Lx/2) * (x<0.75*Lx)

This should also produce the same output as with diag_lo and diag_hi options? Which one is more preferable?

Could you please look at other questions of mine and answer if possible. I am new to WarpX and although the documentation is very good, there are always some subtle points that would be better to clarify for me first otherwise I ended up misunderstanding some features.

@n01r
Copy link
Member

n01r commented Nov 14, 2024

The options diag_lo and diag_hi will actually reduce the size of the output.
The .filter(x,y..) option will give you the full-sized simulation domain but only outputs data that satisfies the filter.

I have not used the particle_fields option together with diag_lo and diag_hi but would assume that they work together.
I would suggest to try this in a minimal example.

I may not fully understand yet why you should not be able to add another ParticleHistogram2D for your p vs x.
Can you explain?

Why not do the following?

pvsx.type                                 = ParticleHistogram2D
pvsx.intervals                            = 100
pvsx.species                              = <your_species>
pvsx.bin_number_abs                       = <number of desired bins in x>
pvsx.bin_number_ord                       = <number of desired bins in p>
pvsx.bin_min_abs                          = Lx/2
pvsx.bin_max_abs                          = 0.75*Lx
pvsx.bin_min_ord                          = <expected lowest limit of p>
pvsx.bin_max_ord                          = <expected highest limit of p>
pvsx.histogram_function_abs(t,x,y,z,ux,uy,uz,w) = "x"
pvsx.histogram_function_ord(t,x,y,z,ux,uy,uz,w) = "sqrt( ux*ux + uy*uy + uz*uz)"
pvsx.value_function(t,x,y,z,ux,uy,uz,w) = "w"
pvsx.filter_function(t,x,y,z,ux,uy,uz,w) = <some other filter you may want to add, like filtering in y>

I can have a look at the other questions but would formulate an answer if I feel like I can be immediately helpful. 😅

@Tissot11
Copy link
Author

I did actually this as you suggested. Since this diagnostics is reduced, I assumed, it is better suited for checking quantities rather than analysing the physics because data in these files is rather small. I am yet to visualise these results though for these reduced diagnostics. I'll let you tomorrow.

@n01r
Copy link
Member

n01r commented Nov 14, 2024

Dunno if I got that right that there may be a misunderstanding, but to clarify: the "reduced" in ReducedDiagnostics does not mean reduced in quality or unfit for detailed analysis. 😁
It rather means that if you already know what kind of analysis you want to do that requires mathematical reduction operations (summation, averaging, ..., etc.), then you can set that up in situ.

The particle fields that you are using are similar to what you can set up with ParticleHistogram2D, only that the <diag_name>.particle_fields_to_plot will always output the data reduced (averaged or summed) over all particles in a simulation cell whereas with the histogram you choose your desired number of bins.

@n01r
Copy link
Member

n01r commented Nov 14, 2024

@Tissot11, feel free to chime in on #5462 in case it's not clear enough, yet. 🙂

@Tissot11
Copy link
Author

Ok. I guess because of the term ReducedDiagnostics and txt file format, I misunderstood the utility of these diagnostics. In fact, I am more interested in these ReducedDiagnostics. I just noticed that it writes data both in txt and h5 files. Though for .h5 format, it creates files for each interval. Can one create a single file as in Full diagnostic using <diagName>.openpmd_encoding = g?. I need to worry about the number of files too for all my simulations.

In fact, now I realized why the file sizes are bigger using Full diagnostic. Essentially it is a checkpoint but with data reductions using filters, coarsening ratio etc. So in the input deck in #5131, I was essentially having a checkpoint minus some data reductions that's why file sizes were really huge as I raised in #5458. I guess I need in Full diagnostic only fields_to_plot and particle tracking using plot_filter_function. For rest, I could use ReducedDiagnostics. Thanks for this discussion and now I understood these subtle points.

Another clarification needed: is <diagName>.particle_fields a data reduction on <diagName>.<species>.variables (for the same variables) in Full diagnostic?

@n01r
Copy link
Member

n01r commented Nov 15, 2024

Can one create a single file as in Full diagnostic using <diagName>.openpmd_encoding = g?

I just checked. Unfortunately, that has not been implemented so far, but I think it is a good change to make.

In fact, now I realized why the file sizes are bigger using Full diagnostic. Essentially it is a checkpoint but with data reductions ...

Yea, in a way. FullDiagnostics give you the raw particle and field data. In addition to that, over time, we implemented filters, reductions, etc. You have some flexibility in creating many diagnostics of the type Full to express which quantities you want at which output frequency, while deselecting the rest. But it is also true that the most memory you save if you can shift some of your analysis to in-situ calculations so that you can use the ReducedDiagnostics.

We are working on gradually making the diagnostics more configurable and comfortable for the user.
I am glad we could clarify your questions. 🙂
If you have new ideas for useful diagnostics, feel free to open issues for them. While we may not be able to immediately allocate time to implement them, another contributor might, or the idea might gather more community interest at which point they will be implemented more quickly.

Another clarification needed: is .particle_fields a data reduction on ..variables (for the same variables) in Full diagnostic?

You can see particle_fields as a spatial histogram with the binning fixed to one value per cell (as opposed to ParticleHistogram2D where you choose the binning). The dimensionality will be based on your simulation dimension and other factors like coarsening_ratio and diag_lo/hi apply. The value inside the bins can be an average or a sum of an arbitrary function based on the input values x, y, z, ux, uy, uz.
It is on our radar to make it possible that the electromagnetic fields interpolated to the particle positions can also be input parameters but so far that is not yet implemented.

@Tissot11
Copy link
Author

Thanks for this clarification and now I have answers to my confusions!

@n01r
Copy link
Member

n01r commented Nov 15, 2024

Cool, thanks for your questions!
Feel free to close the issue if it has been resolved. :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement New feature or request
Projects
None yet
Development

No branches or pull requests

2 participants