# heat2Dadv_solution.py
#
# Subitop modelling course, Edinburgh 2017
# Jeroen van Hunen, March 2017
# purpose: calculates 2D heat diffusion
# method: explicit time integration
import numpy as np
import pylab as plt
import time
# Subfunctions:
def twoDadvdiff (fin,vx,vz,dx,dz,kappa,dt):
# Performs 1 diffusion timestep on array fin of size nx X nz
# and discr dx and dz using timestep dt and diffusion coeff kappa
# Fixed essential b.c. are used.
# Calculate 2nd derivatives in x- & z-dir.:
d2fdx2=np.diff(fin,n=2,axis=1)/dx**2
d2fdz2=np.diff(fin,n=2,axis=0)/dz**2
# Initialize a timestep df/dt vector:
dfdt=np.zeros(np.shape(fin))
# Apply diffusion:
# Internal diffusion:
dfdt[1:-1,1:-1] = kappa*d2fdx2[1:-1,:] + kappa*d2fdz2[:,1:-1]
# Natural b.c.'s at side boundaries:
dfdt[1:-1,0] = dfdt[1:-1, 0] + 2*kappa*(fin[1:-1, 1]-fin[1:-1, 0])/dx**2
dfdt[1:-1,-1] = dfdt[1:-1,-1] + 2*kappa*(fin[1:-1,-2]-fin[1:-1,-1])/dx**2
# Advection: upwind approach:
[nz,nx]=np.shape(fin)
for i in range(1,nx-1):
for j in range(0,nz):
if vx[j,i]>=0:
dfdtx=vx[j,i]*(fin[j,i-1]-fin[j,i])/dx
else:
dfdtx=vx[j,i]*(fin[j,i]-fin[j,i+1])/dx
dfdt[j,i]=dfdt[j,i]+dfdtx
for i in range(0,nx):
for j in range(1,nz-1):
if vz[j,i]>=0:
dfdtz=vz[j,i]*(fin[j-1,i]-fin[j,i])/dz
else:
dfdtz=vz[j,i]*(fin[j,i]-fin[j+1,i])/dz
dfdt[j,i]=dfdt[j,i]+dfdtz
# Add dt * df/dt-vector to old solution:
fout=fin+dt*dfdt
return fout
# Main code:
# Initialisation:
# Heat diffusion variables:
kappa = 1e-6 # thermal diffusivity
# Mesh setup:
h = 1000e3 # height of box: 1000 km
w = 1.0*h # box of aspect ratio 1
dx = 2.5e4 # discretization step in meters
dz = 2.5e4
nx = w/dx+1
nz = h/dz+1
dx = w/(nx-1) # Adjust requested dx & dz to fit in equidistant grid space
dz = h/(nz-1)
x = np.linspace(0,w,nx) # array for the finite difference mesh
z = np.linspace(0,h,nz)
[xx,zz] = np.meshgrid(x,z)
# Time variables:
dt_diff = 0.2*dx**2/kappa # timestep in Myrs
nt = 100 # number of tsteps
secinmyr = 1e6*365*24*3600 # amount of seconds in 1 Myr
t = 0 # set initial time to zero
nplot = 5 # plotting interval: plot every nplot timesteps
# Initial condition:
Tm = 1350 # mantle temperature in degC
Ttop = 0 # surface T
Told=Tm*np.ones(np.shape(xx)) # initial T=Tm everywhere ...
# ... except for a block:
Told[(xx>=400e3) & (xx<=600e3) & (zz>=600e3) & (zz<=800e3)]=0.
# Define solid rotation velocity field:
tau = 100*secinmyr # 100-Myr convection overturn period in seconds
v0 = 2*np.pi/tau # angular velocity (in rad/sec)
xh = xx-0.5*h # x-distances to hor middle of mesh (on xx mesh)
zh = zz-0.5*h # z-distances to vert middle of mesh (on xx mesh)
vx = -v0*zh # define solid rotation velocity field: v=vo*r:
vz = v0*xh
# Determine Courant timestep criterion:
vxmax = (abs(vx)).max()
vzmax = (abs(vz)).max()
dt_adv = min(dx/vxmax,dz/vzmax) # advection timestep
dt = 0.5*min(dt_diff, dt_adv) # total timestep
nplot = 1 # Plot every nplot timesteps
# timestepping
starttime = time.time()
for it in range(1,nt):
# numerical solution
#Tnew = twoDadvdiff(Told, xx, zz, xxMP, zzMP, vx, vz, kappa, dt)
Tnew = twoDadvdiff(Told,vx,vz,dx,dz,kappa,dt)
#update time
t=t+dt # in sec
tmyrs=int(t/secinmyr) # and in Myrs
# plot solution:
if (it%nplot==0):
plt.clf()
plt.figure(1) # T-v plot
plt.imshow(Told,
extent=(0,h,0,h),
clim=(0,Tm),
interpolation='bilinear',
cmap='jet')
plt.quiver(xx, h-zz, vx, -vz, units='width')
plt.title('T after '+str(tmyrs)+' Myrs')
plt.pause(0.00005)
# prepare for next time step:
Told = Tnew