EZspline not giving correct derivatives

[zhucaoxiang]

This is a reminder that @aaroncbader has reported a bug that EZspline doesn't give correct derivatives. More details and reproduction will be provided by @aaroncbader later.

[mdjcole]

Did anything come of this?

[zhucaoxiang]

@mdjcole So far not yet. If you are also using the function, be careful. @aaroncbader stated the bug can be detected by comparing the derivatives with finite difference.

[mdjcole]

EZspline is used for constructing XGC-S equilibria. It would be very useful for me if @aaroncbader could provide his test case demonstrating the bug.

[lazersos]

I'd like @aaroncbader to upload an example of what he detected. Also, I'm a bit confused how one can detect an error in the spline by comparing with finite differences (which should have a larger error). I did some tests with 2D cos function and couldn't find some issues but these were a function of grid resolution.

[lazersos]

Here are some files which exemplify the issue. https://github.com/lazersos/SPLINE_TEST

[lazersos]

Just a note. I think the issue arrises from using the derivatives of the quantities. For the most part in STELLOPT and LIBSTELL I try to avoid using the derivatives of the spline quantities. Since VMEC has a Fourier representation for nearly all quantities the derivatives can be calculated analytically in the poloidal and toroidal directions. Which is what's been done. However, we do use the gradients when calculating radial derivatives. But from what I can remember these are only used in the coordinate inversion routine (not for inversion of the coordinates but for the gradient descent method to calculate the s value for a given R and Z).

Also this is a bigger deal than just STELLOPT. TRANSP uses these routines. In fact they're a direct import of the PSPLINE library found in NTCC

[aaroncbader]

Yes, the issue is specifically that the derivatives are incorrect. I'll try to set up a case demonstrating it, I just have had too many other obligations crowding out working on setting this up. I'll see what I can get done today.

[lazersos]

@aaroncbader Take a look at the repo I shared. It shows the problem without needing to do the finite difference calculation.

[aaroncbader]

Ok, i see. You are able to demonstrate that the u and v derivatives are in error. The radial derivatives though are also a problem, and can't be calculated analytically. The gamma_c metric often requires these radial derivatives, and so do the transport codes.

For how to see this in a comparison with finite differencing, you need to test out multiple finite differences, using the spline itself to get the values. The finite differences will converge to some value (and agree with the quantities calculated by ROSE for the same variable) but these are not the values output by EZSPLINE. If you guys still want I can try to set a test up inside the STELLOPT framework.

[lazersos]

@aaroncbader no need to go through the work. The problem is clear in the example I shared, I believe it exists, the question now is how to fix it.

[zhucaoxiang]

This is what I got from Sam's testing code. Just for a reference.

[lazersos]

