Passing pitch angle requirements

[David Dickinson]

I noticed that at high theta resolution (ntheta=512) the velocity space integration weights errors were getting rather large for a typical ngauss. This got better as I increased ngauss/npassing. Below is a plot of the pitch angle weights errors for low and high npassing for a shaped case with ntheta=512. It can be seen that for the low npassing case the errors can be significant (~1%) and are particularly noticeable near the inboard side.

This error is arising from the passing particle weights. There are two steps to the passing particle weights calculation:

1. Get the Legendre weights, wx, by calling get_legendre_grids_from_cheb (1., 0., xx, wx)
2. Transform to the relevant domain to get the true weights, wl, by doing wl(ig,il) = wx(il)*2.0*sqrt((bmag(ig)/bmax) *((1.0/bmax-lambda_grid(il))/(1.0/bmag(ig)-lambda_grid(il))))

It is clear that wxis independent of theta so alone cannot describe the theta dependence in the error (indeed the error on wx is around 1e-16 for both npassing values). If we define the mapping function as map = 2.0*sqrt((bmag(ig)/bmax) *((1.0/bmax-lambda_grid(il))/(1.0/bmag(ig)-lambda_grid(il)))) we can consider plotting this for different theta values for the two npassing values.

Here is a comparison for theta~0

Here is a comparison for theta~- pi (but not exactly -pi)

The points lie exactly on top of the curve in both cases (as we would expect as it’s an analytic evaluation just on different grids) but we can note that for the inboard case we’re perhaps struggling to really resolve the structure of the mapping function.

It seems that our required npassing resolution is being controlled by our ability to resolve the mapping function at all theta grid values. This may not be particularly obvious so I wanted to flag this.

Note that in the limit of theta==pi the mapping function becomes easy to resolve (it is exactly 2.0 at all passing lambda).

‌

It’s not clear if there is anything we can do to avoid this within our existing approach.

[David Dickinson]

Here is a plot showing the maximum error vs npassing (~2*ngauss) for both values of new_trap_int. We can see that the two curves are ~identical until around npassing ~ 80. This indicates that the passing pitch angle weights dominate this error until this point. After this we see that the new_trap_int = F case hits a minimum error level of ~1e-7 suggesting a single floating point error limit (or that we need to increase the trapped pitch resolution further, however this is a case with ntheta=512), whilst the new_trap_int = T case continues to fall further and will eventually hit the double precision limit of ~1e-16.

‌

‌

The passing region is treated with a spectral approach so we would typically expect errors to drop rapidly with increasing numbers of grid points once we have sufficient points to resolve the appropriate structures. Here we are integrating a constant, which should be easy to resolve, however in effect I think we are integrating the mapping function times our real integrand (1) with the true spectral scheme so we are perhaps seeing that we need a lot of points to appropriately resolve our mapping function.

[David Dickinson]

new_trap_int = T actually levels off at 1e-13.

[David Dickinson]

Here is an expanded form with linear x-scale.

One can see that this does show exponential convergence, albeit at a relatively slow rate.

[David Dickinson]

Here is a comparison with Radau grids enabled.

‌

[David Dickinson]

The mapping function is map = 2.0*sqrt((bmag(ig)/bmax) *((1.0/bmax-lambda_grid(il))/(1.0/bmag(ig)-lambda_grid(il)))) . Suppose we wish to integrate this over the passing pitch angles, 0<lambda<1/bmax, at a particular theta grid point:

I = Int_0^{1/bmax} { map dlambda}

Let us define c = bmax/bmag and v = bmax lambda such that the upper limit of the integral is 1 and dlambda becomes dv / bmax

I = [2.0/sqrt(c)] Int_0^1 { sqrt([1-v]/[c-v]) dv /bmax}

The integral can be written as

Int_0^1{sqrt([1-v]/[c-v]) dv / bmax} = (1/bmax) [ sqrt(c) + (1-c) tanh^-1 (1/sqrt(c))]

so we find our final result:

I = (2/bmax) [1 + {(1-c)/sqrt(c)} tanh^1-(1/sqrt(c))]

We can then plot this as a function of theta and compare to the sum of the passing pitch angle weights (with npassing = 47) as shown below

Clearly these two curves do not agree in structure or magnitude so I’m likely to be doing something wrong here. It can be noted that summing over half of the passing pitch angle weights appears to give fairly good agreement with the analytic result over the central region, but there’s no obvious justification for this.

[David Dickinson]

