tssa_test_plug.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_plug.py (4045B)
---
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
23 H0 = 2000. # //m
24 L = 50.e3 # // 50km half-width
25 dhdx = 0.001 # // pure number, slope of surface & bed
26 tauc0 = 0. # // No basal shear stress
27 B0 = 3.7e8 # // Pa s^{1/3}; hardness
28 # // given on p. 239 of Schoof; why so big?
29 glen_n = 3.
30
31
32 class test_plug(PISM.ssa.SSAExactTestCase):
33
34 def _initGrid(self):
35 self.grid = PISM.IceGrid.Shallow(PISM.Context().ctx, L, L, 0, 0,
36 self.Mx, self.My,
37 PISM.CELL_CORNER,
38 PISM.NOT_PERIODIC)
39
40 def _initPhysics(self):
41 config = self.config
42
43 # // Enthalpy converter is irrelevant for this test.
44 enthalpyconverter = PISM.EnthalpyConverter(config)
45
46 # // Use constant hardness
47 config.set_string("stress_balance.ssa.flow_law", "isothermal_glen")
48 config.set_number("flow_law.isothermal_Glen.ice_softness", pow(B0, -glen_n))
49 config.set_number("stress_balance.ssa.Glen_exponent", glen_n)
50
51 self.modeldata.setPhysics(enthalpyconverter)
52
53 def _initSSACoefficients(self):
54 self._allocStdSSACoefficients()
55 self._allocateBCs()
56
57 vecs = self.modeldata.vecs
58
59 # Set constant coefficients.
60 vecs.land_ice_thickness.set(H0)
61 vecs.tauc.set(tauc0)
62 vecs.mask.set(PISM.MASK_GROUNDED)
63
64 bc_mask = vecs.bc_mask
65 vel_bc = vecs.vel_bc
66 bed = vecs.bedrock_altitude
67 surface = vecs.surface_altitude
68
69 grid = self.grid
70 with PISM.vec.Access(comm=[bc_mask, vel_bc, bed, surface]):
71 for (i, j) in grid.points():
72 x = grid.x(i)
73 y = grid.y(j)
74
75 bed[i, j] = -x * (dhdx)
76 surface[i, j] = bed[i, j] + H0
77
78 edge = ((j == 0) or (j == grid.My() - 1)) or ((i == 0) or (i == grid.Mx() - 1))
79 if edge:
80 bc_mask[i, j] = 1
81 [u, v] = self.exactSolution(i, j, x, y)
82 vel_bc(i, j).u = u
83 vel_bc(i, j).v = v
84
85 def _initSSA(self):
86 # Ensure we never use the strength extension.
87 self.ssa.strength_extension.set_min_thickness(H0 / 2)
88
89 # // The finite difference code uses the following flag to treat the non-periodic grid correctly.
90 # self.config.set_flag("stress_balance.ssa.compute_surface_gradient_inward", True);
91
92 # SSAFEM uses this (even though it has "ssafd" in its name)
93 self.config.set_number("stress_balance.ssa.epsilon", 0.0)
94
95 def exactSolution(self, i, j, x, y):
96 earth_grav = self.config.get_number("constants.standard_gravity")
97 ice_rho = self.config.get_number("constants.ice.density")
98 f = ice_rho * earth_grav * H0 * dhdx
99 ynd = y / L
100
101 u = 0.5 * (f ** 3) * (L ** 4) / ((B0 * H0) ** 3) * (1 - ynd ** 4)
102 return [u, 0]
103
104
105 # The main code for a run follows:
106 if __name__ == '__main__':
107 context = PISM.Context()
108 config = context.config
109
110 PISM.set_abort_on_sigint(True)
111
112 tc = test_plug(int(config.get_number("grid.Mx")), int(config.get_number("grid.My")))
113 tc.run(config.get_string("output.file_name"))