tssa_siple.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_siple.py (24210B) --- 1 # Copyright (C) 2012, 2014, 2015, 2016, 2017, 2018, 2019 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 """Inverse SSA solvers using the siple library.""" 20 21 import PISM 22 from PISM.logging import logError, logMessage 23 from PISM.util import convert, Bunch 24 25 if not PISM.imported_from_sphinx: 26 from petsc4py import PETSc 27 28 import sys 29 import math 30 import traceback 31 32 import siple 33 import PISM.invert.sipletools 34 from PISM.invert.ssa import InvSSASolver 35 from PISM.invert.sipletools import PISMLocalVector 36 from siple.gradient.forward import NonlinearForwardProblem 37 from siple.gradient.nonlinear import InvertNLCG, InvertIGN 38 39 siple.reporting.clear_loggers() 40 siple.reporting.add_logger(PISM.invert.sipletools.pism_logger) 41 siple.reporting.set_pause_callback(PISM.invert.sipletools.pism_pause) 42 43 WIDE_STENCIL = 2 44 45 46 class InvSSASolver_Gradient(InvSSASolver): 47 48 """Inverse SSA solver based on `siple` iterative gradient methods.""" 49 50 def __init__(self, ssarun, method): 51 """ 52 :param ssarun: The :class:`PISM.ssa.SSARun` defining the forward problem. See :class:`PISM.inv_ssa.SSAForwardRunFromInputFile`. 53 :param method: String describing the actual ``siple`` algorithm to use. One of ``sd``, ``nlcg`` or ``ign``.""" 54 InvSSASolver.__init__(self, ssarun, method) 55 56 monitor_adjoint = PISM.OptionBool("-inv_monitor_adjoint", 57 "Track accuracy of the adjoint during computation") 58 morozov_scale_factor = PISM.OptionReal("-inv_morozov_scale", 59 "Scale factor (>=1) for Morozov discrepancy principle", 60 1.0).value() 61 ls_verbose = PISM.OptionBool("-inv_ls_verbose", "Turn on a verbose linesearch.") 62 ign_theta = PISM.OptionReal("-ign_theta", "theta parameter for IGN algorithm", 0.5) 63 max_it = int(self.config.get_number("inverse.max_iterations")) 64 65 target_misfit = PISM.OptionReal("-inv_target_misfit", 66 "m/year; desired root misfit for inversions", 0.0) 67 68 if target_misfit.is_set(): 69 self.target_misfit = target_misfit.value() 70 else: 71 raise RuntimeError("Missing required option -inv_target_misfit") 72 73 velocity_scale = ssarun.grid.ctx().config().get_number("inverse.ssa.velocity_scale", "m/year") 74 self.target_misfit /= velocity_scale 75 76 self.forward_problem = SSAForwardProblem(ssarun) 77 78 # Determine the inversion algorithm, and set up its arguments. 79 if self.method == "ign": 80 Solver = InvertSSAIGN 81 else: 82 Solver = InvertSSANLCG 83 84 self.u_i = None 85 self.zeta_i = None 86 87 params = Solver.defaultParameters() 88 params.ITER_MAX = max_it 89 if self.method == "sd": 90 params.steepest_descent = True 91 params.ITER_MAX = 10000 92 elif self.method == "ign": 93 params.linearsolver.ITER_MAX = 10000 94 params.linearsolver.verbose = True 95 if ls_verbose: 96 params.linesearch.verbose = True 97 params.verbose = True 98 params.thetaMax = ign_theta 99 params.mu = morozov_scale_factor 100 101 params.deriv_eps = 0. 102 103 self.siple_solver = Solver(self.forward_problem, params=params) 104 105 if monitor_adjoint: 106 self.addIterationListener(MonitorAdjoint()) 107 if self.method == 'ign': 108 self.addLinearIterationListener(MonitorAdjointLin()) 109 110 def solveForward(self, zeta, out=None): 111 r"""Given a parameterized design variable value :math:`\zeta`, solve the SSA. 112 See :cpp:class:`IP_TaucParam` for a discussion of parameterizations. 113 114 :param zeta: :cpp:class:`IceModelVec` containing :math:`\zeta`. 115 :param out: optional :cpp:class:`IceModelVec` for storage of the computation result. 116 :returns: An :cpp:class:`IceModelVec` contianing the computation result. 117 """ 118 if out is None: 119 out = self.forward_problem.F(PISMLocalVector(zeta)) 120 else: 121 out = self.forward_problem.F(PISMLocalVector(zeta), out=PISMLocalVector(out)) 122 return out.core() 123 124 def addIterationListener(self, listener): 125 """Add a listener to be called after each iteration. See FIXME.""" 126 self.siple_solver.addIterationListener(SipleIterationListenerAdaptor(self, listener)) 127 128 def addDesignUpdateListener(self, listener): 129 self.siple_solver.addXUpdateListener(SipleDesignUpdateListenerAdaptor(self, listener)) 130 131 def addLinearIterationListener(self, listener): 132 "Add a linear iteration listener." 133 self.siple_solver.addLinearIterationListener(SipleLinearIterationListenerAdaptor(self, listener)) 134 135 def solveInverse(self, zeta0, u_obs, zeta_inv): 136 r"""Executes the inversion algorithm. 137 138 :param zeta0: The best `a-priori` guess for the value of the parameterized design variable value :math:`\zeta`. 139 Ignored by siple gradient algorithms in deference to `zeta_inv`. 140 :param u_obs: :cpp:class:`IceModelVec2V` of observed surface velocities. 141 :param zeta_inv: :cpp:class:`zeta_inv` starting value of :math:`\zeta` for iterative minimization. 142 :returns: A :cpp:class:`TerminationReason`. 143 """ 144 try: 145 vecs = self.ssarun.modeldata.vecs 146 if vecs.has('zeta_fixed_mask'): 147 self.ssarun.ssa.set_tauc_fixed_locations(vecs.zeta_fixed_mask) 148 149 (self.zeta_i, self.u_i) = self.siple_solver.solve(zeta_inv, u_obs, self.target_misfit) 150 except Exception: 151 import traceback 152 exc_type, exc_value, exc_traceback = sys.exc_info() 153 description = "" 154 for l in traceback.format_exception(exc_type, exc_value, exc_traceback): 155 description += l 156 # It would be nice to make siple so that if the inverse solve fails 157 # you can still keep the most recent iteration. 158 self.u_i = None 159 self.zeta_i = None 160 return PISM.GenericTerminationReason(-1, description) 161 return PISM.GenericTerminationReason(1, "Morozov Discrepancy Met") 162 163 def inverseSolution(self): 164 """Returns a tuple ``(zeta, u)`` of :cpp:class:`IceModelVec`'s corresponding to the values 165 of the design and state variables at the end of inversion.""" 166 return (self.zeta_i, self.u_i) 167 168 169 class SSAForwardProblem(NonlinearForwardProblem): 170 171 """Subclass of a :class:`siple.NonlinearForwardProblem` defining the forward problem for ``siple``. 172 The heavy lifting is forwarded to a :cpp:class:`IP_SSATaucForwardProblem` or 173 :cpp:class:`IP_SSAHardnessForwardProblem` living inside a :class:`PISM.invert.ssa.SSATaucForwardRun`.""" 174 175 def __init__(self, ssarun): 176 """ 177 :param ssarun: 178 179 A :class:`PISM.invert.ssa.SSAForwardRun` 180 or :class:`PISM.invert.ssa.SSAForwardRunFromInputFile` defining 181 the forward problem.""" 182 183 self.ssarun = ssarun 184 self.ssa = self.ssarun.ssa 185 self.grid = ssarun.grid 186 187 self.tmpV = PISM.IceModelVec2V() 188 self.tmpV.create(self.grid, "work vector (2V)", PISM.WITHOUT_GHOSTS, WIDE_STENCIL) 189 190 self.tmpS = PISM.IceModelVec2S() 191 self.tmpS.create(self.grid, "work vector (2S)", PISM.WITHOUT_GHOSTS, WIDE_STENCIL) 192 193 self.tmpS2 = PISM.IceModelVec2S() 194 self.tmpS2.create(self.grid, "work vector (2S)", PISM.WITHOUT_GHOSTS, WIDE_STENCIL) 195 196 ksp_rtol = 1e-12 197 self.ksp = PETSc.KSP() 198 self.ksp.create(self.grid.com) 199 self.ksp.setTolerances(ksp_rtol, PETSc.DEFAULT, PETSc.DEFAULT, PETSc.DEFAULT) 200 self.ksp.getPC().setType('bjacobi') 201 self.ksp.setFromOptions() 202 203 (self.designFunctional, self.stateFunctional) = PISM.invert.ssa.createGradientFunctionals(ssarun) 204 205 self.designForm = None 206 207 def F(self, x, out=None, guess=None): 208 """ 209 Returns the value of the forward problem at the design variable ``x``. 210 211 Nonlinear problems often make use of an initial guess; this can be provided in ``guess``. 212 213 Storage in ``out``, if given, is used for the return value. 214 """ 215 if out is None: 216 out = self.rangeVector() 217 reason = self.ssa.linearize_at(x.core()) 218 if reason.failed(): 219 raise PISM.AlgorithmFailureException(reason.description()) 220 out.core().copy_from(self.ssa.solution()) 221 return out 222 223 def T(self, d, out=None): 224 """ 225 Returns the value of the linearization :math:`T` of the forward problem at the design variable :math:`x` 226 specified previously in :meth:`linearizeAt`, in the direction `d`. 227 228 Storage in `out`, if given, is used for the return value. 229 """ 230 if out is None: 231 out = self.rangeVector() 232 self.ssa.apply_linearization(d.core(), out.core()) 233 return out 234 235 def TStar(self, r, out=None): 236 """ 237 Let :math:`T` be the linearization of the forward problem :math:`F` at the design variable `x` (as specified previously in :meth:`linearizeAt`). 238 Its adjoint is :math:`T^*`. This method returns the value of :math:`T^*` in the direction `r`. 239 240 Storage in `out`, if given, is used for the return value. 241 """ 242 if out is None: 243 out = self.domainVector() 244 245 if self.designForm is None: 246 stencil_width = int(self.grid.ctx().config().get_number("grid.max_stencil_width")) 247 da2 = self.grid.get_dm(1, stencil_width) 248 249 if PISM.PETSc.Sys.getVersion() < (3, 5, 0): 250 self.designForm = da2.get().getMatrix("baij") 251 else: 252 da2.get().setMatType("baij") 253 self.designForm = da2.get().getMatrix() 254 255 self.designFunctional.assemble_form(self.designForm) 256 257 # First step 258 self.stateFunctional.gradientAt(r.core(), self.tmpV) 259 self.tmpV.scale(0.5) 260 self.ssa.apply_linearization_transpose(self.tmpV, self.tmpS) 261 262 # Second step 263 if PISM.PETSc.Sys.getVersion() < (3, 5, 0): 264 self.ksp.setOperators(self.designForm, self.designForm, 265 PETSc.Mat.Structure.SAME_NONZERO_PATTERN) 266 else: 267 self.ksp.setOperators(self.designForm, self.designForm) 268 269 self.ksp.solve(self.tmpS.vec(), self.tmpS2.vec()) 270 271 reason = self.ksp.getConvergedReason() 272 if reason < 0: 273 raise RuntimeError('TStarB linear solve failed to converge (KSP reason %s)\n\n' % reason) 274 else: 275 PISM.logging.logPrattle("TStarB converged (KSP reason %s)\n" % reason) 276 277 out.core().copy_from(self.tmpS2) 278 return out 279 280 def linearizeAt(self, x, guess=None): 281 """ 282 Instructs the class that subsequent calls to T and TStar will be conducted for the given value of x. 283 284 Nonlinear problems often make use of an initial guess; this can be provided in 'guess'. 285 """ 286 reason = self.ssa.linearize_at(x.core()) 287 if reason.failed(): 288 raise Exception(reason.description()) 289 290 def evalFandLinearize(self, x, out=None, guess=None): 291 """ 292 Computes the value of F(x) and locks in a linearization. Sometimes there are efficiencies that 293 can be acheived this way. 294 295 Default implementation simply calls F, then linearizeAt. 296 """ 297 if out is None: 298 out = self.rangeVector() 299 self.linearizeAt(x) 300 out.core().copy_from(self.ssa.solution()) 301 return out 302 303 def rangeIP(self, a, b): 304 """ 305 Computes the inner product of two vectors in the range (i.e. state) space. 306 """ 307 return self.stateFunctional.dot(a.core(), b.core()) 308 309 def domainIP(self, a, b): 310 """ 311 Computes the inner product of two vectors in the domain (i.e. design) space. 312 """ 313 return self.designFunctional.dot(a.core(), b.core()) 314 315 def rangeVector(self): 316 """Constructs a brand new vector from the range (i.e. state) vector space""" 317 v = PISM.IceModelVec2V() 318 v.create(self.grid, "", True, WIDE_STENCIL) 319 320 # Add appropriate meta data. 321 intent = "?inverse?" # FIXME 322 desc = "SSA velocity computed by inversion" 323 v.set_attrs(intent, "%s%s" % ("X-component of the ", desc), "m s-1", "m year-1", "", 0) 324 v.set_attrs(intent, "%s%s" % ("Y-component of the ", desc), "m s-1", "m year-1", "", 1) 325 326 huge_vel = convert(1e6, "m/year", "m/second") 327 attrs = [("valid_min", -huge_vel), ("valid_max", huge_vel), ("_FillValue", 2 * huge_vel)] 328 for a in attrs: 329 for component in range(2): 330 v.metadata(component).set_number(a[0], a[1]) 331 332 return PISMLocalVector(v) 333 334 def domainVector(self): 335 """Constructs a brand new vector from the domain (i.e. design) vector space""" 336 v = PISM.IceModelVec2S() 337 v.create(self.grid, "", True, WIDE_STENCIL) 338 return PISMLocalVector(v) 339 340 341 class InvertSSANLCG(InvertNLCG): 342 343 r"""Subclass of :class:`siple.gradient.nonlinear.InvertNLCG` for inversion of SSA velocities 344 from :math:`\tau_c` or hardness.""" 345 346 @staticmethod 347 def defaultParameters(): 348 params = InvertNLCG.defaultParameters() 349 return params 350 351 def __init__(self, forward_problem, params=None): 352 """:param forward_problem: A :class:`SSATaucForwardProblem` that defines the forward problem. 353 :param params: A :class:`siple.params.Parameters` containing run-time parameters. 354 """ 355 InvertNLCG.__init__(self, params) 356 self.forward_problem = forward_problem 357 self.misfit_goal = 0.0 358 359 def forwardProblem(self): 360 """:returns: the associated :class:`SSATaucForwardProblem` provided at construction""" 361 return self.forward_problem 362 363 def stopConditionMet(self, count, x, Fx, y, r): 364 """ 365 Determines if minimization should be halted (based, e.g. on a Morozov discrepancy principle) 366 367 :param count: current iteration count 368 :param x: point in domain of potential minimizer. 369 :param Fx: value of nonlinear function at `x` 370 :param r: current residual, i.e. :math:`y-F(x)` 371 :returns: boolean, `True` if termination condition is met 372 """ 373 374 misfit = math.sqrt(abs(self.forward_problem.rangeIP(r, r))) 375 376 if (misfit < self.misfit_goal): 377 siple.reporting.msg('Stop condition met') 378 return True 379 return False 380 381 def initialize(self, x, y, deltaLInf): 382 """ 383 Hook called at the start of :meth:`solve`. This gives the class a chance to massage the input. 384 385 The remaining arguments are passed directly from solve, and can be used for determining the 386 final stopping criterion. 387 388 Returns vectors corresponding to the initial value of `x` and the desired value of `y` where :math:`y=F(x)`. 389 """ 390 xv = PISMLocalVector(x) 391 yv = PISMLocalVector(y) 392 393 self.misfit_goal = self.params.mu * deltaLInf 394 395 return (xv, yv) 396 397 def finalize(self, x, y): 398 """ 399 Hook called at the end of :meth:`solve`. Gives the chance to massage the return values. 400 """ 401 zeta = x.core() 402 u = y.core() 403 return (zeta, u) 404 405 406 class InvertSSAIGN(InvertIGN): 407 408 r"""Subclass of :class:`siple.gradient.nonlinear.InvertIGN` for inversion of SSA velocities 409 from :math:`\tau_c` or hardness.""" 410 411 @staticmethod 412 def defaultParameters(): 413 params = InvertIGN.defaultParameters() 414 return params 415 416 def __init__(self, forward_problem, params=None): 417 """:param forward_problem: A :class:`SSATaucForwardProblem` that defines the forward problem. 418 :param params: A :class:`siple.params.Parameters` containing run-time parameters. 419 """ 420 InvertIGN.__init__(self, params) 421 self.forward_problem = forward_problem 422 423 def forwardProblem(self): 424 """:returns: the associated :class:`SSATaucForwardProblem` provided at construction""" 425 return self.forward_problem 426 427 def temper_d(self, x, d, y, r): 428 """Method called during iterations in :meth:`solve` to ensure we don't make to large a step. 429 :param x: current design parameter 430 :param d: step in design space 431 :param y: desired value of state paramter 432 :param r: residual of desired and current state parameter values (:math:`y-F(x)`) 433 434 Changes `d` as a side-effect to temper step size. 435 """ 436 dnorm = d.norm('linf') 437 xnorm = x.norm('linf') 438 if dnorm > 2 * xnorm: 439 siple.reporting.msg('wild change predicted by linear step. scaling') 440 d.scale(2 * xnorm / dnorm) 441 442 def initialize(self, x, y, target_misfit): 443 """ 444 Hook called at the start of :meth:`solve`. This gives the class a chance to massage the input. 445 446 The remaining arguments are passed directly from solve, and can be used for determining the 447 final stopping criterion. 448 449 Returns vectors corresponding to the initial value of `x` and the desired value of `y` where :math:`y=F(x)`. 450 """ 451 xv = PISMLocalVector(x) 452 yv = PISMLocalVector(y) 453 454 return (xv, yv, target_misfit) 455 456 def finalize(self, x, y): 457 """ 458 Hook called at the end of :meth:`solve`. Gives the chance to massage the return values. 459 """ 460 461 zeta = x.core() 462 u = y.core() 463 464 return (zeta, u) 465 466 467 class SipleIterationListenerAdaptor(object): 468 469 """Adaptor for passing listening events from `siple`-based solvers to a python object.""" 470 471 def __init__(self, owner, listener): 472 """:param owner: The :class:`InvSSATaucSolver_Tikhonov` that constructed us 473 :param listener: The python-based listener. 474 """ 475 self.owner = owner 476 self.listener = listener 477 478 def __call__(self, siplesolver, it, x, Fx, y, d, r, *args): 479 """Callback from `siple`. Gathers together the long list of arguments 480 into a dictionary and passes it along in a standard form to the python listener.""" 481 data = Bunch(zeta=x.core(), u=Fx.core(), zeta_step=d.core(), residual=r.core(), target_misfit=self.owner.target_misfit) 482 483 if self.owner.method == 'ign': 484 data.update(T_zeta_step=args[0].core()) 485 else: 486 data.update(TStar_residual=args[0].core()) 487 try: 488 self.listener(self.owner, it, data) 489 except Exception: 490 logError("\nERROR: Exception occured during an inverse solver listener callback:\n\n") 491 traceback.print_exc(file=sys.stdout) 492 raise 493 494 class SipleLinearIterationListenerAdaptor(object): 495 496 """Adaptor for passing listening events the linear steps of `siple`-based `ign` solvers to a python object.""" 497 498 def __init__(self, owner, listener): 499 """:param owner: The :class:`InvSSATaucSolver_Tikhonov` that constructed us 500 :param listener: The python-based listener. 501 """ 502 self.owner = owner 503 self.listener = listener 504 505 def __call__(self, siplesolver, it, x, y, d, r, Td, TStarR): 506 """Callback from `siple`. Gathers together the long list of arguments 507 into a dictionary and passes it along in a standard form to the python listener.""" 508 509 data = Bunch(x=x.core(), y=y.core(), r=r.core(), d=d.core(), Td=Td.core(), TSTarR=TStarR.core()) 510 try: 511 self.listener(self.owner, it, data) 512 except Exception: 513 logError("\nERROR: Exception occured during an inverse solver listener callback:\n\n") 514 traceback.print_exc(file=sys.stdout) 515 raise 516 517 class SipleDesignUpdateListenerAdaptor(object): 518 519 """Adaptor for design variable update events of `siple`-based solvers to a python object.""" 520 521 def __init__(self, owner, listener): 522 """:param owner: The :class:`InvSSATaucSolver_Tikhonov` that constructed us 523 :param listener: The python-based listener. 524 """ 525 self.owner = owner 526 self.listener = listener 527 528 def __call__(self, siplesolver, it, zeta, u, u_obs, r): 529 """Callback from `siple`. Gathers together the long list of arguments 530 into a dictionary and passes it along in a standard form to the python listener.""" 531 data = Bunch(zeta=zeta.core(), u=u.core(), r=r.core(), u_obs=u_obs.core()) 532 try: 533 self.listener(self.owner, it, data) 534 except Exception as e: 535 logError("\nWARNING: Exception occured during an inverse solver DesignUpdate listener callback:\n%s\n\n" % str(e)) 536 537 538 class MonitorAdjoint(object): 539 540 r"""Iteration listener that can be used to verify the correctness of the implementation of an adjoint. 541 For adjoint-based interative inverse methods, a residual ``r`` is known in state space and a step direction ``d`` is 542 known in design space. A linearized forward problem :math:`T` maps from design space to state space, and its adjoint :math:`T^*` 543 goes in the opposite direction. The inner products :math:`\left<Td,r\right>_{\rm State}` 544 and :math:`\left<d,T^*r\right>_{\rm Design}` should always be the same; this is a good diagnostic to determine 545 of an adjoint has been coded correctly. The listener prints a comparison of the values of the two inner products 546 at each iteration. 547 """ 548 549 def __init__(self): 550 self.Td = None 551 self.TStarR = None 552 self.didWarning = False 553 554 def __call__(self, inverse_solver, count, data): 555 """ 556 :param inverse_sovler: the solver (e.g. :class:`~InvSolver_Tikhonov`) we are listening to. 557 :param count: the iteration number. 558 :param data: dictionary of data related to the iteration. 559 """ 560 method = inverse_solver.method 561 if method != 'sd' and method != 'nlcg' and method != 'ign': 562 if not self.didWarning: 563 PISM.verbPrintf(1, PISM.Context().com, '\nWarning: unable to monitor adjoint for inverse method: %s\nOption -inv_monitor_adjoint ignored\n' % method) 564 self.didWarning = True 565 return 566 fp = inverse_solver.forward_problem 567 d = PISM.invert.sipletools.PISMLocalVector(data.d) 568 r = PISM.invert.sipletools.PISMLocalVector(data.r) 569 self.Td = fp.T(d, self.Td) 570 self.TStarR = fp.TStar(r, out=self.TStarR) 571 ip1 = fp.domainIP(d, self.TStarR) 572 ip2 = fp.rangeIP(self.Td, r) 573 logMessage("adjoint test: <Td,r>=%g <d,T^*r>=%g (percent error %g)", ip1, ip2, (abs(ip1 - ip2)) / max(abs(ip1), abs(ip2))) 574 575 576 class MonitorAdjointLin(object): 577 578 def __init__(self): 579 self.Td = None 580 self.TStarR = None 581 582 def __call__(self, inverse_solver, count, data): 583 """ 584 :param inverse_sovler: the solver (e.g. :class:`~InvSSASolver_Tikhonov`) we are listening to. 585 :param count: the iteration number. 586 :param data: dictionary of data related to the iteration. 587 """ 588 fp = inverse_solver.forward_problem 589 r = PISM.invert.sipletools.PISMLocalVector(data.r) 590 d = PISM.invert.sipletools.PISMLocalVector(data.d) 591 self.Td = fp.T(d, self.Td) 592 self.TStarR = fp.TStar(r, out=self.TStarR) 593 ip1 = fp.domainIP(d, self.TStarR) 594 ip2 = fp.rangeIP(self.Td, r) 595 logMessage("adjoint test: <Td,r>=%g <d,T^*r>=%g (percent error %g)", ip1, ip2, (abs(ip1 - ip2)) / max(abs(ip1), abs(ip2)))