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)