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

Pauline Vidal, Emily Bourne, Virginie Grandgirard, Michel Mehrenberger,

Eric Sonnendrücker

► To cite this version:

Pauline Vidal, Emily Bourne, Virginie Grandgirard, Michel Mehrenberger, Eric Sonnendrücker. Local cubic spline interpolation for Vlasov-type equations on a multi-patch geometry. 2025. hal-05119228v1

HAL Id: hal-05119228

https://hal.science/hal-05119228v1

Preprint submitted on 18 Jun 2025 (v1), last revised 8 Jan 2026 (v2)

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Local cubic spline interpolation for Vlasov-type equations on a multi-patch geometry.

Pauline Vidal<fnref n="1" />

Emily Bourne<fnref n="2" />

Virginie Grandgirard<fnref n="3" />

Michel Mehrenberger<fnref n="4" />

Eric Sonnendrücker<fnref n="5" />

Last update: June, 18th, 2025

## Abstract

We present a semi-Lagrangian method for the numerical resolution of Vlasov-type equations on multi-patch meshes. Following [10], we employ a local cubic spline interpolation with Hermite boundary conditions between the patches. The derivative reconstruction is adapted to cope with non-uniform meshes as well as non-conforming situations. In the conforming case, the constraint of the number of points for each patch, found in [10], is removed; however, a small global system must now be solved. In that case, the local spline representations coincide with the corresponding global spline reconstruction. Alternatively, we can choose not to apply the global system and the derivatives can be approximated as in [10]. The influence of the most distant points diminishes as the number of points per patch increases. For uniform per patch configurations, a study of the explicit and asymptotic behavior of this influence has been led. The method is validated using a two-dimensional guiding-center model with an O-point [37]. All the numerical results are carried out in the Gyselalib++ library.

**Keywords** Multi-patch approach · Backward semi-Lagrangian method · Cubic B-splines · Hermite polynomials · T-joints.

**Mathematics Subject Classification (2020)** 65M25 · 65D07

**Acknowledgements** This work was carried out using the Gyselalib++ library. The project has received funding from the EoCoE-III project (Energy-oriented Centre of Excellence for Exascale HPC applications). This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Program (Grant Agreement No. 101052200 — EUROfusion). Views and opinions expressed are those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. Emily Bourne's salary was paid for by the EUROfusion Advanced Computing Hub (ACH).

## 1 Introduction

In numerical tokamak simulations with semi-Lagrangian schemes, we align the mesh lines on the flux surfaces of the magnetic field, since the variations of the different physical quantities along the magnetic

<footnote n="1">¹Corresponding author: Pauline Vidal</footnote>

<footnote>Max-Planck-Institut für Plasmaphysik, Garching, Germany Email: pauline.vidal@ipp.mpg.de ORCID: https://orcid.org/0009-0008-7233-6351</footnote>

<footnote n="2">²Emily Bourne</footnote>

<footnote>SCITAS, EPFL, CH-1015 Lausanne, Switzerland ORCID: https://orcid.org/0000-0002-3469-2338</footnote>

<footnote n="3">³Virginie Grandgirard</footnote>

<footnote>CEA, IRFM, 13108 Saint-Paul-lez-Durance Cedex, France ORCID: https://orcid.org/0000-0001-7821-9107</footnote>

<footnote n="4">⁴Michel Mehrenberger</footnote>

<footnote>Aix Marseille Univ, CNRS, I2M, Marseille, France ORCID: https://orcid.org/0000-0001-6376-409X</footnote>

<footnote n="5">⁵Eric Sonnendrücker</footnote>

<footnote>Max-Planck-Institut für Plasmaphysik, Garching, Germany, Technische Universität München, Zentrum Mathematik, Garching, Germany ORCID: https://orcid.org/0000-0002-8340-7230</footnote>
field lines are small [23]. These mesh lines separate the poloidal plane into different zones (see Fig. 1). In the core, the mesh lines converge to a central point called the O-point. On the edge, the separatrix (the boundary between the closed and open magnetic field lines) forms an X-point at its intersection. The mesh lines diverge from the X-point around it. These zones raise different numerical difficulties that are interesting to address using meshes adapted to these specificities. The multi-patch approach is an effective tool for this geometry type. In this approach, the multi-patch geometry corresponds to a global domain split into different subdomains called patches.

The multi-patch approach has interesting optimization properties. For instance, computations on uniform meshes can be more optimized than on non-uniform meshes (see [31] for an example of optimization on uniform meshes). Non-uniform meshes can reduce memory constraints, but they can also potentially increase the calculation requirements [4]. So, if we want to locally refine an area of a global domain, we can split it into multiple patches to be uniform piecewise and remain optimized. Moreover, with a 2D basis formed from the tensor product of 1D function bases, the multi-patch approach allows us to deal with non-conforming meshes (as is the case in Fig. 1). Studies have already been conducted on non-conforming refinement on a multi-patch geometry for finite element methods (FEM) with B-splines [5, 19, 20]. Splitting a domain into independent subdomains is also interesting for parallelization. We can cite the work in [10, 22] on MPI parallelization for Vlasov-Poisson equations with a domain decomposition and B-splines. However, parallelizing operators on a poloidal plane with both an O-point and an X-point may be difficult here to achieve a good load balance. This will therefore not be studied in this article.

