tssa_test_linear.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_test_linear.py (3737B) --- 1 #! /usr/bin/env python3 2 # 3 # Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2018 Ed Bueler and Constantine Khroulev and David Maxwell 4 # 5 # This file is part of PISM. 6 # 7 # PISM is free software; you can redistribute it and/or modify it under the 8 # terms of the GNU General Public License as published by the Free Software 9 # Foundation; either version 3 of the License, or (at your option) any later 10 # version. 11 # 12 # PISM is distributed in the hope that it will be useful, but WITHOUT ANY 13 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14 # FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15 # details. 16 # 17 # You should have received a copy of the GNU General Public License 18 # along with PISM; if not, write to the Free Software 19 # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA 20 21 22 import PISM 23 from PISM.util import convert 24 import math 25 26 context = PISM.Context() 27 28 L = 50.e3 # // 50km half-width 29 H0 = 500 # // m 30 dhdx = 0.005 # // pure number, slope of surface & bed 31 nu0 = convert(30.0, "MPa year", "Pa s") 32 tauc0 = 1.e4 # // 1kPa 33 34 35 class test_linear(PISM.ssa.SSAExactTestCase): 36 37 def _initGrid(self): 38 self.grid = PISM.IceGrid.Shallow(PISM.Context().ctx, L, L, 0, 0, 39 self.Mx, self.My, 40 PISM.CELL_CORNER, 41 PISM.NOT_PERIODIC) 42 43 def _initPhysics(self): 44 config = self.config 45 config.set_flag("basal_resistance.pseudo_plastic.enabled", True) 46 config.set_number("basal_resistance.pseudo_plastic.q", 1.0) 47 48 enthalpyconverter = PISM.EnthalpyConverter(config) 49 50 config.set_string("stress_balance.ssa.flow_law", "isothermal_glen") 51 52 self.modeldata.setPhysics(enthalpyconverter) 53 54 def _initSSACoefficients(self): 55 self._allocStdSSACoefficients() 56 self._allocateBCs() 57 58 vecs = self.modeldata.vecs 59 60 vecs.land_ice_thickness.set(H0) 61 vecs.surface_altitude.set(H0) 62 vecs.bedrock_altitude.set(0.) 63 vecs.tauc.set(tauc0) 64 65 vel_bc = vecs.vel_bc 66 bc_mask = vecs.bc_mask 67 bc_mask.set(0) 68 69 grid = self.grid 70 with PISM.vec.Access(comm=[bc_mask, vel_bc]): 71 for (i, j) in grid.points(): 72 edge = ((j == 0) or (j == grid.My() - 1)) or ((i == 0) or (i == grid.Mx() - 1)) 73 if edge: 74 bc_mask[i, j] = 1 75 x = grid.x(i) 76 y = grid.y(j) 77 [u, v] = self.exactSolution(i, j, x, y) 78 vel_bc[i, j].u = u 79 vel_bc[i, j].v = v 80 81 def _initSSA(self): 82 # The following ensure that the strength extension is used everywhere 83 se = self.ssa.strength_extension 84 se.set_notional_strength(nu0 * H0) 85 se.set_min_thickness(4000 * 10) 86 87 # For the benefit of SSAFD on a non-periodic grid 88 self.config.set_flag("stress_balance.ssa.compute_surface_gradient_inward", True) 89 90 def exactSolution(self, i, j, x, y): 91 tauc_threshold_velocity = self.config.get_number("basal_resistance.pseudo_plastic.u_threshold", 92 "m/second") 93 94 v0 = convert(100, "m/year", "m/second") 95 alpha = math.sqrt((tauc0 / tauc_threshold_velocity) / (4 * nu0 * H0)) 96 return [v0 * math.exp(-alpha * (x - L)), 0] 97 98 99 # The main code for a run follows: 100 if __name__ == '__main__': 101 context = PISM.Context() 102 config = context.config 103 104 PISM.set_abort_on_sigint(True) 105 106 tc = test_linear(int(config.get_number("grid.Mx")), int(config.get_number("grid.My"))) 107 tc.run(config.get_string("output.file_name"))