This is a larger issue than just us. TRANSP uses the PSPLINE package as does the ASCOT code. ASCOT5 implemented the splines in C, and from what I understand did rigorous testing (or that's what was said). Most parts of our codes don't take advantage of the derivatives of the splines so the effect of these errors is minimized. However, the splines should be doing a better job than this.

[lazersos]

These results suggest to me that the issue is somehow related to what's assumed for the boundary condition calculation. I say this because the error grows toward the upward bound of the dimension. I should point out that these error plots are defined on grids with 16 times as many points in U and V as the splines knots themselves. The error in dF/Fu grows in V while the error in dF/dv grows in U, suggesting that this is somehow related to the derivative quantities.

[zhucaoxiang]

@lazersos This might be fixed if one spent some time on debugging. Apparently, it requires non-trivial effort.

@aaroncbader What does ROSE use to have the ability calculating gamma-C?

[lazersos]

@aaroncbader, @zhucaoxiang Looking back at this problem, the thing that bothers me is that the error in F grows at maximum x and y extents. If this was a simple algebraic issue I'd expect the errors to be symmetric. I think the issue may be with how the periodic boundary conditions are enforced.

[lazersos]

OK I looked into it a bit more. It seems that the issue is the underlying grid on which the initial function is sampled. The function I coded up goes like cos(mtheta+nphi) and we evaluate it from -pi to pi in the u and v dimensions. Let me summarize my findings

• 8x8 grid has max error of 0.6% in f(u,v), and 13.6% error in df/du, df/dv, and d^2f/dudv -12x12 grid has max error of 0.1% in f(u,v), and 6% error in df/du, df/dv, and d^2f/dudv -16x16 grid has max error of 0.08% in f(u,v), and 4.7% error in df/du, df/dv, and d^2f/dudv -32x32 grid has max error of 0.01% in f(u,v), and 1.5% error in df/du, df/dv, and d^2f/dudv

If I change the underlying function to be symmetric about (u,v)=(0,0) then the error becomes symmetric. I suspect what we're seeing is just the result of using these types of splines.

I think the takeaway message is that when going from a Fourier space representation to a spline grid we need to be careful to make sure the underlying grid possesses enough gridpoints to represent the Fourier space. Now when it comes to representing the quantity 8xM where m is the highest harmonic gives us <1% accuracy. However, for the derivative quantities we need much higher resolution on the order of 16xM to 32xM.

If we switch to-non-Hermite splines then our derivative terms go as

• 8x8 grid has max error of 0.1% in f(u,v), and 2.7% error in df/du, df/dv, and d^2f/dudv -12x12 grid has max error of 0.02% in f(u,v), and 0.4% error in df/du, df/dv, and d^2f/dudv -16x16 grid has max error of 0.005% in f(u,v), and 0.1% error in df/du, df/dv, and d^2f/dudv -32x32 grid has max error of 0.003% in f(u,v), and 0.01% error in df/du, df/dv, and d^2f/dudv

So my general comment is that we should probably switch to non-Hermite splines in parts of the code like stel_tools.f90 where the quantities we're splining over are defined everywhere and have no strong discontinuities.

[zhucaoxiang]

@lazersos If we use other spline packages to do the same test, will we get better convergence? Like the one from Scipy (https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.CubicHermiteSpline.html)

[lazersos]

@zhucaoxiang There are many ways to do splines, but I personally would like to keep the methods in-house instead of relying on 3rd party libraries. The reason for this is that the call to a spline is often in the RHS function of an ODE integrator. Eventually we'd like to run on GPU's and then I think the splines are the first candidate to 'live' on the GPU, and I'd prefer if we didn't have to rely on 3rd parties to do this.

[lazersos]

Also I'd like to be clear that EZspline can use piecewise linear, C2 cubic splines or C1 cubic splines. The above inaccuracy in the derivative terms comes from using C1 cubic splines. The fix is to use the C2 cubic spline methods (spl_obj%ishermite=0).

[lazersos]

OK I did some testing and there's a complication using C2 over C1. I tested both in the FIELDLINES code using the W7-X standard vacuum configuration as defined from the coils file. If I use C1 splines I clearly resolve the flux surfaces, boundary 5/5 island chain, and the flux surfaces beyond that. For the same grid resolution and construction methods I only get the core flux surfaces and no 5/5 island chain using C2. If I restrict the field generation so we don't evaluate fields larger than 5 T (mask away any spikes caused by the 1/x created by the coils being in the domain), then I get some parts of the 5/5 island chain to reappear.

So here's where we stand. Clearly the C2 splines provide more accurate derivatives for smooth functions. However, for field line tracing the presence of singularities or discontinuities in the domain causes ringing for C2. This results in poor approximation of the underlying function values. For codes like stel_tools.f90 and virtual_casing_mod.f90 we can use (and do) the C2 splines to provide the most accurate representation of the smooth underlying function. For the cylindrical grid quantities like those which are found in BEAMS3D and FIELDLINES we should use the C1 representation (and do).

It should be noted that the FIELDLINES code does not use derivatives of the splines in the ODE it solves but does use the derivative to provide a Jacobian for the ODE integrator. This appears to work incredibly well, better than the C2 splines.

In BEAMS3D the ODE includes both the magnetic field and derivatives of the magnetic field. Clearly the C1 derivatives are not as accurate as the C2 derivatives. However, the representation of the function values is substantially better than those of the C2 given the 'real' nature of the smooth magnetic field. So I believe we're justified in using the C1 representation.

We should entertain new 3D spline routines which can address this issue but for now I'm going to close this issue as I think it's come to a conclusion. That conclusion is that we need to research new spline methods which are both derivative preserving and can handle functions with discontinuities.

[dbindel]

Marking this here only in case someone comes back to this thread to pic up the last point on alternate spline methods. There are common spline families that take advantage of smoothness where present and also allow points of non-differentiability. One of the most common (used in CAD models and in isogeometric finite element analysis) is the non-uniform rational B-spline (NURBS).

[lazersos]

@dbindel I'd be happy to do this but the existing NURBS libraries don't seem particularly open source or easily interfaceable with Fortran. However, it may make a nice summer undergraduate project to write Fortran wrappers one of the NURBS C++ libraries.

[zhucaoxiang]

@dbindel @lazersos I will have a summer student working on splines for coils in FOCUS. After that, I will come back to this problem and see if I can make a contribution.