Unequally spaced element boundaries for grid optimisations

[mrhardman]

This pull request introduces flexibility into the specification of the element boundaries, naturally allowing packing of points near the wall boundaries for the purpose of capturing the square-root shape of the electric potential near the sheath entrance. The code modifications are minor -- in the specification of the grids, weights, and derivative routines there were variables called "scale_factor" and "shift". Using Gauss-Chebyshev-Lobatto points in the spectral element scheme, we can write these scale and shift factors in terms of differences and averages of the element boundary points. All that is required is to specify the desired locations of the element boundaries at the initial time of the simulation. Previously, element boundaries were assumed to be uniformly spaced.

Access the new functionality by specifying

`z_element_spacing_option = "sqrt" '

The default value is

`z_element_spacing_option = "uniform" '.

To see how this optimisation performs consider the output from the file "examples/wall-bc/wall-bc_cheb.toml", using the new element spacing option. Note that the electric field and potential are both very smooth near the sheath entrance. Here z_ngrid = 9 and z_nelement = 32.

wall-bc_cheb_phi(r0,z)_vs_z.pdf
wall-bc_cheb_Ez(r0,z)_vs_z.pdf

We can obtain similar forms for the potential and electric field even down to low resolutions. For z_ngrid = 5 and z_nelement = 8 we can obtain the following plots, keeping all other input parameters fixed. Whilst the curves are not smooth, the resolution is more optimally used so that no large fluctuations are observed where the electric field diverges.

wall-bc_cheb_lowres_Ez(r0,z)_vs_z.pdf
wall-bc_cheb_lowres_phi(r0,z)_vs_z.pdf

Closes #122

[johnomotani]

Looks good to me! Thanks @mrhardman!

There are a few docstrings/comments to tidy up, and I think we can add a reasonable number of tests without too much work.

[mrhardman]

@johnomotani I have just pushed a commit which addresses your comments to the main code. I still have to address your suggestions regarding extending the tests. This appears to be slightly involved as the "sqrt" option appears not to pass the tests to the same tolerances as the "uniform" option. This should be expected in cases where the test function has variation in the centre of the test domain.

[johnomotani]

which tests? I think the derivative tests should be not a problem, as they're differentiating polynomials that should be exact up to machine precision. In the wall-bc test, it may be the opposite problem, that with the sqrt grid the simulation is better resolved than in the benchmark case - anyway it should be fine to just give it however big a tolerance it needs to pass, as there are two parts to the test: one is the comparison to the benchmark results (which might have a fairly big difference); the other, more important I think, is comparing to saved results from the exact same test, which is a regression test to make sure we do not accidentally break the feature in future.

[johnomotani]

About that test failure: the integration code looks right, but I'm wondering if there could be an inconsistency in how element boundaries are treated in derivatives and integrals.

It looks like (from the fact that the test passed on uniform grids, and I guess would pass with one element) on a single element we would have integral(derivative(f)) = f[end] - f[start]. Now what if we have two elements? The integration weights are calculated independently for each element and summed at the element boundary, so if the elements have different widths the total integration weight for the element boundary has unequal contributions from the elements on either side. However, in the derivative routine, the derivative at the element boundary is taken to be the mean of the derivatives in the elements on the left and right. To keep the property that numerically the integral of the derivative is exactly the change in f, I think the derivative at the element boundary would have to be some sort of weighted sum of the derivatives from either element, to match up with the unequal contributions to the integration weights.

The proposal in the previous paragraph would mean giving more (?) weight to the derivative from the wider element, so it's not obvious to me that this is what we want to do - presumably (at least for a non-special function whose steepness is not related to the widths of the elements) the smaller element would give a more accurate estimate of the derivative at the element boundary, so you might think we'd want to give the smaller element a larger weight in the derivative at the element boundary...

[mrhardman]

If this is a matter of precision, and we can pass the tests with tolerances like 1.0e-8 or so, I would be in favour of marking this as an issue to be returned to, and reducing the tolerance in the specific tests using the "sqrt" option. I am not sure that spectral elements naturally satisfy the fundamental theorem of calculus in the general case. To the best of my understanding, no part of our code relies on this feature. To the best of my understanding, as we upwind the distribution function anyway, the fundamental theorem of calculus is broken in the kinetic equation, which means that density is not conserved anyway.

[johnomotani]

Happy to reduce the tolerance for the sqrt spacing tests, and leave anything further as a future todo.

One motivation for 'fixing' this at some point: although at the moment we're hardly using the 'central' derivatives anywhere (mostly using the upwind ones wherever possible), we could use the 'fundamental theorem of calculus' property to conserve particles in moment kinetic modes. If we used a centred element-boundary in the continuity equation (I think @mabarnes did originally and I changed it to an upwinding version when trying to make things more stable) we could conserve particles because we could have numerically that

$$\int_{-L/2}^{L/2}dz \frac{\partial(n u_\parallel)}{\partial z} = (nu_\parallel)|_{L/2} - (nu_\parallel)|_{-L/2}$$

As a bonus, the proposal for doing this would give a kind of partial upwinding - if element sizes get smaller towards the target, then at each element boundary we would be giving more weight to the derivative from the larger (upwind) element.

[mrhardman]

@johnomotani Your comments should now be addressed, with at least some testing support for the "sqrt" `element_spacing_option'. Could you take a quick look and comment whether or not we can merge this into master?

[johnomotani]

Looks good to me. Thanks @mrhardman! We can merge once the tests complete.