tbeddef_lc_restart.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_lc_restart.py (3473B)
       ---
            1 #!/usr/bin/env python3
            2 
            3 """ Sets up and runs the Lingle-Clark bed deformation model.
            4 
            5 Compares results of
            6 
            7 - Two 1000 year steps
            8 - One 1000 year step, followed by
            9   * saving the model state and
           10   * re-initializing
           11   * one more 1000 year step
           12 
           13 Used as a regression test for PISM.LingleClark.
           14 """
           15 
           16 import PISM
           17 from PISM.util import convert
           18 import numpy as np
           19 import os
           20 
           21 ctx = PISM.Context()
           22 
           23 # silence initialization messages
           24 ctx.log.set_threshold(1)
           25 
           26 # disc load parameters
           27 disc_radius = convert(1000, "km", "m")
           28 disc_thickness = 1000.0         # meters
           29 # domain size
           30 Lx = 2 * disc_radius
           31 Ly = Lx
           32 N = 101
           33 
           34 ctx.config.set_number("bed_deformation.lc.grid_size_factor", 2)
           35 
           36 dt = convert(1000.0, "years", "seconds")
           37 
           38 def add_disc_load(ice_thickness, radius, thickness):
           39     "Add a disc load with a given radius and thickness."
           40     grid = ice_thickness.grid()
           41 
           42     with PISM.vec.Access(nocomm=ice_thickness):
           43         for (i, j) in grid.points():
           44             r = PISM.radius(grid, i, j)
           45             if r <= disc_radius:
           46                 ice_thickness[i, j] = disc_thickness
           47 
           48 
           49 def run(dt, restart=False):
           50     "Run the model for 1 time step, stop, save model state, restart, do 1 more step."
           51 
           52     grid = PISM.IceGrid.Shallow(ctx.ctx, Lx, Ly, 0, 0, N, N,
           53                                 PISM.CELL_CORNER, PISM.NOT_PERIODIC)
           54 
           55     model = PISM.LingleClark(grid)
           56 
           57     geometry = PISM.Geometry(grid)
           58 
           59     bed_uplift = PISM.IceModelVec2S(grid, "uplift", PISM.WITHOUT_GHOSTS)
           60 
           61     # start with a flat bed, no ice, and no uplift
           62     geometry.bed_elevation.set(0.0)
           63     geometry.ice_thickness.set(0.0)
           64     geometry.sea_level_elevation.set(0.0)
           65     geometry.ensure_consistency(0.0)
           66 
           67     bed_uplift.set(0.0)
           68 
           69     # initialize the model
           70     model.bootstrap(geometry.bed_elevation, bed_uplift, geometry.ice_thickness, geometry.sea_level_elevation)
           71 
           72     # add the disc load
           73     add_disc_load(geometry.ice_thickness, disc_radius, disc_thickness)
           74 
           75     # do 1 step
           76     model.step(geometry.ice_thickness, geometry.sea_level_elevation, dt)
           77 
           78     if restart:
           79         # save the model state
           80         filename = "lingle_clark_model_state.nc"
           81         try:
           82             PISM.util.prepare_output(filename)
           83             f = PISM.File(grid.com, filename, PISM.PISM_NETCDF3, PISM.PISM_READWRITE)
           84             model.write_model_state(f)
           85             f.close()
           86 
           87             # create a new model
           88             del model
           89             model = PISM.LingleClark(grid)
           90 
           91             # initialize
           92             ctx.config.set_string("input.file", filename)
           93             options = PISM.process_input_options(grid.com, ctx.config)
           94             model.init(options, geometry.ice_thickness, geometry.sea_level_elevation)
           95         finally:
           96             os.remove(filename)
           97 
           98     # do 1 more step
           99     model.step(geometry.ice_thickness, geometry.sea_level_elevation, dt)
          100 
          101     return model
          102 
          103 def compare_vec(v1, v2):
          104     "Compare two vecs."
          105     print("Comparing {}".format(v1.get_name()))
          106     np.testing.assert_equal(v1.numpy(), v2.numpy())
          107 
          108 def compare(model1, model2):
          109     "Compare two models"
          110     compare_vec(model1.bed_elevation(), model2.bed_elevation())
          111     compare_vec(model1.uplift(), model2.uplift())
          112     compare_vec(model1.total_displacement(), model2.total_displacement())
          113     compare_vec(model1.viscous_displacement(), model2.viscous_displacement())
          114     compare_vec(model1.relief(), model2.relief())
          115 
          116 
          117 def lingle_clark_restart_test():
          118     "Compare straight and re-started runs."
          119     compare(run(dt),
          120             run(dt, restart=True))