tbed_smoother.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_smoother.py (3981B) --- 1 #!/usr/bin/env python3 2 3 """Simple testing program for Schoof (2003)-type bed smoothing and 4 roughness- parameterization schemes. Allows comparison of computed 5 theta to result from Matlab/Octave code exampletheta.m in 6 src/base/bedroughplay. Also used in PISM software (regression) test. 7 """ 8 9 10 import PISM 11 from math import sin, pi 12 13 ctx = PISM.Context() 14 config = ctx.config 15 16 17 def grid(): 18 "Create the bed smoother grid." 19 P = PISM.GridParameters(config) 20 21 P.horizontal_size_from_options() 22 P.horizontal_extent_from_options() 23 P.vertical_grid_from_options(config) 24 P.horizontal_extent_from_options() 25 P.ownership_ranges_from_options(ctx.size) 26 27 return PISM.IceGrid(ctx.ctx, P) 28 29 30 def allocate_storage(grid): 31 "Allocate the bed, the smoothed bed, the surface elevation, and theta." 32 topg = PISM.IceModelVec2S() 33 topg.create(grid, "topg", PISM.WITH_GHOSTS, 1) 34 topg.set_attrs("internal", "original topography", 35 "m", "m", "bedrock_altitude", 0) 36 37 topg_smoothed = PISM.IceModelVec2S() 38 topg_smoothed.create(grid, "topg_smoothed", PISM.WITHOUT_GHOSTS, 1) 39 topg_smoothed.set_attrs("internal", "smoothed topography", 40 "m", "m", "bedrock_altitude", 0) 41 42 usurf = PISM.IceModelVec2S() 43 usurf.create(grid, "usurf", PISM.WITH_GHOSTS, 1) 44 usurf.set_attrs("internal", "ice surface elevation", 45 "m", "m", "surface_altitude", 0) 46 47 theta = PISM.IceModelVec2S() 48 theta.create(grid, "theta", PISM.WITH_GHOSTS, 1) 49 theta.set_attrs("internal", 50 "coefficient theta in Schoof (2003) bed roughness parameterization", 51 "", "", "", 0) 52 53 return (topg, topg_smoothed, usurf, theta) 54 55 56 def set_topg(topg): 57 "Initialize the bed topography." 58 grid = topg.grid() 59 60 with PISM.vec.Access(comm=[topg]): 61 for (i, j) in grid.points(): 62 x = grid.x(i) 63 y = grid.y(j) 64 topg[i, j] = (400.0 * sin(2.0 * pi * x / 600.0e3) + 65 100.0 * sin(2.0 * pi * (x + 1.5 * y) / 40.0e3)) 66 67 68 def set_usurf(usurf): 69 "Initialize the surface elevation." 70 usurf.set(1000.0) 71 72 73 def set_config(): 74 "Set configuration parameters." 75 76 config.set_string("grid.periodicity", "none") 77 config.set_string("grid.registration", "corner") 78 79 config.set_number("grid.Mx", 81) 80 config.set_number("grid.My", 81) 81 82 config.set_number("grid.Lx", 1200e3) 83 config.set_number("grid.Ly", 1200e3) 84 85 config.set_number("stress_balance.sia.Glen_exponent", 3.0) 86 config.set_number("stress_balance.sia.bed_smoother.range", 50.0e3) 87 88 89 def smooth(topg, topg_smoothed, usurf, theta): 90 "Smooth the bed topography." 91 grid = topg.grid() 92 93 smoother = PISM.BedSmoother(grid, 1) 94 95 smoother.preprocess_bed(topg) 96 97 smoother.theta(usurf, theta) 98 99 topg_smoothed.copy_from(smoother.smoothed_bed()) 100 101 102 def run(): 103 "Run the bed smoother using synthetic geometry." 104 105 set_config() 106 107 topg, topg_smoothed, usurf, theta = allocate_storage(grid()) 108 109 set_usurf(usurf) 110 111 set_topg(topg) 112 113 smooth(topg, topg_smoothed, usurf, theta) 114 115 return topg, topg_smoothed, usurf, theta 116 117 118 def bed_smoother_test(): 119 "Compare the range of topg, topg_smoothed, and theta to stored values" 120 121 topg, topg_smoothed, usurf, theta = run() 122 123 stored_range = {} 124 stored_range["topg"] = [-500.0, 500.0] 125 stored_range["topg_smoothed"] = [-372.9924735817933, 372.9924735817933] 126 stored_range["theta"] = [0.7147300652935706, 0.9884843647808601] 127 128 computed_range = {} 129 for f in [topg, topg_smoothed, theta]: 130 R = f.range() 131 computed_range[f.get_name()] = [R.min, R.max] 132 133 for name in list(stored_range.keys()): 134 computed = computed_range[name] 135 stored = stored_range[name] 136 137 for k in range(2): 138 assert abs(computed[k] - stored[k]) < 1e-16 139 140 141 if __name__ == "__main__": 142 for field in run(): 143 field.dump("bed_smoother_%s.nc" % field.get_name())