The Characteristic Mapping Method for the Linear Advection of Arbitrary Sets

In this section, we present the mathematical formulation for the Characteristic Mapping (CM) method. First we consider the advection equation, written in its strong form:

where we assume $U\subset\mathbb{R}^{d}$ is a rectangular domain in $d$ dimensions and $\phi:U\times[0,T]\to\mathbb{R}$ is the advected quantity. Throughout the paper, we assume the given velocity field $\vec{v}(\vec{x},t)$

to be smooth and divergence-free.

We apply the method of characteristics to (1) . Let $\vec{x}$ be the solution to the following initial value problem (IVP):

where $\vec{x}$ are the characteristic curves of the equation generated by the velocity $\vec{v}$.

The method of characteristics consists in realizing that (1) simplifies to the following equations (4) and (5). We call (5) the generalized advection equation for some given velocity field $\vec{v}$ generating the characteristics $\vec{x}(t)$. This is more general in that (1) is ill-posed when regularity assumptions on $\phi$ are dropped; we can therefore use (5) to transport arbitrary functions:

or equivalently

for every time $t$ and every initial condition $\vec{x}_{0}$.

The CM method consists in computing a backward-in-time solution operator to the IVP (3). That is, we compute a map $\vec{X}:U\times[0,T]\to U$ such that for any curve $\vec{x}$ satisfying the IVP (3) with some arbitrary initial condition $\vec{x}_{0}\in U$, we have that

Existence and uniqueness for $\vec{X}$ under mild regularity assumptions have been established in [10].

One can check that $\vec{X}$ is the solution to the vector valued advection equation

We call $\vec{X}$ the characteristic map associated to the velocity $\vec{v}$. Under suitable regularity assumptions on $\vec{v}$, $\vec{X}$ is known to be a diffeomorphism.

We see that the map $\vec{X}(\cdot,t)$ takes a point $\vec{x}(t)$ at time $t$ on its characteristic curve and returns the initial position $\vec{x}_{0}$ of this curve. This allows us to obtain the solution of the advection equation (1) as

This framework provides solutions to a generalization of equation (1). Any $\phi$ satisfying equation (5) can be expressed using (8). This allows for the transport of more general sets, since $\phi$ is allowed to be any function defining a transported set. The lack of regularity assumptions on $\phi$ implies that arbitrary sets, even fractal ones, can be advected under this framework. For instance, taking $\phi_{0}$ to be the indicator function of an initial set $S_{0}$, the pullback $\phi(\cdot,t)=\phi_{0}\circ\vec{X}(\cdot,t)$ is the indicator function of

where $\vec{X}^{-1}$ denotes preimage and is valid even without the diffeomorphism assumption on $\vec{X}$.

The characteristic map $\vec{X}$ possesses a semi-group property which we exploit to significantly improve the numerical efficiency of the method. For this purpose, one can solve equation (7) in an arbitrary subinterval $[\tau_{i-1},\tau_{i}]\subseteq[0,T]$. The characteristic map for this subinterval satisfies

We denote by $\vec{X}_{[\tau_{i-1},\tau_{i}]}:U\to U$ the solution of the above map at time $t=\tau_{i}$. That is, $\vec{X}_{[\tau_{i-1},\tau_{i}]}$ represents the characteristic map generated by the velocity in the time interval $[\tau_{i-1},\tau_{i}]$. The semi-group structure allows us to decompose the characteristic map for $[0,T]$ into arbitrarily many maps of subintervals. This presents a numerical advantage since it gives us an extra degree of freedom for adaptivity in time. The decomposition is given by:

We emphasize that the computation of the characteristics map in each individual subinterval $[\tau_{i},\tau_{i+1}]$ start with the identity map as initial condition, and depends only on the velocity in this time interval. Hence, the computation is independent of any previous maps and can be solved individually.

For each subinterval, we compute the map by solving (7) numerically. For this, we employ the Gradient Augmented level-set (GALS) method [23], an unconditionally stable semi-Lagrangian method that couples high order Hermite interpolation with Runge-Kutta time stepping schemes of matching order. Since the deformation map in CM is computed using GALS, and the advected quantity is obtained by pullback whenever a solution is required, it is a direct corollary that the CM method for linear advection is also unconditionally stable.

