ttest_invssaforward.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
---
ttest_invssaforward.py (18703B)
---
1 #! /usr/bin/env python3
2 #
3 # Copyright (C) 2012, 2014, 2015, 2016, 2017, 2018, 2019 David Maxwell
4 #
5 # This file is part of PISM.
6 #
7 # PISM is free software; you can redistribute it and/or modify it under the
8 # terms of the GNU General Public License as published by the Free Software
9 # Foundation; either version 3 of the License, or (at your option) any later
10 # version.
11 #
12 # PISM is distributed in the hope that it will be useful, but WITHOUT ANY
13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14 # FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
15 # details.
16 #
17 # You should have received a copy of the GNU General Public License
18 # along with PISM; if not, write to the Free Software
19 # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
20
21 # Compares the computations of the key methods of a forward problem
22 # with alternative computations (finite differences for derivatives, etc.)
23 # Takes all the arguments of pismi.py, plus -inv_test with the name of
24 # one of the tests in the file. The solver is set with its initial value
25 # of tauc equal to what it would be for pismi, and then the tests occur.
26 #
27 # Tests:
28 # lin -- the linearized forward map, compared to finite difference approximations.
29 # lin_transpose -- the transpose of this map, T^*, compared to the formula <Td,r> = <d,T^*r>
30 # j_design -- the derivative of the SSA residual with respect to the design variable, compared to finite difference approx
31 # J_design_transpose -- the transpose of this map, J^*, compared to the formula <Jd,r> = <d,J^*r>
32
33 import sys
34 import petsc4py
35 petsc4py.init(sys.argv)
36 from petsc4py import PETSc
37 import numpy as np
38 import os
39 import math
40 import PISM.invert.ssa
41 from PISM.logging import logMessage
42
43 import PISM
44
45
46 def view(vec, viewer):
47 if (isinstance(vec, PISM.IceModelVec2S)):
48 v_global = PISM.IceModelVec2S()
49 else:
50 v_global = PISM.IceModelVec2V()
51 v_global.create(vec.grid(), "", PISM.WITHOUT_GHOSTS)
52 v_global.copy_from(vec)
53 v_global.get_vec().view(viewer)
54
55 # def view2(vec,viewer):
56 # v_global = PISM.IceModelVec2V()
57 # v_global.create(vec.grid(),"",PISM.WITHOUT_GHOSTS)
58 # v_global.copy_from(vec)
59 # v_global.get_vec().view(viewer)
60
61
62 def adjustTauc(mask, tauc):
63 """Where ice is floating or land is ice-free, tauc should be adjusted to have some preset default values."""
64
65 grid = mask.grid()
66 high_tauc = grid.ctx().config().get_number("basal_yield_stress.ice_free_bedrock")
67
68 with PISM.vec.Access(comm=tauc, nocomm=mask):
69 for (i, j) in grid.points():
70 if mask.ocean(i, j):
71 tauc[i, j] = 0
72 elif mask.ice_free(i, j):
73 tauc[i, j] = high_tauc
74
75
76 def createDesignVec(grid, design_var, name=None, **kwargs):
77 if name is None:
78 name = design_var
79 if design_var == "tauc":
80 design_vec = PISM.model.createYieldStressVec(grid, name=name, **kwargs)
81 elif design_var == "hardav":
82 design_vec = PISM.model.createAveragedHardnessVec(grid, name=name, **kwargs)
83 else:
84 raise ValueError("Unknown design variable %s" % design_var)
85 return design_vec
86
87
88 WIDE_STENCIL = 2
89
90
91 def test_lin(ssarun):
92 grid = ssarun.grid
93
94 PISM.verbPrintf(1, grid.com, "\nTest Linearization (Comparison with finite differences):\n")
95
96 S = 250
97 d_viewer = PETSc.Viewer().createDraw(title="d", size=S)
98 Td_viewer = PETSc.Viewer().createDraw(title="Td", size=S)
99 Td_fd_viewer = PETSc.Viewer().createDraw(title="Td_fd", size=S)
100 d_Td_viewer = PETSc.Viewer().createDraw(title="d_Td", size=S)
101
102 for (i, j) in grid.points():
103 d = PISM.IceModelVec2S()
104 d.create(grid, "", PISM.WITH_GHOSTS)
105 d.set(0)
106 with PISM.vec.Access(comm=d):
107 d[i, j] = 1
108
109 ssarun.ssa.linearize_at(zeta1)
110 u1 = PISM.IceModelVec2V()
111 u1.create(grid, "", PISM.WITH_GHOSTS)
112 u1.copy_from(ssarun.ssa.solution())
113
114 Td = PISM.IceModelVec2V()
115 Td.create(grid, "", PISM.WITH_GHOSTS)
116 ssarun.ssa.apply_linearization(d, Td)
117
118 eps = 1e-8
119 zeta2 = PISM.IceModelVec2S()
120 zeta2.create(grid, "", PISM.WITH_GHOSTS)
121 zeta2.copy_from(d)
122 zeta2.scale(eps)
123 zeta2.add(1, zeta1)
124 ssarun.ssa.linearize_at(zeta2)
125 u2 = PISM.IceModelVec2V()
126 u2.create(grid, "", PISM.WITH_GHOSTS)
127 u2.copy_from(ssarun.ssa.solution())
128
129 Td_fd = PISM.IceModelVec2V()
130 Td_fd.create(grid, "", PISM.WITH_GHOSTS)
131 Td_fd.copy_from(u2)
132 Td_fd.add(-1, u1)
133 Td_fd.scale(1. / eps)
134
135 d_Td = PISM.IceModelVec2V()
136 d_Td.create(grid, "", PISM.WITH_GHOSTS)
137 d_Td.copy_from(Td_fd)
138 d_Td.add(-1, Td)
139
140 n_Td_fd = Td_fd.norm(PETSc.NormType.NORM_2)
141 n_Td_l2 = Td.norm(PETSc.NormType.NORM_2)
142 n_Td_l1 = Td.norm(PETSc.NormType.NORM_1)
143 n_Td_linf = Td.norm(PETSc.NormType.NORM_INFINITY)
144
145 n_d_Td_l2 = d_Td.norm(PETSc.NormType.NORM_2)
146 n_d_Td_l1 = d_Td.norm(PETSc.NormType.NORM_1)
147 n_d_Td_linf = d_Td.norm(PETSc.NormType.NORM_INFINITY)
148
149 PISM.verbPrintf(1, grid.com, "(i,j)=(%d,%d)\n" % (i, j))
150 PISM.verbPrintf(1, grid.com, "apply_linearization(d): l2 norm %.10g; finite difference %.10g\n" %
151 (n_Td_l2, n_Td_fd))
152
153 r_d_l2 = 0
154 if n_Td_l2 != 0:
155 r_d_l2 = n_d_Td_l2 / n_Td_l2
156 r_d_l1 = 0
157 if n_Td_l1 != 0:
158 r_d_l1 = n_d_Td_l1 / n_Td_l1
159
160 r_d_linf = 0
161 if n_Td_linf != 0:
162 r_d_linf = n_d_Td_linf / n_Td_linf
163 PISM.verbPrintf(1, grid.com, "relative difference: l2 norm %.10g l1 norm %.10g linf norm %.10g\n" %
164 (r_d_l2, r_d_l1, r_d_linf))
165
166 PISM.verbPrintf(1, grid.com, "\n")
167
168 d_global = PISM.IceModelVec2S()
169 d_global.create(grid, "", PISM.WITHOUT_GHOSTS)
170 d_global.copy_from(d)
171 d_global.get_vec().view(d_viewer)
172
173 Td_global = PISM.IceModelVec2V()
174 Td_global.create(grid, "", PISM.WITHOUT_GHOSTS)
175 Td_global.copy_from(Td)
176 # if n_Td_linf > 0:
177 # Td_global.scale(1./n_Td_linf)
178 Td_global.get_vec().view(Td_viewer)
179
180 Td_fd_global = PISM.IceModelVec2V()
181 Td_fd_global.create(grid, "", PISM.WITHOUT_GHOSTS)
182 Td_fd_global.copy_from(Td_fd)
183 # if n_Td_fd_linf > 0:
184 # Td_fd_global.scale(1./n_Td_linf)
185 Td_fd_global.get_vec().view(Td_fd_viewer)
186
187 d_Td_global = PISM.IceModelVec2V()
188 d_Td_global.create(grid, "", PISM.WITHOUT_GHOSTS)
189 d_Td_global.copy_from(d_Td)
190 # if n_Td_linf > 0:
191 # d_Td_global.scale(1./n_Td_linf)
192 d_Td_global.get_vec().view(d_Td_viewer)
193
194 PISM.logging.pause()
195
196 # ######################################################################################################################
197 # Jacobian design
198
199
200 def test_j_design(ssarun):
201 grid = ssarun.grid
202
203 S = 250
204 d_viewer = PETSc.Viewer().createDraw(title="d", size=S)
205 drhs_viewer = PETSc.Viewer().createDraw(title="drhs", size=S)
206 drhs_fd_viewer = PETSc.Viewer().createDraw(title="drhs_fd", size=S)
207 d_drhs_viewer = PETSc.Viewer().createDraw(title="d_drhs", size=S)
208
209 ssarun.ssa.linearize_at(zeta1)
210 u1 = PISM.IceModelVec2V()
211 u1.create(grid, "", PISM.WITH_GHOSTS)
212 u1.copy_from(ssarun.ssa.solution())
213
214 for (i, j) in grid.points():
215 d = PISM.IceModelVec2S()
216 d.create(grid, "", PISM.WITH_GHOSTS)
217 d.set(0)
218 with PISM.vec.Access(comm=d):
219 d[i, j] = 1
220
221 ssarun.ssa.linearize_at(zeta1)
222
223 rhs1 = PISM.IceModelVec2V()
224 rhs1.create(grid, "", PISM.WITHOUT_GHOSTS)
225 ssarun.ssa.assemble_residual(u1, rhs1)
226
227 eps = 1e-8
228 zeta2 = PISM.IceModelVec2S()
229 zeta2.create(grid, "zeta_prior", PISM.WITH_GHOSTS)
230 zeta2.copy_from(d)
231 zeta2.scale(eps)
232 zeta2.add(1, zeta1)
233 ssarun.ssa.set_design(zeta2)
234
235 rhs2 = PISM.IceModelVec2V()
236 rhs2.create(grid, "", PISM.WITHOUT_GHOSTS)
237 ssarun.ssa.assemble_residual(u1, rhs2)
238
239 drhs_fd = PISM.IceModelVec2V()
240 drhs_fd.create(grid, "", PISM.WITHOUT_GHOSTS)
241 drhs_fd.copy_from(rhs2)
242 drhs_fd.add(-1, rhs1)
243 drhs_fd.scale(1. / eps)
244
245 drhs = PISM.IceModelVec2V()
246 drhs.create(grid, "", PISM.WITHOUT_GHOSTS)
247 ssarun.ssa.apply_jacobian_design(u1, d, drhs)
248
249 d_drhs = PISM.IceModelVec2V()
250 d_drhs.create(grid, "", PISM.WITHOUT_GHOSTS)
251
252 d_drhs.copy_from(drhs)
253 d_drhs.add(-1, drhs_fd)
254
255 n_drhs_fd = drhs_fd.norm(PETSc.NormType.NORM_2)
256 n_drhs_l2 = drhs.norm(PETSc.NormType.NORM_2)
257 n_drhs_l1 = drhs.norm(PETSc.NormType.NORM_1)
258 n_drhs_linf = drhs.norm(PETSc.NormType.NORM_INFINITY)
259
260 n_d_drhs_l2 = d_drhs.norm(PETSc.NormType.NORM_2)
261 n_d_drhs_l1 = d_drhs.norm(PETSc.NormType.NORM_1)
262 n_d_drhs_linf = d_drhs.norm(PETSc.NormType.NORM_INFINITY)
263
264 PISM.verbPrintf(1, grid.com, "\nTest Jacobian Design (Comparison with finite differences):\n")
265 PISM.verbPrintf(1, grid.com, "jacobian_design(d): l2 norm %.10g; finite difference %.10g\n" %
266 (n_drhs_l2, n_drhs_fd))
267 if n_drhs_linf == 0:
268 PISM.verbPrintf(1, grid.com, "difference: l2 norm %.10g l1 norm %.10g linf norm %.10g\n" %
269 (n_d_drhs_l2, n_d_drhs_l1, n_d_drhs_linf))
270 else:
271 PISM.verbPrintf(1, grid.com, "relative difference: l2 norm %.10g l1 norm %.10g linf norm %.10g\n" %
272 (n_d_drhs_l2 / n_drhs_l2, n_d_drhs_l1 / n_drhs_l1, n_d_drhs_linf / n_drhs_linf))
273
274 view(d, d_viewer)
275 view(drhs, drhs_viewer)
276 view(drhs_fd, drhs_fd_viewer)
277 view(d_drhs, d_drhs_viewer)
278
279 PISM.logging.pause()
280
281
282 def test_j_design_transpose(ssarun):
283 grid = ssarun.grid
284
285 S = 250
286 r_viewer = PETSc.Viewer().createDraw(title="r", size=S)
287 JStarR_viewer = PETSc.Viewer().createDraw(title="JStarR", size=S)
288 JStarR_indirect_viewer = PETSc.Viewer().createDraw(title="JStarR (ind)", size=S)
289 d_JStarR_viewer = PETSc.Viewer().createDraw(title="d_JStarR_fd", size=S)
290
291 ssarun.ssa.linearize_at(zeta1)
292 u = PISM.IceModelVec2V()
293 u.create(grid, "", PISM.WITH_GHOSTS)
294 u.copy_from(ssarun.ssa.solution())
295
296 Jd = PISM.IceModelVec2V()
297 Jd.create(grid, "", PISM.WITHOUT_GHOSTS)
298
299 JStarR = PISM.IceModelVec2S()
300 JStarR.create(grid, "", PISM.WITHOUT_GHOSTS)
301
302 JStarR_indirect = PISM.IceModelVec2S()
303 JStarR_indirect.create(grid, "", PISM.WITHOUT_GHOSTS)
304
305 for (i, j) in grid.points():
306
307 for k in range(2):
308
309 r = PISM.IceModelVec2V()
310 r.create(grid, "", PISM.WITH_GHOSTS)
311 r.set(0)
312 with PISM.vec.Access(comm=r):
313 if k == 0:
314 r[i, j].u = 1
315 else:
316 r[i, j].v = 1
317
318 ssarun.ssa.apply_jacobian_design_transpose(u, r, JStarR)
319
320 r_global = PISM.IceModelVec2V()
321 r_global.create(grid, "", PISM.WITHOUT_GHOSTS)
322 r_global.copy_from(r)
323
324 for (k, l) in grid.points():
325 with PISM.vec.Access(nocomm=JStarR_indirect):
326 d = PISM.IceModelVec2S()
327 d.create(grid, "", PISM.WITH_GHOSTS)
328 d.set(0)
329 with PISM.vec.Access(comm=d):
330 d[k, l] = 1
331
332 ssarun.ssa.apply_jacobian_design(u, d, Jd)
333
334 JStarR_indirect[k, l] = Jd.get_vec().dot(r_global.get_vec())
335
336 d_JStarR = PISM.IceModelVec2S()
337 d_JStarR.create(grid, "", PISM.WITHOUT_GHOSTS)
338
339 d_JStarR.copy_from(JStarR)
340 d_JStarR.add(-1, JStarR_indirect)
341
342 PISM.verbPrintf(1, grid.com, "\nTest Jacobian Design Transpose (%d,%d):\n" % (i, j))
343
344 view(r_global, r_viewer)
345 view(JStarR, JStarR_viewer)
346 view(JStarR_indirect, JStarR_indirect_viewer)
347 view(d_JStarR, d_JStarR_viewer)
348
349 PISM.logging.pause()
350
351
352 def test_linearization_transpose(ssarun):
353 grid = ssarun.grid
354
355 S = 250
356 r_viewer = PETSc.Viewer().createDraw(title="r", size=S)
357 TStarR_viewer = PETSc.Viewer().createDraw(title="TStarR", size=S)
358 TStarR_indirect_viewer = PETSc.Viewer().createDraw(title="TStarR (ind)", size=S)
359 d_TStarR_viewer = PETSc.Viewer().createDraw(title="d_TStarR_fd", size=S)
360
361 ssarun.ssa.linearize_at(zeta1)
362 u = PISM.IceModelVec2V()
363 u.create(grid, "", PISM.WITH_GHOSTS)
364 u.copy_from(ssarun.ssa.solution())
365
366 Td = PISM.IceModelVec2V()
367 Td.create(grid, "", PISM.WITHOUT_GHOSTS)
368
369 TStarR = PISM.IceModelVec2S()
370 TStarR.create(grid, "", PISM.WITHOUT_GHOSTS)
371
372 TStarR_indirect = PISM.IceModelVec2S()
373 TStarR_indirect.create(grid, "", PISM.WITHOUT_GHOSTS)
374
375 for (i, j) in grid.points():
376
377 for k in range(2):
378
379 r = PISM.IceModelVec2V()
380 r.create(grid, "", PISM.WITH_GHOSTS)
381 r.set(0)
382 with PISM.vec.Access(comm=r):
383 if k == 0:
384 r[i, j].u = 1
385 else:
386 r[i, j].v = 1
387
388 ssarun.ssa.apply_linearization_transpose(r, TStarR)
389
390 r_global = PISM.IceModelVec2V()
391 r_global.create(grid, "", PISM.WITHOUT_GHOSTS)
392 r_global.copy_from(r)
393
394 for (k, l) in grid.points():
395 with PISM.vec.Access(nocomm=TStarR_indirect):
396 d = PISM.IceModelVec2S()
397 d.create(grid, "", PISM.WITH_GHOSTS)
398 d.set(0)
399 with PISM.vec.Access(comm=d):
400 d[k, l] = 1
401
402 ssarun.ssa.apply_linearization(d, Td)
403
404 TStarR_indirect[k, l] = Td.get_vec().dot(r_global.get_vec())
405
406 d_TStarR = PISM.IceModelVec2S()
407 d_TStarR.create(grid, "", PISM.WITHOUT_GHOSTS)
408
409 d_TStarR.copy_from(TStarR)
410 d_TStarR.add(-1, TStarR_indirect)
411
412 PISM.verbPrintf(1, grid.com, "\nTest Linearization Transpose (%d,%d):\n" % (i, j))
413
414 view(r_global, r_viewer)
415 view(TStarR, TStarR_viewer)
416 view(TStarR_indirect, TStarR_indirect_viewer)
417 view(d_TStarR, d_TStarR_viewer)
418
419 PISM.logging.pause()
420
421
422 # Main code starts here
423 if __name__ == "__main__":
424 context = PISM.Context()
425 config = context.config
426 com = context.com
427 PISM.set_abort_on_sigint(True)
428
429 append_mode = False
430
431 input_filename = config.get_string("input.file")
432 inv_data_filename = PISM.OptionString("-inv_data", "inverse data file", input_filename).value()
433 use_design_prior = config.get_flag("inverse.use_design_prior")
434 design_var = PISM.OptionKeyword("-inv_ssa",
435 "design variable for inversion",
436 "tauc,hardav", "tauc").value()
437 using_zeta_fixed_mask = config.get_flag("inverse.use_zeta_fixed_mask")
438
439 ssarun = PISM.invert.ssa.SSAForwardRunFromInputFile(input_filename, inv_data_filename, design_var)
440 ssarun.setup()
441
442 vecs = ssarun.modeldata.vecs
443 grid = ssarun.grid
444
445 # Determine the prior guess for tauc/hardav. This can be one of
446 # a) tauc/hardav from the input file (default)
447 # b) tauc/hardav_prior from the inv_datafile if -inv_use_design_prior is set
448 design_prior = createDesignVec(grid, design_var, '%s_prior' % design_var)
449 long_name = design_prior.string_attr("long_name")
450 units = design_prior.string_attr("units")
451 design_prior.set_attrs("",
452 "best prior estimate for %s (used for inversion)" % long_name,
453 units, units, "", 0)
454 if PISM.util.fileHasVariable(inv_data_filename, "%s_prior" % design_var) and use_design_prior:
455 PISM.logging.logMessage(" Reading '%s_prior' from inverse data file %s.\n" % (design_var, inv_data_filename))
456 design_prior.regrid(inv_data_filename, critical=True)
457 else:
458 if not PISM.util.fileHasVariable(input_filename, design_var):
459 PISM.verbPrintf(1, com, "Initial guess for design variable is not available as '%s' in %s.\nYou can provide an initial guess in the inverse data file.\n" % (
460 design_var, input_filename))
461 exit(1)
462 PISM.logging.logMessage("Reading '%s_prior' from '%s' in input file.\n" % (design_var, design_var))
463 design = createDesignVec(grid, design_var)
464 design.regrid(input_filename, True)
465 design_prior.copy_from(design)
466
467 if using_zeta_fixed_mask:
468 if PISM.util.fileHasVariable(inv_data_filename, "zeta_fixed_mask"):
469 zeta_fixed_mask = PISM.model.createZetaFixedMaskVec(grid)
470 zeta_fixed_mask.regrid(inv_data_filename)
471 vecs.add(zeta_fixed_mask)
472 else:
473 if design_var == 'tauc':
474 logMessage(
475 " Computing 'zeta_fixed_mask' (i.e. locations where design variable '%s' has a fixed value).\n" % design_var)
476 zeta_fixed_mask = PISM.model.createZetaFixedMaskVec(grid)
477 zeta_fixed_mask.set(1)
478 mask = vecs.ice_mask
479 with PISM.vec.Access(comm=zeta_fixed_mask, nocomm=mask):
480 for (i, j) in grid.points():
481 if mask.grounded_ice(i, j):
482 zeta_fixed_mask[i, j] = 0
483 vecs.add(zeta_fixed_mask)
484
485 adjustTauc(vecs.ice_mask, design_prior)
486 elif design_var == 'hardav':
487 pass
488 else:
489 raise NotImplementedError("Unable to build 'zeta_fixed_mask' for design variable %s.", design_var)
490
491 # Convert design_prior -> zeta_prior
492 zeta1 = PISM.IceModelVec2S()
493 zeta1.create(grid, "", PISM.WITH_GHOSTS, WIDE_STENCIL)
494 ssarun.designVariableParameterization().convertFromDesignVariable(design_prior, zeta1)
495
496 ssarun.ssa.linearize_at(zeta1)
497
498 test_type = PISM.OptionKeyword("-inv_test", "",
499 "j_design,j_design_transpose,lin,lin_transpose",
500 "j_design")
501
502 if not test_type.is_set():
503 PISM.verbPrintf(1, com, "Must specify a test type via -inv_test\n")
504 exit(1)
505 else:
506 test_type = test_type.value()
507
508 if test_type == "j_design":
509 test_j_design(ssarun)
510 elif test_type == "j_design_transpose":
511 test_j_design_transpose(ssarun)
512 elif test_type == "lin":
513 test_lin(ssarun)
514 elif test_type == "lin_transpose":
515 test_linearization_transpose(ssarun)