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)