tbedrock_column.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 --- tbedrock_column.py (5160B) --- 1 #!/usr/bin/env python3 2 3 """Simple verification/regression tests for the heat equation in bedrock columns. Tests 4 PISM's "bedrock thermal unit". 5 6 """ 7 8 import PISM 9 from PISM.util import convert 10 import numpy as np 11 12 log10 = np.log10 13 14 ctx = PISM.Context() 15 16 k = ctx.config.get_number("energy.bedrock_thermal.conductivity") 17 c = ctx.config.get_number("energy.bedrock_thermal.specific_heat_capacity") 18 rho = ctx.config.get_number("energy.bedrock_thermal.density") 19 K = k / c 20 # alpha squared 21 alpha2 = k / (c * rho) 22 23 # set to True to plot solutions (for debugging) 24 plot_solutions = False 25 26 def convergence_rate_time(error_func, plot): 27 "Compute the convergence rate with refinement in time." 28 dts = 2.0**np.arange(4, 10) 29 30 max_errors = np.zeros_like(dts) 31 avg_errors = np.zeros_like(dts) 32 for k, dt in enumerate(dts): 33 max_errors[k], avg_errors[k] = error_func(plot_solutions, dt_years=dt, Mz=101) 34 35 p_max = np.polyfit(log10(dts), log10(max_errors), 1) 36 p_avg = np.polyfit(log10(dts), log10(avg_errors), 1) 37 38 if plot: 39 plt.figure() 40 plt.title("Heat equation in the bedrock column\nconvergence as dt -> 0") 41 log_plot(dts, max_errors, 'o', "max errors") 42 log_plot(dts, avg_errors, 'o', "avg errors") 43 log_fit_plot(dts, p_max, "max: dt^{:.3}".format(p_max[0])) 44 log_fit_plot(dts, p_avg, "avg: dt^{:.3}".format(p_avg[0])) 45 plt.axis('tight') 46 plt.grid(True) 47 plt.legend(loc="best") 48 49 return p_max[0], p_avg[0] 50 51 52 def convergence_rate_space(error_func, plot): 53 "Compute the convergence rate with refinement in space." 54 Mz = np.array(2.0**np.arange(2, 7), dtype=int) 55 dzs = 1000.0 / Mz 56 57 max_errors = np.zeros_like(dzs) 58 avg_errors = np.zeros_like(dzs) 59 for k, M in enumerate(Mz): 60 T = 1000.0 61 # time step has to be short enough so that errors due to the time discretization 62 # are smaller than errors due to the spatial discretization 63 dt = 0.001 * T 64 max_errors[k], avg_errors[k] = error_func(plot_solutions, 65 T_final_years=T, 66 dt_years=dt, 67 Mz=M) 68 69 p_max = np.polyfit(log10(dzs), log10(max_errors), 1) 70 p_avg = np.polyfit(log10(dzs), log10(avg_errors), 1) 71 72 if plot: 73 plt.figure() 74 plt.title("Heat equation in the bedrock column\nconvergence as dz -> 0") 75 log_plot(dzs, max_errors, 'o', "max errors") 76 log_plot(dzs, avg_errors, 'o', "avg errors") 77 log_fit_plot(dzs, p_max, "max: dz^{:.3}".format(p_max[0])) 78 log_fit_plot(dzs, p_avg, "avg: dz^{:.3}".format(p_avg[0])) 79 plt.axis('tight') 80 plt.grid(True) 81 plt.legend(loc="best") 82 83 return p_max[0], p_avg[0] 84 85 def exact(L, Q_bottom, U_top): 86 """Exact solution (and an initial state) for the 'Neumann at the base, 87 Dirichlet at the top' setup.""" 88 n = 2 89 lambda_n = 1.0 / L * (-np.pi / 2.0 + n * np.pi) 90 a = L * 25.0 91 92 def f(z, t): 93 v = a * np.exp(-lambda_n**2 * alpha2 * t) * np.sin(lambda_n * z) 94 return v + (U_top + Q_bottom * z) 95 return f 96 97 def errors(plot_results=True, T_final_years=1000.0, dt_years=100, Mz=101): 98 """Test the bedrock temperature solver with Neumann B.C. at the base and 99 Dirichlet B.C. at the top surface. 100 """ 101 T_final = convert(T_final_years, "years", "seconds") 102 dt = convert(dt_years, "years", "seconds") 103 104 Lz = 1000.0 105 dz = Lz / (Mz - 1.0) 106 107 column = PISM.BedrockColumn("btu", ctx.config, dz, int(Mz)) 108 109 z = np.linspace(-Lz, 0, Mz) 110 111 T_base = 240.0 # Kelvin 112 T_surface = 260.0 # Kelvin 113 dT_base = (T_surface - T_base) / Lz 114 115 T_steady = T_base + dT_base * (z - (-Lz)) 116 Q_base = - K * dT_base 117 118 T_exact = exact(Lz, dT_base, T_surface) 119 120 t = 0.0 121 # initial condition 122 x = T_exact(z, t) 123 while t < T_final: 124 x = column.solve(dt, Q_base, T_surface, x) 125 t += dt 126 127 T_exact_final = T_exact(z, t) 128 129 if plot_results: 130 t_years = convert(t, "seconds", "years") 131 132 plt.figure() 133 plt.xlabel("z, meters") 134 plt.ylabel("T, K") 135 plt.step(z, T_exact(z, 0), color="blue", label="initial condition") 136 plt.step(z, T_exact_final, color="green", label="exact solution") 137 plt.step(z, T_steady, "--", color="black", label="steady state profile") 138 plt.grid(True) 139 140 plt.step(z, x, label="T={} years".format(t_years), color="red") 141 142 plt.legend(loc="best") 143 144 errors = T_exact(z, t) - x 145 146 max_error = np.max(np.fabs(errors)) 147 avg_error = np.average(np.fabs(errors)) 148 149 return max_error, avg_error 150 151 def test(plot=False): 152 assert convergence_rate_time(errors, plot)[1] > 0.94 153 assert convergence_rate_space(errors, plot)[1] > 1.89 154 155 if __name__ == "__main__": 156 import pylab as plt 157 158 def log_plot(x, y, style, label): 159 plt.plot(log10(x), log10(y), style, label=label) 160 plt.xticks(log10(x), x) 161 162 def log_fit_plot(x, p, label): 163 plt.plot(log10(x), np.polyval(p, log10(x)), label=label) 164 165 test(plot=True) 166 plt.show()