tcolumn.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 --- tcolumn.py (17056B) --- 1 """Simple verification/regression tests for vertical diffusion and 2 advection of enthalpy in the column. Tests PISM's enthalpy solver. 3 4 """ 5 6 import PISM 7 from PISM.util import convert 8 import numpy as np 9 10 log10 = np.log10 11 12 ctx = PISM.Context() 13 unit_system = ctx.unit_system 14 config = ctx.config 15 16 config.set_string("grid.ice_vertical_spacing", "equal") 17 18 k = config.get_number("constants.ice.thermal_conductivity") 19 c = config.get_number("constants.ice.specific_heat_capacity") 20 rho = config.get_number("constants.ice.density") 21 K = k / c 22 # alpha squared 23 alpha2 = k / (c * rho) 24 25 EC = PISM.EnthalpyConverter(ctx.config) 26 pressure = np.vectorize(EC.pressure) 27 cts = np.vectorize(EC.enthalpy_cts) 28 29 30 class EnthalpyColumn(object): 31 "Set up the grid and arrays needed to run column solvers" 32 33 def __init__(self, Mz, dt): 34 self.Lz = 1000.0 35 self.z = np.linspace(0, self.Lz, Mz) 36 37 param = PISM.GridParameters() 38 param.Lx = 1e5 39 param.Ly = 1e5 40 param.z = PISM.DoubleVector(self.z) 41 param.Mx = 3 42 param.My = 3 43 param.Mz = Mz 44 45 param.ownership_ranges_from_options(1) 46 47 self.dt = dt 48 49 self.grid = PISM.IceGrid(ctx.ctx, param) 50 grid = self.grid 51 52 self.enthalpy = PISM.model.createEnthalpyVec(grid) 53 54 self.strain_heating = PISM.model.createStrainHeatingVec(grid) 55 56 self.u, self.v, self.w = PISM.model.create3DVelocityVecs(grid) 57 58 self.sys = PISM.enthSystemCtx(grid.z(), "energy.enthalpy", 59 grid.dx(), grid.dy(), self.dt, 60 config, 61 self.enthalpy, 62 self.u, self.v, self.w, 63 self.strain_heating, 64 EC) 65 66 # zero ice velocity: 67 self.reset_flow() 68 # no strain heating: 69 self.reset_strain_heating() 70 71 def reset_flow(self): 72 self.u.set(0.0) 73 self.v.set(0.0) 74 self.w.set(0.0) 75 76 def reset_strain_heating(self): 77 self.strain_heating.set(0.0) 78 79 def init_column(self): 80 ice_thickness = self.Lz 81 self.sys.init(1, 1, False, ice_thickness) 82 83 84 def diffusion_convergence_rate_time(title, error_func, plot=False): 85 "Compute the convergence rate with refinement in time." 86 dts = 2.0**np.arange(10) 87 88 max_errors = np.zeros_like(dts) 89 avg_errors = np.zeros_like(dts) 90 for k, dt in enumerate(dts): 91 max_errors[k], avg_errors[k] = error_func(False, dt_years=dt, Mz=101) 92 93 p_max = np.polyfit(log10(dts), log10(max_errors), 1) 94 p_avg = np.polyfit(log10(dts), log10(avg_errors), 1) 95 96 if plot: 97 plt.figure() 98 plt.title(title + "\nTesting convergence as dt -> 0") 99 log_plot(dts, max_errors, 'o', "max errors") 100 log_plot(dts, avg_errors, 'o', "avg errors") 101 log_fit_plot(dts, p_max, "max: dt^{:.2}".format(p_max[0])) 102 log_fit_plot(dts, p_avg, "avg: dt^{:.2}".format(p_avg[0])) 103 plt.axis('tight') 104 plt.grid(True) 105 plt.legend(loc="best") 106 107 return p_max[0], p_avg[0] 108 109 110 def diffusion_convergence_rate_space(title, error_func, plot=False): 111 "Compute the convergence rate with refinement in space." 112 Mz = np.array(2.0**np.arange(3, 10), dtype=int) 113 dzs = 1000.0 / Mz 114 115 max_errors = np.zeros_like(dzs) 116 avg_errors = np.zeros_like(dzs) 117 for k, M in enumerate(Mz): 118 T = 1.0 119 max_errors[k], avg_errors[k] = error_func(False, 120 T_final_years=T, 121 dt_years=T, 122 Mz=M) 123 124 p_max = np.polyfit(log10(dzs), log10(max_errors), 1) 125 p_avg = np.polyfit(log10(dzs), log10(avg_errors), 1) 126 127 if plot: 128 plt.figure() 129 plt.title(title + "\nTesting convergence as dz -> 0") 130 log_plot(dzs, max_errors, 'o', "max errors") 131 log_plot(dzs, avg_errors, 'o', "avg errors") 132 log_fit_plot(dzs, p_max, "max: dz^{:.2}".format(p_max[0])) 133 log_fit_plot(dzs, p_avg, "avg: dz^{:.2}".format(p_avg[0])) 134 plt.axis('tight') 135 plt.grid(True) 136 plt.legend(loc="best") 137 138 return p_max[0], p_avg[0] 139 140 141 def exact_DN(L, U0, QL): 142 """Exact solution (and an initial state) for the 'Dirichlet at the base, 143 Neumann at the top' setup.""" 144 n = 1 145 lambda_n = 1.0 / L * (np.pi / 2.0 + n * np.pi) 146 a = L * 25.0 147 148 def f(z, t): 149 v = a * np.exp(-lambda_n**2 * alpha2 * t) * np.sin(lambda_n * z) 150 return v + (U0 + QL * z) 151 return f 152 153 154 def errors_DN(plot_results=True, T_final_years=1000.0, dt_years=100, Mz=101): 155 """Test the enthalpy solver with Dirichlet B.C. at the base and 156 Neumann at the top surface. 157 """ 158 T_final = convert(T_final_years, "years", "seconds") 159 dt = convert(dt_years, "years", "seconds") 160 161 column = EnthalpyColumn(Mz, dt) 162 163 Lz = column.Lz 164 z = np.array(column.sys.z()) 165 166 E_base = EC.enthalpy(230.0, 0.0, EC.pressure(Lz)) 167 E_surface = EC.enthalpy(270.0, 0.0, 0.0) - 25*Lz 168 dE_surface = (E_surface - E_base) / Lz 169 170 E_steady = E_base + dE_surface * z 171 Q_surface = K * dE_surface 172 173 E_exact = exact_DN(Lz, E_base, dE_surface) 174 175 with PISM.vec.Access(nocomm=[column.enthalpy, 176 column.u, column.v, column.w, 177 column.strain_heating]): 178 column.sys.fine_to_coarse(E_exact(z, 0), 1, 1, column.enthalpy) 179 column.reset_flow() 180 column.reset_strain_heating() 181 182 t = 0.0 183 while t < T_final: 184 column.init_column() 185 186 column.sys.set_surface_heat_flux(Q_surface) 187 column.sys.set_basal_dirichlet_bc(E_base) 188 189 x = column.sys.solve() 190 191 column.sys.fine_to_coarse(x, 1, 1, column.enthalpy) 192 193 t += dt 194 195 E_exact_final = E_exact(z, t) 196 197 if plot_results: 198 t_years = convert(t, "seconds", "years") 199 200 plt.figure() 201 plt.xlabel("z, meters") 202 plt.ylabel("E, J/kg") 203 plt.step(z, E_exact(z, 0), color="blue", label="initial condition") 204 plt.step(z, E_exact_final, color="green", label="exact solution") 205 plt.step(z, cts(pressure(Lz - z)), "--", color="black", label="CTS") 206 plt.step(z, E_steady, "--", color="green", label="steady state profile") 207 plt.grid(True) 208 209 plt.step(z, x, label="T={} years".format(t_years), color="red") 210 211 plt.legend(loc="best") 212 213 errors = E_exact(z, t) - x 214 215 max_error = np.max(np.fabs(errors)) 216 avg_error = np.average(np.fabs(errors)) 217 218 return max_error, avg_error 219 220 221 def exact_ND(L, Q0, UL): 222 """Exact solution (and an initial state) for the 'Dirichlet at the base, 223 Neumann at the top' setup.""" 224 n = 2 225 lambda_n = 1.0 / L * (-np.pi / 2.0 + n * np.pi) 226 a = L * 25.0 227 228 def f(z, t): 229 v = a * np.exp(-lambda_n**2 * alpha2 * t) * np.sin(lambda_n * (L - z)) 230 return v + (UL + Q0 * (z - L)) 231 return f 232 233 234 def errors_ND(plot_results=True, T_final_years=1000.0, dt_years=100, Mz=101): 235 """Test the enthalpy solver with Neumann B.C. at the base and 236 Dirichlet B.C. at the top surface. 237 """ 238 T_final = convert(T_final_years, "years", "seconds") 239 dt = convert(dt_years, "years", "seconds") 240 241 column = EnthalpyColumn(Mz, dt) 242 243 Lz = column.Lz 244 z = np.array(column.sys.z()) 245 246 E_base = EC.enthalpy(240.0, 0.0, EC.pressure(Lz)) 247 E_surface = EC.enthalpy(260.0, 0.0, 0.0) 248 dE_base = (E_surface - E_base) / Lz 249 250 E_steady = E_surface + dE_base * (z - Lz) 251 Q_base = - K * dE_base 252 253 E_exact = exact_ND(Lz, dE_base, E_surface) 254 255 with PISM.vec.Access(nocomm=[column.enthalpy, 256 column.u, column.v, column.w, 257 column.strain_heating]): 258 column.sys.fine_to_coarse(E_exact(z, 0), 1, 1, column.enthalpy) 259 column.reset_flow() 260 column.reset_strain_heating() 261 262 t = 0.0 263 while t < T_final: 264 column.init_column() 265 266 column.sys.set_basal_heat_flux(Q_base) 267 column.sys.set_surface_dirichlet_bc(E_surface) 268 269 x = column.sys.solve() 270 271 column.sys.fine_to_coarse(x, 1, 1, column.enthalpy) 272 273 t += dt 274 275 E_exact_final = E_exact(z, t) 276 277 if plot_results: 278 t_years = convert(t, "seconds", "years") 279 280 plt.figure() 281 plt.xlabel("z, meters") 282 plt.ylabel("E, J/kg") 283 plt.step(z, E_exact(z, 0), color="blue", label="initial condition") 284 plt.step(z, E_exact_final, color="green", label="exact solution") 285 plt.step(z, cts(pressure(Lz - z)), "--", color="black", label="CTS") 286 plt.step(z, E_steady, "--", color="green", label="steady state profile") 287 plt.grid(True) 288 289 plt.step(z, x, label="T={} years".format(t_years), color="red") 290 291 plt.legend(loc="best") 292 293 errors = E_exact(z, t) - x 294 295 max_error = np.max(np.fabs(errors)) 296 avg_error = np.average(np.fabs(errors)) 297 298 return max_error, avg_error 299 300 301 def exact_advection(L, w): 302 "Exact solution of the 'pure advection' problem." 303 C = np.pi / L 304 305 def f(z, t): 306 return np.sin(C * (z - w * t)) 307 308 def df(z, t): 309 return C * np.cos(C * (z - w * t)) 310 return f, df 311 312 313 def errors_advection_up(plot_results=True, T_final=1000.0, dt=100, Mz=101): 314 """Test the enthalpy solver using a 'pure advection' problem with 315 Neumann (in-flow) B.C. at the base. 316 317 We use Dirichlet B.C. at the surface but they are irrelevant due 318 to upwinding. 319 """ 320 w = 1.0 321 322 config.set_number("constants.ice.thermal_conductivity", 0.0) 323 column = EnthalpyColumn(Mz, dt) 324 config.set_number("constants.ice.thermal_conductivity", k) 325 326 Lz = column.Lz 327 z = np.array(column.sys.z()) 328 329 E_exact, dE_exact = exact_advection(Lz, w) 330 331 with PISM.vec.Access(nocomm=[column.enthalpy, 332 column.u, column.v, column.w, 333 column.strain_heating]): 334 column.sys.fine_to_coarse(E_exact(z, 0), 1, 1, column.enthalpy) 335 column.reset_flow() 336 column.w.set(w) 337 column.reset_strain_heating() 338 339 t = 0.0 340 while t < T_final: 341 column.init_column() 342 343 column.sys.set_basal_neumann_bc(dE_exact(0, t+dt)) 344 column.sys.set_surface_dirichlet_bc(E_exact(Lz, t+dt)) 345 346 x = column.sys.solve() 347 348 column.sys.fine_to_coarse(x, 1, 1, column.enthalpy) 349 350 t += dt 351 352 if plot_results: 353 plt.figure() 354 plt.xlabel("z, meters") 355 plt.ylabel("E, J/kg") 356 plt.step(z, E_exact(z, 0), color="blue", label="initial condition") 357 plt.step(z, E_exact(z, t), color="green", label="exact solution") 358 plt.grid(True) 359 360 plt.step(z, x, label="T={} seconds".format(t), color="red") 361 362 plt.legend(loc="best") 363 364 errors = E_exact(z, t) - x 365 366 max_error = np.max(np.fabs(errors)) 367 avg_error = np.average(np.fabs(errors)) 368 369 return max_error, avg_error 370 371 372 def errors_advection_down(plot_results=True, T_final=1000.0, dt=100, Mz=101): 373 """Test the enthalpy solver using a 'pure advection' problem with 374 Neumann (in-flow) B.C. at the surface. 375 376 We use Dirichlet B.C. at the base but they are irrelevant due 377 to upwinding. 378 """ 379 w = -1.0 380 381 config.set_number("constants.ice.thermal_conductivity", 0.0) 382 column = EnthalpyColumn(Mz, dt) 383 config.set_number("constants.ice.thermal_conductivity", k) 384 385 Lz = column.Lz 386 z = np.array(column.sys.z()) 387 388 E_exact, dE_exact = exact_advection(Lz, w) 389 390 with PISM.vec.Access(nocomm=[column.enthalpy, 391 column.u, column.v, column.w, 392 column.strain_heating]): 393 column.sys.fine_to_coarse(E_exact(z, 0), 1, 1, column.enthalpy) 394 column.reset_flow() 395 column.w.set(w) 396 column.reset_strain_heating() 397 398 t = 0.0 399 while t < T_final: 400 column.init_column() 401 402 column.sys.set_basal_dirichlet_bc(E_exact(0, t+dt)) 403 column.sys.set_surface_neumann_bc(dE_exact(Lz, t+dt)) 404 405 x = column.sys.solve() 406 407 column.sys.fine_to_coarse(x, 1, 1, column.enthalpy) 408 409 t += dt 410 411 if plot_results: 412 plt.figure() 413 plt.xlabel("z, meters") 414 plt.ylabel("E, J/kg") 415 plt.step(z, E_exact(z, 0), color="blue", label="initial condition") 416 plt.step(z, E_exact(z, t), color="green", label="exact solution") 417 plt.grid(True) 418 419 plt.step(z, x, label="T={} seconds".format(t), color="red") 420 421 plt.legend(loc="best") 422 423 errors = E_exact(z, t) - x 424 425 max_error = np.max(np.fabs(errors)) 426 avg_error = np.average(np.fabs(errors)) 427 428 return max_error, avg_error 429 430 431 def advection_convergence_rate_time(title, error_func, plot=False): 432 "Compute the convergence rate with refinement in time." 433 dts = np.linspace(1, 101, 11) 434 Mz = 1001 435 T_final = 200 436 437 max_errors = np.zeros_like(dts) 438 avg_errors = np.zeros_like(dts) 439 for k, dt in enumerate(dts): 440 max_errors[k], avg_errors[k] = error_func(False, 441 T_final=T_final, 442 dt=dt, 443 Mz=Mz) 444 445 p_max = np.polyfit(log10(dts), log10(max_errors), 1) 446 p_avg = np.polyfit(log10(dts), log10(avg_errors), 1) 447 448 if plot: 449 plt.figure() 450 plt.title(title + "\nTesting convergence as dt -> 0") 451 log_plot(dts, max_errors, 'o', "max errors") 452 log_plot(dts, avg_errors, 'o', "avg errors") 453 log_fit_plot(dts, p_max, "max: dt^{:.2}".format(p_max[0])) 454 log_fit_plot(dts, p_avg, "avg: dt^{:.2}".format(p_avg[0])) 455 plt.axis('tight') 456 plt.grid(True) 457 plt.legend(loc="best") 458 459 return p_max[0], p_avg[0] 460 461 462 def advection_convergence_rate_space(title, error_func, plot=False): 463 "Compute the convergence rate with refinement in time." 464 dt = 0.4 465 Mzs = np.linspace(500, 5000, 11, dtype="i") 466 T_final = 10 467 468 dzs = 1000.0 / (Mzs - 1) 469 470 max_errors = np.zeros_like(dzs) 471 avg_errors = np.zeros_like(dzs) 472 for k, M in enumerate(Mzs): 473 max_errors[k], avg_errors[k] = error_func(False, 474 T_final=T_final, 475 dt=dt, 476 Mz=M) 477 478 p_max = np.polyfit(log10(dzs), log10(max_errors), 1) 479 p_avg = np.polyfit(log10(dzs), log10(avg_errors), 1) 480 481 if plot: 482 plt.figure() 483 plt.title(title + "\nTesting convergence as dz -> 0") 484 log_plot(dzs, max_errors, 'o', "max errors") 485 log_plot(dzs, avg_errors, 'o', "avg errors") 486 log_fit_plot(dzs, p_max, "max: dz^{:.2}".format(p_max[0])) 487 log_fit_plot(dzs, p_avg, "avg: dz^{:.2}".format(p_avg[0])) 488 plt.axis('tight') 489 plt.grid(True) 490 plt.legend(loc="best") 491 492 return p_max[0], p_avg[0] 493 494 495 def diffusion_DN_test(plot=False): 496 assert diffusion_convergence_rate_time("Diffusion: Dirichlet at the base, Neumann at the surface", 497 errors_DN, plot)[1] > 0.93 498 assert diffusion_convergence_rate_space("Diffusion: Dirichlet at the base, Neumann at the surface", 499 errors_DN, plot)[1] > 2.0 500 501 502 def diffusion_ND_test(plot=False): 503 assert diffusion_convergence_rate_time("Diffusion: Neumann at the base, Dirichlet at the surface", 504 errors_ND, plot)[1] > 0.93 505 assert diffusion_convergence_rate_space("Diffusion: Neumann at the base, Dirichlet at the surface", 506 errors_ND, plot)[1] > 2.0 507 508 509 def advection_up_test(plot=False): 510 assert advection_convergence_rate_time("Advection: Upward flow", 511 errors_advection_up, plot)[1] > 0.87 512 assert advection_convergence_rate_space("Advection: Upward flow", 513 errors_advection_up, plot)[1] > 0.96 514 515 516 def advection_down_test(plot=False): 517 assert advection_convergence_rate_time("Advection: Downward flow", 518 errors_advection_down, plot)[1] > 0.87 519 assert advection_convergence_rate_space("Advection: Downward flow", 520 errors_advection_down, plot)[1] > 0.96 521 522 523 if __name__ == "__main__": 524 import pylab as plt 525 526 def log_plot(x, y, style, label): 527 plt.plot(log10(x), log10(y), style, label=label) 528 plt.xticks(log10(x), ["{:.2f}".format(t) for t in x]) 529 plt.xlabel("log10(spacing)") 530 plt.ylabel("log10(error)") 531 532 def log_fit_plot(x, p, label): 533 plt.plot(log10(x), np.polyval(p, log10(x)), label=label) 534 535 diffusion_ND_test(plot=True) 536 diffusion_DN_test(plot=True) 537 advection_up_test(plot=True) 538 advection_down_test(plot=True) 539 plt.show()