resolved bugs in plotter_pp

minor bugs inside plotter_pp that were stopping plotting of particle surface areas have been resolved

Author: simonom

PyCHAM/__pycache__/accom_coeff_calc.cpython-39.pyc

@@ -24,7 +24,7 @@
 # changes due to gas-phase photochemistry and partitioning are included; 
 # generated in init_conc and treats loss from gas-phase as negative
 
-# File Created at 2024-11-12 15:18:27.878545
+# File Created at 2024-11-12 16:11:40.155400
 
 import numpy as np 
 

PyCHAM/__pycache__/dydt_rec.cpython-39.pyc

@@ -1 +1 @@
-/Users/user/Documents/GitHub/PyCHAM/PyCHAM/output/autoAPRAM_TOY_20241104_scheme/TOY_course_03
+/Users/user/Library/CloudStorage/OneDrive-TheUniversityofManchester/Man_Teach/research_experience/SAPHIR-STAR/PyCHAM_output/mcmAP_scheme/Baker_2024_Exp2p1

PyCHAM/__pycache__/hyst_eq.cpython-39.pyc

@@ -21,7 +21,7 @@
 ##########################################################################################
 '''solution of deliquescence and efflorescence RH, generated by eqn_pars.py in fully functioning mode, or by ui_check.py in testing mode'''
 # module to estimate deliquescence and efflorescence relative humidity as a function of temperature
-# File Created at 2024-11-12 16:00:59.660521
+# File Created at 2024-11-12 16:26:35.051795
 
 # function for deliquescence
 def drh(TEMP):

PyCHAM/__pycache__/ode_solv.cpython-39.pyc

@@ -21,7 +21,7 @@
 ##########################################################################################
 '''solution of ODEs, generated by eqn_pars.py'''
 # module to solve system of ordinary differential equations (ODEs) using solve_ivp of Scipy 
-# File Created at 2024-11-12 15:18:27.869878
+# File Created at 2024-11-12 16:11:40.109821
 
 import numpy as np
 import scipy.sparse as SP
@@ -156,6 +156,56 @@ def dydt(t, y): # define the ODE(s)
 		df_indx = (df_indx == 1) # transform to Boolean array 
 		dd[df_indx, 0] -= y[df_indx, 0]*self.dil_fac_now
 
+		# gas-particle and gas-wall partitioning-----------------
+		# transform component concentrations in particles and walls
+		# into size bins in rows, components in columns
+		ymat = (y[num_comp::, 0]).reshape(num_sb, num_comp)
+		
+		# for particles, force all components in bins with no particle to zero
+		ymat[0:num_asb, :][N_perbin[:, 0] == 0, :] = 0
+		
+		# for particles, calculate total particle-phase concentration per size bin (# molecules/cm3 (air))
+		csum = ((ymat[0:num_asb, :].sum(axis=1)-ymat[0:num_asb, self.seedi].sum(axis=1))+((ymat[0:num_asb, self.seedi]*core_diss).sum(axis=1)).reshape(-1)).reshape(-1, 1)
+		# tile total particle-phase concentration over components (# molecules/cm3 (air))
+		csum = np.tile(csum, [1, num_comp])
+		# concatenate wall bin total concentrations to total particle-phase concentration (# molecules/cm3)
+		csum = np.concatenate((csum, self.Cw), axis=0)
+		
+		# size bins with contents
+		isb = (csum[0:num_asb, 0] > 0.)
+		
+		# wall bins with contents
+		wsb = (self.Cw[:, 0] > 0.)
+		
+		# container for gas-phase concentrations at particle surface and at wall surface
+		Csit = np.zeros((num_sb, num_comp))
+		
+		# mole fractions of components at particle surface
+		Csit[0:num_asb, :][isb, :] = (ymat[0:num_asb, :][isb, :]/csum[0:num_asb, :][isb, :])
+		# mole fraction of components on walls, note that Cw included in csum above
+		Csit[num_asb::, :][wsb, :] = (ymat[num_asb::, :][wsb, :]/csum[num_asb::, :][wsb, :])
+		
+		# gas-phase concentration of components at
+		# particle surface (# molecules/cm3 (air))
+		Csit[0:num_asb, :][isb, :] = Csit[0:num_asb, :][isb, :]*self.Psat[0:num_asb, :][isb, :]*kelv_fac[isb]*act_coeff[0:num_asb, :][isb, :]
+		# partitioning rate (# molecules/cm3/s)
+		dd_all = kimt[0:num_asb, :]*(y[0:num_comp, 0].reshape(1, -1)-Csit[0:num_asb, :])
+		# gas-phase change
+		dd[0:num_comp, 0] -= dd_all.sum(axis=0)
+		# particle change
+		dd[num_comp:num_comp*(num_asb+1), 0] += (dd_all.flatten())
+		
+		if any(wsb):
+			# gas-phase concentration of components at
+			# wall surface (# molecules/cm3 (air))
+			Csit[num_asb::, :][wsb, :] = Csit[num_asb::, :][wsb, :]*self.Psat[num_asb::, :][wsb, :]*act_coeff[num_asb::, :][wsb, :]
+			# partitioning rate (# molecules/cm3/s)
+			dd_all = kimt[num_asb::, :]*(y[0:num_comp, 0].reshape(1, -1)-Csit[num_asb::, :])
+			# gas-phase change (summed over all wall bins)
+			dd[0:num_comp, 0] -= dd_all.sum(axis=0)
+			# wall change
+			dd[num_comp*(num_asb+1)::, 0] += (dd_all.flatten())
+		
 		dd = (dd[:, 0]).reshape(num_sb+1, num_comp)
 		# force all components in size bins with no particle to zero
 		if (num_asb > 0):