In this paper, we use Hermite cubic polynomials for spatial interpolation. We denote by $\mathcal{H}_{\bm{C}}$ the projection operator from $C^{1}(U)$ functions to Hermite cubics on some grid denoted by ${\bm{C}}$ with cell widths ${\delta_{x}}$. To construct the Hermite cubic $\mathcal{H}_{\bm{C}}[f]$, we evaluate $f$, $\nabla f$ and mixed partial derivatives of order at most one in each spacial dimension at grid points of ${\bm{C}}$; $\mathcal{H}_{\bm{C}}[f]$ is then defined in each cell to be the Hermite cubic interpolant with these grid data. For time discretization we use a fixed time step of ${\delta_{t}}$. Time integration is performed using Runge-Kutta 3.

We denote by $\vec{\mathcal{X}}_{i}(\vec{x},t)$ the numerical approximation of $\vec{X}_{i}$ represented as a piecewise Hermite cubic. The Hermite cubic interpolant is updated using GALS stepping: we compute a one-step map $\vec{\Psi}$ which maps a point $\vec{x}$ at time $t$ to its position at time $t-{\delta_{t}}$ following its characteristic curve in $[t-{\delta_{t}},t]$, then we compose $\vec{\Psi}$ with the characteristic map at time $t-{\delta_{t}}$ and project the composition on the space of Hermite cubics. We denote by $\vec{\gamma}_{\vec{x}}(s)$ the characteristic curve satisfying (3) with $\vec{\gamma}_{\vec{x}}(t)=\vec{x}$ and compute $\vec{\Psi}$ by integrating the velocity along this curve:

The time integrations for $\vec{\gamma}_{\vec{x}}$ and $\vec{\Psi}$ are approximated using RK3. The one-step map $\vec{\Psi}$ is evaluated when updating $\vec{\mathcal{X}}_{i}$ in (12b); the spacial derivatives at grid points required for the definition of the Hermite interpolant are obtained by direct chain rule: derivatives of $\vec{\mathcal{X}}_{i}(\vec{x},t-{\delta_{t}})$ are computed from its Hermite interpolant definition and those of $\vec{\Psi}$ are directly evaluated from (12a).

Our strategy is to use the above time-stepping to evolve submaps $\vec{\mathcal{X}}_{i}$ over some time interval $[\tau_{i-1},\tau_{i}]$ which is either selected a priori or determined dynamically. Once time $\tau_{i}$ is reached, we stop the submap evolution for $\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]}$ and inductively store the composition $\vec{\mathcal{X}}_{[0,\tau_{i-1}]}\circ\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]}$ as a single Hermite cubic $\vec{\mathcal{X}}_{[0,\tau_{i}]}$ of high enough resolution; this is called the remapping step. We then repeat with the next submap $\vec{\mathcal{X}}_{[\tau_{i},\tau_{i+1}]}$.

To improve efficiency, we compute the submaps $\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]}$ on a coarse grid, since submap updates happen often, and a low spatial resolution is sufficient to represent the small deformations. On the other hand, the spatial resolution required to represent the global map is high since large deformations can accumulate over time. We therefore store the global map $\mathcal{X}_{[0,\tau_{i}]}$ on a fine grid. Updating the global map is costly, however these updates only happen at the remapping steps. Hence we must select the submap time intervals appropriately as to not incur too many remapping steps and remap only when the coarse grid is no longer able to accurately resolve the submap deformations. The interplay between the submap time intervals and the coarse and fine grid resolutions will be analyzed in the following sections.

To summmarize:

1. Given $0=\tau_{0}<\tau_{1}<\ldots<\tau_{m-1}<\tau_{m}=T$.

2. 1.

   Initialize $\vec{\mathcal{X}}_{[0,\tau_{0}]}(\vec{x})=\vec{x}$. For $i=1$ to $m$, repeat the following 2 steps:

3. 2.

   Compute $\mathcal{X}_{[\tau_{i-1},\tau_{i}]}$ on a coarse grid ${\bm{C}}$ using GALS.

4. 3.

   Compose and project on a fine grid ${\bm{F}}$, as Hermite cubic:

   (13)

   $$\displaystyle\vec{\mathcal{X}}_{[0,\tau_{i}]}=\mathcal{H}_{{\bm{F}}}\left[\vec{\mathcal{X}}_{[0,\tau_{i-1}]}\circ\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]}\right]$$

The computation of $\mathcal{X}_{[\tau_{i-1},\tau_{i}]}$ uses the GALS method. In each GALS step, the method produces two errors, one due to the approximate time integration and the other due to interpolation :

