Coverage for tcvx21/grillix_post/components/normalisation_m.py: 82%

1""" 

2Normalisation parameters, used for converting to SI units. 

3 

4Based on the pint unit-handling library. 

5 

6The typical usage is as follows; 

71. Use read_normalisation_file to read a formatted fortran namelist, typically called 'physical_parameters.nml' 

8 and return a dictionary called 'normalisation' 

92. Pass the 'normalisation' dictionary to dependent paramters 

10 i.e. rho_s0(**normalisation) 

11 which converts the dictionary into keyword arguments 

12 

13Alternatively, one can manually generate the normalisation dictionary 

14""" 

15from pathlib import Path 

16import numpy as np 

17import xarray as xr 

18from tcvx21.units_m import Quantity as Q_, unit_registry as ureg 

19from .namelist_reader_m import read_fortran_namelist 

20 

21 

22class Normalisation: 

23 """ 

24 dict-like class for converting from dimensionless arrays to SI units 

25 """ 

26 

27 @classmethod 

28 def initialise_from_normalisation_file(cls, normalisation_filepath: Path): 

29 """ 

30 Reads from the normalisation namelist, and returns a dictionary of quantities 

31 """ 

32 

33 experiment_parameters = dict(read_fortran_namelist(normalisation_filepath)[ 

34 'physical_parameters']) 

35 

36 physical_parameters = {} 

37 

38 for key in ['b0', 'te0', 'ti0', 'n0', 'r0', 'mi', 'z', 'z_eff']: 

39 assert key in experiment_parameters 

40 

41 # Magnetic field normalisation, usually taken on axis, in Tesla 

42 physical_parameters['B0'] = Q_(experiment_parameters['b0'], 'tesla') 

43 # Electron temperature normalisation, in electron-volts 

44 physical_parameters['Te0'] = Q_( 

45 experiment_parameters['te0'], 

46 'electron_volt') 

47 # Ion temperature normalisation, in electron-volts 

48 physical_parameters['Ti0'] = Q_( 

49 experiment_parameters['ti0'], 

50 'electron_volt') 

51 # Density normalisation, in particles-per-cubic-metres 

52 physical_parameters['n0'] = Q_(experiment_parameters['n0'], 'metre**-3') 

53 # Major radius, in metres (n.b. this is also the scale length for parallel 

54 # quantities) 

55 physical_parameters['R0'] = Q_(experiment_parameters['r0'], 'metre') 

56 # Ion mass, in amu 

57 physical_parameters['Mi'] = Q_( 

58 experiment_parameters['mi'], 

59 'proton_mass') 

60 # Ion charge 

61 physical_parameters['Z'] = experiment_parameters['z'] 

62 # Ion effective charge 

63 physical_parameters['Z_eff'] = experiment_parameters['z_eff'] 

64 

65 return cls(physical_parameters) 

66 

67 def __init__(self, normalisation_dictionary: dict): 

68 """ 

69 Each element of the input normalisation dictionary is copied as an attribute of the object. 

70 

71 i.e. If you define 

72 norm = Normalisation({'a':1, 'b':2}) 

73 then 

74 norm.a == 1 and norm.b == 2 

75 """ 

76 self.B0 = np.nan 

77 self.Te0 = np.nan 

78 self.Ti0 = np.nan 

79 self.n0 = np.nan 

80 self.R0 = np.nan 

81 self.Mi = np.nan 

82 self.Z = np.nan 

83 self.Z_eff = np.nan 

84 self.__dict__ = normalisation_dictionary 

85 

86 def __repr__(self): 

87 """ 

88 Returns a string representation for each element of the namelist 

89 """ 

90 string = f"\t{'Value':<30}\t{'Magnitude':<10}\t{'Units':<20}\n" 

91 

92 for parameter_name in sorted(dir(self), key=str.casefold): 

93 parameter = getattr(self, parameter_name, None) 

94 if isinstance(parameter, Q_): 

95 magnitude = parameter.magnitude 

96 units = str(parameter.units) 

97 

98 string += f"\t{parameter_name:<30}\t{magnitude:<10.4e}\t{units:<30}\n" 

