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