# 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.