FastFlow: GPU Acceleration of Flow and Depression Routing for Landscape Simulation

Input : Terrain cells $\backslash mathcal\{T\}$ with initial discharge (precipitation) $q\; =\; p$ and their donors
Output: Updated accumulated discharge

1 for $i\; \backslash leftarrow\; 1$ to $\backslash log\_2(|\backslash mathcal\{T\}|)$ do
2   foreach cell $c\; \backslash in\; \backslash mathcal\{T\}$ in parallel do
3     foreach donor $d$ of cell $c$ do
4       if $d$ is a leaf then
5         $q\_c\; \backslash leftarrow\; q\_c\; +\; q\_d$
6         remove edge $d\; \backslash rightarrow\; c$
7       end
8     else if $d$ has a single donor then
9       $q\_c\; \backslash leftarrow\; q\_c\; +\; q\_d$
10      donor of $c\; \backslash leftarrow$ donor of $d$
11    end
12  end
13  if all donors are removed then
14    tag $c$ as leaf
15  end
16 end
17 end

With the exception of some cells marked by the outflow mask, stream trees cover all of the terrain, which makes the design of an efficient parallel solution to Equation 1 non-trivial.

Fortunately, this problem can be cast as an instance of Parallel Tree Accumulation, which aims to compute the cumulative sum of quantities in tree nodes, from leaf to root. This can be achieved with straightforward pointer jumping, but this may result in a $\backslash mathcal\{O\}(n)$ span due to the required atomic operations. An alternative is to use an Euler tour [SAF05] along the tree edges combined with pointer jumping, which reduces the span to $\backslash mathcal\{O\}(\backslash log\; n)$, since the number of edges is less than the number of nodes. However, this method typically requires $2\; \backslash log\_2\; n$ iterations, as each edge is traversed twice.

Figure 3: A few iterations of a parallel rake-compress solution to flow routing. On each iteration, leaf nodes (in green) are pruned and recipients of cells with a single donor (in blue) are rerouted.

A still more efficient approach is rake-compress, which combines pointer jumping with leaf pruning to achieve a span of $\backslash mathcal\{O\}(\backslash log\; n)$. This method typically needs only $\backslash log\_2\; n$ iterations, making it faster in practice than the Euler tour. Sevilgen et al. [SAF05] provide a starting point in the form of a general rake-compress solution to graph accumulation problems for distributed architectures. This relies on the observation that leaves constitute half the nodes in a balanced binary tree. So, if we iteratively transfer accumulated discharge from a leaf to its recipient and then prune the leaf and do this in parallel across all leaves, we only need to iterate $\backslash log\_2\; n$ times where $n$ is the number of nodes in the terrain.

In general, river trees might be neither binary nor balanced. In the worst case, where a river tree is composed entirely of single donors, the iterations degrade to $\backslash mathcal\{O\}(n)$. To compensate, we use pointer jumping by re-routing the edges of cells whose donor itself has a single donor (parent). Effectively, the parent (donor) is skipped and an edge is formed directly from the grandparent (donor of the donor). Thereby, worst-case performance is restored to $\backslash log\_2(n)$ parallel iterations.

As shown in Algorithm 1, with an accompanying example in Figure 3, on each iteration and for each terrain cell $c$, we parse each donor $d$ in parallel, with up to 4 donors in our implementation. If $d$ is a leaf, we add its discharge to $c$ and remove it from the donor matrix. Alternatively, if $d$ has only a single donor, we again add its discharge to $c$ but in this case reconfigure the donor matrix so that $c$ instead points at the donor of $d$. Distinguishing between these cases requires a single lookup in the donor count list, e.g., leaves have a donor count of zero. Straightforward extensions of Algorithm 1 can be used for other flow accumulation functions, including those involving maxima or weighted sums of upstream values.

Note that some care must be taken in the implementation of Algorithm 1, which exhibits potential read-after-write (RAW) hazards when updating the discharge and recipients. An immediate solution is to write the updated $q$ and $r$ to temporary arrays and copy them back at the end of each iteration. Unfortunately, this forces the unnecessary transfer on each iteration of a large number of cells that are unmodified. Instead, we propose a per-cell ping pong scheme. We first allocate sufficient space to hold the state of two full terrains, $T\_A$ and $T\_B$. We then introduce an additional source buffer with one 8-bit integer per cell, initialized to 0. The sign of the integer at cell $c$ in the source indicates whether data should be read from $T\_A$ and written to $T\_B$, or vice-versa. When cells are updated,

the sign is flipped. To avoid RAW hazards within the same iteration, the remaining bytes of each source entry store the iteration at which $c$ was last updated: a thread will check the source to negate the read sign if the stored iteration is the current one. At the end of each iteration, the latest state of each cell is either in $T\_A$ or $T\_B$, as indicated by source. Hence, a final merge step is required before returning the result as a single tensor. Incorporating this optimization leads to a $14\backslash times$ speed up on a $1024\; \backslash times\; 1024$ terrain.