@@ -181,7 +231,7 @@ def jac(t, y): # define the Jacobian
 			y = y.reshape(-1, 1)
 
 		# elements of sparse Jacobian matrix
-		data = np.zeros((151))
+		data = np.zeros((4492+jac_mod_len))
 
 		for i in range(self.rindx_g.shape[0]): # gas-phase reaction loop
 			# reaction rate (# molecules/cm3/s)
@@ -195,6 +245,84 @@ def jac(t, y): # define the Jacobian
 			data[self.jac_indx_g[i, 0:self.njac_g[i, 0]]] += jac_coeff
 
 
+		# gas-particle partitioning
+		part_eff = np.zeros((1328))
+		if (sum(N_perbin[:, 0]) > 0.): # if any particles present 
+			part_eff[0:664:2] = -kimt[0:num_asb, :].sum(axis=0) # effect of gas on gas
+		
+		# empty array for any particle-on-gas and particle-on-particle effects on water in the particle-phase for rows of Jacobian
+		part_eff_rw = np.zeros((len(jac_part_hmf_indx)))
+		# empty array for any particle-on-gas and particle-on-particle effects of water in the particle-phase on non-water components in the particle-phase for columns of Jacobian
+		part_eff_cl = np.zeros((len(jac_part_H2O_indx)))
+		# starting index for jacobian row inputs for effect on water
+		sti_rw = 0 
+		
+		# transform particle phase concentrations into
+		# size bins in rows, components in columns
+		ymat = (y[num_comp:num_comp*(num_asb+1), 0]).reshape(num_asb, num_comp)
+		ymat[N_perbin[:, 0] == 0, :] = 0 # ensure zero components where zero particles
+		# total particle-phase concentration per size bin (molecules/cm3 (air))
+		csum = ymat.sum(axis=1)-ymat[:, self.seedi].sum(axis=1)+(ymat[:, self.seedi]*core_diss).sum(axis=1)
+		
+		# effect of particle on gas
+		for isb in range(int(num_asb)): # size bin loop
+			if (csum[isb] > 0): # if components present in this size bin
+				# effect of gas on particle
+				part_eff[1+isb:num_comp*(num_asb+1):num_asb+1] = +kimt[isb, :]
+				# start index
+				sti = int((num_asb+1)*num_comp+isb*(num_comp*2))
+				# diagonal index
+				diag_indxg = sti+np.arange(0, num_comp*2, 2).astype('int')
+				diag_indxp = sti+np.arange(1, num_comp*2, 2).astype('int')
+				# prepare for diagonal (component effect on itself)
+				diag = kimt[isb, :]*self.Psat[0, :]*act_coeff[0, :]*kelv_fac[isb, 0]*(-(csum[isb]-ymat[isb, :])/(csum[isb]**2.)) 
+				# implement to part_eff
+				part_eff[diag_indxg] -= diag
+				part_eff[diag_indxp] += diag
+				
+				if (rw_indx[isb] > -1): # if water in this size bin 
+					# prepare for row(s) (particle-phase non-water component effects on water in particle phase)
+					rw = kimt[isb, rw_indx[isb]]*self.Psat[0, rw_indx[isb]]*act_coeff[0, rw_indx[isb]]*kelv_fac[isb, 0]*(-(-ymat[isb, rw_indx[isb]])/(csum[isb]**2.)) 
+					# indices
+					indxg = sti_rw+np.arange(0, ((num_comp-1)*2), 2).astype('int')
+					indxp = sti_rw+np.arange(1, ((num_comp-1)*2), 2).astype('int')
+					# implement to part_eff_rw
+					part_eff_rw[indxg] -= rw
+					part_eff_rw[indxp] += rw
+					
+					# prepare for column(s) (particle-phase water effect on non-water in particle phase)
+					#cl = kimt[isb, :]*self.Psat[0, :]*act_coeff[0, :]*kelv_fac[isb, 0]*(-(-ymat[isb, :])/(csum[isb]**2.))
+					#cl = np.zeros((num_comp))
+					# remove water
+					#cl = np.concatenate((cl[0:H2Oi], cl[H2Oi+1::]))
+					#indxg = sti_rw+np.arange(0, (num_comp-1)).astype('int')
+					#indxp = sti_rw+np.arange((num_comp-1), (num_comp-1)*2).astype('int')
+					# implement to part_eff_cl
+					#part_eff_cl[indxg] -= cl
+					#part_eff_cl[indxp] += cl
+					
+					# starting index update
+					sti_rw += (num_comp-1)*2
+		
+		data[self.jac_part_indxn] += part_eff # diagonal
+		data[jac_part_hmf_indx] += part_eff_rw # rows
+		#data[jac_part_H2O_indx] += part_eff_cl # columns
+		
+		wsb = 0 # count on wall bins
+		# holder for wall effect
+		wall_eff = np.zeros((1328))
+		# effect of gas on gas 
+		wall_eff[0:664:2] = -np.sum(kimt[num_asb::, :], axis=0) 
+		for wsb in range(int(self.wall_on)): # wall bin loop
+			if (self.Cw[wsb, 0] > 0.):
+				# effect of gas on wall 
+				wall_eff[wsb+1:664:2] = +kimt[num_asb+wsb, :] 
+				# effect of wall on gas
+				wall_eff[wsb*2*num_comp+num_comp*(self.wall_on+1):num_comp*(self.wall_on+1)+(wsb+1)*2*num_comp:2] = +kimt[num_asb::, :][wsb, :]*(self.Psat[num_asb::, :][wsb, :]*act_coeff[num_asb::, :][wsb, :]/self.Cw[wsb, :]) 
+				# effect of wall on wall
+				wall_eff[wsb*2*num_comp+num_comp*(self.wall_on+1)+1:num_comp*(self.wall_on+1)+(wsb+1)*2*num_comp:2] = -kimt[num_asb::, :][wsb, :]*(self.Psat[num_asb::, :][wsb, :]*act_coeff[num_asb::, :][wsb, :]/self.Cw[wsb, :]) 
+		data[self.jac_wall_indxn] += wall_eff
+		
 		data[self.jac_extr_indx] -= 1.*self.dil_fac_now
 
 		# create Jacobian

