A Characteristic Mapping Method with Source Terms:
Applications to Ideal Magnetohydrodynamics

Consider the scalar advection equation (3) with solution given in terms of the initial condition and the characteristic map (4). We note that the evolution of the scalar $\theta(t)$ is given by a linear action of the map in the sense that $\theta(t)={\bm{X}}_{[t,0]}^{*}\theta_{0}=\theta_{0}\circ{\bm{X}}_{[t,0]}$ and, for $a,b$ arbitrary, $a\theta(t)+b\phi(t)={\bm{X}}_{[t,0]}^{*}(a\theta_{0}+b\phi_{0})$ for $\theta,\phi$ solutions to linear homogeneous advection equations with initial conditions $\theta_{0},\phi_{0}$, under the same velocity field ${\bm{u}}$. Solutions to the inhomogeneous equation can therefore be obtained from the homogeneous equations through Duhamel’s principle. This in fact also holds for the Lie-advection of differential $k$-forms where ${\bm{X}}_{[t,0]}^{*}$ denotes the pullback operator, and more generally, for linear (right) actions of the diffeomorphism group.

Therefore, for the inhomogeneous advection equation

we can perform the decomposition

by commuting pullback with time integration. Defining

we can write a recursive formula for the submap decomposition expression

for $0=\tau_{0}<\tau_{1}<\tau_{2}<\ldots<\tau_{n-1}<\tau_{n}=t$ and $\theta(\cdot,0)=\theta_{0}(\cdot)$. Here, we note that within one remapping subinterval, for $t\in[\tau_{i-1},\tau_{i}]$, the source term accumulation $F_{[t,\tau_{i-1}]}$ as a function of $t$ satisfies

Therefore, in the CMM for inhomogeneous advection we must at the same time discretize two quantities, the current submap ${\bm{X}}_{[t,\tau_{i-1}]}$ and the current source term accumulation $F_{[t,\tau_{i-1}]}$.

In the case of equations with source terms, during a given time subinterval $[\tau_{i},\tau_{i+1}]$, the accumulated source term $F_{[t,\tau_{i}]}$ must be computed by discretizing equation (12). Alternatively, given the flow maps ${\bm{X}}_{[t,s]}$ of the velocity field ${\bm{u}}$, we may write $F_{[t_{n+1},\tau_{i}]}$ explicitly using Duhamel’s principle as a function of the map, $f$ and $F_{[t_{n},\tau_{i}]}$:

This gives us a direct semi-Lagrangian discretization formula:

where $a_{j},s_{j}$ are $k$-stage the Runge-Kutta time quadrature weights and loci. In our current implementation, we opted to use the same Runge-Kutta scheme for the one-step map and source term integrations. Therefore, the one-step map $\tilde{\bm{X}}_{[t_{n+1},t_{n}]}$ and its intermediate stage values $\tilde{\bm{X}}_{[t_{n+1},s_{j}]}$ were already computed during the map integration and therefore may be reused at no additional cost.

The source term $\tilde{f}$ is also defined in $[t_{n},t_{n+1}]$ using the time extrapolation operator $E_{\gamma}$ from discrete time data $\tilde{f}_{n}$. At each time $t_{n}$, the source term, in this case, the curl of the Lorentz force, is approximated by sampling the magnetic field using expression (19d) on the grid $A$. Then, the magnetic field is filtered in Fourier space to retain spatial scales well-resolved by the grid, from which we then compute the curl of the Lorentz force, giving us the source term $\tilde{f}_{n}$. We note that this operator is solely a function of ${\bm{B}}_{0}$ and the map $\mathcal{\bm{X}}_{[t,0]}$, hence we denote $\tilde{f}=S_{l}(\mathcal{\bm{X}}_{[t,0]})$ where $S_{l}$ is given by $\nabla\times\left(\nabla\times\tilde{{\bm{B}}}\times\tilde{{\bm{B}}}\right)$ where $\tilde{{\bm{B}}}$ is a filtered version of $\text{adj}(\nabla\mathcal{\bm{X}}_{[t_{n},0]}){\bm{B}}_{0}(\mathcal{\bm{X}}_{[t_{n},0]})$.

Here the method is described using two separate grids $M$ and $A$ for the map and the source term accumulation respectively. In addition to accuracy considerations, these grids’ resolution also controls the remapping frequency. Indeed, since the evolution of both $\mathcal{\bm{X}}$ and $\mathcal{F}$ are discretizations of an advection equation on a fixed Eulerian grid, convergence is lost when the numerical solution reaches the Nyquist frequency of the grid. Therefore a rule of thumb is to trigger a remapping routine whenever either $\mathcal{\bm{X}}$ and $\mathcal{F}$ develops significant small scales ill-resolved by their respective grids $M$ and $A$. The freedom of using different grids then allows for further optimization between the remapping frequencies of the two terms. Generally, using the same grid should be sufficient since both terms are advected by the same velocity and minor tweaking is possible depending on the respective right-hand-sides for the two equations.

