ttest_invssa_gn.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_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))