Text provided under a Creative Commons Attribution license, CC-BY. All code is made available under the FSF-approved BSD-3 license. (c) Lorena A. Barba, Gilbert F. Forsyth 2017. Thanks to NSF for support via CAREER award #1149784. @LorenaABarba
Up to now, all of our work has been in one spatial dimension (Steps 1 to 4). We can learn a lot in just 1D, but let's grow up to flatland: two dimensions.
In the following exercises, you will extend the first four steps to 2D. To extend the 1D finite-difference formulas to partial derivatives in 2D or 3D, just apply the definition: a partial derivative with respect to
In 2D space, a rectangular (uniform) grid is defined by the points with coordinates:
Now, define
Hence, for a first-order partial derivative in the
and similarly in the
The PDE governing 2-D Linear Convection is written as
This is the exact same form as with 1-D Linear Convection, except that we now have two spatial dimensions to account for as we step forward in time.
Again, the timestep will be discretized as a forward difference and both spatial steps will be discretized as backward differences.
With 1-D implementations, we used
Here, we'll again use
With that in mind, our discretization of the PDE should be relatively straightforward.
As before, solve for the only unknown:
We will solve this equation with the following initial conditions:
and boundary conditions:
from mpl_toolkits.mplot3d import Axes3D ##New Library required for projected 3d plots
import numpy
from matplotlib import pyplot, cm
%matplotlib inline
###variable declarations
nx = 81
ny = 81
nt = 100
c = 1
dx = 2 / (nx - 1)
dy = 2 / (ny - 1)
sigma = .2
dt = sigma * dx
x = numpy.linspace(0, 2, nx)
y = numpy.linspace(0, 2, ny)
u = numpy.ones((ny, nx)) ##create a 1xn vector of 1's
un = numpy.ones((ny, nx)) ##
###Assign initial conditions
##set hat function I.C. : u(.5<=x<=1 && .5<=y<=1 ) is 2
u[int(.5 / dy):int(1 / dy + 1),int(.5 / dx):int(1 / dx + 1)] = 2
###Plot Initial Condition
##the figsize parameter can be used to produce different sized images
fig = pyplot.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
X, Y = numpy.meshgrid(x, y)
surf = ax.plot_surface(X, Y, u[:], cmap=cm.viridis)
To plot a projected 3D result, make sure that you have added the Axes3D library.
from mpl_toolkits.mplot3d import Axes3D
The actual plotting commands are a little more involved than with simple 2d plots.
fig = pyplot.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
surf2 = ax.plot_surface(X, Y, u[:])The first line here is initializing a figure window. The figsize and dpi commands are optional and simply specify the size and resolution of the figure being produced. You may omit them, but you will still require the
fig = pyplot.figure()
The next line assigns the plot window the axes label 'ax' and also specifies that it will be a 3d projection plot. The final line uses the command
plot_surface()
which is equivalent to the regular plot command, but it takes a grid of X and Y values for the data point positions.
The X and Y values that you pass to plot_surface are not the 1-D vectors x and y. In order to use matplotlibs 3D plotting functions, you need to generate a grid of x, y values which correspond to each coordinate in the plotting frame. This coordinate grid is generated using the numpy function meshgrid.
X, Y = numpy.meshgrid(x, y)
To evaluate the wave in two dimensions requires the use of several nested for-loops to cover all of the i's and j's. Since Python is not a compiled language there can be noticeable slowdowns in the execution of code with multiple for-loops. First try evaluating the 2D convection code and see what results it produces.
u = numpy.ones((ny, nx))
u[int(.5 / dy):int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2
for n in range(nt + 1): ##loop across number of time steps
un = u.copy()
row, col = u.shape
for j in range(1, row):
for i in range(1, col):
u[j, i] = (un[j, i] - (c * dt / dx * (un[j, i] - un[j, i - 1])) -
(c * dt / dy * (un[j, i] - un[j - 1, i])))
u[0, :] = 1
u[-1, :] = 1
u[:, 0] = 1
u[:, -1] = 1
fig = pyplot.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
surf2 = ax.plot_surface(X, Y, u[:], cmap=cm.viridis)
Here the same 2D convection code is implemented, but instead of using nested for-loops, the same calculations are evaluated using array operations.
u = numpy.ones((ny, nx))
u[int(.5 / dy):int(1 / dy + 1), int(.5 / dx):int(1 / dx + 1)] = 2
for n in range(nt + 1): ##loop across number of time steps
un = u.copy()
u[1:, 1:] = (un[1:, 1:] - (c * dt / dx * (un[1:, 1:] - un[1:, :-1])) -
(c * dt / dy * (un[1:, 1:] - un[:-1, 1:])))
u[0, :] = 1
u[-1, :] = 1
u[:, 0] = 1
u[:, -1] = 1
fig = pyplot.figure(figsize=(11, 7), dpi=100)
ax = fig.gca(projection='3d')
surf2 = ax.plot_surface(X, Y, u[:], cmap=cm.viridis)
The video lesson that walks you through the details for Step 5 (and onwards to Step 8) is Video Lesson 6 on You Tube:
from IPython.display import YouTubeVideo
YouTubeVideo('tUg_dE3NXoY')from IPython.core.display import HTML
def css_styling():
styles = open("../styles/custom.css", "r").read()
return HTML(styles)
css_styling().warning{
color: rgb( 240, 20, 20 )
}
(The cell above executes the style for this notebook.)