ssa.py (22070B)
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 and 20 SSA verification test cases.""" 21 22 import PISM 23 import math 24 from PISM import util, model 25 26 # Conversion from command-line arguments to classes of SSA solver. 27 SSAAlgorithms = {"fem": PISM.SSAFEM, "fd": PISM.SSAFD} 28 29 class SSARun(object): 30 31 """Mediates solving PISM's SSA model from a minimal set of data, without the constrution of an :cpp:class:`iceModel`. 32 It codifies the steps needed to put together the data for an SSA run; subclasses do the work of 33 implementing the steps in :meth:`_setFromOptions`, :meth:`_initGrid`, etc. Uses include: 34 35 * Running SSA test cases. 36 * Running the SSA in standalone mode (e.g. via :command:`ssaforward.py`) 37 * The SSA inversion code. 38 39 Usage: After construction (of a subclass), 40 41 1. Call :meth:`setup` to run through the various 42 steps needed to set up an environment for solving the SSA. 43 2. Solve the SSA with :meth:`solve`. 44 3. Optionally write the the model vectors and solution to a file with :meth:`write`.""" 45 46 def __init__(self): 47 """Do little constructor. Real work is done by :meth:`setup` which should be called prior to :meth:`solve`.""" 48 self.grid = None #: The computation grid; will be set by :meth:`_initGrid` 49 self.config = None #: Placeholder for config dictionary; set indirectly by :meth:`_constructModelData` 50 51 #: Instance of :class:`PISM.model.ModelData` that stores all data needed for solving the SSA. Much of the work of 52 #: the :class:`SSARun` is involved in setting up this object. Tasks include setting up :cpp:class:IceModelVec 53 #: variables as well as model physics (e.g. :cpp:class:`EnthalpyConverter`). 54 self.modeldata = None 55 self.ssa = None #: Subclass of :cpp:class:`SSA` that sovles the SSA. 56 57 def setup(self): 58 """Orchestrates the steps of setting up an environment for running the SSA. The following methods 59 are called in order, and should be impelmeneted by a subclass. 60 61 1. :meth:`_setFromOptions` to set any parameters from command-line options 62 2. :meth:`_initGrid` to determine the computation grid, to be stored as :attr:`grid` 63 3. :meth:`_constructModelData` provide a :class:`ModelData` object (a default implementation is provided) 64 4. :meth:`_initPhysics` to set the non-vec members of the :class:`ModelData`, e.g. the :cpp:class:`EnthalpyConverter`. 65 5. :meth:`_constructSSA` to build the actual subclass of :cpp:class:`SSA` that will be used to solve the SSA 66 6. :meth:`_initSSACoefficients` enter all of the vecs needed for solving the SSA into the :class:`ModelData`. 67 7. :meth:`_initSSA` initialize the :cpp:class:`SSA` returned in step 5 68 """ 69 self._setFromOptions() 70 71 self._initGrid() 72 if self.grid is None: 73 raise RuntimeError("SSARun failed to provide a grid.") 74 75 self.modeldata = self._constructModelData() 76 if self.modeldata is None: 77 raise RuntimeError("SSARun._constructModelData failed to provide a ModelData.") 78 self.config = self.modeldata.config 79 80 self._initPhysics() 81 if self.modeldata.enthalpyconverter is None: 82 raise RuntimeError("SSARun._initPhysics failed to initialize the physics of the underlying SSA solver.") 83 84 self.ssa = self._constructSSA() 85 if self.ssa is None: 86 raise RuntimeError("SSARun._constructSSA failed to provide an SSA.") 87 88 self._initSSACoefficients() 89 # FIXME: is there a reasonable check to do here? 90 91 self._initSSA() 92 93 def solve(self): 94 """Solve the SSA by calling the underlying PISM :cpp:class:`SSA`'s 95 :cpp:member:`update` method. Returns the solution vector (owned by 96 self.ssa, but you should not need to know about ownership). 97 98 """ 99 vecs = self.modeldata.vecs 100 101 # make sure vecs is locked! 102 self.ssa.init() 103 104 melange_back_pressure = PISM.IceModelVec2S() 105 melange_back_pressure.create(self.grid, "melange_back_pressure", PISM.WITHOUT_GHOSTS) 106 melange_back_pressure.set_attrs("diagnostic", 107 "melange back pressure fraction", "1", "1", "", 0) 108 melange_back_pressure.set(0.0) 109 110 PISM.verbPrintf(2, self.grid.com, "* Solving the SSA stress balance ...\n") 111 112 full_update = True 113 114 inputs = PISM.StressBalanceInputs() 115 inputs.melange_back_pressure = melange_back_pressure 116 inputs.geometry = self.geometry 117 inputs.enthalpy = vecs.enthalpy 118 inputs.basal_yield_stress = vecs.tauc 119 if vecs.has('vel_bc'): 120 inputs.bc_mask = vecs.bc_mask 121 inputs.bc_values = vecs.vel_bc 122 123 self.ssa.update(inputs, full_update) 124 125 return self.ssa.velocity() 126 127 def write(self, filename): 128 """Saves all of :attr:`modeldata`'s vecs (and the solution) to an 129 output file.""" 130 grid = self.grid 131 vecs = self.modeldata.vecs 132 133 pio = util.prepare_output(filename) 134 pio.close() 135 136 # Save time & command line 137 util.writeProvenance(filename) 138 139 vel_ssa = self.ssa.velocity() 140 vecs.add(vel_ssa) 141 142 velbar_mag = model.createCBarVec(self.grid) 143 velbar_mag.set_to_magnitude(vel_ssa) 144 velbar_mag.mask_by(vecs.thk, util.convert(-0.01, "m/year", "m/second")) 145 vecs.add(velbar_mag) 146 147 taud = PISM.SSA_taud(self.ssa).compute() 148 vecs.add(taud) 149 150 try: 151 nuH = PISM.SSAFD_nuH(self.ssa).compute() 152 vecs.add(nuH) 153 except: 154 pass 155 156 taud_mag = PISM.SSA_taud_mag(self.ssa).compute() 157 vecs.add(taud_mag) 158 159 vecs.writeall(filename) 160 161 def _setFromOptions(self): 162 """Optionally override to set any data from command line variables.""" 163 pass 164 165 def _constructModelData(self): 166 """Optionally override to return a custom :class:`PISM.model.ModelData` instance.""" 167 return model.ModelData(self.grid) 168 169 def _initGrid(self): 170 """Override to return the computation grid.""" 171 raise NotImplementedError() 172 173 def _initPhysics(self): 174 """Override to set the non-var parts of :attr:`modeldata` (e.g. the basal yeild stress model and the enthalpy converter)""" 175 raise NotImplementedError() 176 177 def _allocStdSSACoefficients(self): 178 """Helper method that allocates the standard :cpp:class:`IceModelVec` variables used to solve the SSA and stores them 179 in :attr:`modeldata```.vecs``: 180 181 * ``surface`` 182 * ``thickness`` 183 * ``bed`` 184 * ``tauc`` 185 * ``enthalpy`` 186 * ``mask`` 187 * ``age`` if -age is given 188 189 Intended to be called from custom implementations of :meth:`_initSSACoefficients` if desired.""" 190 vecs = self.modeldata.vecs 191 grid = self.grid 192 193 self.geometry = PISM.Geometry(grid) 194 geometry = self.geometry 195 196 vecs.add(geometry.ice_surface_elevation) 197 vecs.add(geometry.ice_thickness) 198 vecs.add(geometry.bed_elevation) 199 vecs.add(geometry.sea_level_elevation) 200 vecs.add(geometry.cell_type) 201 vecs.add(model.createYieldStressVec(grid), 'tauc') 202 vecs.add(model.createEnthalpyVec(grid), 'enthalpy') 203 204 # The SIA model might need the "age" field 205 if grid.ctx().config().get_flag("age.enabled"): 206 vecs.add(model.createAgeVec(grid), "age") 207 208 def _allocateBCs(self, velname='_bc', maskname='bc_mask'): 209 """Helper method that allocates standard Dirichlet data 210 :cpp:class:`IceModelVec` variable and stores them in 211 :attr:`modeldata` ``.vecs``: 212 213 * ``vel_bc`` 214 * ``bc_mask`` 215 216 """ 217 vecs = self.modeldata.vecs 218 vecs.add(model.create2dVelocityVec(self.grid, 219 name=velname, 220 desc='SSA velocity boundary condition', 221 intent='intent'), 222 "vel_bc") 223 vecs.add(model.createBCMaskVec(self.grid, name=maskname), 224 "bc_mask") 225 226 def _initSSACoefficients(self): 227 """Override to allocate and initialize all :cpp:class:`IceModelVec` variables in :attr:`modeldata` ``.vecs`` 228 needed for solving the SSA.""" 229 raise NotImplementedError() 230 231 def _constructSSA(self): 232 """Optionally override to return an instance of :cpp:class:`SSA` (e.g. :cpp:class:`SSAFD` or :cpp:class:`SSAFEM`) 233 that will be used for solving the SSA.""" 234 md = self.modeldata 235 return SSAAlgorithms[md.config.get_string("stress_balance.ssa.method")](md.grid) 236 237 def _initSSA(self): 238 """Optionally perform any final initialization of :attr:`ssa`.""" 239 pass 240 241 242 class SSAExactTestCase(SSARun): 243 244 """Base class for implmentation of specific SSA test cases. Provides a mechanism for comparing 245 computed and exact values. Simply construct with a grid size and then call :meth:`run`""" 246 247 def __init__(self, Mx, My): 248 """Initialize with a grid of the specified size.""" 249 SSARun.__init__(self) 250 self.Mx = Mx 251 self.My = My 252 253 # For convenience, provide a grid. It will get initialized later 254 # on when _initGrid is called by our setup method. 255 self.grid = None 256 257 def run(self, output_file): 258 """Main command intended to be called by whatever code executes the test case. 259 Calls :meth:`setup`, :meth:`solve`, :meth:`report`, and :meth:`write`.""" 260 self.setup() 261 self.solve() 262 self.report() 263 self.write(output_file) 264 265 def report(self): 266 """Compares computed and exact solution values and displays a summary report.""" 267 grid = self.grid 268 269 ssa_stdout = self.ssa.stdout_report() 270 PISM.verbPrintf(3, grid.com, ssa_stdout) 271 272 maxvecerr = 0.0 273 avvecerr = 0.0 274 avuerr = 0.0 275 avverr = 0.0 276 maxuerr = 0.0 277 maxverr = 0.0 278 279 if (self.config.get_flag("basal_resistance.pseudo_plastic.enabled") and 280 self.config.get_number("basal_resistance.pseudo_plastic.q") != 1.0): 281 PISM.verbPrintf(1, grid.com, "WARNING: numerical errors not valid for pseudo-plastic till\n") 282 PISM.verbPrintf(1, grid.com, "NUMERICAL ERRORS in velocity relative to exact solution:\n") 283 284 vel_ssa = self.ssa.velocity() 285 286 vel_ssa.begin_access() 287 288 exactvelmax = 0 289 gexactvelmax = 0 290 for (i, j) in self.grid.points(): 291 x = grid.x(i) 292 y = grid.y(j) 293 (uexact, vexact) = self.exactSolution(i, j, x, y) 294 exactnormsq = math.sqrt(uexact * uexact + vexact * vexact) 295 exactvelmax = max(exactnormsq, exactvelmax) 296 solution = vel_ssa[i, j] 297 uerr = abs(solution.u - uexact) 298 verr = abs(solution.v - vexact) 299 avuerr += uerr 300 avverr += verr 301 maxuerr = max(maxuerr, uerr) 302 maxverr = max(maxverr, verr) 303 vecerr = math.sqrt(uerr * uerr + verr * verr) 304 maxvecerr = max(maxvecerr, vecerr) 305 avvecerr = avvecerr + vecerr 306 307 vel_ssa.end_access() 308 309 N = grid.Mx() * grid.My() 310 gexactvelmax = PISM.GlobalMax(grid.com, exactvelmax) 311 gmaxuerr = PISM.GlobalMax(grid.com, maxuerr) 312 gmaxverr = PISM.GlobalMax(grid.com, maxverr) 313 gavuerr = PISM.GlobalSum(grid.com, avuerr) / N 314 gavverr = PISM.GlobalSum(grid.com, avverr) / N 315 gmaxvecerr = PISM.GlobalMax(grid.com, maxvecerr) 316 gavvecerr = PISM.GlobalSum(grid.com, avvecerr) / N 317 318 sys = grid.ctx().unit_system() 319 320 m_year = PISM.UnitConverter(sys, "m / second", "m / year") 321 322 if abs(gexactvelmax) > 0.0: 323 relative_vel_error = (gavvecerr / gexactvelmax) * 100.0 324 else: 325 relative_vel_error = 0.0 326 327 PISM.verbPrintf(1, grid.com, "velocity : maxvector prcntavvec maxu maxv avu avv\n") 328 PISM.verbPrintf(1, grid.com, 329 " %11.4f%13.5f%10.4f%10.4f%10.4f%10.4f\n", 330 m_year(gmaxvecerr), 331 relative_vel_error, 332 m_year(gmaxuerr), 333 m_year(gmaxverr), 334 m_year(gavuerr), 335 m_year(gavverr)) 336 PISM.verbPrintf(1, grid.com, "NUM ERRORS DONE\n") 337 338 def exactSolution(self, i, j, xi, xj): 339 """Override to provide the exact value of the solution at grid index (``i``, ``j``) with 340 coordinates (``xi``, ``xj``).""" 341 raise NotImplementedError() 342 343 def write(self, filename): 344 """Override of :meth:`SSARun.write`. Does all of the above, and saves a copy of the exact solution.""" 345 SSARun.write(self, filename) 346 347 grid = self.grid 348 exact = model.create2dVelocityVec(grid, name="_exact", desc="SSA exact solution", intent="diagnostic") 349 exact.begin_access() 350 for (i, j) in grid.points(): 351 exact[i, j] = self.exactSolution(i, j, grid.x(i), grid.y(j)) 352 exact.end_access() 353 exact.write(filename) 354 355 356 class SSAFromInputFile(SSARun): 357 358 """Class for running the SSA based on data provided in an input file.""" 359 360 def __init__(self, boot_file): 361 SSARun.__init__(self) 362 self.grid = None 363 self.config = PISM.Context().config 364 self.boot_file = boot_file 365 self.phi_to_tauc = False 366 self.is_regional = False 367 368 def _setFromOptions(self): 369 self.phi_to_tauc = PISM.OptionBool("-phi_to_tauc", 370 "Recompute pseudo yield stresses from till friction angles.") 371 self.is_regional = PISM.OptionBool("-regional", "enable 'regional' mode") 372 373 def _initGrid(self): 374 """Override of :meth:`SSARun._initGrid`.""" 375 # FIXME: allow specification of Mx and My different from what's 376 # in the boot_file. 377 378 if self.is_regional and (self.config.get_string("stress_balance.ssa.method") == "fem"): 379 registration = PISM.CELL_CORNER 380 else: 381 registration = PISM.CELL_CENTER 382 383 ctx = PISM.Context().ctx 384 385 pio = PISM.File(ctx.com(), self.boot_file, PISM.PISM_NETCDF3, PISM.PISM_READONLY) 386 self.grid = PISM.IceGrid.FromFile(ctx, pio, "enthalpy", registration) 387 pio.close() 388 389 def _initPhysics(self): 390 """Override of :meth:`SSARun._initPhysics` that sets the physics based on command-line flags.""" 391 config = self.config 392 393 enthalpyconverter = PISM.EnthalpyConverter(config) 394 395 if PISM.OptionString("-ssa_glen", "SSA flow law Glen exponent").is_set(): 396 config.set_string("stress_balance.ssa.flow_law", "isothermal_glen") 397 config.scalar_from_option("flow_law.isothermal_Glen.ice_softness", "ice_softness") 398 else: 399 config.set_string("stress_balance.ssa.flow_law", "gpbld") 400 401 self.modeldata.setPhysics(enthalpyconverter) 402 403 def _allocExtraSSACoefficients(self): 404 """Allocate storage for SSA coefficients.""" 405 vecs = self.modeldata.vecs 406 if util.fileHasVariable(self.boot_file, 'ssa_driving_stress_x'): 407 vecs.add(model.createDrivingStressXVec(self.grid)) 408 409 if util.fileHasVariable(self.boot_file, 'ssa_driving_stress_y'): 410 vecs.add(model.createDrivingStressYVec(self.grid)) 411 412 no_model_mask = None 413 # For a regional run we'll need no_model_mask, usurfstore, thkstore 414 if self.is_regional: 415 no_model_mask = model.createNoModelMaskVec(self.grid) 416 vecs.add(no_model_mask, 'no_model_mask') 417 vecs.add(model.createIceSurfaceStoreVec(self.grid)) 418 vecs.add(model.createIceThicknessStoreVec(self.grid)) 419 420 if self.config.get_flag('stress_balance.ssa.dirichlet_bc'): 421 vecs.add(model.create2dVelocityVec(self.grid, name='_ssa_bc', 422 desc='SSA velocity boundary condition', 423 intent='intent'), 424 "vel_ssa_bc") 425 426 if self.is_regional: 427 vecs.add(no_model_mask, 'bc_mask') 428 else: 429 vecs.add(model.createBCMaskVec(self.grid), 'bc_mask') 430 431 if self.phi_to_tauc: 432 vecs.add(PISM.model.createBasalMeltRateVec(self.grid)) 433 vecs.add(PISM.model.createTillPhiVec(self.grid)) 434 vecs.add(PISM.model.createBasalWaterVec(self.grid)) 435 436 def _initSSACoefficients(self): 437 """Override of :meth:`SSARun._initSSACoefficients` that initializes variables from the 438 contents of the input file.""" 439 # Build the standard thickness, bed, etc 440 self._allocStdSSACoefficients() 441 self._allocExtraSSACoefficients() 442 443 vecs = self.modeldata.vecs 444 445 thickness = vecs.land_ice_thickness 446 bed = vecs.bedrock_altitude 447 enthalpy = vecs.enthalpy 448 mask = vecs.mask 449 surface = vecs.surface_altitude 450 sea_level = vecs.sea_level 451 452 sea_level.set(0.0) 453 454 # Read in the PISM state variables that are used directly in the SSA solver 455 for v in [thickness, bed, enthalpy]: 456 v.regrid(self.boot_file, True) 457 458 # The SIA model might need the age field. 459 if self.config.get_flag("age.enabled"): 460 vecs.age.regrid(self.boot_file, True) 461 462 # variables mask and surface are computed from the geometry previously read 463 464 gc = PISM.GeometryCalculator(self.config) 465 gc.compute(sea_level, bed, thickness, mask, surface) 466 467 if util.fileHasVariable(self.boot_file, 'ssa_driving_stress_x'): 468 vecs.ssa_driving_stress_x.regrid(self.boot_file, critical=True) 469 470 if util.fileHasVariable(self.boot_file, 'ssa_driving_stress_y'): 471 vecs.ssa_driving_stress_y.regrid(self.boot_file, critical=True) 472 473 # For a regional run we'll need no_model_mask, usurfstore, thkstore 474 if self.is_regional: 475 vecs.no_model_mask.regrid(self.boot_file, True) 476 477 if util.fileHasVariable(self.boot_file, 'usurfstore'): 478 vecs.usurfstore.regrid(self.boot_file, True) 479 else: 480 vecs.usurfstore.copy_from(vecs.surface_altitude) 481 482 if util.fileHasVariable(self.boot_file, 'thkstore'): 483 vecs.thkstore.regrid(self.boot_file, True) 484 else: 485 vecs.thkstore.copy_from(vecs.land_ice_thickness) 486 487 # Compute yield stress from PISM state variables 488 # (basal melt rate, tillphi, and basal water height) 489 grid = self.grid 490 491 if self.phi_to_tauc: 492 for v in [vecs.bmr, vecs.tillphi, vecs.bwat]: 493 v.regrid(self.boot_file, True) 494 vecs.add(v) 495 496 if self.is_regional: 497 yieldstress = PISM.RegionalDefaultYieldStress(self.modeldata.grid) 498 else: 499 yieldstress = PISM.MohrCoulombYieldStress(self.modeldata.grid) 500 501 # make sure vecs is locked! 502 yieldstress.init() 503 yieldstress.set_till_friction_angle(vecs.tillphi) 504 yieldstress.update(0, 1) 505 vecs.tauc.copy_from(yieldstress.basal_material_yield_stress()) 506 else: 507 vecs.tauc.regrid(self.boot_file, True) 508 509 if self.config.get_flag('stress_balance.ssa.dirichlet_bc'): 510 has_u_ssa_bc = util.fileHasVariable(self.boot_file, 'u_ssa_bc') 511 has_v_ssa_bc = util.fileHasVariable(self.boot_file, 'v_ssa_bc') 512 513 if (not has_u_ssa_bc) or (not has_v_ssa_bc): 514 PISM.verbPrintf(2, grid.com, 515 "Input file '%s' missing Dirichlet boundary data u/v_ssa_bc;" 516 " using zero default instead." % self.boot_file) 517 vecs.vel_ssa_bc.set(0.0) 518 else: 519 vecs.vel_ssa_bc.regrid(self.boot_file, True) 520 521 if not self.is_regional: 522 bc_mask_name = vecs.bc_mask.metadata().get_string("short_name") 523 if util.fileHasVariable(self.boot_file, bc_mask_name): 524 vecs.bc_mask.regrid(self.boot_file, True) 525 else: 526 PISM.verbPrintf(2, grid.com, 527 "Input file '%s' missing Dirichlet location mask '%s'." 528 " Default to no Dirichlet locations." % (self.boot_file, bc_mask_name)) 529 vecs.bc_mask.set(0) 530 531 def _constructSSA(self): 532 """Constructs an instance of :cpp:class:`SSA` for solving the SSA based on command-line flags ``-regional`` and ``-ssa_method``""" 533 md = self.modeldata 534 if self.is_regional and (md.config.get_string("stress_balance.ssa.method") == "fd"): 535 algorithm = PISM.SSAFD_Regional 536 else: 537 algorithm = SSAAlgorithms[md.config.get_string("stress_balance.ssa.method")] 538 return algorithm(md.grid)