Deriving coefficients for the implicit scheme

The ice sheet energy balance model uses an implicit scheme to solve the heat equation for $N$ layers. It uses the Crank-Nicholson scheme to discretise the equations. For an equation in one space dimension $x$

$\frac{df}{dt} = F$,

the Crank-Nicholson scheme discretises the equation as

$\frac{f^{i+1}_x - f^{i}_x}{\Delta t} = 0.5\left [ F^{i+1}_x + F^{i}_x \right ]$

where the superscript is time and the subscript is space.



In [4]:

    
from sympy import *
init_printing()



In [141]:

    
tnew_x = Symbol('T^{i+1}_x')
tnew_xprev = Symbol('T^{i+1}_{x-1}')
tnew_xafter = Symbol('T^{i+1}_{x+1}')

told_x = Symbol('T^{i}_x')
told_xprev = Symbol('T^{i}_{x-1}')
told_xafter = Symbol('T^{i}_{x+1}')

u_x = Symbol('\kappa_x')
u_xprev = Symbol('\kappa_{x-1}')
u_xafter = Symbol('\kappa_{x+1}')

delta_t = Symbol('\Delta t')
delta_x = Symbol('\Delta x')



In [142]:

    
told_x, u_xprev, tnew_xafter, delta_x









    Out[142]:





$$\left ( T^{i}_x, \quad \kappa_{x-1}, \quad T^{i+1}_{x+1}, \quad \Delta x\right )$$



In [143]:

    
lhs = (tnew_x - told_x)/delta_t



In [144]:

    
lhs # The time derivative









    Out[144]:





$$\frac{T^{i+1}_x - T^{i}_x}{\Delta t}$$



In [158]:

    
rhs_new = 0.5*(u_x*(tnew_xprev - 2*tnew_x + tnew_xafter)/delta_x**2 + 
               ((tnew_x - tnew_xprev)/(delta_x))*((u_x - u_xprev)/(delta_x)))



In [159]:

    
rhs_old = 0.5*(u_x*(told_xprev - 2*told_x + told_xafter)/delta_x**2 + 
               ((told_x - told_xprev)/(delta_x))*((u_x - u_xprev)/(delta_x)))



In [160]:

    
rhs_new, rhs_old # The two parts of the crank-nicholson RHS.









    Out[160]:





$$\left ( \frac{0.5 \kappa_x}{\Delta x^{2}} \left(- 2 T^{i+1}_x + T^{i+1}_{x+1} + T^{i+1}_{x-1}\right) + \frac{0.5}{\Delta x^{2}} \left(T^{i+1}_x - T^{i+1}_{x-1}\right) \left(\kappa_x - \kappa_{x-1}\right), \quad \frac{0.5 \kappa_x}{\Delta x^{2}} \left(- 2 T^{i}_x + T^{i}_{x+1} + T^{i}_{x-1}\right) + \frac{0.5}{\Delta x^{2}} \left(T^{i}_x - T^{i}_{x-1}\right) \left(\kappa_x - \kappa_{x-1}\right)\right )$$



In [161]:

    
expr = lhs - rhs_new - rhs_old



In [162]:

    
expr









    Out[162]:





$$- \frac{0.5 \kappa_x}{\Delta x^{2}} \left(- 2 T^{i+1}_x + T^{i+1}_{x+1} + T^{i+1}_{x-1}\right) - \frac{0.5 \kappa_x}{\Delta x^{2}} \left(- 2 T^{i}_x + T^{i}_{x+1} + T^{i}_{x-1}\right) - \frac{0.5}{\Delta x^{2}} \left(T^{i+1}_x - T^{i+1}_{x-1}\right) \left(\kappa_x - \kappa_{x-1}\right) - \frac{0.5}{\Delta x^{2}} \left(T^{i}_x - T^{i}_{x-1}\right) \left(\kappa_x - \kappa_{x-1}\right) + \frac{T^{i+1}_x - T^{i}_x}{\Delta t}$$



In [163]:

    
poly_form = Poly(expr, tnew_x, tnew_xafter, tnew_xprev, told_x, told_xafter, told_xprev)



In [167]:

    
poly_form









    Out[167]:





$$\operatorname{Poly}{\left( \frac{1.0 T^{i+1}_x}{\Delta t \Delta x^{2}} \left(0.5 \Delta t \kappa_x + 0.5 \Delta t \kappa_{x-1} + 1.0 \Delta x^{2}\right) - \frac{0.5 T^{i+1}_{x+1} \kappa_x}{\Delta x^{2}} - \frac{0.5 T^{i+1}_{x-1} \kappa_{x-1}}{\Delta x^{2}} + \frac{1.0 T^{i}_x}{\Delta t \Delta x^{2}} \left(0.5 \Delta t \kappa_x + 0.5 \Delta t \kappa_{x-1} - 1.0 \Delta x^{2}\right) - \frac{0.5 T^{i}_{x+1} \kappa_x}{\Delta x^{2}} - \frac{0.5 T^{i}_{x-1} \kappa_{x-1}}{\Delta x^{2}}, T^{i+1}_x, T^{i+1}_{x+1}, T^{i+1}_{x-1}, T^{i}_x, T^{i}_{x+1}, T^{i}_{x-1}, domain=\mathbb{R}\left(\Delta t, \Delta x, \kappa_x, \kappa_{x-1}\right) \right)}$$

The coefficients for the $i+1$ temperature (predicted) are



In [165]:

    
(poly_form.coeff_monomial(tnew_xprev)*delta_t).simplify(), (poly_form.coeff_monomial(tnew_x)*delta_t).simplify(), (poly_form.coeff_monomial(tnew_xafter)*delta_t).simplify()









    Out[165]:





$$\left ( - \frac{0.5 \Delta t \kappa_{x-1}}{\Delta x^{2}}, \quad \frac{1}{\Delta x^{2}} \left(0.5 \Delta t \kappa_x + 0.5 \Delta t \kappa_{x-1} + 1.0 \Delta x^{2}\right), \quad - \frac{0.5 \Delta t \kappa_x}{\Delta x^{2}}\right )$$

The coefficients for the $i$ temperature (current) are



In [166]:

    
-(poly_form.coeff_monomial(told_xprev)*delta_t).simplify(), (poly_form.coeff_monomial(told_x)*-delta_t).simplify(), -(poly_form.coeff_monomial(told_xafter)*delta_t).simplify()









    Out[166]:





$$\left ( \frac{0.5 \Delta t \kappa_{x-1}}{\Delta x^{2}}, \quad \frac{1}{\Delta x^{2}} \left(- 0.5 \Delta t \kappa_x - 0.5 \Delta t \kappa_{x-1} + 1.0 \Delta x^{2}\right), \quad \frac{0.5 \Delta t \kappa_x}{\Delta x^{2}}\right )$$

These are the coefficients for the diagonals of the sparse matrix. The need to be multiplied by $\frac{1}{\rho C}$ where $\rho, C$ are the vertically varying density and specific heat to represent snow and ice.

The $\kappa$ values are staggered from the $T$ values, since $T$ represents the interface temperatures and $\kappa$ represents the mid_level values. K_mid and K_bar in the code therefore interpolate $\kappa$ values to the interface.