beddef_iso.py (1707B)
1 #!/usr/bin/env python3 2 3 "Test the pointwise isostasy model." 4 5 import numpy 6 import PISM 7 8 ctx = PISM.Context().ctx 9 10 def exact(ice_thickness_change): 11 "Exact deflection corresponding to a given thickness change." 12 13 config = PISM.Context().config 14 15 ice_density = config.get_number("constants.ice.density") 16 mantle_density = config.get_number("bed_deformation.mantle_density") 17 18 return -(ice_density / mantle_density) * ice_thickness_change 19 20 def bed_def_iso(ice_thickness_change): 21 "Use the pointwise isostasy model to compute plate deflection." 22 23 # grid size and domain size are irrelevant 24 M = 3 25 L = 1000e3 26 27 grid = PISM.IceGrid.Shallow(ctx, L, L, 0, 0, M, M, PISM.CELL_CORNER, PISM.NOT_PERIODIC) 28 29 geometry = PISM.Geometry(grid) 30 geometry.bed_elevation.set(0.0) 31 geometry.sea_level_elevation.set(0.0) 32 geometry.ice_thickness.set(0.0) 33 geometry.ice_area_specific_volume.set(0.0) 34 geometry.ensure_consistency(0.0) 35 36 # uplift is required (but not used) 37 bed_uplift = PISM.IceModelVec2S(grid, "uplift", PISM.WITHOUT_GHOSTS) 38 bed_uplift.set(0.0) 39 40 bed_model = PISM.PointwiseIsostasy(grid) 41 bed_model.bootstrap(geometry.bed_elevation, 42 bed_uplift, 43 geometry.ice_thickness, 44 geometry.sea_level_elevation) 45 46 geometry.ice_thickness.set(ice_thickness_change) 47 48 # time step duration is irrelevant 49 bed_model.update(geometry.ice_thickness, geometry.sea_level_elevation, 0, 1) 50 51 return bed_model.bed_elevation().numpy()[0,0] 52 53 def beddef_iso_test(): 54 "Test the pointwise isostasy model" 55 56 dH = 1000.0 57 58 numpy.testing.assert_almost_equal(bed_def_iso(dH), exact(dH))