Python Exercise 4
Apr 11, 2019
Elliptic and parabolic partial differential equations
Problem 1: Stationary temperature field in beam cross-section
In this problem we shall solve solve problem 1 of last weeks theory exercise, using python.
Figure 1: Beam cross-section.
a) Create a python script that solves for the stationary temperature field in the beam cross-section. Choose \( \Delta x = \Delta y = \frac{1}{N} \), and plot the solution for different values of \( N \). See problem 1 from the theory exercise 4 for details on discretization and treatment of boundary conditions.
Figure 2: Sketch of discretization with \( N=3 \)
Figure 2 shows the discretized grid for \( N=3 \). The nodes with unknown temperature are colored black (\( T_1, \dots T_9 \)), whereas ghost nodes are colored grey and nodes where the temperature is known are colored blue.
If the components of \( \boldsymbol{T} \) is numbered by starting in the lower left part of the grid and the index is increased in the \( x \)-direction first, the function plotSurfaceNeumannDirichlet could be used to plot the solution.
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
def plotSurfaceNeumannDirichlet(Temp, Ttop, Tleft, l, N, nxTicks=4, nyTicks=4):
""" Surface plot of the stationary temperature in quadratic beam cross-section.
Note that the components of T has to be started in the
lower left part of the grid with increasing indexes in the
x-direction first.
Args:
Temp(array): the unknown temperatures, i.e. [T_1 .... T_(NxN)]
Ttop(float): temperature at the top boundary
Tleft(float): temperature at the left boundary
l(float): height/width of the sides
N(int): number of nodes with unknown temperature in x/y direction
nxTicks(int): number of ticks on x-label (default=4)
nyTicks(int): number of ticks on y-label (default=4)
"""
x = np.linspace(0, l, N + 1)
y = np.linspace(0, l, N + 1)
X, Y = np.meshgrid(x, y)
T = np.zeros_like(X)
T[-1,:] = Ttop
T[:,0] = Tleft
k = 1
for j in range(N):
T[j,1:] = Temp[N*(k-1):N*k]
k+=1
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, T, rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
ax.set_zlim(0, Ttop)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('T [$^o$C]')
xticks=np.linspace(0.0, l, nxTicks+1)
ax.set_xticks(xticks)
yticks=np.linspace(0.0, l, nyTicks+1)
ax.set_yticks(yticks)
# python_exercises/pexercise5/problem1.py
import numpy as np
import scipy
import scipy.linalg
import scipy.sparse
import scipy.sparse.linalg
import matplotlib.pylab as plt
import time
from math import sinh
#import matplotlib.pyplot as plt
# Change some default values to make plots more readable on the screen
LNWDT=2; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT
#### start plot
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
def plotSurfaceNeumannDirichlet(Temp, Ttop, Tleft, l, N, nxTicks=4, nyTicks=4):
""" Surface plot of the stationary temperature in quadratic beam cross-section.
Note that the components of T has to be started in the
lower left part of the grid with increasing indexes in the
x-direction first.
Args:
Temp(array): the unknown temperatures, i.e. [T_1 .... T_(NxN)]
Ttop(float): temperature at the top boundary
Tleft(float): temperature at the left boundary
l(float): height/width of the sides
N(int): number of nodes with unknown temperature in x/y direction
nxTicks(int): number of ticks on x-label (default=4)
nyTicks(int): number of ticks on y-label (default=4)
"""
x = np.linspace(0, l, N + 1)
y = np.linspace(0, l, N + 1)
X, Y = np.meshgrid(x, y)
T = np.zeros_like(X)
T[-1,:] = Ttop
T[:,0] = Tleft
k = 1
for j in range(N):
T[j,1:] = Temp[N*(k-1):N*k]
k+=1
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X, Y, T, rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
ax.set_zlim(0, Ttop)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('T [$^o$C]')
xticks=np.linspace(0.0, l, nxTicks+1)
ax.set_xticks(xticks)
yticks=np.linspace(0.0, l, nyTicks+1)
ax.set_yticks(yticks)
#### end plot
def setup_LaplaceNeumann_xy(Ttop, Tleft, N):
mainDiag = -4.*np.ones(N*N)
superDiag1 = np.ones(N*N - 1)
superDiag1[N - 1::N] = 0
subDiag1 = np.ones(N*N - 1)
subDiag1[N - 2::N] = 2
subDiag1[N - 1::N] = 0
superDiagN = np.ones(N*N - N)
superDiagN[0:N] = 2
subDiagN = np.ones(N*N - N)
A = scipy.sparse.diags([mainDiag, superDiag1, subDiag1, superDiagN, subDiagN], [0, 1, -1, N, -N])
b = np.zeros(N*N)
b[0::N] += -Tleft
b[-N:] += -Ttop
return A, b
if __name__ == '__main__':
# Main program
# Set temperature at the top
Ttop=100
Tleft = 50
l=1.0
# Set simulation parameters
#need hx=(1/nx)=hy=(1.5/ny)
N = 20
h=l/N
A, b = setup_LaplaceNeumann_xy(Ttop, Tleft, N)
#print A.toarray()
#print b
Temp = scipy.sparse.linalg.spsolve(A, b)
#print Temp
#print len(Temp)
plotSurfaceNeumannDirichlet(Temp, Ttop, Tleft, l, N)
# figfile='LaPlace_vNeumann.png'
# plt.savefig(figfile, format='png',transparent=True)
plt.show()
Problem 2: Transient temperature field in one-dimensional beam
In this problem we shall solve a parabolic PDE. We are going to find the temperature field along a one-dimensional beam as function of time (transient temperature field). The temperature field is known in the initial state \( t=0 \), and the temperature is known at the boundary.
This problem is goverened by the one-dimensional transient heat conduction equation, $$ \begin{equation} \frac{\partial T}{\partial t}=k \frac{\partial^2 T}{\partial x^2} \label{eq:pde} \end{equation} $$ where \( k \) is the heat conduction coefficient. In this exercise we set \( k=1 \). The beam has length \( L=1 \).
Boundary conditions: \( T(0,t)=100 \) and \( T(1,t)=0 \).
Initial conditions: \( T(x,0)=0 \) Find the temperature field along the beam as function of time. Discretize the length with a suitable number of nodes.
a) Discretize \eqref{eq:pde} with the FTCS scheme and and write a script that solves the problem with given boundary conditions. Explore with different values for \( r \) (the numerical Fourier number), e.g. \( r=0.2 \), \( r=0.5 \) and \( r=1.0 \). Plot the temperature field for \( t=0.001 \), \( t=0.025 \) and \( t=0.4 \).
b) Same as a), but now using the \( \theta \)-scheme instead of the FTCS scheme. Use both \( \theta = \frac{1}{2} \) (Crank) and \( \theta = 1 \) (Laasonen).
This is a suggestion for how to solve the problem. The code is shown below, and the temperature fields are plotted in the following figures. Note that both the FTCS and the \( \theta \)-schemes are coded in two different ways, 1) by using for loops to evaluate expressions for all x-positions and for loops to create RHS/A-matrix, 2) evaluating expressions using arrays and creating Matrices using scipy.sparse.diags. As the number of unknowns grow, the latter is much more efficient.
# python_exercises/pexercise5/problem2.py
import numpy as np
import scipy.sparse.linalg
def FTCS_python(uOld, r):
"""Solve using FTCS looping over all x-positions and evaluating the FTCS-scheme for all points"""
u = np.zeros_like(uOld)
for i in range(1,nx):
u[i] = r*uOld[i-1] + (1.0-2.0*r)*uOld[i]+ r*uOld[i+1]
u[0] = uleft
u[-1] = uright
return u
def FTCS_numpy(uOld, r):
"""Solve using FTCS using a Numpy (array) expression"""
u = np.zeros_like(uOld)
uOld_mid = uOld[1:-1]
uOld_minus = uOld[0:-2]
uOld_plus = uOld[2:]
u[1:-1] = r*uOld_minus + (1.0 - 2.0*r)*uOld_mid + r*uOld_plus
u[0] = uleft
u[-1] = uright
return u
def theta_scheme_python(uOld, r):
"""Solve using the theta scheme, and create matrix/RHS by using for loops"""
u = np.zeros_like(uOld)
A = np.zeros((nx -1, nx-1))
for i in range(nx -1):
for j in range(nx -1):
if i == j:
A[i, j] = (1 + 2.0*r*theta)
if i == j + 1:
A[i, j] = -r*theta
if i == j - 1:
A[i, j] = -r*theta
RHS = np.zeros(nx - 1)
for i in range(1, nx):
RHS[i - 1] = r*(1-theta)*uOld[i - 1] + (1 - 2*r*(1-theta))*uOld[i] + r*(1-theta)*uOld[i + 1]
RHS[0] += r*theta*uOld[0]
RHS[-1] += r*theta*uOld[-1]
# sparse matrix instance
u[1:-1] = np.linalg.solve(A, RHS)
u[0] = uleft
u[-1] = uright
return u
def theta_scheme_numpy(uOld, r):
"""Solve using the theta scheme by creating matrix using scipy.sparse.diags"""
u = np.zeros_like(uOld)
subDiag = -r*theta*np.ones(nx - 2)
mainDiag = (1 + 2.0*r*theta)*np.ones(nx - 1)
superDiag = -r*theta*np.ones(nx - 2)
uOld_mid = uOld[1:-1]
uOld_minus = uOld[0:-2]
uOld_plus = uOld[2:]
RHS = r*(1-theta)*uOld_minus + (1 - 2*r*(1-theta))*uOld_mid + r*(1-theta)*uOld_plus
RHS[0] += r*theta*uOld[0]
RHS[-1] += r*theta*uOld[-1]
A = scipy.sparse.diags([subDiag, mainDiag, superDiag], [-1, 0, 1], format='csc') # sparse matrix instance
u[1:-1] = scipy.sparse.linalg.spsolve(A, RHS)
u[0] = uleft
u[-1] = uright
return u
if __name__ == '__main__':
import matplotlib.pylab as plt
import time as timemodule
LNWDT=2; FNT=15
plt.rcParams['lines.linewidth'] = LNWDT; plt.rcParams['font.size'] = FNT
## Main program starts here
k = 1.0
nx = 100 # number of nodes
L = 1.0 # length of beam
dx = L/nx
x = np.linspace(0, L, nx + 1)
tstart, tend = 0, 0.4 # time period
r = 0.2 # numerical Fourier number
dt = r*dx**2/k # compute timestep based on Fourier number, dx and diffusivity
theta = 0.5 # theta parameter for the theta scheme
uleft, uright, uinitial = 100.0, 0.0, 0.0 # boundary and initial conditions
time = np.arange(tstart, tend + dt, dt)
solvers = [FTCS_numpy, FTCS_python, theta_scheme_numpy, theta_scheme_python]
uSolutions = np.zeros((len (solvers), len(time), nx + 1))
uSolutions[:, 0, 0] = uleft
uSolutions[:, 0, -1] = uright
for k, solver in enumerate(solvers):
cpu_start = timemodule.time()
for n in range(len(time) - 1):
uOld = uSolutions[k, n, :]
uSolutions[k, n + 1, :] = solver(uOld, r)
cpu_end = timemodule.time()
print "cpu-time for " + solver.func_name + " = ", cpu_end - cpu_start
t_plot_list = [0.001, 0.025, 0.4]
l_style=['c', 'ro', 'g^', 'b--']
legendList = []
plt.figure()
for t in t_plot_list:
n = int(round(t/dt))
print time[n]
for k in range(len(solvers)):
plt.plot(x, uSolutions[k, n, :], l_style[k])
for solver in solvers:
legendList.append(solver.func_name)
plt.legend(legendList, frameon=False)
plt.show()