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