Note that since Hermite cubics match the function value and $1^{\text{st}}$ derivative at grid points, the error scales quadratically with the distances between the query point and both the neighbouring grid points in each spacial dimension. For $\vec{x}$ a grid point, the evaluation of $\vec{\mathcal{X}}_{i}(\vec{\Psi}(\vec{x}),t-{\delta_{t}})$ in (12b) occurs at query points located $|\vec{\Psi}(\vec{x})-\vec{x}|=\mathcal{O}({\delta_{t}})$ and $\mathcal{O}({\delta_{x}})$ away from the closest grid points. Hence the $\mathcal{O}({\delta_{t}}^{2}{\delta_{x}}^{2})$ interpolation error.

For each $\mathcal{X}_{[\tau_{i-1},\tau_{i}]}$, we perform $\frac{\tau_{i}-\tau_{i-1}}{{\delta_{t}}}$ GALS steps. With each step accumulating the above error, the overall error for each submap is

After composition of the submaps, the global error due to submap evolution is

we see that for any number of time subdivisions, the composition of the submaps yields a globally third-order accurate method.

In the next section, we will discuss the numerical properties specific to this third-order method. Namely, we will look at the effect of the choice of coarse and fine grids, as well as the dynamic adaptation of time interval length and fine grid resolution.

In practice, the subdivision times $\tau_{i}$ are chosen adaptively and not predetermined. We define a measure of the error of the submap computation $\vec{\mathcal{X}}_{i}$ and take $\tau_{i}$ to be the first time this error exceeds some predetermined tolerance.

To evaluate this error, we use Lagrangian particles $p\in\{1,\dots,m\}$, initially distributed uniformly in $U$ at positions $\vec{y}_{p}^{\,0}$. These particles are independently evolved under $\vec{v}$ by solving the set of ODEs

We use RK3 integration scheme to solve (18) for particles $\vec{y}_{p}$. A measure of the accumulated interpolation error is given by

Equipped with the error measure $M_{1}$, we define $\tau_{i}$ as the first time at which the evolution of $\vec{\mathcal{X}}_{i}$ induces a representation error greater than some predefined tolerance $\mathcal{E}_{1}$, that is

At this time, we stop evolving $\vec{\mathcal{X}}_{i}$ and start computing the next submap. Notice that choosing $\mathcal{E}_{1}$ sufficiently small ensures that the higher frequencies omitted by the coarse grid representation do not become large. Indeed, even though the grid values of $\mathcal{X}_{i}$ may be very precise (for instance, by taking very small ${\delta_{t}}$), the Hermite cubic representation can still be bad due to a lack of resolution. In this case, projecting and composing on a fine grid leads to high error for the fine grid values. This is why the error measure we use takes into account the entire domain and controls error at off-grid points. Consequently, we may oversample $\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]}$ and compose it with $\vec{\mathcal{X}}_{[0,\tau_{i-1}]}$ onto a finer grid with controlled error.

The entire algorithm utilizes two grids: A coarse grid ${\bm{C}}$ of $N_{c}$ points per dimension where the $\vec{\mathcal{X}}_{i}$ submaps are evolved, and a fine grid ${\bm{F}}$ of $N_{f}$ points per dimension where the accumulated global maps $\vec{\mathcal{X}}_{[0,\tau_{i}]}$ are stored.

As seen in the previous subsection, $\mathcal{E}_{1}$ controls the importance of the omitted higher frequencies in the $\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]}$ submaps. These frequencies which cannot be represented on a $N_{c}$ grid arise either from the projection of the velocity onto the coarse grid or from advection effects on existing lower frequencies (e.g. shear effects). Similarly, when we perform the composition (13), advection effects will generate higher frequencies than what is originally present in $\vec{\mathcal{X}}_{[0,\tau_{i-1}]}$. In this case, the original resolution of $N_{f}$ might not be enough to resolve these frequencies. Therefore, we dynamically adjust the fine grid resolution $N_{f}$ to control the error made by omitting higher frequencies.

To do this, we define another error measure. This one is intended to evaluate the capacity of a grid ${\bm{F}}$ to properly represent some $C^{1}$ function $g$:

where the sup-norm is approximated by sampling the error at off-grid points.

We can then decide if a grid ${\bm{F}}$ is fine enough to accurately represent the updated global map $\vec{\mathcal{X}}_{[0,\tau_{i-1}]}\circ\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]}$ by evaluating $M_{2}({\bm{F}},\vec{\mathcal{X}}_{[0,\tau_{i-1}]}\circ\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]})$. If $M_{2}\geq\mathcal{E}_{2}$, for some predefined tolerance $\mathcal{E}_{2}$ we will refine the fine grid, replacing ${\bm{F}}$ by a grid ${\bm{F}}^{f}$ of $2N_{f}$ points per dimension. We will compute $\vec{\mathcal{X}}_{[0,\tau_{i}]}$ from (13) using ${\bm{F}}^{f}$ as ${\bm{F}}$.

