Modular computation of diffusion matrices

[lfarv]

While the tracking in AT is modular (pass methods may be modified or added easily), the computation of the diffusion matrices used in ohmi_envelope is centralised and difficult to maintain.

The computation is done in findmpoleraddiffmatrix.c for all magnets, straight or curved. Wigglers are ignored. This C function combines the tracking and the computation of diffusion matrices and is based on Bnd/StrMPoleSymplectic4RadPass.

This PR is a proof of concept where the computation of diffusion matrices is included as optional in each *RadPass passmethod. This allows the computation to be consistent for each tracking variant, and has the following main advantages:

• it allows the computation for wigglers (to be done),
• it gives a better accuracy for large off-axis orbits with pass methods different from Bnd/StrMPoleSymplectic4RadPass (on-axis, all methods should be equivalent).

As an example, in this branch, the following passmethods include the computation of the diffusion matrix: BndMPoleSymplectic4RadPass, StrMPoleSymplectic4RadPass (with strictly identical results as before) and ExactMultipoleRadPass (with a slightly different result for large off-momentum). A new C function diffusion_matrix.c is added to replace findmpoleraddiffmatrix.c using the new pass methods while keeping a similar interface. ohmi_envelope is modified to call this new function, with a test mode allowing to switch easily between findmpoleraddiffmatrix.c and diffusion_matrix.c.

This approach will require all *RadPass methods to be modified, which is a significant task. As it is now, a non-modified method generates a zero contribution to the diffusion (but still propagates correctly the cumulated matrix): as if the element does not contribute to the emittance.

This PR is open to discuss the potential of this method and identify the difficulties before going further on. It follows the discussion in #715 and integrates the proposal of @carmignani (#715 (comment)) which makes the transition easier.

[lfarv]

The processing of the diffusion matrix requires two things:

1. For each integration step, one needs to propagate the current diffusion matrix. This requires the 6x6 transfer matrix of each of the elementary steps, computed around the closed orbit. Up to now, the transfer matrix is computed analytically for the following steps:

   • “classical” drift
   • “exact” drift
   • kick in straight geometry (same for "classical" and "exact")
   • kick in curved geometry (same for "classical" and "exact")
   • “classical” pole face angle (magnetic wedge + fringe field)
   • Y rotation

   Still missing for implementing the “exact” bends:

   • Curved drift
   • “exact” magnetic wedge
   • “exact” dipole fringe field

   Building the transfer matrix analytically becomes difficult for those: taking the derivatives of E. Forest’s formulae is tedious. It’s possible to build the matrix numerically in the C integrator, but this might be too slow at run time.

2. For the “radiating” steps, one needs to compute the increment to be added to the propagated diffusion matrix. This is already done for:

   • “classical” kick in straight geometry
   • “classical” kick in curved geometry
   • “exact” kick in straight geometry

   Still missing is:

   • “exact” kick in curved geometry

So there is a significant work to provide the “exact” bends. ON the other hand I think that all ingredients are there for the wiggler.

[lfarv]

As it is now, this is already better than the present computation: one could redirect the missing (non-upgraded) passmethods to BndMPoleSymplectic4RadPass, with a result identical to what we have now. But this implies that we keep a hard-coded list of upgraded methods. Otherwise, one can easily check that a passmethod does not modify the diffusion matrix, and in this case revert to BndMPoleSymplectic4RadPass, but this doubles the processing time…

The hard-coded list avoids that, but it has to be maintained… If you think of another way, let us know.

If this is solved, one could think of merging this as it is, or with the GWigSymplecticPass, and have a smooth transition while other passmethods are upgraded. Please advise!

[lfarv]

This study revealed three problems with the present implementation:

1. The rotation around the z-axis at entrance or exit is not taken into account. So a vertical bend is not correctly taken into account. This is easy to correct,
2. The entrance and exit face angles correctly propagate the diffusion matrix, but they do not generate diffusion. For off-axis closed orbit, I think they should: it adds or remove some magnetic field, so it should radiate, it's similar to quadrupole radiation. Implementing this is not obvious,
3. findmpoleraddiffmatrix.c still uses old physical constants. For comparison with the new method, I had to modify it to use the new ones. Results are then identical to machine precision.

I propose to open another pull request to solve points 1 and 3: it will allow vertical bends and make comparisons straightforward. What's your opinion on that?

[swhite2401]

Hello @lfarv, I think we discussed this and came to a conclusion on the strategy did you update the code to switch passmethods when necessary?
Maybe the implementation including vertical bends should be integrated here as well no?

[lfarv]

@ZeusMarti, @joanarenillas
I'm working in the integration of @joanarenillas's work in this branch, which merges tracking and computation of diffusion matrix in a single C file. And I have some suspicion on 2 points:

