Improve Green's Function Performance

[DanShort12]

In GitLab by @NathanUKAEA on Jun 21, 2021, 13:29

In greens_psi, greens_Bx and greens_Bz the same factors, k, E(k) and K(k) are re-computed. I'd suggest placing this into one function that finds these once and uses them as required rather than re-calculating them multiple times for psi, Bx and Bz separately.

Also, the Scipy function uses interpolation for E and K, this may increase errors for values of k near 1, i.e. the field very close to the loop. Two methods to improve this may be to use the AGM calculation (Wikipedia has the simplest method and finds E and K simultaneously) and/or to use a different modulus, k_1 = (1-sqrt(1-k^2))/(1+sqrt(1-k^2)), this approaches 1 much slower for the same distance and so the AGM converges very quickly. Again Wikipedia has the general method for doing this. Maxwell himself seemed to have cracked this so there are plenty of other sources if required.

Lastly, the Green's functions are overloaded to deal with finite-sized wires using the d_xc and d_zc inputs, however, currently, they are only set to take into account the width only in z and not x. I'd suggest including both the d_xc and dz_c but for isolated coils, using Gaussian quadrature points and weights (evaluate at +/- 1/sqrt(3)*d_zc instead of +/- dz_c) so it requires the same 4 evaluations but gives 3rd order accuracy, rather than the trapezoid method which requires the same number of evaluations but gives 2nd order accuracy.

For fields far from the wire, the infinitely thin wire at the central point should be fine and the 3rd order corrections more than enough. However, for the meshed plasma loops, I would use the corner points (or other 'serendipity' nodes) instead of the Gaussian nodes as it will allow 2nd or higher quadrature with few extra function calculations but requires a slightly different method linking the plasma mesh to the green's function calculations instead of calculating them as separate loops i.e. cornerpoints will require 4NxNz but calculating them all together only requires slightly more than Nx*Nz.

More than happy to elaborate on this further if required.

[DanShort12]

In GitLab by @CoronelBuendia on Jun 21, 2021, 16:04

I'm afraid I only follow the first two paragraphs. If you have improvements to compute the elliptic integrals then I'm all for it. I had a look at it recently when I wanted to speed things up, and there are definitely approaches to calculating ellipe and ellipk simultaneously which would no doubt offer some gains.

I have given some thought in the past to just using an external interpolation or fit of ellipk(k) and ellipe(k) over 0 <= k <= 1 which would get rid of scipy altogether here, but I kind of didn't dare because it seemed like an unnecessary introduction of errors. In general though, we should confirm that the greens_.. are actually showing up in bottlenecks and how bad performance is before deciding upon resource allocation.

Also, I'd be very interested to see any improvements to the semi-analytic integral approach taken for a constant current density circular coil of rectangular cross-section.

[DanShort12]

In GitLab by @NathanUKAEA on Jun 21, 2021, 16:31

Yes, apologies the last two paragraphs are rambling walls of text that could be much better described by two labelled diagrams. Happy to add those shortly.

[DanShort12]

In GitLab by @NathanUKAEA on Jun 23, 2021, 11:43

Apologies meant to add these yesterday:

[DanShort12]

In GitLab by @CoronelBuendia on Jun 28, 2021, 11:41

I understand better now, thank you. What we do at present for "meshed"[sic] coils and the plasma is just nr * nz evaluations, so similar to your second graphic but without the external nodes. The d_xcand d_zcif statements I think are no longer directly used, but we should confirm that and remove them.

Would the logarithmic divergence still exist with Gauss-Legendre nodes?

[DanShort12]

In GitLab by @NathanUKAEA on Jun 28, 2021, 14:13

The corner-point or Gauss-Legendre calculations that take into account dx and dz, should greatly increase accuracy for little additional overhead, so I think they should be considered. Especially if the underlying Green's functions can be calculated efficiently.

The Gauss-Legendre nodes aren't exponentially closer to the current source than the corner points, so they shouldn't have any problems. From looking at the convergence of the AGM method for the elliptic integrals, it's only once you get closer than $10^{-3} \times$radius from the loop that accuracy is affected at all and $10^{-6} \times$ radius before more terms in the AGM are really required. For the polynomial approx method used in Scipy, I imagine it'll have similar accuracy.