Conversely, we want to take advantage of situations where the fine grid could be coarsened. To detect such situations, we define a coarser grid ${\bm{F}}^{c}$ of $N_{f}/2$ points per dimension and compute the $M_{2}$ error according to (21). If we have $M_{2}\leq\mathcal{E}_{2}$, it means that the coarser grid is sufficient to represent $\vec{\mathcal{X}}_{[0,\tau_{i-1}]}\circ\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]}$. In this case, we coarsen the fine grid, replacing ${\bm{F}}$ by a grid ${\bm{F}}^{c}$ of $N_{f}/2$ points per dimension and compute $\vec{\mathcal{X}}_{[0,\tau_{i}]}$ from (13) using ${\bm{F}}^{c}$ as ${\bm{F}}$.

Using such a dynamic grid has two main advantages. By refining the grid when needed, we ensure that the transformation is always well represented, and by coarsening it when possible, we reduce the computational time required to do a remapping step. Note that the redefinition of $\vec{\mathcal{X}}_{[0,\tau_{i}]}$ on a different grid is easily performed at low cost due to the Hermite cubic interpolant structure.

We summarize here the CM algorithm with some minor adjustments at the implementation level. We use a single variable $\vec{\mathcal{X}}_{0}$ for all global maps $\vec{\mathcal{X}}_{[0,\tau_{i}]}$, as the old global map is overwritten every time it is updated. We also use a single variable $\vec{\mathcal{X}}$ for all local submaps $\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]}$, as old submaps are no longer useful after the global map is updated. In the following, algorithm 1 describes the GALS method for evolving submaps, algorithm 2 updates the global map by submap composition on a dynamically selected fine grid.

We apply the characteristic mapping method to the following 2D vector field taken from [19]

with $A=16$ in the domain $[0,1]\times[0,1]$. This velocity creates a swirl centered at $(0.5,0.5)$. The swirl reaches its maximal deformation at $t=8$ and then unwinds itself. At $t=16$, the deformation induced in the first half of the period are completely undone and the transformation returns to identity.

The initial set is a circle of radius $0.15$ centered at $(0.5,0.75)$ represented by a level-set function. From times $t=0$ to $t=16$, this flow will stretch the circle into a thin swirl and return it back to its original state.

We use a $N_{c}=32$ grid for the advection and a dynamic-sized fine grid for the remapping. This fine grid starts at $N_{f}=32$ and is refined to a maximum resolution of $N_{f}=512$. The remapping and resizing tolerances were set to $\mathcal{E}_{1}=5\times 10^{-6}$ and $\mathcal{E}_{2}=10^{-4}$. Figure 2 shows the results.

The initial identity transformation is very well represented on the $32\times 32$ grid, but as time advances, the flow creates very fine structures. The transformation at $t=8$ would not be well represented on a $32\times 32$ grid. Indeed, the adaptive grid for $\vec{\mathcal{X}}_{0}$ is refined to $512\times 512$ at $t=8$, which is sufficient to represent correctly the sharp deformations induced by the vector field. Then as the swirl returns to its original shape, the grid coarsens. At $t=16$, the transformation only needs a $16\times 16$ grid to be correctly represented. We see that the initial circle is recovered at the final time, even though the advection was computed on the rather coarse $32\times 32$ grid and only a few remapping steps involved finer grid calculations.

We use a similar test to compare the characteristic mapping method with the standard GALS method [23]. We use the same initial set and transport it in the vector field (22) with $A=8$ until $t=16$, which corresponds to the set being stretched and returned to its original shape twice. We make four different tests by modifying the grid sizes. For the CM method, we use a $N_{c}=32$ for all tests but we set the maximal $N_{f}$ to different values. The maximum allowed fine grid for each test are $N_{f}=32,64,128,256$ respectively. For the GALS method, we fix the grid resolution to be the same as for the finest possible remapping grid of the CM method, i.e. $N_{g}=32,64,128,256$. Results are shown in figure 3.

