ocelot-collab / ocelot

OCELOT is a multiphysics simulation toolkit designed for studying FEL and storage ring-based light sources.
GNU General Public License v3.0
84 stars 56 forks source link

Jitter studies / set fix RF parameters #142

Open jonasbjorklundsvensson opened 2 years ago

jonasbjorklundsvensson commented 2 years ago

Hi,

I want to use Ocelot for a jitter study, mainly of the longitudinal bunch parameters, but it's not clear to me how to implement this jitter in an appropriate way.

If I change for example the on-crest energy gain of a module, then the nominal particle parameters will be reset such that the effective path length is always the same. How can I set "fix" nominal parameters and then examine the variation of the beam parameters around this nominal value by changing RF settings?

I found the p_array.t value, but this does not seem to correspond to the absolute time from the simulation start (does not come even remotely close to p_array.s / c), and p_array.s is anyway constant.

Best regards, Jonas

sergey-tomin commented 2 years ago

HI,

If I understand you correctly we do not have yet a "nice" way to do this. Ideally, we need to specify the reference energy of any magnetic element. Then, when you track the beam which is not matched to this reference energy, it will see different bending angle or quad strength. So we do not have this reference energy for a magnetic element. However, you can use a workaround, if your energy jitter is small. You can shift energy of all particles in the ParticleArray, something like that: p_array.p()[:] -= delta_energy Maybe you need to implement your own Physics Process to automate it.

Cheers, Sergey

jonasbjorklundsvensson commented 2 years ago

Hi Sergey,

All right, that's too bad (it does seem like you understand what I mean more or less correctly!). I want to be able to jitter the amplitude and phase of the RF elements and observe changes in mean energy and chirp, and further down the machine the induced changes to mean energy, chirp(s), arrival time...

What you describe is something like what is implemented in Elegant, where one can choose if the nominal energy of the elements is automatically updated to the mean beam energy or not. Calculating the nominal energy and setting it programmatically should be much easier with Ocelot than with Elegant because of Python, so I need to think about this (not being very experienced with Ocelot or Python), but it seems like it should be relatively simple to do. I just need to understand the programming better :)

Another thing it would be quite useful for is for PWFA, where one might want to set the energy centered on the witness bunch rather than the driver bunch, in the case that their energy difference is significant (can be on the multiple-%-level).

For the jitter study itself, I would expect the actual jitter to be small, (1e-2 %-, 1e-2 degree-level and amplitude and phase, respectively (SRF cavities)), but to be able to match up with some theory I've been working on it would be nice to go to significantly larger values than such to observe "non-linear" shifts of various metrics.

Best regards, Jonas

jonasbjorklundsvensson commented 2 years ago

Hi again Sergey,

As I am thinking about this and working on getting more into the code, it strikes me that there might be a problem which I don't know how to solve: propagation of the bunch jitter downstream in the machine.

I'm considering trying to implement the jitter as a physics process somehow, as you suggested, but there are beam properties missing which I would not know how to get out of the code, either because they are not there or because I'm misundertanding something. For example, mean energy offset from an RF module upstream of a bunch compressor will translate into a timing offset, and therefore phase offset, in RF modules downstream of the compressor. When I look at for example the output tws object, I cannot find, for example, an absolute time of the reference particle from the simulation start. It's also not clear to me if the same information can be extracted from the s-coordinate - the lattice elements should of course have the same nominal s positions, but I wouldn't know how to handle this if the actual path length is longer or shorter than nominal. All the particles in a p_array seem to have the same s too.

I have seen the tau and tautau parameters in the tws object, which I guess are the mean and variance of tau of the bunch, respectively, but I'm not sure how to interpret the mean tau as it seems to start with random values (perhaps because of the random input bunch) and then propagates in a very peculiar way (can both increase and decrease in bunch compressors and only changes on the few-micrometer level).

I hope my problems comes across. If I can't somehow infer the time shift for bunches with different energy offsets, then I will not be able to simulate the jitter propagation at all. Best regards, Jonas

sergey-tomin commented 2 years ago

Hi Jonas,

