Bug Fix & Update: Wiggler diffusion matrix

[joanarenillas]

Dear all,

This is an ongoing Work In Progress, which aims to correct some bugs found in GWigSymplecticPass.c and GWigSymplecticRadPass.c, as well as updating gwigR.c and ohmienvelope to include the computation of wiggler diffusion matrices. To do so, a new function (FDW.c) is currently being developed, following the work of @mashl.

More on current and past issues on wiggler passmethods can be found in #715 and #749.

No need to review yet.

[joanarenillas]

Dear @swhite2401,

I am currently cleaning up the code in FDW.c, which implements the computation of diffusion matrices for wigglers. I am having a little bit of trouble understanding how to code the Matlab mex function and the corresponding Python version. My guess is that it should look similar to the code in findmpoleraddiffmatrix.c, but I am not entirely sure whether I can just use the same lines of code or a more detailed makeover is necessary.

Could I get some advice with this?

Regards,

Joan

[lfarv]

Hello @joanarenillas,

Sorry for my late intervention, but I was away until now and just had a look to your proposal.

There is already a git branch dedicated to the computation of diffusion matrices, called diffusion_matrices, with a corresponding pull request: #742. In its discussion, you will find the motivation and the proposed solution.

In short, the idea is to have the computation of diffusion matrices as modular as the tracking itself, so that introducing a new kind of element, including its contribution to equilibrium emittances, can be made by simply dropping a new C file in the atintegrators directory,

Since the computation of diffusion matrices involves the tracking of the closed orbit through the element, a simple way is to add the diffusion as an optional computation in the integrator C file itself.

For backward compatibility, we cannot change the signature of the integrator, so to indicate the need for the computation of the diffusion matrix, we added a new variable in the struct parameters structure, called bdiff. bdiff is a (double *) pointer. If bdiff == NULL, it means that no computation of diffusion matrix is needed (simple tracking). Otherwise, bdiff points to an array of 36 doubles initialised with zeros. The integrator is expected to fill this array with the diffusion matrix of the element.

Existing integrators ignore this new variable, so are automatically considered as not contributing to diffusion. In the diffusion_matrices branch, most of the *RadPass have been modified to include the computation of diffusion matrices (still missing the so-called “Exact” integrators).

I think you should go directly to this new scheme by basing your working branch on diffusion_matrices instead of master. Then your modifications will concern only a single file: GWigSymplecticRadPass.c.

For an example, look at StrMPoleSymplectic4RadPass,c in the diffusion_matrices branch: it's drift-kick sequence in cartesian geometry, like the wiggler.

From what I can see, your new FDW.c file is a good starting point for this upgraded GWigSymplecticRadPass.c. I does both tracking and computation of the diffusion matrix. You just need to make the computation of diffusion optional by enclosing it in something like:

if (bdiff) {
    wigglerM(...);
    wigglerB(...);
    ATsandwichmmt(...);
    ATaddmm(...);
}

and to add a loop on particles for the "normal" tracking (of course, the computation of diffusion matrices is requested only with a single particle: the closed orbit).

[ZeusMarti]

Hi @joanarenillas, @lfarv and @swhite2401,

Besides Joan contributing to the diffusion_matrices branch, which I think it is a good idea, the present MR solves the #749 bug which I think is already a good reason to merge it, after your supervision of course.

Best,

[lfarv]

@ZeusMarti, @joanarenillas:

Besides Joan contributing to the diffusion_matrices branch, which I think it is a good idea, the present MR solves the https://github.com/atcollab/at/issues/749 bug which I think is already a good reason to merge it, after your supervision of course.

Ok, no objection! I have a few comments:

• it would be nice to add a small description in the documentation of passmethods (at/docs/common/sa_passmethods.md). If you don't want to dive into Markdown, send me a few lines, I can do it…
• Is there any reference to the use of Hessian matrices in this kind of computation?
• It looks difficult to me, but did you think of way to validate the results of this new feature?

Otherwise, for me we could merge this branch.

Also, one should think soon of converting this to the new architecture for computation of diffusion matrices: #742. This branch was designed to integrate easily wigglers but also any future kind of element. Though apparently it did raise much interest since it was opened…

[joanarenillas]

Dear @lfarv,

Regarding your comments:

• I will add a few lines to the documentation of passmethods. If I have any trouble doing so I will let you know.
• The computation of the Hessian matrices was done by @mashl. I am not sure whether he used any reference for the computation or he did it all himself. I will look into this and let you know if I find something.
• Together with @ZeusMarti, we were thinking of doing an slicing of the wiggler into many multipoles to get a very precise model for it, and then compare the emittance contributions of this model with that of FDW.c. We are still not completely sure if this is the right approach.

After this branch is merged, we are planning on adapting the computation of wiggler diffusion matrices of this branch to #742, so that wigglers are easily integrated.

Regards,

Joan

[lfarv]

I have no more comment. I am not able to appreciate the correctness of the result, I'll trust you . If @ZeusMarti and you are confident that everything is now correct, I can approve and merge, but if you have further validation to do, let me know.

[joanarenillas]

Dear @lfarv,

We are currently working on a benchmark to test our results. We hope to complete it soon. Once we have it, we will let you know the results so you can approve and merge.

Regards.

[joanarenillas]

Dear @lfarv,

We are having difficulties in completing the benchmark for the new set of integrators, since the particle dynamics are more complex than we anticipated. We have tried solving the Lorentz equation including energy loss (gamma factor is not constant) and taking into account the Abraham-Lorentz force. Nevertheless, the results we obtain are not 100% consistent.

Regarding the references on the Hessian matrices, we have not managed to get any from @mashl, who initially did the computation. I have however reviewed and summarised his calculation in a report which I am happy to share with you if its is of your interest. I have also included some documentation on the wiggler passmethods as you suggested.

Although the benchmark is not yet finished, I think it is still worth merging this branch, since it fixes many reported bugs with these passmethods, and includes the diffusion matrix computation for the wiggler. Moreover, the tracking in the new passmethods do not introduce significant deviations from the current version (they are smaller than 10^(-8) m).

My internship period at ALBA will finish soon and @ZeusMarti will take over from now.

Kind regards,

Joan

[swhite2401]

I agree with the merge, this is a major improvement. However, before that it would be nice to add simple unit tests. Is your report publicly available on the web? If not adding it, or a summary of it to the documentation pages would be very nice

[joanarenillas]

Hello @swhite2401,

The report will be available online in the following weeks after it is reviewed by the ALBA accelerator's team.

[ZeusMarti]

Hi @swhite2401, by "unit tests" you mean adding an additional test definition in at/pyat/test? What kind of check do you have in mind? a tracking result through the element with and without radiation could be enough?