A Fully Lagrangian Advection Scheme

Abstract

A numerical method for passive scalar and self-advection dynamics, Lagrangian rearrangement, is proposed. This fully Lagrangian advection algorithm introduces no artificial numerical dissipation or interpolation of parcel values. In the zero-viscosity limit, it preserves all of the Casimir invariants associated with parcel rearrangement. In the two-dimensional case presented here, these invariants are arbitrary piecewise continuous functions of the vorticity and concentration fields. The initial parcel centroids are evolved in a Lagrangian frame, using the method of characteristics. At any time this Lagrangian solution may be viewed by projecting it onto an Eulerian grid using a rearrangement map. The resulting rearrangement of initial parcel values is accomplished with a weighted Bresenham algorithm, which identifies quasi-optimal, distributed paths along which chains of parcels are pushed to fill in nearby empty cells. The error introduced by this rearrangement does not propagate to future time steps.

Appendices

Appendix 1: Weighted Bresenham Algorithm

In our weighted version of the Bresenham algorithm, the choice of the next point in the path depends on the weight of the eligible neighbouring cells. Figure 16 shows eight cases depending on the location of the destination point. One seeks a neighbouring cell in the general direction of the path with the lowest weight. If the slope is zero, the cells to the north-east, east, and south-east of the current cell are searched. The north-east and east cells are searched if the slope is between zero and one, the north, north-east, and east cells are searched if the slope is one, and the north and north-east cells are searched if the slope is greater than one. If the slope is infinity, the cells to the north-west, north, and north-east of the current cell are searched. Should the weight of two or more cells be the same, the original Bresenham algorithm will be the tie-breaker. The following is Asymptote code (a vector graphics language for technical drawing [7, 22]), along with eight sample output paths, for our weighted Bresenham algorithm (optimized so that all cases are mapped to the first quadrant).

Fig. 16

Paths of local minimal weight from the center circle to eight edge circles, as determined by the weighted Bresenham algorithm and the indicated numerical parcel weights

Appendix 2: Termination of Weighted Bresenham Algorithm

Theorem 1

(Termination of weighted Bresenham) The weighted Bresenham algorithm produces a finite path between any two points on a regular lattice. For a unit square lattice, at most $\left\lceil 1.82x \right\rceil$ steps are needed to connect two points that are a distance $x$ apart.

Proof

Let $A$ and $B$ be given on the grid. We want to find the desired path between them. For simplicity, we assume that $B$ is inside or on the boundary of the first octant with respect to $A$. Without loss of generality, consider a unit square lattice. As described before, depending on the position of $B$, we have either two or three choices for the next cell to be included in the path. If one of these choices is the grid point $B$, we are done; the algorithm terminates after choosing $B$. Otherwise, in choosing one of the immediate neighbours of the current cell, we will take a step of size 1 or $\sqrt{2}$. We must show that there exists a fixed number $\delta > 0$ such that the distance to the point $B$ in each step is always reduced by at least $\delta$. The algorithm will then terminate in a finite number of steps.

Case (i): Assume $B$ lies on the same horizontal line as $A$, so that the slope of the line from $A$ to $B$ is zero, ($m=0$). In this case (Fig. 17), the next point in the path is one of the points $C$, $D$, or $F$. If we choose $D$, then since $\overline{DB}=\overline{AB}-1$, a step of 1 is taken toward $B$. Suppose instead that we choose $C$. On letting $x=\overline{AB} \ge 2$ and $\delta =x-\overline{CB} =1+\overline{DB}-\overline{CB} < 1$ (since $\overline{DB} < \overline{CB}$), and noting that $\delta =\overline{AD}+\overline{DB}-\overline{CB}= \overline{CD}+\overline{DB}-\overline{CB} > 0$, we find that

$$\begin{aligned} (x-\delta )^2=\overline{CB}^2= \overline{DB}^2+1= (\overline{AB}-1)^2+1 = \overline{AB}^2-2\overline{AB}+2=x^2-2x+2, \end{aligned}$$

so that $-2x\delta +\delta ^2=-2x+2.$ We then deduce from $x\ge 2$ that $2-\delta ^2=2x(1-\delta ) \ge 4(1-\delta ).$ Thus $\delta ^2-4\delta +2 \le 0 \implies \delta \in \left[ 2-\sqrt{2},1\right) .$ The same argument of course also holds for the choice $F$. The distance reduction in this case is thus at least $2-\sqrt{2}$.

Fig. 17

Case(i): $m=0$

Case (ii): Assume that the slope of the line from $A$ to $B$ is 1. In this case (Fig. 18), the next point in the path will be one of the points $C$, $D$, or $F$. Here $\overline{AB} = \overline{DB}-\sqrt{2}$. If we choose $D$, we take a step of size $\sqrt{2}$ toward $B$. Suppose instead that we choose $C$. We see that $\overline{CH}=\overline{AH}=1/\sqrt{2}$. On letting $x=\overline{AB}$, we obtain

$$\begin{aligned} (x-\delta )^2=\overline{CB}^2=\overline{HB}^2+\overline{CH}^2 =\left( x-\frac{1}{\sqrt{2}}\right) ^2+\frac{1}{2}= x^2-\sqrt{2}x+1. \end{aligned}$$