<figure><img src="image_2.png" /><figcaption>Figure 1: On the left, an example of global domain split into 9 patches. Mesh for the ITER tokamak based on data from SOLEDGE3X and generated by K. Obrejan (https://soledge3x.onrender.com, article [32]). At the center, a zoom on the core with a reference mapping with O-point. On the right, a zoom on the separatrix (---) intersection forming the X-point with a reference mapping.</figcaption></figure>

Once the domain is split into multiple patches, we are interested in simulating particle motion. In a tokamak, the time evolution of the particle distribution functions is governed by the Vlasov equation, a 6D kinetic equation (3D in space and 3D in velocity). It can be averaged over the cyclotron motion. This gives the gyrokinetic equation which describes a charge distribution function in a 5D phase space with 3 dimensions in space and 2 dimensions in velocity [15, 16, 24],

$$ \partial_t f(t, \mathbf{x}, v_{\parallel}, \mu) + \mathbf{v}_{GC}(t, \mathbf{x}, v_{\parallel}, \mu) \cdot \nabla f(t, \mathbf{x}, v_{\parallel}, \mu) + a_{GC}(t, \mathbf{x}, v_{\parallel}, \mu) \partial_{v_{\parallel}} f(t, \mathbf{x}, v_{\parallel}, \mu) = 0, \quad (1) $$

with f the gyro-center distribution function of plasma particles, $\mathbf{v}_{GC}$ the gyrokinetic velocity, $a_{GC}$ the gyrokinetic acceleration defined on position $\mathbf{x} \in \mathbb{R}^3$, $v_{\parallel} \in \mathbb{R}$ the parallel velocity, and $\mu \in \mathbb{R}^+$ the modified magnetic moment. To solve the Vlasov equation with the gyrokinetic approximation, we use an operator splitting method to combine different advection equations along each dimension. The main difficulty arises at special points (the O-point and the X-point) which are located in the poloidal cross-section.
So, to study these poloidal cross-section singularities, we focus in this paper on the following 2D model with guiding-center equations where the advection is in the poloidal cross-section,

$$ \left\{ \begin{array}{l} \partial_t \rho + \mathbf{v}_{GC} \cdot \nabla \rho = 0, \\ \Delta \phi = -\rho, \end{array} \right. \qquad (2) $$

with $\rho$ the density of the plasma particles. In the Poisson equation in (2), $\phi$ is the electric potential of the electric field $\mathbf{E} = -\nabla\phi$. Additionally, the magnetic field is supposed to be a unit vector, $\mathbf{B} = \mathbf{e}_z$, so the drift velocity $\mathbf{v}_{GC} = \frac{\mathbf{E} \times \mathbf{B}}{|\mathbf{B}^2|} = \mathbf{E} \times \mathbf{e}_z$.

The numerical method and the results that will be presented are implemented in the Gyselalib++ library [2] which is a C++ rewrite of the GYSELA code [15].

In the library, the Poisson equation on the poloidal cross-section in (2) is solved with a FEM [37]. The advection equation is solved with a backward semi-Lagrangian method (BSL) [34]. The basis functions used in the FEM and the interpolation are cubic B-splines. They have good performance and interpolation qualities [35]. So, we will study the advection on multiple patches solved with a BSL on cubic B-splines to implement the method in Gyselalib++.

In 2019, E. Zoni and Y. Güçlü, introduced methods to solve guiding-center equations on O-point geometries [37]. They present methods to solve Poisson equations with an O-point by modifying the basis functions around the O-point. The advection equation is solved in the physical domain by linearizing the advection field around the O-point. In 2023, E. Bourne et al. studied different Poisson solvers [3] and implemented the Spline FEM solver in Gyselalib++. We use the work of these two articles as a starting point and adapt it to multi-patch geometries with O-point.

The main difficulty in the BSL method on a multi-patch domain lies in the interpolation and dealing with the patch connections. We manage this issue here by using Hermite boundary conditions on the local splines defined on each patch. These conditions confer a $C^1$ regularity to the global piecewise spline representation. However, to impose the interface derivatives, we need to compute them. There are different methods to compute derivatives. For instance, we can apply an epsilon method [33]. It consists of adding two interpolation points at $\pm\epsilon$ of the interfaces and computing the derivatives using the function values at these points. To use this method, we then need to add two interpolation points for each derivative that we want to compute. Another idea specific to the backward semi-Lagrangian method presented by E. Bourne in her PhD thesis [1] is to advect the derivatives simultaneously with the function. Starting from the relation in the backward semi-Lagrangian method, $f(t, \mathbf{X}(t; s, \mathbf{x})) = f(s, \mathbf{x}), \forall s \in \mathbb{R}, \mathbf{x} \in \Omega$, a relation for the derivatives can be directly deduced: $\partial_{\mathbf{x}}(f(t, \mathbf{X}(t; s, \mathbf{x}))) = \partial_{\mathbf{x}}\mathbf{X}(t; s, \mathbf{x})\partial_{\mathbf{x}}f(s, \mathbf{x}), \forall s \in \mathbb{R}, \mathbf{x} \in \Omega$. However, this method additionally requires the derivatives of the characteristics to be determined $\partial_{\mathbf{x}}\mathbf{X}(t; s, \mathbf{x})$.

This article is devoted to the generalization of a formula developed by N. Crouseilles, G. Latu and E. Sonnendrücker, in [10]. It consists of calculating the derivatives from the function values for parallel computation purposes (see also [26, 35]). It is also interesting the multi-patch approach (even without parallelization considerations). This formula was developed for uniform meshes. We use Hermite polynomials to handle non-uniform meshes and non-conforming situations. The generalization introduced in this paper allows us to calculate the derivatives of an equivalent global spline and thus reach $C^2$ continuity in the conforming case. The use of derivatives equal to those of an equivalent global spline comes at the cost of the resolution of a small global system, but the derivatives can also be approximated. The weight of the values at the most distant points in the calculation of the derivative decreases as the number of points per patch increases.

We can draw a parallel between this interface derivatives approach and the hybridizable discontinuous Galerkin (HDG) approach [6–8]. The two approaches add degrees of freedom and split the domain to locally solve a problem. The added interface derivatives in the multi-patch approach can be seen as elements of the numerical traces in the static condensation method.

This paper is organized as follows. Some reminders of the backward semi-Lagrangian method and the B-splines can be found in Section 2. Section 3 presents the generalized method that we suggest. To advect with a backward semi-Lagrangian method on 1D multi-patch geometries, we develop a formulation to calculate the derivatives of an equivalent global spline. We also study an explicit formulation and its
asymptotic behavior in the uniform per patch case. Section 4 describes the application of the method on 2D multi-patch geometries. Section 5 discusses some numerical results on spline interpolations, advections, and the guiding-center equations run with the Gyselalib++ code on multi-patch geometries with O-point. We apply the method only to geometries with an O-point using the work in [1,37], but not to geometries with X-points as further studies are needed.

# 2 Semi-Lagrangian scheme on a global domain

We briefly recall the backward semi-Lagrangian method on a global domain without multi-patch in the first subsection 2.1. For the interpolation step, we use B-splines which are described in the second subsection 2.2. These B-splines are defined on a logical domain and then mapped to a physical domain. These terms are defined in subsection 2.3.

## 2.1 Backward semi-Lagrangian method

The hyperbolic advection equation in (2) is solved with a backward semi-Lagrangian method. One advantage of semi-Lagrangian methods is that they are less restrictive on the time step, as there is no Courant-Friedrichs-Lewy (CFL) condition [9]. The method is detailed in [34]. For a typical advection equation as follows,

$$ \partial_t f(t, \mathbf{x}) + \mathbf{A}(t, \mathbf{x}) \cdot \nabla f(t, \mathbf{x}) = 0, \quad (3) $$

with a given advection field $\mathbf{A}$, the method consists of adopting the Lagrangian point of view. We define the characteristics as the solution $\mathbf{X}$ of the following ordinary differential equation (ODE),

$$ \left\{ \begin{array}{ll} \partial_s \mathbf{X}(s; t, \mathbf{x}) = \mathbf{A}(s, \mathbf{X}(s; t, \mathbf{x})), & s \ge 0, \mathbf{x} \in \Omega, \\ \mathbf{X}(t; t, \mathbf{x}) = \mathbf{x}, & \mathbf{x} \in \Omega. \end{array} \right. $$

Thanks to the property of conservation along the characteristics, we have the following relation for the solution of (3), with fixed $t \in \mathbb{R}$, $\mathbf{x} \in \Omega$,

$$ f(t, \mathbf{X}(s; t, \mathbf{x})) = f(t, \mathbf{x}), \quad \forall t \in \mathbb{R}. $$

Knowing the function at $s = t - \Delta t$, we can deduce it at $t$. We need to solve backward the equations of the characteristics for every interpolation point. The solutions of the ODE are called the feet of the characteristics. We interpolate the function known at $t - \Delta t$ at the feet of the characteristics. Each function value corresponds to the function value at $t$ at the interpolation point $\mathbf{x}$. To adapt the method to a multi-patch geometry, the main difficulty is in the interpolation.

## 2.2 B-splines interpolation

The chosen interpolation method for the function at $t$ is an interpolation on cubic B-splines. These basis functions are well described in [14, 25, 28, 35]. In 1D space, the spline $s$ is given by

$$ s(x) = \sum_{k=0}^{N_c+2} c_k b_k(x), \quad x \in \Omega $$

with $N_c$ the number of cells, $\{b_k\}_{k=0,...,N_c+2}$ the basis functions, called the B-splines, and $\{c_k\}_{k=0,...,N_c+2}$ coefficients to determine (also called control points).

We call break points, $\{x_{b,i}\}_{i=0}^{N_c}$ the boundaries of the cells. There are $N_c + 1$ break points. We add $N_c + 6$ points, $\{k_i\}_{i=0}^{N_c+5}$ called knots to build the B-splines. The central knot sequence is defined on the break points; the additional points are placed at regular intervals outside the domain (this is a common choice in the uniform case) or are defined with repetitions of the first and the last break points (see Appendix A. for more details). The latter case is referred to as an open knots sequence.

Whatever knot sequence we choose, this system has $N_c+3$ degrees of freedom to fix. We then introduce interpolation points $\{x_i\}_{i \in N}$. In the Gyselalib++ code, three interpolation types are implemented.
**Interpolation with Hermite boundary conditions.** If the function values and its derivatives are known on the boundaries, a way to fix the $N_c + 3$ degrees of freedom is to define $N_c + 1$ interpolation points $\{x_i\}_{i=0}^{N_c}$ and fix the two last degrees of freedom with the derivatives at the boundaries as follows.

$$ \left\{ \begin{array}{ll} f_i = \sum_{k=0}^{N_c+2} c_k b_k(x_i), & i = 0, \dots, N_c, \\ (\partial_x f)_i = \sum_{k=0}^{N_c+2} c_k b'_k(x_i), & i = 0, N_c. \end{array} \right. $$

**Interpolation using periodic boundary conditions.** If the derivatives are unknown on the boundaries, but the domain is periodic, then we can define $N_c + 1$ interpolation points and use the periodicity property to fix all the degrees of freedom.

$$ f_i = \sum_{k=0}^{N_c+2} c_k b_k(x_i), \quad i = 0, \dots, N_c. $$

**Interpolation points as closure conditions.** If we do not know the derivatives at the boundaries and the domain is not periodic, we add more interpolation points to fix all the degrees of freedom. We generally use Greville points which are defined from the knot sequence as follows, (here $d = 3$). In this paper, we use open knot sequences.

$$ x_i = \frac{1}{d} \sum_{j=i+1}^{i+d} k_j. \qquad (4) $$

Using the Greville points as interpolation points is usually considered as an optimal solution for accuracy and efficiency, [14,17,18,21]. In the uniform case, the Greville points (except at the boundaries) are the break points for cubic spline. It would have been cell midpoints for even degree splines. With $N_c + 6$ knots, it gives $N_c + 3$ interpolation points $\{x_i\}_{i=0}^{N_c+2}$. This fixes all the degrees of freedom.

$$ f_i = \sum_{k=0}^{N_c+2} c_k b_k(x_i), \quad i = 0, \dots, N_c + 2 $$

Once the degrees of freedom are fixed, we can determine the B-spline coefficients $\{c_k\}_{k=0}^{N_c+2}$. They define a spline representation of the function at $t$ that we can evaluate to determine the function at $t$ in the backward semi-Lagrangian method, 2.1. This spline representation does not depend on the knot sequence used to define the B-splines.

**Remark.** Using a uniform knot sequence provides identical B-splines except for a support shift. It is then easier to optimize the code. So, it is more interesting to split a global non-uniform domain into several uniform patches to optimize it.

## 2.3 Logical and physical domains

Let's introduce some notations. The B-splines are defined on a rectangle-shape domain called the *logical domain* and denoted $\hat{\Omega}$. The 2D B-spline basis is a tensor product of two 1D B-spline bases. This logical domain is mapped to the *physical domain* denoted $\Omega$. An example of a physical domain and its associated logical domain is shown in Fig. 2. The mapping from the logical domain to the physical domain and converting curvilinear coordinates into Cartesian coordinates is noted $\mathcal{F}: \hat{\Omega} \to \Omega$.

As the B-splines are defined on the logical domain, all the interpolations are performed in the logical domain.
<figure><img src="image_3.png" /><figcaption>Figure 2: Logical and physical domains. We note $(r, \theta)$ the curvilinear coordinates on the logical domain and $(x, y)$ the Cartesian coordinates on the physical domain.</figcaption></figure>

# 3 Semi-Lagrangian scheme on a 1D multi-patch domain

When we split a global domain into multiple patches, we need to be aware of the regularity at the patch interfaces. We start by defining the multi-patch domain in 3.1. We explain our choice of a $C^1$ regularity in 3.2. The $C^1$ regularity implies fixing the interface derivatives. The authors in [10] already gave us a formula to compute the derivatives with a precision of $10^{-6}$ for a uniform global mesh. Subsection 3.3 presents a new generalization of this formula to arbitrary meshes and any precision. We end this section with the algorithm for a backward semi-Lagrangian method adapted to a multi-patch mesh, 3.4.

## 3.1 Definition of multi-patch domain

We define a multi-patch domain as a physical domain split into different subdomains $\Omega = \cup_k \Omega^k$ (as the colors show in Fig. 1). The interiors of these subdomains are disjoint. The domain recovery is only on the edges. We build equivalent logical domains $\tilde{\Omega}^k$ mapping to each physical subdomain. These logical domains are called *patches*. The coordinates in the different patches are not necessarily equivalent. The logical domains are not necessarily defined in the same dimensions. On each patch, we define independent spline representations of the function. We called these splines *local splines*. For simplicity, all the local splines will be cubic piecewise polynomials.

## 3.2 Local spline regularity at the interfaces

The local splines are defined independently of each other. However, we want the solution of the guiding-center model, (2) to be continuous in the global domain. We then need to impose a regularity on the interfaces between two patches to connect the local splines one another. The first idea is to impose a $C^0$ regularity.

### 3.2.1 $C^0$ regularity

If we impose a $C^0$ regularity, then we do not need to define derivatives on the boundaries of each local patch. We can define local splines with interpolation points as closure conditions. As mentioned previously in 2.2, optimal interpolation points are the Greville points. However, despite Greville points, we have noticed that using only a $C^0$ regularity can be unstable.

We start from the results in the article [27]. In this article, they study a 1D advection on a periodic domain with a constant advection field. The domain $\tilde{\Omega}$ is split into $N_p$ patches. Unlike this article, we solve the advection with a backward semi-Lagrangian method [34]. We choose to work with a uniform mesh for simplicity. The interpolation is realized with piecewise polynomial basis functions of degree $d = 3$ over $N_c$ cells for each patch,

$$ \varphi_j(x) = \sum_{k=0}^{d} \varphi_{j,k} x^k, \quad j = 0, \dots, N_c - 1. \qquad (5) $$

The $C^{d-1}$ continuity condition of the basis functions at the cell interfaces gives $d(N_c - 1)$ constraints. The coefficients $\{\varphi_{j,k}\}_{j=0, \dots, N_c-1, k=0, \dots, d}$ are fixed with these constraints and the interpolation at the
$N_c + d$ points of the function $f$ we want to advect in (3). To avoid duplication, each patch stores the $N_c+2$ first function values and the last one is taken from the next patch. An illustration of the advection step on these patches is given in Fig. 3. The function values $\mathbf{f}^{n+1}$ at $t^{n+1} = t^n + \Delta t^n$ are deduced from the function values $\mathbf{f}^n$ at $t^n$ for a given time discretization $\{t^n\}_{n \ge 0}$. We study the case of the constant velocity of the advection field $A = v$ such that $v < 0$ and the displacement $|v\Delta t|$ is smaller than a patch length. We note $\mathbf{f}_p^n$ the set of function values on the patch $p$ of the function at $t^n$.

<figure><img src="image_4.png" /><figcaption>Figure 3: Constant advection on a periodic domain split into $N_p$ patches. The function on patch $p$ is advected with a constant velocity $v < 0$. The feet of the characteristics land in the patch $p$ and $p+1$. To determine the function at $t^{n+1}$, we interpolate the function representation at $t^n$ using the information from patch $p, p+1$ and $p+2$.</figcaption></figure>

The function values on patch $p$ at $t^{n+1}$ depend on the values of the function at $t^n$ on patches $p, p+1$ and $p+2$. As the advection is constant and the domain periodic, we can write $\mathbf{f}^{n+1} = \mathbf{A}\mathbf{f}^n$, with $\mathbf{A}$ a block matrix,

$$ \mathbf{A} = \begin{bmatrix} \mathbf{A}_0 & \mathbf{A}_1 & \mathbf{A}_2 & \dots & \mathbf{A}_{N_p-1} \\ 0 & \mathbf{A}_0 & \mathbf{A}_1 & \mathbf{A}_2 & \dots \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & 0 & \mathbf{A}_0 & \mathbf{A}_1 \\ \mathbf{A}_2 & 0 & \dots & 0 & \mathbf{A}_0 \\ \mathbf{A}_1 & \mathbf{A}_2 & 0 & \dots & \mathbf{A}_{N_p-1} \end{bmatrix}. \quad (6) $$

It is a circulant matrix, we can therefore apply the Fourier transform to diagonalize it,

$$ (\mathbb{R}^{N_c+2})^{N_p} \rightarrow (\mathbb{C}^{N_c+2})^{N_p}, \quad (\mathbf{f}_0, \mathbf{f}_1, \dots, \mathbf{f}_{N_p-1}) \mapsto (\hat{\mathbf{f}}_0, \hat{\mathbf{f}}_1, \dots, \hat{\mathbf{f}}_{N_p-1}) \quad (7) $$

with the relations

$$ \hat{\mathbf{f}}_p = \sum_{j=0}^{N_p-1} e^{-2i\pi p j/N_p} \mathbf{f}_j, \quad \hat{\mathbf{A}}_p = \sum_{j=0}^{N_p-1} e^{-2i\pi p j/N_p} \mathbf{A}_j. \quad (8) $$

However, it can be shown that for $k = 1, N_p = 5$ and $N_c = 3$, there exists a displacement $|v\Delta t| = 405/100000$, such that the maximum of the absolute value of the eigenvalues of $\hat{\mathbf{A}}_k$ is strictly superior to 1 with an excess of $10^{-7}$. This counterexample shows that we do not have the eigenvalues strictly inferior to 1. We may encounter other examples with a higher excess, so a bigger instability. So, this case shows that applying a $\mathbb{C}^0$ regularity can be unstable.

Thus, we need a stronger regularity between the patches.

### 3.2.2 $\mathbb{C}^1$ regularity

As a $\mathbb{C}^0$ regularity may sometimes lead to instabilities, we enforce a $\mathbb{C}^1$ regularity at the interfaces. We impose the equality of the derivatives and the spline representations of the function at the interfaces. We then apply an interpolation with Hermite boundary conditions on the local splines to ensure this. However, we do not know the derivatives at the interfaces. To build the spline representations, we need to determine them.
We suggest in this article a new generalization of a formula developed in [10] which consists of calculating the derivatives from twenty function values. This formula is given for uniform meshes. We use Hermite polynomials to extend it to arbitrary meshes. The new generalization introduced in this paper allows us to calculate derivatives equal to those of an equivalent global spline and even reach $C^2$ continuity in the conforming case.

## 3.3 Computation of interface derivatives

### 3.3.1 Assumptions

Before starting the study, we describe the assumptions that are necessary for this work. The study conducted here supposes that the spline representations use cubic B-splines.

Another remark concerns the interpolation points. As we will see, to define Hermite polynomials between two interpolation points, the spline representation between the points needs to be polynomial and not piecewise polynomial. This means that each break point also needs to be an interpolation point. While Greville points are found at break points when uniform cells are used this is not the case generally. In order to fulfil this condition, we cannot use Greville points as interpolation points in the non-uniform case.

### 3.3.2 Application of Hermite polynomials

The authors in [10] wanted to parallelize the computations in the semi-Lagrangian method using a domain decomposition. They split a global spline into local splines and used (9) to be as close as possible to the initial global spline. They developed a formula to compute the derivative $s'$ of the spline representation (34) from the function values {$f_{i+k}$}$_k$ at an interpolation point $x_i$:

$$s'(x_i) = \sum_{j=-10}^{10} \frac{\tilde{\omega}_j}{\Delta x} f_{i+j}, \qquad (9)$$

with given weights {$\tilde{\omega}_j$}$_j$ ($\tilde{\omega}_0 = 0$) (see Table 1 in [10]), and $\Delta x$ the uniform length between the interpolation points.

This formula is interesting for our multi-patch geometry. It is an approximation of the spline derivative with a precision of $10^{-6}$. The coefficients {$\tilde{\omega}_k$}$_k$ are given for a uniform mesh. However, it does not allow us to obtain better precision or deal with non-uniform cases. If we use this formula to compute the derivatives at the interface of two patches, both have to have the same uniformity. In practice, we would like to locally refine some patches. So, some interfaces will connect zones with two different cell lengths. In the following sections, we will generalize this formula to non-uniform meshes and extend it to be able to use it with any precision.

### 3.3.3 Relation between derivatives at three consecutive points

To generalize the formula (9), we define Hermite polynomials. The cubic Hermite polynomials interpolate the function values and the function derivatives at the boundaries. They are defined as follows for $x \in [a, b]$, $a, b \in \mathbb{R}$,

$$P(a + x(b - a)) = f(a)H_0(x) + f(b)H_1(x) + f'(a)(b - a)K_0(x) + f'(b)(b - a)K_1(x) \qquad (10)$$

with the basis functions

$$H_0(x) = (1-x)^2(1+2x), \quad H_1(x) = x^2(3-2x), \quad K_0(x) = (1-x)^2x, \quad K_1(x) = x^2(x-1). \quad (11)$$

We start by defining two Hermite polynomials $P_-$ and $P_+$ on two domains defined using three consecutive interpolation points (see Fig. 4). The polynomials are fixed by the boundary conditions ($s(x_{i-1}^{-})$, $s'(x_{i-1}^{-})$, $s(x_i^{-})$, $s'(x_i^{-})$ and $s(x_{i+1}^{+})$, $s'(x_{i+1}^{+})$, $s(x_{i+1}^{+})$, $s'(x_{i+1}^{+})$). Here, $s$ can represent the same spline on the two cells, or two different local splines. To stay as general as possible, we do not impose that the
interpolation points $x_i^-$ and $x_i^+$ have the same value. This is particularly the case when we define functions on different logical domains which connect via the physical domain, but do not share coordinate systems. The indices of $x_i^-$ and $x_i^+$ could also be different, e.g. $x_i^-$ and $x_{i+}^+$. But to lighten the notation, we will only write $i$.

We also note in this section the $x$-dimension the dimension which we split into multiple patches. This reasoning can be applied equally well to both the $r$-dimension and the $\theta$-dimension of the patches.

<figure><img src="image_5.png" /><figcaption>Figure 4: Two Hermite polynomials: $P_-$ defined on $[x_{i-1}^-, x_i^-]$ and $P_+$ defined on $[x_i^+, x_{i+1}^+]$ with functions connected across patches between $x_i^-$ and $x_i^+$.</figcaption></figure>

We impose a $\mathcal{C}^2$ continuity between the two polynomials. Imposing it enforces the Hermite polynomials to be equal to an equivalent global cubic spline defined on the whole domain $[x_{i-1}^-, x_i^-] \cup [x_i^+, x_{i+1}^+]$. Moreover, it gives us a relation between the derivatives. Indeed, if we impose

$$P''_-(x_i^-) = P''_-(x_i^-),$$

with

$$ \left\{ \begin{aligned} P''_-(x_i^-) &= \frac{1}{(x_i^- - x_{i-1}^-)^2} [s(x_{i-1}^-)H_0''(1) + s(x_i^-)H_1''(1) + s'(x_{i-1}^-)(x_i^- - x_{i-1}^-)K_0''(1) + s'(x_i^-)(x_i^- - x_{i-1}^-)K_1''(1)] \\ P''_+(x_i^+) &= \frac{1}{(x_i^- - x_{i-1}^-)^2} [s(x_i^+)H_0''(0) + s(x_{i+1}^+)H_1''(0) + s'(x_i^+)(x_{i+1}^+ - x_i^+)K_0''(0) + s'(x_{i+1}^+)(x_{i+1}^+ - x_i^+)K_1''(0)] \end{aligned} \right. $$

and write $\Delta x_i^- = x_i^- - x_{i-1}^-$ and $\Delta x_i^+ = x_{i+1}^+ - x_i^+$, the relation $P''_-(x_i^-) = P''_-(x_i^-)$ becomes

$$ \frac{1}{(\Delta x_i^-)^2} [6s(x_{i-1}^-) - 6s(x_i^-) + 2\Delta x_i^- s'(x_{i-1}^-) + 4\Delta x_i^- s'(x_i^-)] \\ = \frac{1}{(\Delta x_i^+)^2} [-6s(x_i^+) + 6s(x_{i+1}^+) - 4\Delta x_i^+ s'(x_i^+) - 2\Delta x_i^+ s'(x_{i+1}^+)] $$

Reordered knowing that $s'(x_i^-) = s'(x_i^+)$ and $s(x_i^-) = s(x_i^+)$ and that the values of the spline at the interpolation points are the function values, the relation is

$$ s'(x_i^-) = \frac{3}{2} \frac{1}{\Delta x_i^+ + \Delta x_i^-} \left[ \frac{\Delta x_i^-}{\Delta x_i^+} f_{i+1}^+ + \left( \frac{\Delta x_i^+}{\Delta x_i^-} - \frac{\Delta x_i^-}{\Delta x_i^+} \right) f_i^+ - \frac{\Delta x_i^+}{\Delta x_i^-} f_{i-1}^- \right] \\ - \frac{1}{2} \frac{1}{\Delta x_i^+ + \Delta x_i^-} \left[ \Delta x_i^- s'(x_{i+1}^+) + \Delta x_i^+ s'(x_{i-1}^-) \right]. $$

We then obtain a relation between the derivatives at three consecutive interpolation points. This relation gives the exact derivative of the equivalent global spline defined on the whole domain $[x_{i-1}^-, x_i^-] \cup [x_i^+, x_{i+1}^+]$. For clarity, we introduce the coefficients $\gamma_i$, $\alpha_i$ and $\beta_i$ such that,

$$ s'(x_i^-) = \gamma_i + \alpha_i s'(x_{i+1}^+) + \beta_i s'(x_{i-1}^-) \quad (12) $$

with

$$ \left\{ \begin{aligned} \gamma_i &= \frac{3}{2} \frac{1}{\Delta x_i^+ + \Delta x_i^-} \left[ \frac{\Delta x_i^-}{\Delta x_i^+} f_{i+1}^+ + \left( \frac{\Delta x_i^+}{\Delta x_i^-} - \frac{\Delta x_i^-}{\Delta x_i^+} \right) f_i^+ - \frac{\Delta x_i^+}{\Delta x_i^-} f_{i-1}^- \right], \\ \alpha_i &= -\frac{1}{2} \frac{\Delta x_i^-}{\Delta x_i^+ + \Delta x_i^-}, \\ \beta_i &= -\frac{1}{2} \frac{\Delta x_i^+}{\Delta x_i^+ + \Delta x_i^-}. \end{aligned} \right. $$
**Remark 1.** It is interesting to remark that the coefficients $\alpha_i$ and $\beta_i$ depend only on the proportionality of the distances between the interpolation points. The coefficient $\gamma_i$ is a linear combination of the function values and is inversely proportional to the distance between the interpolation points.

**Remark 2.** Note also that we wanted $C^1$ regularity, and this formula gives us an improved $C^2$ regularity.

### 3.3.4 Extension of the relation

Let's take the calculations a step further to extend this formula to a relation between derivatives further from each other. We first go forward in the right direction and establish a relation between the derivatives at $x_{i-1}$, $x_i$ and $x_{i+N}$ with $N^R \in \mathbb{N}^*$. To simplify the notation, we remove the superscripts $\cdot$ and $\cdot^+$, knowing that the value of the points can still be different on the left and on the right but the values of the function and its derivatives are the same. We can apply a recursion using the relation (12) at different interpolation points.

The notation of type $c_{n^l, n^r}^i$ corresponds to a relation between the derivative at $s'(x_i)$, the derivative at $n^l$ points to the left $s'(x_{i-n^l})$ (i.e. spaced by $n^l$ cells), and the derivative at $n^r$ points to the right $s'(x_{i+n^r})$ (i.e. spaced by $n^r$ cells).

We can prove the following relation for the derivatives at $x_{i-1}$, $x_i$ and $x_{i+n}$ with $n \in \mathbb{N}^*$,

$$s'(x_i) = c_{1,n}^i + a_{1,n}^i s'(x_{i+n}) + b_{1,n}^i s'(x_{i-1}).$$

with $c_{1,n}^i$, $a_{1,n}^i$ and $b_{1,n}^i$, coefficients depending on the chosen interpolation point $x_i$ and the number of cells $n$.

Indeed, for $n=1$, it corresponds to the relation (12), and the coefficients are

$$c_{1,1}^i = \gamma_i, \quad a_{1,1}^i = \alpha_i, \quad b_{1,1}^i = \beta_i. \qquad (13)$$

For $n=2$, we use the relation (12) at $x_i$ and $x_{i+1}$, and inject the equation (12) at $x_{i+1}$ into the equation (12) at $x_i$ to get

$$s'(x_i) = \frac{1}{1 - \alpha_i \beta_{i+1}} [\gamma_i + \alpha_i \gamma_{i+1}] + \frac{\alpha_i \alpha_{i+1}}{1 - \alpha_i \beta_{i+1}} s'(x_{i+2}) + \frac{\beta_i}{1 - \alpha_i \beta_{i+1}} s'(x_{i-1}),$$

so the coefficients are

$$c_{1,2}^i = \frac{1}{1 - \alpha_i \beta_{i+1}} [\gamma_i + \alpha_i \gamma_{i+1}], \quad a_{1,2}^i = \frac{\alpha_i \alpha_{i+1}}{1 - \alpha_i \beta_{i+1}}, \quad b_{1,2}^i = \frac{\beta_i}{1 - \alpha_i \beta_{i+1}}. \qquad (14)$$

For $n > 2$, we apply the same method with the three following relations

$$\left\{ \begin{array}{l} s'(x_i) = c_{1,n-1}^i + a_{1,n-1}^i s'(x_{i+n-1}) + b_{1,n-1}^i s'(x_{i-1}), \\ s'(x_i) = c_{1,n,1}^i + a_{1,n,1}^i s'(x_{i+n}) + b_{1,n,1}^i s'(x_{i-1}), \\ s'(x_{i+n}) = \gamma_{i+n} + \alpha_{i+n} s'(x_{i+n+1}) + \beta_{i+n} s'(x_{i+n-1}). \end{array} \right.$$

Injecting the third into the second one gives us,

$$s'(x_i) = c_{1,n}^i + a_{1,n}^i \gamma_{i+n} + a_{1,n}^i \alpha_{i+n} s'(x_{i+n+1}) + a_{1,n}^i \beta_{i+n} s'(x_{i+n-1}) + b_{1,n}^i s'(x_{i-1}),$$

and using the first relation to replace $s'(x_{i+n-1})$,

$$s'(x_i) = \frac{1}{1 - \beta_{i+n} \frac{a_{1,n}^i}{a_{1,n-1}^i}} \left[ c_{1,n}^i + a_{1,n}^i \gamma_{i+n} - a_{1,n}^i \beta_{i+n} \frac{1}{a_{1,n-1}^i} c_{1,n-1}^i \right] \\ + \frac{a_{1,n}^i \alpha_{i+n}}{1 - \beta_{i+n} \frac{a_{1,n}^i}{a_{1,n-1}^i}} s'(x_{i+n+1}) + \frac{1}{1 - \beta_{i+n} \frac{a_{1,n}^i}{a_{1,n-1}^i}} \left[ b_{1,n}^i - \beta_{i+n} \frac{a_{1,n}^i}{a_{1,n-1}^i} b_{1,n-1}^i \right] s'(x_{i-1}).$$
The coefficients at $n > 2$ are then given by

$$ egin{cases} c_{1,n+1}^i = \frac{1}{1 - \beta_{i+n} \frac{a_{1,n}^i}{a_{1,n-1}^i}} \left[ c_{1,n}^i + a_{1,n}^i \gamma_{i+n} - \beta_{i+n} \frac{a_{1,n}^i}{a_{1,n-1}^i} c_{1,n-1}^i \right] \\ a_{1,n+1}^i = \frac{a_{1,n}^i \alpha_{i+n}}{1 - \beta_{i+n} \frac{a_{1,n}^i}{a_{1,n-1}^i}}, \\ b_{1,n+1}^i = \frac{1}{1 - \beta_{i+n} \frac{a_{1,n}^i}{a_{1,n-1}^i}} \left[ b_{1,n}^i - \beta_{i+n} \frac{a_{1,n}^i}{a_{1,n-1}^i} b_{1,n-1}^i \right]. \end{cases} \quad (15) $$

Applying this recursion, we find the desired relation between the derivatives at $x_{i-1}$, $x_i$ and $x_{i+N^R}$,

$$ s'(x_i) = c_{1,N^R}^i + a_{1,N^R}^i s'(x_{i+N^R}) + b_{1,N^R}^i s'(x_{i-1}). $$

As a second step, we handle the points to the left and establish a relation between the derivatives at $x_{i-N^L}$, $x_i$ and $x_{i+N^R}$ with $N^L \in \mathbb{N}^*$. We can apply the same method going backward with a new initial condition. We write the relation linking the derivatives at $x_{i-m}$, $x_i$ and $x_{i+N^R}$ with $m \in \mathbb{N}^*$,

$$ s'(x_i) = c_{m,N^R}^i + a_{m,N^R}^i s'(x_{i+N^R}) + b_{m,N^R}^i s'(x_{i-m}). $$

with $c_{m,N^R}^i$, $a_{m,N^R}^i$ and $b_{m,N^R}^i$ coefficients depending on the chosen interpolation point $x_i$, the number of forward cells $N^R$, and the number of backward cells $m$.

For $m = 1$, the coefficients are those given by the previous formula with the coefficients (15). For $m = 2$, using equations (12) and (15), we have

$$ c_{2,N^R}^i = \frac{1}{1 - b_{1,N^R}^i \alpha_{i-1}} \left[ c_{1,N^R}^i + b_{1,N^R}^i \gamma_{i-1} \right], \quad a_{2,N^R}^i = \frac{a_{1,N^R}^i}{1 - b_{1,N^R}^i \alpha_{i-1}}, \quad b_{2,N^R}^i = \frac{b_{1,N^R}^i \beta_{i-1}}{1 - b_{1,N^R}^i \alpha_{i-1}}. \quad (16) $$

By recursion, we can establish that for $m > 2$,

$$ egin{cases} c_{m+1,N^R}^i = \frac{1}{1 - \alpha_{i-m} \frac{b_{m,N^R}^i}{b_{m-1,N^R}^i}} \left[ c_{m,N^R}^i + b_{m,N^R}^i \gamma_{i-m} - \alpha_{i-m} \frac{b_{m,N^R}^i}{b_{m-1,N^R}^i} c_{m-1,N^R}^i \right], \\ a_{m+1,N^R}^i = \frac{1}{1 - \alpha_{i-m} \frac{b_{m,N^R}^i}{b_{m-1,N^R}^i}} \left[ a_{m,N^R}^i - \alpha_{i-m} \frac{b_{m,N^R}^i}{b_{m-1,N^R}^i} a_{m-1,N^R}^i \right], \\ b_{m+1,N^R}^i = \frac{b_{m,N^R}^i \beta_{i-m}}{1 - \alpha_{i-m} \frac{b_{m,N^R}^i}{b_{m-1,N^R}^i}}. \end{cases} \quad (17) $$

With these two recursions (17) for $m = N^L$ and (15) for $n = N^R$, we have now a relation between further derivatives at the interpolating points $x_{i-N^L}$, $x_i$ and $x_{i+N^R}$ for a given number of cells on the right and on the left $N^R, N^L \in \mathbb{N}^*$ of $x_i$,

$$ s'(x_i) = c_{N^L,N^R}^i + a_{N^L,N^R}^i s'(x_{i+N^R}) + b_{N^L,N^R}^i s'(x_{i-N^L}). \quad (18) $$

This equation gives the exact derivatives of an equivalent global spline defined with the same interpolation points and on the same cells.
### 3.3.6 Independence of coefficients from the function

As mentioned previously the coefficients $\gamma_i$ in Equation (12) are a weighted sum of function values. The only coefficient in Equation (19) depending on the coefficients $\gamma_i$ is the coefficient $c_{N_I^L, N_I^R}^I$. It is also a weighted sum of function values. We can write it as follows (with the index $i$ referring to the interpolation points $x_i$ located on the interface $X_I$),

$$c_{N_I^L, N_I^R}^I = \sum_{k=-N^L}^{N^R} \omega_{k, N_I^L, N_I^R} f_{i+k}. \quad (21)$$

The weights $\{\omega_{k, N_I^L, N_I^R}\}_k$ only depend on the interpolation points. They can be determined at the initialization for each interface. In addition, the coefficients $a_{N_I^L, N_I^R}^I$ and $b_{N_I^L, N_I^R}^I$ in the relation (19) only depend on $\alpha_i$ and $\beta_i$ from the relation (12). They then depend only on the proportionality of the distances between the interpolation points, (i.e. rescaling the domain will not affect the values of $a_{N_I^L, N_I^R}^I$ and $b_{N_I^L, N_I^R}^I$).

The computation of the coefficients $a_{N_I^L, N_I^R}^I$ and $b_{N_I^L, N_I^R}^I$ and the weights $\{\omega_k\}_k$ needs only to be done at the initialization for each interface. It is also done pseudo-locally: only the information from two patches is needed to determine them.

### 3.3.7 Different boundary conditions

To establish the matrix system (20), we assumed that the derivatives on the external boundaries (i.e. the boundaries of the global domain) were known. In practice, we do not always know these derivatives. The two main examples are the case where we have defined enough interpolation points to not require additional closure conditions on the global domain, and the case where the global domain is periodic.

**Interpolation points as closure conditions.** If we define $N_c+3$ interpolation points, we do not need the derivatives, and we do not know them. All the degrees of freedom can be fixed. If we use Greville points as interpolation points, the previous figure, Fig. 4 becomes the Fig. 6 for the left boundary. To determine the cubic Hermite polynomials $P_-$, we use the value of the function at this additional interpolation point $x_*^-$ in Fig. 6.

<figure><img src="image_6.png" /><figcaption>Figure 6: Example of the case of interpolation points as closure conditions with two Hermite polynomials: $P_-$ defined on $[x_{i-1}^-, x_i^-]$ and $P_+$ defined on $[x_i^+, x_{i+1}^+]$ with functions connected across patches between $x_i^-$ and $x_i^+$. The left cell is a boundary cell of a global domain. We do not know the derivative at $x_{i-1}^-$. But, we know the function value at $x_*^-$ to provide closure for the system of equations.</figcaption></figure>

The Hermite polynomial $P_-$ has to interpolate the function at the additional interpolation point $x_*^-$, so we have the following relation

$$s(x_*^-) = s(x_{i-1}^-) H_0^{*-,-} + s(x_i^-) H_1^{*-,-} + s'(x_{i-1}^-) \Delta x_i^- K_0^{*-,-} + s'(x_i^-) \Delta x_i^- K_1^{*-,-}$$

with $H_0^{*-,-}$, $H_1^{*-,-}$, $K_0^{*-,-}$ and $K_1^{*-,-}$ the evaluations of $H_0$, $H_1$, $K_0$ and $K_1$ (11) at $\frac{x_*^- - x_{i-1}^-}{\Delta x_i^-}$.

As $s'(x_{i-1}^-)$ is unknown, we define it from the known values, and inject it into the relation 12 established before. This gives,

$$s'(x_i^-) = \frac{1}{1 + \beta_i \frac{K_1^{*-,-}}{K_0^{*-,-}}} \left[ \gamma_i + \frac{\beta_i}{\Delta x_i^- K_0^{*-,-}} s(x_*^-) - \frac{\beta_i H_0^{*-,-}}{\Delta x_i^- K_0^{*-,-}} s(x_{i-1}^-) - \frac{\beta_i H_1^{*-,-}}{\Delta x_i^- K_0^{*-,-}} s(x_i^-) \right] + \frac{\alpha_i}{1 + \beta_i \frac{K_1^{*-,-}}{K_0^{*-,-}}} s'(x_{i+1}^+)$$
So the relation on $[x_{i-1}^{-}, x_i^{-}]$ is

$$s'(x_i^{-}) = \gamma_i^{*,-} + \alpha_i^{*,-} s'(x_{i+1}^{+}). \qquad (22)$$

with

$$\left\{ \begin{aligned} \gamma_i^{*,-} &= \frac{1}{1 + \beta_i \frac{K_0^{*,-}}{K_0^{*,-}}} \left[ \gamma_i + \beta_i \frac{1}{\Delta x_i^{-} K_0^{*,-}} s(x_*^{-}) - \beta_i \frac{H_0^{*,-}}{\Delta x_i^{-} K_0^{*,-}} s(x_{i-1}^{-}) - \beta_i \frac{H_1^{*,-}}{\Delta x_i^{-} K_0^{*,-}} s(x_i^{-}) \right], \\ \alpha_i^{*,-} &= \frac{\alpha_i}{1 + \beta_i \frac{K_0^{*,-}}{K_0^{*,-}}}, \\ \beta_i^{*,-} &= 0. \end{aligned} \right.$$

If the additional Greville point is on $[x_i^{+}, x_{i+1}^{+}]$, the relation is similar

$$s'(x_i^{+}) = \gamma_i^{*,+} + \beta_i^{*,+} s'(x_{i-1}^{-}). \qquad (23)$$

with

$$\left\{ \begin{aligned} \gamma_i^{*,+} &= \frac{1}{1 + \alpha_i \frac{K_0^{*,+}}{K_1^{*,+}}} \left[ \gamma_i + \alpha_i \frac{1}{\Delta x_i^{+} K_1^{*,+}} s(x_*^{+}) - \alpha_i \frac{H_0^{*,+}}{\Delta x_i^{+} K_1^{*,+}} s(x_i^{+}) - \alpha_i \frac{H_1^{*,+}}{\Delta x_i^{+} K_1^{*,+}} s(x_{i+1}^{+}) \right], \\ \alpha_i^{*,+} &= 0, \\ \beta_i^{*,+} &= \frac{\beta_i}{1 + \alpha_i \frac{K_0^{*,+}}{K_1^{*,+}}}. \end{aligned} \right.$$

and $H_0^{*,+}$, $H_1^{*,+}$, $K_0^{*,+}$ and $K_1^{*,+}$ the evaluations of $H_0$, $H_1$, $K_0$ and $K_1$ (11) at $\frac{x_*^{+} - x_i^{+}}{\Delta x_i^{+}}$.

In the recursions (15) and (17), we use equations (22) and/or (23) instead of (12) for the last cells. This change affects only the first and the last interfaces for an equivalent global spline with interpolation points as closure conditions.

**Interpolation with periodic boundary conditions.** If the boundary conditions are periodic, we do not know the derivatives at the boundaries, but we can establish additional relations between the derivatives at $X_0$, $X_1$, $X_{N_p}$ and $X_{N_p-1}$ (see Fig. 7),

$$\left\{ \begin{aligned} s'(X_0) &= s'(X_{N_p}), \\ s'(X_0) &= c_{N_0^L, N_0^R}^0 + a_{N_0^L, N_0^R}^0 s'(X_1) + b_{N_0^L, N_0^R}^0 s'(X_{N_p-1}). \end{aligned} \right.$$

<figure><img src="image_7.png" /><figcaption>Figure 7: Domain split into $N_p$ patches. In the periodic case, the derivatives at $X_0$, $X_1$, $X_{N_p}$ and $X_{N_p-1}$ depend on each other. Moreover, the interfaces at $X_0$ and $X_{N_p-1}$ are physically overlapping.</figcaption></figure>
We add these relations to the matrix system (20). It adds another line, so the matrix is now of size $N_p \times N_p$.

$$ S = \begin{bmatrix} s'(X_0) \\ s'(X_1) \\ \vdots \\ s'(X_I) \\ \vdots \\ s'(X_{N_p-1}) \end{bmatrix} = \begin{bmatrix} 1 & -a_{N_0^L, N_0^R}^0 & & & & -b_{N_0^L, N_0^R}^0 \\ -b_{N_1^L, N_1^R}^1 & 1 & -a_{N_1^L, N_1^R}^1 & & & \\ & -b_{N_2^L, N_2^R}^2 & 1 & -a_{N_2^L, N_2^R}^2 & & \\ & & \ddots & & \ddots & \\ & & & -b_{N_{N_p-1}^L, N_{N_p-1}^R}^{N_p-1} & 1 & \end{bmatrix}^{-1} \begin{bmatrix} c_{N_0^L, N_0^R}^0 \\ c_{N_1^L, N_1^R}^1 \\ \vdots \\ c_{N_I^L, N_I^R}^I \\ \vdots \\ c_{N_{N_p-1}^L, N_{N_p-1}^R}^{N_p-1} \end{bmatrix} $$

### 3.3.8 Explicit formulation for uniform cells per patch

As mentioned before, for optimization reasons, it is more interesting to deal with uniform patches. To locally refine the global domain, some patches may also not have the same cell length. In this situation, we can also optimize the computation by giving an explicit expression for the coefficients $a_{N_I^L, N_I^R}^I, b_{N_I^L, N_I^R}^I$ and $c_{N_I^L, N_I^R}^I$, instead of the recursive relations (19), (20) given for the general case. The coefficients in (12) are simplified and are almost constant in the uniform per patch case. Let's take two uniform patches surrounding the interface $X_I$ at the interpolation point $x_i$. We introduce the following sequence $u$ such that

$$ u_k = (2 + \sqrt{3})^k - (2 - \sqrt{3})^k, \quad \forall k \in \mathbb{N}. $$

We can then prove by recursion that the coefficients in (18) can be explicitly expressed as follows with $\{u_k\}_{k \in \mathbb{N}}$, for $n, m \in \mathbb{N}^*$, (with $n$ representing $N_I^I$ and $m$ representing $N_I^L$ to simplify the notation),

$$ \left\{ \begin{aligned} a_{m,n}^I &= \frac{(-1)^{n-1} u_1 a_{1,1}^I u_m}{u_m u_n + u_m u_{n-1} a_{1,1}^I + u_n u_{m-1} b_{1,1}^I}, \\ b_{m,n}^I &= \frac{(-1)^{m-1} u_1 b_{1,1}^I u_n}{u_m u_n + u_m u_{n-1} a_{1,1}^I + u_n u_{m-1} b_{1,1}^I}, \\ c_{m,n}^I &= \sum_{k=-m}^{n} \omega_{k,m,n}^I f_{i+k}, \end{aligned} \right. \quad (24) $$

with the following weights,

$$ \left\{ \begin{aligned} \omega_{k,m,n}^I &= 3(-1)^k \frac{a_{1,1}^I u_m u_1}{u_m u_n + u_m u_{n-1} a_{1,1}^I + u_n u_{m-1} b_{1,1}^I}, & k = n, \\ \omega_{k,m,n}^I &= 3(-1)^k \frac{a_{1,1}^I u_m(u_{n-k+1} - u_{n-k-1})}{u_m u_n + u_m u_{n-1} a_{1,1}^I + u_n u_{m-1} b_{1,1}^I}, & k = 1, \dots, n-1, \\ \omega_{k,m,n}^I &= 3(-1)^k \frac{a_{1,1}^I u_m(u_n - u_{n-1}) - b_{1,1}^I u_n(u_m - u_{m-1})}{u_m u_n + u_m u_{n-1} a_{1,1}^I + u_n u_{m-1} b_{1,1}^I}, & k = 0, \\ \omega_{k,m,n}^I &= 3(-1)^{k+1} \frac{b_{1,1}^I u_n(u_{m+k+1} - u_{m+k-1})}{u_m u_n + u_m u_{n-1} a_{1,1}^I + u_n u_{m-1} b_{1,1}^I}, & k = -(m-1), \dots, -1, \\ \omega_{k,m,n}^I &= 3(-1)^{k+1} \frac{b_{1,1}^I u_n u_1}{u_m u_n + u_m u_{n-1} a_{1,1}^I + u_n u_{m-1} b_{1,1}^I}, & k = -m, \end{aligned} \right. \quad (25) $$

and $a_{1,1}^I = -\frac{1}{2} \frac{\Delta x^-}{\Delta x^+ + \Delta x^-}$ and $b_{1,1}^I = -\frac{1}{2} \frac{\Delta x^+}{\Delta x^+ + \Delta x^-}$. Some hints of the proof are given in Appendix B.
### 3.3.9 Approximation of the relation

Equation (20) gives derivatives which are exactly equal to the derivatives of an equivalent global spline. To compute them, we need the information at every interface. However, when the number of cells increases, we expect that the derivative at one interface will depend less on the derivatives at the other interfaces. We see here how the coefficients of Equation (20) behave when the number of cells becomes large. If they become negligible, the dependency then also becomes negligible.

In this section, we work with the numbers of cells $N_I^L$ and $N_I^R$ tending to infinity. Thanks to the symmetry of the problem, we can fix $N_I = N_I^R = N_I^L$ and make $N_I$ tend to infinity.

**Behavior of the coefficients $a_{N_I,N_I}^I$ and $b_{N_I,N_I}^I$.** These coefficients describe the dependency on the derivatives at the previous and next interfaces. Their absolute values for an increasing number of interpolation points are shown in Fig. 8.

In the uniform per patch case, the expressions in (24) give us an equivalence. Indeed, it is easy to see that $u_k = (2 + \sqrt{3})^k + o((2 + \sqrt{3})^k)$ when $k$ tends to infinity. So, injecting it in the equations, we find, for $n = N_I$ tending to infinity,

$$ \begin{cases} a_{n,n}^I \approx (-(2 - \sqrt{3}))^n (-4a_{11}^I) \\ b_{n,n}^I \approx (-(2 - \sqrt{3}))^n (-4b_{11}^I) \end{cases} \qquad (26) $$

These results agree with the tests shown in Fig. 8.

<figure><img src="image_8.png" /><figcaption>Figure 8: Absolute values of the coefficients $a_{N_I,N_I}^I$ and $b_{N_I,N_I}^I$ for the same number of cells $N_I$ in the two patches surrounding the interface $I$ with $n = N_I$ varying between 1 and 50. On the left, it is for uniform meshes. On the right several examples for non-uniform meshes generated with random functions. We can observe that for $N_I \simeq 28$, the values of these coefficients reach machine precision. The dependency on the derivatives at the other interfaces becomes negligible. On the left, we compare the coefficients to the equivalence relation (26).</figcaption></figure>

**Behavior of the coefficient $c_{N_I,N_I}^I$.** As mentioned before in (21), the coefficient $c_{N_I,N_I}^I$ can be expressed as a weighted sum of the function values: $c_{N_I,N_I}^I = \sum_{k=-N_I}^{N_I} \omega_{k,N_I,N_I} f_{i+k}$. It is more interesting to study the behavior of these weights than the behavior of the sum.

The absolute values of the weights for an increasing number of interpolation points are shown in the Fig. 9. The support of one patch is of length 1.

In the uniform per patch case, the expressions in (25) give us the following equivalences, for $n = N_I$
tending to infinity.

$$
\left\{
\begin{array}{ll}
\omega_{k,n,n} \approx (-1)^k u_1^2 \frac{a_{1,1}}{\Delta x^+} \left( (2-\sqrt{3})^k + (2-\sqrt{3})^{2n-k} \right), & k=1, \dots, n, \\
\\
\omega_{k,n,n} \approx u_1(1+u_1) \left( \frac{a_{1,1}^I}{\Delta x^+} - \frac{b_{1,1}^I}{\Delta x^-} \right), & k=0, \\
\\
\omega_{k,n,n} \approx (-1)^k u_1^2 \frac{b_{1,1}}{\Delta x^-} \left( (2-\sqrt{3})^{-k} + (2-\sqrt{3})^{2n+k} \right), & k=-n, \dots, -1,
\end{array}
\right.
\quad (27)
$$

These results agree with the tests shown in Fig. 9.

<figure><img src="image_9.png" /><figcaption>Figure 9: Absolute values of the weights $\omega_k$ in (21) for the same number of cells $N_I$ in the two patches surrounding the interface $I$ with $N_I$ varying between 1, 10, 20, 30 and 40. The weights for different $N_I$ are overlaid in the figures. The values of weights close to the interface also change with $N_I$. On the left, it is for uniform meshes. On the right, an example with non-uniform meshes generated with random functions. We can observe that for $N_I \simeq 30$, the values of these weights reach machine precision. The dependency on more distant function values becomes negligible. We can also note that for $N_I = 10$ in the uniform case, we retrieve the precision of $10^{-6}$ mentioned in [10] for Equation (9). On the left, the weights for $N_I = 40$ are compared to $\omega_{k,N_I,N_I} \simeq \omega_{\text{equiv},k,N_I,N_I} = 3N_I[(2-\sqrt{3})^{|k|} + (2-\sqrt{3})^{2N_I-|k|}]$, following the equivalence relation (27).</figcaption></figure>

**Consequences in the matrix system.** If we now look at the matrix system (20), we can see that a lot of the terms are negligible if the number of cells becomes large enough.

$$
\mathbf{S} = \begin{bmatrix} s'(X_1) \\ \vdots \\ s'(X_I) \\ \vdots \\ s'(X_{N_p-1}) \end{bmatrix} = \begin{bmatrix} 1 & -\frac{a_{N_I^L, N_I^R}}{N_I^L} & \cdots & -\frac{a_{N_I^L, N_I^R}}{N_I^L} \\ \frac{b_{N_I^L, N_I^R}}{N_I^L} & 1 & \cdots & \frac{b_{N_I^L, N_I^R}}{N_I^L} \\ \vdots & \vdots & \ddots & \vdots \\ -\frac{b_{N_p-1, N_I^R}}{N_p-1} & \frac{a_{N_p-1, N_I^R}}{N_p-1} & \cdots & 1 \end{bmatrix}^{-1} \begin{bmatrix} c_{N_I^L, N_I^R}^1 + \frac{b_{N_I^L, N_I^R}^1}{N_I^L} s'(X_0) & \cdots & c_{N_I^L, N_I^R}^I \\ \vdots & \ddots & \vdots \\ c_{N_p-1, N_I^R}^1 + \frac{a_{N_p-1, N_I^R}}{N_p-1} s'(X_p) & \cdots & c_{N_p-1, N_I^R}^I \end{bmatrix} \quad \mathbf{S} \to \mathbf{C}_{\text{trunc}}, \quad \text{when } N_I^L, N_I^R \to \infty, \forall I \quad (28)
$$

with $\mathbf{C}_{\text{trunc}}$ the vector of truncated sum,

$$
s'(X_I) \rightarrow \sum_{k=-N_{\text{trunc}}^L}^{N_{\text{trunc}}^R} \omega_k f_{i+k}.
$$

We observe the same behavior for the different boundary conditions mentioned earlier. We can choose different truncation values depending on the desired precision. We refer to Table 1 to select the precision.
<table><caption>Table 1: Different precisions according to the number of cells we take into account. The values in this table are expressed to three significant figures. If we use the approximation, we will have an error of the order of $|a_{N_I,N_I}^I||s'(X_{I+1})|$ because we remove the dependency on the other derivatives. We will also have an error of the order of $|\omega_{N_I,N_I,N_I}||f_{I+N_I}|$, because we truncate the sum in (24). For instance, for a cell length of $10^{-1}$, the error made with the approximation will be of order $(|f_{I+N_I}| + 10^{-2}|s'(X_{I+1})|)|a_{1,1}^I| \cdot 10^{-15}$.</caption><thead><tr><th>N<sub>I,trunc</sub></th><th>5</th><th>10</th><th>15</th><th>20</th><th>25</th><th>30</th></tr></thead><tbody><tr><td>(2 - √3)<sup>N<sub>I</sub></sup></td><td>1.38 · 10<sup>-3</sup></td><td>1.91 · 10<sup>-6</sup></td><td>2.64 · 10<sup>-9</sup></td><td>3.64 · 10<sup>-12</sup></td><td>5.03 · 10<sup>-15</sup></td><td>6.94 · 10<sup>-18</sup></td></tr><tr><td>|a<sup>I</sup><sub>N<sub>I</sub>,N<sub>I</sub></sub>| / |a<sup>I</sup><sub>1,1</sub>|</td><td>5.52 · 10<sup>-3</sup></td><td>7.63 · 10<sup>-6</sup></td><td>1.05 · 10<sup>-8</sup></td><td>1.46 · 10<sup>-11</sup></td><td>2.01 · 10<sup>-14</sup></td><td>2.78 · 10<sup>-17</sup></td></tr><tr><td>|b<sup>I</sup><sub>N<sub>I</sub>,N<sub>I</sub></sub>| / |b<sup>I</sup><sub>1,1</sub>|</td><td>5.52 · 10<sup>-3</sup></td><td>7.63 · 10<sup>-6</sup></td><td>1.05 · 10<sup>-8</sup></td><td>1.46 · 10<sup>-11</sup></td><td>2.01 · 10<sup>-14</sup></td><td>2.78 · 10<sup>-17</sup></td></tr><tr><td>|ω<sub>N<sub>I</sub>,N<sub>I</sub>,N<sub>I</sub></sub>| / |Δx<sup>+</sup>| / |a<sup>I</sup><sub>1,1</sub>|</td><td>1.66 · 10<sup>-2</sup></td><td>2.29 · 10<sup>-5</sup></td><td>3.16 · 10<sup>-8</sup></td><td>4.37 · 10<sup>-11</sup></td><td>6.03 · 10<sup>-14</sup></td><td>8.33 · 10<sup>-17</sup></td></tr></tbody></table>

## 3.4 Semi-Lagrangian advection on a multi-patch domain

**Application of a multi-patch interpolation in BSL.** Now that we know how to interpolate on a multi-patch domain, let's detail the backward semi-Lagrangian algorithm on a multi-patch domain.

1. We first solve the equations of the characteristics to get the feet of the characteristics. The feet may not land on the same patch as the starting point (i.e. their initial condition in the equation of the characteristics).

2. We then identify on which patch the feet are physically located. Different algorithms exist. If the mappings $\mathcal{F}$ from each of the logical domains to the physical domain are analytically invertible, we can identify on which subdomain a physical coordinate can be found. For special geometries such as the concentric rings implemented in our simulations, we can also identify on which patch a coordinate is physically located by computing the distance to the O-point.

3. The function values have been computed at the previous time step. We can now compute the derivatives with the method described in the previous section, (20) from the function values.

4. The function values and the derivatives at the interfaces are then known. We have all the necessary information to build a local spline representation of the function at $t - \Delta t$ on each patch. We interpolate the corresponding local spline at the foot of each characteristic. These values are the values of the function at $t$ at the interpolation points.

# 4 Semi-Lagrangian scheme on a 2D multi-patch domain

Now, we know how to handle 1D multi-patch domains. We may also want to handle 2D multi-patch domains. Indeed, in Fig. 1 we have a 2D multi-patch case.

In this case, the global domain is now split into multiple patches along two dimensions. In the guiding-center model (2), it is the $r$ direction and the $\theta$ direction. Local refinement in a 2D multi-patch domain can be split into two cases: a conforming global mesh or a non-conforming global mesh. A mesh is called non-conforming when it is not the result of a tensor product of two 1D meshes. In multi-patch geometries, this means that there are break points defined on one patch with no counterpart in its neighboring patch. Moreover, the 2D multi-patch case leads to another case which must be studied: T-joints. T-joints designate patch connections involving more than two patches such that only part of an edge of a patch is connected to the edge of another patch. In this section, we present the three cases.

## 4.1 Conforming mesh: a direct application of the 1D case

In the conforming case, the 2D B-spline basis of an equivalent global domain is a tensor product of two 1D B-spline bases of equivalent global 1D domains. The computation of the interface derivatives is then a direct application of the 1D formulae, (19), (20). In the 2D multi-patch case, we apply Hermite boundary conditions in both directions. So, we need to define the derivatives along $r$ and $\theta$ on the edges
and the cross-derivatives on the corners. Let's note $N_r$ the number of cells in the $r$-direction and $N_\theta$ the number of cells in the $\theta$-direction.

$$ 
\left\{
\begin{array}{ll}
\partial_r s(r_i, \theta_j) = \sum_{k=0}^{N_r+2} \sum_{l=0}^{N_\theta+2} c_k c_l b'_k(r_i) b_l(\theta_j), & i=0, N_r, \quad j=0, \dots, N_\theta, \\
\\
\partial_\theta s(r_i, \theta_j) = \sum_{k=0}^{N_r+2} \sum_{l=0}^{N_\theta+2} c_k c_l b_k(r_i) b'_l(\theta_j), & i=0, \dots, N_r, \quad j=0, N_\theta, \\
\\
\partial_{r\theta} s(r_i, \theta_j) = \sum_{k=0}^{N_r+2} \sum_{l=0}^{N_\theta+2} c_k c_l b'_k(r_i) b'_l(\theta_j), & i=0, N_r, \quad j=0, N_\theta.
\end{array}
\right.
$$

For a fixed $j \in [0, N_\theta]$, $\forall i \in [0, N_r]$, we can write the 2D-spline $(r, \theta) \mapsto s(r, \theta)$ as a 1D-spline $r \mapsto s_{\theta_j}(r)$ and its derivatives as,

$$ \partial_r s(r_i, \theta_j) = \left( \sum_{k=0}^{N_r+2} c_k b'_k(r_i) \right) \left( \sum_{l=0}^{N_\theta+2} c_l b_l(\theta_j) \right) = s'_{\theta_j}(r_i) = \gamma_{\theta_j, i} + \alpha_i s'_{\theta_j}(r_{i+1}) + \beta_i s'_{\theta_j}(r_{i-1}), $$

by applying the relation (12) established in the previous section to the function $r \mapsto s_{\theta_j}(r)$. The subscript $j, i$ is used to highlight that here $\gamma_{\theta_j, i}$ is a weighted sum of $s_{\theta_j}(r_{i-1})$, $s_{\theta_j}(r_i)$ and $s_{\theta_j}(r_{i+1})$. Similarly, with the other direction, we then have the following relation,

$$ \partial_{\theta} s(r_i, \theta_j) = \gamma_{r_i, j} + \alpha_j s'_{r_i}(\theta_{j+1}) + \beta_j s'_{r_i}(\theta_{j-1}). $$

By fixing $\theta_j$ or $r_i$, we can proceed similarly for the cross-derivatives. Let's take for example $\theta_j$ fixed:

$$ \partial_{r\theta} s(r_i, \theta_j) = \left( \sum_{k=0}^{N_r+2} c_k b'_k(r_i) \right) \left( \sum_{l=0}^{N_\theta+2} c_l b'_l(\theta_j) \right) = s'_{\theta'_j}(r_i) = \gamma_{\theta'_j, i} + \alpha_i s'_{\theta'_j}(r_{i+1}) + \beta_i s'_{\theta'_j}(r_{i-1}). $$

We then have equivalent formulae for the derivatives along $r$ and $\theta$ and the cross-derivatives. With the method applied in the previous section 3, we can obtain equations linking the derivatives at the interfaces. Let's denote by $\{R_I\}_I$ the interfaces in the $r$-direction and $\{\Theta_J\}_J$ the interfaces in the $\theta$-direction. We note $N_p^R$ the number of patches in the $r$ direction and $N_p^\Theta$ the number of patches in the $\theta$ direction. The equations linking the derivatives are:

$$ 
\left\{
\begin{array}{ll}
\partial_r s(R_I, \theta_j) = c_{N_I^L, N_I^R}^I + a_{N_I^L, N_I^R}^I \partial_r s(R_{I+1}, \theta_j) + b_{N_I^L, N_I^R}^I \partial_r s(R_{I-1}, \theta_j), & \forall j \in [0, N_\theta], I \in [1, N_p^R], \\
\\
\partial_\theta s(r_i, \Theta_J) = c_{N_J^L, N_J^R}^J + a_{N_J^L, N_J^R}^J \partial_\theta s(r_i, \Theta_{J+1}) + b_{N_J^L, N_J^R}^J \partial_\theta s(r_i, \Theta_{J-1}), & \forall i \in [0, N_r], J \in [1, N_p^\Theta], \\
\\
\partial_{r\theta} s(R_I, \Theta_J) = c_{N_I^L, N_I^R}^{I,*} + a_{N_I^L, N_I^R}^I \partial_{r\theta} s(R_{I+1}, \Theta_J) + b_{N_I^L, N_I^R}^I \partial_{r\theta} s(R_{I-1}, \Theta_J), & I \in [1, N_p^R], J \in [1, N_p^\Theta],
\end{array}
\right.
\quad (29) $$

with the coefficients $a_{N_I^L, N_I^R}^I$, $b_{N_I^L, N_I^R}^I$ and $c_{N_I^L, N_I^R}^I$ calculated with the same formula as in (19) for $N_I^R$ cells along the $r$-direction in the right patch and $N_I^L$ cells along the $r$-direction in the left patch. We proceed similarly for the $a_{N_J^L, N_J^R}^J$, $b_{N_J^L, N_J^R}^J$ and $c_{N_J^L, N_J^R}^J$ coefficients in the $\theta$-direction. The $c_{N_I^L, N_I^R}^I$ coefficient is a weighted sum of the function values $\{f_{i+k,j}\}_{k=-N_I^L, \dots, N_I^R}$ in the $r$-direction for a fixed $j$ and $i$ corresponding to the interpolation point located on the interface $R_I$. In the equation describing the cross-derivatives, $c_{N_I^L, N_I^R}^{I,*}$ is a weighted sum of the function derivatives $\{\partial_\theta f_{i+k,j}\}_{k=-N_I^L, \dots, N_I^R}$ along the $\theta$-direction in the $r$-direction for a fixed $j$ and $i$ corresponding to the interpolation point located on the interface $R_I$.

**Remark.** A remark concerning the given relation between the cross-derivatives is that we can also establish it in the other direction,

$$ \partial_{r\theta} s(R_I, \Theta_J) = c_{N_J^L, N_J^R}^{J,*} + a_{N_J^L, N_J^R}^J \partial_{r\theta} s(R_I, \Theta_{J+1}) + b_{N_J^L, N_J^R}^J \partial_{r\theta} s(R_I, \Theta_{J-1}), \quad I \in [1, N_p^R], J \in [1, N_p^\Theta], $$

where $c_{n_J, m_J}^{J,*}$ is the weighted sum of the function derivatives $\{\partial_r f_{i,j+k}\}_{k=-N_J^L, \dots, N_J^R}$ along the $r$-direction in the $\theta$-direction for a fixed $i$ and $j$ corresponding to the interpolation point located on the interface $\Theta_J$. Both formulae are equivalent and the cross-derivatives are equal to the cross-derivatives of an equivalent 2D global spline.
**Algorithm.** As the equations describing the cross-derivatives depend on the derivatives of the function at the interfaces, to compute all the derivatives in the conforming case, the algorithm needs only two steps:

1. From the function values at the interpolation points, compute the derivatives along $r$ and $\theta$ with (29).

2. From the computed derivatives at the interfaces, compute the cross-derivatives at the corners by choosing one of the two directions with (29).

## 4.2 Non-Conforming mesh

Non-conforming multi-patch meshes refer to meshes where break points defined on one patch have no counterpart on its neighboring patch. The global mesh is not the result of a tensor product of two 1D global meshes. So, there is no equivalent global spline. The mesh local to a given patch remains conforming and the result of a tensor product of two 1D local meshes. So we can still define local splines. We suppose here that the non-conforming mesh is not too degenerated. By this we mean that for an interface between two patches, the one-dimensional mesh of one of the two patches is included in the one-dimensional mesh of the other patch.

We cannot apply formulae (29) along all the lines of each direction because they may end at some interfaces. The idea is then to progressively compute the derivatives and the cross-derivatives. As mentioned before, the cross-derivatives can be computed from the derivatives along one direction. From the derivatives and the cross-derivatives, we can build spline representations of the derivatives. If we want to impose a $C^1$ regularity between the patches as mentioned previously, we build a spline on a coarser patch of an interface. We interpolate it at the interpolation points on the refined patch which have no counterpart on the coarser patch. The interpolated values are then assigned to the refined patch to enforce the $C^1$ regularity.

In some cases, we can define an algorithm to compute all the derivatives and cross-derivatives. For other cases, we may be limited and will not use all the information of the patches. For these cases, it is recommended to apply an approximation of the exact formulae (28), if the number of cells per patch is big enough.

To illustrate this point, here are two examples of non-conforming patches. Fig. 10 shows an example of a situation where we can compute all the derivatives and cross-derivatives using all the available information. Fig. 11 shows an example of a situation where we cannot use all the available information if we want to apply the exact formulae (29).
<figure><img src="image_10.png" /></figure>

(a) Compute the derivatives (red lines) along the lines crossing the whole domain in $r$ and $\theta$ directions, the "conforming" lines with (29). From these derivatives, we compute the cross-derivatives ($\times$) on the coarser patches. With the derivatives along $r$ and the cross-derivatives, we can build spline representations of the $r$-derivatives on the coarse patches. We do the same for the $\theta$-derivatives. We interpolate the local splines at the interpolation points of the refined patches with no counterpart in the coarse patches ($\bullet$).

(b) With the derivatives computed at the previous step at the "non-conforming" interpolation points ($\bullet$) and the boundary conditions, we can compute the derivatives along the red lines with (29). From these derivatives, we compute the cross-derivatives at $\times$. As previously, we build a spline representation of the derivatives along $r$ and $\theta$ using the newly computed derivatives and cross-derivatives. The splines have to be interpolated at the $\bullet$ points.

(c) With the derivatives computed at the previous step and the boundary conditions, we can compute the derivatives along the red lines using (29) again. From them, we compute the cross-derivatives at $\times$. The function values are known at every interpolation point and all the derivatives are now determined. We have all the necessary information to build local spline representations of the function on the whole domain.

<caption>Figure 10: Example of an algorithm for non-conforming patches where we can successively compute all the derivatives and cross-derivatives.</caption>
<figure><img src="image_11.png" /></figure>

(a) As in the first example in Fig. 10, we start by computing the derivatives along the "conforming" lines in red on the figure. Then, we can compute the cross-derivatives at the ×. However, by doing that we do not use the fact that patches 15 and 11 are refined. If we build a spline representation of the derivative on the interface 15/11, it will be as coarse as if we build it in the patch 16 (for example). So, we lose information from the patches 15 and 11.

(b) If the option in the sub-figure (a) is a bad idea, we can then think that we can start by computing the cross-derivatives at the × on this sub-figure. However, we will encounter the same problem at the interface 2/6. Because of the dependency of the derivatives at all the interfaces in (29), it will affect the cross-derivatives used to define a spline at the interface 3/7.

(c) The case of the two first subfigures can also be encountered if we compute the cross-derivatives in the θ-direction. This illustrates that if we apply the exact formulae in (29), we may not always use all the available information. A solution is to add enough cells in the problematic patches in order to use an approximation of the formulae (see (28)). Instead of computing the derivatives at ⋄ by interpolating a spline representation of the derivative, we approximate the derivatives with (28).

Figure 11: Example of non-conforming patches where the application of the exact formulae (29) may not use all the available information. In this case, the application of the approximated formulae (28) seems to be a better solution if the number of the cells in the patches is big enough.

In practice, the first case of non-conforming meshes may appear more often. It is for example the case in Fig. 1. The second case was described in this article mainly as a warning.

**Remark.** As mentioned previously, we impose a $C^1$ regularity at the interfaces. It is imposed by building spline representations of the derivatives and evaluating them at the interpolation points with no counterpart in the neighboring patches. Thus, contrary to the conforming case, we lose the $C^2$ regularity that we had with the Hermite polynomials in Section 3.3.3.

## 4.3 Geometries with T-joints

### 4.3.1 Definition of a T-joint

We define a *T-joint* as a connection between at least three patches, where an edge of one patch is not fully connected to the edge of another patch along its entire length. An illustration between three patches is given in Fig. 12. The introduction of T-joints may be necessary for some complex geometries involving an O-point and an X-point. This can be observed in Fig. 1. We will analyze a case with T-joints on a domain with an O-point.
<figure><img src="image_12.png" /><figcaption>Figure 12: Example of a T-joint between three patches. The mapping between the logical domains and the physical domain is the identity. We can see that part of an edge of patch 1 is connected to an edge of patch 2 on a subdomain and to an edge of patch 3 on a complementary subdomain. This type of connection is called a *T-joint*.</figcaption></figure>

### 4.3.2 Application of the formula to the T-joint cases

The computation of the derivatives on geometries with T-joints stays pretty similar to the cases mentioned before. For conforming meshes, we remark that the cross-derivatives can be computed along one of the directions without a difference in the result. On geometries with T-joints, the choice of the direction may be influenced by the geometry. For instance, in the example of Fig. 12, the cross-derivatives will be computed along the *r*-direction. Some *r*-derivatives on patch 1, needed to apply equation (29) along the *θ*-direction, are not available.

**Example applied later in our simulation.** Let's study the geometry in Fig. 13, which will be an application case for the simulations presented in the Section 5 of this article. The global mesh is conforming, so we can build an equivalent global spline on this domain. To compute all the derivatives, we start by computing the derivatives along *r* on each patch using (29). We then have all the necessary information to build spline representations on patches 1 and 5, but not on patches 2, 3 and 4. The derivatives still need to be computed along *θ* on patches 2, 3 and 4. From the derivatives along *r*, we compute the cross-derivatives at the T-joints. All the required derivatives are then computed, and we can build spline representations on patches 2, 3 and 4.

<figure><img src="image_13.png" /><figcaption>Figure 13: Example of the geometry with T-joints and an O-point studied for our simulations in the Section 5. On the left, the patches in the logical domains. On the right, the physical domain. The global mesh is a conforming mesh.</figcaption></figure>
# 5 Numerical results

Here, we present some numerical results of the method presented in the previous sections. The implementation can be found in the Gyselalib++ library [2]. This section is organized as follows. The first results are spline interpolation results on three concentric patches. We check that the approximation formula returns the expected precision. Secondly, we test the interpolation in an advection simulation with a backward semi-Lagrangian method. Finally, results corresponding to more complete simulation on a guiding-center model, (2), involving a coupling with a Poisson equation are shown.

## 5.1 Spline interpolation tests

We compare the results of an interpolation on local splines built with the derivatives computed using the approximation relations and with the exact formula to an equivalent global spline. We test it with the following function

$$f : \hat{\Omega} \to \mathbb{R} : (r, \theta) \mapsto \cos(2\pi x(r, \theta)) \sin(2\pi y(r, \theta)),$$

with $x$ and $y$ given by the Czarny mapping [11],

$$\mathcal{F} : \hat{\Omega} \to \Omega : (r, \theta) \mapsto \begin{cases} x(r, \theta) = \frac{1}{\varepsilon}(1 - \sqrt{1 + \varepsilon(\varepsilon + 2r \cos(\theta))}), \\ y(r, \theta) = \frac{e\xi r \sin(\theta)}{2 - \sqrt{1 + \varepsilon(\varepsilon + 2r \cos(\theta))}}, \end{cases} \quad (30)$$

with $\xi = \frac{1}{\sqrt{1-\frac{1}{4}\varepsilon^2}}$, and $\varepsilon$ and $e$ given constants, that we will fix to $\varepsilon = 0.3$ and $e = 1.4$. The interpolation is carried out on the logical domain. The physical domain is not involved. We introduced a mapping in the test function to have curves which do not follow the mesh lines, as will be the case in the next simulations.

The global domain is split into three patches. The splitting of the global domain is only along the $r$-direction. It is a 1D multi-patch case as only the $r$-grid is split. For the patch connections, we only need to compute the $r$-derivatives. We test two cases: a uniform global mesh and a non-uniform global mesh. In the non-uniform global mesh case, the patches are uniform, but the cell length differs from one patch to another. In the following tests, the number of cells in each patch stays the same in both cases but the position of the interfaces changes. The grid parameters of the patches and the global domain are given in Table 2.

<table><caption>Table 2: Number of cells and domains of the patches and equivalent global domain for the spline interpolation tests.</caption><thead><tr><th>Test</th><th>Patch</th><th>Number of cells</th><th>Domain along <i>r</i></th><th><i>dr</i></th><th>Domain along <i>θ</i></th></tr></thead><tbody><tr><td rowspan="4">1.1 Uniform</td><td>Patch1</td><td>[42 × 256]</td><td>[0, 42/128]</td><td>7.8125 · 10<sup>-3</sup></td><td>[0, 2π]</td></tr><tr><td>Patch2</td><td>[44 × 256]</td><td>[42/128, 86/128]</td><td>7.8125 · 10<sup>-3</sup></td><td>[0, 2π]</td></tr><tr><td>Patch3</td><td>[42 × 256]</td><td>[86/128, 1]</td><td>7.8125 · 10<sup>-3</sup></td><td>[0, 2π]</td></tr><tr><td>Global</td><td>[128 × 256]</td><td>[0, 1]</td><td>7.8125 · 10<sup>-3</sup></td><td>[0, 2π]</td></tr><tr><td rowspan="4">1.2 Non-uniform</td><td>Patch1</td><td>[42 × 256]</td><td>[0, 42/128]</td><td>7.8125 · 10<sup>-3</sup></td><td>[0, 2π]</td></tr><tr><td>Patch2</td><td>[44 × 256]</td><td>[42/128, 60/128]</td><td>3.1960 · 10<sup>-3</sup></td><td>[0, 2π]</td></tr><tr><td>Patch3</td><td>[42 × 256]</td><td>[60/128, 1]</td><td>1.2649 · 10<sup>-2</sup></td><td>[0, 2π]</td></tr><tr><td>Global</td><td>[128 × 256]</td><td>[0, 1]</td><td></td><td>[0, 2π]</td></tr></tbody></table>

The results for different numbers of cells for the approximation (28) and the exact formula (19) are shown in Table 3. We compare the values of the interface derivatives computed with the formulae to the derivatives of an equivalent global spline. We also compare the values of the local and global spline representations on the interpolation points. The given errors are the $L_\infty$-norm of the difference,

$$err_{\text{funct}} = \|s_{\text{local}} - s_{\text{global}}\|_{L_{\infty}(\hat{\Omega})}, \quad err_{\text{deriv}} = \|\partial_r s_{\text{local}} - \partial_r s_{\text{global}}\|_{L_{\infty}(\partial\hat{\Omega}^k)}.$$

We observe that the results in Table 3 agree with the estimations in Table 1. Moreover, as expected, for thirty cells, the computed interface derivative values reach the same precision as if we use all the cells with the exact formula.
<table><thead><tr><th>Number of cells</th><th>Errors on local splines</th><th>Errors on interfaces derivatives</th></tr></thead><tbody><tr><td>5</td><td>5.15e-05</td><td>4.16e-02</td></tr><tr><td>10</td><td>6.37e-08</td><td>5.15e-05</td></tr><tr><td>20</td><td>7.86e-14</td><td>6.35e-11</td></tr><tr><td>30</td><td>8.44e-15</td><td>2.37e-13</td></tr><tr><td>40</td><td>8.44e-15</td><td>2.36e-13</td></tr><tr><td>Exact</td><td>8.44e-15</td><td>2.82e-13</td></tr></tbody></table>

(a) Uniform mesh between the three patches.

<table><thead><tr><th>Number of cells</th><th>Errors on local splines</th><th>Errors on interfaces derivatives</th></tr></thead><tbody><tr><td>5</td><td>4.61e-05</td><td>2.57e-02</td></tr><tr><td>10</td><td>6.08e-08</td><td>3.23e-05</td></tr><tr><td>20</td><td>9.55e-14</td><td>4.77e-11</td></tr><tr><td>30</td><td>8.44e-15</td><td>2.59e-13</td></tr><tr><td>40</td><td>8.44e-15</td><td>3.84e-13</td></tr><tr><td>Exact</td><td>8.44e-15</td><td>3.36e-13</td></tr></tbody></table>

(b) Non-uniform mesh between the three patches.

<caption>Table 3: Interpolation results. The given errors are the difference of the results on local splines defined on patch 1, patch 2 and patch 3 to an equivalent global spline. The indicated number of cells corresponds to the number of cells taken in the approximation of the derivatives. We see that indeed for 30, it reaches the same precision as if we take all the cells with the exact formula. The results agree with the analytical ones given in Table 1 until saturation at $10^{-13}$ (for a space step around $10^{-2}$).</caption>

## 5.2 Advection tests

We test the method on an advection with a test case similar to the one detailed in [37]. The exact formula (20) is applied, so we expect the same results between the local splines and the equivalent global spline. The initial condition on the global domain is given by,

$$ 
\left\{
\begin{array}{l}
\rho(x, y) = \frac{1}{2} [G(r_1(x, y)) + G(r_2(x, y))], \quad (x, y) \in \Omega^2 \\
r_1(x, y) = \sqrt{(x - x_0)^2 + (y - y_0)^2}, \\
r_2(x, y) = \sqrt{8(x - x_0)^2 + (y - y_0)^2}, \\
G(r) = \cos\left(\frac{\pi r}{2a}\right) 1_{r \le a},
\end{array}
\right.
\qquad (31) $$

with $(x_0, y_0) = \mathcal{F}(0.5, 0)$ and $a = 0.3$. The physical domain is defined by the Czarny mapping $\mathcal{F}$ (30), [11]. This function is advected according to the following advection field $\mathbf{A}$,

$$ \mathbf{A} = \mathbf{J}_{\mathcal{F}} \begin{bmatrix} 0 \\ \omega \end{bmatrix}, \qquad (32) $$

with $\mathbf{J}_{\mathcal{F}}$ the Jacobian matrix of the mapping $\mathcal{F}$ (30), and $\omega = 2\pi$.

The trajectory is a rotation following the mesh lines. The results in [37] are given for a uniform global domain defined on $N_r \times N_\theta = 128 \times 256$. The chosen time step is $\Delta t = 0.01$. To find the feet of the characteristics, we apply a third order Runge-Kutta method.

In our tests, we split the domain into three (concentric rings) for (Test 2.1) or five patches (see Fig. 13) for (Test 2.2). We apply the exact formulae (20) in the first test and (29) in the second. The domain discretizations of these two tests are detailed in Table 4.

<table><caption>Table 4: Number of cells and domains of the patches and equivalent global domains for the advection tests. In parentheses, we indicate the dimensions of the multi-patch mesh. If the multi-patch mesh is 1D, we only need to compute the derivatives along $r$. If it is 2D, we compute the derivatives along $r$, along $\theta$ and the cross-derivatives.</caption><thead><tr><th>Test</th><th>Patch</th><th>Number of cells</th><th>Domain along r</th><th>Domain along θ</th></tr></thead><tbody><tr><td rowspan="4">2.1 Non-conforming (1D)</td><td>Patch1</td><td>[21 × 128]</td><td>[0, 42/128]</td><td>[0, 2π]</td></tr><tr><td>Patch2</td><td>[76 × 256]</td><td>[42/128, 80/128]</td><td>[0, 2π]</td></tr><tr><td>Patch3</td><td>[24 × 128]</td><td>[80/128, 1]</td><td>[0, 2π]</td></tr><tr><td>Refined global</td><td>[128 × 256]</td><td>[0, 1]</td><td>[0, 2π]</td></tr><tr><td>Coarse global</td><td>[128 × 128]</td><td>[0, 1]</td><td>[0, 2π]</td></tr><tr><td rowspan="6">2.2 T-joints (2D)</td><td>Patch1</td><td>[42 × 255]</td><td>[0, 42/128]</td><td>[0, 2π]</td></tr><tr><td>Patch2</td><td>[38 × 85]</td><td>[42/128, 80/128]</td><td>[0, 2π/3]</td></tr><tr><td>Patch3</td><td>[38 × 85]</td><td>[42/128, 80/128]</td><td>[2π/3, 4π/3]</td></tr><tr><td>Patch4</td><td>[38 × 85]</td><td>[42/128, 80/128]</td><td>[4π/3, 2π]</td></tr><tr><td>Patch5</td><td>[48 × 255]</td><td>[80/128, 1]</td><td>[0, 2π]</td></tr><tr><td>Global</td><td>[128 × 255]</td><td>[0, 1]</td><td>[0, 2π]</td></tr></tbody></table>
**Test 2.1 Non-conforming patches.** Patch 1 and Patch 3 are made coarser along the $\theta$ direction. The global mesh is then no longer conforming. After computing the derivatives along $r$ which are on the conforming lines, we build $\theta$-spline representations of the $r$-derivatives to compute the $r$-derivatives on the non-conforming lines at the boundaries of Patch 2. Fig. 14 and Fig. 15 show the results. Here, there is no equivalent global spline. So, we cannot expect a machine precision agreement to a global spline. However, we can build a refined global spline and a coarse global spline and compare the three spline representations to the exact solution. We expect the error with respect to the exact solution of the local splines to be bounded between the error of the refined global spline and the coarse global spline. This is what we observe in Fig. 15. The results of this figure also highlight the efficiency of the non-conforming mesh: despite using fewer cells, the local splines tend to be as precise as the refined spline. The non-conforming mesh is, of course, adapted to the solution trajectory presented in Fig. 14. The errors we obtain in Fig. 15 for $t = 0.8$ agree with the error obtained in [37], see Table 2.

<figure><img src="image_14.png" /><figcaption>Figure 14: Advection test 2.1: non-conforming patches. There is no comparison to an equivalent global spline as it does not exist, but we compare here the errors to the exact solution with the errors of a refined spline and a coarse spline. Errors of the spline representations to the exact solution at $t = 2.0$. The dashed line shows the trajectory of the maximum density point.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 15: Advection test 2.1: non-conforming patches. Time variation of the errors of the spline representations with respect to the exact solution at each time steps. The errors of the local splines are bounded by the errors of the refined spline and the errors of the coarse spline as expected. We also observe that the local splines on the non-conforming mesh tend to be as good as the spline on the refined mesh.</figcaption></figure>
**Test 2.2 T-joints.** We apply the method to a case with T-joints. Here the multi-patch splitting is 2D, so we need to compute the r-derivatives, the θ-derivatives and the cross-derivatives as explained in paragraph 4.3.2 for the studied geometry described in Fig. 13. We choose a uniform global mesh. So, we expect a difference between the final local spline representations and the final equivalent global spline representation at machine precision. Fig. 16 shows that the spline representations agree. We obtain an error of 6 · 10⁻¹³ due to machine precision.

<figure><img src="image_16.png" /><figcaption>Figure 16: Advection test 2.3: T-joints. Difference between the results on multiple patches and the results on an equivalent global spline. The dashed line shows the trajectory of the maximum density point. We can see that at t = 2.0 (200 iterations), ||s<sub>local</sub> − s<sub>global</sub>||<sub>L∞</sub> = 5.978 · 10⁻¹³. The two spline representations agree even when the trajectory crosses the interfaces.</figcaption></figure>

## 5.3 Guiding-center model - diocotron instability simulation

The diocotron instability simulation is a reference numerical test for the guiding-center equations (2). Physically, it simulates a non-neutral plasma in a cylinder, whose particles are confined radially by a uniform axial magnetic field with a cylindrical conducting wall located at the outer boundary. The simulation is detailed in [12, 13, 29, 30, 36, 37].

The electric potential φ is null at the boundaries. The initial density profile is given by

$$ \rho(t=0, r, \theta) = \rho_0(r) + \rho_1(t=0, r, \theta) = (1 + \varepsilon \cos(m\theta)) \exp\left[-\left(\frac{r-\bar{r}}{d}\right)^p\right] \mathbf{1}(r^{-} \le r \le r^{+}), \quad (33) $$

with ρ₀ the equilibrium and ρ₁, the perturbation of amplitude ε and mode m. The initial density profile is smoothed by an exponential to avoid discontinuities.

The functions are defined on $\hat{\Omega} = [0, 1] \times [0, 2\pi]$. The physical domain is the arrival space of the following circular mapping,

$$ \mathcal{F}: \hat{\Omega} \to \Omega, (r, \theta) \mapsto \begin{cases} x = r \cos(\theta), \\ y = r \sin(\theta). \end{cases} $$

The reference simulation in [37], is run on a grid of N<sub>r</sub> × N<sub>θ</sub> = 128 × 256 cells, with a time step of dt = 0.1. The initial condition is defined with the amplitude and mode perturbation, ε = 10⁻⁴ and m = 9. It is a ring of diameters r<sup>-</sup> = 0.45 and r<sup>+</sup> = 0.50. The parameters of the smoothing function are: $\bar{r} = \frac{r^{-}+r^{+}}{2}$, d = $\frac{r^{+}-r^{-}}{2}$ and p = 50. The coupling of the two equations is made with an explicit predictor-corrector method described in [37]. Fig. 17 shows the evolution of the density profile for different time steps during the linear regime. The perturbation introduced in the initial condition grows linearly and creates m vortices.
<figure><img src="image_17.png" /><figcaption>Figure 17: Density profile for different time steps. On the right, the initial condition ($t = 0.0$). In the middle, end of the linear regime ($= 50.0$). On the left, after the linear regime ($t = 70.0$), $m = 9$ vortexes are formed. This is also the final time chosen in the simulation in the article [37].</figcaption></figure>

In the article [37], the simulations are run on a global spline. We will compare the results of this simulation to a simulation run with local splines. To compute the derivatives at the interface, we will use the exact formulae developed in (19) and (29). The difference between a simulation on a global spline and a simulation on local splines is computed by taking the infinity norm of the difference of the density profile,

$$ \max_{t \ge 0} \|\rho_{\text{global}} - \rho_{\text{local}}\|_{L_{\infty}} = \max_{t \ge 0} \left( \max_{(r_i, \theta_j)} \left( |\rho_{\text{global}}(r_i, \theta_j) - \rho_{\text{local}}(r_i, \theta_j)| \right) \right). $$

Two tests have been implemented in Gyselalib++ and are presented below: (Test 3.1) an introduction of patches on the same simulation as in the article [37] in 5.3; (Test 3.2) a modification of the simulation mesh to add a local refinement with an interface in the middle of the solution and a cell depopulation of the central patch in a 2D multi-patch geometry with T-joints in 5.3.

**Remark.** The study in this article focuses on the advection on a multi-patch grid. The Poisson equation of the guiding-center equations (2) is solved on an equivalent global spline of the local splines used for the advection. The method used is exactly the same as described in the reference article [37]. It has been implemented in Gyselalib++ [1]. For this reason, the results presented below will be only on conforming global meshes.

**List of tests.** The diocotron instability simulation has been run on different multi-patch grids. Table 5 contains a list of the different parameters we chose for the grids.

<table><caption>Table 5: Diocotron instability simulation. Number of cells and domains of the patches and equivalent global domains for the advection tests. In parentheses, we indicate the type of splitting. If the multi-patch splitting is 1D, we only need to compute the derivatives along $r$. If the splitting is 2D, we compute the derivatives along $r$, along $θ$ and the cross-derivatives.</caption><thead><tr><th>Test</th><th>Patch</th><th>Number of cells</th><th>Domain along r</th><th>Domain along θ</th></tr></thead><tbody><tr><td rowspan="4">3.1 Uniform (1D)</td><td>Patch1</td><td>[42 × 256]</td><td>[0, 42/128]</td><td>[0, 2π]</td></tr><tr><td>Patch2</td><td>[44 × 256]</td><td>[42/128, 86/128]</td><td>[0, 2π]</td></tr><tr><td>Patch3</td><td>[42 × 256]</td><td>[86/128, 1]</td><td>[0, 2π]</td></tr><tr><td>Global</td><td>[128 × 256]</td><td>[0, 1]</td><td>[0, 2π]</td></tr><tr><td rowspan="6">3.2 T-joints with refinement (2D)</td><td>Patch1</td><td>[4 × 255]</td><td>[0, 42/128]</td><td>[0, 2π]</td></tr><tr><td>Patch2</td><td>[72 × 85]</td><td>[42/128, 60/128]</td><td>[0, 2π/3]</td></tr><tr><td>Patch3</td><td>[72 × 85]</td><td>[42/128, 60/128]</td><td>[2π/3, 4π/3]</td></tr><tr><td>Patch4</td><td>[72 × 85]</td><td>[42/128, 60/128]</td><td>[4π/3, 2π]</td></tr><tr><td>Patch5</td><td>[48 × 255]</td><td>[60/128, 1]</td><td>[0, 2π]</td></tr><tr><td>Global</td><td>[124 × 255]</td><td>[0, 1]</td><td>[0, 2π]</td></tr></tbody></table>

### Test 3.1 Uniform global mesh
This first simulation is run on the same mesh as in [37]. The mesh with 128 × 256 cells is split into three patches. The number of total cells and their shape are the same as in the article [37].

The results are presented in Fig. 18. We expect an error of $10^{-14}$ between the two curves at each interpolation (see Table 3). The time step is $\Delta t = 0.1$ and the final time, $T = 70.0$. There are 700 iterations. In the predictor-corrector method applied to solve the guiding-center equations, the density is interpolated twice and the advection field once. So, obtaining an error of $10^{-10}$ in 18 is consistent with the estimations. The results are not as good as for the simple advection. This is due to the errors on the advection field, which is also interpolated using the same method.

Fig. 19 confirms that the $L^2$-norm of the perturbation in the multi-patch case gives rise to the expected slope of 0.18 (for $m = 9$) as calculated in [37]. It does not affect the mass variation. The variation stays below 0.25% of the initial mass.

<figure><img src="image_18.png" /><figcaption>Figure 18: Diocotron instability simulation. On the left, results with advection on a global spline at $t = 70.0$. This result agrees with the result in article [37]. In the middle, results with advection on local splines at $t = 70.0$. The three yellow circles represent the interfaces between the patches. On the right, the difference ($\rho_{\text{global}} - \rho_{\text{local}}$) at the same time.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 19: Diocotron instability simulation. Verification of the results: $L^2$-norm of the perturbed density and relative errors on the total mass. As expected, the growth rate of the $L^2$-norm follows the analytical one given in Fig. 6 of [37]. The relative errors on the mass are small.</figcaption></figure>

**Test 3.2. 2D multi-patch case with T-joints and local refined** We now test the simulation on a domain with a 2D multi-patch case. The patches also form T-joints with the interfaces. The simulation uses the same layout as in Fig. 13, in paragraph 4.3.2.

One advantage of the multi-patch approach is that we can deal with local refinement while working on a uniform B-spline basis, which can be more optimized than the non-uniform ones. Here, the most interesting patch is Patch 2. The density is quasi-null in the other patches. So, we refine Patch 2 and remove some cells in the patches 1 and 3. The global domain is no longer uniform but still conforming.

We depopulated the central patch 1 to the minimum number of cells necessary to define a cubic spline (i.e. four cells).
In addition, we tested shifting one interface to the middle of the solution so as not to be influenced by the fact that the density was quasi-null at the interfaces in the previous simulations. The interface between patch 2 and patch 3 has been shifted to the middle of the ring of the initial density profile.

The patches are defined in Table 5. The results are presented in Fig. 20. The error here is higher than in the previous simulation. This is due to the less adapted mesh. The local spline on patch 1 is not precise due to the depopulation. The location of the interfaces with patch 5 is not optimized.

Fig. 21 shows the admissibility of this simulation despite the depopulation in patch 1 and the ill-placed interfaces at $r = \frac{60}{128}$. The analytically predicted slope of the $L^2$-norm of the perturbation ([37], Fig. 6) is retrieved. The mass variation reaches 7% of the initial mass. This choice of mesh is not recommended for this simulation. However, despite this mesh being poorly adapted, the local splines follow the equivalent global spline, and the results show a good agreement between them.

<figure><img src="image_20.png" /><figcaption>Figure 20: Diocotron instability simulation with a 2D multi-patch case and T-joints. We observe results similar to the previous simulations despite the use of a non-optimal mesh for this simulation (ill-placed interface and depopulation of the central patch). The local splines and the equivalent global spline still agree.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 21: Diocotron instability simulation with a 2D multi-patch case and T-joints. Verification of the results: $L^2$-norm of the perturbed density and relative errors on the total mass. We check that the suggested simulation on T-joints with refinement is an admissible simulation. The growth rate of the $L^2$-norm follows as expected the analytical one given in [37], Fig. 6. The relative errors on the mass stay small.</figcaption></figure>

# 6 Conclusions and perspectives

We presented a method to solve a classical advection equation on a multi-patch domain with the backward semi-Lagrangian method. In the backward semi-Lagrangian method, we mainly focus on interpolating the function at the previous step. The interpolation is made on cubic B-splines. In a multi-patch geometry, the conditions used to enforce patch connections are important to avoid instabilities. In this paper, we showed that a $C^0$ regularity at the patch interfaces may be not sufficient. Instead, we
impose a $C^1$ regularity. We then work with Hermite boundary conditions and determine the derivatives at interfaces every time we want to interpolate a function. Starting from a formula developed in [10], we demonstrate here a new generalization to compute the derivatives. This allows us to work with a variety of meshes (uniform and non-uniform with any cell length, conforming and some non-conforming). The generalized formula is exact but requires to solve a small global linear system. However, as the dependency on points far away from the interface decreases with distance, more local formulas can be derived at any precision. In the uniform per patch case, we also give an explicit formulation and study its asymptotic behavior when the number of cells per patch increases. The ideal case may be an approximation using thirty cells per patch to reconcile locality of the computations and machine precision. Section 4 also describes the application on 2D multi-patch domains, especially the non-conforming and T-joints cases. The numerical results of Section 5 show good agreement between the local splines defined on the patches and an equivalent global spline (if it exists).

In the future, we will consider developing new tools and methods to deal with X-point geometries as described in the Introduction 1. The X-point raises new difficulties that cannot be handled using the same methods as for the O-point. The study in this article with T-joints should allow us to connect O-point geometries to X-point geometries in a next step. So, once interpolations on X-point geometry are handled, we will be able to advect on a complete tokamak-like poloidal cross-section as represented in Fig. 1.

## References

[1] E. Bourne. *Non-Uniform Numerical Schemes for the Modelling of Turbulence in the 5D GYSELA Code*. PhD thesis, ED 184 MATHEMATIQUES ET INFORMATIQUE, December 2022.

[2] E. Bourne, V. Grandgirard, Y. Asahi, J. Bigot, P. Donnel, A. Hoffmann, A. Kara, P. Krah, B. Legouix, E. Malaboeuf, D. Midou, Y. Munschy, M. Peybernes, M. Protais, K. Obrejan, T. Padioleau, and P. Vidal. Gyselalib++.

[3] E. Bourne, P. Leleux, K. Kormann, C. Kruse, V. Grandgirard, Y. Güçlü, M. J. Kühn, U. Rüde, E. Sonnendrücker, and E. Zoni. Solver comparison for poisson-like equations on tokamak geometries. *Journal of Computational Physics*, 488:112249, 2023.

[4] E. Bourne, Y. Munschy, V. Grandgirard, M. Mehrenberger, and P. Ghendrih. Non-uniform splines for semi-lagrangian kinetic simulations of the plasma sheath. *Journal of Computational Physics*, 488:112229, 2023.

[5] M. Campos-Pinto and Y. Güçlü. Broken-feec discretizations and hodge laplace problems, 2024.

[6] J. Carrero, B. Cockburn, and D. Schötzau. Hybridized globally divergence-free ldg methods. part i: The stokes problem. *Mathematics of computation*, 75(254):533–563, 2006.

[7] B. Cockburn. *Static Condensation, Hybridization, and the Devising of the HDG Methods*, pages 129–177. Springer International Publishing, Cham, 2016.

[8] B. Cockburn, J. Gopalakrishnan, and R. Lazarov. Unified hybridization of discontinuous galerkin, mixed, and continuous galerkin methods for second order elliptic problems. *SIAM Journal on Numerical Analysis*, 47(2):1319–1365, 2009.

[9] R. Courant, K. Friedrichs, and H. Lewy. On the partial difference equations of mathematical physics. *IBM Journal of Research and Development*, 11(2):215–234, March 1967.

[10] N. Crouseilles, G. Latu, and E. Sonnendrücker. A parallel vlasov solver based on local cubic spline interpolation on patches. *Journal of Computational Physics*, 228(5):1429–1446, 2009.

[11] O. Czarny and G. Huysmans. Bézier surfaces and finite elements for mhd simulations. *Journal of Computational Physics*, 227(16):7423–7445, 2008.

[12] R. C. Davidson. *Physics of non neutral plasmas*. Clarendon Press, 1990.
[13] Eric Madaule, Sever Adrian Hirstoaga, Michel Mehrenberger and Jérôme Pétri. Semi-Lagrangian simulations of the diocotron instability. *HAL*, July 2013.

[14] G. Farin. *Curves and surfaces for computer-aided geometric design: a practical guide*. Elsevier, 2014.

[15] V. Grandgirard, J. Abiteboul, J. Bigot, T. Cartier-Michaud, N. Crouseilles, G. Dif-Pradalier, C. Ehrlacher, D. Esteve, X. Garbet, P. Ghendrih, G. Latu, M. Mehrenberger, C. Norscini, C. Passeron, F. Rozar, Y. Sarazin, E. Sonnendrücker, A. Strugarek, and D. Zarzoso. A 5d gyrokinetic full-f global semi-lagrangian code for flux-driven ion turbulence simulations. *Computer Physics Communications*, 207:35-68, 2016.

[16] V. Grandgirard, M. Brunetti, P. Bertrand, N. Besse, X. Garbet, P. Ghendrih, G. Manfredi, Y. Sarazin, O. Sauter, E. Sonnendrücker, J. Vaclavik, and L. Villard. A drift-kinetic semi-lagrangian 4d code for ion turbulence simulation. *Journal of Computational Physics*, 217(2):395-423, 2006.

[17] T. N. E. Greville. Numerical procedures for interpolation by spline functions. *Journal of the Society for Industrial and Applied Mathematics Series B Numerical Analysis*, 1(1):53-68, 1964.

[18] T. N. E. Greville. Theory and applications of spline functions. *Theory and applications of spline functions*, 1969.

[19] J. Gu, T. Yu, L. Van Lich, T.-T. Nguyen, and T. Q. Bui. Adaptive multi-patch isogeometric analysis based on locally refined b-splines. *Computer Methods in Applied Mechanics and Engineering*, 339:704-738, 2018.

[20] Y. Güçlü, S. Hadjout, and M. C. Pinto. A broken feec framework for electromagnetic problems on mapped multipatch domains, 2022.

[21] R. W. Johnson. Higher order b-spline collocation at the greville abscissae. *Applied Numerical Mathematics*, 52(1):63-75, 2005.

[22] K. Kormann, K. Reuter, and M. Rampp. A massively parallel semi-lagrangian solver for the six-dimensional vlasov-poisson equation. *The International Journal of High Performance Computing Applications*, 33(5):924-947, 2019.

[23] G. Latu, M. Mehrenberger, Y. Güçlü, M. Ottaviani, and E. Sonnendrücker. Field-aligned interpolation for semi-lagrangian gyrokinetic simulations. *Journal of Scientific Computing*, 74:1601-1650, 2018.

[24] G. Latu, M. Mehrenberger, Y. Güçlü, M. Ottaviani, and E. Sonnendrücker. Field-aligned interpolation for semi-lagrangian gyrokinetic simulations. *Journal of Scientific Computing*, 74:1601-1650, 2018.

[25] W. T. Les Piegl. *The NURBS Book*. Springer Berlin, Heidelberg, November 1996.

[26] L. Mendoza. *A new approach discretising the 2D poloidal plane of fusion devices*. PhD thesis, Technische Universität München, 2017.

[27] N. Crouseilles, M. Mehrenberger, and F. Vecil. Discontinuous galerkin semi-lagrangian method for vlasov-poisson. *ESAIM: Proc.*, 32:211-230, 2011.

[28] N. M. Patrikalakis and T. Maekawa. *Representation of Curves and Surfaces*, pages 1-33. Springer Berlin Heidelberg, Berlin, Heidelberg, 2002.

[29] J. Pétri. The diocotron instability in a pulsar cylindrical electrosphere. *Astronomy & Astrophysics*, 464(1):135-142, 2007.

[30] J. Pétri. Non-linear evolution of the diocotron instability in a pulsar electrosphere: two-dimensional particle-in-cell simulations. *Astronomy & Astrophysics*, 503(1):1-12, 2009.

[31] X. Qian. Topology optimization in b-spline space. *Computer Methods in Applied Mechanics and Engineering*, 265:15-35, 2013.
[32] N. Rivals, P. Tamain, Y. Marandet, X. Bonnin, J.-S. Park, H. Bufferand, R. Pitts, G. Falchetto, H. Yang, and G. Ciraolo. Soledge3x full vessel plasma boundary simulations of iter non-active phase plasmas. *Nuclear Fusion*, 65(2):026038, jan 2025.

[33] B. Seibold, R. R. Rosales, and J.-C. Nave. Jet schemes for advection problems. *Discrete and Continuous Dynamical Systems - B*, 17(4):1229-1259, 2012.

[34] E. Sonnendrücker, J. Roche, P. Bertrand, and A. Ghizzo. The semi-lagrangian method for the numerical resolution of the vlasov equation. *Journal of Computational Physics*, 149(2):201-220, 1999.

[35] M. Unser, A. Aldroubi, and M. Eden. Fast b-spline transforms for continuous image representation and interpolation. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 13(3):277-285, 1991.

[36] E. Zoni. *Theoretical and numerical studies of gyrokinetic models for shaped Tokamak plasmas*. PhD thesis, Technische Universität München, 2019.

[37] E. Zoni and Y. Güçlü. Solving hyperbolic-elliptic problems on singular mapped disk-like domains with the method of characteristics and spline finite elements. *Journal of Computational Physics*, 398:108889, 2019.
# Appendix A. Knot sequences to define B-splines

For $N_c$ cells, the 1D spline $s$ of degree $d = 3$ is given by

$$s(x) = \sum_{k=0}^{N_c+d-1} c_k b_k(x), \quad x \in \Omega \qquad (34)$$

with $\{b_k\}_k$ the B-splines, and $\{c_k\}_k$ coefficients to determine. There are $N_c + 1$ break points $\{x_{b,i}\}_{i=0}^{N_c}$. We add $N_c+2d$ knots $\{k_i\}_{i=0}^{N_c+2d-1}$ to build the B-splines. The knot sequence is generally given as follows (especially in the uniform case),

$$\left\{ \begin{array}{ll} k_i = x_{b,0} - \frac{x_{b,N_c} - x_{b,0}}{N_c} (d-i), & 0 \le i < d, \\ k_i = x_{b,i-d}, & d \le i < N_c + d, \\ k_i = x_{b,N_c} + \frac{x_{b,N_c} - x_{b,0}}{N_c} (i - N_c + d + 1), & N_c + d \le i < N_c + 2d. \end{array} \right.$$

The Fig. 22 shows B-splines we obtain with this type of knot sequence.

<figure><img src="image_22.png" /><figcaption>Figure 22: Cubic B-splines defined on $N_c = 4$ uniform cells with a knot sequence with uniform outside knots (•).</figcaption></figure>

In the $2\pi$-periodic case, they can be given by

$$\left\{ \begin{array}{ll} k_i = x_{b,N_c-1-i} - 2\pi, & 0 \le i < d, \\ k_i = x_{b,i-d}, & d \le i < N_c + d, \\ k_i = x_{b,i-N_c+d} + 2\pi, & N_c + d \le i < N_c + 2d. \end{array} \right.$$

we obtained B-splines basis is similar to the one shown in Fig. 22.

We call open knots sequence, if they are defined with repetitions of the first and the last break points as follows

$$\left\{ \begin{array}{ll} k_i = x_{b,0}, & 0 \le i < d, \\ k_i = x_{b,i-d}, & d \le i < N_c + d, \\ k_i = x_{b,N_c}, & N_c + d \le i < N_c + 2d. \end{array} \right.$$

<figure><img src="image_23.png" /><figcaption>Figure 23: Cubic B-splines defined on $N_c = 4$ uniform cells with an open knot sequence (•).</figcaption></figure>
## Appendix B. Proof hints of the explicit formulation of the coefficients

We detail here some steps to demonstrate the equality between the recursive formulations (15), (17) and the explicit formulations (24), (25) of the coefficients in the uniform per patch case. The index $i$ corresponds to the index of the interpolation point $x_i$ located on the interface $X_I$.

In the uniform per patch case, the coefficients are recursively given by, for $n = 2$ and $m = 1$,

$$c_{1,2}^{i} = \frac{1}{4 + a_{1,1}^{i}} \left[ 4c_{1,1}^{i} + 3 \frac{f_{i+2}^{+} - f_{i}^{+}}{\Delta x^{+}} a_{1,1}^{i} \right], \quad a_{1,2}^{i} = \frac{-a_{1,1}^{i}}{4 + a_{1,1}^{i}}, \quad b_{1,2}^{i} = \frac{4b_{1,1}^{i}}{4 + a_{1,1}^{i}},$$

for $n > 2$ and $m = 1$,

$$\left\{
\begin{aligned}
c_{1,n+1}^{i} &= \frac{1}{4 + \frac{a_{1,n}^{i}}{a_{1,n-1}^{i}}} \left[ 4c_{1,n}^{i} + \frac{a_{1,n}^{i}}{a_{1,n-1}^{i}} c_{1,n-1}^{i} + 3 \frac{f_{i+n+1}^{+} - f_{i+n-1}^{+}}{\Delta x^{+}} a_{1,n}^{i} \right] \\
a_{1,n+1}^{i} &= \frac{-a_{1,n}^{i}}{4 + \frac{a_{1,n}^{i}}{a_{1,n-1}^{i}}}, \\
b_{1,n+1}^{i} &= \frac{1}{4 + \frac{a_{1,n}^{i}}{a_{1,n-1}^{i}}} \left[ 4b_{1,n}^{i} + \frac{a_{1,n}^{i}}{a_{1,n-1}^{i}} b_{1,n-1}^{i} \right],
\end{aligned}
\right.$$

for $n = N_I^R$ and $m = 2$,

$$c_{2,N^R}^{i} = \frac{1}{4 + b_{1,N^R}^{i}} \left[ 4c_{1,N^R}^{i} + 3 \frac{f_{i}^{-} - f_{i-2}^{-}}{\Delta x^{-}} b_{1,N^R}^{i} \right], \quad a_{2,N^R}^{i} = \frac{4a_{1,N^R}^{i}}{4 + b_{1,N^R}^{i}}, \quad b_{2,N^R}^{i} = \frac{-b_{1,N^R}^{i}}{4 + b_{1,N^R}^{i}},$$

for $n = N_I^R$ and $m > 2$,

$$\left\{
\begin{aligned}
c_{m+1,N^R}^{i} &= \frac{1}{4 + \frac{b_{m,N^R}^{i}}{b_{m-1,N^R}^{i}}} \left[ 4c_{m,N^R}^{i} + \frac{b_{m,N^R}^{i}}{b_{m-1,N^R}^{i}} c_{m-1,N^R}^{i} + 3 \frac{f_{i-m+1}^{-} - f_{i-m-1}^{-}}{\Delta x^{-}} b_{m,N^R}^{i} \right], \\
a_{m+1,N^R}^{i} &= \frac{1}{4 + \frac{b_{m,N^R}^{i}}{b_{m-1,N^R}^{i}}} \left[ 4a_{m,N^R}^{i} + \frac{b_{m,N^R}^{i}}{b_{m-1,N^R}^{i}} a_{m-1,N^R}^{i} \right]. \\
b_{m+1,N^R}^{i} &= \frac{-b_{m,N^R}^{i}}{4 + \frac{b_{m,N^R}^{i}}{b_{m-1,N^R}^{i}}}.
\end{aligned}
\right.$$

Looking at the expressions of the {$a_{1,n}^i$}$_{n=2,...,N_I^R}$ and {$b_{m,N^R}^i$}$_{m=2,...,N_I^L}$ sequences, we start by introducing the $w$ sequence such that $\forall k > 1$,

$$w_k = \frac{-1}{4 + w_{k-1}}. \qquad (35)$$

The function $x \mapsto \frac{-1}{4+x}$ has two fixed points: $-(2 - \sqrt{3})$ and $-(2 + \sqrt{3})$. Successively injecting the relation (35) into itself gives us

$$w_k = \frac{-(4 + w_{k-2})}{4(4 - 1) + 4w_{k-2}}, \quad w_k = \frac{-(4(4 - 1) + 4w_{k-2})}{4(4(4 - 1) - 1) + 4(4 - 1)w_{k-2}}.$$

We can see appearing the sequence $u$ such that $u_{k+1} = 4u_k - u_{k-1}$, $\forall k \in \mathbb{Z}$ with $u_0 = 0$ and $u_1 = 1$. It is a linear-recursive sequence whose characteristic polynomial $x^2 + 4x - 1$ has two zeros: $(2 - \sqrt{3})$
and $(2 + \sqrt{3})$. The sequence $u$ is then in the vector space span($\{\{(2 - \sqrt{3})^k\}_{k \in \mathbb{N}}, \{(2 + \sqrt{3})^k\}_{k \in \mathbb{N}}\}$) and can be expressed as follows $\forall k > 1$,

$$u_k = (2 + \sqrt{3})^k - (2 - \sqrt{3})^k,$$

and we can prove that $w$ can be expressed from $u$ as following

$$w_k = -\frac{u_{k-1} + u_{k-2}w_1}{u_k + u_{k-1}w_1}.$$

Indeed, $\{-\frac{u_{k-1}+u_{k-2}w_1}{u_k+u_{k-1}w_1}\}_{k \ge 1}$ verifies the relation (35). For $k = 1$, $\{-\frac{u_{k-1}+u_{k-2}w_1}{u_k+u_{k-1}w_1}\}_{k \ge 1}$ is equal to the same initial condition $w_1$ as $w$.

It is easy to see that by fixing $w_1 = a_{1,1}^i$, the $\{w_n\}_{n=2,...,N_I^R}$ sequence gives the same values as the $\{\frac{a_{1,n}^i}{a_{1,n-1}^i}\}_{n=2,...,N_I^R}$ sequence. By fixing $w_1 = b_{1,N_I^R}^i$, the $\{w_m\}_{m=2,...,N_I^L}$ sequence gives the same values as the $\{\frac{b_{m,N_I^R}^i}{b_{m-1,N_I^R}^i}\}_{m=2,...,N_I^L}$ sequence.

By factoring the two equalities, we deduce explicit formulas for $\{a_{1,n}^i\}_{n=2,...,N_I^R}$ and $\{b_{m,N_I^R}^i\}_{m=2,...,N_I^L}$,

$$a_{1,n}^i = \frac{(-1)^{n-1}u_1 a_{1,1}^i}{u_n + u_{n-1}a_{1,1}^i}, \quad b_{m,N_I^R}^i = \frac{(-1)^{m-1}u_1 b_{1,N_I^R}^i}{u_m + u_{m-1}b_{1,N_I^R}^i}.$$

Manipulating the first terms of the $\{b_{1,n}^i\}_{n=1,...,N_I^R}$ and $\{a_{m,N_I^R}^i\}_{m=1,...,N_I^L}$ sequences helps us to identify the following expressions,

$$b_{1,n}^i = \frac{u_n b_{1,1}^i}{u_n + u_{n-1}a_{1,n}^i}, \quad a_{m,N_I^R}^i = \frac{u_m a_{1,N_I^R}^i}{u_m + u_{m-1}b_{1,N_I^R}^i},$$

which can be shown by recursion.

We then inject the expressions for $a_{1,N_I^R}^i$ and $b_{1,N_I^R}^i$ into $a_{m,n}^i$ and $b_{m,n}^i$ for $n = N_I^R$,

$$a_{m,n}^i = \frac{(-1)^{n-1}u_1 a_{1,1}^i u_m}{u_m u_n + u_m u_{n-1} a_{1,1}^i + u_n u_{m-1} b_{1,1}^i}, \quad b_{m,n}^i = \frac{(-1)^{m-1}u_1 b_{1,1}^i u_n}{u_m u_n + u_m u_{n-1} a_{1,1}^i + u_n u_{m-1} b_{1,1}^i}.$$

The expressions for the weights $\{\omega_{k,N_I^L,N_I^R}\}_{k=-N_I^L,...,N_I^R}$ in the coefficient $c_{1,n-1}^i$ are also demonstrated by recursion. We start by the forward recursion which can be written as follows thanks to the results on $a_{m,n}^i$ and $b_{m,n}^i$ for $n \in \mathbb{N}^*, m = 1$,

$$c_{1,n+1}^i = 4 \frac{u_n + u_{n-1}a_{1,1}}{u_{n+1} + u_n a_{1,1}} c_{1,n}^i - \frac{u_{n-1} + u_{n-2}a_{1,1}}{u_{n+1} + u_n a_{1,1}} c_{1,n-1}^i - \frac{3(-1)^{n-1}u_1 a_{1,1}}{u_{n+1} + u_n a_{1,1}} \frac{f_{i-1+n}^+}{\Delta x^+} + \frac{3(-1)^{n-1}u_1 a_{1,1}}{u_{n+1} + u_n a_{1,1}} \frac{f_{i+1+n}^+}{\Delta x^+}.$$

We identify six cases where $k = n + 1$, $k = n$, $k = n - 1$, $k = 1, ..., n - 2$, $k = 0$ and $k = -1$. Manipulating the first terms helps us to identify the following expressions,

$$\left\{
\begin{array}{ll}
\omega_{k,1,n} = \frac{3(-1)^k}{u_n + u_{n-1}a_{1,1}} \frac{a_{1,1}}{\Delta x^+} u_1, & k = n \\
\\
\omega_{k,1,n} = \frac{3(-1)^k}{u_n + u_{n-1}a_{1,1}} \frac{a_{1,1}}{\Delta x^+} (u_{n-k+1} - u_{n-k-1}), & k = 1, ..., n-1 \\
\\
\omega_{k,1,n} = \frac{3(-1)^k}{u_n + u_{n-1}a_{1,1}} \left( \frac{a_{1,1}}{\Delta x^+} (u_n - u_{n-1}) - \frac{b_{1,1}}{\Delta x^-} u_n \right), & k = 0 \\
\\
\omega_{k,1,n} = \frac{3(-1)^{k+1}}{u_n + u_{n-1}a_{1,1}} \frac{b_{1,1}}{\Delta x^-} u_n, & k = -1
\end{array}
\right.$$

which can be shown by recursion. It is not detailed here, there are no particular difficulties.
We inject these results in the backward recursion. Rewritten, it gives

$$
\left\{
\begin{aligned}
c_{2,n}^i &= \frac{4(u_n + u_{n-1}a_{1,1}^i)}{u_n(4 + b_{1,1}^i) + u_{n-1}4a_{1,1}^i} c_{n,1}^i + \frac{3u_n b_{1,1}^i}{u_n(4 + b_{1,1}^i) + u_{n-1}4a_{1,1}^i} \frac{f_i^+ - f_{i-2}^-}{\Delta x^-} \\
c_{m+1,n}^i &= \frac{1}{u_n u_{m+1} + u_n u_m b_{1,1}^i + u_{n-1} u_{m+1} a_{1,1}^i} \left[
\begin{aligned}
&4(u_n u_m + u_n u_{m-1} b_{1,1}^i + u_{n-1} u_m a_{1,1}^i) c_{m,n}^i \\
&+ 3(-1)^{m-1} u_1 b_{1,1}^i u_n \frac{f_{i+1-m}^- - f_{i-1-m}^-}{\Delta x^-} \\
&- (u_n u_{m-1} + u_n u_{m-2} b_{1,1}^i + u_{n-1} u_{m-1} a_{1,1}^i) c_{m-1,n}^i
\end{aligned}
\right].
\end{aligned}
\right.
$$

We also identify different cases where $k = n, k = 1, \dots, n-1, k = 0, k = -(m-2), \dots, -1, k = -(m-1)$, $k = -m$ and $k = -(m+1)$. It leads us to the followings expressions,

$$
\left\{
\begin{aligned}
\omega_{k,m,n}^i &= 3(-1)^k \frac{\frac{a_{1,1}^i}{\Delta x^+} u_m u_1}{u_m u_n + u_m u_{n-1} a_{1,1}^i + u_n u_{m-1} b_{1,1}^i}, && k = n, \\
\omega_{k,m,n}^i &= 3(-1)^k \frac{\frac{a_{1,1}^i}{\Delta x^+} u_m (u_{n-k+1} - u_{n-k-1})}{u_m u_n + u_m u_{n-1} a_{1,1}^i + u_n u_{m-1} b_{1,1}^i}, && k = 1, \dots, n-1, \\
\omega_{k,m,n}^i &= 3(-1)^k \frac{\frac{a_{1,1}^i}{\Delta x^+} u_m (u_n - u_{n-1}) - \frac{b_{1,1}^i}{\Delta x^-} u_n (u_m - u_{m-1})}{u_m u_n + u_m u_{n-1} a_{1,1} + u_n u_{m-1} b_{1,1}^i}, && k = 0, \\
\omega_{k,m,n}^i &= 3(-1)^{k+1} \frac{\frac{b_{1,1}^i}{\Delta x^-} u_n (u_{m+k+1} - u_{m+k-1})}{u_m u_n + u_m u_{n-1} a_{1,1}^i + u_n u_{m-1} b_{1,1}^i}, && k = -(m-1), \dots, -1, \\
\omega_{k,m,n}^i &= 3(-1)^{k+1} \frac{\frac{b_{1,1}^i}{\Delta x^-} u_n u_1}{u_m u_n + u_m u_{n-1} a_{1,1}^i + u_n u_{m-1} b_{1,1}^i}, && k = -m.
\end{aligned}
\right.
$$

The demonstration is also made by recursion.