you should not use twiss object for this study. ParticleArray is what you need and ParticleArray.tau() is what you are looking for. this is time coordinate for each macroparticle.

Cheers, Sergey.

jonasbjorklundsvensson commented 2 years ago

Hi Sergey,

Thanks for the input, I can try to think about a way to use tau instead.

From what I understand from the definitions of p_array (for example in tutorial N2), tau is the spatial coordinate in the co-moving frame tau = c*t - s/b0? I suppose b0 = v0/c in this definition. I'm not sure I exactly understand this definition. The one I'm used to seeing is something like z = b0*c*t - s, so would this imply that tau itself is also normalized and a factor of b0 needs to be included for it to be the actual spatial coordinate (small differences, I know, but still...).

So, if I were to implement some jitter here with respect to some "global" nominal particle with tau = 0, the "local" nominal particles of the jittering bunches will have tau != 0 ?

Best regards, Jonas

sergey-tomin commented 2 years ago

If I look to the definition I see: tau = ct - s/b0 s is independent variable which is the distance along the beam line (which, in turn, is the path length of the reference particle) and v0 is the velocity of the reference particle, t is the time at which a particle arrives at the position s along the beam line. So tau = ct - s/b0 = c (t - s/v0), s/v0 corresponds to a reference particle and we can call it let's say t0. so tau= cdt as it is written at the beginning.

jonasbjorklundsvensson commented 2 years ago

Hi again @sergey-tomin,

I finally started looking into this in a little bit more detail, and I find myself stuck. I can't figure out in which order things are applied during tracking and where certain parameters are calculated and set. Would appreciate some help/input.

The way I think about the implementation is the following:

I really don't understand the code well enough at the moment to implement these things, could you point me in the right direction?

To note here is that simply doing a "hack" and "artificially" shifting particle energies, etc. will probably not capture the relevant physical effects. If we have a long beam wrt the RF wavelength, for example, we need to accurately model the impact a shift of amplitude and/or phase has on the imposed LPS curvature of the bunch. Also, apart from the RF jitters, shifting the input particle timing with respect to the reference is also of interest. If the reference is just recalculated from the bunch then that will never work, and I can't think of a good workaround.

Best regards, Jonas

jonasbjorklundsvensson commented 2 years ago

Hi,

Just a friendly reminder about my recent post - I would very much appreciate some input on the subject!

Best regards, Jonas

sergey-tomin commented 1 year ago

hi,

I haven't modeled jitters, so I may not have caught the whole problem, but if you have jitter in phase in the cavity, you just put that jitter in the cavity phase. If you have fluctuation in the voltage amplitude of the cavity, it's a little more complicated if you also want to catch the betatron mismatch. I would calculate energy different between reference and actual amplitude and shift the energy of whole beam on that difference. So now your whole beam a bit off from the reference energy. There may be some effects that didn't come to my attention, but that's something you need to figure out, solve, and ideally commit back to ocelot so the other accelerator guys can reuse it :) Cheers, Sergey.

jonasbjorklundsvensson commented 1 year ago

Hi Sergey!

The way I understand how the code works, the beam energy/arrival time is pretty much always centered on the bunch. I guess one of my issues is that I can't find where this reference energy/time is set, so I can't see where or how to introduce a change. Both jitter in phase and amplitude of the cavities give a shift of the particle energies compared to the nominal phase/amplitude value, and after a section with R56 != 0, this offset in energy will induce an arrival time offset which then needs to be included in the following RF sections.

If I could find out where this is set, it is not difficult to calculate the nominal energy based on the settings of the cavities - worst case I could first track a single reference particle through the lattice and use its properties as a reference. I could then apply an offset in amplitude and/or phase to the RF, but this effect is in principle nonlinear for large phase shifts (which, admittedly, are not very realistic to assume) so it would most effectively be added directly to the cavity. Perhaps this is doable, and it's just my programming skills that are lacking (which they are).

I am mostly concerned with the longitudinal (time-energy) effects at the moment, but as you say, an energy jitter is naturally accompanied by a betatron mismatch as well (which I would probably also want to investigate further down the line). To put it in slightly different terms, I want to be able to accurately simulate synchrotron oscillations where the bunch as a whole is basically performing the oscillations.

