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()