Thus $-2x\delta +\delta ^2=-\sqrt{2}x+1.$ We know that $x\ge 2\sqrt{2}$ since $B$ is not one of the choices, so $\delta ^2-1=x(2\delta -\sqrt{2})\ge 2\sqrt{2}(2\delta -\sqrt{2})=4\sqrt{2}\delta -4.$ Now $\delta ^2-4\sqrt{2}\delta +3 \le 0 \implies \delta \ge 2\sqrt{2}-\sqrt{5}$. The same argument of course also holds for the choice $F$. The distance reduction in this case is thus at least $2\sqrt{2}-\sqrt{5}$.

Fig. 18

Case(ii): $m=1$

Case (iii): Assume $B$ lies inside the first octant ($0<m<1$). In this case, Fig. 19, the next point in the path is one of the two points $C_1$ or $C_2$. Let $x=\overline{AB}$. Notice that $x \ge \sqrt{5}$ since $B$ is not one of the points $C_1$ or $C_2$. Let $C$ be the point ($C_1$ or $C_2$) that is selected, drop the perpendicular $\overline{CH}$ to $\overline{AB}$. Let $y=\overline{CB}$ and $z=\overline{AC}$ and note that $z=1$ if $C=C_1$ and $z=\sqrt{2}$ if $C=C_2$. Since $0 < \angle CAH < \pi /4$, we know that $\overline{AH}/z > 1/\sqrt{2}$. On letting $\delta =x-\overline{CB}$, we find that

$$\begin{aligned} (x-\delta )^2&= \overline{CB}^2=\overline{HB}^2+\overline{CH}^2 = (\overline{AB}-\overline{AH})^2+\overline{AC}^2-\overline{AH}^2 \\&= x^2-2x\overline{AH}+z^2 < x^2-2x \frac{z}{\sqrt{2}}+z^2. \end{aligned}$$

Thus $-2x\delta +\delta ^2 < -\sqrt{2}xz+z^2,$ so that $z^2-\delta ^2 > x(\sqrt{2}z-2\delta ) \ge \sqrt{5}(\sqrt{2}z-2\delta ).$ Hence, $\delta ^2 -2\sqrt{5}\delta +\sqrt{10} z-z^2 < 0$. For $z=1$ this implies that $\delta > \sqrt{5}-\sqrt{6-\sqrt{10}}$ and for $z=\sqrt{2}$ this implies that $\delta > \sqrt{5}-\sqrt{7-2\sqrt{5}}$.

Fig. 19

Case(iii): $0<m<1$

In all cases, the distance between the point $B$ and the new included point is thus always less than $\overline{AB}$ by an amount

$$\begin{aligned} \delta&= \mathop {\mathrm{min}}\nolimits \left\{ 1, \sqrt{2}, 2-\sqrt{2}, 2\sqrt{2}-\sqrt{5}, \sqrt{5}-\sqrt{6-\sqrt{10}}, \sqrt{5}-\sqrt{7-2\sqrt{5}}\right\} \\&= \sqrt{5}-\sqrt{6-\sqrt{10}} > 0.551. \end{aligned}$$

That is, at most $\left\lceil 1.82\overline{AB}\,\right\rceil$ steps will be required to reach the point $B$. $\square$

Appendix 3: Multiple-Hole Expected Value

The expected number of holes in a shell $k$ (with $8k$ cells) containing a hole is

$$\begin{aligned} \left<j\right>= \frac{\sum _{j=1}^{8k} j {8k\atopwithdelims ()j} \left( 1-\frac{1}{e} \right) ^{8k-j} \left( \frac{1}{e}\right) ^j}{\sum _{j=1}^{8k} {8k\atopwithdelims ()j} \left( 1-\frac{1}{e} \right) ^{8k-j}\left( \frac{1}{e}\right) ^j}=\frac{\sum _{j=1}^{8k} j {8k\atopwithdelims ()j} r^j}{\sum _{j=1}^{8k} {8k\atopwithdelims ()j} r^j}, \end{aligned}$$

where $r=1/(e-1)$. Since the probability of not having a hole in $8k$ cells is $(1-1/e)^{8k}$, we get

$$\begin{aligned} \sum _{j=1}^{8k} {8k\atopwithdelims ()j} \left( 1-\frac{1}{e} \right) ^{8k-j}\left( \frac{1}{e}\right) ^j=1-\left( 1-\frac{1}{e}\right) ^{8k}, \end{aligned}$$

which on multiplying both sides by $(r+1)^{8k}=(1-1/e)^{-8k}$, can be written as

$$\begin{aligned} \sum _{j=1}^{8k} {8k\atopwithdelims ()j} r^j=(r+1)^{8k}-1. \end{aligned}$$

On differentiating both sides with respect to $r$ and then multiplying by $r$, it follows that

$$\begin{aligned} \sum _{j=1}^{8k} {8k\atopwithdelims ()j} jr^j=8kr(r+1)^{8k-1}. \end{aligned}$$

The expected number of holes thus evaluates to

$$\begin{aligned} \left<j\right>=\frac{\displaystyle 8kr(r+1)^{8k-1}}{(r+1)^{8k}-1} =\frac{\displaystyle 8k\left( 1-\frac{1}{r+1}\right) }{1-(r+1)^{-8k}} =\frac{\displaystyle 8k\left( \frac{1}{e}\right) }{\displaystyle 1-\left( 1-\frac{1}{e}\right) ^{8k}}. \end{aligned}$$

This of course is just the probability $1/e$ of finding a hole times the number of holes, conditioned on the probability that the shell contains at least one hole.