The issue with the above lies with the fact that when we sum the pitch angle weights, this is equivalent to integrating over the mapping function evaluated on the underlying Legendre grid rather than on the lambda grid. To be more explicit, we define lambda = (1-xx^2)/bmax where xx are the Legendre grid points. We can write out mapping function in terms of xx rather than lambda:

map = 2.0*sqrt((bmag(ig)/bmax) *((1.0/bmax-lambda_grid(il))/(1.0/bmag(ig)-lambda_grid(il))))

map = (2.0/sqrt(c)) * sqrt{[(1/bmax)(1-(1-xx^2))] / [(1/bmax)(c-(1-xx^2))]}

Our integral is therefore

I = (2.0/sqrt(c)) Int(sqrt([xx^2]/[(c-1 + xx^2)]) dxx)_0^1

With the final result

I = (2.0/sqrt(c)) [sqrt(c) - sqrt(c-1)] = 2.0 [1 - sqrt(c-1)/sqrt(c)]

This results in the following much better comparison

‌

[David Dickinson]

We’re now able to directly calculate the error on the passing pitch angle integration by simply summing the passing pitch angle weights and comparing to the analytic result from the previous post. This gives us data like

‌

Increasing npassing both increases the rate at which error decays as we move from theta ~ pi to theta~0 and the peak error around theta~pi. Reducing ntheta / reducing the inboard local resolution may allow us to “step over” the region where the passing error is large.

[David Dickinson]

One approach to addressing this issue is to split the passing pitch angle into two separate regions. Each of which are treated spectrally in the same way the current pitch angle grid is. This would allow us to pack points in the challenging parts of the pitch angle space. Of course this will degrade the quality of the treatment elsewhere in the passing pitch angle region. To demonstrate this I have split the passing region such that rather than one Legendre region covering [1,0] we have two covering [1,0.05] and [0.05,0] respectively. I’ve then allocated different numbers of points to each of the two regions such that the total is still equal to npassing. The resulting errors are shown below for a case with npassing = 20.

We can note that splitting into two regions in this way has increased the error relative to the original case in blue for some regions of the theta grid. However, in the two cases where we have two regions we can see that the maximum error is reduced (roughly 7e-4, 3e-5, 8e-6 for the blue, orange and green curves respectively).

[David Dickinson]

It is possible to try to tune how we split this region into two (e.g. by trying to minimise the max or average error). Doing so with npassing=48 gives the following error behaviour (note the two split cases always split the number of points evenly between the two regions). We can see that both tuned cases have substantially lower maximum error, with the case tuned against the average error showing a peak error five orders of magnitude smaller than the original case.

The tuning isn’t always particularly successful as shown by the following npassing=10 case

The maximum error is lower in both of the tuned cases (0.05 for Original, 0.00073 for tuned-max and 0.0016 for tuned-average) but clearly the errors around theta~0 are much larger. This is likely due to the even split of the number of points. We could also consider trialling different distributions of points but this is likely to complicate the optimisation somewhat. Perhaps we just need a better way of choosing a default split rather than even (similar to how we set nesuper and nesub when the user specifies negrid). For example the below shows the case with npassing=10 but with three quarters of the points in the larger lambda region rather than just half. This marginally increases the maximum error in the tuned-max case (0.00078 vs 0.00073) and a bit more for the tuned-average (0.0022 vs 0.0016) but reduces the error substantially around the theta~0 region.

‌

[David Dickinson]

I can’t edit the above, but need to flag that there’s a typo. The maximum error on the original case should be 0.005 rather than 0.05.

[David Dickinson]

Here’s a repeat of the maximum total pitch angle error vs npassing as shown near the start of this issue. We include the original case and then values from the “Tuned-max” case for different point distributions (the value after max, refers to how many points are used in the largest of the pitch angle regions corresponding to the most passing particles).

‌

[David Dickinson]

We now reach the error floor by around npassing = 40 whilst we required around npassing = 90 for the original case.

[David Dickinson]

Here’s a plot of the mean pitch angle error vs npassing for two different points distributions and tuning against either the max or mean error:

This shows that we can get substantially lower mean error at moderate npassing, where we also get somewhat lower max error. At low npassing we find the mean error is comparable to the original case (sometimes a bit higher, sometimes a bit lower) but we do have lower max error.

Overall it looks like splitting into two regions may offer some potential to reduce our velocity space integration errors but this relies somewhat on tuning the position of the split and distribution of the points. We may also find that the error for integrating functions with more structure in the pitch angle may exhibit different error behaviour.