Note

This page was generated from a jupyter notebook.

TVD advection solver and related functions

(Greg Tucker, first version April 2023)

This notebook reviews Landlab tools for solving advection equations numerically. The tools include the AdvectionSolverTVD component and several related functions. AdvectionSolverTVD provides a second-order, total variation diminishing (TVD) numerical finite-volume solution to the 2D advection equation. The component requires a structured grid (raster or hex). The component can advect any node-based field or array. The user needs to provide a link-based advection velocity field.

Overview of new functions

Second-order TVD methods require evaluation of gradients “upwind” of a given location, where “upwind” means in the direction opposite to the advection direction. To facilitate this, Landlab’s components.advection package includes the following functions:

• find_upwind_link_at_link: for each link, identifies the adject parallel link that lies in the upwind direction. This function makes use of the raster- and hex-based property parallel_links_at_link, which is a (number-of-links x 2) data structure that records the IDs of adjacent parallel links that lie in the tail-wise (elements [:,0]) and head-wise (elements [:,1]) directions. The find_upwind_link_at_link function is called directly by AdvectionSolverTVD, but it is also available as a public function for custom-built applications.

• upwind_to_local_grad_ratio: for each link, calculates the ratio of gradient of given node scalar at the upwind link (if it exists) to the gradient at local link. For links that lack an upwind link (e.g., upwind would be beyond the grid boundary), or those where the local gradient is zero, a value of 1 is assigned. (Note that links along a grid boundary are considered not to have an upwind link.) Again, the function is called directly by AdvectionSolverTVD, but it is also available as a public function for custom-built applications.

TVD methods also require the use of a “flux limiter” function, which takes the upwind-to-local gradient ratio as an input, and returns a weighting factor that sets the relative contribution of high-order to low-order terms (see below). Many flux limiter functions have been proposed. As of this writing, Landlab’s advection package includes just one: flux_lim_vanleer, which can be imported from landlab.advection.flux_limiters. As of this writing, AdvectionSolverTVD uses the Van Leer flux limiter (future versions might provide other flux limters as user options).

Theory

The AdvectionSolverTVD and related support functions provide second-order TVD models for 2D advection. Here we briefly review the theory using a 1D advection equation with state variable \(z\):

\[\frac{\partial z}{\partial t} = -u \frac{\partial z}{\partial x}\]

where \(t\) is time, \(x\) is spatial position in 1D, and \(u\) is advection velocity (positive in the \(x\) direction). This type of equation comes in up in numerous contexts in Earth and environmental sciences. The right-hand side is an advection term; often such terms appear in combination with other terms, such as diffusion or source. (See for example Slingerland and Kump Mathematical Modeling of Earth’s Dynamical Systems (2011)).

To learn more about numerical solutions to advection problems, I highly recommend a series of online videos by Dr. Hilary Weller in the Meteorology Department at the University of Reading.

Discretizing the advection equation

A common discretization involves solving \(z\) at a set of nodes, which we’ll indicate with a subscript \(j\), so that \(z_j\) refers to the value of \(z\) at node \(j\). We’ll use superscripts to denote time step, so that \(z_j^n\) means the value of \(z\) at node \(j\) and time step \(n\). Using a first-order forward-difference time discretization, the discretized 1D advection equation can be written as:

\[\frac{z_j^{n+1}-z_j^n}{\Delta t} = -u \frac{z_{j+1/2} - z_{j-1/2}}{\Delta x}\]

where \(\Delta x\) is the spacing between nodes. Note that we’ve omitted the time-step superscript on the right side; for purposes of this example, we’ll assume it is \(n\), indicating an explicit solution, where not otherwise noted.

We can rearrange this equation to solve for the future quantity as:

\[z_j^{n+1} = z_j -c( z_{j+1/2} - z_{j-1/2} )\]

where \(c\) is a dimensionless parameter known as the Courant number:

\[c = \frac{u\Delta t}{\Delta x}\]

Note that \(c\) has the same sign as \(u\).

Defining the “in between” values

