tssa_test_cfbc.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_cfbc.py (4654B) --- 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 help = \ 25 """ 26 SSA_TESTCFBC 27 Testing program for PISM's implementations of the SSA. 28 Does a time-independent calculation. Does not run IceModel or a derived 29 class thereof. Uses the van der Veen flow-line shelf geometry. Also may be 30 used in a PISM software (regression) test. 31 """ 32 33 usage = \ 34 """ 35 usage of SSA_TEST_CFBC: 36 run ssa_test_cfbc -Mx <number> -My <number> 37 """ 38 39 context = PISM.Context() 40 41 H0 = 600. # meters 42 V0 = convert(300, "m/year", "m/second") 43 C = 2.45e-18 # "typical constant ice parameter" 44 T = 400 # time used to compute the calving front location 45 46 Q0 = V0 * H0 47 Hc1 = 4. * C / Q0 48 Hc2 = 1. / (H0 ** 4) 49 50 51 def H_exact(x): 52 return (Hc1 * x + Hc2) ** (-1 / 4.) 53 54 55 def u_exact(x): 56 return Q0 / H_exact(x) 57 58 59 class test_cfbc(PISM.ssa.SSAExactTestCase): 60 61 def _initGrid(self): 62 self.grid = None 63 halfWidth = 250.0e3 # 500.0 km length 64 Lx = halfWidth 65 Ly = halfWidth 66 self.grid = PISM.IceGrid.Shallow(PISM.Context().ctx, Lx, Ly, 0, 0, 67 self.Mx, self.My, 68 PISM.CELL_CENTER, 69 PISM.Y_PERIODIC) 70 71 def _initPhysics(self): 72 config = self.config 73 74 config.set_number("flow_law.isothermal_Glen.ice_softness", 75 pow(1.9e8, -config.get_number("stress_balance.ssa.Glen_exponent"))) 76 config.set_flag("stress_balance.ssa.compute_surface_gradient_inward", False) 77 config.set_flag("stress_balance.calving_front_stress_bc", True) 78 config.set_string("stress_balance.ssa.flow_law", "isothermal_glen") 79 80 enthalpyconverter = PISM.EnthalpyConverter(config) 81 82 self.modeldata.setPhysics(enthalpyconverter) 83 84 def _initSSACoefficients(self): 85 self._allocStdSSACoefficients() 86 self._allocateBCs() 87 88 vecs = self.modeldata.vecs 89 90 vecs.tauc.set(0.0) # irrelevant 91 vecs.bedrock_altitude.set(-1000.0) # assures shelf is floating 92 93 EC = PISM.EnthalpyConverter(PISM.Context().config) 94 enth0 = EC.enthalpy(273.15, 0.01, 0) # 0.01 water fraction 95 vecs.enthalpy.set(enth0) 96 97 grid = self.grid 98 thickness = vecs.land_ice_thickness 99 surface = vecs.surface_altitude 100 bc_mask = vecs.bc_mask 101 vel_bc = vecs.vel_bc 102 ice_mask = vecs.mask 103 104 ocean_rho = self.config.get_number("constants.sea_water.density") 105 ice_rho = self.config.get_number("constants.ice.density") 106 107 x_min = grid.x(0) 108 with PISM.vec.Access(comm=[thickness, surface, bc_mask, vel_bc, ice_mask]): 109 for (i, j) in grid.points(): 110 x = grid.x(i) 111 if i != grid.Mx() - 1: 112 thickness[i, j] = H_exact(x - x_min) 113 ice_mask[i, j] = PISM.MASK_FLOATING 114 else: 115 thickness[i, j] = 0 116 ice_mask[i, j] = PISM.MASK_ICE_FREE_OCEAN 117 118 surface[i, j] = (1.0 - ice_rho / ocean_rho) * thickness[i, j] 119 120 if i == 0: 121 bc_mask[i, j] = 1 122 vel_bc[i, j].u = V0 123 vel_bc[i, j].v = 0. 124 else: 125 bc_mask[i, j] = 0 126 vel_bc[i, j].u = 0. 127 vel_bc[i, j].v = 0. 128 129 def exactSolution(self, i, j, x, y): 130 x_min = self.grid.x(0) 131 if i != self.grid.Mx() - 1: 132 u = u_exact(x - x_min) 133 else: 134 u = 0 135 return [u, 0] 136 137 138 if __name__ == '__main__': 139 140 config = PISM.Context().config 141 142 tc = test_cfbc(int(config.get_number("grid.Mx")), 143 int(config.get_number("grid.My"))) 144 145 tc.run(config.get_string("output.file_name"))