tcolumn.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tcolumn.py (17056B)
---
1 """Simple verification/regression tests for vertical diffusion and
2 advection of enthalpy in the column. Tests PISM's enthalpy solver.
3
4 """
5
6 import PISM
7 from PISM.util import convert
8 import numpy as np
9
10 log10 = np.log10
11
12 ctx = PISM.Context()
13 unit_system = ctx.unit_system
14 config = ctx.config
15
16 config.set_string("grid.ice_vertical_spacing", "equal")
17
18 k = config.get_number("constants.ice.thermal_conductivity")
19 c = config.get_number("constants.ice.specific_heat_capacity")
20 rho = config.get_number("constants.ice.density")
21 K = k / c
22 # alpha squared
23 alpha2 = k / (c * rho)
24
25 EC = PISM.EnthalpyConverter(ctx.config)
26 pressure = np.vectorize(EC.pressure)
27 cts = np.vectorize(EC.enthalpy_cts)
28
29
30 class EnthalpyColumn(object):
31 "Set up the grid and arrays needed to run column solvers"
32
33 def __init__(self, Mz, dt):
34 self.Lz = 1000.0
35 self.z = np.linspace(0, self.Lz, Mz)
36
37 param = PISM.GridParameters()
38 param.Lx = 1e5
39 param.Ly = 1e5
40 param.z = PISM.DoubleVector(self.z)
41 param.Mx = 3
42 param.My = 3
43 param.Mz = Mz
44
45 param.ownership_ranges_from_options(1)
46
47 self.dt = dt
48
49 self.grid = PISM.IceGrid(ctx.ctx, param)
50 grid = self.grid
51
52 self.enthalpy = PISM.model.createEnthalpyVec(grid)
53
54 self.strain_heating = PISM.model.createStrainHeatingVec(grid)
55
56 self.u, self.v, self.w = PISM.model.create3DVelocityVecs(grid)
57
58 self.sys = PISM.enthSystemCtx(grid.z(), "energy.enthalpy",
59 grid.dx(), grid.dy(), self.dt,
60 config,
61 self.enthalpy,
62 self.u, self.v, self.w,
63 self.strain_heating,
64 EC)
65
66 # zero ice velocity:
67 self.reset_flow()
68 # no strain heating:
69 self.reset_strain_heating()
70
71 def reset_flow(self):
72 self.u.set(0.0)
73 self.v.set(0.0)
74 self.w.set(0.0)
75
76 def reset_strain_heating(self):
77 self.strain_heating.set(0.0)
78
79 def init_column(self):
80 ice_thickness = self.Lz
81 self.sys.init(1, 1, False, ice_thickness)
82
83
84 def diffusion_convergence_rate_time(title, error_func, plot=False):
85 "Compute the convergence rate with refinement in time."
86 dts = 2.0**np.arange(10)
87
88 max_errors = np.zeros_like(dts)
89 avg_errors = np.zeros_like(dts)
90 for k, dt in enumerate(dts):
91 max_errors[k], avg_errors[k] = error_func(False, dt_years=dt, Mz=101)
92
93 p_max = np.polyfit(log10(dts), log10(max_errors), 1)
94 p_avg = np.polyfit(log10(dts), log10(avg_errors), 1)
95
96 if plot:
97 plt.figure()
98 plt.title(title + "\nTesting convergence as dt -> 0")
99 log_plot(dts, max_errors, 'o', "max errors")
100 log_plot(dts, avg_errors, 'o', "avg errors")
101 log_fit_plot(dts, p_max, "max: dt^{:.2}".format(p_max[0]))
102 log_fit_plot(dts, p_avg, "avg: dt^{:.2}".format(p_avg[0]))
103 plt.axis('tight')
104 plt.grid(True)
105 plt.legend(loc="best")
106
107 return p_max[0], p_avg[0]
108
109
110 def diffusion_convergence_rate_space(title, error_func, plot=False):
111 "Compute the convergence rate with refinement in space."
112 Mz = np.array(2.0**np.arange(3, 10), dtype=int)
113 dzs = 1000.0 / Mz
114
115 max_errors = np.zeros_like(dzs)
116 avg_errors = np.zeros_like(dzs)
117 for k, M in enumerate(Mz):
118 T = 1.0
119 max_errors[k], avg_errors[k] = error_func(False,
120 T_final_years=T,
121 dt_years=T,
122 Mz=M)
123
124 p_max = np.polyfit(log10(dzs), log10(max_errors), 1)
125 p_avg = np.polyfit(log10(dzs), log10(avg_errors), 1)
126
127 if plot:
128 plt.figure()
129 plt.title(title + "\nTesting convergence as dz -> 0")
130 log_plot(dzs, max_errors, 'o', "max errors")
131 log_plot(dzs, avg_errors, 'o', "avg errors")
132 log_fit_plot(dzs, p_max, "max: dz^{:.2}".format(p_max[0]))
133 log_fit_plot(dzs, p_avg, "avg: dz^{:.2}".format(p_avg[0]))
134 plt.axis('tight')
135 plt.grid(True)
136 plt.legend(loc="best")
137
138 return p_max[0], p_avg[0]
139
140
141 def exact_DN(L, U0, QL):
142 """Exact solution (and an initial state) for the 'Dirichlet at the base,
143 Neumann at the top' setup."""
144 n = 1
145 lambda_n = 1.0 / L * (np.pi / 2.0 + n * np.pi)
146 a = L * 25.0
147
148 def f(z, t):
149 v = a * np.exp(-lambda_n**2 * alpha2 * t) * np.sin(lambda_n * z)
150 return v + (U0 + QL * z)
151 return f
152
153
154 def errors_DN(plot_results=True, T_final_years=1000.0, dt_years=100, Mz=101):
155 """Test the enthalpy solver with Dirichlet B.C. at the base and
156 Neumann at the top surface.
157 """
158 T_final = convert(T_final_years, "years", "seconds")
159 dt = convert(dt_years, "years", "seconds")
160
161 column = EnthalpyColumn(Mz, dt)
162
163 Lz = column.Lz
164 z = np.array(column.sys.z())
165
166 E_base = EC.enthalpy(230.0, 0.0, EC.pressure(Lz))
167 E_surface = EC.enthalpy(270.0, 0.0, 0.0) - 25*Lz
168 dE_surface = (E_surface - E_base) / Lz
169
170 E_steady = E_base + dE_surface * z
171 Q_surface = K * dE_surface
172
173 E_exact = exact_DN(Lz, E_base, dE_surface)
174
175 with PISM.vec.Access(nocomm=[column.enthalpy,
176 column.u, column.v, column.w,
177 column.strain_heating]):
178 column.sys.fine_to_coarse(E_exact(z, 0), 1, 1, column.enthalpy)
179 column.reset_flow()
180 column.reset_strain_heating()
181
182 t = 0.0
183 while t < T_final:
184 column.init_column()
185
186 column.sys.set_surface_heat_flux(Q_surface)
187 column.sys.set_basal_dirichlet_bc(E_base)
188
189 x = column.sys.solve()
190
191 column.sys.fine_to_coarse(x, 1, 1, column.enthalpy)
192
193 t += dt
194
195 E_exact_final = E_exact(z, t)
196
197 if plot_results:
198 t_years = convert(t, "seconds", "years")
199
200 plt.figure()
201 plt.xlabel("z, meters")
202 plt.ylabel("E, J/kg")
203 plt.step(z, E_exact(z, 0), color="blue", label="initial condition")
204 plt.step(z, E_exact_final, color="green", label="exact solution")
205 plt.step(z, cts(pressure(Lz - z)), "--", color="black", label="CTS")
206 plt.step(z, E_steady, "--", color="green", label="steady state profile")
207 plt.grid(True)
208
209 plt.step(z, x, label="T={} years".format(t_years), color="red")
210
211 plt.legend(loc="best")
212
213 errors = E_exact(z, t) - x
214
215 max_error = np.max(np.fabs(errors))
216 avg_error = np.average(np.fabs(errors))
217
218 return max_error, avg_error
219
220
221 def exact_ND(L, Q0, UL):
222 """Exact solution (and an initial state) for the 'Dirichlet at the base,
223 Neumann at the top' setup."""
224 n = 2
225 lambda_n = 1.0 / L * (-np.pi / 2.0 + n * np.pi)
226 a = L * 25.0
227
228 def f(z, t):
229 v = a * np.exp(-lambda_n**2 * alpha2 * t) * np.sin(lambda_n * (L - z))
230 return v + (UL + Q0 * (z - L))
231 return f
232
233
234 def errors_ND(plot_results=True, T_final_years=1000.0, dt_years=100, Mz=101):
235 """Test the enthalpy solver with Neumann B.C. at the base and
236 Dirichlet B.C. at the top surface.
237 """
238 T_final = convert(T_final_years, "years", "seconds")
239 dt = convert(dt_years, "years", "seconds")
240
241 column = EnthalpyColumn(Mz, dt)
242
243 Lz = column.Lz
244 z = np.array(column.sys.z())
245
246 E_base = EC.enthalpy(240.0, 0.0, EC.pressure(Lz))
247 E_surface = EC.enthalpy(260.0, 0.0, 0.0)
248 dE_base = (E_surface - E_base) / Lz
249
250 E_steady = E_surface + dE_base * (z - Lz)
251 Q_base = - K * dE_base
252
253 E_exact = exact_ND(Lz, dE_base, E_surface)
254
255 with PISM.vec.Access(nocomm=[column.enthalpy,
256 column.u, column.v, column.w,
257 column.strain_heating]):
258 column.sys.fine_to_coarse(E_exact(z, 0), 1, 1, column.enthalpy)
259 column.reset_flow()
260 column.reset_strain_heating()
261
262 t = 0.0
263 while t < T_final:
264 column.init_column()
265
266 column.sys.set_basal_heat_flux(Q_base)
267 column.sys.set_surface_dirichlet_bc(E_surface)
268
269 x = column.sys.solve()
270
271 column.sys.fine_to_coarse(x, 1, 1, column.enthalpy)
272
273 t += dt
274
275 E_exact_final = E_exact(z, t)
276
277 if plot_results:
278 t_years = convert(t, "seconds", "years")
279
280 plt.figure()
281 plt.xlabel("z, meters")
282 plt.ylabel("E, J/kg")
283 plt.step(z, E_exact(z, 0), color="blue", label="initial condition")
284 plt.step(z, E_exact_final, color="green", label="exact solution")
285 plt.step(z, cts(pressure(Lz - z)), "--", color="black", label="CTS")
286 plt.step(z, E_steady, "--", color="green", label="steady state profile")
287 plt.grid(True)
288
289 plt.step(z, x, label="T={} years".format(t_years), color="red")
290
291 plt.legend(loc="best")
292
293 errors = E_exact(z, t) - x
294
295 max_error = np.max(np.fabs(errors))
296 avg_error = np.average(np.fabs(errors))
297
298 return max_error, avg_error
299
300
301 def exact_advection(L, w):
302 "Exact solution of the 'pure advection' problem."
303 C = np.pi / L
304
305 def f(z, t):
306 return np.sin(C * (z - w * t))
307
308 def df(z, t):
309 return C * np.cos(C * (z - w * t))
310 return f, df
311
312
313 def errors_advection_up(plot_results=True, T_final=1000.0, dt=100, Mz=101):
314 """Test the enthalpy solver using a 'pure advection' problem with
315 Neumann (in-flow) B.C. at the base.
316
317 We use Dirichlet B.C. at the surface but they are irrelevant due
318 to upwinding.
319 """
320 w = 1.0
321
322 config.set_number("constants.ice.thermal_conductivity", 0.0)
323 column = EnthalpyColumn(Mz, dt)
324 config.set_number("constants.ice.thermal_conductivity", k)
325
326 Lz = column.Lz
327 z = np.array(column.sys.z())
328
329 E_exact, dE_exact = exact_advection(Lz, w)
330
331 with PISM.vec.Access(nocomm=[column.enthalpy,
332 column.u, column.v, column.w,
333 column.strain_heating]):
334 column.sys.fine_to_coarse(E_exact(z, 0), 1, 1, column.enthalpy)
335 column.reset_flow()
336 column.w.set(w)
337 column.reset_strain_heating()
338
339 t = 0.0
340 while t < T_final:
341 column.init_column()
342
343 column.sys.set_basal_neumann_bc(dE_exact(0, t+dt))
344 column.sys.set_surface_dirichlet_bc(E_exact(Lz, t+dt))
345
346 x = column.sys.solve()
347
348 column.sys.fine_to_coarse(x, 1, 1, column.enthalpy)
349
350 t += dt
351
352 if plot_results:
353 plt.figure()
354 plt.xlabel("z, meters")
355 plt.ylabel("E, J/kg")
356 plt.step(z, E_exact(z, 0), color="blue", label="initial condition")
357 plt.step(z, E_exact(z, t), color="green", label="exact solution")
358 plt.grid(True)
359
360 plt.step(z, x, label="T={} seconds".format(t), color="red")
361
362 plt.legend(loc="best")
363
364 errors = E_exact(z, t) - x
365
366 max_error = np.max(np.fabs(errors))
367 avg_error = np.average(np.fabs(errors))
368
369 return max_error, avg_error
370
371
372 def errors_advection_down(plot_results=True, T_final=1000.0, dt=100, Mz=101):
373 """Test the enthalpy solver using a 'pure advection' problem with
374 Neumann (in-flow) B.C. at the surface.
375
376 We use Dirichlet B.C. at the base but they are irrelevant due
377 to upwinding.
378 """
379 w = -1.0
380
381 config.set_number("constants.ice.thermal_conductivity", 0.0)
382 column = EnthalpyColumn(Mz, dt)
383 config.set_number("constants.ice.thermal_conductivity", k)
384
385 Lz = column.Lz
386 z = np.array(column.sys.z())
387
388 E_exact, dE_exact = exact_advection(Lz, w)
389
390 with PISM.vec.Access(nocomm=[column.enthalpy,
391 column.u, column.v, column.w,
392 column.strain_heating]):
393 column.sys.fine_to_coarse(E_exact(z, 0), 1, 1, column.enthalpy)
394 column.reset_flow()
395 column.w.set(w)
396 column.reset_strain_heating()
397
398 t = 0.0
399 while t < T_final:
400 column.init_column()
401
402 column.sys.set_basal_dirichlet_bc(E_exact(0, t+dt))
403 column.sys.set_surface_neumann_bc(dE_exact(Lz, t+dt))
404
405 x = column.sys.solve()
406
407 column.sys.fine_to_coarse(x, 1, 1, column.enthalpy)
408
409 t += dt
410
411 if plot_results:
412 plt.figure()
413 plt.xlabel("z, meters")
414 plt.ylabel("E, J/kg")
415 plt.step(z, E_exact(z, 0), color="blue", label="initial condition")
416 plt.step(z, E_exact(z, t), color="green", label="exact solution")
417 plt.grid(True)
418
419 plt.step(z, x, label="T={} seconds".format(t), color="red")
420
421 plt.legend(loc="best")
422
423 errors = E_exact(z, t) - x
424
425 max_error = np.max(np.fabs(errors))
426 avg_error = np.average(np.fabs(errors))
427
428 return max_error, avg_error
429
430
431 def advection_convergence_rate_time(title, error_func, plot=False):
432 "Compute the convergence rate with refinement in time."
433 dts = np.linspace(1, 101, 11)
434 Mz = 1001
435 T_final = 200
436
437 max_errors = np.zeros_like(dts)
438 avg_errors = np.zeros_like(dts)
439 for k, dt in enumerate(dts):
440 max_errors[k], avg_errors[k] = error_func(False,
441 T_final=T_final,
442 dt=dt,
443 Mz=Mz)
444
445 p_max = np.polyfit(log10(dts), log10(max_errors), 1)
446 p_avg = np.polyfit(log10(dts), log10(avg_errors), 1)
447
448 if plot:
449 plt.figure()
450 plt.title(title + "\nTesting convergence as dt -> 0")
451 log_plot(dts, max_errors, 'o', "max errors")
452 log_plot(dts, avg_errors, 'o', "avg errors")
453 log_fit_plot(dts, p_max, "max: dt^{:.2}".format(p_max[0]))
454 log_fit_plot(dts, p_avg, "avg: dt^{:.2}".format(p_avg[0]))
455 plt.axis('tight')
456 plt.grid(True)
457 plt.legend(loc="best")
458
459 return p_max[0], p_avg[0]
460
461
462 def advection_convergence_rate_space(title, error_func, plot=False):
463 "Compute the convergence rate with refinement in time."
464 dt = 0.4
465 Mzs = np.linspace(500, 5000, 11, dtype="i")
466 T_final = 10
467
468 dzs = 1000.0 / (Mzs - 1)
469
470 max_errors = np.zeros_like(dzs)
471 avg_errors = np.zeros_like(dzs)
472 for k, M in enumerate(Mzs):
473 max_errors[k], avg_errors[k] = error_func(False,
474 T_final=T_final,
475 dt=dt,
476 Mz=M)
477
478 p_max = np.polyfit(log10(dzs), log10(max_errors), 1)
479 p_avg = np.polyfit(log10(dzs), log10(avg_errors), 1)
480
481 if plot:
482 plt.figure()
483 plt.title(title + "\nTesting convergence as dz -> 0")
484 log_plot(dzs, max_errors, 'o', "max errors")
485 log_plot(dzs, avg_errors, 'o', "avg errors")
486 log_fit_plot(dzs, p_max, "max: dz^{:.2}".format(p_max[0]))
487 log_fit_plot(dzs, p_avg, "avg: dz^{:.2}".format(p_avg[0]))
488 plt.axis('tight')
489 plt.grid(True)
490 plt.legend(loc="best")
491
492 return p_max[0], p_avg[0]
493
494
495 def diffusion_DN_test(plot=False):
496 assert diffusion_convergence_rate_time("Diffusion: Dirichlet at the base, Neumann at the surface",
497 errors_DN, plot)[1] > 0.93
498 assert diffusion_convergence_rate_space("Diffusion: Dirichlet at the base, Neumann at the surface",
499 errors_DN, plot)[1] > 2.0
500
501
502 def diffusion_ND_test(plot=False):
503 assert diffusion_convergence_rate_time("Diffusion: Neumann at the base, Dirichlet at the surface",
504 errors_ND, plot)[1] > 0.93
505 assert diffusion_convergence_rate_space("Diffusion: Neumann at the base, Dirichlet at the surface",
506 errors_ND, plot)[1] > 2.0
507
508
509 def advection_up_test(plot=False):
510 assert advection_convergence_rate_time("Advection: Upward flow",
511 errors_advection_up, plot)[1] > 0.87
512 assert advection_convergence_rate_space("Advection: Upward flow",
513 errors_advection_up, plot)[1] > 0.96
514
515
516 def advection_down_test(plot=False):
517 assert advection_convergence_rate_time("Advection: Downward flow",
518 errors_advection_down, plot)[1] > 0.87
519 assert advection_convergence_rate_space("Advection: Downward flow",
520 errors_advection_down, plot)[1] > 0.96
521
522
523 if __name__ == "__main__":
524 import pylab as plt
525
526 def log_plot(x, y, style, label):
527 plt.plot(log10(x), log10(y), style, label=label)
528 plt.xticks(log10(x), ["{:.2f}".format(t) for t in x])
529 plt.xlabel("log10(spacing)")
530 plt.ylabel("log10(error)")
531
532 def log_fit_plot(x, p, label):
533 plt.plot(log10(x), np.polyval(p, log10(x)), label=label)
534
535 diffusion_ND_test(plot=True)
536 diffusion_DN_test(plot=True)
537 advection_up_test(plot=True)
538 advection_down_test(plot=True)
539 plt.show()