In [ ]:
%pylab inline
Se også notebook for den eksplisitte løsningen samt 1D implisitt løsning. Vi ønsker å løse varmelikningen implisitt,
$\frac{\delta u}{\delta t} = \kappa \frac{\delta^2u}{\delta x^2}$,
fortsatt på kartesisk grid, men nå i 2D. Diskretisert implisitt med hensyn på tid får vi
$\frac{1}{\Delta t}\left(u_i^{k} - u_i^{k-1}\right) = \kappa\left[ \frac{1}{\Delta x^2} \left(u^k_{i, j-1} - 2u_{i, j}^k + u_{i, j+1}^k \right) + \frac{1}{\Delta y^2} \left(u^k_{i-1, j} - 2u_{i, j}^k + u_{i+1, j}^k \right) \right]$,
og tilhørende stensil
$u_{i, j}^{k-1} = -s_x u^k_{i, j-1} +2s_x u_{i, j}^k - s_x u_{i, j+1}^k -s_y u^k_{i-1, j} +2s_y u_{i, j}^k - s_y u_{i+1, j}^k + u_{i, j}^k, \quad s_x=\frac{\kappa\Delta t}{\Delta x^2}, s_y=\frac{\kappa\Delta t}{\Delta y^2}$.
Om vi reorganiserer, kan vi skrive dette som
$u_{i, j}^{k-1} =
eller
$u_{i, j}^{k-1} = - \text{LeftFlux} - \text{RightFlux} - \text{DownFlux} - \text{UpFlux} + u_{i, j}^k$,
som illustrerer sammenhengen mellom matrisa for likningssystemet og fluksene samt verdiene $u_i^n$ for løsningen.
In [ ]:
import numpy as np
import time
import matplotlib.pyplot as plt
from IPython.display import clear_output
from mpl_toolkits.mplot3d.axes3d import Axes3D
# Grid spacing
dx = 0.5
dy = 0.5
# Heat diffusion constant
kappa = 0.3
# Maximum time step size (constrained by CFL)
#dt = (dx*dx)/(4*kappa) # ok when dx==dy
dt = (dx*dx*dy*dy)/(dx*dx+dy*dy) / (2*kappa) # if dx!=dy
#dt *= 10
#dt = 1000;
# Number of grid cells
n = 5
# Boundary conditions
ub = np.zeros( (n, n) )
ub[0] = np.linspace(0.5, 1.5, n) # Top row
ub[n-1] = np.linspace(0.5, 1.5, n) # Bottom row
ub[:, 0] = np.linspace(0.5, 0.5, n) # Left column
ub[:, n-1] = np.linspace(1.5, 1.5, n) # Right column
# Grid (fikse denne, bruke dx og dy)
x = np.outer( np.linspace(1, 1, n), np.linspace(0.0, 1.0, n) )
y = np.outer( np.linspace(1.0, 0.0, n), np.linspace(1, 1, n) )
In [ ]:
# Initial temperatures
u = np.zeros( (n, n) )
u[0] = ub[0]
u[n-1] = ub[n-1]
u[:, 0] = ub[:, 0]
u[:, n-1] = ub[:, n-1]
# Data structure containing neigbhour information, nbr[i, j, k] contains the displacement of neigbhour cell number
# k of cell (i, j) (with coordinates x[j], y[i]), i.e. the neighbour cell is cell (i, j) + nbr[i, j, k].
# The number of neighbours of cell (i, j) is stored in num_of_nbr[i, j].
nbr = np.zeros( (n, n, 4, 2), dtype=int )
num_of_nbr = np.zeros( (n, n), dtype=int )
for i in range(n):
for j in range(n):
# left
if (j>0):
nbr[i, j, num_of_nbr[i, j]] = np.array( [0, -1] )
num_of_nbr[i, j] += 1
# right
if (j<n-1):
nbr[i, j, num_of_nbr[i, j]] = np.array( [0, 1] )
num_of_nbr[i, j] += 1
# up (in indices, "down" in the sense that y is less for the neighbour)
if (i>0):
nbr[i, j, num_of_nbr[i, j]] = np.array( [-1, 0] )
num_of_nbr[i, j] += 1
# down
if (i<n-1):
nbr[i, j, num_of_nbr[i, j]] = np.array( [1, 0] )
num_of_nbr[i, j] += 1
In [ ]:
# print or plot?
pl = True
In [ ]:
def assembleA(rx, ry):
A = np.zeros( (n*n, n*n) )
# Assemble the matrix one cell at a time.
for i in range(n):
for j in range(n):
ii = i*n + j # Index of cell (i, j)
A[ii, ii] = 1;
# Looping over neighbours.
for k in range(num_of_nbr[i, j]):
di, dj = nbr[i, j, k]
nbr_ii = (i+di)*n + (j+dj)
# flux = r*(u[i, j] - u[i+di, j+dj])
if ( di == 0 ):
r = rx
else:
r = ry
A[ii, ii] += r
A[ii, nbr_ii] += -r
# Apply Dirichlet boundary conditions, # u_i^n = u_i^n-1, for i=0 and i=n-1.
for i in range(0, n, n-1):
for j in range(n):
ii = i*n + j;
A[ii] = np.zeros(n*n)
A[ii, ii] = 1
for j in range(0, n, n-1):
for i in range(n):
ii = i*n + j;
A[ii] = np.zeros(n*n)
A[ii, ii] = 1
return A;
In [ ]:
sx = kappa * dt / (dx*dx)
sy = kappa * dt / (dy*dy)
In [ ]:
# Rendering a matrix A and its inverse, for visual inspection
fig2 = plt.figure(figsize=(16, 4))
plt.subplot(121); # Denne går inn i current figure, som nå er fig2?
plt.imshow(assembleA(sx, sy), interpolation='nearest');
plt.colorbar();
plt.subplot(122)
plt.imshow(pylab.inv(assembleA(sx, sy)), interpolation='nearest')
qqq=plt.colorbar() # Trenger å putte den i en variabel for å unngå utskrift av "handle"
#plt.sca(ax) # Setting current axis to 'ax' and current figure to the parent figure of 'ax'.
#plt.figure(fig.number); # Setting the current figure
In [ ]:
# Preparing figure for later rendering
if (pl):
fig = plt.figure(figsize=(10, 6))
ax = fig.add_subplot(1, 1, 1, projection='3d')
ax.view_init(50, -110)
plt.close()
In [ ]:
# Initializing with zero in the interior and the Dirichlet boundary conditions. Rendering the initial condition.
u = np.zeros( (n, n) )
u[0] = ub[0]
u[n-1] = ub[n-1]
u[:, 0] = ub[:, 0]
u[:, n-1] = ub[:, n-1]
if (pl):
ax.cla()
clear_output()
#ax.set_title( "Time = " + str(dt*k) )
ax.grid(False) # When 'True', the grid "shines through". How to avoid this?
ax.plot_surface( x, y, u, rstride=1, cstride=1, antialiased=False, linewidth=1, shade=True );
ax.set_xlabel("x")
ax.set_ylabel("y");
ax.set_zlabel("temp");
ax.set_zlim(0.5, 1.5)
display(fig)
In [ ]:
# Setting up the matrix A, and then performing som iterations
A = assembleA(sx, sy)
for k in range(0, 20):
# Indexing convention for the unknowns u^{n+1}_{i, j} ( =v[i, j] ) in the solution vector: v[i*n+j].
v = np.linalg.solve(A, reshape(u, (n*n, 1)));
u, v = reshape(v, (n, n)), u;
if (not(pl)):
print "u =\n", u
else:
ax.cla()
clear_output()
ax.set_title( "Time = " + str(dt*k) )
ax.grid(False) # When 'True', the grid "shines through". How to avoid this?
ax.plot_surface( x, y, u, rstride=1, cstride=1, antialiased=False, linewidth=1, shade=True );
ax.set_xlabel("x")
ax.set_ylabel("y");
ax.set_zlabel("temp");
ax.set_zlim(0.5, 1.5)
display(fig)
plt.close()
In [ ]: