test_invssa_gn.py (5837B)
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 import sys 22 import petsc4py 23 petsc4py.init(sys.argv) 24 from petsc4py import PETSc 25 import numpy as np 26 import os 27 import math 28 29 import PISM 30 31 32 def adjustTauc(mask, tauc): 33 """Where ice is floating or land is ice-free, tauc should be adjusted to have some preset default values.""" 34 35 grid = mask.grid() 36 high_tauc = grid.ctx().config().get_number("basal_yield_stress.ice_free_bedrock") 37 38 with PISM.vec.Access(comm=tauc, nocomm=mask): 39 for (i, j) in grid.points(): 40 if mask.ocean(i, j): 41 tauc[i, j] = 0 42 elif mask.ice_free(i, j): 43 tauc[i, j] = high_tauc 44 45 46 # Main code starts here 47 if __name__ == "__main__": 48 context = PISM.Context() 49 config = context.config 50 com = context.com 51 PISM.set_abort_on_sigint(True) 52 53 append_mode = False 54 55 input_filename = config.get_string("input.file") 56 inv_data_filename = PISM.OptionString("-inv_data", "inverse data file", input_filename).value() 57 use_tauc_prior = PISM.OptionBool("-inv_use_tauc_prior", 58 "Use tauc_prior from inverse data file as initial guess.") 59 60 ssarun = PISM.invert.ssa.SSAForwardRunFromInputFile(input_filename, inv_data_filename, 'tauc') 61 ssarun.setup() 62 63 vecs = ssarun.modeldata.vecs 64 grid = ssarun.grid 65 66 # Determine the prior guess for tauc. This can be one of 67 # a) tauc from the input file (default) 68 # b) tauc_prior from the inv_datafile if -use_tauc_prior is set 69 tauc_prior = PISM.model.createYieldStressVec(grid, 'tauc_prior') 70 tauc_prior.set_attrs("diagnostic", 71 "initial guess for (pseudo-plastic) basal yield stress in an inversion", 72 "Pa", "Pa", "", 0) 73 tauc = PISM.model.createYieldStressVec(grid) 74 if use_tauc_prior: 75 tauc_prior.regrid(inv_data_filename, critical=True) 76 else: 77 if not PISM.util.fileHasVariable(input_filename, "tauc"): 78 PISM.verbPrintf( 79 1, com, "Initial guess for tauc is not available as 'tauc' in %s.\nYou can provide an initial guess as 'tauc_prior' using the command line option -use_tauc_prior." % input_filename) 80 exit(1) 81 tauc.regrid(input_filename, True) 82 tauc_prior.copy_from(tauc) 83 84 adjustTauc(vecs.ice_mask, tauc_prior) 85 86 # Convert tauc_prior -> zeta_prior 87 zeta = PISM.IceModelVec2S() 88 WIDE_STENCIL = int(grid.ctx().config().get_number("grid.max_stencil_width")) 89 zeta.create(grid, "", PISM.WITH_GHOSTS, WIDE_STENCIL) 90 ssarun.tauc_param.convertFromDesignVariable(tauc_prior, zeta) 91 ssarun.ssa.linearize_at(zeta) 92 93 vel_ssa_observed = None 94 vel_ssa_observed = PISM.model.create2dVelocityVec(grid, '_ssa_observed', stencil_width=2) 95 if PISM.util.fileHasVariable(inv_data_filename, "u_ssa_observed"): 96 vel_ssa_observed.regrid(inv_data_filename, True) 97 else: 98 if not PISM.util.fileHasVariable(inv_data_filename, "u_surface_observed"): 99 PISM.verbPrintf( 100 1, context.com, "Neither u/v_ssa_observed nor u/v_surface_observed is available in %s.\nAt least one must be specified.\n" % inv_data_filename) 101 exit(1) 102 vel_surface_observed = PISM.model.create2dVelocityVec(grid, '_surface_observed', stencil_width=2) 103 vel_surface_observed.regrid(inv_data_filename, True) 104 105 sia_solver = PISM.SIAFD 106 if is_regional: 107 sia_solver = PISM.SIAFD_Regional 108 vel_sia_observed = PISM.sia.computeSIASurfaceVelocities(modeldata, sia_solver) 109 110 vel_sia_observed.metadata(0).set_name('u_sia_observed') 111 vel_sia_observed.metadata(0).set_string('long_name', "x-component of the 'observed' SIA velocities") 112 113 vel_sia_observed.metadata(1).set_name('v_sia_observed') 114 vel_sia_observed.metadata(1).set_string('long_name', "y-component of the 'observed' SIA velocities") 115 116 vel_ssa_observed.copy_from(vel_surface_observed) 117 vel_ssa_observed.add(-1, vel_sia_observed) 118 119 (designFunctional, stateFunctional) = PISM.invert.ssa.createTikhonovFunctionals(ssarun) 120 eta = config.get_number("inverse.tikhonov.penalty_weight") 121 122 solver_gn = PISM.InvSSATikhonovGN(ssarun.ssa, zeta, vel_ssa_observed, eta, designFunctional, stateFunctional) 123 124 seed = PISM.OptionInteger("-inv_seed", "random generator seed") 125 if seed.is_set(): 126 np.random.seed(seed.value() + PISM.Context().rank) 127 128 d1 = PISM.vec.randVectorS(grid, 1) 129 d2 = PISM.vec.randVectorS(grid, 1) 130 131 GNd1 = PISM.IceModelVec2S() 132 GNd1.create(grid, "", PISM.WITHOUT_GHOSTS) 133 GNd2 = PISM.IceModelVec2S() 134 GNd2.create(grid, "", PISM.WITHOUT_GHOSTS) 135 136 solver_gn.apply_GN(d1, GNd1) 137 solver_gn.apply_GN(d2, GNd2) 138 139 ip1 = d1.get_vec().dot(GNd2.get_vec()) 140 ip2 = d2.get_vec().dot(GNd1.get_vec()) 141 PISM.verbPrintf(1, grid.com, "Test of Gauss-Newton symmetry (x^t GN y) vs (y^t GN x)\n") 142 PISM.verbPrintf(1, grid.com, "ip1 %.10g ip2 %.10g\n" % (ip1, ip2)) 143 PISM.verbPrintf(1, grid.com, "relative error %.10g\n" % abs((ip1 - ip2) / ip1))