tssa_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
---
tssa_gn.py (3717B)
---
1 # Copyright (C) 2012, 2014, 2015, 2016, 2018 David Maxwell
2 #
3 # This file is part of PISM.
4 #
5 # PISM is free software; you can redistribute it and/or modify it under the
6 # terms of the GNU General Public License as published by the Free Software
7 # Foundation; either version 3 of the License, or (at your option) any later
8 # version.
9 #
10 # PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 # FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 # details.
14 #
15 # You should have received a copy of the GNU General Public License
16 # along with PISM; if not, write to the Free Software
17 # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
18
19 """Tikhonov Gauss-Newton inverse method code."""
20
21 import PISM
22 from PISM.invert.ssa import InvSSASolver
23
24 class InvSSASolver_TikhonovGN(InvSSASolver):
25 "Tikhonov Gauss-Newton inverse solver."
26 def __init__(self, ssarun, method):
27 self.solver = None
28 InvSSASolver.__init__(self, ssarun, method)
29
30 target_misfit = PISM.OptionReal("-inv_target_misfit",
31 "m/year; desired root misfit for inversions", 0.0)
32
33 if target_misfit.is_set():
34 self.target_misfit = target_misfit.value()
35 else:
36 raise RuntimeError("Missing required option -inv_target_misfit")
37
38 # FIXME: m_vel_scale is not defined (what are the units?)
39 self.target_misfit = self.target_misfit / m_vel_scale
40
41 self.listeners = []
42
43 def solveForward(self, zeta, out=None):
44 ssa = self.ssarun.ssa
45
46 reason = ssa.linearize_at(zeta)
47 if reason.failed():
48 raise PISM.AlgorithmFailureException(reason.description())
49 if out is not None:
50 out.copy_from(ssa.solution())
51 else:
52 out = ssa.solution()
53 return out
54
55 def solveInverse(self, zeta0, u_obs, zeta_inv):
56 eta = self.config.get_number("inverse.tikhonov.penalty_weight")
57
58 (designFunctional, stateFunctional) = PISM.invert.ssa.createTikhonovFunctionals(self.ssarun)
59 self.solver = PISM.IP_SSATaucTikhonovGNSolver(self.ssarun.ssa, zeta0, u_obs, eta, designFunctional, stateFunctional)
60
61 vel_scale = self.ssarun.grid.ctx().config().get_number("inverse.ssa.velocity_scale")
62 self.solver.setTargetMisfit(self.target_misfit / vel_scale)
63
64 # pl = [ TikhonovIterationListenerAdaptor(self,l) for l in self.listeners ]
65 pl = []
66 for l in pl:
67 self.solver.addListener(l)
68 self.solver.setInitialGuess(zeta_inv)
69
70 if PISM.OptionBool("-inv_test_adjoint", ""):
71 self.solver.init()
72 grid = self.ssarun.grid
73 d1 = PISM.vec.randVectorS(grid, 1)
74 d2 = PISM.vec.randVectorS(grid, 1)
75 y1 = PISM.IceModelVec2S()
76 y1.create(grid, '', PISM.WITHOUT_GHOSTS)
77 y2 = PISM.IceModelVec2S()
78 y2.create(grid, '', PISM.WITHOUT_GHOSTS)
79 self.solver.apply_GN(d1, y1)
80 self.solver.apply_GN(d2, y2)
81 ip1 = y1.get_vec().dot(d2.get_vec())
82 ip2 = y2.get_vec().dot(d1.get_vec())
83 PISM.logging.logMessage("ip1 %.10g ip2 %.10g\n" % (ip1, ip2))
84 PISM.logging.logMessage("ip1 %g ip2 %g\n" % (ip1, ip2))
85 exit(0)
86
87 vecs = self.ssarun.modeldata.vecs
88 if vecs.has('zeta_fixed_mask'):
89 self.ssarun.ssa.set_tauc_fixed_locations(vecs.zeta_fixed_mask)
90
91 return self.solver.solve()
92
93 def inverseSolution(self):
94 zeta = self.solver.designSolution()
95 u = self.solver.stateSolution()
96 return (zeta, u)