Further towards merging 8 species and CR support

BOLSIGChemistry_NominalRates/detailed_balance.py

@@ -0,0 +1,135 @@
+import numpy as np
+import h5py as h5
+from matplotlib import pyplot as plt
+
+E_lvl_m = 11.548
+E_lvl_r = 11.624
+E_lvl_4p = 12.907
+E_lvl_i = 15.76
+g_m = 5.0
+g_r = 3.0
+g_4p = 3.0
+g_i = 4.0
+n_e = 3.7e15
+
+rates14 = []
+rates21 = []
+rates26 = []
+rates29 = []
+rates35 = []
+data10 = h5.File('./1s-metastable.h5', 'r')["table"]
+data11 = h5.File('./1s-resonance.h5', 'r')["table"]
+data12 = h5.File('./2p-lumped.h5', 'r')["table"]
+data13 = h5.File('./Ionization.h5', 'r')["table"]
+data15 = h5.File('./StepIonization.h5', 'r')["table"]
+
+Te = data10[:,0]
+rates10 = data10[:,1]
+rates11 = data11[:,1]
+rates12 = data12[:,1]
+rates13 = data13[:,1]
+rates15 = data15[:,1]
+
+for j in range(len(Te)):
+    rates14.append(rates10[j] * (1 / g_m) * np.exp(E_lvl_m / (Te[j] / 11604)))
+    rates21.append(rates11[j] * (1 / g_r) * np.exp(E_lvl_r / (Te[j] / 11604)))
+    rates26.append(rates15[j] * (g_m / g_i) * np.exp((E_lvl_i - E_lvl_m) / (Te[j] / 11604)) / n_e * 6.022e23)
+    rates29.append(rates13[j] * (1 / g_i) * np.exp(E_lvl_i / (Te[j] / 11604)) / n_e * 6.022e23)
+    rates35.append(rates12[j] * (1 / g_4p) * np.exp(E_lvl_4p / (Te[j] / 11604)))
+
+
+data14 = np.empty(shape = (len(Te), 2))
+data21 = np.empty(shape = (len(Te), 2))
+data26 = np.empty(shape = (len(Te), 2))
+data29 = np.empty(shape = (len(Te), 2))
+data35 = np.empty(shape = (len(Te), 2))
+
+data14[:,0] = Te
+data14[:,1] = rates14
+data21[:,0] = Te
+data21[:,1] = rates21
+data26[:,0] = Te
+data26[:,1] = rates26
+data29[:,0] = Te
+data29[:,1] = rates29
+data35[:,0] = Te
+data35[:,1] = rates35
+
+with h5.File('./Deexci-metastable.h5', 'w') as f:
+    dataset = f.create_dataset("table", data = data14)
+
+with h5.File('./Deexci-resonance.h5', 'w') as g:
+    dataset = g.create_dataset("table", data = data21)
+
+with h5.File('./3BdyRecomb-ground.h5', 'w') as h:
+    dataset = h.create_dataset("table", data = data29)
+
+with h5.File('./3BdyRecomb-metastable.h5', 'w') as x:
+    dataset = x.create_dataset("table", data = data26)
+
+with h5.File('./Deexci-2p.h5', 'w') as y:
+    dataset = y.create_dataset("table", data = data35)
+
+#rxn10 = h5.File("./BOLSIGChemistry_6SpeciesRates/lumped.metastable.h5", 'r')
+#rxn11 = h5.File("./BOLSIGChemistry_6SpeciesRates/lumped.resonance.h5", 'r')
+
+#data10 = rxn10["table"]
+#data11 = rxn11["table"]
+
+#Te = data10[:,0]
+#Te /= 11604
+#rates10 = data10[:,1]
+#rates11 = data11[:,1]
+#rates10 /= 6.022e23
+#rates11 /= 6.022e23
+
+#rates14 = []
+#rates21 = []
+
+#for i in range(len(rates10)):
+#    rates14.append(rates10[i] * (1 / g_m) * np.exp(E_lvl_m / Te[i]))
+#    rates21.append(rates11[i] * (1 / g_r) * np.exp(E_lvl_r / Te[i]))
+
+
+#arrh_14 = []
+#arrh_21 = []
+
+#for j in range(len(Te)):
+#    arrh_14.append(4.3e-16 * Te[j]**0.74)
+#    arrh_21.append(4.3e-16 * Te[j]**0.74)
+
+
+#fig,ax = plt.subplots()
+#ax.set_xscale('log')
+#ax.set_yscale('log')
+#ax.plot(Te, rates14, color = 'darkblue', label = 'CRSC - 14')
+#ax.plot(Te, arrh_14, color = 'cyan', label = 'Arrhenius - 14')
+#ax.plot(Te, rates21, color = 'darkred', label = 'CRSC - 21')
+#ax.plot(Te, arrh_21, color = 'orange', label = 'Arrhenius - 21')
+#ax.legend()
+#ax.set_xlabel('Te [eV]')
+#ax.set_ylabel('Rate [# / m3-s]')
+#plt.savefig('RatesComparison.png')
+
+
+#for i in range(len(rates14)):
+#    rates14[i] *= 6.022e23
+#    rates21[i] *= 6.022e23
+#    Te[i] *= 11604
+
+#data14 = np.empty(shape = (len(Te), 2))
+#data21 = np.empty(shape = (len(Te), 2))
+
+#data14[:,0] = Te
+#data21[:,0] = Te
+#data14[:,1] = rates14
+#data21[:,1] = rates21
+
+#with h5.File('./BOLSIGChemistry_6SpeciesRates/deexci.metastable.h5', 'w') as f:
+#    dataset = f.create_dataset("table", data = data14)
+
+#with h5.File('./BOLSIGChemistry_6SpeciesRates/deexci.resonance.h5', 'w') as g:
+#    dataset = g.create_dataset("table", data = data21)
+
+
+