One of the challenges in solving advection problems discretized as above is defining what \(j+1/2\) or \(j-1/2\) means. If we imagine that each node sits inside a cell of width \(\Delta x\), then index \(j+1/2\) refers to the interface between the cells that contain nodes \(j\) and \(j+1\). The sketch below shows two nodes (indicated by a dot), labeled \(j\) and \(j+1\), with the ^ mark showing the interface at location \(j+1/2\):

-------------
|  .  |  .  |
-------------
   j  ^ j+1
    j+1/2

To solve the equation, we need a way of mapping values from one or both adjacent nodes onto the interface between them. One approach is to use a first-order upwind method, in which we assign to location \(z_{j+1/2}\) either the value \(z_j\) or \(z_{j+1}\), depending on which is “upwind” (in the direction opposite to that implied by \(u\)). If \(u>0\), we would assign \(z_j\); if \(u<0\), we would assign \(z_{j+1}\). Mathematically, this can be written:

\[z_{j+1/2} = z_j \text{ if } c\ge0, z_{j+1} \text{otherwise}\]

Second-order Total Variation Diminishing scheme

The problem with first-order upwind schemes is that they tend to produce excessive smoothing. One solution is to use a linear second-order advection scheme. One such scheme is the well-known Lax-Wendroff method, which assigns a weighted average of \(z_j\) and \(z_{j+1}\) based on the Courant number:

\[z_{j+1/2} = \frac{1}{2}(1+c)z_j + \frac{1}{2} (1-c)z_{j+1}\]

However, linear second-order methods generate spurious oscillations in the solution. Lax-Wendroff solutions tend to be smooth ahead of a discontinuity but generate oscillations behind it. The Warming and Beam method:

\[z_{j+1/2} = \frac{1}{2}(3-c)z_j - \frac{1}{2} (1-c)z_{j-1}\]

has the opposite problem, with smoothness behind but oscillations ahead.

These oscillation represent an increase in the total variation, which in this discrete formulation is the sum of \(|z_j - z_{j+1}|\), i.e., the sum of all variations in adjacent values of \(z\).

To ensure Total Variation Diminishing, we can calculate \(z_{j+1/2}\) as a weighted average of a high-order flux \(z_H\) and a low-order flux \(z_L\):

\[z_{j+1/2} = \Psi_{j+1/2} z_H + (1-\Psi_{j+1/2}) z_L\]

The low-order flux is just the linear upwind solution:

\[z_{L,j+1/2} = z_j \text{ if } c\ge0, z_{j+1} \text{otherwise}\]

For the high-order flux, Lax-Wendroff is a good choice:

\[z_{H,j+1/2} = \frac{1}{2}(1+c)z_j + \frac{1}{2} (1-c)z_{j+1}\]

\(\Psi\) is called a flux limiter function. The idea is to use \(z_H\) as much as possible. In particular, we want to have \(\Psi\) approach 1 where the solution is smooth, and approach zero where the solution is changing fast. This is a nonlinear method because \(\Psi\) depends on \(z\).

It turns out to be useful to define \(\Psi\) in terms of a variable called \(r\), which represents the ratio of the upwind gradient to the local gradient. For \(u>0\),

\[r_{j+1/2} = \frac{z_j - z_{j-1}}{z_{j+1}-z_j}\]

and for \(u<0\)

\[r_{j+1/2} = \frac{z_{j+2} - z_{j+1}}{z_{j+1}-z_j}\]

(Note: I have not found a source that describes what to do when the denominator is zero. From Campforts et al. (2017) it looks as they set \(r=1\) when that happens.)

One can cast \(\Psi\) as a function of \(r\). There are constraints on doing this to ensure TVD (summarized here). There are lots of these flux limiter functions. As of this writing, the one used in AdvectionSolverTVD is called Van Leer:

\[\Psi(r) = \frac{r+|r|}{1+|r|}\]

Practically speaking, one needs to compute gradients locally and at the upwind locations in order to calculate \(r\), and then calculate the flux limiter \(\Psi\) for each pair of adjacent nodes. After computing \(z_H\) and \(z_L\) for each pair of nodes, you calculate the weighted averages using \(\Psi\). Once that’s done, you calculate the solution to the difference equation above to get the value of \(z\) at each node for the new time step. As described below, the AdvectionSolverTVD takes advantage of Landlab’s node-link geometry to do these calculations.

