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)