{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Devito CFD Tutorial series\n",
"\n",
"The following series of notebook tutorials will demonstrate the use of Devito and it's SymPy-based API to solve a set of classic examples from Computational Fluid Dynamics (CFD). The tutorials are based on the excellent tutorial series _CFD Python: 12 steps to Navier-Stokes_ by Lorena Barba and focus on the implementation with Devito rather than pure CFD or finite difference theory. For a refresher on how to implement 2D finite difference solvers for CFD problems, please see the original tutorial series here: \n",
"http://lorenabarba.com/blog/cfd-python-12-steps-to-navier-stokes/"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Example 1: Linear convection in 2D\n",
"\n",
"Lets start with a simple 2D convection example - step 5 in the original blog. This will already allow us to demonstrate a lot about the use of Devito's symbolic data objects and how to use them to build a simple operator directly from the symbolic notation of the equation. The governing equation we will implement in this tutorial is:\n",
"\n",
"$$\\frac{\\partial u}{\\partial t}+c\\frac{\\partial u}{\\partial x} + c\\frac{\\partial u}{\\partial y} = 0$$\n",
"\n",
"In order to implement this equation we first discretize it using forward differences in time and backward differences in space. Just as the original tutorial, we will use $u_{i,j}^n$ to denote a finite difference stencil point with $i$ and $j$ denoting spatial indices and $n$ denoting the time index. So, after re-arranging the discretized equation for the forward stencil point in time we get \n",
"\n",
"$$u_{i,j}^{n+1} = u_{i,j}^n-c \\frac{\\Delta t}{\\Delta x}(u_{i,j}^n-u_{i-1,j}^n)-c \\frac{\\Delta t}{\\Delta y}(u_{i,j}^n-u_{i,j-1}^n)$$\n",
"\n",
"Using this, we can start deriving the computational stencil for this equation. Let's first look at the original _pure Python_ implementation of the linear convection flow - but first we import our tools and define some parameters:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"scrolled": true
},
"outputs": [
{
"output_type": "stream",
"name": "stdout",
"text": [
"dx 0.025, dy 0.025\n"
]
}
],
"source": [
"from examples.cfd import plot_field, init_hat\n",
"import numpy as np\n",
"%matplotlib inline\n",
"\n",
"# Some variable declarations\n",
"nx = 81\n",
"ny = 81\n",
"nt = 100\n",
"c = 1.\n",
"dx = 2. / (nx - 1)\n",
"dy = 2. / (ny - 1)\n",
"print(\"dx %s, dy %s\" % (dx, dy))\n",
"sigma = .2\n",
"dt = sigma * dx"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**A small note on style:** Throughout this tutorial series we will use utility functions to plot the various 2D functions and data sets we deal with. These are all taken from the original tutorial series, but have been slightly modified for our purposes. One of the differences readers might find is that the original series uses _(y, x)_ indexing for 2d data arrays, whereas many of the examples have been adapted to use _(x, y)_ notation in our tutorials.\n",
"\n",
"So, let's start by creating a simple 2D function and initialising it with a \"hat function\". We will use that initialisation function a lot, so it comes from our utility scripts:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"output_type": "display_data",
"data": {
"text/plain": "",
"image/svg+xml": "\n\n