99 

100 return string 

101 

102 def __getitem__(self, key): 

103 """ 

104 Allows normalisation to be indexed like a dictionary 

105 """ 

106 return getattr(self, key) 

107 

108 def convert_norm(self, input_array: xr.DataArray): 

109 """ 

110 Converts the 'norm' attribute of a xarray from a string to a Quantity 

111 """ 

112 assert isinstance(input_array, xr.DataArray), "convert_norm requires an xr.DataArray input" 

113 assert 'norm' in input_array.attrs, "convert_norm requires that the DataArray has an attribute 'norm'" 

114 

115 if isinstance(input_array.norm, str): 

116 input_array.attrs['norm'] = getattr(self, input_array.norm) 

117 

118 elif not isinstance(input_array.norm, Q_): 

119 raise NotImplementedError(f"Error: norm must be string or Quantity, but type was {type(input_array.norm)}") 

120 

121 return input_array 

122 

123 def as_dict(self): 

124 return_dict = {} 

125 for parameter_name in sorted(dir(self), key=str.casefold): 

126 parameter = getattr(self, parameter_name, None) 

127 if isinstance(parameter, Q_): 

128 return_dict[parameter_name] = parameter 

129 

130 return return_dict 

131 

132 elementary_charge = Q_(1, 'elementary_charge') 

133 electron_mass = Q_(1, 'electron_mass') 

134 proton_mass = Q_(1, 'proton_mass') 

135 atomic_mass_units = Q_(1, 'amu') 

136 speed_of_light = Q_(1, 'speed_of_light') 

137 vacuum_permeability = Q_(1, 'mu0') 

138 vacuum_permittivity = Q_(1, 'vacuum_permittivity') 

139 

140 @property 

141 def rho_s0(self): 

142 """Ion Larmor Radius [m] (n.b. this is also the scale length for perpendicular quantities)""" 

143 return (np.sqrt(self.Ti0 * self.Mi) / (self.elementary_charge * self.B0)).to_base_units() 

144 

145 @property 

146 def c_s0(self): 

147 """Sound speed [m/s]""" 

148 return (np.sqrt(self.Te0 / self.Mi)).to_base_units() 

149 

150 @property 

151 def v_A0(self): 

152 """Alfven speed [m/s]""" 

153 return (self.B0 / np.sqrt(self.vacuum_permeability * self.n0 * self.Mi)).to_base_units() 

154 

155 @property 

156 def tau_0(self): 

157 """Time normalisation [s]""" 

158 return (self.R0 / self.c_s0).to('s') 

159 

160 @staticmethod 

161 @ureg.wraps(ret="", args=['cm**-3', 'eV'], strict=True) 

162 def Coulomb_logarithm_ee(n, Te): 

163 """ 

164 Coulomb logarithm (unitless) for thermal electron-electron collisions 

165 

166 You can pass n and Te in any units, and these will be converted to 

167 cgs+eV such that the following formula applies 

168 """ 

169 if Te > 10: 

170 return 24. - np.log(np.sqrt(n) / Te) 

171 else: 

172 return 23. - np.log(np.sqrt(n) / Te**(1.5)) 

173 

174 

175 @property 

176 def lnLambda0(self): 

177 """ 

178 Thermal electron-electron Coulomb logarithm for reference parameters. 

179 

180 Note that we usually the Coulomb logarithm to be constant to ensure that the 

181 collision operator is bilinear 

182 """ 

183 return self.Coulomb_logarithm_ee(self.n0, self.Te0) 

184 

185 

186 def tau_ee(self, n, Te, Z): 

187 """ 

188 Electron-electron collision time 

189 

190 This is the mean time required for the direction of motion of an individual electron 

191 to be deflected by approximately 90 degrees due to collisions with electrons 

192 

193 Using 3.182 from Fitzpatrick (in SI units), but taking the parametric dependency 

194 from equation 2.5e from Braginskii (n = n_e = Z n_i, so n_i = n/Z.  

195 Therefore Z**2 * n_i = Z n) -- since Braginskii is c.g.s. and conversion to s.i. is a pain 

196 """ 

