tgrounded_cell_fraction.py - 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
---
tgrounded_cell_fraction.py (4193B)
---
1 import PISM
2 import time
3 import numpy as np
4
5 np.set_printoptions(precision=5, suppress=True)
6
7 ctx = PISM.Context()
8
9 ice_density = ctx.config.get_number("constants.ice.density")
10 ocean_density = ctx.config.get_number("constants.sea_water.density")
11
12 mu = ice_density / ocean_density
13
14
15 def allocate_grid(ctx):
16 params = PISM.GridParameters(ctx.config)
17 params.Lx = 1e5
18 params.Ly = 1e5
19 params.Lz = 1000
20 params.Mx = 7
21 params.My = 7
22 params.Mz = 5
23 params.periodicity = PISM.NOT_PERIODIC
24 params.registration = PISM.CELL_CORNER
25 params.ownership_ranges_from_options(ctx.size)
26 return PISM.IceGrid(ctx.ctx, params)
27
28
29 def allocate_storage(grid):
30 ice_thickness = PISM.model.createIceThicknessVec(grid)
31
32 # not used, but needed by GeometryCalculator::compute()
33 surface = PISM.model.createIceSurfaceVec(grid)
34
35 bed_topography = PISM.model.createBedrockElevationVec(grid)
36
37 mask = PISM.model.createIceMaskVec(grid)
38
39 gl_mask = PISM.model.createGroundingLineMask(grid)
40 gl_mask_x = PISM.model.createGroundingLineMask(grid)
41 gl_mask_x.set_name("gl_mask_x")
42 gl_mask_y = PISM.model.createGroundingLineMask(grid)
43 gl_mask_y.set_name("gl_mask_y")
44
45 sea_level = PISM.model.createIceThicknessVec(grid)
46 sea_level.set_name("sea_level")
47
48 return ice_thickness, bed_topography, surface, mask, gl_mask, gl_mask_x, gl_mask_y, sea_level
49
50
51 def compute_mask(sea_level, bed_topography, ice_thickness, mask, surface):
52 gc = PISM.GeometryCalculator(ctx.config)
53 gc.compute(sea_level, bed_topography, ice_thickness, mask, surface)
54
55
56 def print_vec(vec):
57 v0 = vec.allocate_proc0_copy()
58 vec.put_on_proc0(v0.get())
59
60 shape = vec.get_dm().get().sizes
61
62 print(vec.get_name())
63 print(v0.get()[:].reshape(shape, order="f"))
64
65
66 def init(mu, L, sea_level, vec, type="box"):
67 k = {0.0: 8,
68 0.25: 7,
69 0.5: 6,
70 0.75: 5,
71 1.0: 4}
72
73 H0 = (8.0 / k[L]) * (sea_level / mu)
74 H1 = 0.5 * H0
75
76 grid = vec.grid()
77
78 with PISM.vec.Access(nocomm=[vec]):
79 for (i, j) in grid.points():
80 if type == "box" and abs(i - 3) < 2 and abs(j - 3) < 2:
81 vec[i, j] = H0
82 elif type == "cross" and abs(i - 3) < 2 and abs(j - 3) < 2 and (i == 3 or j == 3):
83 vec[i, j] = H0
84 else:
85 vec[i, j] = H1
86
87 if abs(i - 3) >= 3 or abs(j - 3) >= 3:
88 vec[i, j] = 0.0
89
90 vec.update_ghosts()
91
92
93 def grounded_cell_fraction_test():
94
95 # allocation
96 grid = allocate_grid(ctx)
97
98 ice_thickness, bed_topography, surface, mask, gl_mask, gl_mask_x, gl_mask_y, _ = allocate_storage(grid)
99
100 bed_topography.set(0.0)
101
102 # initialization
103 sea_level = 500.0
104 for L in [0.0, 0.25, 0.5, 0.75, 1.0]:
105 init(mu, L, sea_level, ice_thickness, "box")
106
107 compute_mask(sea_level, bed_topography, ice_thickness, mask, surface)
108
109 # computation of gl_mask
110 PISM.compute_grounded_cell_fraction(ice_density, ocean_density, sea_level,
111 ice_thickness, bed_topography, mask, gl_mask,
112 gl_mask_x, gl_mask_y)
113
114 # inspection / comparison
115 print("L = %f" % L)
116 print_vec(mask)
117 print_vec(gl_mask_x)
118 print_vec(gl_mask_y)
119 print_vec(gl_mask)
120
121
122 def new_grounded_cell_fraction_test():
123
124 # allocation
125 grid = allocate_grid(ctx)
126
127 ice_thickness, bed_topography, _, _, gl_mask, _, _, sea_level = allocate_storage(grid)
128
129 # initialization
130 bed_topography.set(0.0)
131 sl = 500.0
132 sea_level.set(sl)
133 for L in [0.0, 0.25, 0.5, 0.75, 1.0]:
134 init(mu, L, sl, ice_thickness, "box")
135
136 # computation of gl_mask
137 PISM.compute_grounded_cell_fraction(ice_density, ocean_density,
138 sea_level,
139 ice_thickness,
140 bed_topography,
141 gl_mask)
142
143 # inspection / comparison
144 print("L = %f" % L)
145 print_vec(gl_mask)
146
147
148 grounded_cell_fraction_test()
149 new_grounded_cell_fraction_test()