PyCHAM/__pycache__/plotter_pp.cpython-39.pyc

@@ -67,6 +67,7 @@ def plotter(caller, dir_path, comp_names_to_plot, self):
 	rbou_rec = np.zeros((self.ro_obj.rad.shape[0], self.ro_obj.rad.shape[1]))
 	rbou_rec[:, :] = self.ro_obj.rad[:, :]
 	group_indx = self.ro_obj.gi
+	y_MV = np.array((self.ro_obj.comp_MV)) # cm3/mol
 
 	# number of actual particle size bins
 	num_asb = (self.ro_obj.nsb-self.ro_obj.wf)
@@ -296,7 +297,7 @@ def plotter(caller, dir_path, comp_names_to_plot, self):
 		ax0.xaxis.set_tick_params(labelsize = 14, direction = 'in')
 		ax0.legend(fontsize = 14)
 
-		# end of particle-phase concentration sub-plot ---------------------------------------
+		# end of particle-phase concentration sub-plot ----------------------------
 
 	# if called by button to plot temporal profile of total particle-phase concentration 
 	# excluding water and seed
@@ -314,7 +315,8 @@ def plotter(caller, dir_path, comp_names_to_plot, self):
 			ppc[:, i::num_comp] = 0.
 
 		# tile molar weights over size bins and times
-		y_mwt = np.tile(np.array((self.ro_obj.comp_MW)).reshape(1, -1), (1, self.ro_obj.nsb-self.ro_obj.wf))
+		y_mwt = np.tile(np.array((self.ro_obj.comp_MW)).reshape(1, -1), 
+			(1, self.ro_obj.nsb-self.ro_obj.wf))
 		y_mwt = np.tile(y_mwt, (ppc.shape[0], 1))
 		# convert from # molecules/cm3 to ug/m3
 		ppc = (ppc/si.N_A)*y_mwt*1.e12