The CMM for the ideal MHD equations in 2d is summarized in the pseudocode algorithm 1. In this algorithm, we have introduced two distinct filtering scales $l_{1}$ and $l_{2}$ for the operators $K$ and $S$ respectively. The filtering $l_{1}$ is chosen to guarantee that the resulting velocity field $\tilde{{\bm{u}}}$ is well resolved on the map grid. Indeed, any subgrid scales of the velocity in the map evolution will only appear as noise since it cannot be captured by the interpolant $\mathcal{I}_{M}$. This filtering allows us to access, for each fixed $l_{1}$, the convergence regime in the limit of fine spatial and temporal grids. Similarly, $l_{2}$ is chosen so that the source term is well resolved by $\mathcal{I}_{A}$. The convergence of the numerical solution to the exact solution then relies on the $l_{1},l_{2}$ limit which one should take after the spatial and temporal grid limits. A diagonalization argument may then be used to select a single subsequence. We give an overview of the error estimate arguments in the following section.

In this section we present a sketch of the error estimates for the CMM with source term, detailed error analysis for the CMM can be found in [31]. One notable difference is that here, instead of the momentum equation, we use the vorticity equation which, in the 2d incompressible case, is simply a scalar advection equation, thereby avoiding derivative loss due to the differential pullback operator. Under suitable assumptions on the spatial and temporal discretization operators, the following error estimate can be reached.

To simplify the analysis, we assume that the spatial discretization $\mathcal{I}_{M}$ for the map and $\mathcal{I}_{A}$ for the source term accumulation, on grids of cell width $N^{-1}$, are smooth projections on Sobolev spaces $H^{s}$ for some $s>n/2+1$ with the contraction and approximation properties

for any $0\leq r\leq s$, for instance, Fourier projections on $N^{2}$ harmonics [3] or mollification to scale $N^{-1}$. We also assume that once the numerical velocity $\tilde{\bm{u}}$ is defined in $[t_{n},t_{n+1}]$ as a Lagrange polynomial, that the numerical integration for the one-step map $\tilde{\bm{X}}_{[t_{n+1},t_{n}]}$ is performed exactly. In practice, the numerical solution obtained using Hermite cubic interpolants and RK4 integrators can then be viewed as a perturbation of the smooth scheme, with perturbation errors given by the accuracy of the interpolation and integration schemes.

We define $\tilde{\mathcal{\bm{X}}}_{[t,0]}=\mathcal{\bm{X}}_{[t_{n},0]}\circ\tilde{\bm{X}}_{[t,t_{n}]}$, thus satisfying

We may then define a time-interpolated map $\mathcal{\bm{X}}_{[t,0]}=\mathcal{I}_{M}[\tilde{\mathcal{\bm{X}}}_{[t,0]}]$ for all $t$, and this map corresponds to the numerical map obtained from the method at discrete times $t=t_{n}$. It follows from the linearity of the interpolation operator that there exists some velocity field ${\bm{v}}_{\mathcal{\bm{X}}}\in L^{1}(0,t,H^{s-1})$ which generates the map $\mathcal{\bm{X}}_{[t,0]}$. Since we may write $\mathcal{\bm{X}}_{[t,0]}$ exactly as the backward flow map of some velocity field, it follows that advected fields obtained by pullback through this map also satisfy an advection equation. Letting $\tilde{\theta}(t)=\theta_{0}\circ\mathcal{\bm{X}}_{[t,0]}$, then $\tilde{\theta}$ solves

Let $\theta(t)$ be the exact solution of the advection of the same initial field $\theta_{0}$ by the velocity field ${\bm{u}}$. Applying Grönwall’s lemma we have that for $r\leq s-1$,

In terms of the velocity error, we have that

where $N_{M}^{-1}$ is the cell width for the map grid $M$. The temporal error $\|{\bm{u}}-\tilde{\bm{u}}\|_{r}$ contains the nonlinearity as the numerical velocities $\tilde{\bm{u}}$ depend on the numerical map. The numerical velocity $\tilde{\bm{u}}$ is defined through an order $\gamma$ extension operator $E_{\gamma}$ which extrapolates a velocity $\tilde{\bm{u}}\in L^{1}(t_{n},t_{n+1},H^{s})$ from given numerical velocity data at steps $n-\gamma+1$ to $n$. We denote by $\tilde{\bm{u}}_{i}$ the numerical velocity at $t_{i}$ and ${\bm{u}}(t_{i})$ the exact velocity. The accuracy and boundedness properties

