pism

[fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
git clone git://src.adamsgaard.dk/pism # fast
git clone https://src.adamsgaard.dk/pism.git # slow
Log | Files | Refs | README | LICENSE Back to index

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))