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)