Numerical implementation on a Landlab grid

In translating the second-order TVD scheme above onto a Landlab structured grid, it’s important to note that the discretized advection equation can be cast as a conservation-law problem. If the horizontal flux of a quantity \(z\) is \(q = uz\), then we can re-write the advection equation as:

\[\frac{\partial z}{\partial t} = -\frac{\partial q}{\partial x}\]

This equation has the form of a 1D mass (or volume) conservation equation with specific mass (or volume) flux \(q\).

For the sake of a concrete example, suppose we are interested in the tectonic advection of topography (with height = \(z(x,t)\)) through a fixed grid. The volume of material inside a grid cell \(j\) with projected surface area \(a_j\) is its average height, \(z_j\), times the area. Along each face \(k\) of the cell, material is being advected in at a rate of \(Q_k\) cubic meters per year. If \(Q_k\) is negative, it means that material is being carried out of the cell at that face. The rate of change of volume is:

\[\frac{\partial a_jz_j}{\partial t} = \sum_{k=1}^{N_j} Q_k\]

where \(N_j\) is the number of faces in cell \(j\) (4 if it is a rectangle, 6 if it is a hexagon, etc.). Let \(\lambda\) represent the length of a cell face. We can use this to express volume flux in terms of specific volume flux, \(q_k\), times \(\lambda\):

\[\frac{\partial a_jz_j}{\partial t} = \sum_{k=1}^{N_j} \lambda_k q_k\]

Assuming the cell areas don’t change with time,

\[\frac{\partial z_j}{\partial t} = \frac{1}{a_j} \sum_{k=1}^{N_j} \lambda_k q_k\]

Landlab already has a built-in function to calculate the right-hand side of this equation: calc_flux_div_at_node. The user just needs to provide an array of \(q\) values defined at links, and the function then uses the grid geometry (\(\lambda\) and \(a\)) to calculate the resulting rate.

(By the way, when we are dealing with a raster grid, the above equation gets even simpler: \(N_j=4\), \(\lambda_k\) is the grid spacing \(\Delta x\), and \(a_j=\Delta x^2\). For a hex grid, \(N_j=6\), and \(\lambda_k\) and \(a_j\) are the length of a hexagon face and the area of a hexagon cell, respectively. But we don’t need to worry about these implementation details, because the calc_flux_div_at_node takes care of them.)

For the topographic advection problem, \(q_k\) for any given cell face is the product of \(z_k\) and the face-perpendicular advection velocity \(u_k\). Because each cell face also has a matching link that connects the two nodes in question, the approach used by AdvectionSolverTVD is to calculate and store \(z_k\), \(u_k\), and \(q_k\) at the links. Among other things, this means that the user needs to provide face-perpendicular (link-parallel) advection velocities \(u_k\).

Setting up face-perpendicular advection velocities

For a RasterModelGrid, assigning face-perpendicular velocities is straightforward. Assuming you know the \(x\) and \(y\) velocity components, you can simply assign the \(x\) component to horizontal links, and the \(y\) component to vertical links. For example, if vel is an at-link velocity array (a required input; see below), grid is a RasterModelGrid, and ux and uy are spatially uniform velocity components, then you can assign the correct link-parallel velocities with:

vel[grid.horizontal_links] = ux
vel[grid.vertical_links] = uy

For a HexModelGrid, it is a bit more complicated, because the faces don’t all align with Cartesian directions. Fortunately, Landlab provides a mapping function called map_vectors_to_links to take care of this, illustrated by this example:

vel = grid.add_zeros('advection__velocity', at='link')
grid.map_vectors_to_links(ux, uy, out=vel)

This function will calculate and assign the vector projection of (ux, uy) onto each link. (The function works with raster grids too: it simply assigns ux to horizontal links and uy to vertical links.)

How AdvectionSolverTVD calculates its solution

To calculate a rate of vertical change at each grid node, the solution algorithm involves these steps:

1. If the advection field is changing in sign/direction, re-identify the link that lies directly upwind of every link. This calculation makes use of the find_upwind_link_at_link function, which in turn relies on the parallel_links_at_link data structure of raster and hex grids (which is only created when needed). A value of -1 is assigned to links that have no upwind neighbor.

2. Use the function map_node_to_link_linear_upwind to assign a “low order” value of the advected quantity at each link (corresponding to \(z_L\) in the preceding math).

3. For each link, use the function map_node_to_link_lax_wendroff to calculate a “high order” value, \(z_H\), of the advected quantity at each link.

4. Use the function upwind_to_local_grad_ratio to calculate the ratio \(r\) of upwind to local gradient in the advected quantity at each link (applying a value of unity where an upwind link does not exist or the local gradient is zero).

5. From the gradient ratio, use the function flux_lim_vanleer to calculate the weighting factor \(\Psi\).

6. Using the weighting factor, calculate a value of the advected quantity at the links using a weighted average of the low-order and high-order values. This corresponds to the “in between” values that were denoted as \(z_{j+1/2}\) in the preceding math.

7. Calculate a specific flux at each active link by multiplying the link-based value of the advected quantity by the link-perpendicular advection velocity. These are the \(q_k\) values.

8. Use the calc_flux_div_at_node function to calculate the flux divergence, and return the (negative) value of this as the time rate of change of the advected variable at each (non-perimeter) node. This corresponds to \(\partial q / \partial x\) in the 1D advection equation above, except that here it is implemented in 2D.

Calling update(dt) or run_one_step(dt) (they are synonyms) computes this rate of change and then updates the advected quantity using a forward-Euler step with the given time-step duration dt. If you want to get the rate of change without immediately updating the state variable \(z\), the component provides a calc_rate_of_change_at_nodes function. The function still requires a time-step duration dt as an input argument, because this is needed to evaluate the Courant number \(c\).

Parameters to AdvectionSolverTVD

Like nearly all Landlab components, AdvectionSolverTVD takes a ModelGrid object as its first argument. As of this writing, it must be a RasterModelGrid or HexModelGrid (unstructured grids don’t have parallel connected links, so there’s no natural “upwind” connectivity … BUT, if someone wants to figure out how to extend this capability to unstructured grids, please go ahead!)

Other parameters include:

• fields_to_advect: specifies which field or fields should be advected. This parameter can be given either as the name of a node-based field (string), or the array itself. Alternatively, if there is more than one field that should be advected with a given velocity, the parameter can be a list of either field names or arrays. For example, in a contaminant-transport problem the field to be advected would be contaminant concentration; for a topographic advection problem it would be the topographic elevation field; for water-wave propagation it would be the water-surface height.

• advection_direction_is_steady: a boolean that defaults to False. If you know that the advection direction is not going to change over time, set this to True to skip the overhead of re-identifying the upwind links at each time step.

Examples

Advection of a step function in 1d

Raster grid

import matplotlib.pyplot as plt
import numpy as np

from landlab import HexModelGrid, RasterModelGrid, imshow_grid
from landlab.components import AdvectionSolverTVD

# Parameters
u = 1.0  # advection velocity in x direction
c = 0.2  # Courant number
nnodes = 100  # number of grid columns
ntimesteps = 100  # number of time steps

# Setup
dx = 1.0 / (nnodes - 1)
dt = c * dx / u

# Create grid
grid = RasterModelGrid((3, nnodes), xy_spacing=dx)

# Create field for the advected quantity (here topographic elevation)
elev = grid.add_zeros("topographic__elevation", at="node")

# Initial condition: a step function from x = 0.4 to 0.6
step_exists = np.logical_and(grid.x_of_node >= 0.4, grid.x_of_node <= 0.6)
elev[step_exists] = 1.0

# Create field for advection velocity
vel = grid.add_zeros("advection__velocity", at="link")
vel[grid.horizontal_links] = u

# Instantiate component
# (Note: set advection_direction_is_stead to False so we can run it
# in reverse too)
advector = AdvectionSolverTVD(
    grid, fields_to_advect="topographic__elevation", advection_direction_is_steady=False
)