are sufficient to control the time approximation error:

It remains to control the discrete step velocity error ${\bm{u}}(t_{i})-\tilde{\bm{u}}_{t_{i}}$ which involves the map pullback, Biot-Savart kernel, and the source term integration.

The numerical accumulated source term $\mathcal{F}_{[t_{n},0]}$ is obtained from semi-Lagrangian methods using velocity $\tilde{\bm{u}}$ and source term $\tilde{f}$. We again try to write an advection equation for $\mathcal{F}_{[t,0]}=\mathcal{I}_{A}[\tilde{\mathcal{F}}_{[t,0]}]$:

where $[\mathcal{I}_{A},\tilde{\bm{u}}\cdot\nabla]$ denotes the commutator between interpolation and the material derivative. This commutator arises from the fact that the source term is computed on a fixed grid rather than directly through the map and therefore subgrid scales generated by the advection are lost. The $L^{1}(0,t,H^{q})$ error of this term is $\mathcal{O}(h^{s-q-1})$, i.e. the order of the scheme, but deteriorates accuracy once $\mathcal{F}_{[t,0]}$ contains significant small-scale features, at which point remapping is needed in order to reset this commutator error to 0.

Applying Grönwall yields the following estimate for the source term accumulation:

Again, defining $\tilde{f}(t)$ through the extension operator $E_{\gamma}$ gives an error

The control one can expect for the numerical source term depends heavily on the equation studied. For the ideal MHD equations and more generally for so-called equations with advected parameters, it turns out that $f$ is solely a function of the backward map ${\bm{X}}_{[t,0]}$ (if we consider initial data as parameters). Indeed, for the curl of the Lorentz force

we may abstractly write a source term operator $S$ as $f(t)=S({\bm{X}}_{[t,0]})$. This can then be approximated by a regularized numerical operator $S_{l_{2}}:\textit{SDiff}^{s}\to H^{q}$ which should be $C^{1}$ so that $\|S_{l_{2}}({\bm{X}}_{[t,0]})-S_{l_{2}}(\mathcal{\bm{X}}_{[t,0]})\|_{q}\leq C\|{\bm{v}}_{\mathcal{\bm{X}}}-{\bm{u}}\|_{L^{1}(0,t,H^{q})}$. The source term error can then be controlled by

where the first term is the error of the approximate source term operator along the exact solution. Therefore,

We now have everything necessary for estimating $\|{\bm{u}}(t_{n})-\tilde{\bm{u}}_{n}\|_{r}$. Denote by $K$ the Biot-Savart kernel $-\nabla\times\Delta^{-1}:H^{r-1}\to H^{r}$ and let $K_{l}$ be a regularized numerical approximation $K_{l}:H^{r-k}\to H^{r}$ for $k\geq 1$. We have the exact and numerical velocities are given by

We obtain the following bound on the error:

where we omitted any constant factors depending on the exact solution and the bounds on the $H^{s}$ norms of the numerical maps and the $H^{s-1}$ norms on the numerical source term accumulation $\mathcal{F}$. Control on the growth of $\|\mathcal{\bm{X}}_{[t,0]}\|_{s}$ is studied in [31], the growth rate of $\|\mathcal{F}_{[t,0]}\|_{s-1}$ relies on the choice of numerical solver for the advection equation. Both of these can be limited using the remapping method. In fact, since all constants in the error estimates are controlled by the $s$ seminorm of the map and $s-1$ seminorm of the accumulated source term, this suggests that the decay rate of the Fourier spectrum of the numerical map and source term data are good criteria for triggering the remapping step.

A last application of Grönwall’s inequality then gives us an order $s-r-1$ in space and $\gamma$ in time convergence for the numerical velocities and by extension the map, vorticity and magnetic fields. We note that in the $(N^{-1},{\Delta t})\to 0$ asymptotic, the leading error is the operator errors $K-K_{l_{1}}$ and $S-S_{l_{2}}$ which depend on the filtering parameters ${l_{1}},{l_{2}}$. This makes the error estimate highly non-trivial since the filtering parameters regularize the solution thereby limiting the presence of small scales in the numerical solution. On the one hand, choosing large ${l_{1}},{l_{2}}$, i.e. a stronger regularization allows the error to enter the asymptotic regime for coarser grids, however, the deviation from the exact solution is also increased by the regularization errors $K-K_{l_{1}}$ and $S-S_{l_{2}}$. On the other hand, choosing smaller filtering scales reduces these errors but also amplifies the discretization errors due to potentially faster growth in the $H^{s}$ norms of the maps. The optimal relation between ${l_{1}},{l_{2}}$ and $N^{-1},{\Delta t}$ is highly dependent on the exact solution. As a general rule of thumb, for a given filter, $N^{-1}$ and ${\Delta t}$ should be chosen sufficiently small so that all scales expected in an ${l_{1}},{l_{2}}$-regularized solution are well resolved. The exact solution can then be approached by taking the ${l_{1}},{l_{2}}\to\infty$ limit with optimal $N^{-1},{\Delta t}$ for each given filter. We note that the convergence of the numerical solution in the ${l_{1}},{l_{2}}\to\infty$ limit might not be in a strong sense with respect to an $H^{r}$ norm. This remains a difficult problem and depends on whether the exact solution can be obtained as some weak limit of strong solutions to a family of regularized problems. The above shows that the regularization employed in this method is close to that of the Langrangian-averaged MHD equations or the related MHD-$\alpha$ equations for which global well-posedness is known [19]. In this sense, the numerical results from the CMM, though fully inviscid, offer a relevant comparison with the Large Eddy Simulation models for which there is an extensive study in terms of spectral scaling and energy dissipation [25, 26, 9].

