High-Order Regularity of the Potentials at the Axis Using High-Order B-Spline Finite Element Discretization

This work presents a rigorous theoretical and numerical framework for the smooth disretization of the potentials at the axis in toroidally confined plasmas. The finite element method (FEM) is widely used in various gyrokinetic codes such as GYGLES, ORB5, and EUTERPE. However, when using FEM to disretize the potentials numerical issues arise at the axis due to the tensor product B-spline basis functions not obeying the regularity condition at the axis. To circumvent this issue, an inner boundary is typically introduced to exclude modes close to the magnetic axis. This is not appropriate for nonlinear simulations involving zonal flows and nonlinear interactions – related to radial spreading of turbulence – when modes near the magnetic axis play a critical role.

While typically the so-called unicity condition is implemented to ensure C^0 regularity of the potentials at the axis, it does not guarantee C^0 regularity of the electromagnetic fields and thus introduces numerical errors.

A new set of basis functions, termed orthonormal smooth polar splines, is introduced to replace the standard tensor-product B-splines at the origin. These new basis functions are linear combinations of tensor-product B-splines and can be incorporated using a extraction method, requiring minimal modifications just limited to the solver part of existing codes.

The method is applicable to standard tensor-product B-splines, phase-factor modified B-splines, and hybrid Fourier–B-spline approaches. The proposed basis exhibits several advantageous properties:

1. It is orthonormal.
   
   It preserves charge conservation.
   
   It significantly reduces the condition number of the system matrix.
   
   It decreases the statistical variance at the origin in Particle-In-Cell (PIC) methods.
   
   It acts as a Fourier filter to suppress nonphysical noise at the origin.

Several numerical tests are presented, including convergence studies for the solution of the Helmholtz equation with both analytical and noisy (PIC-based) source terms. The method has been successfully implemented in the EUTERPE code, where it resolves numerical issues in the ITPA-TAE benchmark case. Additional results are in preparation for the conference.