# For plotting convenience
midrow = np.arange(nnodes, 2 * nnodes, dtype=int)
x = grid.x_of_node[midrow]

plt.plot(x, elev[midrow], label="Step 0")
plt.xlabel("Distance")
plt.ylabel("Height")
for i in range(ntimesteps):
    advector.update(dt)
    if (i + 1) % 20 == 0:
        plt.plot(x, elev[midrow], label="Step " + str(i + 1))
plt.legend()

Now we’ll try reversing the direction and running it for some more time (that’s why we set advection_direction_is_steady=False above):

vel[grid.horizontal_links] = -u  # change direction

plt.plot(x, elev[midrow], label="Step 0")
plt.xlabel("Distance")
plt.ylabel("Height")
for i in range(ntimesteps):
    advector.update(dt)
    if (i + 1) % 20 == 0:
        plt.plot(x, elev[midrow], label="Step " + str(i + 1))
plt.legend()

…and presto, we’re back in the middle again (albeit with some numerical diffusion).

Raster grid, rotated coordinates

This is the same example as above, but with the grid’s long axis and advection direct both aligned with the \(y\) axis.

# Create grid
grid = RasterModelGrid((nnodes, 3), xy_spacing=dx)

# Create field for the advected quantity (here topographic elevation)
elev = grid.add_zeros("topographic__elevation", at="node")

# Initial condition: a step function from x = 0.4 to 0.6
step_exists = np.logical_and(grid.y_of_node >= 0.4, grid.y_of_node <= 0.6)
elev[step_exists] = 1.0

# Create field for advection velocity
vel = grid.add_zeros("advection__velocity", at="link")
vel[grid.vertical_links] = u

# Instantiate component
advector = AdvectionSolverTVD(
    grid, fields_to_advect="topographic__elevation", advection_direction_is_steady=True
)

# For plotting convenience
midrow = np.arange(1, 3 * nnodes - 1, 3, dtype=int)
x = grid.y_of_node[midrow]

plt.plot(x, elev[midrow], label="Step 0")
plt.xlabel("Distance")
plt.ylabel("Height")
for i in range(ntimesteps):
    advector.update(dt)
    if (i + 1) % 20 == 0:
        plt.plot(x, elev[midrow], label="Step " + str(i + 1))
plt.legend()

Hex grid in horizontal orientation

For a hex grid in quasi-1d, we have to be a bit more careful. We have to close the two long sides, because otherwise there would be flux outward to those boundaries along the angling links. We also have to multiply the magnitude the remaining horizontal flux by 1.5x to account for the fact that the face width is only 2/3 of the effective width of the cells. (This illustrates that a hex grid is not a good choice for these quasi-1d problems; it’s included here to show that the code works on a hex grid, and to demonstrate the use of the map_vectors_to_links function.)

# create grid
grid = HexModelGrid((3, nnodes - 1), spacing=dx)
grid.status_at_node[grid.y_of_node == 0.0] = grid.BC_NODE_IS_CLOSED
grid.status_at_node[grid.y_of_node > dx] = grid.BC_NODE_IS_CLOSED

# create fields
elev = grid.add_zeros("topographic__elevation", at="node")

# initialize fields
step_exists = np.logical_and(grid.x_of_node >= 0.4, grid.x_of_node <= 0.6)
elev[step_exists] = 1.0

# Create field for advection velocity
vel = grid.add_zeros("advection__velocity", at="link")
grid.map_vectors_to_links(1.5 * u, 0.0, out=vel)

# Instantiate component
advector = AdvectionSolverTVD(
    grid, fields_to_advect="topographic__elevation", advection_direction_is_steady=True
)

# misc
midrow = np.arange(nnodes - 1, (2 * nnodes) - 1, dtype=int)
x = grid.x_of_node[midrow]

plt.plot(x, elev[midrow], label="Step 0")
plt.xlabel("Distance")
plt.ylabel("Height")
for i in range(ntimesteps):
    advector.update(dt)
    if (i + 1) % 20 == 0:
        plt.plot(x, elev[midrow], label="Step " + str(i + 1))
