texactQ.py - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
texactQ.py (6150B)
---
1 #!/usr/bin/env python3
2
3 # Computes exact similarity solution "test Q" at an initial time and after
4 # a specified runtime. Suggests a PISM run for this runtime. The exact
5 # solution for velocity and for thickness can be compared to PISM output.
6 #
7 # For the similarity solution see equation (3.19) in
8 # [\ref PeglerListerWorster2012] = PLW2012:
9 # S. Pegler, J. Lister, and M. G. Worster, 2012. Release of a viscous
10 # power-law fluid over an inviscid ocean", J. Fluid Mech. 700, 63--76.
11
12 # notes:
13 # don't use -pik because it will do -kill_icebergs and the whole thing is an iceberg
14 # could use -ssa_view_nuh if desired
15
16 # FIXME: need to carefully set hardness
17
18 from pylab import *
19 import sys
20 import time
21 # try different netCDF modules
22 try:
23 from netCDF4 import Dataset as CDF
24 except:
25 print("netCDF4 is not installed!")
26 sys.exit(1)
27
28 from optparse import OptionParser
29
30 parser = OptionParser()
31 parser.usage = \
32 """%prog MX DURATION
33
34 where MX is number of grid points,
35 DURATION is time in years for run,
36
37 Example: Try this diagnostic-only run:
38 $ export MX=101 YY=0
39 $ ./exactQ.py $MX $YY
40 $ pismr -o outQ$MX.nc -y $YY -i initQ$MX.nc -bootstrap -Mx $MX -My $MX -Mz 21 -Lz 1500 -z_spacing equal -surface given -stress_balance ssa -energy none -yield_stress constant -tauc 1e6 -ssa_dirichlet_bc -cfbc -part_grid -o_order zyx -ssa_e 1.0 -ssa_flow_law isothermal_glen
41 """
42 parser.description = "A script which runs Test Q."
43 (options, args) = parser.parse_args()
44 if (len(args) < 2) | (len(args) > 2):
45 print("ERROR; exactQ.py needs two arguments; run with --help to see usage")
46 print("... EXITING")
47 exit(-1)
48
49 SperA = 31556926.0
50
51 Mx = int(args[0])
52 My = Mx
53 runtime = float(args[1]) * SperA
54
55 ncfile = "initQ%d.nc" % Mx
56
57 # basic parameters
58 g = 9.81
59 rho = 910.0 # density of ice; kg/m^3
60 rhow = 1028.0 # density of ocean water; kg/m^3
61 n = 3.0
62 barB = 1.9e8 # strength of shelf; Pa s^(1/3); from MacAyeal et al 1996;
63 # FIXME is this equal to \tilde mu or not?
64
65 # derived parameters
66 m = (1.0 / n) - 1.0
67 gprime = (rhow - rho) * g / rhow # see just after (2.1) in PLW2012
68 nurescale = 3.0 ** (m / 2.0) * barB / rho # see just after (3.19) in PLW2012
69 C0 = 12.0 * nurescale / gprime # constant in (3.19) in PLW2012
70
71
72 def timeQ(H0):
73 # invert the first formula in (3.19) in PLW2012 to get t = f(H0)
74 t = (1.0 / (2.0 * n)) * (C0 / H0) ** n
75 return t
76
77
78 def geomvelQ(t, r, V):
79 # computes (3.19) in PLW2012, assuming t,V are scalar and r is array
80 # returns H,u with same size as r, but R is scalar
81 tnt = 2.0 * n * t
82 H = C0 * tnt ** (-1.0 / n) * ones(shape(r))
83 u = r / tnt
84 R = (V / (C0 * pi)) ** 0.5 * tnt ** (1.0 / (2.0 * n))
85 return H, u, R
86
87
88 # similarity solution: choose dimensions, get told the time and volume
89 H0 = 1000.0 # m
90 R0 = 100.0e3 # m; 100 km
91 V = pi * R0 ** 2 * H0
92 t0 = timeQ(H0)
93
94 t = t0 + runtime
95
96 print('exact Test Q has the following parameters for the start time t=t0:')
97 print(' time t0 = %.3e s = %f a' % (t0, t0 / SperA))
98 print(' thickness H0 = %.3f m' % H0)
99 print(' radius R0 = %.3f km' % (R0 / 1.0e3))
100 print('building PISM bootstrap file %s with Mx = %d and My = %d grid points ...' % (ncfile, Mx, My))
101
102 # set up the grid:
103 Lx = 200.0e3
104 Ly = Lx
105 x = linspace(-Lx, Lx, Mx)
106 y = linspace(-Ly, Ly, My)
107 [xx, yy] = meshgrid(x, y) # if there were "ndgrid" in numpy?
108 rr = sqrt(xx ** 2 + yy ** 2)
109
110 fill_value = nan
111
112 # create initial fields
113 topg = zeros((Mx, My)) - 1200.0 # so it is floating
114 ice_surface_temp = zeros((Mx, My)) + 263.15 # -10 degrees Celsius; just a guess
115 climatic_mass_balance = zeros((Mx, My))
116 thk = zeros((Mx, My))
117 thk[rr <= R0] = H0
118 zerossabc = zeros((Mx, My))
119
120 # create exact solution fields
121 thk_exact, c_exact, R_exact = geomvelQ(t, rr, V)
122 thk_exact[rr > R_exact] = 0.0
123 c_exact *= SperA
124 c_exact[rr > R_exact] = 0.0
125
126 print('exact Test Q at time t=%f years is in these variables:' % (t / SperA))
127 print(' c_exact, with max = %.3e' % c_exact.max())
128 print(' thk_exact, with max = %.3e' % thk_exact.max())
129 print('and R_exact = %.3f km' % (R_exact / 1.0e3))
130
131 # Write the data:
132 nc = CDF(ncfile, "w", format='NETCDF3_CLASSIC') # for netCDF4 module
133
134 # Create dimensions x and y
135 nc.createDimension("x", size=Mx)
136 nc.createDimension("y", size=My)
137
138 x_var = nc.createVariable("x", 'f4', dimensions=("x",))
139 x_var.units = "m"
140 x_var.long_name = "easting"
141 x_var.standard_name = "projection_x_coordinate"
142 x_var[:] = x
143
144 y_var = nc.createVariable("y", 'f4', dimensions=("y",))
145 y_var.units = "m"
146 y_var.long_name = "northing"
147 y_var.standard_name = "projection_y_coordinate"
148 y_var[:] = y
149
150
151 def def_var(nc, name, units, fillvalue):
152 # dimension transpose is standard: "float thk(y, x)" in NetCDF file
153 var = nc.createVariable(name, 'f', dimensions=("y", "x"), fill_value=fillvalue)
154 var.units = units
155 return var
156
157
158 bed_var = def_var(nc, "topg", "m", fill_value)
159 bed_var.standard_name = "bedrock_altitude"
160 bed_var[:] = topg
161
162 thk_var = def_var(nc, "thk", "m", fill_value)
163 thk_var.standard_name = "land_ice_thickness"
164 thk_var[:] = thk
165
166 climatic_mass_balance_var = def_var(nc, "climatic_mass_balance", "kg m-2 s-1", fill_value)
167 climatic_mass_balance_var.standard_name = "land_ice_surface_specific_mass_balance"
168 climatic_mass_balance_var[:] = climatic_mass_balance
169
170 ice_surface_temp_var = def_var(nc, "ice_surface_temp", "K", fill_value)
171 ice_surface_temp_var[:] = ice_surface_temp
172
173 u_ssa_bc_var = def_var(nc, "u_ssa_bc", "m s-1", fill_value)
174 u_ssa_bc_var[:] = zerossabc.copy()
175
176 v_ssa_bc_var = def_var(nc, "v_ssa_bc", "m s-1", fill_value)
177 v_ssa_bc_var[:] = zerossabc.copy()
178
179 bc_mask_var = nc.createVariable("bc_mask", "i", dimensions=("y", "x"))
180 bc_mask_var[:] = ((xx == 0.0) & (yy == 0.0))
181
182 thk_exact_var = def_var(nc, "thk_exact", "m", fill_value)
183 thk_exact_var[:] = thk_exact
184
185 c_exact_var = def_var(nc, "c_exact", "m year-1", fill_value)
186 c_exact_var[:] = c_exact
187
188 # set global attributes
189 nc.Conventions = 'CF-1.4'
190 historysep = ' '
191 historystr = time.asctime() + ': ' + historysep.join(sys.argv) + '\n'
192 setattr(nc, 'history', historystr)
193
194 nc.close()
195 print('file %s written ...' % ncfile)
196 print(' ... now run FIXME')