1. at wiggler entry and exit, FDW.c calls GWigGauge, while GWigSymplecticPass and GWigSymplecticRadPass don't. I don't know what this function does, but the tracking should be identical everywhere, so it should be identical in all 3 functions. What is your opinion?
2. I think there is a bug in the propagation of the cumulated diffusion matrix in FDW. For each slice, I see:

   wigglerM(&pWig, orbit_in, dl1, MKICK);   /* transfer matrix of the slice */
   wigglerB(&pWig, orbit_in, dl1, BKICK);   /* contribution of the slice to the diffusion */
   ATsandwichmmt(MKICK,BKICK);              /* propagation of the slice diffusion through the slice itself */
   ATaddmm(BKICK,BDIFF);                    /* accumulation of all slices in BDIFF */

   Here, BDIFF is the sum of all the slice contributions, but it is not propagated along the wiggler. The correct code should rather be:

   wigglerM(&pWig, orbit_in, dl1, MKICK);   /* transfer matrix of the slice */
   wigglerB(&pWig, orbit_in, dl1, BKICK);   /* contribution of the slice to the diffusion */
   ATsandwichmmt(MKICK,BDIFF);              /* propagation of CUMULATED matrix through the slice */
   ATaddmm(BKICK,BDIFF);                    /* Addition of the slice diffusion to the cumulated matrix*/

   This is how it is done in other integrators, and I think it's how it should be.
   The result is significantly different, and we are back to the question of the validation: how could we cross-check the result?

[ZeusMarti]

Hi @lfarv ,

I'll try to help though I'm not sure I understand all the details.