plt.legend()

Advection of a Gaussian “bump” in 2d

We start by defining a function to make the initial bump, and a function to report the cumulative translation of the peak of the bump.

def make_gaussian_bump(grid, field, width):
    midx = 0.5 * np.amax(grid.x_of_node)
    midy = 0.5 * np.amax(grid.y_of_node)
    dist_sq = (grid.x_of_node - midx) ** 2 + (grid.y_of_node - midy) ** 2
    field[:] = np.exp(-0.5 * dist_sq / bump_width)

def calc_translation_distance(grid, elev, x0, y0, u, dt, nt):
    highest_node = np.argmax(elev)
    distance = np.sqrt(
        (grid.x_of_node[highest_node] - x0) ** 2
        + (grid.y_of_node[highest_node] - y0) ** 2
    )
    print("Velocity magnitude is", u)
    print("Duration is", dt * nt)
    print("Predicted distance is", u * dt * nt)
    print("Modeled distance is", distance)

Raster grid with advection to the east

Note that when we instantiate the component, we will leave its value of advection_direction_is_steady at the default of False. We do this because we are going to reuse the component several times with different advection directions, so we want it to recompute the upwind links.

# Parameters
nrows = 51
ncols = 51
ux = 1.0
uy = 0.0
c = 0.2
nnodes = 100
ntimesteps = 50
bump_width = 0.01  # width of Gaussian bump

# Setup
dx = 1.0 / (nrows - 1)
u = (ux * ux + uy * uy) ** 0.5
dt = c * dx / u

# create grid
grid = RasterModelGrid((nrows, ncols), xy_spacing=dx)

# create fields
elev = grid.add_zeros("topographic__elevation", at="node")
vel = grid.add_zeros("advection__velocity", at="link")

# initialize fields
make_gaussian_bump(grid, elev, bump_width)
maxnode = np.argmax(elev)
x0 = grid.x_of_node[maxnode]
y0 = grid.y_of_node[maxnode]
vel[grid.horizontal_links] = ux
vel[grid.vertical_links] = uy

# Instantiate component
adv = AdvectionSolverTVD(grid, fields_to_advect="topographic__elevation")

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(adv.grid, elev)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Test that we can run it back in the

Raster grid with advection to the north

This time we’ll show the change in value.

# Re-initialize and run again, this time with y-directed motion
ux = 0.0
uy = 1.0
vel[grid.horizontal_links] = ux
vel[grid.vertical_links] = uy

# Re-initialize the Gaussian bump
make_gaussian_bump(grid, elev, bump_width)
elev0 = elev.copy()

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(adv.grid, elev - elev0)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Raster grid with advection to the southwest

# Re-initialize and run again, this time with anti-x-directed motion
ux = -(2.0**0.5)
uy = -(2.0**0.5)
u = (ux * ux + uy * uy) ** 0.5

make_gaussian_bump(grid, elev, bump_width)
vel[grid.horizontal_links] = ux
vel[grid.vertical_links] = uy

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(adv.grid, elev - elev0)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Raster grid with a non-uniform advection field

This example tests the component with a non-uniform advection field. For this test, we will shear the Gaussian bump by having \(x\)-directed motion that accelerates in the \(y\) direction.

# Re-initialize and run again, this time with anti-x-directed motion
ux_max = 2.0
uy = 0.0
u = ux_max

make_gaussian_bump(grid, elev, bump_width)
vel[grid.horizontal_links] = (
    ux_max * grid.y_of_node[grid.node_at_link_tail[grid.horizontal_links]]
)
vel[grid.vertical_links] = 0.0

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(adv.grid, elev)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Case of divergent velocity

Here the velocity field diverges around \(x=1/2\), which produces a rift in the midst of our Gaussian bump: a bit like Iceland might be without the volcanoes …

# Identify places with positive vs. negative x-wise velocity
make_gaussian_bump(grid, elev, bump_width)
vel_is_neg = grid.x_of_node[grid.node_at_link_head] <= 0.5
vel[:] = 1.0
vel[vel_is_neg] = -1.0
vel[grid.vertical_links] = 0.0

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(adv.grid, elev)

