tssa.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.py (27174B)
---
1 # Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 David Maxwell and Constantine Khroulev
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 """Module containing classes managing SSA forward runs, inverse SSA
20 solves, and iteration reporting."""
21
22 import PISM
23 import math
24
25 from PISM.logging import logMessage
26
27 class SSAForwardRun(PISM.ssa.SSARun):
28
29 """Subclass of :class:`PISM.ssa.SSAFromInputFile` where the underlying SSA implementation is an
30 :cpp:class:`IP_SSATaucForwardProblem` or :cpp:class:`IP_SSAHardavForwardProblem`.
31 It is responsible for putting together a :class:`PISM.model.ModelData` containing the auxilliary data
32 needed for solving the SSA (:cpp:class:`IceModelVec`'s, :cpp:class:`EnthalpyConverter`, etc.) as well
33 as the instance of :cpp:class:`IP_SSATaucForwardProblem` that will solve the SSA repeatedly in the course
34 of solving an inverse problem. This class is intended to be subclassed by test cases where the data
35 is not provided from an input file. See also :class:`SSAForwardRunFromInputFile`."""
36
37 def __init__(self, design_var):
38 PISM.ssa.SSARun.__init__(self)
39 assert(design_var in list(ssa_forward_problems.keys()))
40 self.grid = None
41 self.config = PISM.Context().config
42 self.design_var = design_var
43 self.design_var_param = createDesignVariableParam(self.config, self.design_var)
44 self.is_regional = False
45
46 def designVariable(self):
47 """:returns: String description of the design variable of the forward problem (e.g. 'tauc' or 'hardness')"""
48 return self.design_var
49
50 def designVariableParameterization(self):
51 """:returns: Object that performs zeta->design variable transformation."""
52 return self.design_var_param
53
54 def _setFromOptions(self):
55 """Initialize internal parameters based on command-line flags. Called from :meth:`PISM.ssa.SSARun.setup`."""
56 self.is_regional = PISM.OptionBool("-regional", "regional mode")
57
58 def _constructSSA(self):
59 """Returns an instance of :cpp:class:`IP_SSATaucForwardProblem` rather than
60 a basic :cpp:class:`SSAFEM` or :cpp:class:`SSAFD`. Called from :meth:`PISM.ssa.SSARun.setup`."""
61 md = self.modeldata
62 return createSSAForwardProblem(md.grid, md.enthalpyconverter, self.design_var_param, self.design_var)
63
64 def _initSSA(self):
65 """One-time initialization of the :cpp:class:`IP_SSATaucForwardProblem`. Called from :meth:`PISM.ssa.SSARun.setup`."""
66 # init() will cache the values of the coefficeints at
67 # quadrature points once here. Subsequent solves will then not
68 # need to cache these values.
69 self.ssa.init()
70
71
72 class SSAForwardRunFromInputFile(SSAForwardRun):
73
74 """Subclass of :class:`SSAForwardRun` where the vector data
75 for the run is provided in an input :file:`.nc` file."""
76
77 def __init__(self, input_filename, inv_data_filename, design_var):
78 """
79 :param input_filename: :file:`.nc` file containing generic PISM model data.
80 :param inv_data_filename: :file:`.nc` file containing data specific to inversion (e.g. observed SSA velocities).
81 """
82 SSAForwardRun.__init__(self, design_var)
83 self.input_filename = input_filename
84 self.inv_data_filename = inv_data_filename
85
86 def _initGrid(self):
87 """Initialize grid size and periodicity. Called from :meth:`PISM.ssa.SSARun.setup`."""
88
89 if self.is_regional:
90 registration = PISM.CELL_CORNER
91 else:
92 registration = PISM.CELL_CENTER
93
94 ctx = PISM.Context().ctx
95
96 pio = PISM.File(ctx.com(), self.input_filename, PISM.PISM_NETCDF3, PISM.PISM_READONLY)
97 self.grid = PISM.IceGrid.FromFile(ctx, pio, "enthalpy", registration)
98 pio.close()
99
100 def _initPhysics(self):
101 """Override of :meth:`SSARun._initPhysics` that sets the physics based on command-line flags."""
102 config = self.config
103
104 enthalpyconverter = PISM.EnthalpyConverter(config)
105
106 if PISM.OptionBool("-ssa_glen", "SSA flow law Glen exponent"):
107 config.set_string("stress_balance.ssa.flow_law", "isothermal_glen")
108 config.scalar_from_option("flow_law.isothermal_Glen.ice_softness", "ice_softness")
109 else:
110 config.set_string("stress_balance.ssa.flow_law", "gpbld")
111
112 self.modeldata.setPhysics(enthalpyconverter)
113
114 def _initSSACoefficients(self):
115 """Reads SSA coefficients from the input file. Called from :meth:`PISM.ssa.SSARun.setup`."""
116 self._allocStdSSACoefficients()
117
118 # Read PISM SSA related state variables
119 #
120 # Hmmm. A lot of code duplication with SSAFromInputFile._initSSACoefficients.
121
122 vecs = self.modeldata.vecs
123 thickness = vecs.land_ice_thickness
124 bed = vecs.bedrock_altitude
125 enthalpy = vecs.enthalpy
126 mask = vecs.mask
127 surface = vecs.surface_altitude
128 sea_level = vecs.sea_level
129
130 sea_level.set(0.0)
131
132 # Read in the PISM state variables that are used directly in the SSA solver
133 for v in [thickness, bed, enthalpy]:
134 v.regrid(self.input_filename, True)
135
136 # variables mask and surface are computed from the geometry previously read
137
138 gc = PISM.GeometryCalculator(self.config)
139 gc.compute(sea_level, bed, thickness, mask, surface)
140
141 grid = self.grid
142 config = self.modeldata.config
143
144 # Compute yield stress from PISM state variables
145 # (basal melt rate, tillphi, and basal water height) if they are available
146
147 file_has_inputs = (PISM.util.fileHasVariable(self.input_filename, 'bmelt') and
148 PISM.util.fileHasVariable(self.input_filename, 'tillwat') and
149 PISM.util.fileHasVariable(self.input_filename, 'tillphi'))
150
151 if file_has_inputs:
152 bmr = PISM.model.createBasalMeltRateVec(grid)
153 tillphi = PISM.model.createTillPhiVec(grid)
154 tillwat = PISM.model.createBasalWaterVec(grid)
155 for v in [bmr, tillphi, tillwat]:
156 v.regrid(self.input_filename, True)
157 vecs.add(v)
158
159 # The SIA model might need the age field.
160 if self.config.get_flag("age.enabled"):
161 vecs.age.regrid(self.input_filename, True)
162
163 hydrology_model = config.get_string("hydrology.model")
164 if hydrology_model == "null":
165 subglacial_hydrology = PISM.NullTransportHydrology(grid)
166 elif hydrology_model == "routing":
167 subglacial_hydrology = PISM.RoutingHydrology(grid)
168 elif hydrology_model == "distributed":
169 subglacial_hydrology = PISM.DistributedHydrology(grid)
170
171 if self.is_regional:
172 yieldstress = PISM.RegionalDefaultYieldStress(self.modeldata.grid, subglacial_hydrology)
173 else:
174 yieldstress = PISM.MohrCoulombYieldStress(self.modeldata.grid, subglacial_hydrology)
175
176 # make sure vecs is locked!
177 subglacial_hydrology.init()
178 yieldstress.init()
179
180 yieldstress.basal_material_yield_stress(vecs.tauc)
181 elif PISM.util.fileHasVariable(self.input_filename, 'tauc'):
182 vecs.tauc.regrid(self.input_filename, critical=True)
183
184 if PISM.util.fileHasVariable(self.input_filename, 'ssa_driving_stress_x'):
185 vecs.add(PISM.model.createDrivingStressXVec(self.grid))
186 vecs.ssa_driving_stress_x.regrid(self.input_filename, critical=True)
187
188 if PISM.util.fileHasVariable(self.input_filename, 'ssa_driving_stress_y'):
189 vecs.add(PISM.model.createDrivingStressYVec(self.grid))
190 vecs.ssa_driving_stress_y.regrid(self.input_filename, critical=True)
191
192 # read in the fractional floatation mask
193 vecs.add(PISM.model.createGroundingLineMask(self.grid))
194 vecs.gl_mask.regrid(self.input_filename, critical=False, default_value=0.0) # set to zero if not found
195
196 if self.is_regional:
197 vecs.add(PISM.model.createNoModelMaskVec(self.grid), 'no_model_mask')
198 vecs.no_model_mask.regrid(self.input_filename, True)
199 vecs.add(vecs.surface_altitude, 'usurfstore')
200
201 if self.config.get_flag('stress_balance.ssa.dirichlet_bc'):
202 vecs.add(PISM.model.create2dVelocityVec(self.grid, name='_ssa_bc', desc='SSA velocity boundary condition', intent='intent'), "vel_ssa_bc")
203 has_u_ssa_bc = PISM.util.fileHasVariable(self.input_filename, 'u_ssa_bc')
204 has_v_ssa_bc = PISM.util.fileHasVariable(self.input_filename, 'v_ssa_bc')
205 if (not has_u_ssa_bc) or (not has_v_ssa_bc):
206 PISM.verbPrintf(2, self.grid.com, "Input file '%s' missing Dirichlet boundary data u/v_ssa_bc; using zero default instead." % self.input_filename)
207 vecs.vel_ssa_bc.set(0.)
208 else:
209 vecs.vel_ssa_bc.regrid(self.input_filename, True)
210
211 if self.is_regional:
212 vecs.add(vecs.no_model_mask, 'bc_mask')
213 else:
214 vecs.add(PISM.model.createBCMaskVec(self.grid), 'bc_mask')
215 bc_mask_name = vecs.bc_mask.metadata().get_string("short_name")
216 if PISM.util.fileHasVariable(self.input_filename, bc_mask_name):
217 vecs.bc_mask.regrid(self.input_filename, True)
218 else:
219 PISM.verbPrintf(2, self.grid.com, "Input file '%s' missing Dirichlet location mask '%s'. Default to no Dirichlet locations." % (self.input_filename, bc_mask_name))
220 vecs.bc_mask.set(0)
221
222 if PISM.util.fileHasVariable(self.inv_data_filename, 'vel_misfit_weight'):
223 vecs.add(PISM.model.createVelocityMisfitWeightVec(self.grid))
224 vecs.vel_misfit_weight.regrid(self.inv_data_filename, True)
225
226
227 class InvSSASolver(object):
228
229 """Abstract base class for SSA inverse problem solvers."""
230
231 def __init__(self, ssarun, method):
232 """
233 :param ssarun: The :class:`PISM.invert.ssa.SSAForwardRun` defining the forward problem.
234 :param method: String describing the actual algorithm to use. Must be a key in :attr:`tao_types`."""
235
236 self.ssarun = ssarun
237 self.config = ssarun.config
238 self.method = method
239
240 def solveForward(self, zeta, out=None):
241 r"""Given a parameterized design variable value :math:`\zeta`, solve the SSA.
242 See :cpp:class:`IP_TaucParam` for a discussion of parameterizations.
243
244 :param zeta: :cpp:class:`IceModelVec` containing :math:`\zeta`.
245 :param out: optional :cpp:class:`IceModelVec` for storage of the computation result.
246 :returns: An :cpp:class:`IceModelVec` contianing the computation result.
247 """
248 raise NotImplementedError()
249
250 def addIterationListener(self, listener):
251 """Add a listener to be called after each iteration. See :ref:`Listeners`."""
252 raise NotImplementedError()
253
254 def addDesignUpdateListener(self, listener):
255 """Add a listener to be called after each time the design variable is changed."""
256 raise NotImplementedError()
257
258 def solveInverse(self, zeta0, u_obs, zeta_inv):
259 r"""Executes the inversion algorithm.
260
261 :param zeta0: The best `a-priori` guess for the value of the parameterized design variable :math:`\zeta`.
262 :param u_obs: :cpp:class:`IceModelVec2V` of observed surface velocities.
263 :param zeta_inv: :cpp:class:`zeta_inv` starting value of :math:`\zeta` for minimization of the Tikhonov functional.
264 :returns: A :cpp:class:`TerminationReason`.
265 """
266 raise NotImplementedError()
267
268 def inverseSolution(self):
269 """Returns a tuple ``(zeta, u)`` of :cpp:class:`IceModelVec`'s corresponding to the values
270 of the design and state variables at the end of inversion."""
271 raise NotImplementedError()
272
273
274 def createInvSSASolver(ssarun, method=None):
275 """Factory function returning an inverse solver appropriate for the config variable ``inverse.ssa.method``.
276
277 :param ssarun: an instance of :class:`SSAForwardRun:` or :class:`SSAForwardRunFromInputFile`.
278 :param method: a string correpsonding to config variable ``inverse.ssa.method`` describing the inversion method to be used.
279 """
280 if method is None:
281 method = ssarun.config.get_string('inverse.ssa.method')
282 if method == 'tikhonov_gn':
283 from PISM.invert import ssa_gn
284 return ssa_gn.InvSSASolver_TikhonovGN(ssarun, method)
285 elif method.startswith('tikhonov'):
286 try:
287 from PISM.invert import ssa_tao
288 return ssa_tao.InvSSASolver_Tikhonov(ssarun, method)
289 except ImportError:
290 raise RuntimeError("Inversion method '%s' requires the TAO library." % method)
291
292 if method == 'sd' or method == 'nlcg' or method == 'ign':
293 try:
294 from PISM.invert import ssa_siple
295 return ssa_siple.InvSSASolver_Gradient(ssarun, method)
296 except ImportError:
297 raise RuntimeError("Inversion method '%s' requires the siple python library." % method)
298
299 raise Exception("Unknown inverse method '%s'; unable to construct solver.", method)
300
301
302 design_param_types = {"ident": PISM.IPDesignVariableParamIdent,
303 "square": PISM.IPDesignVariableParamSquare,
304 "exp": PISM.IPDesignVariableParamExp,
305 "trunc": PISM.IPDesignVariableParamTruncatedIdent}
306
307
308 def createDesignVariableParam(config, design_var_name, param_name=None):
309 """Factory function for creating subclasses of :cpp:class:`IPDesignVariableParameterization` based on command-line flags."""
310 if param_name is None:
311 param_name = config.get_string("inverse.design.param")
312 design_param = design_param_types[param_name]()
313 design_param.set_scales(config, design_var_name)
314 return design_param
315
316 ssa_forward_problems = {'tauc': PISM.IP_SSATaucForwardProblem,
317 'hardav': PISM.IP_SSAHardavForwardProblem}
318
319
320 def createSSAForwardProblem(grid, ec, design_param, design_var):
321 """Returns an instance of an SSA forward problem (e.g. :cpp:class:`IP_SSATaucForwardProblem`)
322 suitable for the value of `design_var`"""
323 ForwardProblem = ssa_forward_problems[design_var]
324 if ForwardProblem is None:
325 raise RuntimeError("Design variable %s is not yet supported.", design_var)
326
327 return ForwardProblem(grid, design_param)
328
329
330 def createGradientFunctionals(ssarun):
331 """Returns a tuple ``(designFunctional,stateFunctional)`` of :cpp:class:`IP_IPFunctional`'s
332 for gradient-based inversions. The specific functionals are constructed on the basis of
333 command-line parameters ``inverse.state_func`` and ``inverse.design.func``.
334
335 :param ssarun: The instance of :class:`PISM.ssa.SSARun` that encapsulates the forward problem,
336 typically a :class:`SSAForwardRunFromFile`.
337 """
338
339 vecs = ssarun.modeldata.vecs
340 grid = ssarun.grid
341
342 useGroundedIceOnly = PISM.OptionBool("-inv_ssa_grounded_ice_tauc",
343 "Computed norms for tau_c only on elements with all grounded ice.")
344
345 misfit_type = grid.ctx().config().get_string("inverse.state_func")
346 if misfit_type != 'meansquare':
347 inv_method = grid.ctx().config().get_string("inverse.ssa.method")
348 raise Exception("'-inv_state_func %s' is not supported with '-inv_method %s'.\nUse '-inv_state_func meansquare' instead" % (misfit_type, inv_method))
349
350 design_functional = grid.ctx().config().get_string("inverse.design.func")
351 if design_functional != "sobolevH1":
352 inv_method = grid.ctx().config().get_string("inverse.ssa.method")
353 raise Exception("'-inv_design_func %s' is not supported with '-inv_method %s'.\nUse '-inv_design_func sobolevH1' instead" % (design_functional, inv_method))
354
355 designFunctional = createHilbertDesignFunctional(grid, vecs, useGroundedIceOnly)
356
357 stateFunctional = createMeanSquareMisfitFunctional(grid, vecs)
358
359 return (designFunctional, stateFunctional)
360
361
362 def createTikhonovFunctionals(ssarun):
363 """Returns a tuple ``(designFunctional,stateFunctional)`` of :cpp:class:`IP_Functional`'s
364 for Tikhonov-based inversions. The specific functionals are constructed on the basis of
365 command-line parameters ``inv_state_func`` and ``inv_design_func``.
366
367 :param ssarun: The instance of :class:`PISM.ssa.SSARun` that encapsulates the forward problem,
368 typically a :class:`SSATaucForwardRunFromFile`.
369 """
370 vecs = ssarun.modeldata.vecs
371 grid = ssarun.grid
372
373 useGroundedIceOnly = PISM.OptionBool("-inv_ssa_grounded_ice_tauc",
374 "Computed norms for tau_c only on elements with all grounded ice.")
375
376 misfit_type = grid.ctx().config().get_string("inverse.state_func")
377 if misfit_type == "meansquare":
378 stateFunctional = createMeanSquareMisfitFunctional(grid, vecs)
379 elif misfit_type == "log_ratio":
380 vel_ssa_observed = vecs.vel_ssa_observed
381 scale = grid.ctx().config().get_number("inverse.log_ratio_scale")
382 velocity_eps = grid.ctx().config().get_number("inverse.ssa.velocity_eps", "m/second")
383 misfit_weight = None
384 if vecs.has('vel_misfit_weight'):
385 misfit_weight = vecs.vel_misfit_weight
386 stateFunctional = PISM.IPLogRatioFunctional(grid, vel_ssa_observed, velocity_eps, misfit_weight)
387 stateFunctional.normalize(scale)
388 elif misfit_type == "log_relative":
389 vel_ssa_observed = vecs.vel_ssa_observed
390 velocity_scale = grid.ctx().config().get_number("inverse.ssa.velocity_scale", "m/second")
391 velocity_eps = grid.ctx().config().get_number("inverse.ssa.velocity_eps", "m/second")
392 misfit_weight = None
393 if vecs.has('vel_misfit_weight'):
394 misfit_weight = vecs.vel_misfit_weight
395 stateFunctional = PISM.IPLogRelativeFunctional(grid, vel_ssa_observed, velocity_eps, misfit_weight)
396 stateFunctional.normalize(velocity_scale)
397 else:
398 raise RuntimeError("Unknown inv_state_func '%s'; unable to construct solver.", misfit_type)
399
400 design_functional = grid.ctx().config().get_string("inverse.design.func")
401 if design_functional == "sobolevH1":
402 designFunctional = createHilbertDesignFunctional(grid, vecs, useGroundedIceOnly)
403 elif design_functional == "tv":
404 area = 4 * grid.Lx() * grid.Ly()
405 velocity_scale = grid.ctx().config().get_number("inverse.ssa.velocity_scale", "m/second")
406 length_scale = grid.ctx().config().get_number("inverse.ssa.length_scale")
407 lebesgue_exponent = grid.ctx().config().get_number("inverse.ssa.tv_exponent")
408 cTV = 1 / area
409 cTV *= (length_scale) ** (lebesgue_exponent)
410
411 zeta_fixed_mask = None
412 if vecs.has('zeta_fixed_mask'):
413 zeta_fixed_mask = vecs.zeta_fixed_mask
414
415 strain_rate_eps = PISM.OptionReal("-inv_ssa_tv_eps",
416 "regularization constant for 'total variation' functional",
417 0.0)
418 if not strain_rate_eps.is_set():
419 schoofLen = grid.ctx().config().get_number("flow_law.Schoof_regularizing_length", "m")
420 strain_rate_eps = 1 / schoofLen
421 else:
422 strain_rate_eps = strain_rate_eps.value()
423
424 designFunctional = PISM.IPTotalVariationFunctional2S(grid, cTV, lebesgue_exponent, strain_rate_eps, zeta_fixed_mask)
425 else:
426 raise Exception("Unknown inv_design_func '%s'; unable to construct solver." % design_functional)
427
428 return (designFunctional, stateFunctional)
429
430
431 def createMeanSquareMisfitFunctional(grid, vecs):
432 """Creates a :cpp:class:`IPMeanSquareFunctional2V` suitable for use for a
433 state variable function for SSA inversions."""
434
435 misfit_weight = None
436 if vecs.has('vel_misfit_weight'):
437 misfit_weight = vecs.vel_misfit_weight
438
439 velocity_scale = grid.ctx().config().get_number("inverse.ssa.velocity_scale", "m/second")
440 stateFunctional = PISM.IPMeanSquareFunctional2V(grid, misfit_weight)
441 stateFunctional.normalize(velocity_scale)
442 return stateFunctional
443
444
445 def createHilbertDesignFunctional(grid, vecs, useGroundedIceOnly):
446 """Creates a :cpp:class:`IP_H1NormFunctional2S` or a :cpp:class`IPGroundedIceH1NormFunctional2S` suitable
447 for use for a design variable functional.
448
449 :param grid: computation grid
450 :param vecs: model vecs
451 :param useGroundedIceOnly: flag, ``True`` if a :cpp:class`IPGroundedIceH1NormFunctional2S` should be created.
452 """
453 cL2 = grid.ctx().config().get_number("inverse.design.cL2")
454 cH1 = grid.ctx().config().get_number("inverse.design.cH1")
455
456 area = 4 * grid.Lx() * grid.Ly()
457 length_scale = grid.ctx().config().get_number("inverse.ssa.length_scale")
458 cL2 /= area
459 cH1 /= area
460 cH1 *= (length_scale * length_scale)
461
462 zeta_fixed_mask = None
463 if vecs.has('zeta_fixed_mask'):
464 zeta_fixed_mask = vecs.zeta_fixed_mask
465
466 if useGroundedIceOnly:
467 mask = vecs.mask
468 designFunctional = PISM.IPGroundedIceH1NormFunctional2S(grid, cL2, cH1, mask, zeta_fixed_mask)
469 else:
470 designFunctional = PISM.IP_H1NormFunctional2S(grid, cL2, cH1, zeta_fixed_mask)
471
472 return designFunctional
473
474
475 def printIteration(invssa_solver, it, data):
476 "Print a header for an iteration report."
477 logMessage("----------------------------------------------------------\n")
478 logMessage("Iteration %d\n" % it)
479
480
481 def printTikhonovProgress(invssasolver, it, data):
482 "Report on the progress of a Tikhonov iteration."
483 eta = data.tikhonov_penalty
484 stateVal = data.JState
485 designVal = data.JDesign
486 sWeight = 1
487 dWeight = 1.0 / eta
488
489 norm_type = PISM.PETSc.NormType.NORM_2
490
491 logMessage("design objective %.8g; weighted %.8g\n" % (designVal, designVal * dWeight))
492 if 'grad_JTikhonov' in data:
493 logMessage("gradient: design %.8g state %.8g sum %.8g\n" % (data.grad_JDesign.norm(norm_type) * dWeight,
494 data.grad_JState.norm(norm_type) * sWeight,
495 data.grad_JTikhonov.norm(norm_type)))
496 else:
497 logMessage("gradient: design %.8g state %.8g; constraints: %.8g\n" % (data.grad_JDesign.norm(norm_type) * dWeight,
498 data.grad_JState.norm(norm_type) * sWeight,
499 data.constraints.norm(norm_type)))
500 logMessage("tikhonov functional: %.8g\n" % (stateVal * sWeight + designVal * dWeight))
501
502
503 class RMSMisfitReporter(object):
504 "Report RMS misfit."
505 def __init__(self):
506 self.J = None
507
508 def __call__(self, invssa_solver, it, data):
509
510 grid = invssa_solver.ssarun.grid
511
512 if self.J is None:
513 vecs = invssa_solver.ssarun.modeldata.vecs
514 self.J = createMeanSquareMisfitFunctional(grid, vecs)
515
516 Jmisfit = self.J.valueAt(data.residual)
517 rms_misfit = math.sqrt(Jmisfit) * grid.ctx().config().get_number("inverse.ssa.velocity_scale")
518
519 PISM.logging.logMessage("Diagnostic RMS Misfit: %0.8g (m/a)\n" % rms_misfit)
520
521
522 class MisfitLogger(object):
523 "Logger that saves history of misfits to a file."
524 def __init__(self):
525 self.misfit_history = []
526 self.misfit_type = None
527
528 def __call__(self, invssa_solver, it, data):
529 """
530 :param inverse_solver: the solver (e.g. :class:`~InvSSASolver_Tikhonov`) we are listening to.
531 :param count: the iteration number.
532 :param data: dictionary of data related to the iteration.
533 """
534
535 grid = invssa_solver.ssarun.grid
536
537 if self.misfit_type is None:
538 self.misfit_type = grid.ctx().config().get_string("inverse.state_func")
539
540 method = invssa_solver.method
541 if method == 'ign' or method == 'sd' or method == 'nlcg':
542 import PISM.invert.sipletools
543 fp = invssa_solver.forward_problem
544 r = PISM.invert.sipletools.PISMLocalVector(data.residual)
545 Jmisfit = fp.rangeIP(r, r)
546 elif 'JState' in data:
547 Jmisfit = data.JState
548 else:
549 raise RuntimeError("Unable to report misfits for inversion method: %s" % method)
550
551 if self.misfit_type == "meansquare":
552 velScale_m_per_year = grid.ctx().config().get_number("inverse.ssa.velocity_scale")
553
554 rms_misfit = math.sqrt(Jmisfit) * velScale_m_per_year
555
556 logMessage("Misfit: sqrt(J_misfit) = %.8g (m/a)\n" % rms_misfit)
557 self.misfit_history.append(rms_misfit)
558 else:
559 logMessage("Misfit: J_misfit = %.8g (dimensionless)\n" % Jmisfit)
560 self.misfit_history.append(Jmisfit)
561
562 def write(self, output_filename):
563 """Saves a history of misfits as :ncvar:`inv_ssa_misfit`
564
565 :param output_filename: filename to save misfits to."""
566
567 N = len(self.misfit_history)
568
569 ds = PISM.File(PISM.Context().com, output_filename, PISM.PISM_NETCDF3, PISM.PISM_READWRITE)
570
571 ds.redef()
572 ds.define_dimension('inv_ssa_iter', N)
573 ds.define_variable("inv_ssa_misfit", PISM.PISM_DOUBLE, ["inv_ssa_iter"])
574 if self.misfit_type == "meansquare":
575 ds.write_attribute("inv_ssa_misfit", "units", "m/a")
576 ds.write_variable("inv_ssa_misfit", [0], [N], self.misfit_history)
577 ds.close()
578
579
580 class ZetaSaver(object):
581 r"""Iteration listener used to save a copy of the current value
582 of :math:`\zeta` (i.e. a parameterized design variable such as :math:`\tau_c` or hardness)
583 at each iteration during an inversion. The intent is to use a saved value to restart
584 an inversion if need be.
585 """
586
587 def __init__(self, output_filename):
588 """:param output_filename: file to save iterations to."""
589 self.output_filename = output_filename
590
591 def __call__(self, inverse_solver, count, data):
592 zeta = data.zeta
593 # The solver doesn't care what the name of zeta is, and we
594 # want it called 'zeta_inv' in the output file, so we rename it.
595 zeta.metadata().set_name('zeta_inv')
596 zeta.metadata().set_string('long_name',
597 'last iteration of parameterized basal yeild stress computed by inversion')
598 zeta.write(self.output_filename)