{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Finite Difference Operators in Many Dimensions using Kronecker Products\n",
"**Prof. Tony Saad (www.tsaad.net)
Department of Chemical Engineering
University of Utah**\n",
"
"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is very easy to build complex 2D and 3D finite difference operators using Kronecker products. We will start by defining a basic forward difference operator for the first derivative. We can choose other operators as well, but this particular case produces tight stencils for Laplacians.\n",
"\n",
"The forward difference is given by\n",
"\\begin{equation}\n",
"\\frac{\\text{d}u}{\\text{d}x} = \\frac{u_{i+1} - u_{i}}{\\Delta x}\n",
"\\end{equation}\n",
"\n",
"In operator form, and for a periodic case on a 1D grid, this operator looks like\n",
"\\begin{equation}\n",
"\\mathbf{D}^+ =\n",
"\\left[ {\\begin{array}{*{20}{c}}\n",
"{ - 1}&1&{}&{}&{}\\\\\n",
"{}&{ - 1}&1&{}&{}\\\\\n",
"{}&{}&{ - 1}&1&{}\\\\\n",
"{}&{}&{}&{ - 1}&1\\\\\n",
"1&{}&{}&{}&{ - 1}\n",
"\\end{array}} \\right]\n",
"\\end{equation}\n",
"\n",
"For periodic boundaries, this operator can be defined in Python as\n",
"```Python\n",
"def d_plus(n, dx=1):\n",
" # periodic forward difference\n",
" d = np.ones(n, dtype=int)\n",
" ud = np.ones(n-1, dtype=int)\n",
" D = diag(-d) + diag(ud,1)\n",
" D[-1,0] = 1\n",
" return D/dx\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We will also need a backward difference operator, $\\mathbf{D}^-$, to enable us to compute second derivatives and Laplacians. It can be shown that $\\mathbf{D}^-$ is given by the negative of the transpose of $\\mathbf{D}^+$\n",
"\\begin{equation}\n",
"\\mathbf{D}^- = - \\mathbf{D}^{+T}\n",
"\\end{equation}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Two-Dimensions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To construct forward and backward difference operators for an entire 2D grid, we can using the Kronecker product. For a uniform structured grid of size $N_x \\times N_y$ we can construct four operators: $\\mathbf{D}_x^+$, $\\mathbf{D}_x^-$, $\\mathbf{D}_y^+$, $\\mathbf{D}_y^-$ where the subscript refers to the direction in which the derivative is taken, e.g. $\\mathbf{D}_y^+ u= \\frac{\\partial u}{\\partial y}$ represents the forward difference approximation of $\\frac{\\partial u}{\\partial y}$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Kronecker product, defined via the symbol $\\otimes$, is simply a way to produce matrices based on matrix outer products. Given matrices $\\mathbf{A} = a_{ij}$ and $\\mathbf{B} = b_{ij}$, then\n",
"\\begin{equation}\n",
"\\mathbf{A}\\otimes\\mathbf{B} = \\begin{bmatrix} a_{11} \\mathbf{B} & \\cdots & a_{1n}\\mathbf{B} \\\\ \\vdots & \\ddots & \\vdots \\\\ a_{m1} \\mathbf{B} & \\cdots & a_{mn} \\mathbf{B} \\end{bmatrix}\n",
"\\end{equation}\n",
"\n",
"An example from Wikipedia is\n",
"\n",
"\\begin{equation}\n",
" \\begin{bmatrix}\n",
" 1 & 2 \\\\\n",
" 3 & 4 \\\\\n",
" \\end{bmatrix}\n",
" \\otimes\n",
" \\begin{bmatrix}\n",
" 0 & 5 \\\\\n",
" 6 & 7 \\\\\n",
" \\end{bmatrix}=\n",
" \\begin{bmatrix}\n",
" 1 \\cdot \\begin{bmatrix}\n",
" 0 & 5 \\\\\n",
" 6 & 7 \\\\\n",
" \\end{bmatrix} & \n",
" 2 \\cdot \\begin{bmatrix}\n",
" 0 & 5 \\\\\n",
" 6 & 7 \\\\\n",
" \\end{bmatrix} \\\\\n",
" 3 \\cdot \\begin{bmatrix}\n",
" 0 & 5 \\\\\n",
" 6 & 7 \\\\\n",
" \\end{bmatrix} & \n",
" 4 \\cdot \\begin{bmatrix}\n",
" 0 & 5 \\\\\n",
" 6 & 7 \\\\\n",
" \\end{bmatrix} \\\\\n",
" \\end{bmatrix}=\n",
" \\begin{bmatrix}\n",
" 1\\cdot 0 & 1\\cdot 5 & 2\\cdot 0 & 2\\cdot 5 \\\\\n",
" 1\\cdot 6 & 1\\cdot 7 & 2\\cdot 6 & 2\\cdot 7 \\\\\n",
" 3\\cdot 0 & 3\\cdot 5 & 4\\cdot 0 & 4\\cdot 5 \\\\\n",
" 3\\cdot 6 & 3\\cdot 7 & 4\\cdot 6 & 4\\cdot 7 \\\\\n",
" \\end{bmatrix}.\n",
" \\end{equation}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Back to our finite difference operators, on a 2D grid, we have the following\n",
"\\begin{equation}\n",
"\\bf{D}_x^ + = \\left[ \\begin{array}{*{20}{c}}\n",
"\\bf{D}^+ &{}&{}&{}&{}\\\\\n",
"{}&\\bf{D}^+&{}&{}&{}\\\\\n",
"{}&{}& \\ddots &{}&{}\\\\\n",
"{}&{}&{}& \\ddots &{}\\\\\n",
"{}&{}&{}&{}&\\bf{D}^+\n",
"\\end{array} \\right] = \\bf{I}_y \\otimes \\bf{D}^+({N_y})\n",
"\\end{equation}\n",
"where $\\mathbf{I}_y$ is the identity matrix of size $N_y\\times N_y$ and $\\mathbf{D}^+$ is of size $N_x\\times N_x$. Similarly, we have\n",
"\\begin{equation}\n",
"\\mathbf{D}_x^- = \\mathbf{I}_y\\otimes \\mathbf{D}^-(N_x)\n",
"\\end{equation}\n",
"The effect of the Kronecker product $\\mathbf{I}_y\\otimes \\mathbf{D}$ is to replicate $\\mathbf{D}$ by $N_y$ times."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The situation is a bit more different for the $y$-derivatives. In this case, we want to offset the off-diagonal components of $\\mathbf{D}$ by $N_x$, i.e. $(i+1, j) \\to (i, j+1)$. To do this, we have to right-multiply $\\mathbf{D}$ by $\\mathbf{I}_x$.\n",
"\\begin{equation}\n",
"\\mathbf{D}_y^+ = \\mathbf{D}^+(N_y) \\otimes \\mathbf{I}_x\n",
"\\end{equation}\n",
"and\n",
"\\begin{equation}\n",
"\\mathbf{D}_y^- = \\mathbf{D}^-(N_y) \\otimes \\mathbf{I}_x\n",
"\\end{equation}\n",
"where $\\mathbf{I}_x$ is the identity matrix of size $N_x\\times N_x$ and $\\mathbf{D}^{+,-}$ is of size $N_y\\times N_y$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Three Dimensions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In three dimensions, the operators can be constructed as follows\n",
"\\begin{equation}\n",
"\\mathbf{D_x^{+,-}} = \\mathbf{I}_z \\otimes \\mathbf{I}_y \\otimes \\mathbf{D}^{+,-} \\\\\n",
"\\mathbf{D_y^{+,-}} = \\mathbf{I}_z \\otimes \\mathbf{D}^{+,-} \\otimes \\mathbf{I}_x \\\\\n",
"\\mathbf{D_z^{+,-}} = \\mathbf{D}^{+,-} \\otimes \\mathbf{I}_y \\otimes \\mathbf{I}_x \\\\ \n",
"\\end{equation}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Bottom Line**\n",
"1. $\\mathbf{I}(N) \\otimes \\mathbf{B}(M)$ creates a block diagonal matrix of $\\mathbf{B}$'s\n",
"2. $\\mathbf{B}(M) \\otimes \\mathbf{I}(N)$ offsets the off-diagonal terms in $\\mathbf{B}$ by $N$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Implementation\n",
"We will implement the Kronecker product formulation in a single class called `GridOps`"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from numpy import zeros, ones, eye, diag, kron\n",
"\n",
"class GridOps:\n",
" '''\n",
" Creates x, y, and z-direction 2D, or 3D versions of a 1D operator defined by Op1D.\n",
" N : A list or array containing the grid dimensions, [Nx, Ny, Nz].\n",
" Op1D: The name of a basic 1D operator. Op1D must be of the form Op1D(n,h) where n is a single valued integer\n",
" designating the number of points in that one-dimension and h is the spacing.\n",
" Δ : A list or array containing the grid spacing, [dx, dy, dz]. Defaults to [1,1,1]\n",
" '''\n",
" def __init__(self, Op1D, N, Δ=[1,1,1]):\n",
" self.N = N\n",
" self.Δ = Δ\n",
" self.Op1D = Op1D\n",
" \n",
" def OpX(self):\n",
" '''\n",
" Creates x-direction 2D, or 3D versions of a 1D operator defined by Op1D.\n",
" ''' \n",
" N = self.N\n",
" Δ = self.Δ\n",
" Op = self.Op1D\n",
" return kron(eye(N[2]), kron(eye(N[1]), Op (N[0],Δ[0]) ) )\n",
" \n",
" def OpY(self):\n",
" '''\n",
" Creates y-direction 2D, or 3D versions of a 1D operator defined by Op1D.\n",
" '''\n",
" N = self.N\n",
" Δ = self.Δ\n",
" Op = self.Op1D\n",
" return kron( eye(N[2]), kron(Op(N[1],Δ[1]), eye(N[0]) ) ) \n",
"\n",
" def OpZ(self):\n",
" '''\n",
" Creates z-direction 2D, or 3D versions of a 1D operator defined by Op1D.\n",
" ''' \n",
" N = self.N\n",
" Δ = self.Δ\n",
" Op = self.Op1D \n",
" return kron(Op(N[2],Δ[2]), kron(eye(N[1]), eye(N[0]) ) ) \n",
" \n",
" def Sum(self):\n",
" '''\n",
" Returns the sum of all three operators. Useful for returning a Laplacian for example.\n",
" '''\n",
" return self.OpX() + self.OpY() + self.OpZ()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Below are some example 1D finite difference operators"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {},
"outputs": [],
"source": [
"def d_plus1d(n, h=1):\n",
" '''\n",
" Defines a basic forward difference operator in 1D (first derivative, first order)\n",
" '''\n",
" # periodic forward difference\n",
" d = np.ones(n, dtype=int)\n",
" ud = np.ones(n-1, dtype=int)\n",
" D = diag(-d) + diag(ud,1)\n",
" D[-1,0] = 1\n",
" return D/h\n",
"\n",
"def d_minus1d(n, h=1):\n",
" '''\n",
" Defines a basic backward difference operator in 1D (first derivative, first order)\n",
" '''\n",
" return -d_plus1d(n,h).T\n",
"\n",
"def laplacian1d(n, h=1):\n",
" '''\n",
" Defines a basic 1D Laplacian operator (second derivative), cell centered, second order\n",
" ''' \n",
" result = 0 if n <= 1 else d_plus1d(n,h) @ d_minus1d(n,h)\n",
" return result"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from numpy import zeros, ones, eye, diag, kron\n",
"%matplotlib inline\n",
"import matplotlib.pyplot as plt\n",
"%config InlineBackend.figure_format = 'svg'"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[-1. 1. 0. 0. 0.]\n",
" [ 0. -1. 1. 0. 0.]\n",
" [ 0. 0. -1. 1. 0.]\n",
" [ 0. 0. 0. -1. 1.]\n",
" [ 1. 0. 0. 0. -1.]]\n"
]
}
],
"source": [
"print(d_plus1d(5))"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 1. -0. -0. -0. -1.]\n",
" [-1. 1. -0. -0. -0.]\n",
" [-0. -1. 1. -0. -0.]\n",
" [-0. -0. -1. 1. -0.]\n",
" [-0. -0. -0. -1. 1.]]\n"
]
}
],
"source": [
"print(d_minus1d(5))"
]
},
{
"cell_type": "code",
"execution_count": 53,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[-2. 1. 0. 0. 1.]\n",
" [ 1. -2. 1. 0. 0.]\n",
" [ 0. 1. -2. 1. 0.]\n",
" [ 0. 0. 1. -2. 1.]\n",
" [ 1. 0. 0. 1. -2.]]\n"
]
}
],
"source": [
"print(laplacian1d(5))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Two Dimensional Example\n",
"In what follows, we build a few 2D versions of the operators we defined above"
]
},
{
"cell_type": "code",
"execution_count": 54,
"metadata": {},
"outputs": [],
"source": [
"N = [3,3,1] # set grid size to 4x4"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To define 2D versions of $\\mathbf{D}^{+,-}$ and the 1d Laplacian, we simply call the class `GridOps`"
]
},
{
"cell_type": "code",
"execution_count": 55,
"metadata": {},
"outputs": [],
"source": [
"Dplus2D = GridOps(d_plus1d, N)\n",
"Dminus2D = GridOps(d_minus1d, N)\n",
"Lap2D = GridOps(laplacian1d, N)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now visualize these operators by accessing the methods defined by `GridOps`. For example, The forward difference version of $\\partial/\\partial x$ on this 2D Grid is given calling `OpX()` on the corresponding `GridOps`"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[-1. 1. 0. -0. 0. 0. -0. 0. 0.]\n",
" [ 0. -1. 1. 0. -0. 0. 0. -0. 0.]\n",
" [ 1. 0. -1. 0. 0. -0. 0. 0. -0.]\n",
" [-0. 0. 0. -1. 1. 0. -0. 0. 0.]\n",
" [ 0. -0. 0. 0. -1. 1. 0. -0. 0.]\n",
" [ 0. 0. -0. 1. 0. -1. 0. 0. -0.]\n",
" [-0. 0. 0. -0. 0. 0. -1. 1. 0.]\n",
" [ 0. -0. 0. 0. -0. 0. 0. -1. 1.]\n",
" [ 0. 0. -0. 0. 0. -0. 1. 0. -1.]]\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Dxp = Dplus2D.OpX()\n",
"print(Dxp)\n",
"plt.matshow(Dxp)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similary, the backward difference representation of $\\partial / \\partial y$ is"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[ 1. 0. 0. -0. -0. -0. -1. -0. -0.]\n",
" [ 0. 1. 0. -0. -0. -0. -0. -1. -0.]\n",
" [ 0. 0. 1. -0. -0. -0. -0. -0. -1.]\n",
" [-1. -0. -0. 1. 0. 0. -0. -0. -0.]\n",
" [-0. -1. -0. 0. 1. 0. -0. -0. -0.]\n",
" [-0. -0. -1. 0. 0. 1. -0. -0. -0.]\n",
" [-0. -0. -0. -1. -0. -0. 1. 0. 0.]\n",
" [-0. -0. -0. -0. -1. -0. 0. 1. 0.]\n",
" [-0. -0. -0. -0. -0. -1. 0. 0. 1.]]\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Dym = Dminus2D.OpY()\n",
"print(Dym)\n",
"plt.matshow(Dym)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Laplacian can be composed by adding `OpX()` and `OpY()'"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[[-4. 1. 1. 1. 0. 0. 1. 0. 0.]\n",
" [ 1. -4. 1. 0. 1. 0. 0. 1. 0.]\n",
" [ 1. 1. -4. 0. 0. 1. 0. 0. 1.]\n",
" [ 1. 0. 0. -4. 1. 1. 1. 0. 0.]\n",
" [ 0. 1. 0. 1. -4. 1. 0. 1. 0.]\n",
" [ 0. 0. 1. 1. 1. -4. 0. 0. 1.]\n",
" [ 1. 0. 0. 1. 0. 0. -4. 1. 1.]\n",
" [ 0. 1. 0. 0. 1. 0. 1. -4. 1.]\n",
" [ 0. 0. 1. 0. 0. 1. 1. 1. -4.]]\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"DG = Lap2D.OpX() + Lap2D.OpY()\n",
"print(DG)\n",
"plt.matshow(DG)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example: Helmholtz Hodge Decomposition"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We use the operators to conduct a projection\n",
"\\begin{equation}\n",
"\\mathbf{u} = \\mathbf{u}_\\perp + \\nabla \\phi\n",
"\\end{equation}\n",
"such that $\\nabla\\cdot\\mathbf{u}_\\perp = 0$. Then, we take the divergence\n",
"\\begin{equation}\n",
"\\nabla \\cdot \\mathbf{u} = \\nabla^2 \\phi\n",
"\\end{equation}\n",
"In terms of operators, we have\n",
"\\begin{equation}\n",
"D \\mathbf{u} = D G\\phi\n",
"\\end{equation}\n",
"then\n",
"\\begin{equation}\n",
"\\phi = (D G)^{-1} D \\mathbf{u}\n",
"\\end{equation}\n",
"and finally\n",
"\\begin{equation}\n",
"\\mathbf{u}_\\perp = \\mathbf{u} - G \\phi = \\mathbf{u} - G (D G)^{-1} D \\mathbf{u}\n",
"\\end{equation}\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"nx = 32\n",
"ny = 16\n",
"N = [nx, ny, 1]\n",
"\n",
"X0 = 0.0\n",
"Xl = 1.0\n",
"Y0 = 0.0\n",
"Yl = 0.5\n",
"x_ = np.linspace(0,1,nx)\n",
"y_ = np.linspace(0,1,ny)\n",
"\n",
"dx = (Xl-X0)/(nx-1) if nx!=1 else 1\n",
"dy = (Yl-Y0)/(ny-1) if ny!=1 else 1\n",
"Δ = [dx, dy, 1.0]\n",
"\n",
"X,Y = np.meshgrid(x_,y_)\n",
"\n",
"a = 1.0\n",
"b = 0.7\n",
"c = 0.2\n",
"\n",
"ux = a*np.exp((-(X-b+0.25)**2)/(2*c**2) -((Y-b+0.25)**2)/(2*c**2))\n",
"uy = a*np.exp((-(X-b-0.1)**2)/(2*c**2) -((Y-b+0.1)**2)/(2*c**2))\n",
"\n",
"plt.quiver(ux,uy)\n",
"# plt.contourf(X,Y,ux,10,cmap=\"RdBu_r\")"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"ux = np.sin(X*Y)\n",
"uy = np.cos(4.0*X*Y)\n",
"plt.quiver(ux,uy)"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"168393.32145065966\n"
]
}
],
"source": [
"Dplus2D = GridOps(d_plus1d, N, Δ)\n",
"Dminus2D = GridOps(d_minus1d, N, Δ)\n",
"Lap2D = GridOps(laplacian1d, N, Δ)\n",
"\n",
"div_x = Dplus2D.OpX()\n",
"div_y = Dplus2D.OpY()\n",
"\n",
"grad_x = Dminus2D.OpX()\n",
"grad_y = Dminus2D.OpY()\n",
"\n",
"DG = div_x@grad_x + div_y@grad_y\n",
"DG[0]= 0\n",
"DG[0,0] = 1.0\n",
"print(np.linalg.cond(DG))\n",
"DG_inv = np.linalg.inv(DG)\n",
"Du = div_x@ux.ravel() + div_y@uy.ravel()\n",
"phi = DG_inv@Du"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"uxperp = ux.ravel() - grad_x@phi\n",
"uyperp = uy.ravel() - grad_y@phi\n",
"plt.quiver(X,Y,uxperp.reshape(nx,ny), uyperp.reshape(nx,ny),pivot='middle')"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"execution_count": 63,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"div_uperp = div_x@uxperp.ravel() + div_y@uyperp.ravel()\n",
"plt.matshow(div_uperp.reshape(ny,nx))\n",
"plt.colorbar()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Three Dimensions"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"Text(0.5, 1.05, '$\\nabla^2$, central diff \\n')"
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"N=[4,4,4]\n",
"\n",
"Dplus3D = GridOps(d_plus1d, N)\n",
"Dminus3D = GridOps(d_minus1d, N)\n",
"Lap3D = GridOps(laplacian1d, N)\n",
"\n",
"plt.matshow(Dplus3D.OpX())\n",
"plt.title('$\\partial / \\partial x$, forward diff \\n')\n",
"\n",
"plt.matshow(Dplus3D.OpY())\n",
"plt.title('$\\partial / \\partial y$, forward diff \\n')\n",
"\n",
"plt.matshow(Dplus3D.OpZ())\n",
"plt.title('$\\partial / \\partial z$, forward diff \\n')\n",
"\n",
"plt.matshow(Dminus3D.OpX())\n",
"plt.title('$\\partial / \\partial x$, backward diff \\n')\n",
"\n",
"plt.matshow(Dminus3D.OpY())\n",
"plt.title('$\\partial / \\partial y$, backward diff \\n')\n",
"\n",
"plt.matshow(Dminus3D.OpZ())\n",
"plt.title('$\\partial / \\partial z$, backward diff \\n')\n",
"\n",
"DG = Lap3D.Sum()\n",
"plt.matshow(DG)\n",
"plt.title('$\\nabla^2$, central diff \\n')"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"CSS style adapted from https://github.com/barbagroup/CFDPython. Copyright (c) Barba group\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import urllib\n",
"import requests\n",
"from IPython.core.display import HTML\n",
"def css_styling():\n",
" styles = requests.get(\"https://raw.githubusercontent.com/saadtony/NumericalMethods/master/styles/custom.css\")\n",
" return HTML(styles.text)\n",
"css_styling()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"anaconda-cloud": {},
"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.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}