Hex grid with tests of advection in several different directions

The last few examples use a hex grid. Note the use of map_vectors_to_links to assign link-parallel advection speeds to the links on the basis of their orientation relative to the advection direction and magnitude.

Case of right-ward motion

# Parameters
ux = 1.0
uy = 0.0
u = (ux * ux + uy * uy) ** 0.5

ncols = 51
nrows = int(ncols / (0.5 * 3.0**0.5))

# Setup
dx = 1.0 / (ncols - 1)
dt = c * dx / u

# Create grid
grid = HexModelGrid((nrows, ncols), spacing=dx, node_layout="rect")

# Create fields
elev = grid.add_zeros("topographic__elevation", at="node")
vel = grid.add_zeros("advection__velocity", at="link")

# Initialize fields
make_gaussian_bump(grid, elev, bump_width)
maxnode = np.argmax(elev)
x0 = grid.x_of_node[maxnode]
y0 = grid.y_of_node[maxnode]
grid.map_vectors_to_links(ux, uy, out=vel)

# Instantiate component
adv = AdvectionSolverTVD(grid, fields_to_advect="topographic__elevation")

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(grid, elev)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Case of upward (northward) motion

This time showing the differences after advection.

# Parameters
ux = 0.0
uy = 1.0

# Initialize fields
make_gaussian_bump(grid, elev, bump_width)
grid.map_vectors_to_links(ux, uy, out=vel)

# Remember the starting elevation field
elev0 = elev.copy()

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(grid, elev - elev0)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Case of leftward (westward) motion

# Parameters
ux = -1.0
uy = 0.0

# Initialize fields
make_gaussian_bump(grid, elev, bump_width)
grid.map_vectors_to_links(ux, uy, out=vel)

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(grid, elev - elev0)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Case of downward (southward) motion

# Parameters
ux = 0.0
uy = -1.0

# Initialize fields
make_gaussian_bump(grid, elev, bump_width)
grid.map_vectors_to_links(ux, uy, out=vel)

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(grid, elev - elev0)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Case of nne motion

# Parameters
ux = np.cos(np.degrees(60))
uy = np.sin(np.degrees(60))
u = np.sqrt(ux * ux + uy * uy)

# Initialize fields
make_gaussian_bump(grid, elev, bump_width)
grid.map_vectors_to_links(ux, uy, out=vel)

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(grid, elev - elev0)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Case of nnw motion

# Parameters
ux = -np.cos(np.degrees(60))
uy = np.sin(np.degrees(60))
u = np.sqrt(ux * ux + uy * uy)

# Initialize fields
make_gaussian_bump(grid, elev, bump_width)
grid.map_vectors_to_links(ux, uy, out=vel)

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(grid, elev - elev0)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Case of ssw motion

# Parameters
ux = -np.cos(np.degrees(60))
uy = -np.sin(np.degrees(60))
u = np.sqrt(ux * ux + uy * uy)

# Initialize fields
make_gaussian_bump(grid, elev, bump_width)
grid.map_vectors_to_links(ux, uy, out=vel)

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(grid, elev - elev0)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Case of sse motion

# Parameters
ux = np.cos(np.degrees(60))
uy = -np.sin(np.degrees(60))
u = np.sqrt(ux * ux + uy * uy)

# Initialize fields
make_gaussian_bump(grid, elev, bump_width)
grid.map_vectors_to_links(ux, uy, out=vel)

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(grid, elev - elev0)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Case with vertically oriented grid, rightward (eastward) motion

# Parameters
ux = 1.0
uy = 0.0
u = np.sqrt(ux * ux + uy * uy)
nrows = 51
ncols = int(nrows / (0.5 * 3.0**0.5))

# Setup
dx = 1.0 / (nrows - 1)
dt = c * dx / u

# Create grid
grid = HexModelGrid(
    (nrows, ncols), spacing=dx, node_layout="rect", orientation="vertical"
)