BOLSIGChemistry_NominalRates/plot_rates.py

@@ -0,0 +1,157 @@
+import numpy as np
+import csv
+from matplotlib import pyplot as plt
+from scipy.interpolate import CubicSpline
+import matplotlib.colors as mcolors
+
+tau = (1./13.56e6)
+nAr = 3.22e22
+np0 = 8e16
+Nr = 9
+iSample = 0
+
+reactionExpressionTypelist =  np.array([True,True,True,False,False,False,False,False,False])
+
+reactionsList = []
+	# import h5py as h5
+	# rxnName = ["Ionization", ...] <- dictionary containing reaction name root string
+	# for r in range(Nr):
+	#	if reactionExpressionTypeList[r]:
+	#		fileName = "{0:s}.{1:08d}.h5".format(rxnName[r], iSample)
+	#		f = h5.File(filename, "r")
+	#		D = f["table"]
+
+for i in range(Nr):
+        if reactionExpressionTypelist[i]:
+            Nsample = 1
+            N300 = 200
+
+            root_dir = "./"
+            rate_file = open("{0:s}reaction300K_{1:d}.txt".format(root_dir, i), 'r', encoding='utf-8-sig')
+            temp_file = open('reaction300K_Te.txt', 'r', encoding='utf-8-sig')
+
+            rateCoeff = np.genfromtxt(rate_file, dtype='float')
+            rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample]
+            rate_file.close()
+
+            Te = np.genfromtxt(temp_file, dtype='float')
+            print(Te)
+            Te = np.reshape(Te,[Nsample, N300]).T[:,iSample]
+            print(Te)
+            temp_file.close()
+
+            # Sorting mean energy array and rate coefficient array based on
+            # the mean energy array.
+            Teinds = Te.argsort()
+            rateCoeff = rateCoeff[Teinds]
+            Te = Te[Teinds]
+
+            # Find duplicates
+            TeDuplicateinds = np.where(np.abs(np.diff(Te, axis=0)) > 0.0)
+            TeDuplicateindsForLog = np.where(np.abs(np.diff(Te, axis=0)) == 0.0)
+            rateCoeff = rateCoeff[TeDuplicateinds]
+            Te = Te[TeDuplicateinds]
+
+            # Nondimensionalization of mean energy.
+            # Te *= 1.5
+
+            # Find first non-zero value of the coefficient rate.
+            I = np.nonzero(rateCoeff)
+
+            diffRateCoeff = [j-i for i, j in zip(rateCoeff[:-1], rateCoeff[1:])]
+            diffTe = [j-i for i, j in zip(Te[:-1], Te[1:])]
+
+            Monotonicity = np.asarray([j/i for i, j in zip(diffTe, diffRateCoeff)])
+            Monotonicity = np.insert(Monotonicity, 0, 0.0, axis=0)
+
+            Nan = np.isnan(Monotonicity)
+            Inf = np.isinf(Monotonicity)
+            indexPositive = np.where(Monotonicity>0.0)
+            Positive = np.full(Monotonicity.shape, False, dtype=bool)
+            Positive[indexPositive] = True
+
+            indices = Nan + Inf + Positive
+
+            lastFalse = np.where(indices==False)[-1][-1] + 2
+
+            # Transformation to log scale.
+            TeLog = np.log(Te)
+
+            # Compute the slope of the rate coefficient between its first two non-zero values.
+            # Finite differences are used.
+            dydx = (rateCoeff[lastFalse + 1] - rateCoeff[lastFalse]) \
+                 / (Te[lastFalse + 1] - Te[lastFalse])
+
+            # Arrhenius form: kf = A * exp(-C / Te)
+            C = Te[lastFalse]**2.0*dydx / rateCoeff[lastFalse]
+
+            # Compute pre-exponential coefficient, A, in log scale.
+            ALog = np.log(rateCoeff[lastFalse]) + C / Te[lastFalse]
+
+            # Transform rate coefficient in log scale.
+            rateCoeffLog = np.zeros(rateCoeff.shape)
+            rateCoeffLog[lastFalse:] = np.log(rateCoeff[lastFalse:])
+            # For the troublesome values, we use the Arrhenius form.
+            rateCoeffLog[0:lastFalse] = ALog - C / Te[0:lastFalse]
+            # Nondimensionalization in log scale.
+            #if i < 2:
+            #    rateCoeffLog += - np.log(1.0/tau) + np.log(nAr)
+            #else:
+            #    rateCoeffLog += - np.log(1.0/tau) + np.log(np0)
+
+            # Interpolation in log scale.
+            reactionExpressionsLog = CubicSpline(TeLog, rateCoeffLog)
+            # Gradient in log scale
+            reactionTExpressionsLog = CubicSpline.derivative(reactionExpressionsLog)
+
+            #reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]],
+            #                    rxnBolsig = reactionExpressionTypelist[i],
+            #                    kf_log = reactionExpressionsLog,
+            #                    kf_T_log = reactionTExpressionsLog)
+            #reactionsList.append(reaction)
+
+
+            # rxn   = eval("lambda energy :" + reactionExpressionslist[i])
+            # rxn_T = eval("lambda energy :" + reactionTExpressionslist[i])
+            # setting the axes at the centre
+            fig ,ax = plt.subplots(figsize=(9, 6))
+            ax.spines["top"].set_visible(True)
+            ax.spines["right"].set_visible(True)
+            ax.set_yscale('log')
+            ax.set_xscale('log')
+
+            # plot the function
+            #plt.plot(rateCoeffXFiner, np.exp(reactionExpressions_cubicSplineDerivative_log(rateCoeffXFiner)),
+            #   color='salmon', linestyle='--', label='interBolsig')
+            #plt.plot(rateCoeffXFine, np.exp(reactionExpressions_cubicSpline_log(rateCoeffXFine)),
+            #   color='lightgreen', linestyle='--', label='interBolsig')
+            #plt.plot(Te[:,0], reactionTExpressionsLogFiltered(TeLog[:]) * np.exp(reactionExpressionsLog(TeLog[:])) / Te[:,0],
+            #   color='blue', linestyle='-', label='interBolsig')
+            plt.plot(Te, np.exp(reactionExpressionsLog(TeLog[:])),
+               color='green', linestyle='-', label='Bolsig')
+            #plt.plot(Te, rxn(Te),
+            #   color='salmon', linestyle='--', label='Liu')
+            #plt.plot(Te[:,0], rxn_T(Te[:,0]),
+            #   color='red', linestyle='--', label='interBolsig')
+            plt.xlim((0.05,100))
+            plt.ylim((1e-50,300))
+            plt.legend()
+            plt.savefig("./PlotRates_%s.pdf" %str(i), dpi=300)
+            #plt.xlim((-0.0001,0.0255))
+            plt.show()
+
+            # Export rates to CSV Files
+            header = ['Te', 'kf']
+            row_temp = []
+            data = []
+            for j in range(len(Te)):
+                row_temp.append(Te[j])
+                row_temp.append(np.exp(reactionExpressionsLog(TeLog[j])))
+                data.append(row_temp)
+                row_temp = []
+
+            f = open("./PlotRates_Rxn%s.csv" %str(i), 'w')
+            writer = csv.writer(f)
+            writer.writerow(header)
+            writer.writerows(data)
+            f.close()