Lateral advection and diffusion

The second example is adapted from Cushman-Roisin and Beckers (2011), Introduction to Geophysical Fluid Dynamics p.196. It is an advection-diffusion problem that describes the transport and dispersal of a tracer in a 2D flow field. The problem is simple enough to have an analytical solution which makes it possible to directly compare results from the analytical model with those from the particle tracking model. The sensitivity of the result to the number of simulated particles is explored.

Background theory

The model describes the temporal evolution of a tracer in a 2D flow field with uniform velocity components \(u\) and \(v\), and uniform constant diffusion \(A\). Given these conditions, the 2D advection-diffusion equation:

\[\begin{equation} \frac{\partial c}{\partial t} + u \frac{\partial c}{\partial x} + v \frac{\partial c}{\partial y} = \frac{\partial}{\partial x} \left(A \frac{\partial c}{\partial x}\right) + \frac{\partial}{\partial y} \left(A \frac{\partial c}{\partial y}\right) \end{equation}\]

has the following analytical solution:

\[\begin{equation} c\left(x, y, t\right) = \frac{M}{4 \pi A t} e^{\frac{\left(x - u t\right)^{2} + \left(y - v t\right)^{2}}{4 A t}} \end{equation}\]

where \(c\) (g m\(^{\mathrm{-3}}\)) is the tracer concentration and \(M\) (g) is the amount of tracer released at time \(t = t_0\) and position \((x,y) = (0,0)\). For testing purposes we set \(u\) and \(v\) to 1 m s\(^{\mathrm{-1}}\), \(A\) = 10 m\(^{\mathrm{2}}\) s\(^{\mathrm{-1}}\) and \(M\) = 1 x 10\(^{\mathrm{4}}\) g. The solution is shown below, where it can be seen the tracer slowly disperses away from the centre of the release site, while the centre of mass is advected toward the top-left hand corner of the grid:

import warnings
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from matplotlib import animation, rc
from IPython.display import HTML

from pylag.mock import MockVelocityEddyDiffusivityDataReader

from pylag.processing.plot import create_figure, colourmap
from pylag.processing.utils import get_grid_bands


warnings.filterwarnings('ignore')

# Ensure inline plotting
%matplotlib inline

# Set up the spatial grid
xmin = -20.
ymin = -20.
xmax = 60.
ymax = 60.
n_x_points = 41
n_y_points = 41
x_grid = np.linspace(xmin, xmax, n_x_points, dtype=float)
y_grid = np.linspace(ymin, ymax, n_y_points, dtype=float)

# Grid edges for plotting
x_grid_edges = get_grid_bands(x_grid)
y_grid_edges = get_grid_bands(y_grid)

# Set up the time grid
t_min = 0.1
t_max = 25.0
time_step = 0.1
t_grid = time = np.arange(t_min, t_max + time_step, time_step)
n_times = len(t_grid)

# Compute the analytic solution
data_reader = MockVelocityEddyDiffusivityDataReader()
C = np.empty((n_times, n_x_points, n_y_points), dtype=float)
for t_idx, t in enumerate(t_grid):
    for x_idx, x in enumerate(x_grid):
        for y_idx, y in enumerate(y_grid):
            C[t_idx, x_idx, y_idx] = data_reader.get_concentration_analytic(t, x, y)


# Animate the concentration field
# -------------------------------

font_size = 15
fig, ax = create_figure(figure_size=(16., 16.), font_size=font_size)
ax.set_xlim((x_grid_edges.min(), x_grid_edges.max()))
ax.set_ylim((y_grid_edges.min(), y_grid_edges.max()))
ax.set_xlabel('x', fontsize=font_size)
ax.set_ylabel('y', fontsize=font_size)
ax.set_aspect('equal')

# Set colour map and limits
cmap = colourmap('DYE')
vmin = 0.0
vmax = 10.0

# Set up the plot
pcolormesh= ax.pcolormesh(x_grid_edges, y_grid_edges, C[0,:,:], vmin=vmin, vmax=vmax, cmap=cmap)
pos = ax.get_position().get_points()
cax = fig.add_axes([pos[1,0]+0.02, pos[0,1], 0.02, pos[1,1] - pos[0,1]])
cbar = fig.colorbar(pcolormesh, cax=cax)
cbar.ax.tick_params(labelsize=font_size)
cbar.set_label('C (kg m$\mathrm{^{-3}}$)', fontsize=font_size)

