tconverter.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
---
tconverter.py (7947B)
---
1 import PISM
2 import numpy as np
3
4 config = PISM.Context().config
5
6 # list of converters
7 converters = {"Default": PISM.EnthalpyConverter(config),
8 "verification (cold)": PISM.ColdEnthalpyConverter(config)}
9
10
11 def try_all_converters(test):
12 print("")
13 for name, converter in list(converters.items()):
14 print("Testing '%s' converter..." % name)
15 test(name, converter)
16 print("done")
17
18
19 def reversibility_test():
20 "Converting from (E, P) to (T, omega, P)"
21
22 def run(name, EC):
23 # for a fixed pressure...
24 H = 1000.0
25 P = EC.pressure(H)
26
27 # cold ice
28 # enthalpy form
29 omega_prescribed = 0.0
30 E = EC.enthalpy(250.0, omega_prescribed, P)
31 # temperature form
32 T = EC.temperature(E, P)
33 omega = EC.water_fraction(E, P)
34 # we should get the same E back
35 assert E == EC.enthalpy(T, omega, P)
36 assert omega == omega_prescribed
37
38 # temperate ice
39 T_m = EC.melting_temperature(P)
40 # enthalpy form
41 omega_prescribed = 0.1
42 E = EC.enthalpy(T_m, omega_prescribed, P)
43 # temperature form
44 T = EC.temperature(E, P)
45 omega = EC.water_fraction(E, P)
46 # we should get the same E back
47 assert E == EC.enthalpy(T, omega, P)
48 assert np.fabs(omega - omega_prescribed) < 1e-16
49
50 try_all_converters(run)
51
52
53 def temperate_temperature_test():
54 "For temperate ice, an increase of E should not change T."
55
56 def run(name, EC):
57 H = 1000.0
58 P = EC.pressure(H)
59 T_m = EC.melting_temperature(P)
60 E = EC.enthalpy(T_m, 0.005, P)
61
62 if not EC.is_temperate(E, P):
63 # skip the test if our converter still treats this ice as
64 # cold
65 return
66
67 for delta_E in [0, 100, 1000]:
68 assert EC.temperature(E + delta_E, P) == T_m
69
70 try_all_converters(run)
71
72
73 def cts_computation_test():
74 "E_cts should be the same no matter how you compute it."
75
76 def run(name, EC):
77 H = 1000.0
78 P = EC.pressure(H)
79 T_m = EC.melting_temperature(P)
80
81 assert EC.enthalpy(T_m, 0.0, P) == EC.enthalpy_cts(P)
82
83 E_cts = EC.enthalpy_cts(P)
84 assert EC.pressure_adjusted_temperature(E_cts, P) == EC.melting_temperature(0)
85
86 try_all_converters(run)
87
88
89 def water_fraction_at_cts_test():
90 "Water fraction at CTS is zero."
91
92 def run(name, EC):
93 H = 1000.0
94 P = EC.pressure(H)
95 E = EC.enthalpy_cts(P)
96
97 assert EC.water_fraction(E, P) == 0
98
99 try_all_converters(run)
100
101
102 def enthalpy_of_water_test():
103 """Test the dependence of the enthalpy of water at T_m(p) on p."""
104
105 config = PISM.Context().config
106 c_w = config.get_number("constants.fresh_water.specific_heat_capacity")
107
108 EC = converters["Default"]
109
110 depth0 = 0.0
111 p0 = EC.pressure(depth0)
112 T0 = EC.melting_temperature(p0)
113 omega0 = 1.0
114 E0 = EC.enthalpy(T0, omega0, p0)
115
116 depth1 = 1000.0
117 p1 = EC.pressure(depth1)
118 T1 = EC.melting_temperature(p1)
119 omega1 = 1.0
120 E1 = EC.enthalpy(T1, omega1, p1)
121
122 # if we change the pressure of water from p0 to p1 while keeping
123 # it at T_m(p), its enthalpy should change and this change should
124 # be equal to c_w * (T_m(p1) - T_m(p0))
125 assert np.fabs((E1 - E0) - c_w * (T1 - T0)) < 1e-9
126
127
128 def invalid_inputs_test():
129 "Test that invalid inputs trigger errors."
130 def run(name, EC):
131 depth = 1000
132 pressure = EC.pressure(depth)
133 E_cts = EC.enthalpy_cts(pressure)
134 E_liquid = EC.enthalpy_liquid(pressure)
135 T_melting = EC.melting_temperature(pressure)
136
137 # don't test the converter that thinks this is cold:
138 if not EC.is_temperate(E_cts, pressure):
139 print("skipped...")
140 return
141
142 # temperature exceeds pressure melting
143 E = EC.enthalpy_permissive(T_melting + 1.0, 0.0, pressure)
144 # omega exceeds 1
145 E = EC.enthalpy_permissive(T_melting, 1.1, pressure)
146 # non-zero omega even though the ice is cold
147 E = EC.enthalpy_permissive(T_melting - 1.0, 0.1, pressure)
148 # negative omega
149 E = EC.enthalpy_permissive(T_melting, -0.1, pressure)
150
151 if not PISM.Pism_DEBUG:
152 # Skip remaining tests if PISM was built with Pism_DEBUG set to "OFF". (These
153 # checks are disabled in optimized builds, although we do ensure that
154 # EnthalpyConverter produces reasonable outputs even if it's fed garbage.)
155 return
156
157 try:
158 E = EC.temperature(E_liquid + 1.0, pressure)
159 raise AssertionError("failed to catch E > E_liquid in temperature()")
160 except RuntimeError:
161 pass
162
163 try:
164 E = EC.water_fraction(E_liquid + 1.0, pressure)
165 raise AssertionError("failed to catch E > E_liquid in water_fraction()")
166 except RuntimeError:
167 pass
168
169 try:
170 E = EC.enthalpy(T_melting + 1.0, 0.0, pressure)
171 raise AssertionError("failed to catch T > T_melting in enthalpy()")
172 except RuntimeError:
173 pass
174
175 try:
176 E = EC.enthalpy(T_melting, -0.1, pressure)
177 raise AssertionError("failed to catch omega < 0 in enthalpy()")
178 except RuntimeError:
179 pass
180
181 try:
182 E = EC.enthalpy(T_melting, 1.1, pressure)
183 raise AssertionError("failed to catch omega > 1 in enthalpy()")
184 except RuntimeError:
185 pass
186
187 try:
188 E = EC.enthalpy(-1.0, 0.0, pressure)
189 raise AssertionError("failed to catch T < 0 in enthalpy()")
190 except RuntimeError:
191 pass
192
193 try:
194 E = EC.enthalpy(T_melting - 1.0, 0.1, pressure)
195 raise AssertionError("failed to catch T < T_melting and omega > 0 in enthalpy()")
196 except RuntimeError:
197 pass
198
199
200 try_all_converters(run)
201
202
203 def plot_converter(name, EC):
204 """Test an enthalpy converter passed as the argument."""
205
206 H = 5000.0 # ice thickness
207
208 Z = np.linspace(0, H, int(H / 10)) # vertical levels
209
210 p = np.zeros_like(Z)
211 T_melting = np.zeros_like(Z)
212 E_cts = np.zeros_like(Z)
213 E = np.zeros_like(Z)
214 T = np.zeros_like(Z)
215 E_wet = np.zeros_like(Z)
216 omega = np.zeros_like(Z)
217
218 E[:] = 97000.0
219 delta_E = 3000.0
220
221 for i, z in enumerate(Z):
222 depth = H - z
223 p[i] = EC.pressure(depth)
224 T_melting[i] = EC.melting_temperature(p[i])
225 E_cts[i] = EC.enthalpy_cts(p[i])
226 T[i] = EC.temperature(E[i], p[i])
227 # dependence on pressure for high omega and T_melting
228 E_wet[i] = EC.enthalpy(T_melting[i], 0.95, p[i])
229 omega[i] = EC.water_fraction(E_cts[i] + delta_E, p[i]) * 100.0
230
231 plt.figure(figsize=(15, 8))
232 plt.subplot(2, 2, 1)
233 plt.title("%s enthalpy converter" % name)
234 plt.plot(Z, E, label="constant enthalpy", lw=2)
235 plt.plot(Z, E_cts, label="enthalpy corresponding to CTS", lw=2)
236 plt.legend(loc='best')
237 plt.ylabel("J/kg")
238 plt.grid()
239
240 plt.subplot(2, 2, 2)
241 plt.title("%s enthalpy converter" % name)
242 plt.plot(Z, omega, label="water fraction for E = E_cts + C", lw=2)
243 plt.legend(loc='best')
244 plt.ylabel("percent")
245 plt.grid()
246
247 plt.subplot(2, 2, 3)
248 plt.plot(Z, T_melting, label="melting temperature", lw=2)
249 plt.plot(Z, T, label="temperature corresponding to constant E", lw=2)
250 plt.legend(loc='best')
251 plt.ylabel("Kelvin")
252 plt.grid()
253 plt.xlabel("height above the base of the ice, m")
254
255 plt.subplot(2, 2, 4)
256 plt.plot(Z, E_wet, label="temperate ice enthalpy with high omega", lw=2)
257 plt.legend(loc='best')
258 plt.xlabel("height above the base of the ice, m")
259 plt.ylabel("J/kg")
260 plt.grid()
261
262 if __name__ == "__main__":
263
264 import pylab as plt
265
266 for name, converter in list(converters.items()):
267 plot_converter(name, converter)
268
269 plt.show()