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()