In this paper, we employ a Hermite cubic spatial interpolation with third-order Lagrange extrapolation in time. This corresponds to a perturbation of a smooth scheme which is stable as long as the numerical solution $\mathcal{\bm{X}}_{[t,0]}$ is not highly oscillatory. More precisely, the error of this perturbation is controlled by the difference of $\mathcal{\bm{X}}_{[t,0]}$ with its Fourier projection $\mathcal{I}_{M,f}[\mathcal{\bm{X}}_{[t,0]}]$ on the same grid $M$. This suggests that, in the case of Hermite cubic interpolation, we should take the $H^{4}$ norm of the Fourier transform of the map data as a remapping criterion. This holds similarly for $\mathcal{F}_{[t,0]}$.

In the first numerical example, we study linear advection with the source term

using a divergence-free velocity field that generates a swirl:

on a periodic $[0,2\pi]^{2}$ domain and for $t\in[0,2]$. To test our framework, we impose a manufactured solution:

with $\epsilon=0.1$ and where the corresponding source term $f$ is determined analytically by Eq. 45. This test case is set up such that the solution has a constant spatial profile with narrowing crossing lines towards $t=0.5$ and increasing thickness after $t>0.5$. Therefore, the source term $f$ needs to compensate for the deformation created by the swirl velocity field. The resulting evolution of the source term and solution is shown in Fig. 1. Additionally, we show the subintegrals $F_{[t,\tau_{i}]}$ for the different time instances and their Fourier spectra.

At time $t=0.4$ we see that the sub-integral and source are identical since no remapping was initiated. Furthermore, the Fourier spectra of the source term/sub-integral and analytical solution $\theta$ are similar because of the small-scale structures present in the solution at time $t=0.4$. In a source term free case, the scale in the solution would not affect the numerical precision because the evolved quantity is the map, not the state itself. Since in the presence of source terms, we have to store quantities that depend on the solution itself, the scales of the solution play a role in the numerical precision. Hence, to preserve these fine scales, storing the sub-integrals on finer grids might be necessary. As soon as the first remapping is performed, we see that the spectra of the source term and sub-integrals differ. As the sub-integrals capture only local information of the solution, they have in general a faster decay than $f(x,t)$ itself. However, note that the spectrum of the submap decays for all time instances the fastest.

Finally, we show the achieved third-order convergence in space and time in Fig. 2.

Convergence in space and time is computed for $\Delta_{x/t}\theta=\|\theta^{\text{ref}}(\cdot,T_{\text{end}})-\theta^{N}\|_{\infty}$ using the manufactured solution Eq. 47. Space convergence is achieved by increasing the number of grid points of the map $N=64,128,\dots,512$ in each dimension while keeping a constant time step of $\Delta t=1/512$ and integrating to $T_{\text{end}}=1$. Convergence in time is achieved by decreasing the step size $\Delta t=1/64,1/128,\dots,1/512$, while keeping the number of grid points of the map at $N=512$ in each dimension. For both convergence tests, we keep the sampling grid of the velocity and the grid for evaluating all quantities at $512^{2}$. For all convergence tests, remapping was switched off. We remark that as mentioned before the fine-scale structure of the solution can not be decoupled from the numerical evolved quantities, as would be the case for the source-free swirl test case. Hence, we observe that for $\epsilon\ll 0.1$ the convergence rate is reduced, especially for the coarser resolutions.

Our numerical experiments are based on the classical Orszag–Tang MHD test case that was implemented to study singularities in MHD turbulence [24].

In this case, we expect $3^{\text{rd}}$ order convergence in the $L^{2}$ norm in space and time. In fact, our numerical experiments show $3^{\text{rd}}$ convergence in the $L^{\infty}$ norm. We define the following errors: