Skip to content
Snippets Groups Projects
Commit 9491fb54 authored by Miro Pütz's avatar Miro Pütz
Browse files

Replace 2d_heat_conduction_v2.py 

# no animation anymore
# magma chambers are easier to place and to adjust
# temperature gradient gets calculated now (but needs some work)
# better layout
# final temperature array gets saved in a file
parent 1515e1d7
Branches master
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
2d_heat_conduction_v2.py
+ 86
− 115
View file @ 9491fb54
# 2d time dependent heat equation
# for a non-moving medium with constant conductivity
# and no radiogenic heat production
import numpy as np
import matplotlib
matplotlib.use('GTKAgg')
import gobject
from pylab import *
# 1 for animation, 0 for plot
animation = 0
### Numerical model parameters ###
import matplotlib.pyplot as plt
# -------------------- Numerical model parameters -----------------------------#
# Boundary conditions = constant temperature (gradient)
# Model size [m]
xsize = 1500000 # Profile length [m]
ysize = 35000 # Depth [m]
X_SIZE = 100000 # Profile length [m]
Y_SIZE = 150000 # Depth [m]
# Number of nodes
xnum = 201
ynum = 201
X_NUM = 201
Y_NUM = 201
# Grid step
xstep = float(xsize)/(xnum-1) # [m]
ystep = float(ysize)/(ynum-1) # [m]
ystep2 = ystep**2 # to save cpu
xstep2 = xstep**2
timesteps = 1000 # Number of timesteps
X_STEP = float(X_SIZE)/(X_NUM-1) # [m]
Y_STEP = float(Y_SIZE)/(Y_NUM-1) # [m]
Y_STEP2 = Y_STEP**2 # to save cpu
X_STEP2 = X_STEP**2
TIMESTEPS = 1000 # Number of timesteps
# Material properties
k = 2 # Thermal conductivity [W/(m*K)]
cp = 1000 # Heat cpacity [J/(kg*K)]
rho = 3200 # Density [kg/m**3]
rhocp = rho*cp # Volumetric heat capacity [J/m**3]
kappa = float(k)/rhocp # Thermal diffusivity [m**2/s]
K_uc = 2.75 # Thermal conductivity: upper crust [W/(m*K)]
K_m = 2.75 # mantle
CP = 1000 # Heat cpacity [J/(kg*K)]
RHO = 3200 # Density [kg/m**3]
RHOCP = RHO*CP # Volumetric heat capacity [J/m**3]
KAPPA = float(K_uc)/RHOCP # Thermal diffusivity [m**2/s]
S_uc = 0.2/1.e6 # Heat Production: upper crust [W/m**3]
S_m = 0.05/1.e6 # mantle
QMAX = 20./1000 # Basal heat flow [W/m**2]
ZC = 35000. # Crustal thickness [m]
ZL = 125000. # Lithospheric thickness [m]
# Vectors for nodal points positions
cols = np.arange(0,xsize+1,xstep)
rows = np.arange(0,ysize+1,ystep)
cols = np.linspace(0, X_SIZE, X_NUM)
rows = np.linspace(0, Y_SIZE, Y_NUM)
# Timestep
dtmax = min(xstep,ystep)**2/(3*kappa) # timestep limitation for stability
dt = 0.5*dtmax # timestep [s]
time = (timesteps*dt)/(10**6*365.25*24*3600) # time in Myr
### initial conditions (temperature profile) ###
tmantle = 1200 # mantle temperature in grad or Kelvin???
dtmax = min(X_STEP, Y_STEP)**2/(3*KAPPA) # timestep limitation for stability
dt = 0.9*dtmax # timestep [s]
time = (TIMESTEPS*dt)/(10**6*365.25*24*3600) # time in Myr
# ----------------- initial conditions (temperature profile) ------------------#
# create temperature arrays
# temperature gradient of 25.71 grad/km
temp_gradient = np.arange(0,900+1,((25.71)*ystep)/1000)[:,newaxis]
temp_gradient_2d = temp_gradient*np.ones((len(temp_gradient),xnum))
told = np.ones((ynum,xnum))*tmantle # old temperature (initial condition)
told[:len(temp_gradient),:] = temp_gradient_2d # old temperature (initial condition = gradient)
tnew = told # new temperature
told[60:100,70:130] = tmantle # anomaly
tplot_initial = np.copy(told) # for plotting the initial temperature profile
# Explicit solving
# for the plot
def explicit_2D(told):
for i in xrange(timesteps):
tnew[1:-1, 1:-1] = dt*kappa*( (told[2:, 1:-1]-2*told[1:-1, 1:-1]+told[:-2, 1:-1])/ystep2 +
(told[1:-1, 2:]-2*told[1:-1, 1:-1]+told[1:-1,:-2])/xstep2 )+told[1:-1, 1:-1]
told = np.copy(tnew)
def plots():
plt.figure(1)
X, Y = np.meshgrid(cols,rows)
im = plt.contourf(X/1000, Y/1000, tplot_initial, 100, cmap=cm.hot, interpolation='nearest', origin='lower')
plt.title('Initial Temperature Profile')
plt.xlabel('Profile Distance [km]')
plt.ylabel('Depth [km]')
plt.gca().invert_yaxis()
plt.colorbar(clim(0,1200))
plt.show(block=False)
# Final plot
plt.figure(2)
myplot = plt.contourf(X/1000., Y/1000., told, 100, cmap=cm.hot, interpolation='nearest', origin='lower')
plt.title('After %.3f Myr' %time)
plt.xlabel('Profile Distance [km]')
plt.ylabel('Depth [km]')
plt.gca().invert_yaxis()
plt.colorbar(clim(0,1200))
plt.show()
# for the animation
def explicit_animation_2D(*args):
global told, tnew, t
im.set_array(told)
manager.canvas.draw()
# Uncomment the next two lines to save images as png
# filename ='animation'+str(t)+'.png'
# fig.savefig(filename)
tnew[1:-1, 1:-1] = dt*kappa*( (told[2:, 1:-1]-2*told[1:-1, 1:-1]+told[:-2, 1:-1])/ystep2 +
(told[1:-1, 2:]-2*told[1:-1, 1:-1]+told[1:-1,:-2])/xstep2 )+told[1:-1, 1:-1]
temp_gradient = np.zeros((Y_NUM, X_NUM))
temp_gradient[:47] = (-(S_uc * rows[:47]**2.)/(2.*K_uc) + (QMAX + S_uc * ZL)/K_uc * rows[:47])[:,np.newaxis]
temp_gradient[47:] = (-(S_m * rows[47:]**2.)/(2.*K_m) + (QMAX + S_m * ZL)/K_m * rows[47:])[:,np.newaxis]
told = temp_gradient*np.ones((len(temp_gradient), X_NUM))
tnew = told # new temperature
# create logical array for magma chambers
magma_chamber = 900 # temperature in grad
mask = np.zeros((Y_NUM, X_NUM), dtype=bool)
for i in range(X_NUM):
for j in range(Y_NUM):
if ( ( (i*X_STEP-70000)**2+(j*Y_STEP-40000)**2 <= 13000**2)
or ((i*X_STEP-25000)**2+(j*Y_STEP-50000)**2 <= 10000**2) ):
mask[j,i] = True # True if there is anomaly (magma chamber)
told[mask] = magma_chamber
tplot_initial = np.copy(told) # for plotting the initial temperature profile
# ----------------------------- Explicit solving ------------------------------#
for i in xrange(TIMESTEPS):
tnew[1:-1, 1:-1] = dt*KAPPA*( (told[2:, 1:-1]-2*told[1:-1, 1:-1]+told[:-2, 1:-1])/Y_STEP2 +
(told[1:-1, 2:]-2*told[1:-1, 1:-1]+told[1:-1,:-2])/X_STEP2)+told[1:-1, 1:-1]+S_uc
told = np.copy(tnew)
t+=1
time = (t*dt)/(10**6*365.25*24*3600) # time in Myr)
plt.title('After %.3f Myr' %time)
if t > timesteps:
return False
return True
##### MAIN PROGRAM #####
if animation:
fig = plt.figure(1)
im = plt.imshow(told, cmap=cm.hot, interpolation='nearest', origin='lower')
plt.gca().invert_yaxis()
manager = get_current_fig_manager()
t=0
fig.colorbar(im) # Show the colorbar along the side
# call function explicit_animation_2D until it returns false.
gobject.idle_add(explicit_animation_2D)
show()
else:
explicit_2D(told)
plots()
# Function for colorbar (Not in use)
def fmt(x, pos):
a, b = '{:.2e}'.format(x).split('e')
b = int(b)
return r'${} \times 10^{{{}}}$'.format(a, b)
\ No newline at end of file
if (i*dt)/(10**6*365.25*24*3600) < 2.5: # time after magma chamber is gone
told[mask] = magma_chamber
np.save('temp_array', told) # save the final temp. array in a numpy format
# for the strength calculation
# ------------------------------- Plots ---------------------------------------#
# Initital conditions
plt.figure(1)
X, Y = np.meshgrid(cols,rows)
im = plt.contourf(X/1000, Y/1000, tplot_initial, 100, cmap=plt.cm.hot, interpolation='nearest', origin='lower')
plt.title('Initial Temperature Profile')
plt.xlabel('Profile Distance [km]')
plt.ylabel('Depth [km]')
plt.gca().invert_yaxis()
plt.colorbar(plt.clim(0,np.max(told)))
plt.show(block=False)
# Final plot
plt.figure(2)
myplot = plt.contourf(X/1000., Y/1000., told, 100, cmap=plt.cm.hot, interpolation='nearest', origin='lower')
plt.title('After %.3f Myr' %time)
plt.xlabel('Profile Distance [km]')
plt.ylabel('Depth [km]')
plt.gca().invert_yaxis()
plt.colorbar(plt.clim(0,np.max(told)))
plt.show(block=False)
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment