Welcome to mesh-sweeper!

This readme will give a brief introduction into what mesh sweeping is as well as the potential applications of it. A basic understanding of linear algebra is useful here!

What is Mesh Sweeping?

Solving Linear Systems with Mesh Sweeping

Fundamentally we want to solve a linear system; that is how to find $x\in\mathbb{R}^n$ such that $Ax=b$, given we know both $A\in\mathbb{R}^{n\times n}$ and $b\in\mathbb{R}^n$. There are many ways to solve this including exact solutions using Gaussian elimination for example or inexactly using Krylov methods for example. There are many advantages to both approaches with exact solutions generally having higher computational complexity requirements. For example, Gaussian elimination requires $\mathcal{O}(n^3)$ operations whereas Krylov methods can reach $\mathcal{O}(n)$ operations with specific matrix structures. Mesh sweeping is a variation of Gaussian elimination called Forward substitution. This method relies on the matrix $A$ being lower triangular but reduces the computational complexity to $\mathcal{O}(n^2)$.

Sorting Linear Equations + Transport

The difficulty of mesh sweeping however is not with the matrix solver but rather the ordering of the linear equations that form the matrix $A$. Generally, when you apply a numerical method to a partial differential equation (PDE), you end up with a system of equations to solve. In lots of cases, this system of equations depends on some spatial-temporal decomposition. For example if you consider transport over a spatial-temporal domain, $\Omega\times[0, T]\subset\mathbb{R}\times\mathbb{R}$, then you can decompose $\Omega$ and $[0, T] into intervals (potentially uniformly) and model the change over time of your modelled quantity in each of the spatial regions. We can generate an approximate linear equation for the change in this quantity and hence can combine all such equations for each spatial regions together to form $A$. As mentioned, this is where the magic happens, how do we determine the order in which we place these equations?

We can do this by taking advantage of the structure of the PDE. Each region in the transport setting has flow into and out of itself on its two boundaries. At each time in the temporal decomposition, flow on a boundary is either travelling into the region or out of the region. This is something we can harness to generate an algorithm. If the transport direction is fixed over the whole spatial-temporal domain, we get flow either to the left or the right - there is a natural ordering of the regions. This allows us to say that each region solely depends on its upstream neighbour. In terms of a linear equation, we state that the change in a region's modelled quantity is equal to what is already there plus what flows in from the upstream neighbour minus what leaves to the downstream neighbour. If the downstream flow is proportional to the current quantity, then the result is a matrix $A$ with entries on the diagonal and sub-diagonal. This is lower triangular and we can sweep.

Limitations + Current Best Practices

The downside here is that our example is a little contrived. In one spatial dimension, this is a trivial consequence with the given setup. In two- or more dimensions this is a hard problem. For example, in two dimensions, if a region is non-convex (maybe c-shaped) then flow will leave this region and re-enter at a later time. This generates a cycle in the mesh with the effect that our matrix $A$ can't ever be lower triangular. There are certain meshes with nice structures that one can use to prevent this as demonstrated in algorithms such as KBA for regular cuboidal meshes but for general polygonal/polytopal meshes, this is a little more difficult. The current best algorithm is the "number of upstreams" algorithm, but we further present a novel algorithm for a new class of structured yet irregular meshes.

Example usage here...