tnode_types.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 --- tnode_types.py (2052B) --- 1 import PISM 2 import numpy as np 3 import pylab as plt 4 5 # This script tests compute_node_types() using a small 11*11 grid. 6 # Note that icy "tongues" of width 1 cannot be resolved by the Q1 FEM 7 # grid and so nodes corresponding to these are considered "exterior". 8 9 np.set_printoptions(precision=5, suppress=True, linewidth=100) 10 11 ctx = PISM.Context() 12 13 def allocate_grid(ctx): 14 params = PISM.GridParameters(ctx.config) 15 params.Lx = 1e5 16 params.Ly = 1e5 17 params.Lz = 1000 18 params.Mx = 11 19 params.My = 11 20 params.Mz = 5 21 params.registration = PISM.CELL_CORNER 22 params.periodicity = PISM.NOT_PERIODIC 23 params.ownership_ranges_from_options(ctx.size) 24 return PISM.IceGrid(ctx.ctx, params) 25 26 def allocate_storage(grid): 27 ice_thickness = PISM.model.createIceThicknessVec(grid) 28 29 mask = PISM.model.createIceMaskVec(grid) 30 mask.set_name("node_type") 31 32 return ice_thickness, mask 33 34 def spy_vec(vec, value): 35 plt.title(vec.get_name()) 36 plt.imshow(vec.numpy(), interpolation="nearest") 37 38 def init(H, vec): 39 40 grid = vec.grid() 41 42 K = 5 43 R = 2 44 45 with PISM.vec.Access(nocomm=[vec]): 46 for (i, j) in grid.points(): 47 if abs(i - K) < R or abs(j - K) < R: 48 vec[i, j] = H 49 else: 50 vec[i, j] = 0.0 51 52 if abs(i - K) > R or abs(j - K) > R: 53 vec[i, j] = 0.0 54 55 if abs(i - K) < 1 or abs(j - K) < 1: 56 vec[i, j] = H 57 58 if abs(i - K) > R + 2 or abs(j - K) > R + 2: 59 vec[i, j] = 0.0 60 61 vec.update_ghosts() 62 63 def node_type_test(): 64 65 # allocation 66 grid = allocate_grid(ctx) 67 68 ice_thickness, mask = allocate_storage(grid) 69 70 H = 1.0 71 thickness_threshold = 0.5 72 73 # initialization 74 sea_level = 500.0 75 76 init(H, ice_thickness) 77 78 PISM.compute_node_types(ice_thickness, thickness_threshold, mask) 79 80 # inspection / comparison 81 plt.figure() 82 spy_vec(ice_thickness, H) 83 plt.figure() 84 spy_vec(mask, 1.0) 85 plt.show() 86 87 if __name__ == "__main__": 88 node_type_test()