def animate(i):
    pcolormesh.set_array(C[i,:,:].ravel())
    return (pcolormesh,)

# Prevent the basic figure from being plotted too
plt.close(fig)

anim = animation.FuncAnimation(fig, animate, frames=n_times, interval=40, blit=True)

HTML(anim.to_html5_video())

In PyLag, particle trajectories are again calculated using a Random Displacement Model with the general form:

\[\begin{equation} dX_{j} = \left[u_{j}\left(\mathbf{x},t\right) + \dfrac{\partial D_{jk}\left(\mathbf{x},t\right)}{\partial x_{k}}\right]dt + \left(2 D_{jk}\left(\mathbf{x},t\right)\right)^{1/2}dW_{k}(t) \end{equation}\]

In this particular problem, the diffusion term is non-zero. However, as it is both constant in time and uniform in space, the pseudo-advection term that arises due to inhomogeneities in the diffusivity field drops out. This allows the above equation to be simplified. Written out in component form we have:

\[\begin{equation} \label{eqn:dx} dX = u\left(\mathbf{x},t\right) dt + \left(2 A\left(\mathbf{x},t\right)\right)^{1/2}dW_X(t), \end{equation}\]

\[\begin{equation} \label{eqn:dy} dY = v\left(\mathbf{x},t\right) dt + \left(2 A\left(\mathbf{x},t\right)\right)^{1/2}dW_Y(t), \end{equation}\]

\[\begin{equation} \label{eqn:dz} dZ = 0, \end{equation}\]

PyLag computes solutions to the above equations using a numerical method for solving stochastic differential equations which can be selected at run time by the user. In contrast to the previous example, particle motion now involves both a deterministic and stochastic component, and the particle tracking model must combine both to correctly model the temporal evolution of the tracer field.

Running the particle tracking model

As before, we first subclass PyLag’s DataReader. The example has been implemented in pylag.mock, and is called MockVelocityEddyDiffusivityDataReader. Again, we use an object of type StdNumIntegrator to manage the integration. Updates to particle positions are calculated using a 2D Millstein method.

import warnings
import numpy as np
from configparser import SafeConfigParser

import pylag.random as random
from pylag.numerics import StdNumMethod
from pylag.mock import MockVelocityEddyDiffusivityDataReader, MockTwoDNumMethod

warnings.filterwarnings('ignore')

# Seed the random number generator
random.seed(10)

# Start/stop times
time_start = 0.0 # Start time (s)
time_end = 25.0  # End time (s)
time_step = 0.1  # Time step (s)

# PyLag reads many of its configuration parameters from a config which we create here
config = SafeConfigParser()

# Set the coordinate system, which here is a simple cartesian coordinate system
config.add_section('SIMULATION')
config.set('SIMULATION', 'coordinate_system', 'cartesian')

# Parameters controlling the type of numerical integration scheme used
config.add_section("NUMERICS")
config.set("NUMERICS", "num_method", "standard")
config.set("NUMERICS", "iterative_method", "AdvDiff_Milstein_3D")
config.set("NUMERICS", "time_step_diff", str(time_step))

# Parameters controlling particle vertical motion, which we set to default values
config.set("SIMULATION", "depth_restoring", "False")
config.set("SIMULATION", "fixed_depth", "0.0")

# The data reader
data_reader = MockVelocityEddyDiffusivityDataReader()

# Numerical method for doing the actual integration
num_method_wrapper = MockTwoDNumMethod(config)

# Time points at which to compute particle positions
time = np.arange(time_start, time_end + time_step, time_step)

# Total number of time points at which to calcualte particle positions
n_times = len(time)

# The number of particles
n_particles = 1000

# Temporary arrays in which to store particle positions. All particles start at the origin.
x_positions = [0.0]*n_particles
y_positions = [0.0]*n_particles

# Arrays in which to save particle x and y positions
particle_x_positions = np.empty((n_times, n_particles), dtype=float)
particle_y_positions = np.empty((n_times, n_particles), dtype=float)

