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