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 (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)