tssa_coeffs.mac - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tssa_coeffs.mac (2535B)
---
1 /* -*- mode: maxima -*- */
2 isolate_wrt_times:true$
3
4 /* Shifts in x and y directions. */
5 shift_x(expr, d) := op(expr)[args(expr)[1]+d, args(expr)[2]]$
6 shift_y(expr, d) := op(expr)[args(expr)[1], args(expr)[2]+d]$
7 shift[x] : shift_x; shift[y] : shift_y;
8
9 /* weights */
10 weight(var, direction, d) := shift[direction](subst(w, op(var), var), d)$
11
12 /* one-sided differences with weights */
13 d_px(var) := weight(var, x, 1/2) * (shift[x](var, 1) - var)$
14 d_mx(var) := weight(var, x, -1/2) * (var - shift[x](var, -1))$
15
16 d_py(var) := weight(var, y, 1/2) * (shift[y](var, 1) - var)$
17 d_my(var) := weight(var, y, -1/2) * (var - shift[y](var, -1))$
18
19 /* centered differences defined as sums of one-sided differences with weights */
20 D_x(foo) := d_px(foo) + d_mx(foo)$
21 D_y(foo) := d_py(foo) + d_my(foo)$
22
23 load("ssa.mac")$
24
25 denominator : (4 * dx**2 * dy**2)$
26
27 /* Clear the denominator */
28 e1 : ''lhs1 * denominator, expand;
29 e2 : ''lhs2 * denominator, expand;
30
31 /* define weights to be equal to 1 in the interior; give them PIK names otherwise */
32 a_ones : [w[i+1/2,j] = 1,
33 w[i-1/2,j] = 1,
34 w[i,j+1/2] = 1,
35 w[i,j-1/2] = 1,
36 N[i+1/2,j] = c_e,
37 N[i-1/2,j] = c_w,
38 N[i,j+1/2] = c_n,
39 N[i,j-1/2] = c_s,
40 w[i-1/2,j+1] = 1,
41 w[i+1/2,j+1] = 1,
42 w[i+1,j+1/2] = 1,
43 w[i+1,j-1/2] = 1,
44 w[i+1/2,j-1] = 1,
45 w[i-1/2,j-1] = 1,
46 w[i-1,j-1/2] = 1,
47 w[i-1,j+1/2] = 1]$
48
49 a_names : [w[i+1/2,j] = aPP,
50 w[i-1/2,j] = aMM,
51 w[i,j+1/2] = bPP,
52 w[i,j-1/2] = bMM,
53 N[i+1/2,j] = c_e,
54 N[i-1/2,j] = c_w,
55 N[i,j+1/2] = c_n,
56 N[i,j-1/2] = c_s,
57 w[i-1/2,j+1] = aMn,
58 w[i+1/2,j+1] = aPn,
59 w[i+1,j+1/2] = bPe,
60 w[i+1,j-1/2] = bMe,
61 w[i+1/2,j-1] = aPs,
62 w[i-1/2,j-1] = aMs,
63 w[i-1,j-1/2] = bMw,
64 w[i-1,j+1/2] = bPw]$
65
66 if interior then (e1 : at(e1, a_ones), e2 : at(e2, a_ones))
67 else (e1 : at(e1, a_names), e2 : at(e2, a_names))$
68
69 /* Finally compute coefficients: */
70 for m: -1 thru 1 do (for n: -1 thru 1 do (c1u[m,n] : combine(expand(coeff(e1, u[i+m,j+n]) / denominator))));
71 for m: -1 thru 1 do (for n: -1 thru 1 do (c1v[m,n] : combine(expand(coeff(e1, v[i+m,j+n]) / denominator))));
72
73 for m: -1 thru 1 do (for n: -1 thru 1 do (c2u[m,n] : combine(expand(coeff(e2, u[i+m,j+n]) / denominator))));
74 for m: -1 thru 1 do (for n: -1 thru 1 do (c2v[m,n] : combine(expand(coeff(e2, v[i+m,j+n]) / denominator))));
75