We see that for a given grid size, the characteristic mapping method gives a better solution than the GALS method. This is in part due to the fact that the CM method starts from a linear initial condition $\vec{\mathcal{X}}_{0}=\vec{x}$ which can be exactly represented by Hermite cubics. Any feature of complexity higher than linear arise from the velocity field and accumulates at the rate of $\mathcal{O}({\delta_{t}})$ per step, hence it takes longer for significant high frequencies to appear. In contrast, the GALS method has a more complicated initial condition. The level-set function already contains $\mathcal{O}(1)$ linear, quadratic and cubic features. Under advection effects, these can quickly generate sharp features that the grid cannot represent. Therefore, when the maximum resolution level $N_{f}=N_{g}$ is fixed, CM will be more accurate than GALS, as shown in figure 3.

A particularly attractive feature of the CM method is that computation efficiency is optimized by the adaptive separation between coarse grid advection calculations and fine grid deformation representation. Indeed, in a short time interval, deformations caused by the advection are small, and errors are dominated by velocity integration in time. It is therefore most efficient to perform these steps on a coarse grid: they happen frequently, but computation is cheap on a coarse grid, and the coarse grid is sufficient to represent the small deformation. In fact, doing submap updates on a fine grid would be wasteful since fine grid representation errors would be negligible compared to error from velocity integration. For the global time map, the error is dominated by difficulties in representing sharp features which requires a high grid resolution. However, global map updates do not happen often, and we can afford to use a fine grid representation.

Table 1 compares the computational times for the GALS and CM method. For very coarse grids, the CM method is slower due to the overhead from processing the adaptivity. For the cases where the representation grid is finer than the advection grid, the CM method is clearly faster. The efficiency aspect is studied in more detail in section 4.5.

The CM method generalizes directly to any number of dimentions. We apply the CM method to a 3D surface evolution problem. The initial surface is a sphere of radius $0.15$ centered at $(0.35,0.35,0.35)$ in the domain $[0,1]^{3}$. The deformation velocity is given by the following vector field taken from [19]

The advection of $\vec{\mathcal{X}}$ is done on a $N_{c}=16$ grid, and the remapping of $\vec{\mathcal{X}}_{0}$ is done on a fixed $N_{f}=128$ grid with a remapping tolerance $\mathcal{E}_{1}=10^{-4}$. We compare the results with and without remapping to demonstrate the benefits of the remapping step. We also compare the CM method to the GALS method computed on a $N_{g}=128$ grid. Results for $t=0$, $t=1$ and $t=2$ for the three cases are shown in figure 4.

The surface at $t=2$ is expected to be identical to the surface at $t=0$ (sphere) due to symmetries in the velocity (23). This vector field causes the sphere to be stretched along the $y=1-x$ plane and creates very fine structures that cannot be represented on the coarse $16^{3}$ grid. We see from figure 4(d) that without remapping, those structures are indeed not well represented and oscillations are observed on the scale of the coarse grid. These oscillations distort the surface for all subsequent times, and the final shape in figure 4(g) is significantly different from the initial sphere.

When using the remapping on the fine grid, the fine structures caused by the deformation field can be accurately represented as a result of storing $\vec{\mathcal{X}}_{0}$ on a $N_{f}=128$ grid. In figure 4(e), the width of the stretched surface is of the order of the fine grid’s cell width, therefore causing no major representation issues. At $t=2$ (figure 4(h)), the surface is visually identical to the initial sphere.

The GALS method also uses a $N_{g}=128$ grid and therefore uses comparable computational resources for spacial resolution as CM with remapping. The main differences between the two methods is that CM method is initialized with a trivial identity map and does not discretize the more complicated initial condition of the level-set function. The effects of this can be seen later in the simulation in figures 4(h) and 4(i), where GALS accumulates more error than CM.

In addition to testing the accuracy of the CM method with remapping, this experiment also demonstrates the computational efficiency of the method. Table 2 compares the time taken to compute the three tests of figure 4. As expected, the time taken by the CM method without remapping is small, but the results are not precise. The time taken by the CM method with remapping is approximately 5 times larger, but achieves the best result out of the three cases. The GALS method takes about 20 times longer than the CM method with remapping and fails to provide better results. This indicates that the separation of scales using the semi-group property of the characteristic maps provides significant improvements to computational efficiency.

This section examines the application of the CM method for the advection of arbitrary sets. We take the deformation field $\vec{v}(\vec{x},t)$ given by

with A=16 and where $L$ and $R$ are smooth weight functions defined by

This vector field causes the left region to swirl into two vortices and the right region to expand around the center of the domain. We apply this vector field to the Mandelbrot set and compute the advection of the set with the CM method. We use $N_{c}=32$ for the coarse grid, a fixed $N_{f}=1024$ for the fine grid and a remapping tolerance $\mathcal{E}_{1}=10^{-7}$. The results at times $t=0$ to $t=16$ are shown in figure 5.

