tfrontal_melt_models.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
---
tfrontal_melt_models.py (6470B)
---
1 #!/usr/bin/env python3
2 """
3 Tests of PISM's frontal melt models.
4 """
5
6 import PISM
7 from PISM.util import convert
8 import sys, os, numpy
9 from unittest import TestCase
10 import netCDF4
11
12 config = PISM.Context().config
13
14 # reduce the grid size to speed this up
15 config.set_number("grid.Mx", 3)
16 config.set_number("grid.My", 3)
17 config.set_number("grid.Mz", 5)
18
19 log = PISM.Context().log
20 # silence models' initialization messages
21 log.set_threshold(1)
22
23 options = PISM.PETSc.Options()
24
25 seconds_per_day = 86400
26
27 def create_geometry(grid):
28 geometry = PISM.Geometry(grid)
29
30 geometry.latitude.set(0.0)
31 geometry.longitude.set(0.0)
32
33 geometry.bed_elevation.set(0.0)
34 geometry.sea_level_elevation.set(0.0)
35
36 geometry.ice_thickness.set(1000.0)
37 geometry.ice_area_specific_volume.set(0.0)
38
39 geometry.ensure_consistency(0.0)
40
41 return geometry
42
43 def sample(vec):
44 return vec.numpy()[0,0]
45
46 def create_dummy_forcing_file(filename, variable_name, units, value):
47 f = netCDF4.Dataset(filename, "w")
48 f.createDimension("time", 1)
49 t = f.createVariable("time", "d", ("time",))
50 t.units = "seconds"
51 delta_T = f.createVariable(variable_name, "d", ("time",))
52 delta_T.units = units
53 t[0] = 0.0
54 delta_T[0] = value
55 f.close()
56
57 def create_grid():
58 "Create a dummy grid"
59 ctx = PISM.Context()
60 params = PISM.GridParameters(ctx.config)
61 params.ownership_ranges_from_options(ctx.size)
62 return PISM.IceGrid(ctx.ctx, params)
63
64 def create_given_input_file(filename, grid, temperature, mass_flux):
65 PISM.util.prepare_output(filename)
66
67 T = PISM.IceModelVec2S(grid, "shelfbtemp", PISM.WITHOUT_GHOSTS)
68 T.set_attrs("climate", "shelf base temperature", "Kelvin", "Kelvin", "", 0)
69 T.set(temperature)
70 T.write(filename)
71
72 M = PISM.IceModelVec2S(grid, "shelfbmassflux", PISM.WITHOUT_GHOSTS)
73 M.set_attrs("climate", "shelf base mass flux", "kg m-2 s-1", "kg m-2 s-1", "", 0)
74 M.set(mass_flux)
75 M.write(filename)
76
77 def check(vec, value):
78 "Check if values of vec are almost equal to value."
79 numpy.testing.assert_almost_equal(sample(vec), value)
80
81 def check_model(model, melt_rate):
82 check(model.frontal_melt_rate(), melt_rate)
83
84 def constant_test():
85 "Model Constant"
86
87 # compute mass flux
88 melt_rate = config.get_number("frontal_melt.constant.melt_rate", "m second-1")
89
90 grid = create_grid()
91 geometry = create_geometry(grid)
92
93 inputs = PISM.FrontalMeltInputs()
94 water_flux = PISM.IceModelVec2S(grid, "water_flux", PISM.WITHOUT_GHOSTS)
95 water_flux.set(0.0)
96 inputs.geometry = geometry
97 inputs.subglacial_water_flux = water_flux
98
99 model = PISM.FrontalMeltConstant(grid)
100
101 model.init(geometry)
102 model.update(inputs, 0, 1)
103
104 check_model(model, melt_rate)
105
106 assert model.max_timestep(0).infinite()
107
108 class DischargeRoutingTest(TestCase):
109
110 def frontal_melt(self, h, q_sg, TF):
111 """
112 h: water depth, meters
113 q_sg: subglacial water flux, m / day
114 TF: thermal forcing, Celsius
115
116 Returns the melt rate in m / day
117 """
118 alpha = config.get_number("frontal_melt.routing.power_alpha")
119 beta = config.get_number("frontal_melt.routing.power_beta")
120 A = config.get_number("frontal_melt.routing.parameter_a")
121 B = config.get_number("frontal_melt.routing.parameter_b")
122
123 return (A * h * q_sg ** alpha + B) * TF ** beta
124
125 def setUp(self):
126 self.depth = 1000.0 # meters
127 self.potential_temperature = 4.0 # Celsius
128 self.water_flux = 10.0 # m / day
129
130 self.grid = create_grid()
131
132 self.theta = PISM.IceModelVec2S(self.grid, "theta_ocean", PISM.WITHOUT_GHOSTS)
133 self.theta.set(self.potential_temperature)
134
135 self.Qsg = PISM.IceModelVec2S(self.grid, "subglacial_water_flux",
136 PISM.WITHOUT_GHOSTS)
137 self.Qsg.set_attrs("climate", "subglacial water flux", "m2 / s", "m2 / s", "", 0)
138
139 grid_spacing = 0.5 * (self.grid.dx() + self.grid.dy())
140 cross_section_area = self.depth * grid_spacing
141
142 self.Qsg.set(self.water_flux * cross_section_area / (grid_spacing * seconds_per_day))
143
144 self.geometry = create_geometry(self.grid)
145 self.geometry.ice_thickness.set(self.depth)
146 self.geometry.sea_level_elevation.set(self.depth)
147
148 self.geometry.ensure_consistency(config.get_number("geometry.ice_free_thickness_standard"))
149
150 self.inputs = PISM.FrontalMeltInputs()
151 self.inputs.geometry = self.geometry
152 self.inputs.subglacial_water_flux = self.Qsg
153
154 def runTest(self):
155 "Model DischargeRouting"
156
157 model = PISM.FrontalMeltDischargeRouting(self.grid)
158
159 model.initialize(self.theta)
160
161 model.update(self.inputs, 0, 1)
162
163 melt_rate = self.frontal_melt(self.depth, self.water_flux, self.potential_temperature)
164 # convert from m / day to m / s
165 melt_rate /= seconds_per_day
166
167 check_model(model, melt_rate)
168
169 assert model.max_timestep(0).infinite()
170
171 def tearDown(self):
172 pass
173
174 class GivenTest(TestCase):
175 def create_input(self, filename, melt_rate):
176 PISM.util.prepare_output(filename)
177
178 Fmr = PISM.IceModelVec2S(self.grid, "frontal_melt_rate", PISM.WITHOUT_GHOSTS)
179 Fmr.set_attrs("climate", "frontal melt rate", "m / s", "m / s", "", 0)
180
181 Fmr.set(melt_rate)
182
183 Fmr.write(filename)
184
185 def setUp(self):
186
187 self.frontal_melt_rate = 100.0
188
189 self.grid = create_grid()
190 self.geometry = create_geometry(self.grid)
191
192 self.filename = "given_input.nc"
193
194 self.create_input(self.filename, self.frontal_melt_rate)
195
196 config.set_string("frontal_melt.given.file", self.filename)
197
198 self.water_flux = PISM.IceModelVec2S(self.grid, "water_flux", PISM.WITHOUT_GHOSTS)
199 self.water_flux.set(0.0)
200
201 self.inputs = PISM.FrontalMeltInputs()
202 self.inputs.geometry = self.geometry
203 self.inputs.subglacial_water_flux = self.water_flux
204
205 def runTest(self):
206 "Model Given"
207
208 model = PISM.FrontalMeltGiven(self.grid)
209 model.init(self.geometry)
210
211 model.update(self.inputs, 0, 1)
212
213 assert model.max_timestep(0).infinite()
214
215 check_model(model, self.frontal_melt_rate)
216
217 def tearDown(self):
218 os.remove(self.filename)
219
220 if __name__ == "__main__":
221
222 t = DischargeRoutingTest()
223
224 t.setUp()
225 t.runTest()
226 t.tearDown()