@@ -323,7 +325,8 @@ def plotter(caller, dir_path, comp_names_to_plot, self):
 		ppc = np.sum(ppc, axis=1)
 
 		# plot
-		ax0.plot(self.ro_obj.thr, ppc, '+', linewidth = 4., label = 'total particle-phase excluding seed and water')
+		ax0.plot(self.ro_obj.thr, ppc, '+', linewidth = 4., 
+			label = 'total particle-phase excluding seed and water')
 		ax0.set_ylabel(r'Concentration ($\rm{\mu}$g$\,$m$\rm{^{-3}}$)', fontsize = 14)
 		ax0.set_xlabel(r'Time through simulation (hours)', fontsize = 14)
 		ax0.yaxis.set_tick_params(labelsize = 14, direction = 'in')
@@ -340,7 +343,8 @@ def plotter(caller, dir_path, comp_names_to_plot, self):
 			ppc = yrec[:, self.ro_obj.nc::]
 
 		# tile molar weights over size bins and times
-		y_mwt = np.tile(np.array((self.ro_obj.comp_MW)).reshape(1, -1), (1, (self.ro_obj.nsb-self.ro_obj.wf)))
+		y_mwt = np.tile(np.array((self.ro_obj.comp_MW)).reshape(1, -1), 
+			(1, (self.ro_obj.nsb-self.ro_obj.wf)))
 		y_mwt = np.tile(y_mwt, (ppc.shape[0], 1))
 
 		# convert to ug/m3 from # molecules/cm3
@@ -351,7 +355,8 @@ def plotter(caller, dir_path, comp_names_to_plot, self):
 
 		# sum particle-phase concentrations over size bins (ug/m3), 
 		# but keeping components separate
-		ppc_t = np.sum(ppc_t.reshape(self.ro_obj.nsb-self.ro_obj.wf, self.ro_obj.nc), axis=0)
+		ppc_t = np.sum(ppc_t.reshape(self.ro_obj.nsb-self.ro_obj.wf, 
+			self.ro_obj.nc), axis=0)
 
 		# convert to mass contributions
 		ppc_t = ((ppc_t/np.sum(ppc_t))*100.).reshape(-1, 1)
@@ -377,15 +382,18 @@ def plotter(caller, dir_path, comp_names_to_plot, self):
 			# get contribution (%)
 			ppci = (ppci[ppc_sbc>0.]/ppc_sbc[ppc_sbc>0.])*100.
 
-			ax0.plot(self.ro_obj.thr[ppc_sbc>0.], ppci, '-+', linewidth = 4., label = namei)
+			ax0.plot(self.ro_obj.thr[ppc_sbc>0.], ppci, '-+', linewidth = 4., 
+				label = namei)
 
-		ax0.set_ylabel(r'Contribution to particle-phase mass concentration (%)', fontsize = 14)
+		ax0.set_ylabel(r'Contribution to particle-phase mass concentration (%)', 
+			fontsize = 14)
 		ax0.set_xlabel(r'Time through simulation (hours)', fontsize = 14)
 		ax0.yaxis.set_tick_params(labelsize = 14, direction = 'in')
 		ax0.xaxis.set_tick_params(labelsize = 14, direction = 'in')
 		ax0.legend(fontsize = 14)
 
-	# if called by button to plot top contributors to particle-phase excluding seed and water
+	# if called by button to plot top contributors to particle-phase excluding 
+	# seed and water
 	if (caller == 7):
 
 		import scipy.constants as si		
@@ -402,7 +410,8 @@ def plotter(caller, dir_path, comp_names_to_plot, self):
 		ppc[:, self.ro_obj.H2O_ind::self.ro_obj.nc] = 0. # zero water
 
 		# tile molar weights over size bins and times
-		y_mwt = np.tile(np.array((self.ro_obj.comp_MW)).reshape(1, -1), (1, (self.ro_obj.nsb-self.ro_obj.wf)))
+		y_mwt = np.tile(np.array((self.ro_obj.comp_MW)).reshape(1, -1), 
+			(1, (self.ro_obj.nsb-self.ro_obj.wf)))
 		y_mwt = np.tile(y_mwt, (ppc.shape[0], 1))
 
 		# convert to ug/m3 from # molecules/cm3