Traditional advection schemes would have difficulties dealing with complex geometries such as the Mandelbrot set. With the CM method, the complexity of the initial set does not translate to difficulties in numerical simulation. Since the advected quantity is the pullback of the initial function by the deformation map and since we know the initial set to arbitrary precision, there is no difficulty capturing all the fine structures. This is a property of the interpolation structure of the method; since the characteristic map is evolved as a function defined everywhere in the domain, $\phi(\vec{x},t)=\phi_{0}(\vec{\mathcal{X}}(\vec{x},t))$ can be evaluated directly for any $\vec{x}$. Therefore, we achieve arbitrary subgrid resolution and preserve the information contained in the initial condition $\phi_{0}$ throughout the simulation. We observe indeed in figures 5(a) and 5(i) that the same resolution is available at the beginning and at the end of the simulation. In both figures, we show a zoomed view of a single cell in the fine grid to show that at the end we recuperate, with arbitrary subgrid resolution, the same set as the initial condition, as expected. This is a remarkable and unique property of the CM framework and has important implications in real-life applications. As illustrated by the above example, the spacial resolution of the solution to the advection equation is not affected by the resolution of the characteristic map. Indeed, the Mandelbrot set requires infinite resolution to represent whereas the map is computed on a $N_{f}=1024$ grid; $N_{f}$ is selected to provide a good approximation to the characteristic map $\vec{X}$, it only needs to resolve the dynamics of the flow. In practice, one might have a very finely sampled discretized initial condition. Traditional methods such as GALS need to numerically represent the initial condition accurately in order to maintain resolution throughout the simulation. If the dynamics of the flow contain less spacial features, this would be a waste. The CM method is more efficient in that the characteristic map uses only the necessary computational resources to represent the deformation. The initial condition, no matter how finely resolved, will be pulled-back by the deformation map to generate solutions of the sample resolution as the initial condition.

For a given velocity field $\vec{v}$, the corresponding characteristic map $\vec{\mathcal{X}}$ is a solution operator which maps any initial condition to the solution of the advection under $\vec{v}$. This implies that multiple initial conditions can be advected simultaneously using a single map $\vec{\mathcal{X}}$. Furthermore, these solutions remain mutually consistent throughout the entire simulation because they arise from the same deformation map. We present two tests showing this feature.

The first test involves open curves. These are represented implicitly using two functions: a level-set function and a mask function. For instance, a segment of a straight line can be represented using an affine function whose zero level-set is the line containing the segment, and a mask function which represents a region whose intersection with the line gives the desired segment. Figure 6 shows an example where we advect a circle along with three independent line segments forming a triple point. These objects are transported in the swirling velocity field (22) with $A=4$. Normally, these seven functions used to define the curves (three lines, three masks and one circle) need to be advected separately; one would need to solve seven advection equations independently, giving potentially mutually incoherent results. The CM method however, decouples the advection from the initial set, we can advect those seven functions at the same time using the same characteristic map. When plotting is required, the functions are evaluated at $\vec{\mathcal{X}}(\vec{x},t)$, incurring only the extra cost of seven function evaluations and an evaluation of the map on the plotting grid. Since the solutions are obtained from pullback by the same characteristic map, the coherence of the results is guaranteed (e.g. triple points are guaranteed to stay triple points).

The second test shows an example of the transport of a Mosaic pattern. We take $[0,1]^{2}$ as a periodic domain and subdivide it into multiple regions. This kind of initial condition arises for instance in the simulation of multiphase flows. The difficulty of these simulations resides in the multiple intersection points caused by the junction of more than two different fluids. These points are hard to represent using a single level-set function, but can be represented piecewise as multiple level-sets. Saye and Sethian [28, 27] propose an approach to multiphase flows using level-sets, while Da et al. [7] suggest to use multimaterial front tracking. In the context of the CM method, this is not a problem since we transform the whole domain and can deal with functions of any level of complexity afterwards. The different regions can simply be represented by a single function by multiplying a different constant to their respective indicator functions.

The velocity field used for this test is

This velocity generates a periodic flow map with period 2. We used $N_{c}=32$ and $N_{f}=512$. It took $65$ seconds to compute the $2048$ steps of this simulation on a single $3.0$GHz CPU, taking $0.008$ seconds for a regular step and $0.45$ seconds when remapping, which was necessary on average approximately every $12$ step. We also highlight three regions in the domain using dashed circles. These regions represent two triple points and one quadruple point. The dashed circles are transported in the flow as passive particles with the initial intersection points as starting position. This shows that the intersections obtained from CM do follow the right path.

