tssa_testj.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_testj.py (4091B) --- 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 import PISM 22 from PISM.util import convert 23 24 class testj(PISM.ssa.SSAExactTestCase): 25 26 def _initGrid(self): 27 halfWidth = 300.0e3 28 Lx = halfWidth 29 Ly = halfWidth 30 ctx = PISM.Context().ctx 31 self.grid = PISM.IceGrid.Shallow(ctx, Lx, Ly, 0, 0, 32 self.Mx, self.My, 33 PISM.CELL_CENTER, 34 PISM.XY_PERIODIC) 35 36 def _initPhysics(self): 37 config = self.modeldata.config 38 config.set_flag("basal_resistance.pseudo_plastic.enabled", False) 39 40 enthalpyconverter = PISM.EnthalpyConverter(config) 41 42 config.set_string("stress_balance.ssa.flow_law", "isothermal_glen") 43 44 self.modeldata.setPhysics(enthalpyconverter) 45 46 def _initSSACoefficients(self): 47 self._allocStdSSACoefficients() 48 self._allocateBCs() 49 50 vecs = self.modeldata.vecs 51 52 vecs.tauc.set(0.0) # irrelevant for test J 53 # ensures that the ice is floating (max. thickness if 770 m) 54 vecs.bedrock_altitude.set(-1000.0) 55 vecs.mask.set(PISM.MASK_FLOATING) 56 vecs.bc_mask.set(0) # No dirichlet data. 57 58 EC = PISM.EnthalpyConverter(PISM.Context().config) 59 enth0 = EC.enthalpy(273.15, 0.01, 0) # 0.01 water fraction 60 vecs.enthalpy.set(enth0) 61 62 ocean_rho = self.config.get_number("constants.sea_water.density") 63 ice_rho = self.config.get_number("constants.ice.density") 64 65 # The PISM.vec.Access object ensures that we call beginAccess for each 66 # variable in 'vars', and that endAccess is called for each one on exiting 67 # the 'with' block. 68 69 with PISM.vec.Access(comm=[vecs.land_ice_thickness, 70 vecs.surface_altitude, 71 vecs.bc_mask, 72 vecs.vel_bc]): 73 grid = self.grid 74 for (i, j) in grid.points(): 75 p = PISM.exactJ(grid.x(i), grid.y(j)) 76 vecs.land_ice_thickness[i, j] = p.H 77 vecs.surface_altitude[i, j] = (1.0 - ice_rho / ocean_rho) * p.H # // FIXME task #7297 78 79 # special case at center point (Dirichlet BC) 80 if (i == grid.Mx() // 2) and (j == grid.My() // 2): 81 vecs.bc_mask[i, j] = 1 82 vecs.vel_bc[i, j] = [p.u, p.v] 83 84 def _initSSA(self): 85 # Test J has a viscosity that is independent of velocity. So we force a 86 # constant viscosity by settting the strength_extension 87 # thickness larger than the given ice thickness. (max = 770m). 88 89 nu0 = convert(30.0, "MPa year", "Pa s") 90 H0 = 500.0 # 500 m typical thickness 91 92 ssa = self.ssa 93 ssa.strength_extension.set_notional_strength(nu0 * H0) 94 ssa.strength_extension.set_min_thickness(800.) 95 96 def exactSolution(self, i, j, x, y): 97 p = PISM.exactJ(x, y) 98 return [p.u, p.v] 99 100 101 # The main code for a run follows: 102 if __name__ == '__main__': 103 context = PISM.Context() 104 config = context.config 105 106 tc = testj(int(config.get_number("grid.Mx")), int(config.get_number("grid.My"))) 107 tc.run(config.get_string("output.file_name"))