The above are only first-order effects - going through for example a compressor off-energy will naturally cause higher orders to play a larger role as well. This is something that my analytical approach is nowhere near covering (yet?).

Anyway, to summarize, if you could help me out with where to look in the code, and maybe give me some ideas as to what to try and fiddle with, it would help me tremendously. If I eventually figure out a solution that is not extremely janky, I would definitely be happy to commit it.

Best regards, Jonas

sergey-tomin commented 1 year ago

hi,

I introduced in ocelot what can help with your problem a7aaa15c786a8ab030f57562b4f47e2d30fba102

example how to use it see here:

https://nbviewer.org/github/ocelot-collab/ocelot/blob/dev/demos/ipython_tutorials/small_useful_features.ipynb#jitter

jonasbjorklundsvensson commented 1 year ago

Hi Sergey,

Thank you very much for the meeting and for doing this! It will be immensely useful and I look forward to playing around with it :)

I checked out both the code and the jupyter notebook and I think there might be a small correction that needs to be made. I don't think it's quite correct to set p_array.E = self.Eref at the second-to-last line in the LatticeEnergyProfile physics process. The energy offset should definitely propagate through the entire machine downstream of the source of the offset - only the relative offset (parray.p()) should be reset according to the reference energy as this is what the magnets will "see". I need to think more about the calculation of p_new, I currently don't understand it perfectly.

Further, this approach would probably be easy to extend to a function for shift in reference time as well, right? I forgot during our meeting that timing jitter is of course also of interest, for example estimating the impact of timing jitter from the gun laser which then also propagates downstream from the gun and can therefore not be simulated simply by an energy/chirp shift in the input parameters. Perhaps this is easy enough to do in the input parameters and doesn't need its own function, I'll check it later.

Best regards, Jonas

sergey-tomin commented 1 year ago

Hi Jonas,

I don't think it's quite correct to set p_array.E = self.Eref ...

I think as it is now it is correct. this procedure should not change the mean energy of the beam. It should just shift p_array.p() relative to new p_array.E

jonasbjorklundsvensson commented 1 year ago

Hi Sergey,

Maybe I need to think more about it, but from the jupyter notebook it also looks wrong. Below the last figures where the beam energy is printed, it has exactly the reference energies, which is not what I would expect from a jittering RF. There is clearly a shift in arrival time, which one would expect from a shift in energy w.r.t the reference, but the relative energy looks identical to the case without jitter (but is hard to tell by eye).

I think the absolute energy should be the same as the jittering second case, because obviously the RF parameters are a bit shifted, but the relative energy offset, delta, should be shifted to be centered around the reference energy. I think this will happen if p_array.E = self.Eref is removed or commented out - this line just sets the beam energy to the reference energy if I understand things correctly.

sergey-tomin commented 1 year ago

Energy coordinate in the OCELOT is p=(E - E0)/p0c if you calculate for each particle energy E before this procedure and after this procedure it should be the same. Particle should not gain or lose energy. But how you can achieve this if you shift "p" but E0 will be the same? It is obvious that E will be different in that case.

That is what I meant in my previous message but instead of calculating E for each particle one can calculate np.mean(E).

jonasbjorklundsvensson commented 1 year ago

Hi,

All right, I tested it out a bit now and it seems to just be me being confused about the different quantities in Ocelot. When I take the mean energy as Emean = np.average( p_array_out.p()*p_array_out.p0c + p_array_out.E ), there is only a very small difference between having the energy profile on or off - that small difference (something like 10-20 keV) I think I can safely attribute to randomness because I generate a new bunch with quite few particles each time. For reference, I'm testing with phase shifts of 1 degree, which for my lattice and parameters gives a shift in mean energy of ~0.6 MeV in the mean energy.

I also did a simple double-check of the expected shift in time/space (to first order), which seems to make sense. Good to know.

Could the mean energies wiht/without energy profile be added to the python notebook, to avoid someone else being confused by the same thing as me?

I haven't done more than very basic testing, but I would say that the issue can be closed. Thanks for all the help and patience!

Best regards, Jonas