tmake_synth_ssa.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
---
tmake_synth_ssa.py (10023B)
---
1 #! /usr/bin/env python3
2 #
3 # Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2018, 2019, 2020 David Maxwell and Constantine Khroulev
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 sys, time, math
25
26 design_prior_scale = 0.2
27 design_prior_const = None
28
29
30 def groundedIceMisfitWeight(modeldata):
31 grid = modeldata.grid
32 weight = PISM.model.createVelocityMisfitWeightVec(grid)
33 mask = modeldata.vecs.mask
34 with PISM.vec.Access(comm=weight, nocomm=mask):
35 weight.set(0.)
36 grounded = PISM.MASK_GROUNDED
37 for (i, j) in grid.points():
38 if mask[i, j] == grounded:
39 weight[i, j] = 1
40 return weight
41
42
43 def fastIceMisfitWeight(modeldata, vel_ssa, threshold):
44 grid = modeldata.grid
45 weight = PISM.model.createVelocityMisfitWeightVec(grid)
46 mask = modeldata.vecs.ice_mask
47 threshold = threshold * threshold
48 with PISM.vec.Access(comm=weight, nocomm=[vel_ssa, mask]):
49 weight.set(0.)
50 grounded = PISM.MASK_GROUNDED
51 for (i, j) in grid.points():
52 u = vel_ssa[i, j].u
53 v = vel_ssa[i, j].v
54 if mask[i, j] == grounded:
55 if u * u + v * v > threshold:
56 weight[i, j] = 1
57 return weight
58
59
60 # The main code for a run follows:
61 if __name__ == '__main__':
62 context = PISM.Context()
63 com = context.com
64
65 PISM.set_abort_on_sigint(True)
66
67 PISM.verbPrintf(2, PISM.Context().com, "SSA forward model.\n")
68 if PISM.OptionBool("-version", "stop after printing PISM version"):
69 sys.exit(0)
70
71 usage = \
72 """ %s -i IN.nc -Mx number -My number [-o file.nc]
73 or (at python prompt)
74 run %s -i IN.nc -Mx number -My number [-o file.nc]
75 where:
76 -i IN.nc is input file in NetCDF format: contains PISM-written model state
77 -Mx number of grid points in the x direction
78 -My number of grid points in the y direction
79 notes:
80 * -i is required
81 """ % (sys.argv[0], sys.argv[0])
82
83 PISM.show_usage_check_req_opts(context.log, sys.argv[0], ["-i"], usage)
84
85 config = context.config
86 if not PISM.OptionString("-ssa_method", "").is_set():
87 config.set_string("stress_balance.ssa.method", "fem")
88
89 input_file_name = config.get_string("input.file")
90
91 config.set_string("output.file_name", "make_synth_ssa.nc", PISM.CONFIG_DEFAULT)
92 output_file_name = config.get_string("output.file_name")
93
94 design_prior_scale = PISM.OptionReal("-design_prior_scale",
95 "initial guess for design variable to be this factor of the true value",
96 design_prior_scale)
97
98 design_prior_const = PISM.OptionReal("-design_prior_const",
99 "initial guess for design variable to be this constant",
100 0.0)
101 design_prior_const = design_prior_const.value() if design_prior_const.is_set() else None
102
103 noise = PISM.OptionReal("-rms_noise", "pointwise rms noise to add (in m/a)", 0.0)
104 noise = noise.value() if noise.is_set() else None
105
106 misfit_weight_type = PISM.OptionKeyword("-misfit_type",
107 "Choice of misfit weight function",
108 "grounded,fast",
109 "grounded").value()
110
111 fast_ice_speed = PISM.OptionReal("-fast_ice_speed",
112 "Threshold in m/a for determining if ice is fast",
113 500.0)
114
115 generate_ssa_observed = PISM.OptionBool("-generate_ssa_observed",
116 "generate observed SSA velocities")
117
118 is_regional = PISM.OptionBool("-regional",
119 "Compute SIA/SSA using regional model semantics")
120
121 design_var = PISM.OptionKeyword("-inv_ssa",
122 "design variable for inversion",
123 "tauc,hardav",
124 "tauc").value()
125
126 ssa_run = PISM.ssa.SSAFromInputFile(input_file_name)
127
128 ssa_run.setup()
129
130 modeldata = ssa_run.modeldata
131 grid = modeldata.grid
132 vecs = modeldata.vecs
133
134 # add everything we need to "vecs" *before* it is locked in ssa_run.setup()
135 if design_var == 'tauc':
136 # Generate a prior guess for tauc
137 tauc_prior = PISM.model.createYieldStressVec(grid, name='tauc_prior',
138 desc="initial guess for (pseudo-plastic)"
139 " basal yield stress in an inversion")
140 vecs.add(tauc_prior, writing=True)
141 elif design_var == 'hardav':
142 # Generate a prior guess for hardav
143 vecs.add(PISM.model.createAveragedHardnessVec(grid))
144
145 hardav_prior = PISM.model.createAveragedHardnessVec(grid,
146 name='hardav_prior',
147 desc="initial guess for vertically averaged"
148 " ice hardness in an inversion")
149 vecs.add(hardav_prior, writing=True)
150
151 solve_t0 = time.time()
152 vel_ssa = ssa_run.solve()
153 solve_t = time.time() - solve_t0
154
155 PISM.verbPrintf(2, context.com, "Solve time %g seconds.\n", solve_t)
156
157 if design_var == 'tauc':
158 if design_prior_const is not None:
159 vecs.tauc_prior.set(design_prior_const)
160 else:
161 vecs.tauc_prior.copy_from(modeldata.vecs.tauc)
162 vecs.tauc_prior.scale(design_prior_scale)
163
164 tauc_true = modeldata.vecs.tauc
165 tauc_true.metadata(0).set_name('tauc_true')
166 tauc_true.set_attrs("diagnostic",
167 "value of basal yield stress used to generate synthetic SSA velocities",
168 "Pa", "Pa", "", 0)
169 vecs.markForWriting(tauc_true)
170 elif design_var == 'hardav':
171 # Generate a prior guess for hardav
172
173 EC = PISM.EnthalpyConverter(config)
174 ice_factory = PISM.IceFlowLawFactory(grid.com, "stress_balance.ssa.", config, EC)
175 ice_factory.removeType(PISM.ICE_GOLDSBY_KOHLSTEDT)
176 ice_factory.setType(config.get_string("stress_balance.ssa.flow_law"))
177 ice_factory.setFromOptions()
178 flow_law = ice_factory.create()
179 averaged_hardness_vec(flow_law, vecs.land_ice_thickness, vecs.enthalpy, vecs.hardav)
180
181 if design_prior_const is not None:
182 vecs.hardav_prior.set(design_prior_const)
183 else:
184 vecs.hardav_prior.copy_from(vecs.hardav)
185 vecs.hardav_prior.scale(hardav_prior_scale)
186
187 hardav_true = vecs.hardav
188 hardav_true.metadata(0).set_name('hardav_true')
189 hardav_true.set_attrs("diagnostic",
190 "vertically averaged ice hardness used to generate synthetic SSA velocities",
191 "Pa s^0.33333", "Pa s^0.33333", "", 0)
192 vecs.markForWriting(hardav_true)
193
194 vel_ssa_observed = vel_ssa # vel_ssa = ssa_run.solve() earlier
195
196 vel_ssa_observed.metadata(0).set_name("u_ssa_observed")
197 vel_ssa_observed.metadata(0).set_string("long_name", "x-component of 'observed' SSA velocities")
198
199 vel_ssa_observed.metadata(1).set_name("v_ssa_observed")
200 vel_ssa_observed.metadata(1).set_string("long_name", "y-component of 'observed' SSA velocities")
201
202 if generate_ssa_observed:
203 vecs.markForWriting(vel_ssa_observed)
204 final_velocity = vel_ssa_observed
205 else:
206 sia_solver = PISM.SIAFD
207 if is_regional:
208 sia_solver = PISM.SIAFD_Regional
209 vel_sia_observed = PISM.sia.computeSIASurfaceVelocities(modeldata, sia_solver)
210
211 vel_sia_observed.metadata(0).set_name('u_sia_observed')
212 vel_sia_observed.metadata(0).set_string('long_name', "x-component of the 'observed' SIA velocities")
213
214 vel_sia_observed.metadata(1).set_name('v_sia_observed')
215 vel_sia_observed.metadata(1).set_string('long_name', "y-component of the 'observed' SIA velocities")
216
217 vel_surface_observed = PISM.model.create2dVelocityVec(grid, "_surface_observed",
218 "observed surface velocities",
219 stencil_width=1)
220 vel_surface_observed.copy_from(vel_sia_observed)
221 vel_surface_observed.add(1., vel_ssa_observed)
222 vecs.markForWriting(vel_surface_observed)
223 final_velocity = vel_surface_observed
224
225 # Add the misfit weight.
226 if misfit_weight_type == "fast":
227 misfit_weight = fastIceMisfitWeight(modeldata, vel_ssa,
228 convert(fast_ice_speed, "m/year", "m/second"))
229 else:
230 misfit_weight = groundedIceMisfitWeight(modeldata)
231 modeldata.vecs.add(misfit_weight, writing=True)
232
233 if not noise is None:
234 u_noise = PISM.vec.randVectorV(grid, noise / math.sqrt(2), final_velocity.stencil_width())
235 final_velocity.add(convert(1.0, "m/year", "m/second"),
236 u_noise)
237
238 pio = PISM.util.prepare_output(output_file_name)
239 pio.close()
240
241 vecs.write(output_file_name)
242
243 # Save time & command line
244 PISM.util.writeProvenance(output_file_name)