tssa_testi.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_testi.py (4355B)
---
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 import math
24
25 m_schoof = 10 # (pure number)
26 L_schoof = 40e3 # meters
27 aspect_schoof = 0.05 # (pure)
28 H0_schoof = aspect_schoof * L_schoof
29 # = 2000 m THICKNESS
30 B_schoof = 3.7e8 # Pa s^{1/3}; hardness
31 # given on p. 239 of Schoof; why so big?
32 p_schoof = 4.0 / 3.0 # = 1 + 1/n
33
34
35 class testi(PISM.ssa.SSAExactTestCase):
36
37 def _initGrid(self):
38 Mx = self.Mx
39 My = self.My
40 Ly = 3 * L_schoof # 300.0 km half-width (L=40.0km in Schoof's choice of variables)
41 Lx = max(60.0e3, ((Mx - 1) / 2.) * (2.0 * Ly / (My - 1)))
42 self.grid = PISM.IceGrid.Shallow(PISM.Context().ctx, Lx, Ly, 0, 0,
43 Mx, My,
44 PISM.CELL_CORNER,
45 PISM.NOT_PERIODIC)
46
47 def _initPhysics(self):
48 config = self.config
49 config.set_flag("basal_resistance.pseudo_plastic.enabled", False)
50
51 # irrelevant
52 enthalpyconverter = PISM.EnthalpyConverter(config)
53
54 config.set_string("stress_balance.ssa.flow_law", "isothermal_glen")
55 config.set_number("flow_law.isothermal_Glen.ice_softness", pow(
56 B_schoof, -config.get_number("stress_balance.ssa.Glen_exponent")))
57
58 self.modeldata.setPhysics(enthalpyconverter)
59
60 def _initSSACoefficients(self):
61 self._allocStdSSACoefficients()
62 self._allocateBCs()
63 vecs = self.modeldata.vecs
64
65 vecs.bc_mask.set(0)
66 vecs.thk.set(H0_schoof)
67 vecs.mask.set(PISM.MASK_GROUNDED)
68
69 # The finite difference code uses the following flag to treat
70 # the non-periodic grid correctly.
71 self.config.set_flag("stress_balance.ssa.compute_surface_gradient_inward", True)
72 self.config.set_number("stress_balance.ssa.epsilon", 0.0) # don't use this lower bound
73
74 standard_gravity = self.config.get_number("constants.standard_gravity")
75 ice_rho = self.config.get_number("constants.ice.density")
76 theta = math.atan(0.001)
77 f = ice_rho * standard_gravity * H0_schoof * math.tan(theta)
78 grid = self.grid
79 with PISM.vec.Access(comm=[vecs.tauc]):
80 for (i, j) in grid.points():
81 y = grid.y(j)
82 vecs.tauc[i, j] = f * (abs(y / L_schoof) ** m_schoof)
83
84 bc_mask = vecs.bc_mask
85 vel_bc = vecs.vel_bc
86 surface = vecs.surface_altitude
87 bed = vecs.bedrock_altitude
88 grid = self.grid
89 with PISM.vec.Access(comm=[surface, bed, vel_bc, bc_mask]):
90 for (i, j) in grid.points():
91 p = PISM.exactI(m_schoof, grid.x(i), grid.y(j))
92 bed[i, j] = p.bed
93 surface[i, j] = p.bed + H0_schoof
94
95 edge = ((j == 0) or (j == grid.My() - 1)) or ((i == 0) or (i == grid.Mx() - 1))
96 if edge:
97 bc_mask[i, j] = 1
98 vel_bc[i, j].u = p.u
99 vel_bc[i, j].v = p.v
100
101 def exactSolution(self, i, j, x, y):
102 p = PISM.exactI(m_schoof, x, y)
103 return [p.u, p.v]
104
105
106 # The main code for a run follows:
107 if __name__ == '__main__':
108 context = PISM.Context()
109 config = context.config
110
111 PISM.set_abort_on_sigint(True)
112
113 config.set_number("grid.Mx", 11, PISM.CONFIG_DEFAULT)
114 config.set_number("grid.My", 61, PISM.CONFIG_DEFAULT)
115
116 tc = testi(int(config.get_number("grid.Mx")),
117 int(config.get_number("grid.My")))
118 tc.run(config.get_string("output.file_name"))