tdiscretization.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
---
tdiscretization.mac (1708B)
---
1 linel : 200$
2
3 /* The definition of mu. */
4 mu_eq : mu = dt / (2 * dz);
5 /* The definition of R. */
6 R_eq[k] := R[k] = K[k] * dt / (rho * dz^2);
7 /* Equations used to get rid of K in the discretization. */
8 Rplus : solve(R_eq[k+half], K[k+half])[1]$
9 Rminus : solve(R_eq[k-half], K[k-half])[1]$
10 /* Define one over dz to make the equation prettier. */
11 one_over_dz : 1 / dz$
12
13 /* Upwinding */
14 Up(a, b) := if b >= 0 then op(a)[k] - op(a)[k - 1] else op(a)[k + 1] - op(a)[k];
15
16 /* The discretization itself. */
17 arg(x) := args(x)[1];
18 Delta(x) := (op(x)[arg(x)+1] - op(x)[arg(x)-1]);
19 delta_p(x) := (op(x)[arg(x)+1] - op(x)[arg(x)]);
20 delta_m(x) := (op(x)[arg(x)] - op(x)[arg(x)-1]);
21
22 eq :
23 rho * ((E[k] - Eold[k]) / dt
24 + lambda * w[k] * 'Delta(E[k]) / (2 * dz)
25 + (1 - lambda) * w[k] * 'Up(E[k], w[k]) / dz)
26 - 'one_over_dz * (K[k+half] * 'delta_p(E[k]) / dz - K[k-half] * 'delta_m(E[k]) / dz) = Phi[k];
27
28 /* The discretization with the upwinding "operator" expanded according
29 to the sign of w */
30 eq1 : eq, nouns;
31
32 /* Multiply by dt / rho and expand. */
33 eq2 : eq * dt / rho, nouns, expand$
34
35 /* Simplify */
36 eq3 : ev(lhs(eq2), Rplus, Rminus, solve(mu_eq, dt)[1]) = rhs(eq2)$
37
38 /* Factor out E[k-1], E[k], E[k+1] */
39 eq4 : facsum(eq3, E[k-1], E[k], E[k+1]);
40
41 l : facsum(coeff(lhs(eq4), E[k-1]), w[k])$
42 d : facsum(coeff(lhs(eq4), E[k]), w[k])$
43 u : facsum(coeff(lhs(eq4), E[k+1]), w[k])$
44 b : rhs(eq4) - (lhs(eq4) - E[k-1]*l - E[k]*d - E[k+1]*u), expand$
45
46 /* Test that the lower, diagonal, upper diagonal entries and the
47 b give an equation equivalent to eq4. test should simplify to 0. */
48 test : lhs(eq4) - rhs(eq4) - (E[k-1]*l + E[k]*d + E[k+1]*u - b), expand;
49
50 l_eq : L[k] = l;
51 d_eq : D[k] = d;
52 u_eq : U[k] = u;
53 b_eq : B[k] = b;