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"))