Implicit descretization

import numpy as np
from scipy.sparse.linalg import factorized
from scipy.sparse import diags, csc_matrix
import matplotlib.pyplot as plt
from matplotlib import animation
from IPython.display import HTML
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

/Users/lruepke/miniconda3/envs/py312_fem_class/lib/python3.12/site-packages/scipy/sparse/_index.py:210: SparseEfficiencyWarning: Changing the sparsity structure of a csc_matrix is expensive. lil and dok are more efficient.
  self._set_arrayXarray(i, j, x)

# Build the matrix
# define diagonals
#lower = -beta * np.ones(nx - 1)
# diag  = ???
# upper = ???^

# puzzle into sparse matrix
#A = diags([?, ?, ?], offsets=[-1, 0, 1], format='csc')

# apply Dirichlet BCs
#A[0, :] = ?; A[0, 0] = ?
#A[-1, :] = ?; A[-1, -1] = ?

# factorize matrix A for fast solves
#solve_A = factorized(A)

print(A[-11:,-12:].toarray())

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

def fdm_solve(Told: np.ndarray) -> np.ndarray:
    # fast per-timestep solve
    return solve_A(Told)

# make initial conditions
T_init  = np.ones(nx)*Trock;              # everything is cold initially
T_init[np.abs(Xvec) <= W/2] = Tmagma   # and hot where the dike is
time    = 0                                     # track the run time

# 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 (°C)')

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

# Initialize Tnew
T=T_init
# Animation function which updates figure data.  This is called sequentially
def animate(i):
    timeLabel.set_text('Time: %.1f day'%(i*dt/day))

    # use global keyword to store the latest solution and update it using fdm_solve function
    global T
    T=fdm_solve(T)
    line.set_data(Xvec, T)

    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