Local cubic spline interpolation for Vlasov-type equations on a multi-patch geometry

Technical Summary: Local Cubic Spline Interpolation for Vlasov-type Equations on a Multi-patch Geometry

Problem Statement
The paper addresses the numerical resolution of Vlasov-type equations, specifically guiding-center models, on complex multi-patch geometries typical of tokamak simulations. In these simulations, mesh lines are aligned with magnetic flux surfaces to minimize variations along field lines. This geometry creates specific challenges: the core features an O-point where mesh lines converge, and the edge features an X-point (separatrix) where lines diverge. While splitting the global domain into subdomains (patches) allows for local refinement and optimized uniform meshes within each patch, it introduces difficulties in maintaining solution regularity at patch interfaces. Previous studies using $C^0$ continuity (matching only function values) at interfaces were found to be potentially unstable for semi-Lagrangian advection schemes. The core problem is how to enforce higher-order continuity ($C^1$) and accurately compute interface derivatives for local cubic B-spline interpolations on arbitrary (uniform and non-uniform) multi-patch meshes without sacrificing the efficiency of local computations.

Methodology
The authors propose a backward semi-Lagrangian (BSL) method utilizing local cubic B-splines defined on each patch. To ensure stability and accuracy, the method enforces $C^1$ regularity at patch interfaces by matching both function values and their first derivatives.

1. Derivative Reconstruction: The central innovation is a generalized formula for computing interface derivatives. Building on a formula by Crouseilles et al. [11] designed for uniform meshes, the authors extend the approach to arbitrary non-uniform meshes using Hermite polynomials.

   • Exact Formulation: By imposing $C^2$ continuity between local Hermite polynomials defined on adjacent cells, the authors derive a recursive relation linking the derivative at an interface to the derivatives at neighboring interfaces and the function values. This results in a small global linear system (size $(N_p-1) \times (N_p-1)$ for $N_p$ patches) that must be solved to obtain derivatives equivalent to those of a global spline.
   • Approximate Formulation: The authors analyze the asymptotic behavior of the coefficients in the global system. They demonstrate that the influence of distant interfaces decays exponentially as the number of cells per patch increases. Consequently, for sufficiently large patches (e.g., $\ge 30$ cells), the dependency on distant derivatives becomes negligible (reaching machine precision). This allows for a truncated, local approximation where derivatives are computed using only a limited number of neighboring cells, avoiding the global solve while maintaining high precision.
   • Explicit Formulas: For the specific case of uniform patches, the authors provide explicit closed-form expressions for the coefficients, further optimizing computation.

2. 2D Extension: The method is extended to 2D multi-patch geometries, handling three cases:

   • Conforming Meshes: A direct tensor-product application of the 1D formulas to compute $r$-derivatives, $\theta$-derivatives, and cross-derivatives.
   • Non-Conforming Meshes: Where break points do not align across patches. The authors propose an algorithm that progressively computes derivatives on coarser interfaces and interpolates them onto refined patches to enforce $C^1$ continuity.
   • T-Joints: Connections involving three or more patches where an edge is only partially connected. The method adapts the order of derivative computation based on the geometry to ensure all necessary derivatives are available.

3. Implementation: The algorithm is implemented in the Gyselalib++ library. The Poisson equation is solved using a Finite Element Method (FEM) with cubic B-splines on an equivalent global spline, while the advection equation is solved using the proposed BSL method on local splines.

Key Contributions

• Generalization of Interface Derivative Formulas: The paper generalizes the derivative computation formula from uniform to arbitrary non-uniform meshes using Hermite polynomials, enabling the handling of non-conforming situations and local refinement.
• Stability Analysis: The authors demonstrate that $C^0$ regularity is insufficient for stability in BSL schemes on multi-patch domains and that $C^1$ regularity is required.
• Local vs. Global Trade-off: The work establishes that while an exact global system ensures $C^2$ continuity, a local approximation using a truncated stencil (e.g., 30 cells) achieves machine precision equivalence to the global solution, offering a computationally efficient alternative.
• Handling Complex Geometries: The method successfully handles T-joints and non-conforming meshes, which are necessary for modeling complex tokamak cross-sections involving both O-points and X-points.

Numerical Results
The method was validated using the Gyselalib++ library with the following tests:

1. Spline Interpolation: Tests on concentric patches (uniform and non-uniform) confirmed that the local spline representations converge to the equivalent global spline with order 4. The error in interface derivatives using the truncated approximation matched the exact formula once the number of cells exceeded $\approx 30$, reaching machine precision ($10^{-13}$).
2. Advection Tests:
   • Non-conforming patches: The local splines on a non-conforming mesh achieved accuracy comparable to a refined global mesh, despite using fewer total cells.
   • T-joints: In a 2D T-joint configuration, the local spline results agreed with the global spline results up to machine precision ($\approx 6 \times 10^{-13}$).
3. Guiding-Center Model (Diocotron Instability):
   • Simulations on conforming multi-patch grids reproduced the results of global spline simulations (growth rates and mass conservation) with errors consistent with the interpolation precision.
   • A simulation with a non-optimal mesh (depopulated central patch and shifted interfaces) showed that while the mesh choice affected the absolute error magnitude, the local spline method still tracked the global solution behavior correctly, validating the robustness of the approach.

Significance and Claims
The paper claims to provide a robust framework for solving Vlasov-type equations on multi-patch geometries, which are essential for realistic tokamak simulations. The primary significance lies in enabling the use of local cubic spline interpolation with Hermite boundary conditions to maintain stability and high-order accuracy across patch interfaces.

The authors emphasize that their method removes the constraint found in previous studies regarding the number of points per patch in conforming cases, replacing it with the requirement to solve a small global system (or a local approximation). They note that the method is currently validated for geometries with O-points. While the introduction of T-joints is a step toward connecting O-point and X-point geometries, the authors explicitly state that handling X-point singularities requires further study and is not yet fully resolved in this work. The method is presented as a significant step toward efficient, high-order semi-Lagrangian solvers for complex magnetic confinement geometries.