Implicit descretization

import numpy as np
import numpy.matlib
from scipy.linalg import solve
from scipy.sparse import spdiags
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML
import matplotlib as mpl
from matplotlib.patches import Rectangle

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      = 500*day         # Timestep [s]

nx      = 101 # 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

Now we need an array to store the coefficient matrix. We use the scipy spdiags command for this, which puts the coefficients in “data” onto the diagonals defined in “diags”. The result is a sparse matrix, which only stores the coefficients and not all the zeros. Boundary conditions are done as outlined in the script. Note how the setting of boundary conditions “destroys” the symmetry of the matrix. This is numerically not very smart; later we will learn how to do this better.

# hide: the code in this cell is hidden by the author

A[-10:nx-1,-10:nx-1]

array([[ 87.4, -43.2,   0. ,   0. ,   0. ,   0. ,   0. ,   0. ,   0. ],
       [-43.2,  87.4, -43.2,   0. ,   0. ,   0. ,   0. ,   0. ,   0. ],
       [  0. , -43.2,  87.4, -43.2,   0. ,   0. ,   0. ,   0. ,   0. ],
       [  0. ,   0. , -43.2,  87.4, -43.2,   0. ,   0. ,   0. ,   0. ],
       [  0. ,   0. ,   0. , -43.2,  87.4, -43.2,   0. ,   0. ,   0. ],
       [  0. ,   0. ,   0. ,   0. , -43.2,  87.4, -43.2,   0. ,   0. ],
       [  0. ,   0. ,   0. ,   0. ,   0. , -43.2,  87.4, -43.2,   0. ],
       [  0. ,   0. ,   0. ,   0. ,   0. ,   0. , -43.2,  87.4, -43.2],
       [  0. ,   0. ,   0. ,   0. ,   0. ,   0. ,   0. , -43.2,  87.4]])

# Now we need to setup a coefficient matrix!

# build the coefficient matrix
# data  = (np.ones((nx,1))*np.array([-beta, (1+2*beta), -beta ])).T
# diags = np.array([-1, 0, 1])
# A     = spdiags(???).toarray()

# and add boundary conditions
# A[0,0]       = ???
# A[0,1]       = ???
# A[nx-1,nx-1] = ???
# A[nx-1,nx-2] = ???

# make 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                                     # track the run time

# We only store the latest time step
def fdm_solve(Told):
    Tnew=solve(A,Told)
    return Tnew

# First set up the figure, the axis, and the plot element we want to animate
fig = plt.figure(figsize=(10,5))
ax = plt.axes(xlim=(-L/2, L/2), ylim=(0, Tmagma))
ax.add_patch(Rectangle((-W/2, 0), W, 1200, alpha=0.5, color='tab:red'))
line, = ax.plot([], [], lw=1)
timeLabel=ax.text(0.02,0.98,'Time: ',transform=ax.transAxes,va='top')
ax.set_xlabel('X (m)')
ax.set_ylabel('Temperature ($^{\circ}$C)')

# Initialization function: plot the background of each frame
def init():
    line.set_data(Xvec, T_init)
    return line,

# Initialize Tnew
Tnew=T_init
# 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)

    return line,

# 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

B     = spdiags(data, diags, nx, nx).toarray()

print(B)

[[ 87.4 -43.2   0.  ...   0.    0.    0. ]
 [-43.2  87.4 -43.2 ...   0.    0.    0. ]
 [  0.  -43.2  87.4 ...   0.    0.    0. ]
 ...
 [  0.    0.    0.  ...  87.4 -43.2   0. ]
 [  0.    0.    0.  ... -43.2  87.4 -43.2]
 [  0.    0.    0.  ...   0.  -43.2  87.4]]