Resolution of thin filaments and accurate tracking of junction points are both important elements of a multiphase flow simulation. From figure 7, we see that the CM method manages to track very thin filaments. If each level-set was advected separately using traditional methods, the errors might not be consistent across individual level-sets. This could result in important qualitative errors on the presence and location of thin filaments, persisting throughout the entire simulation. However, since all level-sets in the CM method are transported using the same deformation map, relationships between different regions are naturally preserved, therefore resolving and maintaining the thin filaments. This also leads to an accurate and consistent tracking of triple and quadruple junctions. As seen in figure 7, even though we can see some error in the position of the quadruple junction when deformations are large, this error is not amplified and does not persist throughout the simulation. We see indeed that the junctions return to the correct positions at time $t=2$. This also indicates that the time symmetry present in the velocity is better preserved in the numerical solution: the characteristic map deforms until $t=1$, then returns to the identity map at $t=2$. This would be difficult to achieve with traditional methods, working directly with the level-set functions as the sharp and convoluted shapes of the regions will incur high numerical dissipation, thereby moving or destroying junction points.

Several elements factor into the efficiency of an algorithm, notably the number of computations per step and memory usage necessary to reach a certain accuracy, and the parallelizability of the method. In the following, we give some estimates on the error and running time of the method. For simplicity, we will restrict ourselves to fixed fine grid $N_{f}$ and try to examine the effect on efficiency from varying $N_{c}$, $N_{f}$ and $\mathcal{E}_{1}$.

For GALS, the error is know to be $\mathcal{O}({\delta_{t}}N_{g}^{-2}+N_{g}^{-4})$. For the CM method, we proceed as follows.

First, we look at the error from each submap computation, $\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]}$. Let $L_{i}=\tau_{i}-\tau_{i-1}$. Each time step, the error from time integration and evaluation at foot point is of order $\mathcal{O}({\delta_{t}}^{2}N_{c}^{-2}+{\delta_{t}}^{4})$. Taking ${\delta_{t}}\leq N_{c}^{-1}$, we simplify this error to roughly $(\mu_{i}+\sigma){\delta_{t}}^{2}N_{c}^{-2}$, where $\mu_{i}$ is a constant of the order of $D^{4}_{x}\vec{X}_{[\tau_{i-1},\tau_{i}]}$ coming from interpolation errors and $\sigma\approx D^{4}_{t}\vec{v}$ from integration errors. The total submap grid errors is then

There is also a submap representation error when evaluating the Hermite cubic $\vec{\mathcal{X}}_{[\tau_{i-1},\tau_{i}]}$ at the remapping step. This error exists even if grid values are exact and is due to the limited resolution available for a Hermite cubic on a coarse grid:

Finally, at each remapping step, we incur an error from evaluating the global map $\vec{\mathcal{X}}_{[0,\tau_{i-1}]}$:

where $\nu_{i}$ is on the order of $D_{x}^{4}\vec{X}_{[0,\tau_{i-1}]}$ (note that $\nu_{1}=0$).

We let $\mu$ and $\nu$ be the average of $\mu_{i}$ and $\nu_{i}$ over all remapping steps. Assuming there are $m$ remapping steps between $t=0$ and $t=T$, the global error is roughly

Both $\mu$ and $m$ are influenced by the choice of $\mathcal{E}_{1}$. If $\mathcal{E}_{1}$ is small, remapping happens more often, making $L_{i}$ shorter. Since $m$ correlates inversely with $L_{i}$, it would increase. On the other hand, $\mu$ can be assumed to correlate directly with $L_{i}$, and hence would decrease. We make the following assumptions and estimates.

We suppose that $\mu_{i}\approx BL_{i}$ for some constant $B$, that is, for short times, higher frequencies not resolved by the coarse grid accumulate linearly in time. It is also implied from this that $\nu\approx BT$, although this may be an overestimate. For small times, we can assume that sharp features in the deformation grow linearly with time. However, in global time, this might not always be true as seen in the swirl test, where the periodic velocity undoes the sharp features it created in the first half of the period. We will however assume that the estimate is generally of correct order of magnitude. From these assumptions, we reach the following estimate for the global error:

where $L$ is a function which scales directly with $\mathcal{E}_{1}$ and inversely with ${\delta_{t}}^{3}+N_{c}^{-4}$. Details of this calculation can be found in Appendix A.