# Create fields
elev = grid.add_zeros("topographic__elevation", at="node")
vel = grid.add_zeros("advection__velocity", at="link")

# Initialize fields
make_gaussian_bump(grid, elev, bump_width)
elev0 = elev.copy()
maxnode = np.argmax(elev)
x0 = grid.x_of_node[maxnode]
y0 = grid.y_of_node[maxnode]
grid.map_vectors_to_links(ux, uy, out=vel)

# Instantiate component
adv = AdvectionSolverTVD(
    grid, fields_to_advect="topographic__elevation", advection_direction_is_steady=True
)

for _ in range(ntimesteps):
    adv.update(dt)
imshow_grid(grid, elev - elev0)
calc_translation_distance(grid, elev, x0, y0, u, dt, ntimesteps)

Tiny unit tests

It’s helpful to have small, simple unit tests. Here is a description of such a test, applied to single and multiple fields. A quantity advects from left to right, such that in the wake of the advection front the values should be equal (or close to) to the starting value at the left side of the grid. In the test, there are 5 core nodes; the idea is that after advection over a total distance of 3 nodes, the value at the middle node should be within a percent or so of the maximum starting value (unity in the first two tests, and for the 3rd test, which advects two different quantities called C1 and C2, it should be unity for C2 and 10 for C1).

grid = RasterModelGrid((3, 7))
C = grid.add_zeros("advected__quantity", at="node")
C[grid.x_of_node < 2.5] = 1.0
u = grid.add_zeros("advection__velocity", at="link")
u[grid.horizontal_links] = 1.0
print("Before:", C[grid.core_nodes])

advector = AdvectionSolverTVD(grid, advection_direction_is_steady=True)

for _ in range(15):
    advector.run_one_step(0.2)

print("After:", C[grid.core_nodes])

grid = RasterModelGrid((3, 7))
C = grid.add_zeros("only_advected__quantity", at="node")
C[grid.x_of_node < 2.5] = 1.0
u = grid.add_zeros("advection__velocity", at="link")
u[grid.horizontal_links] = 1.0
print("Before:", C[grid.core_nodes])

advector = AdvectionSolverTVD(
    grid, fields_to_advect="only_advected__quantity", advection_direction_is_steady=True
)

for _ in range(15):
    advector.run_one_step(0.2)

print("After:", C[grid.core_nodes])

grid = RasterModelGrid((3, 7))
C1 = grid.add_zeros("first_advected__quantity", at="node")
C1[grid.x_of_node < 2.5] = 10.0
C2 = grid.add_zeros("second_advected__quantity", at="node")
C2[grid.x_of_node < 2.5] = 1.0
u = grid.add_zeros("advection__velocity", at="link")
u[grid.horizontal_links] = 1.0
print("Before:")
print(" C1:", C1[grid.core_nodes])
print(" C2:", C2[grid.core_nodes])

advector = AdvectionSolverTVD(
    grid,
    fields_to_advect=["first_advected__quantity", "second_advected__quantity"],
    advection_direction_is_steady=True,
)

for _ in range(15):
    advector.run_one_step(0.2)

print("After:")
print(" C1:", C1[grid.core_nodes])
print(" C2:", C2[grid.core_nodes])

Summary

The AdvectionSolverTVD component provides a second-order, nonlinear, Total Variation Diminishing (TVD) solution to the 2D advection equation. The component can be used with any user-supplied node field as the quantity to be advected. The advection field can be unsteady and/or nonuniform. The component works with either a raster or hex grid.

References and further reading

Campforts, B., Schwanghart, W., & Govers, G. (2017). Accurate simulation of transient landscape evolution by eliminating numerical diffusion: the TTLEM 1.0 model. Earth Surface Dynamics, 5(1), 47-66.

Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. (2007). Numerical recipes 3rd edition: The art of scientific computing. Cambridge university press.

Slingerland, R., & Kump, L. (2011). Mathematical modeling of Earth’s dynamical systems: A primer. Princeton University Press.

Weller, Hilary (various dates). Short videos on numerical methods for atmospheric modelling. Available on YouTube.

---

Generated by nbsphinx from a Jupyter notebook.