tmiscellaneous.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
---
tmiscellaneous.py (36128B)
---
1 """This file contains very superficial tests of the PISM Python
2 wrappers. The goal is to be able to detect changes in the API
3 (function signatures, etc), not to test correctness.
4
5 Use with nose (https://pypi.python.org/pypi/nose/) and coverage.py
6 (https://pypi.python.org/pypi/coverage)
7
8 Run this to get a coverage report:
9
10 nosetests --with-coverage --cover-branches --cover-html --cover-package=PISM test/miscellaneous.py
11 """
12
13 import PISM
14 import PISM.testing
15 import sys
16 import os
17 import numpy as np
18 from unittest import TestCase
19
20 ctx = PISM.Context()
21 ctx.log.set_threshold(0)
22
23 def create_dummy_grid():
24 "Create a dummy grid"
25 ctx = PISM.Context()
26 params = PISM.GridParameters(ctx.config)
27 params.ownership_ranges_from_options(ctx.size)
28 return PISM.IceGrid(ctx.ctx, params)
29
30
31 def context_test():
32 "Test creating a new PISM context"
33 ctx = PISM.Context()
34 config = ctx.config
35 us = ctx.unit_system
36 EC = ctx.enthalpy_converter
37
38
39 def context_missing_attribute_test():
40 "Test the handling of missing attributes"
41 ctx = PISM.Context()
42 try:
43 config = ctx.foo # there is no "foo", this should fail
44 return False
45 except AttributeError:
46 return True
47
48
49 def create_grid_test():
50 "Test the creation of the IceGrid object"
51 grid1 = create_dummy_grid()
52
53 grid2 = PISM.model.initGrid(PISM.Context(), 100e3, 100e3, 4000, 11, 11, 21, PISM.CELL_CORNER)
54
55
56 def algorithm_failure_exception_test():
57 "Test the AlgorithmFailureException class"
58 try:
59 raise PISM.AlgorithmFailureException("no good reason")
60 return False # should not be reached
61 except PISM.AlgorithmFailureException as e:
62 print("calling e.reason(): ", e.reason())
63 print("{}".format(e))
64 return True
65
66
67 def printing_test():
68 "Test verbPrintf"
69 ctx = PISM.Context()
70 PISM.verbPrintf(1, ctx.com, "hello %s!\n", "world")
71
72
73 def random_vec_test():
74 "Test methods creating random fields"
75 grid = PISM.IceGrid_Shallow(PISM.Context().ctx, 1e6, 1e6, 0, 0, 61, 31,
76 PISM.NOT_PERIODIC, PISM.CELL_CENTER)
77
78 vec_scalar = PISM.vec.randVectorS(grid, 1.0)
79 vec_vector = PISM.vec.randVectorV(grid, 2.0)
80
81 vec_scalar_ghosted = PISM.vec.randVectorS(grid, 1.0, 2)
82 vec_vector_ghosted = PISM.vec.randVectorV(grid, 2.0, 2)
83
84
85 def vec_metadata_test():
86 "Test accessing IceModelVec metadata"
87 grid = create_dummy_grid()
88
89 vec_scalar = PISM.vec.randVectorS(grid, 1.0)
90
91 m = vec_scalar.metadata()
92
93 m.set_string("units", "kg")
94
95 print(m.get_string("units"))
96
97
98 def vars_ownership_test():
99 "Test passing IceModelVec ownership from Python to C++ (i.e. PISM)."
100 grid = create_dummy_grid()
101 variables = PISM.Vars()
102
103 print("Adding 'thk'...")
104 variables.add(PISM.model.createIceThicknessVec(grid))
105 print("Returned from add_thk()...")
106
107 print("Getting 'thk' from variables...")
108 thk = variables.get("thk")
109 print(thk)
110 thk.begin_access()
111 print("thickness at 0,0 is", thk[0, 0])
112 thk.end_access()
113
114
115 def vec_access_test():
116 "Test the PISM.vec.Access class and IceGrid::points, points_with_ghosts, coords"
117 grid = create_dummy_grid()
118
119 vec_scalar = PISM.vec.randVectorS(grid, 1.0)
120 vec_scalar_ghosted = PISM.vec.randVectorS(grid, 1.0, 2)
121
122 with PISM.vec.Access(comm=[vec_scalar_ghosted], nocomm=vec_scalar):
123 for (i, j) in grid.points_with_ghosts():
124 pass
125
126 with PISM.vec.Access(comm=vec_scalar_ghosted, nocomm=[vec_scalar]):
127 for (i, j) in grid.points():
128 # do something
129 pass
130
131 for (i, j, x, y) in grid.coords():
132 # do something with coordinates
133 pass
134
135 # try with nocomm=None
136 with PISM.vec.Access(comm=vec_scalar_ghosted):
137 pass
138
139
140 def create_modeldata_test():
141 "Test creating the ModelData class"
142 grid = create_dummy_grid()
143 md = PISM.model.ModelData(grid)
144
145 md2 = PISM.model.ModelData(grid, config=grid.ctx().config())
146
147
148 def grid_from_file_test():
149 "Intiialize a grid from a file"
150 grid = create_dummy_grid()
151
152 enthalpy = PISM.model.createEnthalpyVec(grid)
153 enthalpy.set(80e3)
154
155 output_file = "test_grid_from_file.nc"
156 try:
157 pio = PISM.util.prepare_output(output_file)
158
159 enthalpy.write(pio)
160
161 pio = PISM.File(grid.com, output_file, PISM.PISM_NETCDF3, PISM.PISM_READONLY)
162
163 grid2 = PISM.IceGrid.FromFile(ctx.ctx, pio, "enthalpy", PISM.CELL_CORNER)
164 finally:
165 os.remove(output_file)
166
167 def create_special_vecs_test():
168 "Test helpers used to create standard PISM fields"
169 grid = create_dummy_grid()
170
171 usurf = PISM.model.createIceSurfaceVec(grid)
172
173 thk = PISM.model.createIceThicknessVec(grid)
174
175 sea_level = PISM.model.createSeaLevelVec(grid)
176
177 usurfstore = PISM.model.createIceSurfaceStoreVec(grid)
178
179 thkstore = PISM.model.createIceThicknessStoreVec(grid)
180
181 bed = PISM.model.createBedrockElevationVec(grid)
182
183 tauc = PISM.model.createYieldStressVec(grid)
184
185 strainheat = PISM.model.createStrainHeatingVec(grid)
186
187 u, v, w = PISM.model.create3DVelocityVecs(grid)
188
189 hardav = PISM.model.createAveragedHardnessVec(grid)
190
191 enthalpy = PISM.model.createEnthalpyVec(grid)
192
193 age = PISM.model.createAgeVec(grid)
194
195 bmr = PISM.model.createBasalMeltRateVec(grid)
196
197 tillphi = PISM.model.createTillPhiVec(grid)
198
199 basal_water = PISM.model.createBasalWaterVec(grid)
200
201 gl_mask = PISM.model.createGroundingLineMask(grid)
202
203 vel = PISM.model.create2dVelocityVec(grid)
204
205 taudx = PISM.model.createDrivingStressXVec(grid)
206
207 taudy = PISM.model.createDrivingStressYVec(grid)
208
209 vel_misfit_weight = PISM.model.createVelocityMisfitWeightVec(grid)
210
211 cbar = PISM.model.createCBarVec(grid)
212
213 mask = PISM.model.createIceMaskVec(grid)
214
215 bcmask = PISM.model.createBCMaskVec(grid)
216
217 no_model_mask = PISM.model.createNoModelMaskVec(grid)
218
219 zeta_fixed_mask = PISM.model.createZetaFixedMaskVec(grid)
220
221 lon = PISM.model.createLongitudeVec(grid)
222
223 lat = PISM.model.createLatitudeVec(grid)
224
225 # test ModelVecs.add()
226 modeldata = PISM.model.ModelData(grid)
227 vecs = modeldata.vecs
228
229 vecs.add(mask)
230
231 print(vecs)
232 # test getattr
233 vecs.mask
234
235 return True
236
237
238 def pism_vars_test():
239 """Test adding fields to and getting them from pism::Vars."""
240 grid = create_dummy_grid()
241
242 v = grid.variables()
243
244 v.add(PISM.model.createIceThicknessVec(grid))
245
246 # test getting by short name
247 print(v.get("thk").metadata().get_string("units"))
248
249 # test getting by standard name
250 print(v.get("land_ice_thickness").metadata().get_string("units"))
251
252
253 def modelvecs_test():
254 "Test the ModelVecs class"
255
256 grid = create_dummy_grid()
257
258 mask = PISM.model.createIceMaskVec(grid)
259 mask.set(PISM.MASK_GROUNDED)
260
261 modeldata = PISM.model.ModelData(grid)
262 vecs = modeldata.vecs
263
264 vecs.add(mask, "ice_mask", writing=True)
265
266 # use the default name, no writing
267 vecs.add(PISM.model.createIceThicknessVec(grid))
268
269 try:
270 vecs.add(mask, "ice_mask")
271 return False
272 except RuntimeError:
273 # should fail: mask was added already
274 pass
275
276 # get a field:
277 print("get() method: ice mask: ", vecs.get("ice_mask").metadata().get_string("long_name"))
278
279 print("dot notation: ice mask: ", vecs.ice_mask.metadata().get_string("long_name"))
280
281 try:
282 vecs.invalid
283 return False
284 except AttributeError:
285 # should fail
286 pass
287
288 try:
289 vecs.get("invalid")
290 return False
291 except RuntimeError:
292 # should fail
293 pass
294
295 # test __repr__
296 print(vecs)
297
298 # test has()
299 print("Has thickness?", vecs.has("thickness"))
300
301 # test markForWriting
302 vecs.markForWriting("ice_mask")
303
304 vecs.markForWriting(mask)
305
306 vecs.markForWriting("thk")
307
308 # test write()
309 output_file = "test_ModelVecs.nc"
310 try:
311 pio = PISM.util.prepare_output(output_file)
312 pio.close()
313
314 vecs.write(output_file)
315
316 # test writeall()
317 vecs.writeall(output_file)
318 finally:
319 os.remove(output_file)
320
321 def sia_test():
322 "Test the PISM.sia module"
323 ctx = PISM.Context()
324 params = PISM.GridParameters(ctx.config)
325 params.Lx = 1e5
326 params.Ly = 1e5
327 params.Lz = 1000
328 params.Mx = 100
329 params.My = 100
330 params.Mz = 11
331 params.registration = PISM.CELL_CORNER
332 params.periodicity = PISM.NOT_PERIODIC
333 params.ownership_ranges_from_options(ctx.size)
334 grid = PISM.IceGrid(ctx.ctx, params)
335
336 enthalpyconverter = PISM.EnthalpyConverter(ctx.config)
337
338 mask = PISM.model.createIceMaskVec(grid)
339 mask.set(PISM.MASK_GROUNDED)
340
341 thk = PISM.model.createIceThicknessVec(grid)
342 thk.set(1000.0)
343
344 surface = PISM.model.createIceSurfaceVec(grid)
345 surface.set(1000.0)
346
347 bed = PISM.model.createBedrockElevationVec(grid)
348 bed.set(0.0)
349
350 enthalpy = PISM.model.createEnthalpyVec(grid)
351 enthalpy.set(enthalpyconverter.enthalpy(270.0, 0.0, 0.0))
352
353 modeldata = PISM.model.ModelData(grid)
354 modeldata.setPhysics(enthalpyconverter)
355
356 vecs = grid.variables()
357
358 fields = [thk, surface, mask, bed, enthalpy]
359
360 for field in fields:
361 vecs.add(field)
362
363 vel_sia = PISM.sia.computeSIASurfaceVelocities(modeldata)
364
365
366 def util_test():
367 "Test the PISM.util module"
368 grid = create_dummy_grid()
369
370 output_file = "test_pism_util.nc"
371 try:
372 pio = PISM.File(grid.com, output_file, PISM.PISM_NETCDF3, PISM.PISM_READWRITE_MOVE)
373 pio.close()
374
375 PISM.util.writeProvenance(output_file)
376 PISM.util.writeProvenance(output_file, message="history string")
377
378 PISM.util.fileHasVariable(output_file, "data")
379 finally:
380 os.remove(output_file)
381
382 # Test PISM.util.Bunch
383 b = PISM.util.Bunch(a=1, b="string")
384 b.update(c=3.0)
385
386 print(b.a, b["b"], "b" in b, b)
387
388
389 def logging_test():
390 "Test the PISM.logging module"
391 grid = create_dummy_grid()
392
393 import PISM.logging as L
394
395 log_filename = "log.nc"
396 try:
397 PISM.File(grid.com, log_filename, PISM.PISM_NETCDF3, PISM.PISM_READWRITE_MOVE)
398 c = L.CaptureLogger(log_filename)
399
400 L.clear_loggers()
401
402 L.add_logger(L.print_logger)
403 L.add_logger(c)
404
405 L.log("log message\n", L.kError)
406
407 L.logError("error message\n")
408
409 L.logWarning("warning message\n")
410
411 L.logMessage("log message (again)\n")
412
413 L.logDebug("debug message\n")
414
415 L.logPrattle("prattle message\n")
416
417 c.write() # default arguments
418 c.readOldLog()
419 finally:
420 os.remove(log_filename)
421
422 log_filename = "other_log.nc"
423 try:
424 PISM.File(grid.com, log_filename, PISM.PISM_NETCDF3, PISM.PISM_READWRITE_MOVE)
425 c.write(log_filename, "other_log") # non-default arguments
426 finally:
427 os.remove(log_filename)
428
429
430 def column_interpolation_test(plot=False):
431 """Test ColumnInterpolation by interpolating from the coarse grid to the
432 fine grid and back."""
433 import numpy as np
434
435 Lz = 1000.0
436 Mz = 41
437
438 def z_quadratic(Mz, Lz):
439 "Compute levels of a quadratic coarse grid."
440 result = np.zeros(Mz)
441 z_lambda = 4.0
442 for k in range(Mz - 1):
443 zeta = float(k) / (Mz - 1)
444 result[k] = Lz * ((zeta / z_lambda) * (1.0 + (z_lambda - 1.0) * zeta))
445 result[Mz - 1] = Lz
446 return result
447
448 def fine_grid(z_coarse):
449 "Compute levels of the fine grid corresponding to a given coarse grid."
450 Lz = z_coarse[-1]
451 dz = np.min(np.diff(z_coarse))
452 Mz = int(np.ceil(Lz / dz) + 1)
453 dz = Lz / (Mz - 1.0)
454 result = np.zeros(Mz)
455 for k in range(1, Mz):
456 result[k] = z_coarse[0] + k * dz
457
458 return result
459
460 def test_quadratic_interp():
461 z_coarse = z_quadratic(Mz, Lz)
462 f_coarse = (z_coarse / Lz) ** 2
463 z_fine = fine_grid(z_coarse)
464
465 print("Testing quadratic interpolation")
466 return test_interp(z_coarse, f_coarse, z_fine, "Quadratic interpolation")
467
468 def test_linear_interp():
469 z_coarse = np.linspace(0, Lz, Mz)
470 f_coarse = (z_coarse / Lz) ** 2
471 z_fine = fine_grid(z_coarse)
472
473 print("Testing linear interpolation")
474 return test_interp(z_coarse, f_coarse, z_fine, "Linear interpolation")
475
476 def test_interp(z, f, z_fine, title):
477 interp = PISM.ColumnInterpolation(z, z_fine)
478
479 f_fine = interp.coarse_to_fine(f, interp.Mz_fine())
480
481 f_fine_numpy = np.interp(z_fine, z, f)
482
483 f_roundtrip = interp.fine_to_coarse(f_fine)
484
485 if plot:
486 plt.figure()
487 plt.plot(z, f, 'o-', label="original coarse-grid data")
488 plt.plot(z_fine, f_fine, 'o-', label="interpolated onto the fine grid")
489 plt.plot(z, f_roundtrip, 'o-', label="interpolated back onto the coarse grid")
490 plt.plot(z, f_roundtrip - f, 'o-', label="difference after the roundtrip")
491 plt.legend(loc="best")
492 plt.title(title)
493 plt.grid(True)
494
495
496 delta = np.linalg.norm(f - f_roundtrip, ord=1)
497 delta_numpy = np.linalg.norm(f_fine - f_fine_numpy, ord=1)
498 print("norm1(fine_to_coarse(coarse_to_fine(f)) - f) = %f" % delta)
499 print("norm1(PISM - NumPy) = %f" % delta_numpy)
500
501 return delta, delta_numpy
502
503 if plot:
504 import pylab as plt
505
506 linear_delta, linear_delta_numpy = test_linear_interp()
507
508 quadratic_delta, _ = test_quadratic_interp()
509
510 if plot:
511 plt.show()
512
513 if (linear_delta > 1e-12 or linear_delta_numpy > 1e-12 or quadratic_delta > 1e-3):
514 return False
515 return True
516
517
518 def pism_join_test():
519 "Test PISM.join()"
520 assert PISM.join(["one", "two"], ':') == "one:two"
521
522
523 def pism_split_test():
524 "Test PISM.split()"
525 assert PISM.split("one,two,three", ',') == ("one", "two", "three")
526
527
528 def pism_ends_with_test():
529 "Test PISM.ends_with()"
530 assert PISM.ends_with("foo.nc", ".nc") == True
531 assert PISM.ends_with("foo.nc and more text", ".nc") == False
532 assert PISM.ends_with("short_string", "longer_suffix") == False
533
534
535 def linear_interpolation_test(plot=False):
536 "Test linear interpolation code used to regrid fields"
537 import numpy as np
538
539 M_in = 11
540 M_out = 101
541 a = 0.0
542 b = 10.0
543 padding = 1.0
544 x_input = np.linspace(a, b, M_in)
545 x_output = np.sort(((b + padding) - (a - padding)) * np.random.rand(M_out) + (a - padding))
546
547 def F(x):
548 return x * 2.0 + 5.0
549
550 values = F(x_input)
551
552 i = PISM.Interpolation(PISM.LINEAR, x_input, x_output)
553
554 F_interpolated = i.interpolate(values)
555
556 F_desired = F(x_output)
557 F_desired[x_output < a] = F(a)
558 F_desired[x_output > b] = F(b)
559
560 if plot:
561 import pylab as plt
562
563 plt.plot(x_output, F_interpolated, 'o-', color='blue', label="interpolated result")
564 plt.plot(x_output, F_desired, 'x-', color='green', label="desired result")
565 plt.plot(x_input, values, 'o-', color='red', label="input")
566 plt.grid(True)
567 plt.legend(loc="best")
568 plt.show()
569
570 assert np.max(np.fabs(F_desired - F_interpolated)) < 1e-16
571
572
573 def pism_context_test():
574 "Test creating and using a C++-level Context"
575
576 com = PISM.PETSc.COMM_WORLD
577 system = PISM.UnitSystem("")
578
579 logger = PISM.Logger(com, 2)
580
581 config = PISM.DefaultConfig(com, "pism_config", "-config", system)
582 config.init_with_default(logger)
583
584 EC = PISM.EnthalpyConverter(config)
585
586 time = PISM.Time(config, "360_day", system)
587
588 ctx = PISM.cpp.Context(com, system, config, EC, time, logger, "greenland")
589
590 print(ctx.com().Get_size())
591 print(ctx.config().get_number("constants.standard_gravity"))
592 print(ctx.enthalpy_converter().L(273.15))
593 print(ctx.time().current())
594 print(PISM.convert(ctx.unit_system(), 1, "km", "m"))
595 print(ctx.prefix())
596
597
598 def check_flow_law(factory, flow_law_name, EC, stored_data):
599 factory.set_default(flow_law_name)
600 law = factory.create()
601
602 depth = 2000
603 gs = 1e-3
604 sigma = [1e4, 5e4, 1e5, 1.5e5]
605
606 T_pa = [-30, -5, 0, 0]
607 omega = [0.0, 0.0, 0.0, 0.005]
608
609 assert len(T_pa) == len(omega)
610
611 p = EC.pressure(depth)
612 Tm = EC.melting_temperature(p)
613
614 data = []
615 print(" Flow table for %s" % law.name())
616 print("| Sigma | Temperature | Omega | Flow factor |")
617 print("|--------------+--------------+--------------+--------------|")
618 for S in sigma:
619 for Tpa, O in zip(T_pa, omega):
620 T = Tm + Tpa
621 E = EC.enthalpy(T, O, p)
622 F = law.flow(S, E, p, gs)
623 data.append(F)
624
625 print("| %e | %e | %e | %e |" % (S, T, O, F))
626 print("|--------------+--------------+--------------+--------------|")
627 print("")
628
629 data = np.array(data)
630
631 assert np.max(np.fabs(data - stored_data)) < 1e-16
632
633
634 def flowlaw_test():
635 data = {}
636 data["arr"] = [3.91729503e-18, 6.42803396e-17, 1.05746828e-16, 1.05746828e-16,
637 9.79323757e-17, 1.60700849e-15, 2.64367070e-15, 2.64367070e-15,
638 3.91729503e-16, 6.42803396e-15, 1.05746828e-14, 1.05746828e-14,
639 8.81391381e-16, 1.44630764e-14, 2.37930363e-14, 2.37930363e-14]
640 data["arrwarm"] = [1.59798478e-19, 1.04360343e-16, 3.30653997e-16, 3.30653997e-16,
641 3.99496194e-18, 2.60900856e-15, 8.26634991e-15, 8.26634991e-15,
642 1.59798478e-17, 1.04360343e-14, 3.30653997e-14, 3.30653997e-14,
643 3.59546574e-17, 2.34810771e-14, 7.43971492e-14, 7.43971492e-14]
644 data["gk"] = [1.1636334595808724e-16, 6.217445758362754e-15, 2.5309103327753672e-14,
645 2.5309103327753672e-14, 2.5947947614616463e-16, 2.0065832524499375e-14,
646 9.158056141786197e-14, 9.158056141786197e-14, 4.493111202368685e-16,
647 3.469816186746473e-14, 1.6171243121742907e-13, 1.6171243121742907e-13,
648 7.12096200221403e-16, 4.879162291119208e-14, 2.2895389865988545e-13, 2.2895389865988545e-13]
649 data["gpbld"] = [4.65791754e-18, 1.45114704e-16, 4.54299921e-16, 8.66009225e-16,
650 1.16447938e-16, 3.62786761e-15, 1.13574980e-14, 2.16502306e-14,
651 4.65791754e-16, 1.45114704e-14, 4.54299921e-14, 8.66009225e-14,
652 1.04803145e-15, 3.26508084e-14, 1.02217482e-13, 1.94852076e-13]
653 data["hooke"] = [5.26775897e-18, 2.12325906e-16, 5.32397091e-15, 5.32397091e-15,
654 1.31693974e-16, 5.30814764e-15, 1.33099273e-13, 1.33099273e-13,
655 5.26775897e-16, 2.12325906e-14, 5.32397091e-13, 5.32397091e-13,
656 1.18524577e-15, 4.77733287e-14, 1.19789346e-12, 1.19789346e-12]
657 data["isothermal_glen"] = [3.16890000e-16, 3.16890000e-16, 3.16890000e-16, 3.16890000e-16,
658 7.92225000e-15, 7.92225000e-15, 7.92225000e-15, 7.92225000e-15,
659 3.16890000e-14, 3.16890000e-14, 3.16890000e-14, 3.16890000e-14,
660 7.13002500e-14, 7.13002500e-14, 7.13002500e-14, 7.13002500e-14]
661 data["pb"] = [4.65791754e-18, 1.45114704e-16, 4.54299921e-16, 4.54299921e-16,
662 1.16447938e-16, 3.62786761e-15, 1.13574980e-14, 1.13574980e-14,
663 4.65791754e-16, 1.45114704e-14, 4.54299921e-14, 4.54299921e-14,
664 1.04803145e-15, 3.26508084e-14, 1.02217482e-13, 1.02217482e-13]
665
666 ctx = PISM.context_from_options(PISM.PETSc.COMM_WORLD, "flowlaw_test")
667 EC = ctx.enthalpy_converter()
668 factory = PISM.FlowLawFactory("stress_balance.sia.", ctx.config(), EC)
669
670 for flow_law_name, data in data.items():
671 check_flow_law(factory, flow_law_name, EC, np.array(data))
672
673
674 def ssa_trivial_test():
675 "Test the SSA solver using a trivial setup."
676
677 context = PISM.Context()
678
679 L = 50.e3 # // 50km half-width
680 H0 = 500 # // m
681 dhdx = 0.005 # // pure number, slope of surface & bed
682 nu0 = PISM.util.convert(30.0, "MPa year", "Pa s")
683 tauc0 = 1.e4 # // 1kPa
684
685 class TrivialSSARun(PISM.ssa.SSAExactTestCase):
686 def _initGrid(self):
687 self.grid = PISM.IceGrid.Shallow(context.ctx, L, L, 0, 0,
688 self.Mx, self.My, PISM.CELL_CORNER, PISM.NOT_PERIODIC)
689
690 def _initPhysics(self):
691 self.modeldata.setPhysics(context.enthalpy_converter)
692
693 def _initSSACoefficients(self):
694 self._allocStdSSACoefficients()
695 self._allocateBCs()
696
697 vecs = self.modeldata.vecs
698
699 vecs.land_ice_thickness.set(H0)
700 vecs.surface_altitude.set(H0)
701 vecs.bedrock_altitude.set(0.0)
702 vecs.tauc.set(tauc0)
703
704 # zero Dirichler B.C. everywhere
705 vecs.vel_bc.set(0.0)
706 vecs.bc_mask.set(1.0)
707
708 def _initSSA(self):
709 # The following ensure that the strength extension is used everywhere
710 se = self.ssa.strength_extension
711 se.set_notional_strength(nu0 * H0)
712 se.set_min_thickness(4000 * 10)
713
714 # For the benefit of SSAFD on a non-periodic grid
715 self.config.set_flag("ssa.compute_surface_gradient_inward", True)
716
717 def exactSolution(self, i, j, x, y):
718 return [0, 0]
719
720 output_file = "ssa_trivial.nc"
721 try:
722 Mx = 11
723 My = 11
724 test_case = TrivialSSARun(Mx, My)
725 test_case.run(output_file)
726 finally:
727 os.remove(output_file)
728
729 def epsg_test():
730 "Test EPSG to CF conversion."
731 l = PISM.StringLogger(PISM.PETSc.COMM_WORLD, 2)
732 system = PISM.Context().unit_system
733
734 # test supported formats
735 for template in ["{epsg}:{code}",
736 "+init={epsg}:{code}",
737 "+units=m +init={epsg}:{code}",
738 "+init={epsg}:{code} +units=m"]:
739 for epsg in ["EPSG", "epsg"]:
740 for code in [3413, 3031, 3057, 5936, 26710]:
741 string = template.format(epsg=epsg, code=code)
742 print("Trying {}".format(string))
743 l.reset()
744 v = PISM.epsg_to_cf(system, string)
745 v.report_to_stdout(l, 2)
746 print(l.get())
747 print("done.")
748
749 # test that unsupported codes trigger an exception
750 try:
751 v = PISM.epsg_to_cf(system, "+init=epsg:3032")
752 raise AssertionError("should fail with 3032: only 3413 and 3031 are supported")
753 except RuntimeError as e:
754 print("unsupported codes trigger exceptions: {}".format(e))
755
756 # test that an invalid PROJ string (e.g. an EPSG code is not a
757 # number) triggers an exception
758 try:
759 v = PISM.epsg_to_cf(system, "+init=epsg:not-a-number +units=m")
760 raise AssertionError("an invalid PROJ string failed to trigger an exception")
761 except RuntimeError as e:
762 print("invalid codes trigger exceptions: {}".format(e))
763
764 def regridding_test():
765 "Test 2D regridding: same input and target grids."
766 import numpy as np
767
768 ctx = PISM.Context()
769 params = PISM.GridParameters(ctx.config)
770 params.Mx = 3
771 params.My = 3
772 params.ownership_ranges_from_options(1)
773
774 grid = PISM.IceGrid(ctx.ctx, params)
775
776 thk1 = PISM.model.createIceThicknessVec(grid)
777 thk2 = PISM.model.createIceThicknessVec(grid)
778 x = grid.x()
779 x_min = np.min(x)
780 x_max = np.max(x)
781 y = grid.y()
782 y_min = np.min(y)
783 y_max = np.max(y)
784 with PISM.vec.Access(nocomm=[thk1]):
785 for (i, j) in grid.points():
786 F_x = (x[i] - x_min) / (x_max - x_min)
787 F_y = (y[j] - y_min) / (y_max - y_min)
788 thk1[i, j] = (F_x + F_y) / 2.0
789 thk1.dump("thk1.nc")
790
791 thk2.regrid("thk1.nc", critical=True)
792
793 with PISM.vec.Access(nocomm=[thk1, thk2]):
794 for (i, j) in grid.points():
795 v1 = thk1[i, j]
796 v2 = thk2[i, j]
797 if np.abs(v1 - v2) > 1e-12:
798 raise AssertionError("mismatch at {},{}: {} != {}".format(i, j, v1, v2))
799
800 import os
801 os.remove("thk1.nc")
802
803
804
805 def interpolation_weights_test():
806 "Test 2D interpolation weights."
807
808 def interp2d(grid, F, x, y):
809 i_left, i_right, j_bottom, j_top = grid.compute_point_neighbors(x, y)
810 w = grid.compute_interp_weights(x, y)
811
812 i = [i_left, i_right, i_right, i_left]
813 j = [j_bottom, j_bottom, j_top, j_top]
814
815 result = 0.0
816 for k in range(4):
817 result += w[k] * F[j[k], i[k]]
818
819 return result
820
821 Mx = 100
822 My = 200
823 Lx = 20
824 Ly = 10
825
826 grid = PISM.IceGrid_Shallow(PISM.Context().ctx,
827 Lx, Ly, 0, 0, Mx, My,
828 PISM.CELL_CORNER,
829 PISM.NOT_PERIODIC)
830
831 x = grid.x()
832 y = grid.y()
833 X, Y = np.meshgrid(x, y)
834 Z = 2 * X + 3 * Y
835
836 N = 1000
837 np.random.seed(1)
838 x_pts = np.random.rand(N) * (2 * Lx) - Lx
839 y_pts = np.random.rand(N) * (2 * Ly) - Ly
840 # a linear function should be recovered perfectly
841 exact = 2 * x_pts + 3 * y_pts
842
843 result = np.array([interp2d(grid, Z, x_pts[k], y_pts[k]) for k in range(N)])
844
845 np.testing.assert_almost_equal(result, exact)
846
847
848 def vertical_extrapolation_during_regridding_test():
849 "Test extrapolation in the vertical direction"
850 # create a grid with 11 levels, 1000m thick
851 ctx = PISM.Context()
852 params = PISM.GridParameters(ctx.config)
853 params.Lx = 1e5
854 params.Ly = 1e5
855 params.Mx = 3
856 params.My = 3
857 params.Mz = 11
858 params.Lz = 1000
859 params.registration = PISM.CELL_CORNER
860 params.periodicity = PISM.NOT_PERIODIC
861 params.ownership_ranges_from_options(ctx.size)
862
863 z = np.linspace(0, params.Lz, params.Mz)
864 params.z[:] = z
865
866 grid = PISM.IceGrid(ctx.ctx, params)
867
868 # create an IceModelVec that uses this grid
869 v = PISM.IceModelVec3()
870 v.create(grid, "test", PISM.WITHOUT_GHOSTS)
871 v.set(0.0)
872
873 # set a column
874 with PISM.vec.Access(nocomm=[v]):
875 v.set_column(1, 1, z)
876
877 # save to a file
878 v.dump("test.nc")
879
880 # create a taller grid (to 2000m):
881 params.Lz = 2000
882 params.Mz = 41
883 z_tall = np.linspace(0, params.Lz, params.Mz)
884 params.z[:] = z_tall
885
886 tall_grid = PISM.IceGrid(ctx.ctx, params)
887
888 # create an IceModelVec that uses this grid
889 v_tall = PISM.IceModelVec3()
890 v_tall.create(tall_grid, "test", PISM.WITHOUT_GHOSTS)
891
892 # Try regridding without extrapolation. This should fail.
893 try:
894 ctx.ctx.log().disable()
895 v_tall.regrid("test.nc", PISM.CRITICAL)
896 ctx.ctx.log().enable()
897 raise AssertionError("Should not be able to regrid without extrapolation")
898 except RuntimeError as e:
899 pass
900
901 # allow extrapolation during regridding
902 ctx.config.set_flag("grid.allow_extrapolation", True)
903
904 # regrid from test.nc
905 ctx.ctx.log().disable()
906 v_tall.regrid("test.nc", PISM.CRITICAL)
907 ctx.ctx.log().enable()
908
909 # get a column
910 with PISM.vec.Access(nocomm=[v_tall]):
911 column = np.array(v_tall.get_column_vector(1, 1))
912
913 # compute the desired result
914 desired = np.r_[np.linspace(0, 1000, 21), np.zeros(20) + 1000]
915
916 # compare
917 np.testing.assert_almost_equal(column, desired)
918
919 # clean up
920 import os
921 os.remove("test.nc")
922
923 class PrincipalStrainRates(TestCase):
924 def u_exact(self, x, y):
925 "Velocity field for testing"
926 return (np.cos(y) + np.sin(x),
927 np.sin(y) + np.cos(x))
928
929 def eps_exact(self, x, y):
930 "Principal strain rates corresponding to u_exact."
931 u_x = np.cos(x)
932 u_y = - np.sin(y)
933
934 v_x = - np.sin(x)
935 v_y = np.cos(y)
936
937 A = 0.5 * (u_x + v_y)
938 B = 0.5 * (u_x - v_y)
939 Dxy = 0.5 * (v_x + u_y)
940 q = np.sqrt(B**2 + Dxy**2);
941
942 return (A + q, A - q)
943
944 def create_grid(self, Mx):
945 My = Mx + 10 # a non-square grid
946 return PISM.IceGrid.Shallow(self.ctx.ctx,
947 2*np.pi, 2*np.pi,
948 0, 0, int(Mx), int(My), PISM.CELL_CENTER, PISM.NOT_PERIODIC)
949
950 def create_velocity(self, grid):
951 velocity = PISM.IceModelVec2V(grid, "bar", PISM.WITH_GHOSTS)
952 with PISM.vec.Access(nocomm=velocity):
953 for (i, j) in grid.points():
954 u, v = self.u_exact(grid.x(i), grid.y(j))
955 velocity[i, j].u = u
956 velocity[i, j].v = v
957 velocity.update_ghosts()
958
959 return velocity
960
961 def create_cell_type(self, grid):
962 cell_type = PISM.IceModelVec2CellType(grid, "cell_type", PISM.WITH_GHOSTS)
963 cell_type.set(PISM.MASK_GROUNDED)
964 cell_type.update_ghosts()
965
966 return cell_type
967
968 def error(self, Mx):
969 grid = self.create_grid(Mx)
970
971 velocity = self.create_velocity(grid)
972 cell_type = self.create_cell_type(grid)
973 strain_rates = PISM.IceModelVec2(grid, "strain_rates", PISM.WITHOUT_GHOSTS, 0, 2)
974
975 PISM.compute_2D_principal_strain_rates(velocity, cell_type, strain_rates)
976 rates = strain_rates.numpy()
977
978 e1 = rates[:,:,0]
979 e2 = rates[:,:,1]
980
981 # compute e1, e2 using formulas
982 xx, yy = np.meshgrid(grid.x(), grid.y())
983 e1_exact, e2_exact = self.eps_exact(xx, yy)
984
985 delta_e1 = np.max(np.fabs(e1 - e1_exact))
986 delta_e2 = np.max(np.fabs(e2 - e2_exact))
987
988 return Mx, delta_e1, delta_e2
989
990 def setUp(self):
991 self.ctx = PISM.Context()
992
993 def test_principal_strain_rates(self):
994 "Test principal strain rate computation"
995 errors = np.array([self.error(M) for M in [10, 20, 40]])
996
997 Ms = errors[:,0]
998 p_e1 = np.polyfit(np.log(1.0 / Ms), np.log(errors[:,1]), 1)
999 p_e2 = np.polyfit(np.log(1.0 / Ms), np.log(errors[:,2]), 1)
1000
1001 assert p_e1[0] > 1.8
1002 assert p_e2[0] > 1.8
1003
1004 def tearDown(self):
1005 pass
1006
1007 def test_constant_interpolation():
1008 "Piecewise-constant interpolation in Timeseries"
1009
1010 ctx = PISM.Context()
1011 ts = PISM.Timeseries(ctx.com, ctx.unit_system, "test", "time")
1012
1013 t = [0, 1, 2, 3]
1014 y = [2, 3, 1]
1015 N = len(y)
1016
1017 for k in range(N):
1018 ts.append(y[k], t[k], t[k + 1])
1019
1020 # extrapolation on the left
1021 assert ts(t[0] - 1) == y[0]
1022 # extrapolation on the right
1023 assert ts(t[-1] + 1) == y[-1]
1024 # interior intervals
1025 for k in range(N):
1026 T = 0.5 * (t[k] + t[k + 1])
1027
1028 assert ts(T) == y[k], (ts(T), y[k])
1029
1030 def test_timeseries_linear_interpolation():
1031 "Linear interpolation in Timeseries"
1032
1033 ctx = PISM.Context()
1034 ts = PISM.Timeseries(ctx.com, ctx.unit_system, "test", "time")
1035
1036 def f(x):
1037 "a linear function that can be reproduced exactly"
1038 return 2.0 * x - 3.0
1039
1040 N = 4
1041 t = np.arange(N + 1, dtype=np.float64)
1042
1043 for k in range(N):
1044 # use right end point
1045 ts.append(f(t[k + 1]), t[k], t[k + 1])
1046
1047 ts.set_use_bounds(False)
1048
1049 # extrapolation on the left (note that t[0] is not used)
1050 assert ts(t[1] - 1) == f(t[1])
1051
1052 # extrapolation on the right
1053 assert ts(t[-1] + 1) == f(t[-1])
1054
1055 # interior intervals
1056 for k in range(1, N):
1057 dt = t[k + 1] - t[k]
1058 T = t[k] + 0.25 * dt # *not* the midpoint of the interval
1059
1060 assert ts(T) == f(T), (T, ts(T), f(T))
1061
1062 def test_trapezoid_integral():
1063 "Linear integration weights"
1064 x = [0.0, 0.5, 1.0, 2.0]
1065 y = [2.0, 2.5, 3.0, 2.0]
1066
1067 def f(a, b):
1068 return PISM.Interpolation(PISM.LINEAR, x, [a, b]).integral(y)
1069
1070 assert f(0, 2) == 5.0
1071 assert f(0, 0.5) + f(0.5, 1) + f(1, 1.5) + f(1.5, 2) == f(0, 2)
1072 assert f(0, 0.5) + f(0.5, 1.5) + f(1.5, 2) == f(0, 2)
1073 assert f(0, 0.25) + f(0.25, 1.75) + f(1.75, 2) == f(0, 2)
1074
1075 # constant extrapolation:
1076 assert f(-2, -1) == 1 * y[0]
1077 assert f(-2, 0) == 2 * y[0]
1078 assert f(2, 5) == 3 * y[-1]
1079 assert f(3, 4) == 1 * y[-1]
1080
1081 def test_linear_periodic():
1082 "Linear (periodic) interpolation"
1083
1084 period = 1.0
1085 x_p = np.linspace(0, 0.9, 10) + 0.05
1086 y_p = np.sin(2 * np.pi * x_p)
1087
1088 # grid with end points added (this makes periodic interpolation unnecessary)
1089 x = np.r_[0, x_p, 1]
1090 y = np.sin(2 * np.pi * x)
1091
1092 # target grid
1093 xx = np.linspace(0, 1, 21)
1094
1095 yy_p = PISM.Interpolation(PISM.LINEAR_PERIODIC, x_p, xx, period).interpolate(y_p)
1096
1097 yy = PISM.Interpolation(PISM.LINEAR, x, xx).interpolate(y)
1098
1099 np.testing.assert_almost_equal(yy, yy_p)
1100
1101 def test_nearest_neighbor():
1102 "Nearest neighbor interpolation"
1103 x = [-1, 1]
1104 y = [0, 1]
1105
1106 xx = np.linspace(-0.9, 0.9, 10)
1107 yy = np.ones_like(xx) * (xx > 0)
1108
1109 zz = PISM.Interpolation(PISM.NEAREST, x, xx).interpolate(y)
1110
1111 np.testing.assert_almost_equal(yy, zz)
1112
1113 class AgeModel(TestCase):
1114 def setUp(self):
1115 self.output_file = "age.nc"
1116 self.grid = self.create_dummy_grid()
1117
1118 def create_dummy_grid(self):
1119 "Create a dummy grid"
1120 params = PISM.GridParameters(ctx.config)
1121 params.ownership_ranges_from_options(ctx.size)
1122 return PISM.IceGrid(ctx.ctx, params)
1123
1124 def test_age_model_runs(self):
1125 "Check if AgeModel runs"
1126 ice_thickness = PISM.model.createIceThicknessVec(self.grid)
1127
1128 u = PISM.IceModelVec3(self.grid, "u", PISM.WITHOUT_GHOSTS)
1129 v = PISM.IceModelVec3(self.grid, "v", PISM.WITHOUT_GHOSTS)
1130 w = PISM.IceModelVec3(self.grid, "w", PISM.WITHOUT_GHOSTS)
1131
1132 ice_thickness.set(4000.0)
1133 u.set(0.0)
1134 v.set(0.0)
1135 w.set(0.0)
1136
1137 model = PISM.AgeModel(self.grid, None)
1138 input_options = PISM.process_input_options(ctx.com, ctx.config)
1139 model.init(input_options)
1140
1141 inputs = PISM.AgeModelInputs(ice_thickness, u, v, w)
1142
1143 dt = PISM.util.convert(1, "years", "seconds")
1144
1145 model.update(0, dt, inputs)
1146
1147 model.age().dump(self.output_file)
1148
1149 def tearDown(self):
1150 os.remove(self.output_file)
1151
1152 def checksum_test():
1153 "Check if a small change in an IceModelVec affects checksum() output"
1154 grid = PISM.testing.shallow_grid(Mx=101, My=201)
1155
1156 v = PISM.IceModelVec2S(grid, "dummy", PISM.WITHOUT_GHOSTS)
1157 v.set(1e15)
1158
1159 old_checksum = v.checksum()
1160
1161 with PISM.vec.Access(nocomm=v):
1162 for (i, j) in grid.points():
1163 if i == 0 and j == 0:
1164 v[i, j] += 1
1165
1166 assert old_checksum != v.checksum()
1167
1168 class ForcingOptions(TestCase):
1169 def setUp(self):
1170 # store current configuration parameters
1171 self.config = PISM.DefaultConfig(ctx.com, "pism_config", "-config", ctx.unit_system)
1172 self.config.init_with_default(ctx.log)
1173 self.config.import_from(ctx.config)
1174
1175 self.filename = "forcing-options-input.nc"
1176 PISM.util.prepare_output(self.filename)
1177
1178 def test_without_file(self):
1179 "ForcingOptions: xxx.file is not set"
1180 ctx.config.set_string("input.file", self.filename)
1181
1182 opt = PISM.ForcingOptions(ctx.ctx, "surface.given")
1183
1184 assert opt.filename == self.filename
1185 assert opt.period == 0
1186 assert opt.reference_time == 0
1187
1188 def test_with_file(self):
1189 "ForcingOptions: xxx.file is set"
1190 ctx.config.set_string("surface.given.file", self.filename)
1191
1192 opt = PISM.ForcingOptions(ctx.ctx, "surface.given")
1193
1194 assert opt.filename == self.filename
1195 assert opt.period == 0
1196 assert opt.reference_time == 0
1197
1198 def test_without_file_and_without_input_file(self):
1199 "ForcingOptions: xxx.file is not set and -i is not set"
1200 try:
1201 opt = PISM.ForcingOptions(ctx.ctx, "surface.given")
1202 assert False, "failed to stop with an error message"
1203 except RuntimeError:
1204 pass
1205
1206 def test_negative_period(self):
1207 "ForcingOptions: negative period"
1208 ctx.config.set_string("input.file", self.filename)
1209 ctx.config.set_number("surface.given.period", -1)
1210
1211 try:
1212 opt = PISM.ForcingOptions(ctx.ctx, "surface.given")
1213 assert False, "failed to catch negative xxx.period"
1214 except RuntimeError:
1215 pass
1216
1217 def tearDown(self):
1218 # reset configuration parameters
1219 ctx.config.import_from(self.config)
1220
1221 os.remove(self.filename)