Script for solving the 1D diffusion equation - T(t) at x meters from dike
We will solve for the cooling of hot dike that was emplaced into cooler country rock. The method we are using is a simple explicit 1D finite differences scheme. This script also shows how we can record the temperature throught time at a certain distance from the dike country rock interface
First we load some useful libraries
[1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation
from matplotlib.patches import Rectangle
from IPython.display import HTML
Next we define some physics constants
[2]:
L = 100 # Length of modeled domain [m]
Tmagma = 1200 # Temperature of magma [C]
Trock = 300 # Temperature of country rock [C]
kappa = 1e-6 # Thermal diffusivity of rock [m2/s]
W = 5; # Width of dike [m]
day = 60*60*24 # # seconds per day
dt = 1*day # Timestep [s]
and some numerical constants constants
[10]:
nx = 201 # Number of gridpoints in x-direction
nt = 500 # Number of timesteps to compute
Xvec,dx = np.linspace(-L/2, L/2, nx, retstep=True) # X coordinate vector, constant spacing
beta = kappa*dt/dx**2
Time = np.arange(0,nt*dt, dt)
Set the initial conditions
[11]:
T_init = np.ones(nx)*Trock; # everything is cold initially
T_init[np.nonzero(np.abs(Xvec) <= W/2)] = Tmagma # and hot where the dike is
time = 0
Make a function that computes temperature
[12]:
# hide: the code in this cell is hidden by the author
[13]:
# write a function that solves for temperature!
#def fdm_solve(Told):
# Tnew=Told.copy()
# for i in range(1,len(Xvec)-1):
# Tnew[i] = ???
# return Tnew
[14]:
# temperature close to dike
# x_test = 5 + W/2
# idx = ??? find index closest to 5m from dike/country rock interface
# T_x = np.empty(Time.shape) # initialize T_x array for storing temp at this point through time
# T_x[:] = np.nan
[15]:
# hide: the code in this cell is hidden by the author
Compute and visualize everything
[16]:
# First set up the figure, the axis, and the plot element we want to animate
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10,7.5))
# prepare axes
# temperature (space)
ax1.set_xlim(-L/2, L/2)
ax1.set_ylim(0, Tmagma)
ax1.add_patch(Rectangle((-W/2, 0), W, 1200, alpha=0.5, color='tab:red'))
line, = ax1.plot([], [], lw=1)
marker, = ax1.plot(Xvec[idx], Trock, marker='o', markersize=5, color="red", alpha=0.5)
timeLabel=ax1.text(0.02,0.98,'Time: ',transform=ax1.transAxes,va='top')
ax1.set_xlabel('X (m)')
ax1.set_ylabel('Temperature ($^{\circ}$C)')
#temperature (time)
ax2.set_xlim(0, Time[-1]/day)
ax2.set_ylim(0, Tmagma)
evo, = ax2.plot([],[],lw=1)
ax2.set_xlabel('Time (day)')
ax2.set_ylabel('Temperature ($^{\circ}$C)')
tempLabel=ax2.text(0.02,0.98,'Marker temperature: ',transform=ax2.transAxes,va='top')
# Initialization function: plot the background of each frame
def init():
line.set_data([], [])
evo.set_data([],[])
marker.set_data([],[])
return line,evo,marker
# Initialize Tnew
Tnew=T_init
T_x[0] = T_init[idx]
# Animation function which updates figure data. This is called sequentially
def animate(i):
timeLabel._text='Time: %.1f day'%(i*dt/day)
# use global keyword to store the latest solution and update it using fdm_solve function
global Tnew
Tnew=fdm_solve(Tnew)
line.set_data(Xvec, Tnew)
T_x[i] = Tnew[idx]
evo.set_data(Time/day, T_x)
marker.set_data(Xvec[idx], T_x[i])
tempLabel._text='Marker temperature: %.1f ($^{\circ}$C)'%(T_x[i])
return line,evo,marker
# Call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(fig, animate, init_func=init,
frames=nt, interval=30, blit=True)
plt.close(anim._fig)
# Call function to display the animation
#HTML(anim.to_html5_video()) # lower resolution
HTML(anim.to_jshtml()) # higher resolution
[16]: