import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import time
from IPython.display import display, clear_output
Lxy = 1.0 # domain length Lx=Ly
nTauRun = 0.2 # number of intrinsic timescale to run for
nxy = 22 # number of grid points in x and y directions
α = 1.0 # thermal diffusivity
tFac = 0.5 # timestep size factor.
tau = Lxy**2/α # diffusive timescale based on the whole domain.
tend = nTauRun*tau # run time
dxy = Lxy/(nxy-1) # grid spacing in x and y directions
dt = tFac*dxy**2*α/4 # timestep size
nt = int(np.ceil(tend/dt)) # number of timesteps
nPlt = np.ceil(1/tFac)*10 # how often to plot
f = np.zeros((nxy,nxy)) # solution array
S = np.ones((nxy,nxy)) # source term array
X,Y = np.meshgrid(np.linspace(0,Lxy,nxy), np.linspace(0,Lxy,nxy))
i = np.arange(1,nxy-1)
j = np.arange(1,nxy-1)
ix = np.ix_
for it in range(nt):
f[ix(i,j)] += α*dt/dxy**2*(f[ix(i-1,j)] - 2*f[ix(i,j)] + f[ix(i+1,j)]) + \
α*dt/dxy**2*(f[ix(i,j-1)] - 2*f[ix(i,j)] + f[ix(i,j+1)]) + \
dt*S[ix(i,j)]
if it%nPlt==0:
plt.clf()
plt.contourf(X,Y,f,levels=np.linspace(0,dt*260,21))
plt.colorbar()
plt.gca().set_aspect('equal', adjustable='box')
clear_output(wait=True)
display(plt.gcf())
time.sleep(0.1)
plt.clf()