tagetwo.m - 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
---
tagetwo.m (3035B)
---
1 function [x,a1,a2,v,t,a1sum,a2sum] = agetwo(N,Tf)
2 % AGETWO Compare two methods for solving a pure advection problem, which
3 % can be thought of as the "conservative form" age equation
4 % (1) a_t + (v a)_x = 1
5 % on the interval 0 <= x <= 10. Here a = a(x,t) is the age and v = v(x)
6 % is the scalar velocity.
7 %
8 % Note that in 3D ice problems we are solving a_t + div((u,v,w) a) = 1
9 % which is equivalent to a_t + (u,v,w) . grad(a) = 1 because the velocity
10 % vector field (u,v,w) is incompressible. In this example code we are
11 % paying attention to conservative properties of tracer schemes. The age
12 % interpretation is incidental.
13 %
14 % Both schemes are upwinded explicit methods. The first scheme applies to
15 % the equivalent equation
16 % (2) a_t + v a_x = 1 - v' a
17 % and uses classical upwinding for the left side. The second scheme uses
18 % the conservative form (1) stated above. Both schemes use CFL to determine
19 % the time step.
20 %
21 % The particular example problem is on the interval [0,L] = [0,10] and has
22 % v(x) = cos(pi x / (2L)) (1 - 0.7 exp(-(x-3)^2))
23 % which is positive but has v(L) = 0. Roughly, the left end x = 0 is
24 % analogous to the ice surface, and the right end x = L is analogous to
25 % the ice base, and v(L) = 0 is analogous to a frozen base. Thus
26 % PDE equation (1) above
27 % BC a(0,t) = 0
28 % IC a(x,0) = 0
29 %
30 % Try: COMPAREAGE
31
32 L = 10;
33 dx = L/N;
34 x = 0:dx:L;
35 v = vv(x,L);
36 a1 = zeros(size(x)); % zero initial condition
37 a2 = a1;
38
39 dt = dx / max(vv(x,L));
40 M = ceil(Tf / dt);
41 dt = Tf / M;
42 fprintf(' M=%d and dt=%.4f\n',M,dt)
43 t = 0:dt:Tf; % length M+1
44 a1sum = zeros(size(t)); a2sum = a1sum;
45
46 nu = dt / dx;
47
48 for m = 1:M
49 % method one: classic first-order upwinding
50 a1new(1) = 0.0; % left b.c. a(0,t) = 0
51 for j = 2:N+1
52 vj = vv(x(j),L);
53 if vj > 0
54 da = a1(j) - a1(j-1);
55 else
56 da = a1(j+1) - a1(j);
57 end
58 a1new(j) = a1(j) - nu * vj * da + dt - dt * dvv(x(j),L) * a1(j);
59 end
60 a1 = a1new;
61 a1sum(m+1) = (dx/2) * [1 repmat([2],1,N-1) 1] * a1'; % trap rule (finite vol)
62
63 % method two: conservative first-order upwinding
64 % precompute cell-boundary fluxes q = v a
65 q = zeros(1,N);
66 for j = 1:N
67 vmid = 0.5 * (vv(x(j),L) + vv(x(j+1),L));
68 if vmid > 0
69 q(j) = vmid * a2(j);
70 else
71 q(j) = vmid * a2(j+1);
72 end
73 end
74 % update a2
75 % left b.c. a(0,t) = 0 so q(0,t) = 0
76 % q(x=0) = 0 but there *is* a dx/2 cell:
77 a2new(1) = a2(1) - 2 * nu * (q(1) - 0) + dt;
78 for j = 2:N
79 a2new(j) = a2(j) - nu * (q(j) - q(j-1)) + dt;
80 end
81 % need q(x=L); there *is* a dx/2 cell:
82 a2new(N+1) = a2(N+1) - 2 * nu * (vv(x(N+1),L) * a2(N+1) - q(N)) + dt;
83 a2 = a2new;
84 a2sum(m+1) = (dx/2) * [1 repmat([2],1,N-1) 1] * a2'; % trap rule (finite vol)
85 end
86
87 function y = vv(x,L)
88 y = cos(pi * x / (2*L)) .* (1 - 0.7 * exp(-(x-3).^2));
89 end
90
91 function dy = dvv(x,L)
92 ee = exp(-(x-3).^2);
93 cc = cos(pi * x / (2*L));
94 ss = sin(pi * x / (2*L));
95 dy = - (pi/(2*L)) * ss .* (1 - 0.7 * ee) - 0.7 * cc .* (-2 * (x-3)) .* ee;
96 end
97
98 end