#!/usr/bin/env python
# coding: utf-8
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](https://twitter.com/LorenaABarba)
# 12 steps to Navierâ€“Stokes
# =====
# ***
# Up to now, all of our work has been in one spatial dimension (Steps [1](./01_Step_1.ipynb) to [4](./05_Step_4.ipynb)). 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 $x$ is the variation in the $x$ direction *at constant* $y$.
#
# In 2D space, a rectangular (uniform) grid is defined by the points with coordinates:
#
# $$x_i = x_0 +i \Delta x$$
#
# $$y_i = y_0 +i \Delta y$$
#
# Now, define $u_{i,j} = u(x_i,y_j)$ and apply the finite-difference formulas on either variable $x,y$ *acting separately* on the $i$ and $j$ indices. All derivatives are based on the 2D Taylor expansion of a mesh point value around $u_{i,j}$.
#
# Hence, for a first-order partial derivative in the $x$-direction, a finite-difference formula is:
#
# $$ \frac{\partial u}{\partial x}\biggr\rvert_{i,j} = \frac{u_{i+1,j}-u_{i,j}}{\Delta x}+\mathcal{O}(\Delta x)$$
#
# and similarly in the $y$ direction. Thus, we can write backward-difference, forward-difference or central-difference formulas for Steps 5 to 12. Let's get started!
# Step 5: 2-D Linear Convection
# ----
# ***
# The PDE governing 2-D Linear Convection is written as
#
# $$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x} + c\frac{\partial u}{\partial y} = 0$$
#
# 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 $i$ subscripts to denote movement in space (e.g. $u_{i}^n-u_{i-1}^n$). Now that we have two dimensions to account for, we need to add a second subscript, $j$, to account for all the information in the regime.
#
# Here, we'll again use $i$ as the index for our $x$ values, and we'll add the $j$ subscript to track our $y$ values.
# With that in mind, our discretization of the PDE should be relatively straightforward.
#
# $$\frac{u_{i,j}^{n+1}-u_{i,j}^n}{\Delta t} + c\frac{u_{i, j}^n-u_{i-1,j}^n}{\Delta x} + c\frac{u_{i,j}^n-u_{i,j-1}^n}{\Delta y}=0$$
#
# As before, solve for the only unknown:
#
# $$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)$$
#
# We will solve this equation with the following initial conditions:
#
# $$u(x,y) = \begin{cases}
# \begin{matrix}
# 2\ \text{for} & 0.5 \leq x, y \leq 1 \cr
# 1\ \text{for} & \text{everywhere else}\end{matrix}\end{cases}$$
#
# and boundary conditions:
#
# $$u = 1\ \text{for } \begin{cases}
# \begin{matrix}
# x = 0,\ 2 \cr
# y = 0,\ 2 \end{matrix}\end{cases}$$
# In[1]:
from mpl_toolkits.mplot3d import Axes3D ##New Library required for projected 3d plots
import numpy
from matplotlib import pyplot, cm
get_ipython().run_line_magic('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)
# ### 3D Plotting Notes
# 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.
# ```python
# 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.
#
# ##### Note
#
#
# 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)
#
#
# ### Iterating in two dimensions
# 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.
# In[2]:
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)
# Array Operations
# ----------------
#
# Here the same 2D convection code is implemented, but instead of using nested for-loops, the same calculations are evaluated using array operations.
# In[3]:
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)
# ## Learn More
# The video lesson that walks you through the details for Step 5 (and onwards to Step 8) is **Video Lesson 6** on You Tube:
# In[4]:
from IPython.display import YouTubeVideo
YouTubeVideo('tUg_dE3NXoY')
# In[5]:
from IPython.core.display import HTML
def css_styling():
styles = open("../styles/custom.css", "r").read()
return HTML(styles)
css_styling()
# > (The cell above executes the style for this notebook.)