tbeddef_iso.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
       ---
       tbeddef_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))