@@ -439,7 +448,8 @@ def plotter(caller, dir_path, comp_names_to_plot, self):
 			# get contribution (%)
 			ppci = (ppci[ppc_sbc>0.]/ppc_sbc[ppc_sbc>0.])*100.
 
-			ax0.plot(self.ro_obj.thr[ppc_sbc>0.], ppci, '-+', linewidth = 4., label = namei)
+			ax0.plot(self.ro_obj.thr[ppc_sbc>0.], ppci, '-+', 
+				linewidth = 4., label = namei)
 
 		ax0.set_ylabel(r'Contribution to particle-phase mass concentration (%)', fontsize = 14)
 		ax0.set_xlabel(r'Time through simulation (hours)', fontsize = 14)
@@ -460,7 +470,8 @@ def plotter(caller, dir_path, comp_names_to_plot, self):
 		asp = np.sum(asp, axis=1)
 		# plot
 		ax0.plot(self.ro_obj.thr, asp, '-+', linewidth = 4.)
-		ax0.set_ylabel(r'Total particle-phase surface area concentration ($\rm{m^{2}\,m^{-3}}$)', fontsize = 14)
+		ax0.set_ylabel(str(r'Total particle-phase surface area concentration ' +
+		r'($\rm{m^{2}\,m^{-3}}$)'), fontsize = 14)
 		ax0.set_xlabel(r'Time through simulation (hours)', fontsize = 14)
 		ax0.yaxis.set_tick_params(labelsize = 14, direction = 'in')
 		ax0.xaxis.set_tick_params(labelsize = 14, direction = 'in')
@@ -480,9 +491,10 @@ def plotter(caller, dir_path, comp_names_to_plot, self):
 		ppcs = np.zeros((ppc.shape[0], (num_sb-wall_on)))
 
 		for i in seedi: # loop through seed indices
+
 			# get just seed component particle-phase volume concentration 
 			# (cm3/cm3)
-			ppcs[:, :] += (ppc[:, i::num_comp]/si.N_A)*(np.array((y_MV))[i])
+			ppcs[:, :] += (ppc[:, i::num_comp]/si.N_A)*y_MV[i]
 
 		# convert total volume to volume per particle (cm3)
 		ppcs[Nwet>0] = ppcs[Nwet>0]/Nwet[Nwet>0]

PyCHAM/__pycache__/rate_coeffs.cpython-39.pyc

@@ -21,7 +21,7 @@
 ##########################################################################################
 '''module for calculating reaction rate coefficients (automatically generated)'''
 # module to hold expressions for calculating rate coefficients # 
-# created at 2024-11-12 16:00:59.641214
+# created at 2024-11-12 16:26:35.038128
 
 import numpy
 import photolysisRates

PyCHAM/__pycache__/retr_out.cpython-39.pyc

@@ -45,7 +45,8 @@ def retr_out(self):
 	try: # try opening file
 		const_in = open(fname)
 	except:
-		err_mess = str('Error - no such file ' + fname + ', please check it still exists')
+		err_mess = str('Error - no such file ' + fname + 
+			', please check it still exists')
 		self.l203a.setText(err_mess)
 		# set border around error message
 		if (self.bd_pl == 1):

PyCHAM/dydt_rec.py

@@ -442,14 +442,16 @@ def saving(y_mat, Nresult_dry, Nresult_wet, t_out, num_comp,
 		Nresult_dry = (np.array((Nresult_dry)))
 		rbou_rec = (np.array((rbou_rec)))	
 
-	# to ensure quick use of output, store results in an object in an identical way to retr_out
+	# to ensure quick use of output, store results in an object in an identical 
+	# way to retr_out
 	# create a class to hold outputs
 	class ro_outputs:
 		sp = output_by_sim_sch_ext # chemical scheme path
 		vp = output_by_sim_mv_ext # model variables path
 		gi = group_indx # indices of groups of components
 		gen_numbers = self.gen_num # for each component, the generation number
-		HyC = self.HC.tolist() # hydrogen:carbon ratios for each component, this output added on 31/05/2022
+		# hydrogen:carbon ratios for each component, this output added on 31/05/2022
+		HyC = self.HC.tolist()
 		nominal_mass = self.nom_mass.tolist()
 		nsb = numsb
 		nc = num_comp

PyCHAM/err_log.txt

@@ -8,7 +8,7 @@
 
 setup(
     name="PyCHAM",
-    version="5.2.5",
+    version="5.2.6",
     author="Simon P. O'Meara, Shuxuan Xu and Ademipo Onanuga",
     author_email="simon.omeara@manchester.ac.uk",
     description="PyCHAM: CHemistry with Aerosol Microphysics in Python",