# Time integration
# ----------------

for t_idx, t in enumerate(time):
    # Save x and y positions for the last time point
    particle_x_positions[t_idx, :] = x_positions[:]
    particle_y_positions[t_idx, :] = y_positions[:]

    # Compute new x and y positions
    x_positions[:], y_positions[:] = num_method_wrapper.step(data_reader, t, x_positions, y_positions)

Visualising the result

We can plot changes in particles positions as a function of time, as simulated by PyLag, in the following way:

# Trim arrays to match analytic solution
particle_x_positions = particle_x_positions[1:]
particle_y_positions = particle_y_positions[1:]

# Set up the plot, as before
fig, ax = create_figure(figure_size=(16., 16.), font_size=font_size)
ax.set_xlim((x_grid_edges.min(), x_grid_edges.max()))
ax.set_ylim((y_grid_edges.min(), y_grid_edges.max()))
ax.set_xlabel('x', fontsize=font_size)
ax.set_ylabel('y', fontsize=font_size)
ax.set_aspect('equal')

# Animate particle trajectories
scatter = ax.scatter([], [], c='r', s=20, zorder=2)

def animate(i):
    data = np.array([particle_x_positions[i], particle_y_positions[i]])
    scatter.set_offsets(data.transpose())
    return (scatter,)

# Prevent the basic figure from being plotted too
plt.close(fig)

anim = animation.FuncAnimation(fig, animate, frames=n_times-1, interval=40, blit=True)

HTML(anim.to_html5_video())

It is difficult to discern a correspondance with the analytical solution when looking at the movemen of just a few particles. However, as the number of particles is increased, the global pattern becomes clearer. This can be tried by changing the number of particles above, rerunning the model and plotting the new result. Try this with, say, 10, 100, 1000 and 10000 particles.

To make a more direct comparison with the analytical solution, we can use kernel density estimators to compute the probability density field from individual particle positions. Here, we accomplish this in the following way using the scipy stats package:

from scipy import stats


# Define the grid
X, Y = np.meshgrid(x_grid, y_grid)
XY = np.vstack([X.ravel(), Y.ravel()])

# Compute the particle density field
density = np.empty((n_times-1, n_x_points, n_y_points), dtype=float)
for i in range(n_times-1):
    values = np.vstack([particle_x_positions[i, :], particle_y_positions[i, :]])
    kernel = stats.gaussian_kde(values)
    density[i, :, :] = np.reshape(kernel(XY).T, X.shape)

# Convert this into a concentration field
conc = density * data_reader.M

# Animate the result
# ------------------

font_size = 15
fig, ax = create_figure(figure_size=(16., 16.), font_size=font_size)
ax.set_xlim((x_grid_edges.min(), x_grid_edges.max()))
ax.set_ylim((y_grid_edges.min(), y_grid_edges.max()))
ax.set_xlabel('x', fontsize=font_size)
ax.set_ylabel('y', fontsize=font_size)
ax.set_aspect('equal')

# Set colour map and limits
cmap = colourmap('DYE')
vmin = 0.0
vmax = 10.0

# Set up the plot
pcolormesh= ax.pcolormesh(x_grid_edges, y_grid_edges, conc[0,:,:], vmin=vmin, vmax=vmax, cmap=cmap)
pos = ax.get_position().get_points()
cax = fig.add_axes([pos[1,0]+0.02, pos[0,1], 0.02, pos[1,1] - pos[0,1]])
cbar = fig.colorbar(pcolormesh, cax=cax)
cbar.ax.tick_params(labelsize=font_size)
cbar.set_label('C (kg m$\mathrm{^{-3}}$)', fontsize=font_size)

def animate(i):
    pcolormesh.set_array(conc[i,:,:].ravel())
    return (pcolormesh,)

# Prevent the basic figure from being plotted too
plt.close(fig)

anim = animation.FuncAnimation(fig, animate, frames=n_times-1, interval=40, blit=True)

HTML(anim.to_html5_video())

To take this comparison a step further, one can repeat the above steps for different numbers of particles. By quantifying differences between the two solutions at different times, one can also estimate the error in the particle field, and how this varies with the number of particles.