tbed_deformation.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 --- tbed_deformation.py (5232B) --- 1 #!/usr/bin/env python3 2 3 import PISM 4 from PISM.util import convert 5 from math import cos, pi 6 7 # Simple testing program for Lingle & Clark bed deformation model. 8 # Runs go for 150,000 years on 63.5km grid with 100a time steps and Z=2 in L&C model. 9 # SCENARIOS: run 'python bed_deformation.py -scenario N' where N=1,2,3,4 as follows 10 # (1) dump ice disc on initially level, non-uplifting land, use only viscous 11 # half-space model: 12 # include_elastic = FALSE, do_uplift = FALSE, H0 = 1000.0 13 # center depth b(0,0) should eventually equilibriate to near 14 # -1000 * (910/3300) = -275.76 m 15 # (2) dump ice disc on initially level, non-uplifting land, use both viscous 16 # half-space model and elastic model 17 # include_elastic = TRUE, do_uplift = FALSE, H0 = 1000.0 18 # (3) never loaded, initially level, uplifting land, use only viscous 19 # half-space model (because elastic model gives no additional when no load): 20 # include_elastic = FALSE, do_uplift = TRUE, H0 = 0.0 21 # (4) dump ice disc on initially level, uplifting land, use both viscous 22 # half-space model and elastic model: 23 # include_elastic = TRUE, do_uplift = TRUE, H0 = 1000.0; 24 25 ctx = PISM.Context() 26 config = ctx.config 27 28 R0 = 1000e3 29 30 31 def initialize_uplift(uplift): 32 "Initialize the uplift field." 33 grid = uplift.grid() 34 peak_uplift = convert(10, "mm/year", "m/second") 35 with PISM.vec.Access(nocomm=[uplift]): 36 for (i, j) in grid.points(): 37 r = PISM.radius(grid, i, j) 38 if r < 1.5 * R0: 39 uplift[i, j] = peak_uplift * (cos(pi * (r / (1.5 * R0))) + 1.0) / 2.0 40 else: 41 uplift[i, j] = 0.0 42 43 44 def initialize_thickness(thickness, H0): 45 grid = thickness.grid() 46 with PISM.vec.Access(nocomm=[thickness]): 47 for (i, j) in grid.points(): 48 r = PISM.radius(grid, i, j) 49 if r < R0: 50 thickness[i, j] = H0 51 else: 52 thickness[i, j] = 0.0 53 54 55 def allocate(grid): 56 H = PISM.model.createIceThicknessVec(grid) 57 bed = PISM.model.createBedrockElevationVec(grid) 58 uplift = PISM.IceModelVec2S() 59 uplift.create(grid, "uplift", PISM.WITHOUT_GHOSTS) 60 uplift.set_attrs("internal", "bed uplift", "m / second", "m / second", "", 0) 61 62 sea_level = PISM.IceModelVec2S(grid, "sea_level", PISM.WITHOUT_GHOSTS) 63 64 return H, bed, uplift, sea_level 65 66 67 def create_grid(): 68 P = PISM.GridParameters(config) 69 P.horizontal_size_from_options() 70 P.horizontal_extent_from_options() 71 P.vertical_grid_from_options(config) 72 P.ownership_ranges_from_options(ctx.size) 73 74 return PISM.IceGrid(ctx.ctx, P) 75 76 77 def run(scenario, plot, pause, save): 78 79 # set grid defaults 80 config.set_number("grid.Mx", 193) 81 config.set_number("grid.My", 129) 82 83 config.set_number("grid.Lx", 3000e3) 84 config.set_number("grid.Ly", 2000e3) 85 86 config.set_number("grid.Mz", 2) 87 config.set_number("grid.Lz", 1000) 88 89 scenarios = {"1": (False, False, 1000.0), 90 "2": (True, False, 1000.0), 91 "3": (False, True, 0.0), 92 "4": (True, True, 1000.0)} 93 94 elastic, use_uplift, H0 = scenarios[scenario] 95 96 print("Using scenario %s: elastic model = %s, use uplift = %s, H0 = %f m" % (scenario, elastic, use_uplift, H0)) 97 98 config.set_flag("bed_deformation.lc.elastic_model", elastic) 99 100 grid = create_grid() 101 102 thickness, bed, uplift, sea_level = allocate(grid) 103 104 # set initial geometry and uplift 105 bed.set(0.0) 106 thickness.set(0.0) 107 sea_level.set(0.0) 108 if use_uplift: 109 initialize_uplift(uplift) 110 111 time = ctx.ctx.time() 112 113 time.init(ctx.ctx.log()) 114 115 model = PISM.LingleClark(grid) 116 117 model.bootstrap(bed, uplift, thickness, sea_level) 118 119 # now add the disc load 120 initialize_thickness(thickness, H0) 121 122 dt = convert(100, "365 day", "seconds") 123 124 # the time-stepping loop 125 while time.current() < time.end(): 126 # don't go past the end of the run 127 dt_current = min(dt, time.end() - time.current()) 128 129 model.update(thickness, sea_level, time.current(), dt_current) 130 131 if plot: 132 model.bed_elevation().view(400) 133 model.uplift().view(400) 134 135 print("t = %s years, dt = %s years" % (time.date(), time.convert_time_interval(dt_current, "years"))) 136 time.step(dt_current) 137 138 print("Reached t = %s years" % time.date()) 139 140 if pause: 141 print("Pausing for 5 seconds...") 142 PISM.PETSc.Sys.sleep(5) 143 144 if save: 145 model.bed_elevation().dump("bed_elevation.nc") 146 model.uplift().dump("bed_uplift.nc") 147 148 149 if __name__ == "__main__": 150 scenario = PISM.OptionKeyword("-scenario", "choose one of 4 scenarios", "1,2,3,4", "1") 151 plot = PISM.OptionBool("-plot", "Plot bed elevation and uplift.") 152 save = PISM.OptionBool("-save", "Save final states of the bed elevation and uplift.") 153 pause = PISM.OptionBool("-pause", "Pause for 5 seconds to look at runtime 2D plots.") 154 155 run(scenario.value(), plot, pause, save) 156 157 158 def scenario1_test(): 159 "Test if scenario 1 runs" 160 run("1", False, False, False) 161 162 163 def scenario3_test(): 164 "Test if scenario 3 runs" 165 run("3", False, False, False)