tocean_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
---
tocean_models.py (14481B)
---
1 #!/usr/bin/env python3
2 """
3 Tests of PISM's ocean models and modifiers.
4 """
5
6 import PISM
7 from PISM.testing import *
8 import os
9 import numpy as np
10 from unittest import TestCase
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", 5)
17 config.set_number("grid.Mz", 5)
18
19 config.set_string("ocean.delta_sl_2d.file", "ocean_delta_SL_input.nc")
20
21 seconds_per_year = 365 * 86400
22 # ensure that this is the correct year length
23 config.set_string("time.calendar", "365_day")
24
25 # silence models' initialization messages
26 PISM.Context().log.set_threshold(1)
27
28 def create_given_input_file(filename, grid, temperature, mass_flux):
29 PISM.util.prepare_output(filename)
30
31 T = PISM.IceModelVec2S(grid, "shelfbtemp", PISM.WITHOUT_GHOSTS)
32 T.set_attrs("climate", "shelf base temperature", "Kelvin", "Kelvin", "", 0)
33 T.set(temperature)
34 T.write(filename)
35
36 M = PISM.IceModelVec2S(grid, "shelfbmassflux", PISM.WITHOUT_GHOSTS)
37 M.set_attrs("climate", "shelf base mass flux", "kg m-2 s-1", "kg m-2 s-1", "", 0)
38 M.set(mass_flux)
39 M.write(filename)
40
41 def check_model(model, T, SMB, MBP):
42 "Check values returned by an ocean model"
43 check(model.shelf_base_temperature(), T)
44 check(model.shelf_base_mass_flux(), SMB)
45 check(model.melange_back_pressure_fraction(), MBP)
46
47 def check_modifier(model, modifier, dT, dSMB, dMBP):
48 check_difference(modifier.shelf_base_temperature(),
49 model.shelf_base_temperature(),
50 dT)
51
52 check_difference(modifier.shelf_base_mass_flux(),
53 model.shelf_base_mass_flux(),
54 dSMB)
55
56 check_difference(modifier.melange_back_pressure_fraction(),
57 model.melange_back_pressure_fraction(),
58 dMBP)
59
60 def constant_test():
61 "Model Constant"
62
63 depth = 1000.0 # meters
64
65 # compute mass flux
66 melt_rate = config.get_number("ocean.constant.melt_rate", "m second-1")
67 ice_density = config.get_number("constants.ice.density")
68 mass_flux = melt_rate * ice_density
69
70 # compute pressure melting temperature
71 T0 = config.get_number("constants.fresh_water.melting_point_temperature")
72 beta_CC = config.get_number("constants.ice.beta_Clausius_Clapeyron")
73 g = config.get_number("constants.standard_gravity")
74
75 pressure = ice_density * g * depth
76 T_melting = T0 - beta_CC * pressure
77
78 melange_back_pressure = 0.0
79
80 grid = shallow_grid()
81 geometry = PISM.Geometry(grid)
82 geometry.ice_thickness.set(depth)
83
84 model = PISM.OceanConstant(grid)
85
86 model.init(geometry)
87 model.update(geometry, 0, 1)
88
89 check_model(model, T_melting, mass_flux, melange_back_pressure)
90
91 assert model.max_timestep(0).infinite() == True
92
93
94 def pik_test():
95 "Model PIK"
96 grid = shallow_grid()
97 geometry = PISM.Geometry(grid)
98
99 depth = 1000.0 # meters
100
101 # compute pressure melting temperature
102 ice_density = config.get_number("constants.ice.density")
103 T0 = config.get_number("constants.fresh_water.melting_point_temperature")
104 beta_CC = config.get_number("constants.ice.beta_Clausius_Clapeyron")
105 g = config.get_number("constants.standard_gravity")
106
107 pressure = ice_density * g * depth
108 T_melting = T0 - beta_CC * pressure
109
110 melange_back_pressure = 0.0
111
112 mass_flux = 5.36591610659e-06 # stored mass flux value returned by the model
113
114 # create the model
115 geometry.ice_thickness.set(depth)
116
117 model = PISM.OceanPIK(grid)
118
119 model.init(geometry)
120 model.update(geometry, 0, 1)
121
122 check_model(model, T_melting, mass_flux, melange_back_pressure)
123
124 assert model.max_timestep(0).infinite() == True
125
126 class GivenTest(TestCase):
127 "Test the Given class"
128
129 def setUp(self):
130 self.grid = shallow_grid()
131 self.geometry = PISM.Geometry(self.grid)
132 self.filename = "ocean_given_input.nc"
133
134 self.temperature = 263.0
135 self.mass_flux = 3e-3
136 self.melange_back_pressure = 0.0
137
138 create_given_input_file(self.filename, self.grid, self.temperature, self.mass_flux)
139
140 config.set_string("ocean.given.file", self.filename)
141
142 def test_ocean_given(self):
143 "Model Given"
144
145 model = PISM.OceanGiven(self.grid)
146 model.init(self.geometry)
147 model.update(self.geometry, 0, 1)
148
149 assert model.max_timestep(0).infinite() == True
150
151 check_model(model, self.temperature, self.mass_flux, self.melange_back_pressure)
152
153 def tearDown(self):
154 os.remove(self.filename)
155
156 class GivenTHTest(TestCase):
157 def setUp(self):
158
159 depth = 1000.0
160 salinity = 35.0
161 potential_temperature = 270.0
162 self.melange_back_pressure = 0.0
163 self.temperature = 270.17909999999995
164 self.mass_flux = -6.489250000000001e-05
165
166 self.grid = shallow_grid()
167 self.geometry = PISM.Geometry(self.grid)
168
169 self.geometry.ice_thickness.set(depth)
170
171 filename = "ocean_given_th_input.nc"
172 self.filename = filename
173
174 PISM.util.prepare_output(filename)
175
176 Th = PISM.IceModelVec2S(self.grid, "theta_ocean", PISM.WITHOUT_GHOSTS)
177 Th.set_attrs("climate", "potential temperature", "Kelvin", "Kelvin", "", 0)
178 Th.set(potential_temperature)
179 Th.write(filename)
180
181 S = PISM.IceModelVec2S(self.grid, "salinity_ocean", PISM.WITHOUT_GHOSTS)
182 S.set_attrs("climate", "ocean salinity", "g/kg", "g/kg", "", 0)
183 S.set(salinity)
184 S.write(filename)
185
186 config.set_string("ocean.th.file", self.filename)
187
188 def test_ocean_th(self):
189 "Model GivenTH"
190
191 model = PISM.OceanGivenTH(self.grid)
192 model.init(self.geometry)
193 model.update(self.geometry, 0, 1)
194
195 assert model.max_timestep(0).infinite() == True
196
197 check_model(model, self.temperature, self.mass_flux, self.melange_back_pressure)
198
199 def tearDown(self):
200 os.remove(self.filename)
201
202 class DeltaT(TestCase):
203 def setUp(self):
204 self.filename = "ocean_delta_T_input.nc"
205 self.grid = shallow_grid()
206 self.geometry = PISM.Geometry(self.grid)
207 self.geometry.ice_thickness.set(1000.0)
208 self.model = PISM.OceanConstant(self.grid)
209 self.dT = -5.0
210
211 PISM.testing.create_scalar_forcing(self.filename, "delta_T", "Kelvin", [self.dT], [0])
212
213 def test_ocean_delta_t(self):
214 "Modifier Delta_T"
215
216 modifier = PISM.OceanDeltaT(self.grid, self.model)
217
218 config.set_string("ocean.delta_T.file", self.filename)
219
220 modifier.init(self.geometry)
221 modifier.update(self.geometry, 0, 1)
222
223 check_modifier(self.model, modifier, self.dT, 0.0, 0.0)
224
225 def tearDown(self):
226 os.remove(self.filename)
227
228
229 class DeltaSMB(TestCase):
230 def setUp(self):
231 self.filename = "ocean_delta_SMB_input.nc"
232 self.grid = shallow_grid()
233 self.geometry = PISM.Geometry(self.grid)
234 self.model = PISM.OceanConstant(self.grid)
235 self.dSMB = -5.0
236
237 create_scalar_forcing(self.filename, "delta_mass_flux", "kg m-2 s-1", [self.dSMB], [0])
238
239 def test_ocean_delta_smb(self):
240 "Modifier Delta_SMB"
241
242 modifier = PISM.OceanDeltaSMB(self.grid, self.model)
243
244 config.set_string("ocean.delta_mass_flux.file", self.filename)
245
246 modifier.init(self.geometry)
247 modifier.update(self.geometry, 0, 1)
248
249 check_modifier(self.model, modifier, 0.0, self.dSMB, 0.0)
250
251 def tearDown(self):
252 os.remove(self.filename)
253
254 class AnomalyBMB(TestCase):
255 def setUp(self):
256 self.filename = "ocean_delta_BMB_input.nc"
257 self.grid = shallow_grid()
258 self.geometry = PISM.Geometry(self.grid)
259 self.model = PISM.OceanConstant(self.grid)
260 self.dBMB = -5.0
261
262 delta_BMB = PISM.IceModelVec2S(self.grid, "shelf_base_mass_flux_anomaly",
263 PISM.WITHOUT_GHOSTS)
264 delta_BMB.set_attrs("climate_forcing",
265 "2D shelf base mass flux anomaly", "kg m-2 s-1", "kg m-2 s-1", "", 0)
266 delta_BMB.set(self.dBMB)
267
268 delta_BMB.dump(self.filename)
269
270 def test_ocean_anomaly(self):
271 "Modifier Anomaly"
272
273 config.set_string("ocean.anomaly.file", self.filename)
274
275 modifier = PISM.OceanAnomaly(self.grid, self.model)
276
277 modifier.init(self.geometry)
278 modifier.update(self.geometry, 0, 1)
279
280 check_modifier(self.model, modifier, 0.0, self.dBMB, 0.0)
281
282 def tearDown(self):
283 os.remove(self.filename)
284
285 class FracMBP(TestCase):
286 def setUp(self):
287 self.filename = "ocean_frac_MBP_input.nc"
288 self.grid = shallow_grid()
289 self.geometry = PISM.Geometry(self.grid)
290 self.model = PISM.OceanConstant(self.grid)
291 self.dMBP = 0.5
292
293 create_scalar_forcing(self.filename, "frac_MBP", "1", [self.dMBP], [0])
294 config.set_number("ocean.melange_back_pressure_fraction", 1.0)
295
296 def test_ocean_frac_mpb(self):
297 "Modifier Frac_MBP"
298
299 modifier = PISM.OceanFracMBP(self.grid, self.model)
300
301 config.set_string("ocean.frac_MBP.file", self.filename)
302
303 modifier.init(self.geometry)
304 modifier.update(self.geometry, 0, 1)
305
306 model = self.model
307
308 check_difference(modifier.shelf_base_temperature(),
309 model.shelf_base_temperature(),
310 0.0)
311
312 check_difference(modifier.shelf_base_mass_flux(),
313 model.shelf_base_mass_flux(),
314 0.0)
315
316 check_ratio(modifier.melange_back_pressure_fraction(),
317 model.melange_back_pressure_fraction(),
318 self.dMBP)
319
320 def tearDown(self):
321 os.remove(self.filename)
322 config.set_number("ocean.melange_back_pressure_fraction", 0.0)
323
324
325 class FracSMB(TestCase):
326 def setUp(self):
327 self.filename = "ocean_frac_SMB_input.nc"
328 self.grid = shallow_grid()
329 self.geometry = PISM.Geometry(self.grid)
330 self.model = PISM.OceanConstant(self.grid)
331 self.dSMB = 0.5
332
333 create_scalar_forcing(self.filename, "frac_mass_flux", "1", [self.dSMB], [0])
334
335 def test_ocean_frac_smb(self):
336 "Modifier Frac_SMB"
337
338 modifier = PISM.OceanFracSMB(self.grid, self.model)
339
340 config.set_string("ocean.frac_mass_flux.file", self.filename)
341
342 modifier.init(self.geometry)
343 modifier.update(self.geometry, 0, 1)
344
345 model = self.model
346
347 check_difference(modifier.shelf_base_temperature(),
348 model.shelf_base_temperature(),
349 0.0)
350
351 check_ratio(modifier.shelf_base_mass_flux(),
352 model.shelf_base_mass_flux(),
353 self.dSMB)
354
355 check_difference(modifier.melange_back_pressure_fraction(),
356 model.melange_back_pressure_fraction(),
357 0.0)
358
359 def tearDown(self):
360 os.remove(self.filename)
361
362 class Cache(TestCase):
363 def setUp(self):
364 self.filename = "ocean_dT.nc"
365 self.grid = shallow_grid()
366 self.geometry = PISM.Geometry(self.grid)
367
368 self.constant = PISM.OceanConstant(self.grid)
369 self.delta_T = PISM.OceanDeltaT(self.grid, self.constant)
370
371 time_bounds = np.array([0, 1, 1, 2, 2, 3, 3, 4]) * seconds_per_year
372 create_scalar_forcing(self.filename, "delta_T", "Kelvin", [1, 2, 3, 4],
373 times=None, time_bounds=time_bounds)
374
375 config.set_string("ocean.delta_T.file", self.filename)
376 config.set_number("ocean.cache.update_interval", 2.0)
377
378 def test_ocean_cache(self):
379 "Modifier Cache"
380
381 modifier = PISM.OceanCache(self.grid, self.delta_T)
382
383 modifier.init(self.geometry)
384
385 dt = seconds_per_year
386
387 N = 4
388 ts = np.arange(float(N)) * dt
389 diff = []
390 for t in ts:
391 modifier.update(self.geometry, t, dt)
392
393 original = sample(self.constant.shelf_base_temperature())
394 cached = sample(modifier.shelf_base_temperature())
395
396 diff.append(cached - original)
397
398 np.testing.assert_almost_equal(diff, [1, 1, 3, 3])
399
400 def tearDown(self):
401 os.remove(self.filename)
402
403 class DeltaSL(TestCase):
404 def setUp(self):
405 self.filename = "ocean_delta_SL_input.nc"
406 self.grid = shallow_grid()
407 self.geometry = PISM.Geometry(self.grid)
408 self.model = PISM.SeaLevel(self.grid)
409 self.dSL = -5.0
410
411 create_scalar_forcing(self.filename, "delta_SL", "meters", [self.dSL], [0])
412
413 def test_ocean_delta_sl(self):
414 "Modifier Delta_SL"
415
416 modifier = PISM.SeaLevelDelta(self.grid, self.model)
417
418 config.set_string("ocean.delta_sl.file", self.filename)
419
420 modifier.init(self.geometry)
421 modifier.update(self.geometry, 0, 1)
422
423 check_difference(modifier.elevation(),
424 self.model.elevation(),
425 self.dSL)
426
427 def tearDown(self):
428 os.remove(self.filename)
429
430 class DeltaSL2D(TestCase):
431 def create_delta_SL_file(self, filename, times, sea_level_offsets):
432 output = PISM.util.prepare_output(filename, append_time=False)
433
434 SL = PISM.IceModelVec2S(self.grid, "delta_SL", PISM.WITHOUT_GHOSTS)
435 SL.set_attrs("forcing", "sea level forcing", "meters", "meters", "", 0)
436
437 for k in range(len(times)):
438 PISM.append_time(output, config, times[k])
439 SL.set(sea_level_offsets[k])
440 SL.write(output)
441
442 def setUp(self):
443 self.filename = "ocean_delta_SL_input.nc"
444 self.grid = shallow_grid()
445 self.geometry = PISM.Geometry(self.grid)
446 self.model = PISM.SeaLevel(self.grid)
447 self.dSL = -5.0
448
449 self.create_delta_SL_file(self.filename, [0, seconds_per_year], [0, 2 * self.dSL])
450
451 config.set_string("ocean.delta_sl_2d.file", self.filename)
452
453 def test_ocean_delta_sl_2d(self):
454 "Modifier Delta_SL_2D"
455
456 modifier = PISM.SeaLevelDelta2D(self.grid, self.model)
457
458 modifier.init(self.geometry)
459 # Use a one second time step to try to sample sea level forcing midway through the
460 # interval from 0 to 1 year.
461 modifier.update(self.geometry, 0.5 * seconds_per_year, 1)
462
463 check_difference(modifier.elevation(),
464 self.model.elevation(),
465 self.dSL)
466
467 def tearDown(self):
468 os.remove(self.filename)