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]

radindex = radindex & ~Wig

is more efficient and easier to understand

[joanarenillas]

Thanks for the comment, I'll change this.

[lfarv]

Has it really been initiated in 2018 ? You could mention your contribution…

[joanarenillas]

I will modify it and add my contribution to all files in the next commit.

[lfarv]

@ZeusMarti, @joanarenillas:

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.

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…

[lfarv]

You should propagate the diffusion matrix through the rotation matrix. You are right that at input, propagating a zero matrix is useless, but I think here it should be done.

[joanarenillas]

Ok, I see there is a diff_yrot.c function in your branch that does the job. I will update that when adapting it to #742.

[lfarv]

Be careful, Yrot is a rotation around the Y axis. Here we are talking of a rotation around the Z axis. I think that the propagation consists just in applying the "ATsandwich" to the rotation matrix (R2).

Note that this is also missing in the present findmpoleraddifmatrix. I found that while working on the new branch.

[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]

T2 is not a 6x6 transfer matrix, it's a vector ! You cannot send it to ATsandwichmmt. Did you try and check the result ? It should crash…

A translation has no effect on beam matrices. There is nothing to do to BDIFF.

[joanarenillas]

That is right. It is now corrected.

[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.