1. My guess is that GWigGauge redefines the momentum variables within the wiggler (see Ohmi's paper eq.35). We have been discussing a similar issue with @oscarxblanco recently, I do not understand it well but we have seen that the effect is quite quite small.
2. I can't help much here, it could be a typo, no idea, I should spend some time looking into it to say something well informed.

Regarding the tests I can think about some possibilities:

1. Lattice emitance calculation tests (FDW.c) against the results from integrals method.
2. GWigSymplecticPass: Single element tracking tests against a kickmap method. We did such tests for a single harmonic where the kick map is easy to calculate and the agreement was quite good but no idea how to interpret the discrepancies since the kick map is an approximation. We could compare it to a Runge Kutta, it should be easy since we did it also with radiation.
3. GWigSymplecticRadPass: Single element tracking with radiation tests. We used a Runge Kutta model, the problem is that it is very slow (half an hour in my low budget desktop PC) and since the effect is so faint, to have a good accuracy the field must be unrealistic above 30 T or so, but it is just for test.

Not sure if any of this sounds interesting to you, I should be able to help but our MAC is in too weeks so I could really start looking into after Nov 12...

[lfarv]

Hi all,

This branch is now ready. To summarise:

• the computation of diffusion matrices is now part of the *RadPass methods, so that we can get rid of the specialised functions findmpoleraddiffmatrix.c and FDW.c. This opens the possibility of adding new radiating elements by simple providing a new *RadPass method.
• up to now, the converted methods are BndMPoleSymplectic4RadPass, StrMPoleSymplectic4RadPass, GWigSymplecticRadPass, ExactMultipoleRadPass. Other RadPass methods are converted to BndMPoleSymplectic4RadPass, which is similar to what happens now, and gives a good approximation for bending magnets.
• Work has started for other exact methods, but require very long analytical developments…
• This branch also includes the contents of Access to the lattice energy in integrators #816, which allows energy ramping.

[lfarv]

The number of modified files is very large. Here is a summary of the modifications:

• BndMPoleSymplectic4RadPass, StrMPoleSymplectic4RadPass, GWigSymplecticRadPass, ExactMultipoleRadPass have been heavily modified to include the computation of diffusion matrices,
• All other *RadPass methods are modified to handle the ring energy, always taken from the lattice in python, and preferably from the lattice in Matlab. Some other methods also using the energy are modified in the same way,
• All passmethod are slightly modified in the Matlab-only part to allow additional input parameters in the MexFunction (necessary to provide the energy and particle to the MexFunction),
• python and Matlab function dealing with tracking or diffusion are adapted to match the new passmethods.

Concerning the wigglers, I keep the modification discussed above for the propagation of the diffusion matrix (differing from FDW.c). This has to be validated.

[oscarxblanco]

@ZeusMarti, @joanarenillas I'm working in the integration of @joanarenillas's work in this branch, which merges tracking and computation of diffusion matrix in a single C file. And I have some suspicion on 2 points:
1. at wiggler entry and exit, FDW.c calls GWigGauge, while GWigSymplecticPass and GWigSymplecticRadPass don't. I don't know what this function does, but the tracking should be identical everywhere, so it should be identical in all 3 functions. What is your opinion?

@lfarv I think you found the function I was looking for the last months but was unable to find because I was going through the diffussion matrices and not the tracking.

I agree this is an error. For me GWigSymplecticRadPass should contain GWigGauge because it will calculate the vector potential components $A_x$ and $A_y$ and add them to the kinetic momenta $\gamma_r m_odx/dt$, $\gamma_r m_o dy/dt$, as in Eq. (35) of Ohmi's article

('From the beam-envelope matrix to synchrotron-radiation integrals', Phys. Rev. E Vol. 49, Number 1, jan 1994, https://journals.aps.org/pre/abstract/10.1103/PhysRevE.49.751) in order to keep the canonical momentum $p_{x,y}$ vector in the second and fourth components of the 6D-variable x. From there at every tracking step the transverse components of the vector potential are substracted before tracking and added back after the tracking.

I think that GWigGauge should be used twice. First time at the entry point of the element before tracking and the second time at the exit point of the element after tracking. The transverse components of the vector potential are added only when inside the element.

Apart from this issue in the transverse plane, there might be discrepancy in the longitudinal plane due to the difference between $z$ in Ohmi's coordinate vector and $c\tau$ in AT as longitudinal coordinates. This difference do not affect directly the diffusion matrix $B$ in Eq (44) because the relevant columns are zero

but, could affect the constant in front of the matrix explicitely written in Eq~(45)

and also the tracking. Unfortunately, I cannot tell immediately what is correct because I was lost on the transverse plane issue that now appears to me clearer. Thank you for bringing GWigGauge up @lfarv .

[lfarv]

@oscarxblanco: thanks for your analysis.

As it is now, in this branch GWigGauge is called, so that should be correct.

While you are on it, I have another concern: when computing the radiation kick (GWigRadiationKicks in wigrad.c, line 117), the transverse magnetic field is used (longitudinal component ignored):

  /* B^2 in T^2 */
  B2 = (Bxy[0]*Bxy[0]) + (Bxy[1]*Bxy[1]);
 ...
  /* Beam rigidity in T*m */
  H = 1.0e9 * pWig->Po * __E0 / C0;

  /* 1/rho^2 */
  irho2 = B2/(H*H);
  ...
  /* Classical radiation loss */
    dDelta = -(pWig->srCoef)*dFactor*irho2*dl;

Normally, one should use the field component perpendicular to the trajectory, as it is done in other RadPass methods (ex. diff_str_kick.c, line 44):

  B2P = B2perp(ImSum, ReSum, x, xpr, y, ypr);
  ...
  /* Momentum loss */
  r6[4] -= rad_const * B2P * factor;

where B2perp make the vector product of B and a unit vector following the momentum.

This should not make a big difference, but maybe one should do it ?

[oscarxblanco]

@oscarxblanco: thanks for your analysis.

As it is now, in this branch GWigGauge is called, so that should be correct.

While you are on it, I have another concern: when computing the radiation kick (GWigRadiationKicks in wigrad.c, line 117), the transverse magnetic field is used (longitudinal component ignored):

  /* B^2 in T^2 */
  B2 = (Bxy[0]*Bxy[0]) + (Bxy[1]*Bxy[1]);
 ...
  /* Beam rigidity in T*m */
  H = 1.0e9 * pWig->Po * __E0 / C0;

  /* 1/rho^2 */
  irho2 = B2/(H*H);
  ...
  /* Classical radiation loss */
    dDelta = -(pWig->srCoef)*dFactor*irho2*dl;

Normally, one should use the field component perpendicular to the trajectory, as it is done in other RadPass methods (ex. diff_str_kick.c, line 44):

  B2P = B2perp(ImSum, ReSum, x, xpr, y, ypr);
  ...
  /* Momentum loss */
  r6[4] -= rad_const * B2P * factor;

where B2perp make the vector product of B and a unit vector following the momentum.

This should not make a big difference, but maybe one should do it ?

@lfarv right, I agree with you. One should calculate the B field perpendicular to the particle trajectory.

I cannot yet find it anywhere for the wiggler, this implementation assumes the particle is parallel to the longitudinal axis.

[ZeusMarti]

Hi @lfarv and @oscarxblanco ,

I'm not sure this is exactly what you meant but I remember Joan implemented something similar: see commit before merge:

at/atintegrators/gwigR.c

Line 607 in 57f9003

static double vectorialnB2(struct gwigR *pWig, double *X, double *Bxyz)

I remember removing this new formula because we could not make the Runge-Kutta precise enough to validate such an small effect (not even for an unrealistic 30T 1m ID). The same argument can be used to actually include it so, as you prefer!

Cheers,

[oscarxblanco]

Dear @ZeusMarti , do you remember what were the longitudinal harmonics you tested ? In AT matlab it is called B5 on the GWigSymplecticPass help message.

[oscarxblanco]

Dear @ZeusMarti and @lfarv , I have checked the function vectorialnB2 and I have the following comments:

1. The following line is an approximation to the particle energy, but it is not the total energy because you are ignoring the transverse momenta. It might be ok as it is
   

   at/atintegrators/gwigR.c

   Line 615 in 57f9003

   E=(pWig->E0)*(1+X[4]);

2. The following line is wrong. The correct expression for coef should have the nominal energy $E_0$ on the numerator, not the particle energy $E$.
   

   at/atintegrators/gwigR.c

   Line 618 in 57f9003

   coef=(clight*sqrt(E*E-XMC2*XMC2))/(XMC2*gamma);

   
   Here is the expression I did to confirm it :
   
3. The normalization of $B$, called b as described in the comments, is ambiguous. It is not clear if he intends to normalize with respect to the reference momentum $P_0$ or the particle momentum $P$.
   

   at/atintegrators/gwigR.c

   Lines 608 to 609 in 57f9003

   	/*Calculates SQR(|n x b|), where n is a unit vector in the direction of velocity and
   	b is the normalised magnetic field*/

   
   therefore, the following expression the total energy $E$ but naming it $P_0$ is hard to decipher unless you know exactly where and why the output of this function is used.
   

   at/atintegrators/gwigR.c

   Line 626 in 57f9003

   P0=sqrt(E*E-XMC2*XMC2)*1e09*q_e/clight;

[ZeusMarti]

Thanks @oscarxblanco ,

Sorry guys, I should have checked Joan's work more thoughtfully...

[oscarxblanco]

Dear @lfarv, how would you like to proceed ?
would you rewrite it ? use part of the code fixing it in some parts ? or leave it as it is considering that the error is small ?

I have just talked with @ZeusMarti , and @ZeusMarti and me still miss to answer your second question.
I am still going through the article and checking the code, and it is tipically a slow process while I understand it. Therefore, it will take me a while to give an answer.

@ZeusMarti, @joanarenillas I'm working in the integration of @joanarenillas's work in this branch, which merges tracking and computation of diffusion matrix in a single C file. And I have some suspicion on 2 points:
...
2. I think there is a bug in the propagation of the cumulated diffusion matrix in FDW. For each slice, I see:
c wigglerM(&pWig, orbit_in, dl1, MKICK); /* transfer matrix of the slice */ wigglerB(&pWig, orbit_in, dl1, BKICK); /* contribution of the slice to the diffusion */ ATsandwichmmt(MKICK,BKICK); /* propagation of the slice diffusion through the slice itself */ ATaddmm(BKICK,BDIFF); /* accumulation of all slices in BDIFF */
Here, BDIFF is the sum of all the slice contributions, but it is not propagated along the wiggler. The correct code should rather be:
c wigglerM(&pWig, orbit_in, dl1, MKICK); /* transfer matrix of the slice */ wigglerB(&pWig, orbit_in, dl1, BKICK); /* contribution of the slice to the diffusion */ ATsandwichmmt(MKICK,BDIFF); /* propagation of CUMULATED matrix through the slice */ ATaddmm(BKICK,BDIFF); /* Addition of the slice diffusion to the cumulated matrix*/
This is how it is done in other integrators, and I think it's how it should be.
The result is significantly different, and we are back to the question of the validation: how could we cross-check the result?

[lfarv]

@oscarxblanco, @ZeusMarti:

Concerning "question 2", it is not specific to wiggler, it's just about the proper way of cumulating and propagating diffusion matrices. I am still convinced that the modified version is correct, and I would keep for the moment.

For the rest: this version is exactly as merged by @joanarenillas in #759 (which does not include the vectorialnB2 function mentioned by @ZeusMarti). The only difference is that tracking and computation of diffusion are now merged in a single file, as for all the other *RadPass methods. The wiggler still misses some validation, but it was considered anyway better that what was there before. I would not like to block this PR for wiggler-specific problems, considering that the question of "B2perp" should not make a big difference. So I would propose to go on with this branch as it is, and tackle the wiggler questions in another PR (based on this one!).

For this one, coming back to the topic of "Modular computation of diffusion matrices", I would like to have feedback from @swhite2401, @simoneliuzzo, @carmignani and others. As long as the PR is not approved, improvements on GWigSymplecticRadPass are of course still welcome!

[lfarv]

I just created #860 to continue the work on wiggler radiation, and I would like to merge this one. So please @swhite2401, @carmignani, @simoneliuzzo, @oscarxblanco, give you comments or approval so that we can go on…