tvfnow.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 --- tvfnow.py (13454B) --- 1 #!/usr/bin/env python3 2 3 # @package vfnow 4 # \author Ed Bueler and Constantine Khroulev, University of Alaska Fairbanks, USA 5 # \brief A script for verification of numerical schemes in PISM. 6 # \details It specifies a refinement path for each of Tests ABCDEFGIJKL and runs 7 # pismv accordingly. 8 # Copyright (C) 2007--2013, 2015, 2016, 2018 Ed Bueler and Constantine Khroulev 9 ## 10 # Organizes the process of verifying PISM. It specifies standard refinement paths for each of the tests described in the user manual. It runs the tests, times them, and summarizes the numerical errors reported at the end. 11 ## 12 # Examples: 13 # - \verbatim vfnow.py \endverbatim use one processor and do three levels of refinement; this command is equivalent to \verbatim vfnow.py -n 2 -l 2 -t CGIJ \endverbatim, 14 # - \verbatim vfnow.py -n 8 -l 5 -t J --prefix=bin/ --mpido='aprun -n' \endverbatim will use \verbatim aprun -n 8 bin/pismv \endverbatim as the command and do five levels (the maximum) of refinement only on test J, 15 # - \verbatim vfnow.py -n 2 -l 3 -t CEIJGKLO \endverbatim uses two processers (cores) and runs in about an hour, 16 # - \verbatim vfnow.py -n 40 -l 5 -t ABCDEFGIJKLO \endverbatim will use forty processors to do all possible verification as managed by \c vfnow.py; don't run this unless you have a big computer and you are prepared to wait. 17 # For a list of options do \verbatim test/vfnow.py --help \endverbatim. 18 # Timing information is given in the \c vfnow.py output so performance, including parallel performance, can be assessed along with accuracy. 19 20 import sys 21 import time 22 import subprocess 23 from numpy import array 24 25 # A class describing a refinement path and command-line options 26 # for a particular PISM verification test. 27 28 29 class PISMVerificationTest: 30 31 # max number of levels that will work with 32 N = 50 33 # one-letter test name 34 name = "" 35 # test description 36 test = "" 37 # description of the refinement path 38 path = "" 39 Mx = [] 40 My = [] 41 # 31 levels in the ice 42 Mz = [31] * N 43 # no bedrock by default 44 Mbz = [1] * N 45 # extra options (such as -y, -ys, -ssafd_picard_rtol) 46 opts = "" 47 executable = "pismv" 48 49 def build_command(self, exec_prefix, level): 50 M = list(zip(self.Mx, self.My, self.Mz, self.Mbz)) 51 52 if level > len(M): 53 print("Test %s: Invalid refinement level: %d (only %d are available)" % ( 54 self.name, level, len(M))) 55 return "" 56 57 grid_options = "-Mx %d -My %d -Mz %d -Mbz %d" % M[level - 1] 58 return "%s%s -test %s %s %s" % (exec_prefix, self.executable, self.name, grid_options, self.opts) 59 60 61 def run_test(executable, name, level, extra_options="", debug=False): 62 try: 63 test = tests[name] 64 except: 65 print("Test %s not found." % name) 66 return 67 68 if level == 1: 69 print("# ++++ TEST %s: verifying with %s exact solution ++++\n# %s" % ( 70 test.name, test.test, test.path)) 71 else: 72 extra_options += " -append" 73 74 command = test.build_command(executable, level) + " " + extra_options 75 76 if debug: 77 print('# L%d\n%s' % (level, command)) 78 return 79 else: 80 print(' L%d: trying "%s"' % (level, command)) 81 82 # run PISM: 83 try: 84 lasttime = time.time() 85 (status, output) = subprocess.getstatusoutput(command) 86 elapsetime = time.time() - lasttime 87 except KeyboardInterrupt: 88 sys.exit(2) 89 if status: 90 sys.exit(status) 91 print(' finished in %7.4f seconds; reported numerical errors as follows:' % elapsetime) 92 93 # process the output: 94 position = output.find('NUMERICAL ERRORS') 95 if position >= 0: 96 report = output[position:output.find('NUM ERRORS DONE')] 97 endline = report.find('\n') 98 print(' ' + report[0:endline]) 99 report = report[endline + 1:] 100 while (len(report) > 1) and (endline > 0): 101 endline = report.find('\n') 102 if endline == -1: 103 endline = len(report) 104 print(' #' + report[0:endline]) 105 report = report[endline + 1:] 106 endline = report.find('\n') 107 if endline == -1: 108 endline = len(report) 109 print(' |' + report[0:endline]) 110 report = report[endline + 1:] 111 else: 112 print(" ERROR: can't find reported numerical error") 113 sys.exit(99) 114 115 116 def define_refinement_paths(KSPRTOL, SSARTOL): 117 # Define all the supported refinement paths: 118 tests = {} 119 # A 120 A = PISMVerificationTest() 121 A.name = "A" 122 A.test = "steady, marine margin isothermal SIA" 123 A.path = "(refine dx=53.33,40,26.67,20,13.33,km, dx=dy and Mx=My=31,41,61,81,121)" 124 A.Mx = [31, 41, 61, 81, 121] 125 A.My = A.Mx 126 A.opts = "-y 25000.0" 127 tests['A'] = A 128 # B 129 B = PISMVerificationTest() 130 B.name = "B" 131 B.test = "moving margin isothermal SIA (Halfar)" 132 B.path = "(refine dx=80,60,40,30,20,km, dx=dy and Mx=My=31,41,61,81,121)" 133 B.Mx = [31, 41, 61, 81, 121] 134 B.My = B.Mx 135 B.opts = "-ys 422.45 -y 25000.0" 136 tests['B'] = B 137 # C 138 C = PISMVerificationTest() 139 C.name = "C" 140 C.test = "non-zero accumulation moving margin isothermal SIA" 141 C.path = "(refine dx=50,33.33,25,20,16,km, dx=dy and Mx=My=41,61,81,101,121)" 142 C.Mx = [41, 61, 81, 101, 121] 143 C.My = C.Mx 144 C.opts = "-y 15208.0" 145 tests['C'] = C 146 # D 147 D = PISMVerificationTest() 148 D.name = "D" 149 D.test = "time-dependent isothermal SIA" 150 D.path = "(refine dx=50,33.33,25,20,16.67,km, dx=dy and Mx=My=41,61,81,101,121)" 151 D.Mx = [41, 61, 81, 101, 121] 152 D.My = D.Mx 153 D.opts = "-y 25000.0" 154 tests['D'] = D 155 # E 156 E = PISMVerificationTest() 157 E.name = "E" 158 E.test = "steady sliding marine margin isothermal SIA" 159 E.path = "(refine dx=53.33,40,26.67,20,13.33,km, dx=dy and Mx=My=31,41,61,81,121)" 160 E.Mx = [31, 41, 61, 81, 121] 161 E.My = E.Mx 162 E.opts = "-y 25000.0" 163 tests['E'] = E 164 # F 165 F = PISMVerificationTest() 166 F.name = "F" 167 F.test = "steady thermomechanically-coupled SIA" 168 F.path = "(refine dx=30,20,15,10,7.5,km, dx=dy, dz=66.67,44.44,33.33,22.22,16.67 m and Mx=My=Mz=61,91,121,181,241)" 169 F.Mx = [61, 91, 121, 181, 241] 170 F.My = F.Mx 171 F.Mz = F.Mx 172 F.opts = "-y 25000.0" 173 tests['F'] = F 174 # G 175 G = PISMVerificationTest() 176 G.name = "G" 177 G.test = "time-dependent thermomechanically-coupled SIA" 178 G.path = "(refine dx=30,20,15,10,7.5,km, dx=dy, dz=66.67,44.44,33.33,22.22,16.67 m and Mx=My=Mz=61,91,121,181,241)" 179 G.Mx = [61, 91, 121, 181, 241] 180 G.My = G.Mx 181 G.Mz = G.Mx 182 G.opts = "-y 25000.0" 183 tests['G'] = G 184 # H 185 H = PISMVerificationTest() 186 H.name = "H" 187 H.test = "moving margin, isostatic bed, isothermal SIA" 188 H.path = "(refine dx=80,60,40,30,20,km, dx=dy and Mx=My=31,41,61,81,121)" 189 H.Mx = [31, 41, 61, 81, 121] 190 H.My = H.Mx 191 H.opts = "-bed_def iso -y 60000.0" 192 tests['H'] = H 193 # I 194 I = PISMVerificationTest() 195 I.executable = "ssa_testi" 196 I.name = "I" 197 I.test = "plastic till ice stream (SSA)" 198 I.path = "(refine dy=5000,1250,312.5,78.13,19.53,m, My=49,193,769,3073,12289)" 199 I.Mx = [5] * 5 200 I.My = [49, 193, 769, 3073, 12289] 201 I.executable = "ssa_testi" 202 I.opts = "-ssa_method fd -ssafd_picard_rtol %1.e -ssafd_ksp_rtol %1.e" % (SSARTOL, KSPRTOL) 203 tests['I'] = I 204 # J 205 J = PISMVerificationTest() 206 J.executable = "ssa_testj" 207 J.name = "J" 208 J.test = "periodic ice shelf (linearized SSA)" 209 J.Mx = [49, 98, 196, 392, 784] 210 J.My = J.Mx 211 J.path = "(refine Mx={})".format(J.Mx) 212 J.Mz = [11] * 5 213 J.executable = "ssa_testj" 214 J.opts = "-ssa_method fd -ssafd_pc_type asm -ssafd_sub_pc_type lu -ssafd_ksp_rtol %1.e" % KSPRTOL 215 tests['J'] = J 216 # K 217 K = PISMVerificationTest() 218 K.name = "K" 219 K.test = "pure conduction problem in ice and bedrock" 220 K.path = "(refine dz=100,50,25,12.5,6.25,m, Mz=41,81,161,321,641)" 221 K.Mx = [8] * 5 222 K.My = K.Mx 223 K.Mz = array([41, 81, 161, 321, 641]) 224 K.Mbz = (K.Mz - 1) / 4 + 1 225 K.opts = "-y 130000.0 -Lbz 1000 -z_spacing equal" 226 tests['K'] = K 227 # L 228 L = PISMVerificationTest() 229 L.name = "L" 230 L.test = "non-flat bed stead isothermal SIA" 231 L.path = "(refine dx=60,30,20,15,10,km, dx=dy and Mx=My=31,61,91,121,181)" 232 L.Mx = [31, 61, 91, 121, 181] 233 L.My = L.Mx 234 L.opts = "-y 25000.0" 235 tests['L'] = L 236 # M 237 M = PISMVerificationTest() 238 M.name = "M" 239 M.test = "annular ice shelf with a calving front (SSA)" 240 M.path = "(refine dx=50,25,16.666,12.5,8.333 km; dx=dy and My=31,61,91,121,181)" 241 M.Mx = [31, 61, 91, 121, 181] 242 M.My = M.Mx 243 M.Mz = [11] * 5 244 M.opts = "-ssafd_picard_rtol %1.e -ssafd_ksp_rtol %1.e" % (SSARTOL, KSPRTOL) 245 tests['M'] = M 246 # O 247 O = PISMVerificationTest() 248 O.name = "O" 249 O.test = "basal melt rate from conduction problem in ice and bedrock" 250 O.path = "(refine dz=100,50,25,12.5,6.25,m, Mz=41,81,161,321,641)" 251 O.Mx = [8] * 5 252 O.My = O.Mx 253 O.Mz = array([41, 81, 161, 321, 641]) 254 O.Mbz = (O.Mz - 1) / 4 + 1 255 O.opts = "-z_spacing equal -zb_spacing equal -Lbz 1000 -y 1000 -no_mass" 256 tests['O'] = O 257 258 # test K (for a figure in the User's Manual) 259 K = PISMVerificationTest() 260 K.name = "K" 261 K.test = "pure conduction problem in ice and bedrock" 262 K.path = "(lots of levels)" 263 K.Mz = array([101, 121, 141, 161, 181, 201, 221, 241, 261, 281, 301, 321]) 264 K.Mbz = (K.Mz - 1) / 4 + 1 265 K.Mx = [8] * len(K.Mz) 266 K.My = K.Mx 267 K.opts = "-y 130000.0 -Lbz 1000" 268 tests['K_manual'] = K 269 270 # test B (for a figure in the User's Manual) 271 B = PISMVerificationTest() 272 B.name = "B" 273 B.test = "moving margin isothermal SIA (Halfar)" 274 B.path = "(lots of levels)" 275 B.Mx = [31, 41, 51, 61, 71, 81, 91, 101, 111, 121] 276 B.My = B.Mx 277 B.Mz = [31] * len(B.Mx) 278 B.Mbz = [1] * len(B.Mx) 279 B.opts = "-ys 422.45 -y 25000.0" 280 tests['B_manual'] = B 281 282 # test G (for a figure in the User's Manual) 283 G = PISMVerificationTest() 284 G.name = "G" 285 G.test = "time-dependent thermomechanically-coupled SIA" 286 G.path = "(lots of levels)" 287 G.Mx = [61, 71, 81, 91, 101, 111, 121, 151, 181] 288 G.My = G.Mx 289 G.Mz = G.Mx 290 G.opts = "-y 25000.0" 291 tests['G_manual'] = G 292 293 # test I (for a figure in the User's Manual) 294 I = PISMVerificationTest() 295 I.executable = "ssa_testi" 296 I.name = "I" 297 I.test = "plastic till ice stream (SSA)" 298 I.path = "(lots of levels)" 299 I.My = [51, 101, 151, 201, 401, 601, 801, 1001, 1501, 2001, 2501, 3073] 300 I.Mx = [5] * len(I.My) 301 I.opts = "-ssa_method fd -ssafd_picard_rtol %1.e -ssafd_ksp_rtol %1.e" % (SSARTOL, KSPRTOL) 302 tests['I_manual'] = I 303 304 return tests 305 306 307 from argparse import ArgumentParser 308 309 parser = ArgumentParser() 310 parser.description = """PISM verification script""" 311 parser.add_argument("--eta", dest="eta", action="store_true", 312 help="to add '-eta' option to pismv call") 313 parser.add_argument("-l", dest="levels", type=int, default=2, 314 help="number of levels of verification; '-l 1' fast, '-l 5' slowest") 315 parser.add_argument("--mpido", dest="mpido", default="mpiexec -np", 316 help="specify MPI executable (e.g. 'mpirun -np' or 'aprun -n')") 317 parser.add_argument("-n", dest="n", type=int, default=2, 318 help="number of processors to use") 319 parser.add_argument("--prefix", dest="prefix", default="", 320 help="path prefix to pismv executable") 321 parser.add_argument("-r", dest="report_file", default="", 322 help="name of the NetCDF error report file") 323 parser.add_argument("-t", dest="tests", nargs="+", 324 help="verification tests to use (A,B,C,D,E,F,G,H,I,J,K,L,M,O); specify a space-separated list", default=['C', 'G', 'I', 'J']) 325 parser.add_argument("-u", dest="unequal", action="store_true", 326 help="use quadratic vertical grid spacing") 327 parser.add_argument("--debug", dest="debug", action="store_true", 328 help="just print commands in sequence (do not run pismv)") 329 parser.add_argument("--manual", dest="manual", action="store_true", 330 help="run tests necessary to produce figures in the User's Manual") 331 332 options = parser.parse_args() 333 extra_options = "" 334 335 if options.eta: 336 extra_options += " -eta" 337 338 if options.unequal: 339 extra_options += " -z_spacing quadratic" 340 341 if options.report_file: 342 extra_options += " -report_file %s" % options.report_file 343 344 predo = "" 345 if options.n > 1: 346 predo = "%s %d " % (options.mpido, options.n) 347 348 exec_prefix = predo + options.prefix 349 350 KSPRTOL = 1e-12 # for tests I, J, M 351 SSARTOL = 5e-7 # ditto 352 tests = define_refinement_paths(KSPRTOL, SSARTOL) 353 354 manual_tests = ["B_manual", "G_manual", "K_manual", "I_manual"] 355 if options.manual: 356 print("# VFNOW.PY: test(s) %s, using '%s...'\n" % (manual_tests, exec_prefix) + 357 "# and ignoring options -t and -l") 358 for test in manual_tests: 359 N = len(tests[test].Mx) 360 for j in range(1, N + 1): 361 run_test(exec_prefix, test, j, extra_options, 362 options.debug) 363 else: 364 print("# VFNOW.PY: test(s) %s, %d refinement level(s), using '%s...'" % ( 365 options.tests, options.levels, exec_prefix)) 366 367 for test in options.tests: 368 for j in range(1, options.levels + 1): 369 run_test(exec_prefix, test, j, extra_options, 370 options.debug)