Context Navigation

photochem_postproc.py @ 3537

Last change on this file since 3537 was 3528, checked in by mmaurice, 5 weeks ago

Generic PCM:

Python postprocessing: cleaner implementation, add automatic network
import and export from and to various format (Generic PCM, vulcan) and
more practical tools to handle networks. Updated the python notebook
with additional examples of use.

File size: 51.2 KB

Line
1	"""
2	Generic PCM Photochemistry post-processing library
3	Written by Maxime Maurice in 2024 anno domini
4	"""
5
6	import numpy as np
7	import matplotlib.pyplot as plt
8	import matplotlib.ticker as tk
9	import xarray as xr
10	from scipy.constants import R, N_A
11
12	import warnings
13	warnings.filterwarnings("ignore", message="The following kwargs were not used by contour: 't'")
14	warnings.filterwarnings("ignore", message="The following kwargs were not used by contour: 'lon'")
15	warnings.filterwarnings("ignore", message="The following kwargs were not used by contour: 'lat'")
16
17	background = 'co2' # background gas of the atmosphere
18
19	class GPCM_simu:
20	""" Generic PCM Simulation
21
22	Stores the netCDF file and path to the simulation
23
24	Attributes
25	----------
26	data : xr.Dataset
27	NetCDF file of the simulation
28	path : str
29	Path to simulation directory
30
31	Methods
32	-------
33	get_subset(field,**kw)
34	Get pa subset at fixed given coordinate of the data
35	plot_meridional_slice(field,logcb,labelcb,**kw)
36	Plot a meridional slice of a field
37	plot_time_evolution(field,logcb,labelcb,**kw)
38	Plot time evolution of a field (as a time vs altitude contour plot)
39	plot_time_series(field,lat,lon,alt,logy)
40	Plot time series of a field (at a specific location (lon, lat, alt), averged, or a combination thereof)
41	plot_atlas(field,t,alt,logcb,labelcb)
42	Plot atlas of a field
43	plot_profile(field,**kw)
44	Plot a profile of a field
45
46	"""
47	def __init__(self,path,filename='diagfi.c',verbose=False):
48	"""
49	Parameters
50	----------
51	path : str
52	Path to simulation directory
53	datafilename : str (optional)
54	Name of the netCDF file (by default: diagfi)
55	works with start, startfi, concat etc.
56	Do not add .nc type suffix
57	"""
58	self.path = path
59	try:
60	self.data = xr.open_dataset(path+'/'+filename,decode_times=False)
61	print(path+'/'+filename,'loaded, simulations lasts',self.data['Time'].values[-1],'sols')
62	except:
63	raise Exception('Data not found')
64
65	def __getitem__(self,field):
66	return self.data[field]
67
68	def __setitem__(self,field,value):
69	self.data[field] = value
70
71	def __area_weight__(self,data_array):
72	return self['aire']*data_array/self['aire'].mean('latitude')
73
74	def get_subset(self,field='all',**kw):
75	""" Get a subset at fixed given coordinate of the data
76
77	Can also average over a given dimension. In this case,
78	the meridional average is area-weighted.
79
80	Parameters
81	----------
82	field : str (optional, default = 'all')
83	Field name. If nothing or 'all'
84	specified, return all fields.
85	t : float (optional)
86	Time of the slice. If nothing
87	specified, use time average
88	lon : float (optional)
89	Longitude of the slice. If nothing
90	specified, use meridional average
91	lat : float (optional)
92	Latitude of the slice. If nothing
93	specified, use area-weighted zonal average
94	alt : float (optional)
95	Altitude of the slice. If nothing
96	specified, use time-average
97	"""
98	if field == 'all':
99	data_subset = self.data
100	else:
101	data_subset = self[field]
102
103	if 't' in kw and 'Time' in data_subset.dims:
104	if kw['t'] == 'avg':
105	data_subset = data_subset.mean(dim='Time')
106	else:
107	data_subset = data_subset.sel(Time=kw['t'],method='nearest')
108	if 'lat' in kw and 'latitude' in data_subset.dims:
109	if kw['lat'] == 'avg':
110	data_subset = self.__area_weight__(data_subset).mean(dim='latitude')
111	else:
112	data_subset = data_subset.sel(latitude=kw['lat'],method='nearest')
113	if 'lon' in kw and 'longitude' in data_subset.dims:
114	if kw['lon'] == 'avg':
115	data_subset = data_subset.mean(dim='longitude')
116	else:
117	data_subset = data_subset.sel(longitude=kw['lon'],method='nearest')
118	if 'alt' in kw and 'altitude' in data_subset.dims:
119	if kw['alt'] == 'avg':
120	data_subset = data_subset.mean(dim='altitude')
121	else:
122	data_subset = data_subset.sel(altitude=kw['alt'],method='nearest')
123
124	return data_subset
125
126	def plot_meridional_slice(self,field,t='avg',lon='avg',logcb=False,labelcb=None,**plt_kw):
127	""" Plot a meridional slice of a field
128
129	Parameters
130	----------
131	field : str
132	Field name to plot
133	logcb : bool (optional)
134	Use logarithmic colorscale
135	labelcb : str (optional)
136	Use custom colorbar label
137	t : float (keyword)
138	Time at which to plot (if nothing specified use time average)
139	lon : float (keyword)
140	Longitude at which to plot (if nothing specified use zonal average)
141	"""
142
143	if self['latitude'].size == 1:
144	# safety check
145	raise Exception('Trying to plot a meridional slice of a 1D simulation')
146
147	meridional_slice = self.get_subset(field,t=t,lon=lon)
148
149	if logcb:
150	plt.contourf(meridional_slice['latitude'],meridional_slice['altitude'],meridional_slice,locator=tk.LogLocator(),**plt_kw)
151	else:
152	plt.contourf(meridional_slice['latitude'],meridional_slice['altitude'],meridional_slice,**plt_kw)
153
154	plt.colorbar(label=field if labelcb==None else labelcb)
155	plt.xlabel('latitude [°]')
156	plt.ylabel('altitude [km]')
157
158	def plot_time_evolution(self,field,lat='avg',lon='avg',logcb=False,labelcb=None,**plt_kw):
159	""" Plot time evolution of a field (as a time vs altitude contour plot)
160
161	Parameters
162	----------
163	field : str
164	Field name to plot
165	lat : float (optional)
166	Latitude at which to plot (if nothing specified use area-weighted meridional average)
167	lon : float (optional)
168	Longitude at which to plot (if nothing specified use zonal average)
169	logcb : bool (optional)
170	Use logarithmic colorscale
171	labelcb : str (optional)
172	Use custom colorbar label
173	matplotlib contourf keyword arguments
174	"""
175
176	time_evolution = self.get_subset(field,lon=lon,lat=lat)
177
178	if logcb:
179	plt.contourf(time_evolution['Time'],time_evolution['altitude'],time_evolution.T,locator=tk.LogLocator(),**plt_kw)
180	else:
181	plt.contourf(time_evolution['Time'],time_evolution['altitude'],time_evolution.T,**plt_kw)
182
183	plt.colorbar(label=field if labelcb==None else labelcb)
184	plt.xlabel('time [day]')
185	plt.ylabel('altitude [km]')
186
187	def plot_time_series(self,field,lat='avg',lon='avg',alt='avg',logy=False,**plt_kw):
188	""" Plot time series of a field (at a specific location (lon, lat, alt), averged, or a combination thereof)
189
190	Parameters
191	----------
192	field : str
193	Field name to plot
194	lat : float (optional)
195	Latitude at which to plot (if nothing specified use area-weighted meridional average)
196	lon : float (optional)
197	Longitude at which to plot (if nothing specified use zonal average)
198	logy : bool (optional)
199	Use logarithmic y-axis
200	matplotlib plot keyword arguments
201	"""
202
203	time_series = self.get_subset(field,lon=lon,lat=lat,alt=alt)
204
205	if not 'label' in plt_kw:
206	plt_kw['label'] = self[field].units
207	if logy:
208	plt.semilogy(time_series['Time'],time_series,**plt_kw)
209	else:
210	plt.plot(time_series['Time'],time_series,**plt_kw)
211
212	plt.xlabel('time [day]')
213	plt.ylabel(field+' ['+self[field].units+']')
214
215	def plot_atlas(self,field,t='avg',alt='avg',logcb=False,labelcb=None,**plt_kw):
216	""" Plot atlas of a field
217
218	Parameters
219	----------
220	field : str
221	Field name to plot
222	t : float (optional)
223	Time at which to pot (if nothing specified, use time average)
224	alt : float (optional)
225	Altitude at which to plot (if nothing specified use vertical average)
226	logcb : bool (optional)
227	Use logarithmic colorscale
228	labelcb : str (optional)
229	Use custom colorbar label
230	matplotlib contourf keyword arguments
231	"""
232
233	if 'altitude' in self[field].dims:
234	atlas = self.get_subset(field,t=t,alt=alt)
235	else:
236	atlas = self.get_subset(field,t=t)
237
238	if logcb:
239	plt.contourf(atlas['longitude'],atlas['latitude'],atlas,locator=tk.LogLocator(),**plt_kw)
240	else:
241	plt.contourf(atlas['longitude'],atlas['latitude'],atlas,**plt_kw)
242
243	plt.colorbar(label=field if labelcb==None else labelcb)
244	plt.xlabel('longitude [°]')
245	plt.ylabel('matitude [°]')
246
247	def plot_profile(self,field,t='avg',lon='avg',lat='avg',logx=False,**plt_kw):
248	""" Plot a profile of a field
249
250	Parameters
251	----------
252	field : str
253	Field name to plot
254	logx : bool (optional)
255	Use logarithmic x axis
256	t : float (optional)
257	Time at which to select (if nothing specified use time-average)
258	lat : float (optional)
259	Latitude at which to plot (if nothing specified use area-weighted meridional average)
260	lon : float (optional)
261	Longitude at which to plot (if nothing specified use zonal average)
262	matplotlib's plot / semilogx keyword arguments
263	"""
264
265	profile = self.get_subset(field,t=t,lon=lon,lat=lat)
266
267	if logx:
268	plt.semilogx(profile,profile['altitude'],**plt_kw)
269	else:
270	plt.plot(profile,profile['altitude'],**plt_kw)
271
272	plt.xlabel(field+' ['+self[field].units+']')
273	plt.ylabel('altitude [km]')
274
275	def read_tracfile(self,filename='traceur.def'):
276	""" Read the traceurs of a simulation
277
278	Parameters
279	----------
280	filename : string (optional)
281	Name of the tracer file. Default: 'traceur.def'
282	"""
283	self.tracers = {}
284	self.M = {}
285	with open(self.path+'/'+filename) as tracfile:
286
287	for iline,line in enumerate(tracfile):
288
289	# First line
290	if iline == 0:
291	if not '#ModernTrac-v1' in line:
292	raise Exception('Can only read modern traceur.def')
293	continue
294
295	# Second line (number of tracers)
296	elif iline == 1:
297	continue
298
299	# Empty line
300	elif len(line.split()) == 0:
301	continue
302
303	# Commented line
304	elif line[0] == '!':
305	continue
306
307	# Regular entry
308	else:
309	line = line.replace('=',' ').split()
310	tracparams = {line[2i+1]:float(line[2i+2]) for i in range(int(len(line)/2))}
311	self.tracers = self.tracers \| {line[0]:tracparams}
312
313	def compute_rates(self):
314	""" Computes reaction rates for a simulation
315
316	Parameters
317	----------
318	s : GPCM_simu
319	Simulation object
320	reactions : dict (optional)
321	Dictionnary of reactions whose rate to compute as returned by read_reactfile
322	If nothing passed, call read_reactfile to identify reactions
323
324	Returns
325	-------
326	GPCM_simu
327	Simulation object with reactions rates, rates constants, species vmr and densities as new fields
328	"""
329
330	# self.network = network.from_file(self.path+'/chemnetwork/reactfile')
331	self.read_tracfile() # we need to read the traceur.def to know the molar mass of species
332	reactions = self.network.reactions
333	if not set(self.network.species) <= set(list(self.tracers.keys())):
334	raise Exception('Chemical network contains species that are not in the traceur.def file')
335
336	densities = {}
337
338	# Background density
339	self['total density'] = (self['p'] / R / self['temp'] / 1e6 * N_A).assign_attrs({'units':'cm^-3.s^-1'}) # 1e6 converts m³ to cm^3
340
341	for sp in self.network.species:
342	# volume mixing ratios
343	self[sp+' vmr'] = (self[sp] * self.tracers[background]['mmol'] / self.tracers[sp]['mmol']).assign_attrs({'units':'m^3/m^3'})
344	# molecular densities
345	self[sp+' density'] = (self[sp+' vmr'] * self['p'] / R / self['temp'] / 1e6 * N_A).assign_attrs({'units':'cm^-3.s^-1'}) # 1e6 converts m³ to cm^3
346	densities[sp] = self[sp+' density']
347
348	for r in reactions:
349
350	# Photolysis
351	if type(reactions[r]) == photolysis:
352
353	# Cases with branching ratios
354	if reactions[r].reactants[0] == 'co2':
355	if 'o1d' in reactions[r].products:
356	self['rate ('+reactions[r].formula+')'] = reactions[r].rate(densities,j=self['jco2_o1d'])
357	else:
358	self['rate ('+reactions[r].formula+')'] = reactions[r].rate(densities,j=self['jco2_o'])
359	elif reactions[r].reactants[0] == 'o2':
360	if 'o1d' in reactions[r].products:
361	self['rate ('+reactions[r].formula+')'] = reactions[r].rate(densities,j=self['jo2_o1d'])
362	else:
363	self['rate ('+reactions[r].formula+')'] = reactions[r].rate(densities,j=self['jo2_o'])
364	elif reactions[r].reactants[0] == 'o3':
365	if 'o1d' in reactions[r].products:
366	self['rate ('+reactions[r].formula+')'] = reactions[r].rate(densities,j=self['jo3_o1d'])
367	else:
368	self['rate ('+reactions[r].formula+')'] = reactions[r].rate(densities,j=self['jo3_o'])
369	elif reactions[r].reactants[0] == 'ch2o':
370	if 'cho' in reactions[r].products:
371	self['rate ('+reactions[r].formula+')'] = reactions[r].rate(densities,j=self['jch2o_cho'])
372	else:
373	self['rate ('+reactions[r].formula+')'] = reactions[r].rate(densities,j=self['jch2o_co'])
374	elif reactions[r].reactants[0] == 'h2o_vap':
375	self['rate ('+reactions[r].formula+')'] = reactions[r].rate(densities,j=self['jh2o'])
376	else:
377	# General case
378	self['rate ('+reactions[r].formula+')'] = reactions[r].rate(densities,j=self['j'+reactions[r].reactants[0]])
379	else:
380	self['k ('+reactions[r].formula+')'] = reactions[r].constant(self['temp'],densities[background])
381	self['rate ('+reactions[r].formula+')'] = reactions[r].rate(self['temp'],densities)
382
383	# 3-body reaction
384	if type(reactions[r]) == termolecular_reaction:
385	self['k ('+reactions[r].formula+')'] = self['k ('+reactions[r].formula+')'].assign_attrs({'units':'cm^6.s^-1'})
386
387	# 2-body reaction
388	else:
389	self['k ('+reactions[r].formula+')'] = self['k ('+reactions[r].formula+')'].assign_attrs({'units':'cm^3.s^-1'})
390
391	self['rate ('+reactions[r].formula+')'] = self['rate ('+reactions[r].formula+')'].assign_attrs({'units':'cm^-3.s^-1'})
392
393	def to_chempath(self,t,dt,filename_suffix='_chempath',lon='avg',lat='avg',alt='avg'):
394	""" Create files redable by the chemical path analyzer chempath (DOI: 10.5194/gmd-2024-163)
395
396	Parameters
397	----------
398	t : float
399	Initial time
400	dt : float
401	Timestep
402	filename_suffix : str (optional)
403	Suffix to append to the files being created (species.txt,
404	reactions.txt, model_time.dat, rates.dat, concentrations.dat)
405	lat : float (optional)
406	Latitude at which to plot (if nothing specified use area-weighted meridional average)
407	lon : float (optional)
408	Longitude at which to plot (if nothing specified use zonal average)
409	alt : float (optional)
410	Altitude of the slice. If nothing
411	matplotlib's plot / semilogx keyword arguments
412	"""
413
414	# Save pecies list in chempath format
415	with open(self.path+'/species'+filename_suffix+'.txt', 'w') as fp:
416	for sp in self.network.species:
417	fp.write("%s\n" % sp)
418
419	# Save reactions list in chempath format
420	with open(self.path+'/reactions'+filename_suffix+'.txt', 'w') as fp:
421	for r in self.network:
422	fp.write("%s\n" % r.to_string(format='chempath'))
423
424	# Save rates, concentrations and times in chempath format
425	rates = np.array([],dtype=np.float128)
426	conc = np.array([],dtype=np.float128)
427	conc_dt = np.array([],dtype=np.float128)
428	times = np.array([t,t+dt],dtype=np.float128)
429
430	out = self.get_subset(t=t,lon=lon,lat=lat,alt=alt)
431	for r in self.reactions:
432	rates = np.append(rates,out[f'rate ({r})'])
433	for sp in self.species:
434	conc = np.append(conc,out[f'{sp} vmr'])
435
436	out = self.get_subset(t=t+dt,lon=lon,lat=lat,alt=alt)
437	for sp in self.species:
438	conc_dt = np.append(conc_dt,out[f'{sp} vmr'])
439	conc = np.vstack((conc,conc_dt))
440
441	times.tofile(self.path+'/model_time'+filename_suffix+'.dat')
442	rates.tofile(self.path+'/rates'+filename_suffix+'.dat')
443	conc.tofile(self.path+'/concentrations'+filename_suffix+'.dat')
444
445	class reaction:
446	""" Instantiates a basic two-body reaction
447
448	Attributes
449	----------
450	formula : str
451	Reaction formula (e.g. "A + B -> C + D")
452	reactants : list
453	Reactanting molecules formulae (e.g. ["A", "B"])
454	products : list
455	Produced molecules formulae (e.g. ["C", "D"])
456	constant : callable
457	Reaction rate constant, function of potentially temperature and densities
458
459	Methods
460	-------
461	rate(T,densities)
462	Reaction rate for given temperature and densities
463	from_string(line,format)
464	Set up from an ASCII string
465	to_string(format)
466	Return an ASCII line readable by a photochemical model
467	"""
468	def __init__(self,reactants,products,constant):
469	"""
470	Parameters
471	----------
472	reactants : list(string)
473	Reacting molecules formulae (e.g. ["A", "B"])
474	products : list(string)
475	Produced molecules formulae (e.g. ["C", "D"])
476	constant : fun
477	Reaction rate constant, function of potentially temperature and densities
478	"""
479
480	self.formula = ''.join([r_+' + ' for r_ in reactants[:-1]])+reactants[-1]+' -> '+''.join([r_+' + ' for r_ in products[:-1]])+products[-1]
481	self.products = products
482	self.reactants = reactants
483	self.constant = constant
484
485	def rate(self,T,densities):
486	""" Computes reaction rate
487
488	Parameters
489	----------
490	T : float
491	Temperature [K]
492	densities : dict
493	Molecular densities [cm^-3]
494
495	Returns
496	-------
497	float
498	Value of the reaction rate [cm^-3.s^-1]
499	"""
500
501	return self.constant(T,densities[background])densities[self.reactants[0]]densities[self.reactants[1]]
502
503	@classmethod
504	def from_string(cls,line,format='GPCM',high_pressure_term=False):
505	""" Set up from an ASCII string
506
507	Format
508	------
509	GPCM
510	A B (... to col. 50) B C (... to col. 100) + cst string
511	vulcan
512	[ A + B -> C + D (... to col. 40) ] + cst string
513
514	Parameter
515	---------
516	line : string
517	ASCII string (usually formula and rate constant parameter)
518	format : string (optional)
519	Model format to write in (default: GPCM, options: GPCM, vulcan)
520	high_pressure_term : bool (optional)
521	Does the rate constant include a high-pressure term? (default: False)
522	"""
523	if format == 'GPCM':
524
525	reactants = line[:50].split()
526	products = line[50:100].split()
527	cst_params = line[101:]
528
529	if 'hv' in reactants:
530	reactants.pop(reactants.index('hv'))
531	if int(line[100]) == 0:
532	# photolysis calculated with cross sections
533	return cls(reactants,products)
534	else:
535	high_pressure_term = int(line[100]) == 2
536
537	elif format == 'vulcan':
538
539	reactants = line[line.index('[')+1:line.index('->')].split()[::2]
540	products = line[line.index('->')+2:line.index(']')].split()[::2]
541	cst_params = line[line.index(']')+1:]
542
543	if cst_params.split()[0][0].isalpha():
544	# photolysis calculated with cross sections
545	return cls(reactants,products,0.)
546
547	if 'M' in reactants:
548	reactants.pop(reactants.index('M'))
549	if 'M' in products:
550	products.pop(products.index('M'))
551
552	if high_pressure_term:
553	return cls(reactants,products,reaction_constant_dens_dep.from_string(cst_params,format))
554	else:
555	return cls(reactants,products,reaction_constant.from_string(cst_params,format))
556
557	def to_string(self,format='GPCM'):
558	""" Return an ASCII line readable by a photochemical model
559
560	Format
561	------
562	GPCM
563	A B (... to col. 50) B C (... to col. 100) + cst string
564	vulcan
565	[ A + B -> C + D (... to col. 40) ] + cst string
566	chempath
567	A+B=C+D
568
569	Parameter
570	---------
571	format : string (optional)
572	Model format to write in (default: GPCM, options: GPCM, vulcan, chempath)
573
574	Returns
575	-------
576	string
577	ASCII line readable by a photochemical model
578	"""
579
580	if format == 'GPCM':
581	line = ''
582	# reactants (characters 1 to 50)
583	for molecule in self.reactants:
584	line += molecule.lower().ljust(16,' ')
585	line = line.ljust(50,' ')
586
587	# products (characters 51 to 100)
588	for molecule in self.products:
589	line += molecule.lower().ljust(16,' ')
590	line = line.ljust(100,' ')
591
592	elif format == 'vulcan':
593	# formula
594	line = '[ '
595	for molecule in self.reactants:
596	line = line + molecule[:-4].upper() + ' + ' if '_vap' in molecule else line + molecule.upper() + ' + '
597	line = line[:-2] + '-> '
598	for molecule in self.products:
599	line = line + molecule[:-4].upper() + ' + ' if '_vap' in molecule else line + molecule.upper() + ' + '
600	line = line[:-2]
601	line = line.ljust(40,' ') + ' ] '
602
603	elif format == 'chempath':
604	line = r.formula
605	line = line.replace(' ','')
606	line = line.replace('->','=')
607	return line
608
609	# constant
610	if type(self.constant ) in [reaction_constant,reaction_constant_dens_dep]:
611	line += self.constant.to_string(format)
612	else:
613	print(self.formula,'has a custom reaction constant: add it manually')
614	return line
615
616
617	class termolecular_reaction(reaction):
618	""" Instantiates a three-body reaction
619
620	Attributes
621	----------
622	formula : str
623	Reaction formula (e.g. "A + B -> C + D")
624	reactants : list
625	Reactanting molecules formulae (e.g. ["A", "B"])
626	products : list
627	Produced molecules formulae (e.g. ["C", "D"])
628	constant : callable
629	Reaction rate constant, function of potentially temperature and densities
630
631	Methods
632	-------
633	rate(T,densities)
634	Reaction rate for given temperature and densities
635	from_string(line,format)
636	Set up from an ASCII string
637	to_string(format)
638	Return an ASCII line readable by a photochemical model
639	"""
640	def __init__(self,reactants,products,constant):
641
642	self.formula = ''.join([r_+' + ' for r_ in reactants[:-1]])+reactants[-1]+' + M -> '+''.join([r_+' + ' for r_ in products[:-1]])+products[-1]+' + M'
643	self.products = products
644	self.reactants = reactants
645	self.constant = constant
646
647	@classmethod
648	def from_string(cls,line,format='GPCM',high_pressure_term=False):
649	""" Set up from an ASCII string
650
651	Format
652	------
653	GPCM
654	A B (... to col. 50) B C (... to col. 100) + cst string
655	vulcan
656	[ A + B -> C + D (... to col. 40) ] + cst string
657
658	Parameter
659	---------
660	line : string
661	ASCII string (usually formula and rate constant parameter)
662	format : string (optional)
663	Model format to write in (default: GPCM, options: GPCM, vulcan)
664	high_pressure_term : bool (optional)
665	Does the rate constant include a high-pressure term? (default: False)
666	"""
667	new_instance = super().from_string(line,format,high_pressure_term)
668	if not high_pressure_term:
669	# In case we have a 3-body reaction without
670	# a high pressure term, enforce density dependence
671	new_instance.constant.params['d'] = 1.
672	return new_instance
673
674	def to_string(self,format='GPCM'):
675	""" Return an ASCII line readable by a photochemical model
676
677	Format
678	------
679	GPCM
680	A B M (... to col. 50) B C (... to col. 100) + cst string
681	vulcan
682	[ A + B + M -> C + D + M (... to col. 40) ] + cst string
683	chempath
684	A+B=C+D
685
686	Parameter
687	---------
688	format : string (optional)
689	Model format to write in (default: GPCM, options: GPCM, vulcan, chempath)
690
691	Returns
692	-------
693	string
694	ASCII line readable by a photochemical model
695	"""
696	if format == 'GPCM':
697	line = ''
698	# reactants (characters 1 to 50)
699	for molecule in self.reactants:
700	line += molecule.lower().ljust(16,' ')
701	line += 'M'
702	line = line.ljust(50,' ')
703
704	# products (characters 51 to 100)
705	for molecule in self.products:
706	line += molecule.lower().ljust(16,' ')
707	line = line.ljust(100,' ')
708
709	# constant (characters 101 to end of line)
710	line += self.constant.to_string(format)
711	return line
712
713	elif format == 'vulcan':
714	# formula
715	line = '[ '
716	for molecule in self.reactants:
717	line += molecule[:-4].upper() + ' + ' if '_vap' in molecule else molecule.upper() + ' + '
718	line += 'M -> '
719	for molecule in self.products:
720	line += molecule[:-4].upper() + ' + ' if '_vap' in molecule else molecule.upper() + ' + '
721	line += 'M'
722	line = line.ljust(40,' ') + ' ] '
723
724	# constant
725	line += self.constant.to_string(format)
726	return line
727
728	elif format == 'chempath':
729	line = r.formula
730	line = line.replace(' ','')
731	line = line.replace('->','=')
732	return line
733
734	class photolysis(reaction):
735
736	def __init__(self,reactants,products,constant=None):
737
738	self.formula = ''.join([r_+' + ' for r_ in reactants[:-1]])+reactants[-1]+' + hv -> '+''.join([r_+' + ' for r_ in products[:-1]])+products[-1]
739	self.products = products
740	self.reactants = reactants
741	self.constant = constant
742
743	def rate(self,densities,**kw):
744	""" Computes reaction rate
745
746	Parameters
747	----------
748	j : float
749	Photolysis rate [s^-1]
750	densities : dict
751	Molecular densities [cm^-3]
752
753	Returns
754	-------
755	float
756	Value of the reaction rate [cm^-3.s^-1]
757	"""
758
759	if 'j' in kw:
760	return kw['j']*densities[self.reactants[0]]
761	else:
762	# if a photolysis is prescribed with an Arrhenius constant, it is a trick to give
763	# it a constant rate. We can simply call the constant with an arbitrary temperature.
764	return self.constant(1.,densities[background])*densities[self.reactants[0]]
765
766	def to_string(self,format='GPCM'):
767	""" Return an ASCII line readable by a photochemical model
768
769	Format
770	------
771	GPCM
772	A hv (... to col. 50) B C (... to col. 100) + cst string
773	vulcan
774	[ A -> C + D (... to col. 40) ] A br (to add manually)
775	chemapth
776	A=C+D (to check)
777
778	Parameter
779	---------
780	format : string (optional)
781	Model format to write in (default: GPCM, options: GPCM, vulcan)
782
783	Returns
784	-------
785	string
786	ASCII line readable by a photochemical model
787	"""
788
789	if format == 'GPCM':
790	line = ''
791	# reactants (characters 1 to 50)
792	for molecule in self.reactants:
793	line = line + molecule.lower().ljust(16,' ')
794	line += 'hv'
795	line = line.ljust(50,' ')
796
797	# products (characters 51 to 100)
798	for molecule in self.products:
799	line += molecule.lower().ljust(16,' ')
800	line = line.ljust(100,' ')
801
802	# constant (characters 101 to end of line)
803	if self.constant == None:
804	line += '0'
805	print(self.formula,'is a photolysis: you need to add the cross section files manually')
806	else:
807	line += self.constant.to_string(format)
808	return line
809
810	elif format == 'vulcan':
811	# formula
812	line = '[ '
813	for molecule in self.reactants:
814	line += molecule[:-4].upper() + ' + ' if '_vap' in molecule else molecule.upper() + ' + '
815	line = line[:-2] + '->'
816	for molecule in self.products:
817	line += molecule[:-4].upper() + ' + ' if '_vap' in molecule else molecule.upper() + ' + '
818	line = line[:-2]
819	line = line.ljust(40,' ') + ' ] '
820	molecule = self.reactants[0]
821	line += molecule[:-4].upper() + ' + ' if '_vap' in molecule else molecule.upper()
822	print(self.formula,'is a photolysis: you need to add the branching ratio manually')
823	return line
824
825	elif format == 'chempath':
826	line = r.formula
827	line = line.replace(' ','')
828	line = line.replace('hv','')
829	line = line.replace('->','=')
830	return line
831
832	class reaction_constant:
833	""" Basic (Arrhenius) reaction rate constant
834
835	Instantiates type 1 rate constant for a particular reaction
836	(https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Photochemistry#Reaction_rate_formulae)
837
838	Attributes
839	----------
840	params : dict
841	Reaction-specific set of parameters for the rate constant: a,T0,c,b,d
842
843	Methods
844	-------
845	call(T,density)
846	Compute the reaction rate for given temperature and density
847	"""
848
849	def __init__(self,params):
850	"""
851	Parameters
852	----------
853	params : dict
854	Reaction-specific set of parameters for the rate constant: a,T0,c,b,d
855	"""
856	self.params = params
857
858	def __call__(self,T,density):
859	"""
860	Parameters
861	----------
862	T : float
863	Temperature [K]
864	density : float
865	Background gas density [cm^-3]
866
867	Returns
868	-------
869	float
870	Value of the reaction rate constant [cm^3.s^-1]
871	"""
872	return self.params['a'](T/self.params['T0'])self.params['c']np.exp(-self.params['b']/T)density*self.params['d']
873
874	@classmethod
875	def from_string(cls,line,format='GPCM'):
876	""" Creates an instance from an ASCII string in a variety of formats
877
878	Currently read formats: Generic PDM, vulcan
879
880	Parameters
881	----------
882	line : string
883	Rate constant parameters
884	format : string (optional)
885	Format in which parameters are writtenn (default: Generic PCM)
886
887	Returns
888	-------
889	reaction_constant
890	The instance of reaction_constant created
891
892	"""
893	if format == 'GPCM':
894	cst_param = line.split()
895	return cls({'a':float(cst_param[0]),'b':float(cst_param[1]),
896	'c':float(cst_param[2]),'T0':float(cst_param[3]),
897	'd':float(cst_param[4])})
898
899	elif format == 'vulcan':
900	cst_param = line.split()
901	T0 = 300.
902	return cls({'a':float(cst_param[0])T0*float(cst_param[1]),'b':float(cst_param[1]),
903	'c':float(cst_param[2]),'T0':T0,'d':0.})
904
905	def to_string(self,format='GPCM'):
906	""" Return an ASCII line readable by a photochemical model
907
908	Format
909	------
910	GPCM
911	1 a b c T0 d
912	vulcan
913	A B C
914
915	Parameter
916	---------
917	format : string (optional)
918	Model format to write in (default: GPCM, options: GPCM, vulcan)
919
920	Returns
921	-------
922	string
923	ASCII line readable by a photochemical model
924	"""
925
926	if format == 'GPCM':
927	return '1 '+'{:1.2e}'.format(self.params['a']).ljust(12,' ')+str(self.params['T0']).ljust(12,' ') \
928	+str(self.params['c']).ljust(12,' ')+str(self.params['b']).ljust(12,' ') \
929	+str(self.params['d'])
930	elif format == 'vulcan':
931
932	return '{:1.2e}'.format(self.params['a']/self.params['T0']**self.params['c']).ljust(12,' ') \
933	+str(self.params['c']).ljust(12,' ')+str(self.params['b'])
934
935	class reaction_constant_dens_dep(reaction_constant):
936	""" Type 2 reaction rate constant (Arrhenius with high pressure term)
937
938	Instantiates type 2 rate constant for a particular reaction
939	(https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Photochemistry#Reaction_rate_formulae)
940
941	Attributes
942	----------
943	params : dict
944	Reaction-specific set of parameters for the rate constant: k0,T0,n,a0,kinf,m,b0,g,h,dup,ddown,fc
945
946	Methods
947	-------
948	call(T,density)
949	Computes the reaction rate for given temperature and density
950	"""
951	def __call__(self,T,density):
952	"""
953	Parameters
954	----------
955	T : float
956	Temperature [K]
957	density : float
958	Background gas density [cm^-3]
959
960	Returns
961	-------
962	float
963	The value of the reaction rate constant [cm^3.s^-1]
964	"""
965	num = self.params['k0'](T/self.params['T0'])self.params['n']np.exp(-self.params['a0']/T)
966	den = self.params['kinf'](T/self.params['T0'])self.params['m']np.exp(-self.params['b0']/T)
967
968	return self.params['g']np.exp(-self.params['h']/T)+numdensity*self.params['dup']/(1+num/dendensity*self.params['ddown'])self.params['fc']*(1/(1+np.log10(num/dendensity)**2))
969
970	@classmethod
971	def from_string(cls,line,format='GPCM'):
972	""" Creates an instance from an ASCII string in a variety of formats
973
974	Currently read formats: Generic PDM, vulcan
975
976	Parameters
977	----------
978	line : string
979	Rate constant parameters
980	format : string (optional)
981	Format in which parameters are writtenn (default: Generic PCM)
982
983	Returns
984	-------
985	reaction_constant
986	The instance of reaction_constant created
987
988	"""
989	if format == 'GPCM':
990	cst_param = line.split()
991	return cls({'k0':float(cst_param[0]), 'n':float(cst_param[1]), 'a0':float(cst_param[2]),
992	'kinf':float(cst_param[3]),'m':float(cst_param[4]), 'b0':float(cst_param[5]),
993	'T0':float(cst_param[6]), 'fc':float(cst_param[7]), 'g':float(cst_param[8]),
994	'h':float(cst_param[9]), 'dup':float(cst_param[10]),'ddown':float(cst_param[10])})
995
996	elif format == 'vulcan':
997	cst_param = line.split()
998	T0 = 300.
999	return cls({'k0':float(cst_param[0])T0*float(cst_param[1]), 'n':float(cst_param[1]),'a0':float(cst_param[2]),
1000	'kinf':float(cst_param[3])T0*float(cst_param[4]),'m':float(cst_param[4]),'b0':float(cst_param[5]),
1001	'T0':T0,'fc':0.6,'g':0.,'h':0.,'dup':1,'ddown':1})
1002
1003	def to_string(self,format='GPCM'):
1004	""" Return an ASCII line readable by a photochemical model
1005
1006	Format
1007	------
1008	GPCM
1009	1 k0 n a0 kinf m b0 T0 fc g h dup ddown
1010	vulcan
1011	A_0 B_0 C_0 A_inf B_int C_inf
1012
1013	Parameter
1014	---------
1015	format : string (optional)
1016	Model format to write in (default: GPCM, options: GPCM, vulcan)
1017
1018	Returns
1019	-------
1020	string
1021	ASCII line readable by a photochemical model
1022	"""
1023	if format == 'GPCM':
1024	return '2 '+'{:1.2e}'.format(self.params['k0']).ljust(12,' ') +str(self.params['T0']).ljust(12,' ') +str(self.params['n']).ljust(12,' ')\
1025	+str(self.params['a0']).ljust(12,' ') +'{:1.2e}'.format(self.params['kinf']).ljust(12,' ')+str(self.params['m']).ljust(12,' ') \
1026	+str(self.params['b0']).ljust(12,' ') +str(self.params['g']).ljust(12,' ') +str(self.params['h']).ljust(12,' ') \
1027	+str(self.params['dup']).ljust(12,' ')+str(self.params['ddown']).ljust(12,' ')+str(self.params['fc'])
1028	elif format == 'vulcan':
1029	return '{:1.2e}'.format(self.params['k0']/self.params['T0']**self.params['n']).ljust(12,' ')+str(self.params['n']).ljust(12,' ')\
1030	+str(self.params['a0']).ljust(12,' ')+'{:1.2e}'.format(self.params['kinf']/self.params['T0']**self.params['m']).ljust(12,' ') \
1031	+str(self.params['m']).ljust(12,' ')+str(self.params['b0'])
1032
1033	class network:
1034	""" Reaction network object
1035
1036	Attributes
1037	----------
1038	reactions : dict
1039	Chemical reactions in the network
1040	species : list(string)
1041	Chemical species in the network
1042
1043	Methods
1044	-------
1045	append(to_append)
1046	Computes the reaction rate for given temperature and density
1047	update_species()
1048	Update list of species from list of reactions
1049	from_file(path,format)
1050	Instantiates a network from a file
1051	to_file(path,format)
1052	Save the network into a file
1053	get_subnetwork(criteria)
1054	Generate a subnetwork based on a dictionnary of given criteria
1055	"""
1056	def __init__(self,reactions={}):
1057
1058	self.reactions = reactions
1059	self.species = []
1060	if reactions != {}:
1061	self.update_species()
1062
1063	def __getitem__(self,formula):
1064
1065	return self.reactions[formula]
1066
1067	def __iter__(self):
1068
1069	self.current = list(self.reactions.keys())[0]
1070	return self
1071
1072	def __next__(self):
1073
1074	if self.current == 'finished':
1075	raise StopIteration
1076	elif self.current == list(self.reactions.keys())[-1]:
1077	current = self.current
1078	self.current = 'finished'
1079	else:
1080	current = self.current
1081	self.current = list(self.reactions.keys())[list(self.reactions.keys()).index(self.current)+1]
1082
1083	return self.reactions[current]
1084
1085	def append(self,to_append):
1086	""" Append a reaction to the network
1087
1088	Updates list of species to account for new species
1089
1090	Parameter
1091	---------
1092	to_append : reaction or network
1093	Reaction or network of reactions to append
1094	"""
1095	if type(to_append) == network:
1096	for r in to_append:
1097	self.append(r)
1098	else:
1099	self.reactions = self.reactions \| {to_append.formula:to_append}
1100
1101	self.update_species()
1102
1103	def update_species(self):
1104	""" Update list of species from list of reactions
1105
1106	"""
1107	if self.reactions == {}:
1108	raise Exception('Network empty')
1109
1110	for r in self:
1111	for sp in r.reactants:
1112	if not sp in self.species:
1113	self.species.append(sp)
1114	for sp in r.products:
1115	if not sp in self.species:
1116	self.species.append(sp)
1117
1118	@classmethod
1119	def from_file(cls,path,format='GPCM'):
1120	""" Instantiates a network from a file
1121
1122	Currently read formats: Generic PCM, VULCAN
1123
1124	Parameters
1125	----------
1126	path : str
1127	Path to the network file
1128	format : string (optional)
1129	Format of the network file to read. Default: GPCM
1130
1131	Returns
1132	-------
1133	network
1134	The network instance created
1135	"""
1136	reactions = {}
1137
1138	if format == 'GPCM':
1139	with open(path) as reactfile:
1140	for line in reactfile:
1141	# Commented line
1142	if line[0] == '!':
1143	# Hard-coded reaction
1144	if 'hard' in line and 'coded' in line:
1145	hard_coded_reaction = reaction(line[1:51].split(),line[51:101].split(),None)
1146	print('reaction ',hard_coded_reaction.formula,'seems to be hard-coded. Add it manually if needed.')
1147	continue
1148	else:
1149	# Photolysis
1150	if 'hv' in line:
1151	new_reaction = photolysis.from_string(line,format)
1152	# Other reactions
1153	else:
1154	# three-body reaction
1155	if 'M' in line:
1156	if line[line.index('M')+2] != ' ':
1157	# if third body is not the background gas, treat it as a simple reaction
1158	new_reaction = reaction.from_string(line.replace('M',' '),format)
1159	else:
1160	new_reaction = termolecular_reaction.from_string(line,format)
1161	# two-body reaction
1162	else:
1163	new_reaction = reaction.from_string(line,format)
1164	reactions[new_reaction.formula] = new_reaction
1165
1166	elif format == 'vulcan':
1167
1168	with open(path) as reactfile:
1169	re_tri = False
1170	re_tri_k0 = False
1171	special_re = False
1172	conden_re = False
1173	recomb_re = False
1174	photo_re = False
1175	ion_re = False
1176	for line in reactfile:
1177	if line.startswith("# 3-body"):
1178	re_tri = True
1179
1180	if line.startswith("# 3-body reactions without high-pressure rates"):
1181	re_tri_k0 = True
1182
1183	elif line.startswith("# special"):
1184	re_tri = False
1185	re_tri_k0 = False
1186	special_re = True # switch to reactions with special forms (hard coded)
1187
1188	elif line.startswith("# condensation"):
1189	re_tri = False
1190	re_tri_k0 = False
1191	special_re = False
1192	conden_re = True
1193
1194	elif line.startswith("# radiative"):
1195	re_tri = False
1196	re_tri_k0 = False
1197	special_re = False
1198	conden_re = False
1199	recomb_re = True
1200
1201	elif line.startswith("# photo"):
1202	re_tri = False
1203	re_tri_k0 = False
1204	special_re = False # turn off reading in the special form
1205	conden_re = False
1206	recomb_re = False
1207	photo_re = True
1208
1209	elif line.startswith("# ionisation"):
1210	re_tri = False
1211	re_tri_k0 = False
1212	special_re = False # turn off reading in the special form
1213	conden_re = False
1214	recomb_re = False
1215	photo_re = False
1216	ion_re = True
1217
1218	if line.startswith("#") or line.split() == []:
1219	continue
1220	else:
1221	line = line[line.index('['):]
1222
1223	if re_tri:
1224	new_reaction = termolecular_reaction.from_string(line,format,False)
1225	elif re_tri_k0:
1226	new_reaction = termolecular_reaction.from_string(line,format,True)
1227	elif special_re:
1228	print('special reaction :',line[line.index('[')+1:line.index(']')])
1229	elif conden_re:
1230	print('condensation reaction :',line[line.index('[')+1:line.index(']')])
1231	elif recomb_re:
1232	print('recombination reaction :',line[line.index('[')+1:line.index(']')])
1233	elif photo_re:
1234	new_reaction = photolysis.from_string(line,format)
1235	elif ion_re:
1236	print('ionisation reaction :',line[line.index('[')+1:line.index(']')])
1237	else:
1238	new_reaction = reaction.from_string(line,format)
1239
1240	reactions[new_reaction.formula] = new_reaction
1241
1242	return cls(reactions)
1243
1244	def to_file(self,path,format='GPCM'):
1245	""" Save the network into a file
1246
1247	Currently read formats: Generic PCM, VULCAN
1248
1249	Parameters
1250	----------
1251	path : str
1252	Path to the network file to create
1253	format : string (optional)
1254	Format of the network file to read. Default: GPCM
1255	"""
1256	if format == 'GPCM':
1257	with open(path, 'w') as reactfile:
1258	for i,r in enumerate(self):
1259	reactfile.write(f'# Reaction {str(i+1)}: {r.formula}\n')
1260	reactfile.write(r.to_string(format)+'\n')
1261
1262	elif format == 'vulcan':
1263	with open(path, 'w') as reactfile:
1264
1265	# header
1266	reactfile.write('# VULCAN photochemical network\n')
1267	reactfile.write('# Generated by the photochemical tools of the Generic PCM\n')
1268	reactfile.write('#########################################################\n')
1269	reactfile.write('# in the form of k = A T^B exp(-C/T)\n')
1270	# 2-body reactions
1271	reactfile.write('# Two-body Reactions\n')
1272	reactfile.write('# id Reactions A B C\n')
1273	reactfile.write('\n')
1274	n = 1
1275	for r in self.reactions:
1276	if type(self.reactions[r]) == reaction:
1277	reactfile.write(' '+str(n).ljust(7,' ')+self.reactions[r].to_string(format)+'\n')
1278	n += 1
1279
1280	# 3-body reactions with high-pressure term
1281	reactfile.write('\n')
1282	reactfile.write('# 3-body and Disscoiation Reactions\n')
1283	reactfile.write('# id # Reactions A_0 B_0 C_0 A_inf B_inf C_inf\n')
1284	reactfile.write('\n')
1285	for r in self.reactions:
1286	if type(self.reactions[r]) == termolecular_reaction and type(self.reactions[r].constant) == reaction_constant_dens_dep:
1287	reactfile.write(' '+str(n).ljust(7,' ')+self.reactions[r].to_string(format)+'\n')
1288	n += 1
1289
1290	# 3-body reactions without high-pressure term
1291	reactfile.write('\n')
1292	reactfile.write('# 3-body reactions without high-pressure rates\n')
1293	reactfile.write('# id # Reactions A_0 B_0 C_0\n')
1294	reactfile.write('\n')
1295	for r in self.reactions:
1296	if type(self.reactions[r]) == termolecular_reaction and type(self.reactions[r].constant) == reaction_constant:
1297	reactfile.write(' '+str(n).ljust(7,' ')+self.reactions[r].to_string(format)+'\n')
1298	n += 1
1299
1300	# photolysis
1301	reactfile.write('\n')
1302	reactfile.write('# reverse stops\n')
1303	reactfile.write('# photo disscoiation (no reversals) # use sp to link br_index to RXXX\n')
1304	reactfile.write('# id # Reactions sp br_index #(starting from 1)\n')
1305	reactfile.write('\n')
1306	for r in self.reactions:
1307	if type(self.reactions[r]) == photolysis:
1308	reactfile.write(' '+str(n).ljust(7,' ')+self.reactions[r].to_string(format)+'\n')
1309	n += 1
1310
1311	def save_traceur_file(self,path):
1312
1313	with open(path, 'w') as tracfile:
1314	tracfile.write('#ModernTrac-v1\n')
1315	tracfile.write(str(len(self.species))+'\n')
1316	for sp in self.species:
1317	tracfile.write(sp.ljust(24,' ')+'mmol='.ljust(10,' ')+'is_chim=1\n')
1318
1319	def get_subnetwork(self,criteria):
1320	""" Generate a subnetwork based on a dictionnary of given criteria
1321
1322	Parameter
1323	---------
1324	criteria : dict
1325	Selection criteria to include reactions in the subnetwork.
1326	Criteria are: - 'species' : list of any of the species from the network
1327	- 'element' : list of elements to include
1328	- 'type' : reaction, termolecular_reaction, photolysis
1329
1330	"""
1331	subnetwork = network()
1332
1333	for r in self.reactions:
1334
1335	keep = False
1336
1337	if 'type' in criteria:
1338	if type(self.reactions[r]) in criteria['type']:
1339	keep = True
1340
1341	if 'species' in criteria:
1342	for sp in criteria['species']:
1343	if sp in self.reactions[r].reactants + self.reactions[r].products:
1344	keep = True
1345
1346	if 'elements' in criteria:
1347	for elem in criteria['elements']:
1348	for sp in self.reactions[r].reactants + self.reactions[r].products:
1349	if elem in sp:
1350	keep = True
1351
1352	if keep:
1353	subnetwork.append(self.reactions[r])
1354
1355	return subnetwork

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: trunk/LMDZ.GENERIC/utilities/photochemistry/photochem_postproc.py @ 3537

Download in other formats: