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

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

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

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

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

# hide
# only store the latest time step solution to Told
def fdm_solve(Told):
    Tnew=Told.copy()
    for i in range(1,len(Xvec)-1):
        Tnew[i] = beta * (Told[i+1] - 2*Told[i] + Told[i-1]) + Told[i]
    return Tnew

# 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

# 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

# hide
# temperature close to dike
x_test = 5 + W/2
idx = (np.abs(Xvec - x_test)).argmin()
T_x     = np.empty(Time.shape)
T_x[:]  = np.nan

Compute and visualize everything

# 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

Once Loop Reflect