{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Lateral advection\n",
"\n",
"PyLag can be easily used with analytical models of different flow fields. This provides a convenient means of testing\n",
"the model implementation. The example described here is taken from Kreyszig, E. (2006) Advanced Engineering\n",
"Mathematics, Ch. 18. P762.\n",
"\n",
"## Background theory\n",
"\n",
"In the example, the complex potential with $F(z) = z^{2} = x^{2} - y^{2} + 2ixy$ models a flow with equipotential lines and streamlines given by:\n",
"\n",
"$$\n",
"\\begin{eqnarray}\n",
" \\Phi &=& x^2 - y^2 &=& const \\\\\n",
" \\Psi &=& 2xy &=& const\n",
"\\end{eqnarray}\n",
"$$ \n",
"\n",
"The velocity vector is:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" \\mathbf{u} = 2(x - iy)\n",
"\\end{equation}\n",
"$$\n",
"\n",
"with components:\n",
"\n",
"$$\n",
"\\begin{eqnarray}\n",
" u &=& &2&x \\\\\n",
" v &=& -&2&y\n",
"\\end{eqnarray}\n",
"$$\n",
"\n",
"The speed of the flow is given by:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" |\\mathbf{u}| = \\sqrt{x^2 + y^2}\n",
"\\end{equation}\n",
"$$\n",
"\n",
"The flow may be interpreted as the flow in a channel bounded by the positive coordinate axes and a hyperbola, given by $xy = const$. To calculate particle trajectiories, PyLag solves the equation:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\label{eqn:partmotion}\n",
" \\dfrac{\\partial}{\\partial t} \\mathbf{X}_{i}(t, \\mathbf{r}_{i}) = \\mathbf{U}_{i}(t, \\mathbf{X}_{i})\n",
"\\end{equation}\n",
"$$\n",
"\n",
"where $\\mathbf{r}_{i} = \\mathbf{X}_i(t=t_0)$ is the position vector of particle $i$ at time $t=t_{0}$, and $\\mathbf{U}_{i}$ is the particle's velocity vector. Here, $\\mathbf{U}_{i} = [\\mathbf{u}(t, \\mathbf{x})]_{\\mathbf{x} = \\mathbf{X}_{i}}$ where $\\mathbf{u}$ is the fluid velocity vector, as given above.\n",
"\n",
"In practise, *PyLag* computes solutions to the above equation using a Random Displacement Model of the 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",
"Here, $dX_{j} = d\\mathbf{X}$ is the incremental change in a particle's position and $D_{jk}$ is the diffusion coefficient or tensor. The first term on the right-hand side (RHS) of the equation is a deterministic drift term; the second is a stochastic term that models the action of diffusion. $dW_{j}(t)$ is an incremental Wiener process that builds stochasticity into the model.\n",
"\n",
"In this particular problem, the diffusion coefficient is zero. In component form, the model then reduces to:\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\label{eqn:dx}\n",
" dX = u\\left(\\mathbf{x},t\\right) dt,\n",
"\\end{equation}\n",
"$$\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\label{eqn:dy}\n",
" dY = v\\left(\\mathbf{x},t\\right) dt,\n",
"\\end{equation}\n",
"$$\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\label{eqn:dz}\n",
" dZ = 0,\n",
"\\end{equation}\n",
"$$\n",
"\n",
"which can be solved using one of the standard numerical integration schemes supported by *PyLAg*."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Running the particle tracking model\n",
"\n",
"To compute particle trajectories within the given flow field using *PyLag*, we first subclass *PyLag's* `DataReader`. The new subclass must be implemented in *Cython*, since several of its functions have arguments without a direct counterpart in Python. The above model has been implemented in `pylag.mock`, and is called `MockVelocityDataReader`.\n",
"\n",
"
\n",
"\n",
"**Note:**\n",
"\n",
"It is possible to extend *PyLag* by implementing alternative concrete `DataReader` subclasses in `pylag.mock`. Alternatively, you may choose to implement these in a separate module. Either way, each module must first by cythonized and compiled before it can be used.\n",
"\n",
"
\n",
"\n",
"We will start three particles off with initial coordinates $X_1(t = t_0) = (0.125, 1.)$, $X_2(t = t_0) = (0.3525, 1.)$ and $X_3(t = t_0) = (0.5, 1.)$. The corresponding streamlines are given by $\\Psi = \\{0.5, 1, 2\\}$. Given the initial particle coordinates, we can use PyLag to compute particle pathlines, which can be compared directly with the above streamlines. To do this, we use one of PyLag's numerical integration schemes. The integration is managed by an object of type `pylag.numerics.NumIntegrator`, which in turn farms out most of the hard work to an object of type `pylag.numerics.ItMethod`. Given the type of flow field, we will use `StdNumIntegrator` to manage the integration, and the `Adv_RK4_2D` iterative method, which encodes a fourth order runga kutta scheme. The following code uses *PyLag* to compute particle positions as a function of time:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import warnings\n",
"import numpy as np\n",
"import matplotlib\n",
"from configparser import SafeConfigParser\n",
"\n",
"from pylag.numerics import StdNumMethod\n",
"from pylag.particle_cpp_wrapper import ParticleSmartPtr\n",
"from pylag.mock import MockVelocityDataReader\n",
"\n",
"# Ensure inline plotting\n",
"%matplotlib inline\n",
"\n",
"warnings.filterwarnings('ignore')\n",
"\n",
"# Set the time step to 0.1s\n",
"time_step = 0.01\n",
"\n",
"# PyLag reads many of its parameters from a run config, which we create here\n",
"config = SafeConfigParser()\n",
"config.add_section('SIMULATION')\n",
"config.add_section('NUMERICS')\n",
"\n",
"# Set the coordinate system, which here is a simple cartesian coordinate system\n",
"config.set('SIMULATION', 'coordinate_system', 'cartesian')\n",
"\n",
"# Set the type of iterative method to a 2D fourth order Runga Kutta method\n",
"config.set('NUMERICS', 'iterative_method', 'Adv_RK4_2D')\n",
"\n",
"# Set the time step\n",
"config.set('NUMERICS', 'time_step_adv', str(time_step))\n",
"\n",
"# The data reader\n",
"data_reader = MockVelocityDataReader()\n",
"\n",
"# The numerical integration scheme - this uses the config to initialise its state\n",
"num_method = StdNumMethod(config)\n",
"\n",
"# The three particles, saved in a dictionary. \n",
"particles = {'Particle1': ParticleSmartPtr(x1=0.125, x2=1.0),\n",
" 'Particle2': ParticleSmartPtr(x1=0.25, x2=1.0),\n",
" 'Particle3': ParticleSmartPtr(x1=0.375, x2=1.0)}\n",
"\n",
"# Dictionary in which particle position data will be saved. x and y positions are stored\n",
"# in lists. Initialise this with the particle's position at t = t_0.\n",
"positions_out = {}\n",
"for key in particles.keys():\n",
" positions_out[key] = {'x1': [particles[key].x1], 'x2': [particles[key].x2]}\n",
"\n",
"# We will start the model off at t = 0 s and integrate it forward for 2 s\n",
"time_grid = np.arange(0.0, 2.0, time_step).tolist()\n",
"\n",
"# Perform the integration\n",
"for time in time_grid:\n",
" for key in particles.keys():\n",
" # Udate partcle positions.\n",
" num_method.step_wrapper(data_reader, time, particles[key])\n",
"\n",
" # Save particle positions\n",
" positions_out[key]['x1'].append(particles[key].x1)\n",
" positions_out[key]['x2'].append(particles[key].x2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Visualising the result\n",
"\n",
"Below we animate the particle positions on top of the streamlines. A quiver plot of the velocity field is also shown."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"%matplotlib inline\n",
"\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from matplotlib import animation, rc\n",
"from IPython.display import HTML\n",
"\n",
"font_size=15\n",
"\n",
"fig, ax = plt.subplots(figsize = (8, 8))\n",
"\n",
"# Plot the line x = y\n",
"ax.plot([0,1], [0,1], 'k--')\n",
"\n",
"# Add a quiver plot of the velocity field\n",
"x_grid, y_grid = np.meshgrid(np.linspace(0,1,20),np.linspace(0,1,20)) \n",
"quiver_plot = ax.quiver(x_grid, y_grid, 2.*x_grid, -2.*y_grid)\n",
"\n",
"# Add streamlines\n",
"y = np.linspace(0.01, 5., 100, dtype=float)\n",
"for psi in [0.25, 0.5, 0.75]:\n",
" x = 0.5 * psi / y\n",
" \n",
" ax.plot(x, y, c='b', zorder=1)\n",
"\n",
"# Configure the plot\n",
"ax.set_xlim((0, 1))\n",
"ax.set_ylim((0, 1))\n",
"\n",
"ax.set_xlabel('x', fontsize=font_size)\n",
"ax.set_ylabel('y', fontsize=font_size)\n",
"\n",
"ax.tick_params(axis='both', which='major', labelsize=font_size)\n",
"ax.set_aspect('equal')\n",
"\n",
"# Animate particle trajectories\n",
"scatter = ax.scatter([], [], c='r', s=100, zorder=2)\n",
"\n",
"def animate(i):\n",
" xpos = []\n",
" ypos = []\n",
" for key in particles.keys():\n",
" xpos.append(positions_out[key]['x1'][i])\n",
" ypos.append(positions_out[key]['x2'][i])\n",
"\n",
" data = np.array([xpos, ypos])\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=110, interval=20, blit=True)\n",
"\n",
"HTML(anim.to_html5_video())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is clear that all three particles closely follow their respective streamlines, as they should. The dashed line marks the point along each streamline where the speed is at a minimum. This can be seen in the particle trajectories, with particles decelerating as they approach the point where the dashed line intersects each streamline, then accelerating away from it."
]
}
],
"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
}