Closer examination of this formula reveals three main regimes of the magnitude of the error depending on the choice of $\mathcal{E}_{1}$. When $\mathcal{E}_{1}$ is taken on the order of the local truncation error for submap evolution, i.e. $\mathcal{E}_{1}=\mathcal{O}({\delta_{t}}^{4}+{\delta_{t}}N_{c}^{-4})$, the remapping step happens every or almost every submap update step. The resulting global error is $\mathcal{O}({\delta_{t}}^{2}N_{c}^{-2}+{\delta_{t}}^{3}+N_{c}^{-4}+{\delta_{t}}^{-1}N_{f}^{-4})$. The other extreme is when $\mathcal{E}_{1}$ is of the order of global trunction error, i.e. $\mathcal{E}_{1}=\mathcal{O}({\delta_{t}}^{3}+{\delta_{t}}N_{c}^{-2}+N_{c}^{-4})$. In this case, barely any remapping occurs. The resulting error is $\mathcal{O}({\delta_{t}}^{3}+{\delta_{t}}N_{c}^{-2}+N_{c}^{-4}+N_{f}^{-4})$, with the $N_{f}$ term appearing if at least one remapping step occurred. Lastly, in the intermediate regime, with $\mathcal{E}_{1}$ roughly $\mathcal{O}({\delta_{t}}^{2}N_{c}^{-2})$, we should observe errors of order $\mathcal{O}({\delta_{t}}^{3}+N_{c}^{-4})$.

In order to compare the efficiency of GALS and CM, we need the following estimates on computational time.

For the GALS method, the cost consists of tracing back the footpoints, and then evaluating the interpolant at those locations.

for some constants $C_{1}$ and $C_{2}$, and where $d$ denotes the number of dimensions and $N_{g}$ is the number of grid cells for each spatial dimensions.

For the CM method, we also need to trace back the footpoints and interpolate the function, but that is done on the coarse $N_{c}^{d}$ grid, once for each of the $d$ coordinate functions. Additionally, we need to advect the test particles and do remapping steps when required. The computational cost is

where $\rho$ is the number of particles per cell. Recall that $L$ scales directly with respect to $\mathcal{E}_{1}$, so picking small $\mathcal{E}_{1}$ increases the frequency of remapping steps. This is the most costly operation of the method and hence strongly affects the efficiency. Roughly, the cost of CM is

In order to illustrate the behavior of error, computational time and efficiency of the algorithms, we ran several numerical experiments using the swirl test presented in section 4.1 with $A=8$ and a final time $T=16$. We do not use a dynamic fine grid and set $N_{f}=N_{g}$, ${\delta_{t}}=N_{f}^{-1}$, and fix $N_{c}=32$ for all computations. For values of $\mathcal{E}_{1}$ ranging from $10^{-4}$ to $10^{-7}$, we test the error and running times for different values of $N_{f}$. The results are show in figure 8.

In figure 8(a), we observe that when $\mathcal{E}_{1}$ is much smaller than the global error, we see third order convergence with respect to ${\delta_{t}}$ as described earlier ($N_{c}^{-2}$ was small enough for $N_{c}^{-2}{\delta_{t}}^{2}$ terms not to appear). However, when $\mathcal{E}_{1}$ is close to the global error, almost no remapping step occurs and the error scales linearly in ${\delta_{t}}$.

In figure 8(b), we can observe the effect of $\mathcal{E}_{1}$ on running time. When $\mathcal{E}_{1}$ is relatively large, the running time is dominated by submap time stepping: computational time scales linearly with ${\delta_{t}}^{-1}$. However, when $\mathcal{E}_{1}$ is picked very small, the cost of remapping becomes more obvious. Indeed, we see that when picking $\mathcal{E}_{1}=10^{-7}$, for large $N_{f}$, running time scales like ${\delta_{t}}^{-1}N_{f}^{2}$.

These results suggest that when a prescribed resolution on the velocity is chosen by fixing $N_{c}$, and a target error $\mathcal{E}_{1}$ for the advection by the approximated velocity is given, one can find an optimal fine grid size $N_{f}$ minimizing both error and computational time.

Lastly, the CM method is easily parallelizable. Multiple transported interfaces can be advected by evolving a single characteristic map. Moreover, since the evolution of the characteristic map is done using the local GALS scheme, grid value updates are fully explicit and use compact localized stencils. The same applies to the remapping step due to it being essentially composed of two Hermite interpolant evaluations. It is not the aim of this paper to investigate the parallel implementation of the CM method, however work in this direction has the potential for creating a fast and accurate parallelized solver for the linear advection equation.