{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lateral advection and diffusion\n", "\n", "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.\n", "\n", "## Background theory\n", "\n", "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:\n", "\n", "$$\n", "\\begin{equation}\n", " \\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)\n", "\\end{equation}\n", "$$\n", "\n", "has the following analytical solution:\n", "\n", "$$\n", "\\begin{equation}\n", " 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}}\n", "\\end{equation}\n", "$$\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import warnings\n", "import numpy as np\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "from matplotlib import animation, rc\n", "from IPython.display import HTML\n", "\n", "from pylag.mock import MockVelocityEddyDiffusivityDataReader\n", "\n", "from pylag.processing.plot import create_figure, colourmap\n", "from pylag.processing.utils import get_grid_bands\n", "\n", "\n", "warnings.filterwarnings('ignore')\n", "\n", "# Ensure inline plotting\n", "%matplotlib inline\n", "\n", "# Set up the spatial grid\n", "xmin = -20.\n", "ymin = -20.\n", "xmax = 60.\n", "ymax = 60.\n", "n_x_points = 41\n", "n_y_points = 41\n", "x_grid = np.linspace(xmin, xmax, n_x_points, dtype=float)\n", "y_grid = np.linspace(ymin, ymax, n_y_points, dtype=float)\n", "\n", "# Grid edges for plotting\n", "x_grid_edges = get_grid_bands(x_grid)\n", "y_grid_edges = get_grid_bands(y_grid)\n", "\n", "# Set up the time grid\n", "t_min = 0.1\n", "t_max = 25.0\n", "time_step = 0.1\n", "t_grid = time = np.arange(t_min, t_max + time_step, time_step)\n", "n_times = len(t_grid)\n", "\n", "# Compute the analytic solution\n", "data_reader = MockVelocityEddyDiffusivityDataReader()\n", "C = np.empty((n_times, n_x_points, n_y_points), dtype=float)\n", "for t_idx, t in enumerate(t_grid):\n", " for x_idx, x in enumerate(x_grid):\n", " for y_idx, y in enumerate(y_grid):\n", " C[t_idx, x_idx, y_idx] = data_reader.get_concentration_analytic(t, x, y)\n", "\n", "\n", "# Animate the concentration field\n", "# -------------------------------\n", " \n", "font_size = 15\n", "fig, ax = create_figure(figure_size=(16., 16.), font_size=font_size)\n", "ax.set_xlim((x_grid_edges.min(), x_grid_edges.max()))\n", "ax.set_ylim((y_grid_edges.min(), y_grid_edges.max()))\n", "ax.set_xlabel('x', fontsize=font_size)\n", "ax.set_ylabel('y', fontsize=font_size)\n", "ax.set_aspect('equal')\n", "\n", "# Set colour map and limits\n", "cmap = colourmap('DYE')\n", "vmin = 0.0\n", "vmax = 10.0\n", "\n", "# Set up the plot\n", "pcolormesh= ax.pcolormesh(x_grid_edges, y_grid_edges, C[0,:,:], vmin=vmin, vmax=vmax, cmap=cmap)\n", "pos = ax.get_position().get_points()\n", "cax = fig.add_axes([pos[1,0]+0.02, pos[0,1], 0.02, pos[1,1] - pos[0,1]])\n", "cbar = fig.colorbar(pcolormesh, cax=cax)\n", "cbar.ax.tick_params(labelsize=font_size)\n", "cbar.set_label('C (kg m$\\mathrm{^{-3}}$)', fontsize=font_size)\n", "\n", "def animate(i):\n", " pcolormesh.set_array(C[i,:,:].ravel())\n", " return (pcolormesh,)\n", "\n", "# Prevent the basic figure from being plotted too\n", "plt.close(fig)\n", "\n", "anim = animation.FuncAnimation(fig, animate, frames=n_times, interval=40, blit=True)\n", "\n", "HTML(anim.to_html5_video())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In *PyLag*, particle trajectories are again calculated using a Random Displacement Model with the general form:\n", "\n", "$$\n", "\\begin{equation}\n", " 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)\n", "\\end{equation}\n", "$$\n", "\n", "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:\n", "\n", "$$\n", "\\begin{equation}\n", "\\label{eqn:dx}\n", " dX = u\\left(\\mathbf{x},t\\right) dt + \\left(2 A\\left(\\mathbf{x},t\\right)\\right)^{1/2}dW_X(t),\n", "\\end{equation}\n", "$$\n", "\n", "$$\n", "\\begin{equation}\n", "\\label{eqn:dy}\n", " dY = v\\left(\\mathbf{x},t\\right) dt + \\left(2 A\\left(\\mathbf{x},t\\right)\\right)^{1/2}dW_Y(t),\n", "\\end{equation}\n", "$$\n", "\n", "$$\n", "\\begin{equation}\n", "\\label{eqn:dz}\n", " dZ = 0,\n", "\\end{equation}\n", "$$\n", "\n", "*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." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running the particle tracking model\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import warnings\n", "import numpy as np\n", "from configparser import SafeConfigParser\n", "\n", "import pylag.random as random\n", "from pylag.numerics import StdNumMethod\n", "from pylag.mock import MockVelocityEddyDiffusivityDataReader, MockTwoDNumMethod\n", "\n", "warnings.filterwarnings('ignore')\n", "\n", "# Seed the random number generator\n", "random.seed(10)\n", "\n", "# Start/stop times\n", "time_start = 0.0 # Start time (s)\n", "time_end = 25.0 # End time (s)\n", "time_step = 0.1 # Time step (s)\n", "\n", "# PyLag reads many of its configuration parameters from a config which we create here\n", "config = SafeConfigParser()\n", "\n", "# Set the coordinate system, which here is a simple cartesian coordinate system\n", "config.add_section('SIMULATION')\n", "config.set('SIMULATION', 'coordinate_system', 'cartesian')\n", "\n", "# Parameters controlling the type of numerical integration scheme used\n", "config.add_section(\"NUMERICS\")\n", "config.set(\"NUMERICS\", \"num_method\", \"standard\")\n", "config.set(\"NUMERICS\", \"iterative_method\", \"AdvDiff_Milstein_3D\")\n", "config.set(\"NUMERICS\", \"time_step_diff\", str(time_step))\n", "\n", "# Parameters controlling particle vertical motion, which we set to default values\n", "config.set(\"SIMULATION\", \"depth_restoring\", \"False\")\n", "config.set(\"SIMULATION\", \"fixed_depth\", \"0.0\")\n", "\n", "# Specify the source of ocean current data which ensures data will be read\n", "config.add_section(\"OCEAN_DATA\")\n", "config.set(\"OCEAN_DATA\", \"name\", \"mock\")\n", "\n", "# The data reader\n", "data_reader = MockVelocityEddyDiffusivityDataReader()\n", "\n", "# Numerical method for doing the actual integration\n", "num_method_wrapper = MockTwoDNumMethod(config)\n", "\n", "# Time points at which to compute particle positions\n", "time = np.arange(time_start, time_end + time_step, time_step)\n", "\n", "# Total number of time points at which to calcualte particle positions\n", "n_times = len(time)\n", "\n", "# The number of particles\n", "n_particles = 1000\n", "\n", "# Temporary arrays in which to store particle positions. All particles start at the origin.\n", "x_positions = [0.0]*n_particles\n", "y_positions = [0.0]*n_particles\n", "\n", "# Arrays in which to save particle x and y positions\n", "particle_x_positions = np.empty((n_times, n_particles), dtype=float)\n", "particle_y_positions = np.empty((n_times, n_particles), dtype=float)\n", "\n", "# Time integration\n", "# ----------------\n", "\n", "for t_idx, t in enumerate(time):\n", " # Save x and y positions for the last time point\n", " particle_x_positions[t_idx, :] = x_positions[:]\n", " particle_y_positions[t_idx, :] = y_positions[:]\n", "\n", " # Compute new x and y positions\n", " x_positions[:], y_positions[:] = num_method_wrapper.step(data_reader, t, x_positions, y_positions)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Visualising the result\n", "\n", "We can plot changes in particles positions as a function of time, as simulated by *PyLag*, in the following way:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Trim arrays to match analytic solution\n", "particle_x_positions = particle_x_positions[1:]\n", "particle_y_positions = particle_y_positions[1:]\n", "\n", "# Set up the plot, as before\n", "fig, ax = create_figure(figure_size=(16., 16.), font_size=font_size)\n", "ax.set_xlim((x_grid_edges.min(), x_grid_edges.max()))\n", "ax.set_ylim((y_grid_edges.min(), y_grid_edges.max()))\n", "ax.set_xlabel('x', fontsize=font_size)\n", "ax.set_ylabel('y', fontsize=font_size)\n", "ax.set_aspect('equal')\n", "\n", "# Animate particle trajectories\n", "scatter = ax.scatter([], [], c='r', s=20, zorder=2)\n", "\n", "def animate(i):\n", " data = np.array([particle_x_positions[i], particle_y_positions[i]])\n", " scatter.set_offsets(data.transpose())\n", " return (scatter,)\n", "\n", "# Prevent the basic figure from being plotted too\n", "plt.close(fig)\n", "\n", "anim = animation.FuncAnimation(fig, animate, frames=n_times-1, interval=40, blit=True)\n", "\n", "HTML(anim.to_html5_video())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "\n", "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:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scipy import stats\n", "\n", "\n", "# Define the grid\n", "X, Y = np.meshgrid(x_grid, y_grid)\n", "XY = np.vstack([X.ravel(), Y.ravel()])\n", "\n", "# Compute the particle density field\n", "density = np.empty((n_times-1, n_x_points, n_y_points), dtype=float)\n", "for i in range(n_times-1):\n", " values = np.vstack([particle_x_positions[i, :], particle_y_positions[i, :]])\n", " kernel = stats.gaussian_kde(values)\n", " density[i, :, :] = np.reshape(kernel(XY).T, X.shape)\n", "\n", "# Convert this into a concentration field\n", "conc = density * data_reader.M\n", "\n", "# Animate the result\n", "# ------------------\n", " \n", "font_size = 15\n", "fig, ax = create_figure(figure_size=(16., 16.), font_size=font_size)\n", "ax.set_xlim((x_grid_edges.min(), x_grid_edges.max()))\n", "ax.set_ylim((y_grid_edges.min(), y_grid_edges.max()))\n", "ax.set_xlabel('x', fontsize=font_size)\n", "ax.set_ylabel('y', fontsize=font_size)\n", "ax.set_aspect('equal')\n", "\n", "# Set colour map and limits\n", "cmap = colourmap('DYE')\n", "vmin = 0.0\n", "vmax = 10.0\n", "\n", "# Set up the plot\n", "pcolormesh= ax.pcolormesh(x_grid_edges, y_grid_edges, conc[0,:,:], vmin=vmin, vmax=vmax, cmap=cmap)\n", "pos = ax.get_position().get_points()\n", "cax = fig.add_axes([pos[1,0]+0.02, pos[0,1], 0.02, pos[1,1] - pos[0,1]])\n", "cbar = fig.colorbar(pcolormesh, cax=cax)\n", "cbar.ax.tick_params(labelsize=font_size)\n", "cbar.set_label('C (kg m$\\mathrm{^{-3}}$)', fontsize=font_size)\n", "\n", "def animate(i):\n", " pcolormesh.set_array(conc[i,:,:].ravel())\n", " return (pcolormesh,)\n", "\n", "# Prevent the basic figure from being plotted too\n", "plt.close(fig)\n", "\n", "anim = animation.FuncAnimation(fig, animate, frames=n_times-1, interval=40, blit=True)\n", "\n", "HTML(anim.to_html5_video())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 2 }