tbed_smoother.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
---
tbed_smoother.py (3981B)
---
1 #!/usr/bin/env python3
2
3 """Simple testing program for Schoof (2003)-type bed smoothing and
4 roughness- parameterization schemes. Allows comparison of computed
5 theta to result from Matlab/Octave code exampletheta.m in
6 src/base/bedroughplay. Also used in PISM software (regression) test.
7 """
8
9
10 import PISM
11 from math import sin, pi
12
13 ctx = PISM.Context()
14 config = ctx.config
15
16
17 def grid():
18 "Create the bed smoother grid."
19 P = PISM.GridParameters(config)
20
21 P.horizontal_size_from_options()
22 P.horizontal_extent_from_options()
23 P.vertical_grid_from_options(config)
24 P.horizontal_extent_from_options()
25 P.ownership_ranges_from_options(ctx.size)
26
27 return PISM.IceGrid(ctx.ctx, P)
28
29
30 def allocate_storage(grid):
31 "Allocate the bed, the smoothed bed, the surface elevation, and theta."
32 topg = PISM.IceModelVec2S()
33 topg.create(grid, "topg", PISM.WITH_GHOSTS, 1)
34 topg.set_attrs("internal", "original topography",
35 "m", "m", "bedrock_altitude", 0)
36
37 topg_smoothed = PISM.IceModelVec2S()
38 topg_smoothed.create(grid, "topg_smoothed", PISM.WITHOUT_GHOSTS, 1)
39 topg_smoothed.set_attrs("internal", "smoothed topography",
40 "m", "m", "bedrock_altitude", 0)
41
42 usurf = PISM.IceModelVec2S()
43 usurf.create(grid, "usurf", PISM.WITH_GHOSTS, 1)
44 usurf.set_attrs("internal", "ice surface elevation",
45 "m", "m", "surface_altitude", 0)
46
47 theta = PISM.IceModelVec2S()
48 theta.create(grid, "theta", PISM.WITH_GHOSTS, 1)
49 theta.set_attrs("internal",
50 "coefficient theta in Schoof (2003) bed roughness parameterization",
51 "", "", "", 0)
52
53 return (topg, topg_smoothed, usurf, theta)
54
55
56 def set_topg(topg):
57 "Initialize the bed topography."
58 grid = topg.grid()
59
60 with PISM.vec.Access(comm=[topg]):
61 for (i, j) in grid.points():
62 x = grid.x(i)
63 y = grid.y(j)
64 topg[i, j] = (400.0 * sin(2.0 * pi * x / 600.0e3) +
65 100.0 * sin(2.0 * pi * (x + 1.5 * y) / 40.0e3))
66
67
68 def set_usurf(usurf):
69 "Initialize the surface elevation."
70 usurf.set(1000.0)
71
72
73 def set_config():
74 "Set configuration parameters."
75
76 config.set_string("grid.periodicity", "none")
77 config.set_string("grid.registration", "corner")
78
79 config.set_number("grid.Mx", 81)
80 config.set_number("grid.My", 81)
81
82 config.set_number("grid.Lx", 1200e3)
83 config.set_number("grid.Ly", 1200e3)
84
85 config.set_number("stress_balance.sia.Glen_exponent", 3.0)
86 config.set_number("stress_balance.sia.bed_smoother.range", 50.0e3)
87
88
89 def smooth(topg, topg_smoothed, usurf, theta):
90 "Smooth the bed topography."
91 grid = topg.grid()
92
93 smoother = PISM.BedSmoother(grid, 1)
94
95 smoother.preprocess_bed(topg)
96
97 smoother.theta(usurf, theta)
98
99 topg_smoothed.copy_from(smoother.smoothed_bed())
100
101
102 def run():
103 "Run the bed smoother using synthetic geometry."
104
105 set_config()
106
107 topg, topg_smoothed, usurf, theta = allocate_storage(grid())
108
109 set_usurf(usurf)
110
111 set_topg(topg)
112
113 smooth(topg, topg_smoothed, usurf, theta)
114
115 return topg, topg_smoothed, usurf, theta
116
117
118 def bed_smoother_test():
119 "Compare the range of topg, topg_smoothed, and theta to stored values"
120
121 topg, topg_smoothed, usurf, theta = run()
122
123 stored_range = {}
124 stored_range["topg"] = [-500.0, 500.0]
125 stored_range["topg_smoothed"] = [-372.9924735817933, 372.9924735817933]
126 stored_range["theta"] = [0.7147300652935706, 0.9884843647808601]
127
128 computed_range = {}
129 for f in [topg, topg_smoothed, theta]:
130 R = f.range()
131 computed_range[f.get_name()] = [R.min, R.max]
132
133 for name in list(stored_range.keys()):
134 computed = computed_range[name]
135 stored = stored_range[name]
136
137 for k in range(2):
138 assert abs(computed[k] - stored[k]) < 1e-16
139
140
141 if __name__ == "__main__":
142 for field in run():
143 field.dump("bed_smoother_%s.nc" % field.get_name())