tbed_deformation.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_deformation.py (5232B)
---
1 #!/usr/bin/env python3
2
3 import PISM
4 from PISM.util import convert
5 from math import cos, pi
6
7 # Simple testing program for Lingle & Clark bed deformation model.
8 # Runs go for 150,000 years on 63.5km grid with 100a time steps and Z=2 in L&C model.
9 # SCENARIOS: run 'python bed_deformation.py -scenario N' where N=1,2,3,4 as follows
10 # (1) dump ice disc on initially level, non-uplifting land, use only viscous
11 # half-space model:
12 # include_elastic = FALSE, do_uplift = FALSE, H0 = 1000.0
13 # center depth b(0,0) should eventually equilibriate to near
14 # -1000 * (910/3300) = -275.76 m
15 # (2) dump ice disc on initially level, non-uplifting land, use both viscous
16 # half-space model and elastic model
17 # include_elastic = TRUE, do_uplift = FALSE, H0 = 1000.0
18 # (3) never loaded, initially level, uplifting land, use only viscous
19 # half-space model (because elastic model gives no additional when no load):
20 # include_elastic = FALSE, do_uplift = TRUE, H0 = 0.0
21 # (4) dump ice disc on initially level, uplifting land, use both viscous
22 # half-space model and elastic model:
23 # include_elastic = TRUE, do_uplift = TRUE, H0 = 1000.0;
24
25 ctx = PISM.Context()
26 config = ctx.config
27
28 R0 = 1000e3
29
30
31 def initialize_uplift(uplift):
32 "Initialize the uplift field."
33 grid = uplift.grid()
34 peak_uplift = convert(10, "mm/year", "m/second")
35 with PISM.vec.Access(nocomm=[uplift]):
36 for (i, j) in grid.points():
37 r = PISM.radius(grid, i, j)
38 if r < 1.5 * R0:
39 uplift[i, j] = peak_uplift * (cos(pi * (r / (1.5 * R0))) + 1.0) / 2.0
40 else:
41 uplift[i, j] = 0.0
42
43
44 def initialize_thickness(thickness, H0):
45 grid = thickness.grid()
46 with PISM.vec.Access(nocomm=[thickness]):
47 for (i, j) in grid.points():
48 r = PISM.radius(grid, i, j)
49 if r < R0:
50 thickness[i, j] = H0
51 else:
52 thickness[i, j] = 0.0
53
54
55 def allocate(grid):
56 H = PISM.model.createIceThicknessVec(grid)
57 bed = PISM.model.createBedrockElevationVec(grid)
58 uplift = PISM.IceModelVec2S()
59 uplift.create(grid, "uplift", PISM.WITHOUT_GHOSTS)
60 uplift.set_attrs("internal", "bed uplift", "m / second", "m / second", "", 0)
61
62 sea_level = PISM.IceModelVec2S(grid, "sea_level", PISM.WITHOUT_GHOSTS)
63
64 return H, bed, uplift, sea_level
65
66
67 def create_grid():
68 P = PISM.GridParameters(config)
69 P.horizontal_size_from_options()
70 P.horizontal_extent_from_options()
71 P.vertical_grid_from_options(config)
72 P.ownership_ranges_from_options(ctx.size)
73
74 return PISM.IceGrid(ctx.ctx, P)
75
76
77 def run(scenario, plot, pause, save):
78
79 # set grid defaults
80 config.set_number("grid.Mx", 193)
81 config.set_number("grid.My", 129)
82
83 config.set_number("grid.Lx", 3000e3)
84 config.set_number("grid.Ly", 2000e3)
85
86 config.set_number("grid.Mz", 2)
87 config.set_number("grid.Lz", 1000)
88
89 scenarios = {"1": (False, False, 1000.0),
90 "2": (True, False, 1000.0),
91 "3": (False, True, 0.0),
92 "4": (True, True, 1000.0)}
93
94 elastic, use_uplift, H0 = scenarios[scenario]
95
96 print("Using scenario %s: elastic model = %s, use uplift = %s, H0 = %f m" % (scenario, elastic, use_uplift, H0))
97
98 config.set_flag("bed_deformation.lc.elastic_model", elastic)
99
100 grid = create_grid()
101
102 thickness, bed, uplift, sea_level = allocate(grid)
103
104 # set initial geometry and uplift
105 bed.set(0.0)
106 thickness.set(0.0)
107 sea_level.set(0.0)
108 if use_uplift:
109 initialize_uplift(uplift)
110
111 time = ctx.ctx.time()
112
113 time.init(ctx.ctx.log())
114
115 model = PISM.LingleClark(grid)
116
117 model.bootstrap(bed, uplift, thickness, sea_level)
118
119 # now add the disc load
120 initialize_thickness(thickness, H0)
121
122 dt = convert(100, "365 day", "seconds")
123
124 # the time-stepping loop
125 while time.current() < time.end():
126 # don't go past the end of the run
127 dt_current = min(dt, time.end() - time.current())
128
129 model.update(thickness, sea_level, time.current(), dt_current)
130
131 if plot:
132 model.bed_elevation().view(400)
133 model.uplift().view(400)
134
135 print("t = %s years, dt = %s years" % (time.date(), time.convert_time_interval(dt_current, "years")))
136 time.step(dt_current)
137
138 print("Reached t = %s years" % time.date())
139
140 if pause:
141 print("Pausing for 5 seconds...")
142 PISM.PETSc.Sys.sleep(5)
143
144 if save:
145 model.bed_elevation().dump("bed_elevation.nc")
146 model.uplift().dump("bed_uplift.nc")
147
148
149 if __name__ == "__main__":
150 scenario = PISM.OptionKeyword("-scenario", "choose one of 4 scenarios", "1,2,3,4", "1")
151 plot = PISM.OptionBool("-plot", "Plot bed elevation and uplift.")
152 save = PISM.OptionBool("-save", "Save final states of the bed elevation and uplift.")
153 pause = PISM.OptionBool("-pause", "Pause for 5 seconds to look at runtime 2D plots.")
154
155 run(scenario.value(), plot, pause, save)
156
157
158 def scenario1_test():
159 "Test if scenario 1 runs"
160 run("1", False, False, False)
161
162
163 def scenario3_test():
164 "Test if scenario 3 runs"
165 run("3", False, False, False)