197 lnLambda = self.Coulomb_logarithm_ee(n, Te) 

198 

199 return ((6. * np.sqrt(2.) * np.pi**1.5 * self.vacuum_permittivity**2 * np.sqrt(self.electron_mass) * Te**1.5)/ 

200 (lnLambda * self.elementary_charge**4 * Z * n) 

201 ).to('s') 

202 

203 

204 def tau_ii(self, n, Te, Ti, Mi, Z): 

205 """ 

206 Ion-ion collision time 

207 

208 This is the mean time required for the direction of motion of an individual ion 

209 to be deflected by approximately 90 degrees due to collisions with ions 

210  

211 Using 3.184 from Fitzpatrick (in SI units), but taking the parametric dependency from 

212 equation 2.5i from Braginskii (n = n_e = Z n_i, so n_i = n/Z.  

213 Therefore Z**4 * n_i = Z**3 n) 

214 """ 

215 lnLambda = self.Coulomb_logarithm_ee(n, Te) 

216 

217 return ((12. * np.pi**1.5 * self.vacuum_permittivity**2 * np.sqrt(Mi) * Ti**1.5)/ 

218 (lnLambda * self.elementary_charge**4 * Z**3 * n) 

219 ).to('s') 

220 

221 

222 @property 

223 def tau_e(self): 

224 """ 

225 Reference electron-electron collision time, calculated with the effective charge Z_eff 

226 """ 

227 return self.tau_ee(self.n0, self.Te0, self.Z_eff) 

228 

229 

230 @property 

231 def tau_i(self): 

232 """ 

233 Reference ion-ion collision time 

234 """ 

235 return self.tau_ii(self.n0, self.Te0, self.Ti0, self.Mi, self.Z) 

236 

237 

238 @property 

239 def delta(self): 

240 """" 

241 Ratio of major radius to ion larmor radius. 

242 Conversion from perpendicular length scale to parallel length scale (>>1) 

243 """ 

244 return self.R0 / self.rho_s0 

245 

246 @property 

247 def rhostar(self): 

248 """ 

249 Ratio of ion larmor radius to major radius (inverse delta) 

250 Conversion from parallel length scale to perpendicular length scale (<<1) 

251 """ 

252 return self.rho_s0 / self.R0 

253 

254 @property 

255 def beta_0(self): 

256 """Electron dynamical beta""" 

257 return self.c_s0 ** 2 / (self.v_A0 ** 2) 

258 

259 @property 

260 def mu(self): 

261 """Electron to ion mass ratio""" 

262 return (self.electron_mass / self.Mi).to('') 

263 

264 @property 

265 def zeta(self): 

266 """Ti0/Te0: Ion to electron temperature ratio""" 

267 return self.Ti0 / self.Te0 

268 

269 @property 

270 def tau_e_norm(self): 

271 """Normalised electron collision time""" 

272 return self.tau_e * self.c_s0 / self.R0 

273 

274 @property 

275 def tau_i_norm(self): 

276 """Normalised ion collision time""" 

277 return self.tau_i * self.c_s0 / self.R0 

278 

279 @property 

280 def nu_e0(self): 

281 """Normalised electron collision frequency""" 

282 return self.tau_e_norm ** -1 

283 

284 @property 

285 def nu_i0(self): 

286 """Normalised ion collision frequency""" 

287 return self.tau_i_norm ** -1 

288 

289 @property 

290 def chipar_e(self): 

291 """Parallel electron heat conductivity""" 

292 return 3.16 * self.tau_e_norm * self.mu**-1 

293 

294 @property 

295 def chipar_i(self): 

296 """Parallel ion heat conductivity""" 

297 return 3.90 * self.zeta * self.tau_i_norm 

298 

299 @property 

300 def etapar(self): 

301 """Parallel resistivity""" 

302 return 0.51 * self.nu_e0 * self.mu 

303 

304 @property 

305 def ion_viscosity(self): 

306 """Ion viscosity coefficient""" 

307 return 0.96 * self.tau_i_norm