diff --git a/.github/workflows/regression_test.yaml b/.github/workflows/regression_test.yaml index 71d55abb8..a21f1d20c 100644 --- a/.github/workflows/regression_test.yaml +++ b/.github/workflows/regression_test.yaml @@ -24,7 +24,7 @@ jobs: - name: Install dependencies run: | python -m pip install --upgrade pip - pip install numpy scipy matplotlib + pip install numpy scipy matplotlib h5py - name: Run Liu case test run: ./testRunBase.sh - name: Run constant background CN test @@ -33,3 +33,5 @@ jobs: run: ./testInterpTrans.sh - name: Run variable background CN test run: ./testRunBackground.sh + - name: Run six species test + run: ./test6Species.sh diff --git a/BOLSIGChemistry_6SpeciesRates/StepwiseExcitations.nominal.h5 b/BOLSIGChemistry_6SpeciesRates/StepwiseExcitations.nominal.h5 new file mode 100644 index 000000000..c02b8063c Binary files /dev/null and b/BOLSIGChemistry_6SpeciesRates/StepwiseExcitations.nominal.h5 differ diff --git a/BOLSIGChemistry_6SpeciesRates/ionization.h5 b/BOLSIGChemistry_6SpeciesRates/ionization.h5 new file mode 100644 index 000000000..1468ca181 Binary files /dev/null and b/BOLSIGChemistry_6SpeciesRates/ionization.h5 differ diff --git a/BOLSIGChemistry_6SpeciesRates/lumped.2p.h5 b/BOLSIGChemistry_6SpeciesRates/lumped.2p.h5 new file mode 100644 index 000000000..b8c33b178 Binary files /dev/null and b/BOLSIGChemistry_6SpeciesRates/lumped.2p.h5 differ diff --git a/BOLSIGChemistry_6SpeciesRates/lumped.metastable.h5 b/BOLSIGChemistry_6SpeciesRates/lumped.metastable.h5 new file mode 100644 index 000000000..a2fdd27e4 Binary files /dev/null and b/BOLSIGChemistry_6SpeciesRates/lumped.metastable.h5 differ diff --git a/BOLSIGChemistry_6SpeciesRates/lumped.resonance.h5 b/BOLSIGChemistry_6SpeciesRates/lumped.resonance.h5 new file mode 100644 index 000000000..da36bca5d Binary files /dev/null and b/BOLSIGChemistry_6SpeciesRates/lumped.resonance.h5 differ diff --git a/BOLSIGChemistry_6SpeciesRates/plot_rates_h5.py b/BOLSIGChemistry_6SpeciesRates/plot_rates_h5.py new file mode 100644 index 000000000..d580735ba --- /dev/null +++ b/BOLSIGChemistry_6SpeciesRates/plot_rates_h5.py @@ -0,0 +1,211 @@ +import numpy as np +import csv +from matplotlib import pyplot as plt +from scipy.interpolate import CubicSpline +import matplotlib.colors as mcolors +import h5py as h5 + +tau = (1./13.56e6) +nAr = 3.22e22 +np0 = 8e16 +Nr = 23 +iSample = 0 + +reactionExpressionTypelist = np.array([False, False, False, False, False, False, False, False, False, + True, True, True, True, False, True, True, True, False, True, True, False, True, + True]) + +rxnNameDict = {0: "2AR(m) => 2AR", + 1: "AR(m) + AR(r) => E + AR+ + AR", + 2: "2AR(4p) => E + AR+ + AR", + 3: "2AR(m) => E + AR+ + AR", + 4: "AR(m) + AR => 2AR", + 5: "AR(r) => AR", + 6: "AR(4p) => AR", + 7: "AR(4p) => AR(m)", + 8: "AR(4p) => AR(r)", + 9: "lumped.metastable", + 10: "lumped.resonance", + 11: "lumped.2p", + 12: "ionization", + 13: "E + AR(m) => E + AR", + 14: "step_ionization", + 15: "E + Ar(m) => E + Ar(r)", + 16: "E + Ar(m) => E + Ar(4p)", + 17: "E + Ar(4p) => 2E + Ar+", + 18: "E + Ar(4p) => E + Ar(r)", + 19: "E + Ar(4p) => E + Ar(m)", + 20: "E + Ar(r) => E + Ar", + 21: "E + Ar(r) => E + Ar(m)", + 22: "E + Ar(r) => E + Ar(4p)"} + +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 = "./" + if (i < 15): + f = h5.File("{0:s}.h5".format(rxnNameDict[i]), 'r') + data = f["table"] + else: + f = h5.File("StepwiseExcitations.nominal.h5", 'r') + data = f[rxnNameDict[i]] + + Te = data[:,0] + rateCoeff = data[:,1] + rateCoeff /= 6.022e23 + + row_temp = [] + data = [] + header = ["Te(K)", "kf(m3/s)"] + for j in range(len(Te)): + row_temp.append(Te[j]) + row_temp.append(rateCoeff[j]) + data.append(row_temp) + row_temp = [] + + g = open("./RateProfile_{0:s}.csv".format(rxnNameDict[i]), 'w') + writer = csv.writer(g) + writer.writerow(header) + writer.writerows(data) + g.close() + + + #rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + #rate_file.close() + + #Te = data[:,0] + + Te /= 11604 + #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 + for k in range(len(Te)): + if (Te[k] < 4.5 and indices[k] == False): + lastFalse = k + 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() diff --git a/BOLSIGChemistry_6SpeciesRates/step_ionization.h5 b/BOLSIGChemistry_6SpeciesRates/step_ionization.h5 new file mode 100644 index 000000000..beabf8706 Binary files /dev/null and b/BOLSIGChemistry_6SpeciesRates/step_ionization.h5 differ diff --git a/Liu2014Properties.py b/Liu2014Properties.py index 5b30da9b2..3ffdaa39d 100644 --- a/Liu2014Properties.py +++ b/Liu2014Properties.py @@ -9,7 +9,7 @@ def __init__(self, *initial_data, **kwargs): for key in kwargs: setattr(self, key, kwargs[key]) -def setLiu2014Properties(gam, inputV0, inputVDC, params, Nr): +def setLiu2014Properties(gam, inputV0, inputVDC, params, Nr, iSample): """Sets non-dimensional properties corresponding to Liu 2014 paper. Inputs: @@ -19,6 +19,8 @@ def setLiu2014Properties(gam, inputV0, inputVDC, params, Nr): Outputs: None params data is overwritten using values from Liu 2014. """ + assert iSample == 0, "This scenario is not set up for sampling" + ################################################################### # User specified parameters (you may change these if you wish to # run a different scenario from Liu 2014) @@ -164,7 +166,7 @@ def setLiu2014Properties(gam, inputV0, inputVDC, params, Nr): rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) reaction = Reaction(rxnAlfa = params.alfa, rxnBeta = params.beta, - kf = rxn, kf_T = rxn_T) + kf = rxn, kf_T = rxn_T, rxnBolsig = False) reactionsList.append(reaction) params.reactionsList = reactionsList diff --git a/chebSolver.py b/chebSolver.py index 2968e6975..6e654858f 100644 --- a/chebSolver.py +++ b/chebSolver.py @@ -1,6 +1,5 @@ import numpy as np import numpy.polynomial.chebyshev as cheb -import matplotlib.pyplot as plt import time from Liu2014Properties import setLiu2014Properties @@ -9,6 +8,13 @@ from psaapPropertiesTestArmInterpTrans import setPsaapPropertiesTestArmInterpTrans from psaapPropertiesCurrentTestCase import setPsaapPropertiesCurrentTestCase from psaapPropertiesCurrentTestCase100mTorr import setPsaapPropertiesCurrentTestCase100mTorr +from psaapPropertiesWithSampling import setPsaapPropertiesWithSampling +from psaapPropertiesTestJP import setPsaapPropertiesTestJP +from psaapPropertiesTestJP_Nominal import setPsaapPropertiesTestJP_Nominal +from psaapPropertiesTestJP_Arrhenius import setPsaapPropertiesTestJP_Arrhenius +from psaapProperties_6Species import setPsaapProperties_6Species +from psaapProperties_6Species_Sampling import setPsaapProperties_6Species_Sampling +from psaapProperties_6Species_Nominal import setPsaapProperties_6Species_Nominal class modelClosures: """Class providing model parameters.""" @@ -89,13 +95,13 @@ def __init__(self, Ns, Nr): self.alfa = np.zeros((Ns,Nr),dtype=np.int64) # reactants # this represents a single rxn: Ar+e -> Ar+ + e + e - self.beta[0,0] = 2 - self.beta[1,0] = 1 - self.beta[2,0] = 0 + #self.beta[0,0] = 2 + #self.beta[1,0] = 1 + #self.beta[2,0] = 0 - self.alfa[0,0] = 1 - self.alfa[1,0] = 0 - self.alfa[2,0] = 1 + #self.alfa[0,0] = 1 + #self.alfa[1,0] = 0 + #self.alfa[2,0] = 1 # other non-dimensional parameters self.qStar = 100.0 @@ -187,10 +193,20 @@ def diffusivity(self, i, energy, mu, nb, EinsteinForm): DEf[:,0] = self.diffusivityList[i].D_expression((2./3)*energy[:,i]) / nb - elif EinsteinForm: + elif EinsteinForm and self.Z[i] == -1: V0 = self.qStar * 1.0 # V0 = qStar * 1eV DEf = 2.0 / 3.0 * np.multiply(energy[:,[i]], mu[:,[i]]) / V0 + # Einstein for electrons only (for testing purposes) + # if EinsteinForm: + # #V0 = self.qStar * 1.0 # V0 = qStar * 1eV + # #DEf = 2.0 / 3.0 * np.multiply(energy[:,[i]], mu[:,[i]]) / V0 + # if self.Z[i] == -1: ## Added this to use Einstein Relation only with electrons and constant values for heavies + # V0 = self.qStar * 1.0 + # DEf = 2.0 / 3.0 * np.multiply(energy[:,[i]], mu[:,[i]]) / V0 + # else: + # DEf[:,0] = self.D[i] / nb + else: DEf[:,0] = self.D[i] / nb @@ -219,6 +235,12 @@ def diffusivity_U(self, i, j, energy, energy_U, mu, D, nb, EinsteinForm): if (j == self.Ns - 1): D_U[:,:] -= np.diag(D[:,i] / nb) + # Einstein for electrons only (for testing purposes) + # if EinsteinForm: + # if self.Z[i] == -1: + # V0 = self.qStar * 1.0 # V0 = qStar * 1eV + # D_U = 2.0 / 3.0 * np.multiply(mu[:,[i]], energy_U[i,j,:,:]) / V0 + return D_U @@ -292,7 +314,10 @@ def rxnRateCoefficient(self, energy, i): indFix = (energy[:,0]<=0.0) energy[indFix,0] = 1.0 - kf = self.reactionsList[i].kf(energy) + if self.reactionsList[i].rxnBolsig: + kf = np.exp(self.reactionsList[i].kf_log(np.log(energy))) + else: + kf = self.reactionsList[i].kf(energy) kf[indFix,0] = 0 return kf #a * (energy**b) * np.exp(-Ea/energy) @@ -304,7 +329,11 @@ def rxnRateCoefficientJac(self, energy, i): indFix = (energy[:,0]<=0.0) energy[indFix,0] = 1.0 - kf_T = self.reactionsList[i].kf_T(energy) + if self.reactionsList[i].rxnBolsig: + kf_T = self.reactionsList[i].kf_T_log(np.log(energy)) \ + * np.exp(self.reactionsList[i].kf_log(np.log(energy))) / energy + else: + kf_T = self.reactionsList[i].kf_T(energy) kf_T[indFix,0] = 0 return kf_T #a * (energy**(b-1)) * np.exp(-Ea/energy) * (b + Ea/energy) @@ -359,7 +388,7 @@ class timeDomainCollocationSolver: def __init__(self, Ns, NT, Np, elasticCollisionActivationFactor, backgroundSpecieActivationFactor, EinsteinForm, gam=0.01, V0 = 100.0, VDC = 0.0, - scenario=0, scheme="BE"): + scenario=0, scheme="BE", iSample = 0): """Initializes storage and operaters required for solve.""" # parameters of the time marching scheme @@ -399,6 +428,20 @@ def __init__(self, Ns, NT, Np, elasticCollisionActivationFactor, Nr = 8 elif(scenario==4): Nr = 8 + elif(scenario==5): + Nr = 7 + elif(scenario==6): + Nr = 9 + elif(scenario==7): + Nr = 9 + elif(scenario==8): + Nr = 9 + elif(scenario==9): + Nr = 23 + elif(scenario==10): + Nr = 23 + elif(scenario==12): + Nr = 23 elif(scenario==21): Nr = 8 else: @@ -412,18 +455,31 @@ def __init__(self, Ns, NT, Np, elasticCollisionActivationFactor, self.params = modelClosures(self.Ns, Nr) if(scenario==0): - setLiu2014Properties(gam, V0, VDC, self.params, Nr) + setLiu2014Properties(gam, V0, VDC, self.params, Nr, iSample) elif(scenario==1): - setPsaapProperties(gam, V0, VDC, self.params, Nr) + setPsaapProperties(gam, V0, VDC, self.params, Nr, iSample) elif(scenario==2): - setPsaapPropertiesTestArm(gam, V0, VDC, self.params, Nr) + setPsaapPropertiesTestArm(gam, V0, VDC, self.params, Nr, iSample) elif(scenario==3): - setPsaapPropertiesCurrentTestCase(gam, V0, VDC, self.params, Nr) + setPsaapPropertiesCurrentTestCase(gam, V0, VDC, self.params, Nr, iSample) elif(scenario==4): - setPsaapPropertiesCurrentTestCase100mTorr(gam, V0, VDC, self.params, Nr) + setPsaapPropertiesCurrentTestCase100mTorr(gam, V0, VDC, self.params, Nr, iSample) + elif(scenario==5): + setPsaapPropertiesWithSampling(gam, V0, VDC, self.params, Nr, iSample) + elif(scenario==6): + setPsaapPropertiesTestJP(gam, V0, VDC, self.params, Nr, iSample) + elif(scenario==7): + setPsaapPropertiesTestJP_Nominal(gam, V0, VDC, self.params, Nr, iSample) + elif(scenario==8): + setPsaapPropertiesTestJP_Arrhenius(gam, V0, VDC, self.params, Nr, iSample) + elif(scenario==9): + setPsaapProperties_6Species(gam, V0, VDC, self.params, Nr, iSample) + elif(scenario==10): + setPsaapProperties_6Species_Sampling(gam, V0, VDC, self.params, Nr, iSample) + elif(scenario==12): + setPsaapProperties_6Species_Nominal(gam, V0, VDC, self.params, Nr, iSample) elif(scenario==21): - setPsaapPropertiesTestArmInterpTrans(gam, V0, VDC, self.params, Nr) - + setPsaapPropertiesTestArmInterpTrans(gam, V0, VDC, self.params, Nr, iSample) # Points used to define state and collocation # (Gauss-Lobatto-Chebyshev points) @@ -822,13 +878,6 @@ def residualCN(self, Uin, time, dt, weak_bc=False): res[ self.Ns*self.Np-1 ] = dens[ -1,self.Ns-1] \ - ((self.params.p0 - nT[ -1]) / self.params.Tg0 - ntot[-1]) / self.params.nAronp0 - # enforce Dirichlet condition on heavy species temperature - ntot = np.zeros(self.Np) - - # add all heavies but background - for i in range(1, self.Ns-1): - ntot += dens[:,i] - # electron temperature res[self.Ns*self.Np ] = (nT[ 0] - self.params.EeBC*dens[0,iele]) res[(self.Ns+1)*self.Np-1] = (nT[-1] - self.params.EeBC*dens[-1,iele]) @@ -869,9 +918,22 @@ def residualLCN(self, Uin, time, dt, weak_bc=False): # change in the associated variables is zero. However, this # *only* works if the IC satisfies the BC. Any errors # introduced by the IC will never be eliminated. - if (self.Ns>2): - res[2*self.Np ] = 0.0 #dens[ 0,2] - 0.0 - res[3*self.Np-1] = 0.0 #dens[-1,2] - 0.0 + #if (self.Ns>2): + # res[2*self.Np ] = 0.0 #dens[ 0,2] - 0.0 + # res[3*self.Np-1] = 0.0 #dens[-1,2] - 0.0 + + if (self.Ns > 2): + for i in range(2, self.Ns-1): + res[i*self.Np ] = 0.0 + res[(i+1)*self.Np-1] = 0.0 + + #if (self.Ns == 6): + #res[2*self.Np ] = 0.0 + #res[3*self.Np-1] = 0.0 + # res[3*self.Np ] = 0.0 + # res[4*self.Np-1] = 0.0 + # res[4*self.Np ] = 0.0 + # res[5*self.Np-1] = 0.0 # electron temperature res[self.Ns*self.Np ] = 0.0 #(nT[ 0] - 0.75*dens[0,iele]) @@ -1364,12 +1426,20 @@ def jacobianLCN(self, Uin, time, dt, weak_bc=False): print("Error: Only weak electron flux BCs supported for linearized CN.") exit(-1) - if (self.Ns>2): - self.jac[2*self.Np,:] = np.zeros((1,self.Nv*self.Np)) - self.jac[2*self.Np,2*self.Np] = 1.0 + #if (self.Ns>2): + # self.jac[2*self.Np,:] = np.zeros((1,self.Nv*self.Np)) + # self.jac[2*self.Np,2*self.Np] = 1.0 + + # self.jac[3*self.Np-1,:] = np.zeros((1,self.Nv*self.Np)) + # self.jac[3*self.Np-1,3*self.Np-1] = 1.0 + + if (self.Ns > 2): + for i in range(2,self.Ns-1): + self.jac[i*self.Np,:] = np.zeros((1,self.Nv*self.Np)) + self.jac[i*self.Np,i*self.Np] = 1.0 + self.jac[(i+1)*self.Np-1,:] = np.zeros((1,self.Nv*self.Np)) + self.jac[(i+1)*self.Np-1,(i+1)*self.Np-1] = 1.0 - self.jac[3*self.Np-1,:] = np.zeros((1,self.Nv*self.Np)) - self.jac[3*self.Np-1,3*self.Np-1] = 1.0 self.jac[self.Ns*self.Np,:] = np.zeros((1,self.Nv*self.Np)) self.jac[self.Ns*self.Np,self.Ns*self.Np] = 1.0 @@ -1412,9 +1482,10 @@ def jacobian0(self, time, dt, weak_bc=False): self.jac0[0 ,:] = np.zeros((1,self.Nv*self.Np)) self.jac0[self.Np-1,:] = np.zeros((1,self.Nv*self.Np)) - if (self.Ns>2): - self.jac0[2*self.Np,:] = np.zeros((1,self.Nv*self.Np)) - self.jac0[3*self.Np-1,:] = np.zeros((1,self.Nv*self.Np)) + for i in range(2,self.Ns-1): + self.jac0[i*self.Np,:] = np.zeros((1,self.Nv*self.Np)) + self.jac0[(i+1)*self.Np-1,:] = np.zeros((1,self.Nv*self.Np)) + # Dirichlet on heavy species temperature if (self.backgroundSpecieActivationFactor > 0): @@ -1792,6 +1863,8 @@ def plot(self, col, create=True): action='store_true', help="Activate the background specie density equation.") parser.add_argument('--EinsteinForm', default=False, action='store_true', help="Activate Einstein's form for diffusion coefficient.") + parser.add_argument('--iSample', metavar='iSample', default=0, + type=int, help='Sample index, if BOLSIG chemistry is used.') args = parser.parse_args() # Dump inputs to the screen for posterity @@ -1836,11 +1909,29 @@ def plot(self, col, create=True): elif(args.scenario==4): print("# Running scenario = 4 (4 species, 8 rxn, Liu 2017)") Ns = 4 + elif(args.scenario==5): + print("# Running scenario = 5 (4 species, 7 rxn, Bolsing and Lay, Moss et al, 2003)") + Ns = 4 + elif(args.scenario==6): + print("# Running scenario = 6 (4 species, 9 rxn, Juan's mechanism)") + Ns = 4 + elif(args.scenario==7): + print("# Running scenario = 7 (4 species, 9 rxn, Nominal reaction rates)") + Ns = 4 + elif(args.scenario==8): + Ns = 4 + elif(args.scenario==9): + print("# Running scenario = 9 (6 species, 23 rxn)") + Ns = 6 + elif(args.scenario==10): + Ns = 6 + elif(args.scenario==12): + Ns = 6 elif(args.scenario==21): print("# Running scenario = 21 (4 species, 8 rxn, Liu 2017, interpolated transport)") Ns = 4 else: - print("ERROR: Scenario = {0:d} not recognized. Use --scenario i with i=0, 1, 2, or 3. Exiting.".format(args.scenario)) + print("ERROR: Scenario = {0:d} not recognized. Exiting.".format(args.scenario)) exit(-1) elasticCollisionActivationFactor = 1.0 @@ -1880,7 +1971,8 @@ def plot(self, col, create=True): tds = timeDomainCollocationSolver(Ns, 1, args.Np, elasticCollisionActivationFactor, backgroundSpecieActivationFactor, EinsteinForm, gam=0.01, V0 = args.V0, VDC = args.VDC, - scenario=args.scenario, scheme=args.tscheme) + scenario=args.scenario, scheme=args.tscheme, + iSample = args.iSample) # Default IC (overwritten below if we are restarting) #tds.U1[0:tds.Ns*tds.Np] = 1e-4 diff --git a/generate_sample.py b/generate_sample.py new file mode 100644 index 000000000..4c33dc885 --- /dev/null +++ b/generate_sample.py @@ -0,0 +1,57 @@ +import numpy as np +import h5py +import yaml +import argparse +import re + +# reaction units based on the number of reactants. +unitDict = {1: '1/s', 2: 'm3/mol/s', 3: 'm6/mol2/s'} +parser = argparse.ArgumentParser(description = "", formatter_class = argparse.RawTextHelpFormatter) +parser.add_argument('input_file', metavar = 'string', type = str, help = 'filename for an input file.\n') + + +def getNumberOfReactants(equation): + reactants, direction, products = re.split("( => | <=> )", equation) + reactants = re.split(" \+ ", reactants) + return len(reactants) + + +if __name__ == '__main__': + args = parser.parse_args() + with open(args.input_file) as c: + dict_ = yaml.safe_load(c) + + # parse arrhenius-type reactions from the input file. + rxnList = dict_['reactions'] + for rxn in rxnList: + if (rxn['rate_type'] == 'arrhenius'): + print("Equation: %s" % rxn['equation']) + print("Nominal coefficients: %s" % rxn['arrhenius']['coefficients']) + print("%d-body reaction" % getNumberOfReactants(rxn['equation'])) + + # create the directory for samples. + import os + outputDir = dict_['directory'] + if (outputDir != '.'): + os.makedirs(outputDir, exist_ok = True) + + prefix = dict_['prefix'] + error = 0.01* dict_['relative_error'] + nSample = dict_['number_of_samples'] + for k in range(nSample): + filename = '%s/%s.%08d.h5' % (outputDir, prefix, k) + with h5py.File(filename, 'w') as f: + for rxn in rxnList: + if (rxn['rate_type'] =='arrhenius'): + rateUnit = unitDict[getNumberOfReactants(rxn['equation'])] + # Take the nominal value from the input file + coeffs0 = np.array(rxn['arrhenius']['coefficients'], dtype = np.double) + # Add a relative error + coeffs0[0] *= (1.0 + error) ** np.random.normal() + # Create a dataset for the reaction + dset = f.create_dataset(rxn['equation'], (3,), data = coeffs0) + dset.attrs['rate_unit'] = rateUnit + dset.attrs['temperature_unit'] = 'K' + if (k % 10 == 0): + print("%d-th sample generated." % k) + diff --git a/plot_results.py b/plot_results.py new file mode 100644 index 000000000..3e6d197b9 --- /dev/null +++ b/plot_results.py @@ -0,0 +1,268 @@ +import numpy as np +import matplotlib.pyplot as plt +import csv +import sys +#import chebSolver as cs +import numpy.polynomial.chebyshev as cheb +import matplotlib +import matplotlib.font_manager +from matplotlib import rc +rc('font',**{'family':'serif','serif':['Times']}) +rc('text', usetex=True) + +sys.path.insert(0, '../') +import chebSolver as cs + +# these values are required to "redimensionalize" the results +# they must be consistent with the scenario input file +ne0 = 8e16 # "nominal" electron density [1/m^3] +L = 2.00*0.005 # half-gap-width [m] (gap width is 1in) +tau = (1./13.56e6) # period of driving voltage [s] +V0 = 100.0 + +# Number of Chebyshev modes +Np=150 + +# load solution file +D = np.load('newton_4spec_CN_Np150_fullsoln.npy') +D = np.transpose(D) + +# T = number of time steps +T=128 + +# indices to plot (0, 0.25, 0.5, 0.75)*period +t0 = 0 +t1 = 32 +t2 = 64 +t3 = 96 + +# time +tp = np.linspace(0,1,T+1) +tr = tp*tau + +# spatial grid +xp = -np.cos(np.pi*np.linspace(0,Np-1,Np)/(Np-1)) +xr = (xp+1)*L + +# pull solution out of D +ne = D[0:Np,:] # electron density +ni = D[Np:2*Np,:] # ion density +nm = D[2*Np:3*Np,:] # metastable density +nb = D[3*Np:4*Np,:] # "background" (argon neutral) density +nee = D[4*Np:,:] # electron energy (ne * ee) + + +# (3/2)*electron temp +Te = D[4*Np:,:]/ne + + +# evaluate electric potential (this requires a poisson solve, which is +# why the timeDomainCollocationSolver is loaded) +#phi = np.zeros((Np,T+1)) +#tds = cs.timeDomainCollocationSolver(4,1,Np,scenario=6) +#print("Solving poisson...") +#for i in range(0,T+1): +# tds.solve_poisson(ne[:,i],ni[:,i],tp[i]) +# phi[:,i] = tds.phi[:] + +# make some plots + +# ne +fig,ax = plt.subplots() +h = ax.contourf(xr, tr, np.transpose(ne)*ne0, levels=np.linspace(0,ne0,129)) +cbar = plt.colorbar(h, ticks=np.linspace(0,ne0,9)) +cbar.ax.set_ylabel(r"$n_e$ [m$^{-3}$]", labelpad=6, fontsize=18) +plt.setp(cbar.ax.get_yticklabels(), rotation='horizontal', fontsize=12) +ax.contour(xr, tr, np.transpose(ne)*ne0, levels=np.linspace(0,ne0,129)) +ax.set_xlabel(r"$x$ [m]", fontsize=18) +plt.setp(ax.get_xticklabels(), fontsize=12) +ax.set_ylabel(r"$t$ [s]", fontsize=18) +plt.setp(ax.get_yticklabels(), fontsize=12) +plt.savefig('ne_4spec_contour.pdf') + +# ni +fig,ax = plt.subplots() +h = ax.contourf(xr, tr, np.transpose(ni)*ne0, levels=np.linspace(0,ne0,129)) +cbar = plt.colorbar(h, ticks=np.linspace(0,ne0,9)) +cbar.ax.set_ylabel(r"$n_i$ [m$^{-3}$]", labelpad=6, fontsize=18) +plt.setp(cbar.ax.get_yticklabels(), rotation='horizontal', fontsize=12) +ax.contour(xr, tr, np.transpose(ni)*ne0, levels=np.linspace(0,ne0,129)) +ax.set_xlabel(r"$x$ [m]", fontsize=18) +plt.setp(ax.get_xticklabels(), fontsize=12) +ax.set_ylabel(r"$t$ [s]", fontsize=18) +plt.setp(ax.get_yticklabels(), fontsize=12) +plt.savefig('ni_4spec_contour.pdf') + +# nm +fig,ax = plt.subplots() +h = ax.contourf(xr, tr, np.transpose(nm)*ne0, levels=np.linspace(0,5*ne0,129)) +cbar = plt.colorbar(h, ticks=np.linspace(0,ne0,9)) +cbar.ax.set_ylabel(r"$n^{\ast}$ [m$^{-3}$]", labelpad=6, fontsize=18) +plt.setp(cbar.ax.get_yticklabels(), rotation='horizontal', fontsize=12) +ax.contour(xr, tr, np.transpose(ni)*ne0, levels=np.linspace(0,ne0,129)) +ax.set_xlabel(r"$x$ [m]", fontsize=18) +plt.setp(ax.get_xticklabels(), fontsize=12) +ax.set_ylabel(r"$t$ [s]", fontsize=18) +plt.setp(ax.get_yticklabels(), fontsize=12) +plt.savefig('nm_4spec_contour.pdf') + +# Te +fig,ax = plt.subplots() +h = ax.contourf(xr, tr, np.transpose(Te)*(2./3.), levels=np.linspace(0,16,129)) +cbar = plt.colorbar(h, ticks=np.linspace(0,16,5)) +cbar.ax.set_ylabel(r"$T_e$ [eV]", labelpad=6, fontsize=18) +plt.setp(cbar.ax.get_yticklabels(), rotation='horizontal', fontsize=12) +ax.contour(xr, tr, np.transpose(Te)*(2./3.), levels=np.linspace(0,16,129)) +ax.set_xlabel(r"$x$ [m]", fontsize=18) +plt.setp(ax.get_xticklabels(), fontsize=12) +ax.set_ylabel(r"$t$ [s]", fontsize=18) +plt.setp(ax.get_yticklabels(), fontsize=12) +plt.savefig('Te_4spec_contour.pdf') + +# nb +fig,ax = plt.subplots() +h = ax.contourf(xr, tr, np.transpose(nb)*ne0, levels=np.linspace(0,5*ne0,129)) +cbar = plt.colorbar(h, ticks=np.linspace(0,ne0,9)) +cbar.ax.set_ylabel(r"$n_b$ [m$^{-3}$]", labelpad=6,fontsize=18) +plt.setp(cbar.ax.get_yticklabels(), rotation='horizontal', fontsize=12) +ax.contour(xr, tr, np.transpose(nb)*ne0, levels=np.linspace(0,ne0,129)) +ax.set_xlabel(r"$x$ [m]", fontsize=18) +plt.setp(ax.get_xticklabels(), fontsize=12) +ax.set_ylabel(r"$t$ [s]", fontsize=18) +plt.setp(ax.get_yticklabels(), fontsize=12) +plt.savefig('nb_4spec_contour.pdf') + +# phi +#fig,ax = plt.subplots() +#h = ax.contourf(xr, tr, np.transpose(phi)*V0, levels=np.linspace(-120,120,129)) +#cbar = plt.colorbar(h, ticks=np.linspace(-120,120,13)) +#cbar.ax.set_ylabel(r"$\phi$ [V]", labelpad=6, fontsize=18) +#plt.setp(cbar.ax.get_yticklabels(), rotation='horizontal', fontsize=12) +#ax.contour(xr, tr, np.transpose(phi)*V0, levels=np.linspace(-120,120,129)) +#ax.set_xlabel(r"$x$ [m]", fontsize=18) +#plt.setp(ax.get_xticklabels(), fontsize=12) +#ax.set_ylabel(r"$t$ [s]", fontsize=18) +#plt.setp(ax.get_yticklabels(), fontsize=12) +#plt.savefig('phi_4spec_contour.pdf') + +print("Plotting lines...") + +# ne +fig,ax = plt.subplots() +ax.plot(xr, ne0*ne[:,t0], 'b-', lw=2, label=r"t={0:.2f}T".format(tp[t0])) +ax.plot(xr, ne0*ne[:,t1], 'r-', lw=2, label=r"t={0:.2f}T".format(tp[t1])) +ax.plot(xr, ne0*ne[:,t2], 'g-', lw=2, label=r"t={0:.2f}T".format(tp[t2])) +ax.plot(xr, ne0*ne[:,t3], 'c-', lw=2, label=r"t={0:.2f}T".format(tp[t3])) +ax.plot(xr, ne0*ni[:,t0], 'b--', lw=2) +ax.plot(xr, ne0*ni[:,t1], 'r--', lw=2) +ax.plot(xr, ne0*ni[:,t2], 'g--', lw=2) +ax.plot(xr, ne0*ni[:,t3], 'c--', lw=2) +ax.plot(xr, ne0*nm[:,t0], 'b:', lw=2) +ax.plot(xr, ne0*nm[:,t1], 'r:', lw=2) +ax.plot(xr, ne0*nm[:,t2], 'g:', lw=2) +ax.plot(xr, ne0*nm[:,t3], 'c:', lw=2) +ax.legend() +ax.set_xlim((0,0.0254)) +ax.set_xlabel(r"$x$ [m]", fontsize=18) +plt.setp(ax.get_xticklabels(), fontsize=12) +ax.set_ylabel(r"Density [m$^{-3}$]", fontsize=18) +#ax.set_ylim((0,8e16)) +plt.setp(ax.get_yticklabels(), fontsize=12) +plt.savefig('density_4spec_line.pdf') + +# Te/phi +fig,ax = plt.subplots() +ax.plot(xr, (2./3.)*Te[:,t0], 'b-', lw=2, label=r"t={0:.2f}T".format(tp[t0])) +ax.plot(xr, (2./3.)*Te[:,t1], 'r-', lw=2, label=r"t={0:.2f}T".format(tp[t1])) +ax.plot(xr, (2./3.)*Te[:,t2], 'g-', lw=2, label=r"t={0:.2f}T".format(tp[t2])) +ax.plot(xr, (2./3.)*Te[:,t3], 'c-', lw=2, label=r"t={0:.2f}T".format(tp[t3])) +ax.set_xlim((0,0.0254)) +ax.set_xlabel(r"$x$ [m]", fontsize=18) +plt.setp(ax.get_xticklabels(), fontsize=12) +ax.set_ylim((0,30)) +ax.set_ylabel(r"$T_e$ [eV]", fontsize=18) +plt.setp(ax.get_yticklabels(), fontsize=12) + +#ax2 = ax.twinx() +#ax2.plot(xr, V0*phi[:,t0], 'b--', lw=2) +#ax2.plot(xr, V0*phi[:,t1], 'r--', lw=2) +#ax2.plot(xr, V0*phi[:,t2], 'g--', lw=2) +#ax2.plot(xr, V0*phi[:,t3], 'c--', lw=2) +#ax2.set_ylim((-120,120)) +#ax2.set_ylabel(r"$\phi$ [V]", labelpad=-5, fontsize=18) +#plt.setp(ax2.get_yticklabels(), fontsize=12) +plt.savefig('Te_4spec_line.pdf') + +print("Plotting means...") + +# ne +fig,ax = plt.subplots() +ax.semilogy(xr, ne0*np.mean(ne,axis=1), 'b-', lw=2, label=r"$n_e$") +ax.semilogy(xr, ne0*np.mean(ni,axis=1), 'r--', lw=2, label=r"$n_i$") +#ax.plot(xr, ne0*np.mean(nm,axis=1), 'g:', lw=2, label=r"$n^{\ast}$") +ax.legend(fontsize=12) +#ax.set_xlim((0,0.0254)) +ax.set_xlim((0,0.02)) +ax.set_xlabel(r"$x$ [m]", fontsize=18) +plt.setp(ax.get_xticklabels(), fontsize=12) +#ax.set_ylim((0,8e16)) +ax.set_ylabel(r"Density [m$^{-3}$]", fontsize=18) +plt.setp(ax.get_yticklabels(), fontsize=12) +plt.savefig('density_4spec_mean.pdf') + +# nm +fig,ax = plt.subplots() +ax.plot(xr, ne0*np.mean(nm,axis=1), 'g-', lw=2, label=r"$n^{\ast}$") +ax.legend(fontsize=12) +ax.set_xlim((0,0.0254)) +ax.set_xlabel(r"$x$ [m]", fontsize=18) +plt.setp(ax.get_xticklabels(), fontsize=12) +#ax.set_ylim((0,4*8e16)) +ax.set_ylabel(r"Density [m$^{-3}$]", fontsize=18) +plt.setp(ax.get_yticklabels(), fontsize=12) +plt.savefig('metastable_4spec_mean.pdf') + +# Te/phi +fig,ax = plt.subplots() +ax.plot(xr, (2./3.)*np.mean(Te,axis=1), 'b-', lw=2, label=r"$T_e$") +ax.legend(fontsize=12,loc=2) +ax.set_xlim((0,0.0254)) +ax.set_xlabel(r"$x$ [m]", fontsize=18) +plt.setp(ax.get_xticklabels(), fontsize=12) +ax.set_ylim((0,5.5)) +ax.set_ylabel(r"$T_e$ [eV]", fontsize=18) +plt.setp(ax.get_yticklabels(), fontsize=12) +# +#ax2 = ax.twinx() +#ax2.plot(xr, V0*np.mean(phi,axis=1), 'r--', lw=2, label=r"$\phi$") +#ax2.legend(fontsize=12,loc=1) +##ax2.set_ylim((0,70)) +#ax2.set_ylabel(r"$\phi$ [V]", fontsize=18) +#plt.setp(ax2.get_yticklabels(), fontsize=12) +plt.savefig('Te_4spec_mean.pdf') + +# Write CSV File +header = ['x', 'ne', 'ni', 'nm', 'nb', 'Te'] +data = [] +row = [] +for i in range(len(xr)): + x_temp = xr[i] + ne_temp = ne0*np.mean(ne,axis=1)[i] + ni_temp = ne0*np.mean(ni,axis=1)[i] + nm_temp = ne0*np.mean(nm,axis=1)[i] + nb_temp = ne0*np.mean(nb,axis=1)[i] + Te_temp = (2./3.)*11604.*np.mean(Te,axis=1)[i] + row.append(x_temp) + row.append(ne_temp) + row.append(ni_temp) + row.append(nm_temp) + row.append(nb_temp) + row.append(Te_temp) + data.append(row) + row = [] + +f = open('GD_mean_results.csv', 'w') +writer = csv.writer(f) +writer.writerow(header) +writer.writerows(data) +f.close() diff --git a/psaapProperties.py b/psaapProperties.py index bac59d238..549ce9335 100644 --- a/psaapProperties.py +++ b/psaapProperties.py @@ -1,4 +1,5 @@ import numpy as np +from scipy.interpolate import CubicSpline class Reaction(object): @@ -9,7 +10,7 @@ def __init__(self, *initial_data, **kwargs): for key in kwargs: setattr(self, key, kwargs[key]) -def setPsaapProperties(gam, inputV0, inputVDC, params, Nr): +def setPsaapProperties(gam, inputV0, inputVDC, params, Nr, iSample): """Sets non-dimensional properties corresponding to Liu 2014 paper. Inputs: @@ -158,14 +159,73 @@ def setPsaapProperties(gam, inputV0, inputVDC, params, Nr): reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)"] + reactionExpressionTypelist = [False] + reactionsList = [] for i in range(Nr): - rxn = eval("lambda energy :" + reactionExpressionslist[i]) - rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) - - reaction = Reaction(rxnAlfa = params.alfa, rxnBeta = params.beta, - kf = rxn, kf_T = rxn_T) - reactionsList.append(reaction) + if reactionExpressionTypelist[i]: + Nsample = 72 + N300 = 200 + + rateCoeff = np.fromfile('./BOLSIGChemistry/reaction300K_%s.dat' %str(i)) + rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + + Te = np.fromfile('./BOLSIGChemistry/reaction300K.Te.dat') + Te = np.reshape(Te,[Nsample, N300]).T[:,iSample] + + # Nondimensionalization of mean energy and transformation to log scale. + # Te *= 1.5 + TeLog = np.log(Te) + + # Find first non-zero value of the coefficient rate. + I = np.nonzero(rateCoeff) + + # Compute the slope of the rate coefficient between its first two non-zero values. + # Finite differences are used. + dydx = (rateCoeff[I[0][0] + 1] - rateCoeff[I[0][0]]) \ + / (Te[I[0][0] + 1] - Te[I[0][0]]) + + # Arrhenius form: kf = A * exp(-C / Te) + # C = (dkf/dTe) / kf * Te**2.0 + # A = kf / exp(-C / Te) + C = Te[I[0][0]]**2.0*dydx / rateCoeff[I[0][0]] + # A = rateCoeff[I[0][0]] / np.exp(-C/Te[I[0][0]]) + + # Compute pre-exponential coefficient, A, in log scale. + ALog = np.log(rateCoeff[I[0][0]]) + C / Te[I[0][0]] + + # Transform rate coefficient in log scale. + rateCoeffLog = np.zeros(rateCoeff.shape) + rateCoeffLog[I[0][0]:] = np.log(rateCoeff[I[0][0]:]) + # For the troublesome values, we use the Arrhenius form. + rateCoeffLog[0:I[0][0]] = ALog - C / Te[0:I[0][0]] + # 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) + # Nondimensionalization of the original rate, used for the plot and comparison. + rateCoeff *= tau * nAr + + # Interpolation in log scale. + reactionExpressionsLog = CubicSpline(TeLog, rateCoeffLog) + # Gradient in log scale + reactionTExpressionsLog = CubicSpline.derivative(reactionExpressionsLog) + + reaction = Reaction(rxnAlfa = params.alfa, rxnBeta = params.beta, + rxnBolsig = reactionExpressionTypelist[i], + kf_log = reactionExpressionsLog, + kf_T_log = reactionTExpressionsLog) + reactionsList.append(reaction) + + else: + rxn = eval("lambda energy :" + reactionExpressionslist[i]) + rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) + + reaction = Reaction(rxnAlfa = params.alfa, rxnBeta = params.beta, + rxnBolsig = reactionExpressionTypelist[i], + kf = rxn, kf_T = rxn_T) + reactionsList.append(reaction) params.reactionsList = reactionsList diff --git a/psaapPropertiesCurrentTestCase.py b/psaapPropertiesCurrentTestCase.py index ddf9701eb..b372b994b 100644 --- a/psaapPropertiesCurrentTestCase.py +++ b/psaapPropertiesCurrentTestCase.py @@ -9,7 +9,7 @@ def __init__(self, *initial_data, **kwargs): for key in kwargs: setattr(self, key, kwargs[key]) -def setPsaapPropertiesCurrentTestCase(gam, inputV0, inputVDC, params, Nr): +def setPsaapPropertiesCurrentTestCase(gam, inputV0, inputVDC, params, Nr, iSample): """Sets non-dimensional properties corresponding to Liu 2014 paper. Inputs: diff --git a/psaapPropertiesCurrentTestCase100mTorr.py b/psaapPropertiesCurrentTestCase100mTorr.py index 24cfa7a58..60bf77c14 100644 --- a/psaapPropertiesCurrentTestCase100mTorr.py +++ b/psaapPropertiesCurrentTestCase100mTorr.py @@ -9,7 +9,7 @@ def __init__(self, *initial_data, **kwargs): for key in kwargs: setattr(self, key, kwargs[key]) -def setPsaapPropertiesCurrentTestCase100mTorr(gam, inputV0, inputVDC, params, Nr): +def setPsaapPropertiesCurrentTestCase100mTorr(gam, inputV0, inputVDC, params, Nr, iSample): """Sets non-dimensional properties corresponding to Liu 2014 paper. Inputs: diff --git a/psaapPropertiesTestArm.py b/psaapPropertiesTestArm.py index 1ed7fdb8d..b10553372 100644 --- a/psaapPropertiesTestArm.py +++ b/psaapPropertiesTestArm.py @@ -25,7 +25,7 @@ def __init__(self, *initial_data, **kwargs): for key in kwargs: setattr(self, key, kwargs[key]) -def setPsaapPropertiesTestArm(gam, inputV0, inputVDC, params, Nr): +def setPsaapPropertiesTestArm(gam, inputV0, inputVDC, params, Nr, iSample): """Sets non-dimensional properties corresponding to Liu 2014 paper. Inputs: @@ -35,6 +35,8 @@ def setPsaapPropertiesTestArm(gam, inputV0, inputVDC, params, Nr): Outputs: None params data is overwritten using values from Liu 2014. """ + assert iSample == 0, "This scenario is not set up for sampling" + ################################################################### # User specified parameters (you may change these if you wish to # run a different scenario from Liu 2014) @@ -239,7 +241,7 @@ def setPsaapPropertiesTestArm(gam, inputV0, inputVDC, params, Nr): rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]], - kf = rxn, kf_T = rxn_T) + kf = rxn, kf_T = rxn_T, rxnBolsig = False ) reactionsList.append(reaction) params.reactionsList = reactionsList diff --git a/psaapPropertiesTestArmInterpTrans.py b/psaapPropertiesTestArmInterpTrans.py index db88f0e2c..c1c6f0626 100644 --- a/psaapPropertiesTestArmInterpTrans.py +++ b/psaapPropertiesTestArmInterpTrans.py @@ -25,7 +25,7 @@ def __init__(self, *initial_data, **kwargs): for key in kwargs: setattr(self, key, kwargs[key]) -def setPsaapPropertiesTestArmInterpTrans(gam, inputV0, inputVDC, params, Nr): +def setPsaapPropertiesTestArmInterpTrans(gam, inputV0, inputVDC, params, Nr, iSample): """Sets non-dimensional properties corresponding to Liu 2014 paper. Inputs: @@ -35,6 +35,7 @@ def setPsaapPropertiesTestArmInterpTrans(gam, inputV0, inputVDC, params, Nr): Outputs: None params data is overwritten using values from Liu 2014. """ + assert iSample == 0, "This scenario is not set up for sampling" ################################################################### # User specified parameters (you may change these if you wish to # run a different scenario from Liu 2014) @@ -241,7 +242,7 @@ def setPsaapPropertiesTestArmInterpTrans(gam, inputV0, inputVDC, params, Nr): rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]], - kf = rxn, kf_T = rxn_T) + kf = rxn, kf_T = rxn_T, rxnBolsig = False) reactionsList.append(reaction) params.reactionsList = reactionsList diff --git a/psaapPropertiesTestArmVioleta.py b/psaapPropertiesTestArmVioleta.py new file mode 100644 index 000000000..9c3d55670 --- /dev/null +++ b/psaapPropertiesTestArmVioleta.py @@ -0,0 +1,431 @@ +import numpy as np +from scipy.interpolate import CubicSpline + +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors +import logging + +class Reaction(object): + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) + +class Diffusivity(object): + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) + +class Mobility(object): + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) + +def setPsaapPropertiesTestArm(gam, inputV0, inputVDC, params, Nr, iSample): + """Sets non-dimensional properties corresponding to Liu 2014 paper. + + Inputs: + gam : Secondary electron emission coefficient + params : chebSolver.modelParams class + + Outputs: None + params data is overwritten using values from Liu 2014. + """ + ################################################################### + # User specified parameters (you may change these if you wish to + # run a different scenario from Liu 2014) + ################################################################### + + # densities + nAr = 3.22e22 # background number density of Ar [1/m^3] (corresponds to p=100 mTorr) + np0 = 8e16 # "nominal" electron density [1/m^3] + + # masses + # me = 9.10938356e-31 # mass of an electron [kg] + # me = 5.489e-4 # mass of an electron [u] + me = 0.511e6 # mass of an electron [eV/c2] + # mAr = 39.948 # mass of an argon atom [u] + # mAr = 39.948 * 1.66054e-27 # mass of an argon atom [kg] + mAr = 37.2158e9 # mass of an electron [eV/c2] + # u = 931.4941e6 # eV/c2 + c = 299792458 # speed of light [m/s] + se = 40 # momentum cross section [m^2] + + # nominal electron energy + e0 = 1.0 # [eV] + + # pressure + p = 133.3224*1.5 # [J/m^3] *1.5 to convert it to energy (1 Torr) + + # gas energy at the wall + Tg0 = 0.038778 # 3/2*300K*kB ~ (p0 - nT[:,0])/ntot + + # characteristics of driving voltage + V0 = inputV0 # amplitude of driving voltage [V] + verticalShift = inputVDC # DC voltage (vertical shift in driving voltage) + tau = (1./13.6e6) # period of driving voltage [s] + L = 2.00*0.005 # half-gap-width [m] (gap width is 2.54cm) + electrodeArea = np.pi*0.05**2 # electrode area [m^2] (electrode diameter = 0.1 m) + + # transport parameters + nmue = 9.66e21 # argon number density times electron mobility [1/(V*cm*s)] + nmui = 4.65e19 # argon number density times ion mobility [1/(V*cm*s)] + nmum = 1 / (np.sqrt(16.0 * (mAr + mAr) * 300 * 8.62e-5 * c**2 + / (3.0 * np.pi * mAr * mAr)) * se * mAr * 1.6e-19 / c**2) + nDe = 3.86e22 # argon number density times electron diffusivity [1/(cm*s)] + nDi = 2.07e18 # argon number density times ion diffusivity [1/(cm*s)] + nDm = 2.42e18 # argon number density times metastable diffusivity [1/(cm*s)] + + # reaction parameters (NB: k_i = Ck*exp(-A/Te)) + #Ck = np.array([1.235e-7,3.712e-8,2.05e-7,1.818e-9,2e-7]) # pre-exponential factors [cm^3/s] + #A = np.array([18.687,15.06,4.95,2.14,0.0]) # activation temperature [eV] + #dH = np.array([15.7,11.56,4.14,-11.56,0.0]) # energy lost per electron due to ionization rxn [eV] + + #Ck = np.array([1.235e-7,3.712e-8,2.05e-7,1.818e-9]) # pre-exponential factors [cm^3/s] + #A = np.array([18.687,15.06,4.95,2.14]) # activation temperature [eV] + #dH = np.array([15.7,11.56,4.14,-11.56]) # energy lost per electron due to ionization rxn [eV] + + #Ck = np.array([1.235e-7,0.0,0.0,0.0,0.0]) # pre-exponential factors [cm^3/s] + #A = np.array([18.687,15.06,4.95,2.14,0.0]) # activation temperature [eV] + #dH = np.array([15.7,0.0,0.0,0.0,0.0]) # energy lost per electron due to ionization rxn [eV] + + # nominal + # Ck = np.array([1.235e-7,3.712e-8,2.05e-7,1.818e-9,2e-7,6.2e-10,3.0e-15,1.1e-31]) # pre-exponential factors [cm^3/s] + Ck = np.array([1.235e-7,3.712e-8,2.05e-7,4.0e-13,2e-7,5.0e-10,2.5e-15]) # pre-exponential factors [cm^3/s] + + # slow excitation rate + #Ck = np.array([1.235e-7,0.5*3.712e-8,2.05e-7,1.818e-9,2e-7,6.2e-10,3.0e-15,1.1e-31]) # pre-exponential factors [cm^3/s] + + ## fast excitation rate + #Ck = np.array([1.235e-7,2.0*3.712e-8,2.05e-7,1.818e-9,2e-7,6.2e-10,3.0e-15,1.1e-31]) # pre-exponential factors [cm^3/s] + + # A = np.array([18.687,15.06,4.95,2.14,0.0,0.0,0.0,0.0]) # activation temperature [eV] + # dH = np.array([15.7,11.56,4.14,-11.56,0.0,0.0,0.0,0.0]) # energy lost per electron due to ionization rxn [eV] + # dEps = np.array([0.0,15.7,11.56,0.0]) + + A = np.array([18.687,15.06,4.95,0.0,0.0,0.0,0.0]) # activation temperature [eV] + dH = np.array([15.76,11.56,4.2,0.0,-11.56,-7.56,-11.56]) # energy lost per electron due to ionization rxn [eV] + dEps = np.array([0.0,15.76,11.56,0.0]) + + # BC parameters + # ks = 1.19e7 # electron recombination rate [cm/s] + ks = 1.366109824889323e7 # electron recombination rate [cm/s/eV] + + ################################################################### + # Constants of nature (probably shouldn't change unless you have + # root privileges on universe) + ################################################################### + qe = 1.6e-19 # unit charge [C] + eps0 = 8.86e-12 # permittivity of free space [F/m] + kB = 1.38e-23 # Boltzmann constant [J/K] + # kB = 8.62e−5 # Boltzmann constant [eV/K] + + + ################################################################### + # Calculate non-dimensional parameters + ################################################################### + + # 1) Convert input units to base SI (except eV) + nDe *= 100. # 1/(m*s) + nDi *= 100. # 1/(m*s) + nDm *= 100. + nmue *= 100. # 1/(V*m*s) + nmui *= 100. # 1/(V*m*s) + Ck *= 1e-6 # m^3/s + # Ck[7] *= 1e-6 # Ck[7] is now in m^6/s + ks *= 0.01 # m/s + se *= 1.0e-20 # m^2 + + # 2) Compute "raw" transport parameters + De = nDe/nAr + Di = nDi/nAr + Dm = nDm/nAr + + mue = nmue/nAr + mui = nmui/nAr + mum = nmum/nAr + + # 3) Compute non-dimensional properties required by solver + De = De*tau/(L*L) + Di = Di*tau/(L*L) + Dm = Dm*tau/(L*L) + + mue = mue*V0*tau/(L*L) + mui = mui*V0*tau/(L*L) + mum = mum*V0*tau/(L*L) + + Ck[0:2] = Ck[0:2]*tau*nAr + #Ck[2:5] = Ck[2:5]*tau*np0 + Ck[2:6] = Ck[2:6]*tau*np0 + Ck[6] *= tau*nAr + # Ck[7] *= tau*nAr*nAr + A = A*1.5/e0 # 1.5 to convert from temperature to energy + dH = dH/e0 + qStar = V0/e0 # qe*V0/e0, since e0 in eV, need qe*V0 in eV, which is just V0 in V + alpha = qe*np0*L*L/(V0*eps0) + ks = ks*tau/L + p0 = p/qe/np0 + kappaB = 4.878171165833662*1.6129 # non-dimensional thermal conductivity of background specie + # (2/3)*tau/L**2*Kb/np0/kB, + # where Kb is the thermal conductivity of background specie + + #params.beta = np.array([[2,2,2,1,1],[1,0,1,0,0],[0,1,0,0,0],[0,0,0,1,1]]) + #params.alfa = np.array([[1,1,1,1,1],[0,0,0,0,0],[0,0,1,1,1],[1,1,0,0,0]]) + #params.beta = np.array([[2,1,2,1],[1,0,1,0],[0,1,0,0],[0,0,0,1]], dtype=np.int) + #params.alfa = np.array([[1,1,1,1],[0,0,0,0],[0,0,1,1],[1,1,0,0]], dtype=np.int) + # params.beta = np.array([[2,1,2,1,1,1,0,0],[1,0,1,0,0,1,0,0],[0,1,0,0,0,0,0,0],[0,0,0,1,0,1,2,1]], dtype=np.int64) + # params.alfa = np.array([[1,1,1,1,1,0,0,0],[0,0,0,0,0,0,0,0],[0,0,1,1,1,2,1,1],[1,1,0,0,0,0,1,2]], dtype=np.int64) + params.beta = np.array([[2,1,2,0,1,1,0],[1,0,1,0,0,1,0],[0,1,0,1,0,0,0],[0,0,0,0,1,1,2]], dtype=np.int64) + params.alfa = np.array([[1,1,1,1,1,0,0],[0,0,0,1,0,0,0],[0,0,1,0,1,2,1],[1,1,0,0,0,0,1]], dtype=np.int64) + + # 4) Set values in params class + params.D[0] = De + params.D[1] = Di + params.D[2] = Dm + + params.mu[0] = mue + params.mu[1] = mui + params.mu[2] = mum + + params.A[:] = Ck[:] + params.B[:] = 0.0 + params.C[:] = A[:] + + params.dH[:] = dH[:] + params.dEps[:] = dEps[:] + params.qStar = qStar + params.alpha = alpha + params.ks = ks + params.gam = gam + params.kappaB = kappaB + params.nAronp0 = nAr / np0 + params.p0 = p0 + params.Tg0 = Tg0 + params.EC = 2.0 * me / mAr \ + * np.sqrt(16.0 * (me + mAr) * e0 * c**2 + / (3.0 * np.pi * me * mAr)) * se * nAr * tau + # params.EC = 2.0 * me / mAr * 3.8e9 * tau + + params.verticalShift = verticalShift / V0 + + # Parameters needed to compute the current with dimensions + params.V0Ltau = V0 / (L * tau) + params.V0L = V0 / L + params.LLV0tau = (L*L) / (V0*tau) + params.tauL = L / tau + params.np0 = np0 # "nominal" electron density [1/m^3] + params.qe = qe # unit charge [C] + params.eps0 = eps0 # unit charge [C] + params.eArea = electrodeArea # electrode area [m^2] + + # reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + # f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + # f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + # f"{params.A[3]} * energy**{params.B[3]} * np.exp(-{params.C[3]} / energy)", + # f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + # f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + # f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)", + # f"{params.A[7]} * energy**{params.B[7]} * np.exp(-{params.C[7]} / energy)"] + + # reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + # f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + # f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + # f"{params.A[3]} * (energy**({params.B[3]}-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]} + {params.C[3]}/energy)", + # f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + # f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + # f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)", + # f"{params.A[7]} * (energy**({params.B[7]}-1)) * np.exp(-{params.C[7]}/energy) * ({params.B[7]} + {params.C[7]}/energy)"] + + reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + f"{params.A[3]} * energy**({params.B[3]}-0.5)", + f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)"] + + reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + f"{params.A[3]} * (energy**({params.B[3]}-0.5-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]}-0.5 + {params.C[3]}/energy)", + f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)"] + + reactionExpressionTypelist = np.array([True,True,True,False,False,False,False]) + + reactionsList = [] + LOGFilename = 'interpolationSample%s.log'%str(iSample) + f = open(LOGFilename, 'w') + for i in range(Nr): + if reactionExpressionTypelist[i]: + Nsample = 72 + N300 = 200 + + root_dir = "/g/g92/jbarbere/glowDischarge/toyProblems/timeDomain" + rate_file = "{0:s}/BOLSIGChemistry/reaction300K_{1:d}.dat".format(root_dir, i) + temp_file = "{0:s}/BOLSIGChemistry/reaction300K.Te.dat".format(root_dir) + + #rateCoeff = np.fromfile('/usr/workspace/violetak/BOLSIGSamples/glowDischarge/toyProblems/timeDomain/BOLSIGChemistry/reaction300K_%s.dat' %str(i)) + rateCoeff = np.fromfile(rate_file) + rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + + #Te = np.fromfile('/usr/workspace/violetak/BOLSIGSamples/glowDischarge/toyProblems/timeDomain/BOLSIGChemistry/reaction300K.Te.dat') + Te = np.fromfile(temp_file) + Te = np.reshape(Te,[Nsample, N300]).T[:,iSample] + + # 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) + + logging.basicConfig(filename=LOGFilename) + logging.warning('Interpolation info (Reaction %s):', i) + logging.warning('First non-zero entry in the rate coefficient: %s', I[0][0]) + logging.warning('Monotonicity of the rate coefficient start from entry: %s', lastFalse) + logging.warning('Position of possible duplicates in mean energy array: %s', + TeDuplicateindsForLog[0]) + + # 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("./Rates_%s.pdf" %str(i), dpi=300) + # # plt.xlim((-0.0001,0.0255)) + # plt.show() + else: + rxn = eval("lambda energy :" + reactionExpressionslist[i]) + rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) + + reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]], + rxnBolsig = reactionExpressionTypelist[i], + kf = rxn, kf_T = rxn_T) + reactionsList.append(reaction) + + params.reactionsList = reactionsList + #params.Nr = 1 + + diffList = [] + Te = np.linspace(0, 1000, 10) + De_interp = params.D[0]*np.ones(10) + De_spline = CubicSpline(Te, De_interp) + De_Te_spline = CubicSpline.derivative(De_spline) + diffusivity = Diffusivity(interpolate = False, D_expression = De_spline, D_T_expression = De_Te_spline) + diffList.append(diffusivity) + + Ns = 4 + for i in range(1, Ns): + diffList.append(Diffusivity(interpolate = False)) + + params.diffusivityList = diffList + + muList = [] + Te = np.linspace(0, 1000, 10) + mue_interp = params.mu[0]*np.ones(10) + mue_spline = CubicSpline(Te, mue_interp) + mue_Te_spline = CubicSpline.derivative(mue_spline) + mobility = Mobility(interpolate = False, mu_expression = mue_spline, mu_T_expression = mue_Te_spline) + muList.append(mobility) + + Ns = 4 + for i in range(1, Ns): + muList.append(Mobility(interpolate = False)) + + params.mobilityList = muList + + # 5) Dump to screen + params.print() diff --git a/psaapPropertiesTestJP.py b/psaapPropertiesTestJP.py new file mode 100644 index 000000000..c826f9b0d --- /dev/null +++ b/psaapPropertiesTestJP.py @@ -0,0 +1,396 @@ +import numpy as np +from scipy.interpolate import CubicSpline + +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors + +import logging + +class Reaction(object): + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) + +def setPsaapPropertiesTestJP(gam, inputV0, inputVDC, params, Nr, iSample): + """Sets non-dimensional properties corresponding to Liu 2014 paper. + + Inputs: + gam : Secondary electron emission coefficient + params : chebSolver.modelParams class + + Outputs: None + params data is overwritten using values from Liu 2014. + """ + ################################################################### + # User specified parameters (you may change these if you wish to + # run a different scenario from Liu 2014) + ################################################################### + + # densities + nAr = 3.22e22 # background number density of Ar [1/m^3] (corresponds to p=100 mTorr) + np0 = 8e16 # "nominal" electron density [1/m^3] + + # masses + # me = 9.10938356e-31 # mass of an electron [kg] + # me = 5.489e-4 # mass of an electron [u] + me = 0.511e6 # mass of an electron [eV/c2] + # mAr = 39.948 # mass of an argon atom [u] + # mAr = 39.948 * 1.66054e-27 # mass of an argon atom [kg] + mAr = 37.2158e9 # mass of an electron [eV/c2] + # u = 931.4941e6 # eV/c2 + c = 299792458 # speed of light [m/s] + se = 40 # momentum cross section [A^2] + + # nominal electron energy + e0 = 1.0 # [eV] + + # pressure + p = 133.3224*1.5 # [J/m^3] *1.5 to convert it to energy (1 Torr) + + # gas energy at the wall + Tg0 = 0.038778 # 3/2*300K*kB ~ (p0 - nT[:,0])/ntot + + # characteristics of driving voltage + V0 = inputV0 # amplitude of driving voltage [V] + verticalShift = inputVDC # DC voltage (vertical shift in driving voltage) + tau = (1./13.6e6) # period of driving voltage [s] + L = 2.00*0.005 # half-gap-width [m] (gap width is 2.54cm) + electrodeArea = np.pi*0.05**2 # electrode area [m^2] (electrode diameter = 0.1 m) + + # transport parameters + nmue = 9.66e21 # argon number density times electron mobility [1/(V*cm*s)] + nmui = 4.65e19 # argon number density times ion mobility [1/(V*cm*s)] + #nmum = 1 / (np.sqrt(16.0 * (mAr + mAr) * 300 * 8.62e-5 * c**2 + # / (3.0 * np.pi * mAr * mAr)) * se * mAr * 1.6e-19 / c**2) + nmum = 0.0 + nDe = 3.86e22 # argon number density times electron diffusivity [1/(cm*s)] + nDi = 2.07e18 # argon number density times ion diffusivity [1/(cm*s)] + nDm = 2.42e18 # argon number density times metastable diffusivity [1/(cm*s)] + + # reaction parameters (NB: k_i = Ck*Ee^B*exp(-A/Ee)) + # Ee = 3/2*Te (Te in eV) + # -> k_i = [Ck*(2/3)^B] * Ee^B * exp[-(3/2)*A/Ee] + # nominal + Ck = np.array([1.235e-7,3.712e-8,2.05e-7,4.0e-13,5.0e-10,4.3e-10,2.1e-15,10.0,5.0e-27]) # pre-exponential factors [cm^3/s] + B = np.array([0.0,0.0,0.0,-0.5,0,0.74,0,0,-4.5]) + A = np.array([18.687,15.06,4.95,0.0,0.0,0.0,0.0,0.0,0.0]) # activation temperature [eV] + dH = np.array([15.76,11.56,4.2,0.0,-7.36,-11.56,-11.56,0.0,-4.2]) # energy lost per electron due to ionization rxn [eV] + dEps = np.array([0.0,15.76,11.56,0.0]) + + # BC parameters + # ks = 1.19e7 # electron recombination rate [cm/s] + ks = 1.366109824889323e7 # electron recombination rate [cm/s/eV] + + ################################################################### + # Constants of nature (probably shouldn't change unless you have + # root privileges on universe) + ################################################################### + qe = 1.6e-19 # unit charge [C] + eps0 = 8.86e-12 # permittivity of free space [F/m] + kB = 1.38e-23 # Boltzmann constant [J/K] + # kB = 8.62e−5 # Boltzmann constant [eV/K] + + + ################################################################### + # Calculate non-dimensional parameters + ################################################################### + + # 1) Convert input units to base SI (except eV) + nDe *= 100. # 1/(m*s) + nDi *= 100. # 1/(m*s) + nDm *= 100. + nmue *= 100. # 1/(V*m*s) + nmui *= 100. # 1/(V*m*s) + Ck[0:7] *= 1e-6 # m^3/s + Ck[7] *= 1 # 1/s + Ck[8] *= 1e-12 # m^6/s + ks *= 0.01 # m/s + se *= 1.0e-20 # m^2 + + # 2) Compute "raw" transport parameters + De = nDe/nAr + Di = nDi/nAr + Dm = nDm/nAr + + mue = nmue/nAr + mui = nmui/nAr + mum = nmum/nAr + + # 3) Compute non-dimensional properties required by solver + De = De*tau/(L*L) + Di = Di*tau/(L*L) + Dm = Dm*tau/(L*L) + + mue = mue*V0*tau/(L*L) + mui = mui*V0*tau/(L*L) + mum = mum*V0*tau/(L*L) + + Ck[0:2] = Ck[0:2]*tau*nAr + #Ck[2:5] = Ck[2:5]*tau*np0 + Ck[2:6] = Ck[2:6]*tau*np0 + Ck[6] *= tau*nAr + Ck[7] *= tau + Ck[8] *= tau*np0*np0 + # Ck[7] *= tau*nAr*nAr + A = A*1.5/e0 # 1.5 to convert from temperature to energy + dH = dH/e0 + qStar = V0/e0 # qe*V0/e0, since e0 in eV, need qe*V0 in eV, which is just V0 in V + alpha = qe*np0*L*L/(V0*eps0) + ks = ks*tau/L + p0 = p/qe/np0 + kappaB = 4.878171165833662*1.6129 # non-dimensional thermal conductivity of background specie + # (2/3)*tau/L**2*Kb/np0/kB, + # where Kb is the thermal conductivity of background specie + + #params.beta = np.array([[2,2,2,1,1],[1,0,1,0,0],[0,1,0,0,0],[0,0,0,1,1]]) + #params.alfa = np.array([[1,1,1,1,1],[0,0,0,0,0],[0,0,1,1,1],[1,1,0,0,0]]) + #params.beta = np.array([[2,1,2,1],[1,0,1,0],[0,1,0,0],[0,0,0,1]], dtype=np.int) + #params.alfa = np.array([[1,1,1,1],[0,0,0,0],[0,0,1,1],[1,1,0,0]], dtype=np.int) + # params.beta = np.array([[2,1,2,1,1,1,0,0],[1,0,1,0,0,1,0,0],[0,1,0,0,0,0,0,0],[0,0,0,1,0,1,2,1]], dtype=np.int64) + # params.alfa = np.array([[1,1,1,1,1,0,0,0],[0,0,0,0,0,0,0,0],[0,0,1,1,1,2,1,1],[1,1,0,0,0,0,1,2]], dtype=np.int64) + params.beta = np.array([[2,1,2,0,1,1,0,0,1],[1,0,1,0,1,0,0,0,0],[0,1,0,1,0,0,0,0,1],[0,0,0,0,1,1,2,1,0]], dtype=np.int64) + params.alfa = np.array([[1,1,1,1,0,1,0,0,2],[0,0,0,1,0,0,0,0,1],[0,0,1,0,2,1,1,1,0],[1,1,0,0,0,0,1,0,0]], dtype=np.int64) + # Rxn1: E + AR -> 2E + AR+ + # Rxn2: E + AR -> E + AR* + # Rxn3: E + AR* -> 2E + AR+ + # Rxn4: E + AR+ -> AR* + # Rxn5: 2AR* -> E + AR+ + AR + # Rxn6: E + AR* -> E + AR + # Rxn7: AR* + AR -> 2AR + # Rxn8: AR* -> AR + # Rxn9: 2E + AR+ -> E + AR* + + # 4) Set values in params class + params.D[0] = De + params.D[1] = Di + params.D[2] = Dm + + params.mu[0] = mue + params.mu[1] = mui + params.mu[2] = mum + + params.A[:] = Ck[:] + params.B[:] = B[:] + params.C[:] = A[:] + + # Account for the 2/3 term to convert from electron temperature to electron energy + for i in range(len(params.A)): + params.A[i] *= (2/3)**(params.B[i]) + + params.dH[:] = dH[:] + params.dEps[:] = dEps[:] + params.qStar = qStar + params.alpha = alpha + params.ks = ks + params.gam = gam + params.kappaB = kappaB + params.nAronp0 = nAr / np0 + params.p0 = p0 + params.Tg0 = Tg0 + params.EC = 2.0 * me / mAr \ + * np.sqrt(16.0 * (me + mAr) * e0 * c**2 + / (3.0 * np.pi * me * mAr)) * se * nAr * tau + # params.EC = 2.0 * me / mAr * 3.8e9 * tau + + params.verticalShift = verticalShift / V0 + + # Parameters needed to compute the current with dimensions + params.V0Ltau = V0 / (L * tau) + params.V0L = V0 / L + params.LLV0tau = (L*L) / (V0*tau) + params.tauL = L / tau + params.np0 = np0 # "nominal" electron density [1/m^3] + params.qe = qe # unit charge [C] + params.eps0 = eps0 # unit charge [C] + params.eArea = electrodeArea # electrode area [m^2] + + # reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + # f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + # f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + # f"{params.A[3]} * energy**{params.B[3]} * np.exp(-{params.C[3]} / energy)", + # f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + # f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + # f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)", + # f"{params.A[7]} * energy**{params.B[7]} * np.exp(-{params.C[7]} / energy)"] + + # reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + # f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + # f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + # f"{params.A[3]} * (energy**({params.B[3]}-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]} + {params.C[3]}/energy)", + # f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + # f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + # f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)", + # f"{params.A[7]} * (energy**({params.B[7]}-1)) * np.exp(-{params.C[7]}/energy) * ({params.B[7]} + {params.C[7]}/energy)"] + + reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + f"{params.A[3]} * energy**{params.B[3]} * np.exp(-{params.C[3]} / energy)", + f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)", + f"{params.A[7]} * energy**{params.B[7]} * np.exp(-{params.C[7]} / energy)", + f"{params.A[8]} * energy**{params.B[8]} * np.exp(-{params.C[8]} / energy)"] + + reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + f"{params.A[3]} * (energy**({params.B[3]}-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]} + {params.C[3]}/energy)", + f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)", + f"{params.A[7]} * (energy**({params.B[7]}-1)) * np.exp(-{params.C[7]}/energy) * ({params.B[7]} + {params.C[7]}/energy)", + f"{params.A[8]} * (energy**({params.B[8]}-1)) * np.exp(-{params.C[8]}/energy) * ({params.B[8]} + {params.C[8]}/energy)"] + + reactionExpressionTypelist = np.array([True,True,True,False,False,False,False,False,False]) + + reactionsList = [] + LOGFilename = 'interpolationSample%s.log'%str(iSample) + f = open(LOGFilename, 'w') + # 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 = 72 + N300 = 200 + + root_dir = "/g/g92/jbarbere/glowDischarge/toyProblems/timeDomain" + rate_file = "{0:s}/BOLSIGChemistry/reaction300K_{1:d}.dat".format(root_dir, i) + temp_file = "{0:s}/BOLSIGChemistry/reaction300K.Te.dat".format(root_dir) + + #rateCoeff = np.fromfile('/usr/workspace/violetak/BOLSIGSamples/glowDischarge/toyProblems/timeDomain/BOLSIGChemistry/reaction300K_%s.dat' %str(i)) + rateCoeff = np.fromfile(rate_file) + rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + + #Te = np.fromfile('/usr/workspace/violetak/BOLSIGSamples/glowDischarge/toyProblems/timeDomain/BOLSIGChemistry/reaction300K.Te.dat') + Te = np.fromfile(temp_file) + Te = np.reshape(Te,[Nsample, N300]).T[:,iSample] + + # 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) + + logging.basicConfig(filename=LOGFilename) + logging.warning('Interpolation info (Reaction %s):', i) + logging.warning('First non-zero entry in the rate coefficient: %s', I[0][0]) + logging.warning('Monotonicity of the rate coefficient start from entry: %s', lastFalse) + logging.warning('Position of possible duplicates in mean energy array: %s', + TeDuplicateindsForLog[0]) + + # 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("./Rates_%s.pdf" %str(i), dpi=300) + # # plt.xlim((-0.0001,0.0255)) + # plt.show() + else: + rxn = eval("lambda energy :" + reactionExpressionslist[i]) + rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) + + reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]], + rxnBolsig = reactionExpressionTypelist[i], + kf = rxn, kf_T = rxn_T) + reactionsList.append(reaction) + + params.reactionsList = reactionsList + #params.Nr = 1 + + # 5) Dump to screen + params.print() diff --git a/psaapPropertiesTestJP_Arrhenius.py b/psaapPropertiesTestJP_Arrhenius.py new file mode 100755 index 000000000..84b12f6e5 --- /dev/null +++ b/psaapPropertiesTestJP_Arrhenius.py @@ -0,0 +1,373 @@ +import numpy as np +from scipy.interpolate import CubicSpline +import csv +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors + +import logging + +class Reaction(object): + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) + +def setPsaapPropertiesTestJP_Arrhenius(gam, inputV0, inputVDC, params, Nr, iSample): + """Sets non-dimensional properties corresponding to Liu 2014 paper. + + Inputs: + gam : Secondary electron emission coefficient + params : chebSolver.modelParams class + + Outputs: None + params data is overwritten using values from Liu 2014. + """ + ################################################################### + # User specified parameters (you may change these if you wish to + # run a different scenario from Liu 2014) + ################################################################### + + # densities + nAr = 3.22e22 # background number density of Ar [1/m^3] (corresponds to p=100 mTorr) + np0 = 8e16 # "nominal" electron density [1/m^3] + + # masses + # me = 9.10938356e-31 # mass of an electron [kg] + # me = 5.489e-4 # mass of an electron [u] + me = 0.511e6 # mass of an electron [eV/c2] + # mAr = 39.948 # mass of an argon atom [u] + # mAr = 39.948 * 1.66054e-27 # mass of an argon atom [kg] + mAr = 37.2158e9 # mass of an electron [eV/c2] + # u = 931.4941e6 # eV/c2 + c = 299792458 # speed of light [m/s] + se = 40 # momentum cross section [A^2] + + # nominal electron energy + e0 = 1.0 # [eV] + + # pressure + p = 133.3224*1.5 # [J/m^3] *1.5 to convert it to energy (1 Torr) + + # gas energy at the wall + Tg0 = 0.038778 # 3/2*300K*kB ~ (p0 - nT[:,0])/ntot + + # characteristics of driving voltage + V0 = inputV0 # amplitude of driving voltage [V] + verticalShift = inputVDC # DC voltage (vertical shift in driving voltage) + tau = (1./13.56e6) # period of driving voltage [s] + L = 2.00*0.005 # half-gap-width [m] (gap width is 2 cm) + electrodeArea = np.pi*0.05**2 # electrode area [m^2] (electrode diameter = 0.1 m) + + # transport parameters + nmue = 9.66e21 # argon number density times electron mobility [1/(V*cm*s)] + #nmui = 4.65e19 # argon number density times ion mobility [1/(V*cm*s)] + #nmum = 1 / (np.sqrt(16.0 * (mAr + mAr) * 300 * 8.62e-5 * c**2 + # / (3.0 * np.pi * mAr * mAr)) * se * mAr * 1.6e-19 / c**2) + nmum = 0.0 + nmui = 8.0e19 + #nmum = 9.35e19 + nDe = 3.86e22 # argon number density times electron diffusivity [1/(cm*s)] + nDi = 2.07e18 # argon number density times ion diffusivity [1/(cm*s)] + nDm = 2.42e18 # argon number density times metastable diffusivity [1/(cm*s)] + #nDm = 3.914e20 + + # reaction parameters (NB: k_i = Ck*Ee^B*exp(-A/Ee)) + # Ee = 3/2*Te (Te in eV) + # -> k_i = [Ck*(2/3)^B] * Ee^B * exp[-(3/2)*A/Ee] + # nominal + Ck = np.array([1.235e-7,3.712e-8,2.05e-7,4.0e-13,5.0e-10,4.3e-10,2.1e-15,10.0,5.0e-27]) # pre-exponential factors [cm^3/s] + B = np.array([0.0,0.0,0.0,-0.5,0.0,0.74,0.0,0.0,-4.5]) + A = np.array([18.687,15.06,4.95,0.0,0.0,0.0,0.0,0.0,0.0]) # activation temperature [eV] + dH = np.array([15.76,11.56,4.2,0.0,-7.56,-11.56,-11.56,0.0,-4.2]) # energy lost per electron due to ionization rxn [eV] + dEps = np.array([0.0,15.76,11.56,0.0]) + + # BC parameters + # ks = 1.19e7 # electron recombination rate [cm/s] + ks = 1.366109824889323e7 # electron recombination rate [cm/s/eV] + + ################################################################### + # Constants of nature (probably shouldn't change unless you have + # root privileges on universe) + ################################################################### + qe = 1.6e-19 # unit charge [C] + eps0 = 8.86e-12 # permittivity of free space [F/m] + kB = 1.38e-23 # Boltzmann constant [J/K] + # kB = 8.62e−5 # Boltzmann constant [eV/K] + + + ################################################################### + # Calculate non-dimensional parameters + ################################################################### + + # 1) Convert input units to base SI (except eV) + nDe *= 100. # 1/(m*s) + nDi *= 100. # 1/(m*s) + nDm *= 100. + nmue *= 100. # 1/(V*m*s) + nmui *= 100. # 1/(V*m*s) + Ck[0:7] *= 1e-6 # m^3/s + Ck[7] *= 1 # 1/s + Ck[8] *= 1e-12 # m^6/s + # Ck[7] *= 1e-6 # Ck[7] is now in m^6/s + ks *= 0.01 # m/s + se *= 1.0e-20 # m^2 + + # 2) Compute "raw" transport parameters + De = nDe/nAr + Di = nDi/nAr + Dm = nDm/nAr + + mue = nmue/nAr + mui = nmui/nAr + mum = nmum/nAr + + # 3) Compute non-dimensional properties required by solver + De = De*tau/(L*L) + Di = Di*tau/(L*L) + Dm = Dm*tau/(L*L) + + mue = mue*V0*tau/(L*L) + mui = mui*V0*tau/(L*L) + mum = mum*V0*tau/(L*L) + + Ck[0:2] = Ck[0:2]*tau*nAr + #Ck[2:5] = Ck[2:5]*tau*np0 + Ck[2:6] = Ck[2:6]*tau*np0 + Ck[6] *= tau*nAr + Ck[7] *= tau + Ck[8] *= tau*np0*np0 + # Ck[7] *= tau*nAr*nAr + A = A*1.5/e0 # 1.5 to convert from temperature to energy + dH = dH/e0 + qStar = V0/e0 # qe*V0/e0, since e0 in eV, need qe*V0 in eV, which is just V0 in V + alpha = qe*np0*L*L/(V0*eps0) + ks = ks*tau/L + p0 = p/qe/np0 + kappaB = 4.878171165833662*1.6129 # non-dimensional thermal conductivity of background specie + # (2/3)*tau/L**2*Kb/np0/kB, + # where Kb is the thermal conductivity of background specie + + #params.beta = np.array([[2,2,2,1,1],[1,0,1,0,0],[0,1,0,0,0],[0,0,0,1,1]]) + #params.alfa = np.array([[1,1,1,1,1],[0,0,0,0,0],[0,0,1,1,1],[1,1,0,0,0]]) + #params.beta = np.array([[2,1,2,1],[1,0,1,0],[0,1,0,0],[0,0,0,1]], dtype=np.int) + #params.alfa = np.array([[1,1,1,1],[0,0,0,0],[0,0,1,1],[1,1,0,0]], dtype=np.int) + # params.beta = np.array([[2,1,2,1,1,1,0,0],[1,0,1,0,0,1,0,0],[0,1,0,0,0,0,0,0],[0,0,0,1,0,1,2,1]], dtype=np.int64) + # params.alfa = np.array([[1,1,1,1,1,0,0,0],[0,0,0,0,0,0,0,0],[0,0,1,1,1,2,1,1],[1,1,0,0,0,0,1,2]], dtype=np.int64) + params.beta = np.array([[2,1,2,0,1,1,0,0,1],[1,0,1,0,1,0,0,0,0],[0,1,0,1,0,0,0,0,1],[0,0,0,0,1,1,2,1,0]], dtype=np.int64) + params.alfa = np.array([[1,1,1,1,0,1,0,0,2],[0,0,0,1,0,0,0,0,1],[0,0,1,0,2,1,1,1,0],[1,1,0,0,0,0,1,0,0]], dtype=np.int64) + # Rxn1: E + AR -> 2E + AR+ + # Rxn2: E + AR -> E + AR* + # Rxn3: E + AR* -> 2E + AR+ + # Rxn4: E + AR+ -> AR* + # Rxn5: 2AR* -> E + AR+ + AR + # Rxn6: E + AR* -> E + AR + # Rxn7: AR* + AR -> 2AR + # Rxn8: AR* -> AR + # Rxn9: 2E + AR+ -> E + AR* + + # 4) Set values in params class + params.D[0] = De + params.D[1] = Di + params.D[2] = Dm + + params.mu[0] = mue + params.mu[1] = mui + params.mu[2] = mum + + params.A[:] = Ck[:] + params.B[:] = B[:] + params.C[:] = A[:] + + # Account for the 2/3 term to convert from electron temperature to electron energy + for i in range(len(params.A)): + params.A[i] *= (2/3)**(params.B[i]) + + params.dH[:] = dH[:] + params.dEps[:] = dEps[:] + params.qStar = qStar + params.alpha = alpha + params.ks = ks + params.gam = gam + params.kappaB = kappaB + params.nAronp0 = nAr / np0 + params.p0 = p0 + params.Tg0 = Tg0 + params.EC = 2.0 * me / mAr \ + * np.sqrt(16.0 * (me + mAr) * e0 * c**2 + / (3.0 * np.pi * me * mAr)) * se * nAr * tau + # params.EC = 2.0 * me / mAr * 3.8e9 * tau + + params.verticalShift = verticalShift / V0 + + # Parameters needed to compute the current with dimensions + params.V0Ltau = V0 / (L * tau) + params.V0L = V0 / L + params.LLV0tau = (L*L) / (V0*tau) + params.tauL = L / tau + params.np0 = np0 # "nominal" electron density [1/m^3] + params.qe = qe # unit charge [C] + params.eps0 = eps0 # unit charge [C] + params.eArea = electrodeArea # electrode area [m^2] + + # reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + # f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + # f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + # f"{params.A[3]} * energy**{params.B[3]} * np.exp(-{params.C[3]} / energy)", + # f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + # f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + # f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)", + # f"{params.A[7]} * energy**{params.B[7]} * np.exp(-{params.C[7]} / energy)"] + + # reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + # f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + # f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + # f"{params.A[3]} * (energy**({params.B[3]}-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]} + {params.C[3]}/energy)", + # f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + # f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + # f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)", + # f"{params.A[7]} * (energy**({params.B[7]}-1)) * np.exp(-{params.C[7]}/energy) * ({params.B[7]} + {params.C[7]}/energy)"] + + reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + f"{params.A[3]} * energy**{params.B[3]} * np.exp(-{params.C[3]} / energy)", + f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)", + f"{params.A[7]} * energy**{params.B[7]} * np.exp(-{params.C[7]} / energy)", + f"{params.A[8]} * energy**{params.B[8]} * np.exp(-{params.C[8]} / energy)"] + + reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + f"{params.A[3]} * (energy**({params.B[3]}-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]} + {params.C[3]}/energy)", + f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)", + f"{params.A[7]} * (energy**({params.B[7]}-1)) * np.exp(-{params.C[7]}/energy) * ({params.B[7]} + {params.C[7]}/energy)", + f"{params.A[8]} * (energy**({params.B[8]}-1)) * np.exp(-{params.C[8]}/energy) * ({params.B[8]} + {params.C[8]}/energy)"] + + reactionExpressionTypelist = np.array([False,False,False,False,False,False,False,False,False]) + + reactionsList = [] + LOGFilename = 'interpolationSample%s.log'%str(iSample) + f = open(LOGFilename, 'w') + # 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}/BOLSIGChemistry_NominalRates/reaction300K_{1:d}.txt".format(root_dir, i), 'r') + temp_file = open("{0:s}/BOLSIGChemistry_NominalRates/reaction300K_Te.txt".format(root_dir), 'r') + + #rateCoeff = np.fromfile(rate_file) + rateCoeff = np.genfromtxt(rate_file) + rate_file.close() + rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + + #Te = np.fromfile(temp_file) + Te = np.genfromtxt(temp_file) + temp_file.close() + Te = np.reshape(Te,[Nsample, N300]).T[:,iSample] + + # 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) + + logging.basicConfig(filename=LOGFilename) + logging.warning('Interpolation info (Reaction %s):', i) + logging.warning('First non-zero entry in the rate coefficient: %s', I[0][0]) + logging.warning('Monotonicity of the rate coefficient start from entry: %s', lastFalse) + logging.warning('Position of possible duplicates in mean energy array: %s', + TeDuplicateindsForLog[0]) + + else: + rxn = eval("lambda energy :" + reactionExpressionslist[i]) + rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) + + reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]], + rxnBolsig = reactionExpressionTypelist[i], + kf = rxn, kf_T = rxn_T) + reactionsList.append(reaction) + + params.reactionsList = reactionsList + + # 5) Dump to screen + params.print() diff --git a/psaapPropertiesTestJP_Nominal.py b/psaapPropertiesTestJP_Nominal.py new file mode 100755 index 000000000..8b8ce532b --- /dev/null +++ b/psaapPropertiesTestJP_Nominal.py @@ -0,0 +1,371 @@ +import numpy as np +from scipy.interpolate import CubicSpline +import csv +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors + +import logging + +class Reaction(object): + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) + +def setPsaapPropertiesTestJP_Nominal(gam, inputV0, inputVDC, params, Nr, iSample): + """Sets non-dimensional properties corresponding to Liu 2014 paper. + + Inputs: + gam : Secondary electron emission coefficient + params : chebSolver.modelParams class + + Outputs: None + params data is overwritten using values from Liu 2014. + """ + ################################################################### + # User specified parameters (you may change these if you wish to + # run a different scenario from Liu 2014) + ################################################################### + + # densities + nAr = 3.22e22 # background number density of Ar [1/m^3] (corresponds to p=100 mTorr) + np0 = 8e16 # "nominal" electron density [1/m^3] + + # masses + # me = 9.10938356e-31 # mass of an electron [kg] + # me = 5.489e-4 # mass of an electron [u] + me = 0.511e6 # mass of an electron [eV/c2] + # mAr = 39.948 # mass of an argon atom [u] + # mAr = 39.948 * 1.66054e-27 # mass of an argon atom [kg] + mAr = 37.2158e9 # mass of an electron [eV/c2] + # u = 931.4941e6 # eV/c2 + c = 299792458 # speed of light [m/s] + se = 40 # momentum cross section [A^2] + + # nominal electron energy + e0 = 1.0 # [eV] + + # pressure + p = 133.3224*1.5 # [J/m^3] *1.5 to convert it to energy (1 Torr) + + # gas energy at the wall + Tg0 = 0.038778 # 3/2*300K*kB ~ (p0 - nT[:,0])/ntot + + # characteristics of driving voltage + V0 = inputV0 # amplitude of driving voltage [V] + verticalShift = inputVDC # DC voltage (vertical shift in driving voltage) + tau = (1./13.56e6) # period of driving voltage [s] + L = 2.00*0.005 # half-gap-width [m] (gap width is 2 cm) + electrodeArea = np.pi*0.05**2 # electrode area [m^2] (electrode diameter = 0.1 m) + + # transport parameters + nmue = 9.66e21 # argon number density times electron mobility [1/(V*cm*s)] + #nmui = 4.65e19 # argon number density times ion mobility [1/(V*cm*s)] + #nmum = 1 / (np.sqrt(16.0 * (mAr + mAr) * 300 * 8.62e-5 * c**2 + # / (3.0 * np.pi * mAr * mAr)) * se * mAr * 1.6e-19 / c**2) + nmum = 0.0 + nmui = 8.0e19 + #nmum = 9.35e19 + nDe = 3.86e22 # argon number density times electron diffusivity [1/(cm*s)] + nDi = 2.07e18 # argon number density times ion diffusivity [1/(cm*s)] + nDm = 2.42e18 # argon number density times metastable diffusivity [1/(cm*s)] + #nDm = 3.914e20 + + # reaction parameters (NB: k_i = Ck*Ee^B*exp(-A/Ee)) + # Ee = 3/2*Te (Te in eV) + # -> k_i = [Ck*(2/3)^B] * Ee^B * exp[-(3/2)*A/Ee] + # nominal + Ck = np.array([1.235e-7,3.712e-8,2.05e-7,4.0e-13,5.0e-10,4.3e-10,2.1e-15,10.0,5.0e-27]) # pre-exponential factors [cm^3/s] + B = np.array([0.0,0.0,0.0,-0.5,0.0,0.74,0.0,0.0,-4.5]) + A = np.array([18.687,15.06,4.95,0.0,0.0,0.0,0.0,0.0,0.0]) # activation temperature [eV] + dH = np.array([15.76,11.56,4.2,0.0,-7.56,-11.56,-11.56,0.0,-4.2]) # energy lost per electron due to ionization rxn [eV] + dEps = np.array([0.0,15.76,11.56,0.0]) + + # BC parameters + # ks = 1.19e7 # electron recombination rate [cm/s] + ks = 1.366109824889323e7 # electron recombination rate [cm/s/eV] + + ################################################################### + # Constants of nature (probably shouldn't change unless you have + # root privileges on universe) + ################################################################### + qe = 1.6e-19 # unit charge [C] + eps0 = 8.86e-12 # permittivity of free space [F/m] + kB = 1.38e-23 # Boltzmann constant [J/K] + # kB = 8.62e−5 # Boltzmann constant [eV/K] + + + ################################################################### + # Calculate non-dimensional parameters + ################################################################### + + # 1) Convert input units to base SI (except eV) + nDe *= 100. # 1/(m*s) + nDi *= 100. # 1/(m*s) + nDm *= 100. + nmue *= 100. # 1/(V*m*s) + nmui *= 100. # 1/(V*m*s) + Ck[0:7] *= 1e-6 # m^3/s + Ck[7] *= 1 # 1/s + Ck[8] *= 1e-12 # m^6/s + # Ck[7] *= 1e-6 # Ck[7] is now in m^6/s + ks *= 0.01 # m/s + se *= 1.0e-20 # m^2 + + # 2) Compute "raw" transport parameters + De = nDe/nAr + Di = nDi/nAr + Dm = nDm/nAr + + mue = nmue/nAr + mui = nmui/nAr + mum = nmum/nAr + + # 3) Compute non-dimensional properties required by solver + De = De*tau/(L*L) + Di = Di*tau/(L*L) + Dm = Dm*tau/(L*L) + + mue = mue*V0*tau/(L*L) + mui = mui*V0*tau/(L*L) + mum = mum*V0*tau/(L*L) + + Ck[0:2] = Ck[0:2]*tau*nAr + Ck[2:6] = Ck[2:6]*tau*np0 + Ck[6] *= tau*nAr + Ck[7] *= tau + Ck[8] *= tau*np0*np0 + A = A*1.5/e0 # 1.5 to convert from temperature to energy + dH = dH/e0 + qStar = V0/e0 # qe*V0/e0, since e0 in eV, need qe*V0 in eV, which is just V0 in V + alpha = qe*np0*L*L/(V0*eps0) + ks = ks*tau/L + p0 = p/qe/np0 + kappaB = 4.878171165833662*1.6129 # non-dimensional thermal conductivity of background specie + # (2/3)*tau/L**2*Kb/np0/kB, + # where Kb is the thermal conductivity of background specie + + #params.beta = np.array([[2,2,2,1,1],[1,0,1,0,0],[0,1,0,0,0],[0,0,0,1,1]]) + #params.alfa = np.array([[1,1,1,1,1],[0,0,0,0,0],[0,0,1,1,1],[1,1,0,0,0]]) + #params.beta = np.array([[2,1,2,1],[1,0,1,0],[0,1,0,0],[0,0,0,1]], dtype=np.int) + #params.alfa = np.array([[1,1,1,1],[0,0,0,0],[0,0,1,1],[1,1,0,0]], dtype=np.int) + # params.beta = np.array([[2,1,2,1,1,1,0,0],[1,0,1,0,0,1,0,0],[0,1,0,0,0,0,0,0],[0,0,0,1,0,1,2,1]], dtype=np.int64) + # params.alfa = np.array([[1,1,1,1,1,0,0,0],[0,0,0,0,0,0,0,0],[0,0,1,1,1,2,1,1],[1,1,0,0,0,0,1,2]], dtype=np.int64) + params.beta = np.array([[2,1,2,0,1,1,0,0,1],[1,0,1,0,1,0,0,0,0],[0,1,0,1,0,0,0,0,1],[0,0,0,0,1,1,2,1,0]], dtype=np.int64) + params.alfa = np.array([[1,1,1,1,0,1,0,0,2],[0,0,0,1,0,0,0,0,1],[0,0,1,0,2,1,1,1,0],[1,1,0,0,0,0,1,0,0]], dtype=np.int64) + # Rxn1: E + AR -> 2E + AR+ + # Rxn2: E + AR -> E + AR* + # Rxn3: E + AR* -> 2E + AR+ + # Rxn4: E + AR+ -> AR* + # Rxn5: 2AR* -> E + AR+ + AR + # Rxn6: E + AR* -> E + AR + # Rxn7: AR* + AR -> 2AR + # Rxn8: AR* -> AR + # Rxn9: 2E + AR+ -> E + AR* + + # 4) Set values in params class + params.D[0] = De + params.D[1] = Di + params.D[2] = Dm + + params.mu[0] = mue + params.mu[1] = mui + params.mu[2] = mum + + params.A[:] = Ck[:] + params.B[:] = B[:] + params.C[:] = A[:] + + # Account for the 2/3 term to convert from electron temperature to electron energy + for i in range(len(params.A)): + params.A[i] *= (2/3)**(params.B[i]) + + params.dH[:] = dH[:] + params.dEps[:] = dEps[:] + params.qStar = qStar + params.alpha = alpha + params.ks = ks + params.gam = gam + params.kappaB = kappaB + params.nAronp0 = nAr / np0 + params.p0 = p0 + params.Tg0 = Tg0 + params.EC = 2.0 * me / mAr \ + * np.sqrt(16.0 * (me + mAr) * e0 * c**2 + / (3.0 * np.pi * me * mAr)) * se * nAr * tau + # params.EC = 2.0 * me / mAr * 3.8e9 * tau + + params.verticalShift = verticalShift / V0 + + # Parameters needed to compute the current with dimensions + params.V0Ltau = V0 / (L * tau) + params.V0L = V0 / L + params.LLV0tau = (L*L) / (V0*tau) + params.tauL = L / tau + params.np0 = np0 # "nominal" electron density [1/m^3] + params.qe = qe # unit charge [C] + params.eps0 = eps0 # unit charge [C] + params.eArea = electrodeArea # electrode area [m^2] + + # reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + # f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + # f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + # f"{params.A[3]} * energy**{params.B[3]} * np.exp(-{params.C[3]} / energy)", + # f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + # f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + # f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)", + # f"{params.A[7]} * energy**{params.B[7]} * np.exp(-{params.C[7]} / energy)"] + + # reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + # f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + # f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + # f"{params.A[3]} * (energy**({params.B[3]}-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]} + {params.C[3]}/energy)", + # f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + # f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + # f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)", + # f"{params.A[7]} * (energy**({params.B[7]}-1)) * np.exp(-{params.C[7]}/energy) * ({params.B[7]} + {params.C[7]}/energy)"] + + reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + f"{params.A[3]} * energy**{params.B[3]} * np.exp(-{params.C[3]} / energy)", + f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)", + f"{params.A[7]} * energy**{params.B[7]} * np.exp(-{params.C[7]} / energy)", + f"{params.A[8]} * energy**{params.B[8]} * np.exp(-{params.C[8]} / energy)"] + + reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + f"{params.A[3]} * (energy**({params.B[3]}-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]} + {params.C[3]}/energy)", + f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)", + f"{params.A[7]} * (energy**({params.B[7]}-1)) * np.exp(-{params.C[7]}/energy) * ({params.B[7]} + {params.C[7]}/energy)", + f"{params.A[8]} * (energy**({params.B[8]}-1)) * np.exp(-{params.C[8]}/energy) * ({params.B[8]} + {params.C[8]}/energy)"] + + reactionExpressionTypelist = np.array([True,True,True,False,False,False,False,False,False]) + + reactionsList = [] + LOGFilename = 'interpolationSample%s.log'%str(iSample) + f = open(LOGFilename, 'w') + # 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}/BOLSIGChemistry_NominalRates/reaction300K_{1:d}.txt".format(root_dir, i), 'r') + temp_file = open("{0:s}/BOLSIGChemistry_NominalRates/reaction300K_Te.txt".format(root_dir), 'r') + + #rateCoeff = np.fromfile(rate_file) + rateCoeff = np.genfromtxt(rate_file) + rate_file.close() + rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + + #Te = np.fromfile(temp_file) + Te = np.genfromtxt(temp_file) + temp_file.close() + Te = np.reshape(Te,[Nsample, N300]).T[:,iSample] + + # 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) + + logging.basicConfig(filename=LOGFilename) + logging.warning('Interpolation info (Reaction %s):', i) + logging.warning('First non-zero entry in the rate coefficient: %s', I[0][0]) + logging.warning('Monotonicity of the rate coefficient start from entry: %s', lastFalse) + logging.warning('Position of possible duplicates in mean energy array: %s', + TeDuplicateindsForLog[0]) + + else: + rxn = eval("lambda energy :" + reactionExpressionslist[i]) + rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) + + reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]], + rxnBolsig = reactionExpressionTypelist[i], + kf = rxn, kf_T = rxn_T) + reactionsList.append(reaction) + + params.reactionsList = reactionsList + + # 5) Dump to screen + params.print() diff --git a/psaapPropertiesWithSampling.py b/psaapPropertiesWithSampling.py new file mode 100644 index 000000000..ce6a6977f --- /dev/null +++ b/psaapPropertiesWithSampling.py @@ -0,0 +1,382 @@ +import numpy as np +from scipy.interpolate import CubicSpline + +#import matplotlib.pyplot as plt +#import matplotlib.colors as mcolors + +import logging + +class Reaction(object): + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) + +def setPsaapPropertiesWithSampling(gam, inputV0, inputVDC, params, Nr, iSample): + """Sets non-dimensional properties corresponding to Liu 2014 paper. + + Inputs: + gam : Secondary electron emission coefficient + params : chebSolver.modelParams class + + Outputs: None + params data is overwritten using values from Liu 2014. + """ + ################################################################### + # User specified parameters (you may change these if you wish to + # run a different scenario from Liu 2014) + ################################################################### + + # densities + nAr = 3.22e22 # background number density of Ar [1/m^3] (corresponds to p=100 mTorr) + np0 = 8e16 # "nominal" electron density [1/m^3] + + # masses + # me = 9.10938356e-31 # mass of an electron [kg] + # me = 5.489e-4 # mass of an electron [u] + me = 0.511e6 # mass of an electron [eV/c2] + # mAr = 39.948 # mass of an argon atom [u] + # mAr = 39.948 * 1.66054e-27 # mass of an argon atom [kg] + mAr = 37.2158e9 # mass of an electron [eV/c2] + # u = 931.4941e6 # eV/c2 + c = 299792458 # speed of light [m/s] + se = 40 # momentum cross section [m^2] + + # nominal electron energy + e0 = 1.0 # [eV] + + # pressure + p = 133.3224*1.5 # [J/m^3] *1.5 to convert it to energy (1 Torr) + + # gas energy at the wall + Tg0 = 0.038778 # 3/2*300K*kB ~ (p0 - nT[:,0])/ntot + + # characteristics of driving voltage + V0 = inputV0 # amplitude of driving voltage [V] + verticalShift = inputVDC # DC voltage (vertical shift in driving voltage) + tau = (1./13.6e6) # period of driving voltage [s] + L = 2.00*0.005 # half-gap-width [m] (gap width is 2.54cm) + electrodeArea = np.pi*0.05**2 # electrode area [m^2] (electrode diameter = 0.1 m) + + # transport parameters + nmue = 9.66e21 # argon number density times electron mobility [1/(V*cm*s)] + nmui = 4.65e19 # argon number density times ion mobility [1/(V*cm*s)] + nmum = 1 / (np.sqrt(16.0 * (mAr + mAr) * 300 * 8.62e-5 * c**2 + / (3.0 * np.pi * mAr * mAr)) * se * mAr * 1.6e-19 / c**2) + nDe = 3.86e22 # argon number density times electron diffusivity [1/(cm*s)] + nDi = 2.07e18 # argon number density times ion diffusivity [1/(cm*s)] + nDm = 2.42e18 # argon number density times metastable diffusivity [1/(cm*s)] + + # reaction parameters (NB: k_i = Ck*exp(-A/Te)) + #Ck = np.array([1.235e-7,3.712e-8,2.05e-7,1.818e-9,2e-7]) # pre-exponential factors [cm^3/s] + #A = np.array([18.687,15.06,4.95,2.14,0.0]) # activation temperature [eV] + #dH = np.array([15.7,11.56,4.14,-11.56,0.0]) # energy lost per electron due to ionization rxn [eV] + + #Ck = np.array([1.235e-7,3.712e-8,2.05e-7,1.818e-9]) # pre-exponential factors [cm^3/s] + #A = np.array([18.687,15.06,4.95,2.14]) # activation temperature [eV] + #dH = np.array([15.7,11.56,4.14,-11.56]) # energy lost per electron due to ionization rxn [eV] + + #Ck = np.array([1.235e-7,0.0,0.0,0.0,0.0]) # pre-exponential factors [cm^3/s] + #A = np.array([18.687,15.06,4.95,2.14,0.0]) # activation temperature [eV] + #dH = np.array([15.7,0.0,0.0,0.0,0.0]) # energy lost per electron due to ionization rxn [eV] + + # nominal + # Ck = np.array([1.235e-7,3.712e-8,2.05e-7,1.818e-9,2e-7,6.2e-10,3.0e-15,1.1e-31]) # pre-exponential factors [cm^3/s] + Ck = np.array([1.235e-7,3.712e-8,2.05e-7,4.0e-13,2e-7,5.0e-10,2.5e-15]) # pre-exponential factors [cm^3/s] + + # slow excitation rate + #Ck = np.array([1.235e-7,0.5*3.712e-8,2.05e-7,1.818e-9,2e-7,6.2e-10,3.0e-15,1.1e-31]) # pre-exponential factors [cm^3/s] + + ## fast excitation rate + #Ck = np.array([1.235e-7,2.0*3.712e-8,2.05e-7,1.818e-9,2e-7,6.2e-10,3.0e-15,1.1e-31]) # pre-exponential factors [cm^3/s] + + # A = np.array([18.687,15.06,4.95,2.14,0.0,0.0,0.0,0.0]) # activation temperature [eV] + # dH = np.array([15.7,11.56,4.14,-11.56,0.0,0.0,0.0,0.0]) # energy lost per electron due to ionization rxn [eV] + # dEps = np.array([0.0,15.7,11.56,0.0]) + + A = np.array([18.687,15.06,4.95,0.0,0.0,0.0,0.0]) # activation temperature [eV] + dH = np.array([15.76,11.56,4.2,0.0,-11.56,-7.56,-11.56]) # energy lost per electron due to ionization rxn [eV] + dEps = np.array([0.0,15.76,11.56,0.0]) + + # BC parameters + # ks = 1.19e7 # electron recombination rate [cm/s] + ks = 1.366109824889323e7 # electron recombination rate [cm/s/eV] + + ################################################################### + # Constants of nature (probably shouldn't change unless you have + # root privileges on universe) + ################################################################### + qe = 1.6e-19 # unit charge [C] + eps0 = 8.86e-12 # permittivity of free space [F/m] + kB = 1.38e-23 # Boltzmann constant [J/K] + # kB = 8.62e−5 # Boltzmann constant [eV/K] + + + ################################################################### + # Calculate non-dimensional parameters + ################################################################### + + # 1) Convert input units to base SI (except eV) + nDe *= 100. # 1/(m*s) + nDi *= 100. # 1/(m*s) + nDm *= 100. + nmue *= 100. # 1/(V*m*s) + nmui *= 100. # 1/(V*m*s) + Ck *= 1e-6 # m^3/s + # Ck[7] *= 1e-6 # Ck[7] is now in m^6/s + ks *= 0.01 # m/s + se *= 1.0e-20 # m^2 + + # 2) Compute "raw" transport parameters + De = nDe/nAr + Di = nDi/nAr + Dm = nDm/nAr + + mue = nmue/nAr + mui = nmui/nAr + mum = nmum/nAr + + # 3) Compute non-dimensional properties required by solver + De = De*tau/(L*L) + Di = Di*tau/(L*L) + Dm = Dm*tau/(L*L) + + mue = mue*V0*tau/(L*L) + mui = mui*V0*tau/(L*L) + mum = mum*V0*tau/(L*L) + + Ck[0:2] = Ck[0:2]*tau*nAr + #Ck[2:5] = Ck[2:5]*tau*np0 + Ck[2:6] = Ck[2:6]*tau*np0 + Ck[6] *= tau*nAr + # Ck[7] *= tau*nAr*nAr + A = A*1.5/e0 # 1.5 to convert from temperature to energy + dH = dH/e0 + qStar = V0/e0 # qe*V0/e0, since e0 in eV, need qe*V0 in eV, which is just V0 in V + alpha = qe*np0*L*L/(V0*eps0) + ks = ks*tau/L + p0 = p/qe/np0 + kappaB = 4.878171165833662*1.6129 # non-dimensional thermal conductivity of background specie + # (2/3)*tau/L**2*Kb/np0/kB, + # where Kb is the thermal conductivity of background specie + + #params.beta = np.array([[2,2,2,1,1],[1,0,1,0,0],[0,1,0,0,0],[0,0,0,1,1]]) + #params.alfa = np.array([[1,1,1,1,1],[0,0,0,0,0],[0,0,1,1,1],[1,1,0,0,0]]) + #params.beta = np.array([[2,1,2,1],[1,0,1,0],[0,1,0,0],[0,0,0,1]], dtype=np.int) + #params.alfa = np.array([[1,1,1,1],[0,0,0,0],[0,0,1,1],[1,1,0,0]], dtype=np.int) + # params.beta = np.array([[2,1,2,1,1,1,0,0],[1,0,1,0,0,1,0,0],[0,1,0,0,0,0,0,0],[0,0,0,1,0,1,2,1]], dtype=np.int64) + # params.alfa = np.array([[1,1,1,1,1,0,0,0],[0,0,0,0,0,0,0,0],[0,0,1,1,1,2,1,1],[1,1,0,0,0,0,1,2]], dtype=np.int64) + params.beta = np.array([[2,1,2,0,1,1,0],[1,0,1,0,0,1,0],[0,1,0,1,0,0,0],[0,0,0,0,1,1,2]], dtype=np.int64) + params.alfa = np.array([[1,1,1,1,1,0,0],[0,0,0,1,0,0,0],[0,0,1,0,1,2,1],[1,1,0,0,0,0,1]], dtype=np.int64) + + # 4) Set values in params class + params.D[0] = De + params.D[1] = Di + params.D[2] = Dm + + params.mu[0] = mue + params.mu[1] = mui + params.mu[2] = mum + + params.A[:] = Ck[:] + params.B[:] = 0.0 + params.C[:] = A[:] + + params.dH[:] = dH[:] + params.dEps[:] = dEps[:] + params.qStar = qStar + params.alpha = alpha + params.ks = ks + params.gam = gam + params.kappaB = kappaB + params.nAronp0 = nAr / np0 + params.p0 = p0 + params.Tg0 = Tg0 + params.EC = 2.0 * me / mAr \ + * np.sqrt(16.0 * (me + mAr) * e0 * c**2 + / (3.0 * np.pi * me * mAr)) * se * nAr * tau + # params.EC = 2.0 * me / mAr * 3.8e9 * tau + + params.verticalShift = verticalShift / V0 + + # Parameters needed to compute the current with dimensions + params.V0Ltau = V0 / (L * tau) + params.V0L = V0 / L + params.LLV0tau = (L*L) / (V0*tau) + params.tauL = L / tau + params.np0 = np0 # "nominal" electron density [1/m^3] + params.qe = qe # unit charge [C] + params.eps0 = eps0 # unit charge [C] + params.eArea = electrodeArea # electrode area [m^2] + + reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + f"{params.A[3]} * energy**({params.B[3]}-0.5)", + f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)"] + + reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + f"{params.A[3]} * (energy**({params.B[3]}-0.5-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]}-0.5 + {params.C[3]}/energy)", + f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)"] + + reactionExpressionTypelist = np.array([True,True,True,False,False,False,False]) + + reactionsList = [] + LOGFilename = 'interpolationSample%s.log'%str(iSample) + f = open(LOGFilename, 'w') + for i in range(Nr): + if reactionExpressionTypelist[i]: + Nsample = 7200 + N300 = 200 + + #rateCoeff = np.fromfile("/g/g92/jbarbere/glowDischarge/toyProblems/timeDomain/BOLSIGChemistryZeroIonDeg7200Samples/reaction300K_%s.dat" %str(i)) + #rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + + #Te = np.fromfile("/g/92/jbarbere/glowDischarge/toyProblems/timeDomain/BOLSIGChemistryZeroIonDeg7200Samples/reaction300K.Te.dat") + #Te = np.reshape(Te,[Nsample, N300]).T[:,iSample] + + root_dir = "/g/g92/jbarbere/glowDischarge/toyProblems/timeDomain" + rate_file = "{0:s}/BOLSIGChemistryZeroIonDeg7200Samples/reaction300K_{1:d}.dat".format(root_dir, i) + temp_file = "{0:s}/BOLSIGChemistryZeroIonDeg7200Samples/reaction300K.Te.dat".format(root_dir) + + rateCoeff = np.fromfile(rate_file) + rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + + Te = np.fromfile(temp_file) + Te = np.reshape(Te,[Nsample, N300]).T[:,iSample] + + # 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] + 1 + + # Transformation to log scale. + TeLog = np.log(Te) + + # Find first non-zero value of the coefficient rate. + # I = np.nonzero(rateCoeff) + + # 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] + C = (np.log(rateCoeff[lastFalse + 1]) - np.log(rateCoeff[lastFalse])) \ + / (1.0 / Te[lastFalse] - 1.0 / Te[lastFalse + 1]) + + # 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) + + logging.basicConfig(filename=LOGFilename) + logging.warning('Interpolation info (Reaction %s):', i) + logging.warning('First non-zero entry in the rate coefficient: %s', I[0][0]) + logging.warning('Monotonicity of the rate coefficient start from entry: %s', lastFalse) + logging.warning('Position of possible duplicates in mean energy array: %s', + TeDuplicateindsForLog[0]) + + # 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') + + # TeFinerResolution = np.linspace(5.0e-4, 100, 20000) + # TeFinerResolutionLog = np.log(TeFinerResolution) + + # # plot the function + # # plt.plot(rateCoeffXFiner, np.exp(reactionExpressions_cubicSplineDerivative_log(rateCoeffXFiner)), + # # color='salmon', linestyle='--', label='interBolsig') + # plt.plot(TeFinerResolution, np.exp(reactionExpressionsLog(TeFinerResolutionLog)), + # 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("./Rates_%s.pdf" %str(i), dpi=300) + # # plt.xlim((-0.0001,0.0255)) + # plt.show() + else: + rxn = eval("lambda energy :" + reactionExpressionslist[i]) + rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) + + reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]], + rxnBolsig = reactionExpressionTypelist[i], + kf = rxn, kf_T = rxn_T) + reactionsList.append(reaction) + + params.reactionsList = reactionsList + #params.Nr = 1 + + # 5) Dump to screen + params.print() diff --git a/psaapProperties_4plus2Species.py b/psaapProperties_4plus2Species.py new file mode 100755 index 000000000..e9f41a1b5 --- /dev/null +++ b/psaapProperties_4plus2Species.py @@ -0,0 +1,374 @@ +import numpy as np +from scipy.interpolate import CubicSpline +import csv +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors + +import logging + +class Reaction(object): + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) + +def setPsaapProperties_4plus2Species(gam, inputV0, inputVDC, params, Nr, iSample): + """Sets non-dimensional properties corresponding to Liu 2014 paper. + + Inputs: + gam : Secondary electron emission coefficient + params : chebSolver.modelParams class + + Outputs: None + params data is overwritten using values from Liu 2014. + """ + ################################################################### + # User specified parameters (you may change these if you wish to + # run a different scenario from Liu 2014) + ################################################################### + + # densities + nAr = 3.22e22 # background number density of Ar [1/m^3] (corresponds to p=100 mTorr) + np0 = 8e16 # "nominal" electron density [1/m^3] + + # masses + # me = 9.10938356e-31 # mass of an electron [kg] + # me = 5.489e-4 # mass of an electron [u] + me = 0.511e6 # mass of an electron [eV/c2] + # mAr = 39.948 # mass of an argon atom [u] + # mAr = 39.948 * 1.66054e-27 # mass of an argon atom [kg] + mAr = 37.2158e9 # mass of an electron [eV/c2] + # u = 931.4941e6 # eV/c2 + c = 299792458 # speed of light [m/s] + se = 40 # momentum cross section [A^2] + + # nominal electron energy + e0 = 1.0 # [eV] + + # pressure + p = 133.3224*1.5 # [J/m^3] *1.5 to convert it to energy (1 Torr) + + # gas energy at the wall + Tg0 = 0.038778 # 3/2*300K*kB ~ (p0 - nT[:,0])/ntot + + # characteristics of driving voltage + V0 = inputV0 # amplitude of driving voltage [V] + verticalShift = inputVDC # DC voltage (vertical shift in driving voltage) + tau = (1./13.56e6) # period of driving voltage [s] + L = 2.00*0.005 # half-gap-width [m] (gap width is 2 cm) + electrodeArea = np.pi*0.05**2 # electrode area [m^2] (electrode diameter = 0.1 m) + + # transport parameters + nmue = 9.66e21 # argon number density times electron mobility [1/(V*cm*s)] + #nmui = 4.65e19 # argon number density times ion mobility [1/(V*cm*s)] + #nmum = 1 / (np.sqrt(16.0 * (mAr + mAr) * 300 * 8.62e-5 * c**2 + # / (3.0 * np.pi * mAr * mAr)) * se * mAr * 1.6e-19 / c**2) + nmum = 0.0 + nmur = 0.0 + nmu4p = 0.0 + nmui = 8.0e19 + nDe = 3.86e22 # argon number density times electron diffusivity [1/(cm*s)] + nDi = 2.07e18 # argon number density times ion diffusivity [1/(cm*s)] + nDm = 2.42e18 # argon number density times metastable diffusivity [1/(cm*s)] + nDr = 2.42e18 + nD4p = 2.42e18 + #nDm = 3.914e20 + + # reaction parameters (NB: k_i = Ck*Ee^B*exp(-A/Ee)) + # Ee = 3/2*Te (Te in eV) + # -> k_i = [Ck*(2/3)^B] * Ee^B * exp[-(3/2)*A/Ee] + # nominal + Ck = np.array([1.235e-7,3.712e-8,2.05e-7,4.0e-13,5.0e-10,4.3e-10,2.1e-15,10.0,5.0e-27]) # pre-exponential factors [cm^3/s] + B = np.array([0.0,0.0,0.0,-0.5,0.0,0.74,0.0,0.0,-4.5]) + A = np.array([18.687,15.06,4.95,0.0,0.0,0.0,0.0,0.0,0.0]) # activation temperature [eV] + dH = np.array([15.76,11.56,4.2,0.0,-7.56,-11.56,-11.56,0.0,-4.2]) # energy lost per electron due to ionization rxn [eV] + dEps = np.array([0.0,15.76,11.56,0.0,0.0,0.0]) + + # BC parameters + # ks = 1.19e7 # electron recombination rate [cm/s] + ks = 1.366109824889323e7 # electron recombination rate [cm/s/eV] + + ################################################################### + # Constants of nature (probably shouldn't change unless you have + # root privileges on universe) + ################################################################### + qe = 1.6e-19 # unit charge [C] + eps0 = 8.86e-12 # permittivity of free space [F/m] + kB = 1.38e-23 # Boltzmann constant [J/K] + # kB = 8.62e-5 # Boltzmann constant [eV/K] + + + ################################################################### + # Calculate non-dimensional parameters + ################################################################### + + # 1) Convert input units to base SI (except eV) + nDe *= 100. # 1/(m*s) + nDi *= 100. # 1/(m*s) + nDm *= 100. + nDr *= 100. + nD4p *= 100. + nmue *= 100. # 1/(V*m*s) + nmui *= 100. # 1/(V*m*s) + Ck[0:7] *= 1e-6 # m^3/s + Ck[7] *= 1 # 1/s + Ck[8] *= 1e-12 # m^6/s + ks *= 0.01 # m/s + se *= 1.0e-20 # m^2 + + # 2) Compute "raw" transport parameters + De = nDe/nAr + Di = nDi/nAr + Dm = nDm/nAr + Dr = nDr/nAr + D4p = nD4p/nAr + + mue = nmue/nAr + mui = nmui/nAr + mum = nmum/nAr + mur = nmur/nAr + mu4p = nmu4p/nAr + + # 3) Compute non-dimensional properties required by solver + De = De*tau/(L*L) + Di = Di*tau/(L*L) + Dm = Dm*tau/(L*L) + Dr = Dr*tau/(L*L) + D4p = D4p*tau/(L*L) + + mue = mue*V0*tau/(L*L) + mui = mui*V0*tau/(L*L) + mum = mum*V0*tau/(L*L) + mur = mur*V0*tau/(L*L) + mu4p = mu4p*V0*tau/(L*L) + + Ck[0:2] = Ck[0:2]*tau*nAr + Ck[2:6] = Ck[2:6]*tau*np0 + Ck[6] *= tau*nAr + Ck[7] *= tau + Ck[8] *= tau*np0*np0 + A = A*1.5/e0 # 1.5 to convert from temperature to energy + dH = dH/e0 + qStar = V0/e0 # qe*V0/e0, since e0 in eV, need qe*V0 in eV, which is just V0 in V + alpha = qe*np0*L*L/(V0*eps0) + ks = ks*tau/L + p0 = p/qe/np0 + kappaB = 4.878171165833662*1.6129 # non-dimensional thermal conductivity of background specie + # (2/3)*tau/L**2*Kb/np0/kB, + # where Kb is the thermal conductivity of background specie + + params.beta = np.array([[2,1,2,0,1,1,0,0,1], + [1,0,1,0,1,0,0,0,0], + [0,1,0,1,0,0,0,0,1], + [0,0,0,0,0,0,0,0,0], + [0,0,0,0,0,0,0,0,0], + [0,0,0,0,1,1,2,1,0]], dtype=np.int64) + params.alfa = np.array([[1,1,1,1,0,1,0,0,2], + [0,0,0,1,0,0,0,0,1], + [0,0,1,0,2,1,1,1,0], + [0,0,0,0,0,0,0,0,0], + [0,0,0,0,0,0,0,0,0], + [1,1,0,0,0,0,1,0,0]], dtype=np.int64) + # Rxn1: E + AR -> 2E + AR+ + # Rxn2: E + AR -> E + AR* + # Rxn3: E + AR* -> 2E + AR+ + # Rxn4: E + AR+ -> AR* + # Rxn5: 2AR* -> E + AR+ + AR + # Rxn6: E + AR* -> E + AR + # Rxn7: AR* + AR -> 2AR + # Rxn8: AR* -> AR + # Rxn9: 2E + AR+ -> E + AR* + + # 4) Set values in params class + params.D[0] = De + params.D[1] = Di + params.D[2] = Dm + params.D[3] = Dr + params.D[4] = D4p + + params.mu[0] = mue + params.mu[1] = mui + params.mu[2] = mum + params.mu[3] = mur + params.mu[4] = mu4p + + params.A[:] = Ck[:] + params.B[:] = B[:] + params.C[:] = A[:] + + # Account for the 2/3 term to convert from electron temperature to electron energy + for i in range(len(params.A)): + params.A[i] *= (2/3)**(params.B[i]) + + params.dH[:] = dH[:] + params.dEps[:] = dEps[:] + params.qStar = qStar + params.alpha = alpha + params.ks = ks + params.gam = gam + params.kappaB = kappaB + params.nAronp0 = nAr / np0 + params.p0 = p0 + params.Tg0 = Tg0 + params.EC = 2.0 * me / mAr \ + * np.sqrt(16.0 * (me + mAr) * e0 * c**2 + / (3.0 * np.pi * me * mAr)) * se * nAr * tau + # params.EC = 2.0 * me / mAr * 3.8e9 * tau + + params.verticalShift = verticalShift / V0 + + # Parameters needed to compute the current with dimensions + params.V0Ltau = V0 / (L * tau) + params.V0L = V0 / L + params.LLV0tau = (L*L) / (V0*tau) + params.tauL = L / tau + params.np0 = np0 # "nominal" electron density [1/m^3] + params.qe = qe # unit charge [C] + params.eps0 = eps0 # unit charge [C] + params.eArea = electrodeArea # electrode area [m^2] + + + reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + f"{params.A[3]} * energy**{params.B[3]} * np.exp(-{params.C[3]} / energy)", + f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)", + f"{params.A[7]} * energy**{params.B[7]} * np.exp(-{params.C[7]} / energy)", + f"{params.A[8]} * energy**{params.B[8]} * np.exp(-{params.C[8]} / energy)"] + + reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + f"{params.A[3]} * (energy**({params.B[3]}-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]} + {params.C[3]}/energy)", + f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)", + f"{params.A[7]} * (energy**({params.B[7]}-1)) * np.exp(-{params.C[7]}/energy) * ({params.B[7]} + {params.C[7]}/energy)", + f"{params.A[8]} * (energy**({params.B[8]}-1)) * np.exp(-{params.C[8]}/energy) * ({params.B[8]} + {params.C[8]}/energy)"] + + reactionExpressionTypelist = np.array([True,True,True,False,False,False,False,False,False]) + + reactionsList = [] + LOGFilename = 'interpolationSample%s.log'%str(iSample) + f = open(LOGFilename, 'w') + # 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}/BOLSIGChemistry_NominalRates/reaction300K_{1:d}.txt".format(root_dir, i), 'r') + temp_file = open("{0:s}/BOLSIGChemistry_NominalRates/reaction300K_Te.txt".format(root_dir), 'r') + + #rateCoeff = np.fromfile(rate_file) + rateCoeff = np.genfromtxt(rate_file) + rate_file.close() + rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + + #Te = np.fromfile(temp_file) + Te = np.genfromtxt(temp_file) + temp_file.close() + Te = np.reshape(Te,[Nsample, N300]).T[:,iSample] + + # 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) + + logging.basicConfig(filename=LOGFilename) + logging.warning('Interpolation info (Reaction %s):', i) + logging.warning('First non-zero entry in the rate coefficient: %s', I[0][0]) + logging.warning('Monotonicity of the rate coefficient start from entry: %s', lastFalse) + logging.warning('Position of possible duplicates in mean energy array: %s', + TeDuplicateindsForLog[0]) + + else: + rxn = eval("lambda energy :" + reactionExpressionslist[i]) + rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) + + reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]], + rxnBolsig = reactionExpressionTypelist[i], + kf = rxn, kf_T = rxn_T) + reactionsList.append(reaction) + + params.reactionsList = reactionsList + + # 5) Dump to screen + params.print() diff --git a/psaapProperties_6Species.py b/psaapProperties_6Species.py new file mode 100755 index 000000000..e02c158af --- /dev/null +++ b/psaapProperties_6Species.py @@ -0,0 +1,430 @@ +import numpy as np +from scipy.interpolate import CubicSpline +import csv +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors + +import logging + +class Reaction(object): + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) + +def setPsaapProperties_6Species(gam, inputV0, inputVDC, params, Nr, iSample): + """Sets non-dimensional properties corresponding to Liu 2014 paper. + + Inputs: + gam : Secondary electron emission coefficient + params : chebSolver.modelParams class + + Outputs: None + params data is overwritten using values from Liu 2014. + """ + ################################################################### + # User specified parameters (you may change these if you wish to + # run a different scenario from Liu 2014) + ################################################################### + + # densities + nAr = 3.22e22 # background number density of Ar [1/m^3] (corresponds to p=100 mTorr) + np0 = 8e16 # "nominal" electron density [1/m^3] + + # masses + # me = 9.10938356e-31 # mass of an electron [kg] + # me = 5.489e-4 # mass of an electron [u] + me = 0.511e6 # mass of an electron [eV/c2] + # mAr = 39.948 # mass of an argon atom [u] + # mAr = 39.948 * 1.66054e-27 # mass of an argon atom [kg] + mAr = 37.2158e9 # mass of an argon atom [eV/c2] + # u = 931.4941e6 # eV/c2 + c = 299792458 # speed of light [m/s] + se = 40 # momentum cross section [A^2] + + # nominal electron energy + e0 = 1.0 # [eV] + + # pressure + p = 133.3224*1.5 # [J/m^3] *1.5 to convert it to energy (1 Torr) + + # gas energy at the wall + Tg0 = 0.038778 # 3/2*300K*kB ~ (p0 - nT[:,0])/ntot + + # characteristics of driving voltage + V0 = inputV0 # amplitude of driving voltage [V] + verticalShift = inputVDC # DC voltage (vertical shift in driving voltage) + tau = (1./13.56e6) # period of driving voltage [s] + L = 2.00*0.005 # half-gap-width [m] (gap width is 2 cm) + electrodeArea = np.pi*0.05**2 # electrode area [m^2] (electrode diameter = 0.1 m) + + # transport parameters + nmue = 9.66e21 # argon number density times electron mobility [1/(V*cm*s)] + #nmui = 4.65e19 # argon number density times ion mobility [1/(V*cm*s)] + #nmum = 1 / (np.sqrt(16.0 * (mAr + mAr) * 300 * 8.62e-5 * c**2 + # / (3.0 * np.pi * mAr * mAr)) * se * mAr * 1.6e-19 / c**2) + nmum = 0.0 + nmur = 0.0 + nmu4p = 0.0 + nmui = 8.0e19 + #nmum = 9.35e19 + nDe = 3.86e22 # argon number density times electron diffusivity [1/(cm*s)] + nDi = 2.07e18 # argon number density times ion diffusivity [1/(cm*s)] + nDm = 2.42e18 # argon number density times AR(m) diffusivity [1/(cm*s)] + nDr = 2.42e18 + nD4p = 2.42e18 + #nDm = 3.914e20 + + # reaction parameters (NB: k_i = Ck*Ee^B*exp(-A/Ee)) + # Ee = 3/2*Te (Te in eV) + # -> k_i = [Ck*(2/3)^B] * Ee^B * exp[-(3/2)*A/Ee] + # nominal + Ck = np.array([2.0e-7,2.1e-9,5.0e-10,6.4e-10,2.1e-15,1.0e5,3.2e7,3.0e7,3.0e7,0.0,0.0,0.0,0.0,4.3e-10,0.0,3.7e-8,8.9e-7,1.8e-7,3.0e-7,3.0e-7,4.3e-10,9.1e-7,8.9e-7]) # pre-exponential factors [cm^3/s] + B = np.array([0,0,0,0,0,0,0,0,0,0,0,0,0,0.74,0,0,0.51,0.61,0.51,0.51,0.74,0,0.51]) + A = np.array([0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.59,2.61,0,0,0,0,1.59]) # activation temperature [eV] + dH = np.array([0.0,-7.412,-10.054,-7.336,0.0,0.0,0.0,0.0,0.0,11.548,11.624,12.907,15.76,-11.548,4.212,0.076,1.359,2.853,-1.283,-1.359,-11.624,-0.076,0.983]) # energy lost per electron due to ionization rxn [eV] + dEps = np.array([0.0,15.76,11.548,11.624,12.907,0.0]) # E, AR+, AR(m), AR(r), AR(4p), AR + + # BC parameters + # ks = 1.19e7 # electron recombination rate [cm/s] + ks = 1.366109824889323e7 # electron recombination rate [cm/s/eV] + + ################################################################### + # Constants of nature (probably shouldn't change unless you have + # root privileges on universe) + ################################################################### + qe = 1.6e-19 # unit charge [C] + eps0 = 8.86e-12 # permittivity of free space [F/m] + kB = 1.38e-23 # Boltzmann constant [J/K] + # kB = 8.62e-5 # Boltzmann constant [eV/K] + + + ################################################################### + # Calculate non-dimensional parameters + ################################################################### + + # 1) Convert input units to base SI (except eV) + nDe *= 100. # 1/(m*s) + nDi *= 100. # 1/(m*s) + nDm *= 100. + nDr *= 100. + nD4p *= 100. + nmue *= 100. # 1/(V*m*s) + nmui *= 100. # 1/(V*m*s) + Ck[0:5] *= 1e-6 # m^3/s + Ck[5:9] *= 1 # 1/s + Ck[9:22] *= 1e-6 # m^3/s + ks *= 0.01 # m/s + se *= 1.0e-20 # m^2 + + # 2) Compute "raw" transport parameters + De = nDe/nAr + Di = nDi/nAr + Dm = nDm/nAr + Dr = nDr/nAr + D4p = nD4p/nAr + + mue = nmue/nAr + mui = nmui/nAr + mum = nmum/nAr + mur = nmur/nAr + mu4p = nmu4p/nAr + + # 3) Compute non-dimensional properties required by solver + De = De*tau/(L*L) + Di = Di*tau/(L*L) + Dm = Dm*tau/(L*L) + Dr = Dr*tau/(L*L) + D4p = D4p*tau/(L*L) + + mue = mue*V0*tau/(L*L) + mui = mui*V0*tau/(L*L) + mum = mum*V0*tau/(L*L) + mur = mur*V0*tau/(L*L) + mu4p = mu4p*V0*tau/(L*L) + + Ck[0:4] *= tau*np0 + Ck[4] *= tau*nAr + Ck[5:9] *= tau + Ck[9:13] *= tau*nAr + Ck[13:] *= tau*np0 + A = A*1.5/e0 # 1.5 to convert from temperature to energy + dH = dH/e0 + qStar = V0/e0 # qe*V0/e0, since e0 in eV, need qe*V0 in eV, which is just V0 in V + alpha = qe*np0*L*L/(V0*eps0) + ks = ks*tau/L + p0 = p/qe/np0 + kappaB = 4.878171165833662*1.6129 # non-dimensional thermal conductivity of background specie + # (2/3)*tau/L**2*Kb/np0/kB, + # where Kb is the thermal conductivity of background specie + + params.beta = np.array([[0,1,1,1,0,0,0,0,0,1,1,1,2,1,2,1,1,2,1,1,1,1,1], # E + [0,1,1,1,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,0,0], # AR+ + [0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,1,0,1,0], # AR(m) + [0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,1,0,0,1,0,0,0,0], # AR(r) + [0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1], # AR(4p) + [2,1,1,1,2,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0]], dtype=np.int64) # AR + + params.alfa = np.array([[0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1], # E + [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], # AR+ + [2,1,0,2,1,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0], # AR(m) + [0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1], # AR(r) + [0,0,2,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,1,1,0,0,0], # AR(4p) + [0,0,0,0,1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0]], dtype=np.int64) # AR + # Rxn1: 2AR(m) -> 2AR + # Rxn2: AR(m) + AR(r) -> E + AR+ + AR + # Rxn3: 2AR(4p) -> E + AR+ + AR + # Rxn4: 2AR(m) -> E + AR+ + AR + # Rxn5: AR(m) + AR -> 2AR + # Rxn6: AR(r) -> AR + # Rxn7: AR(4p) -> AR + # Rxn8: AR(4p) -> AR(m) + # Rxn9: AR(4p) -> AR(r) + # Rxn10: E + AR -> E + AR(m) + # Rxn11: E + AR -> E + AR(r) + # Rxn12: E + AR -> E + AR(4p) + # Rxn13: E + AR -> 2E + AR+ + # Rxn14: E + AR(m) -> E + AR + # Rxn15: E + AR(m) -> 2E + AR+ + # Rxn16: E + AR(m) -> E + AR(r) + # Rxn17: E + AR(m) -> E + AR(4p) + # Rxn18: E + AR(4p) -> 2E + AR+ + # Rxn19: E + AR(4p) -> E + AR(r) + # Rxn20: E + AR(4p) -> E + AR(m) + # Rxn21: E + AR(r) -> E + AR + # Rxn22: E + AR(r) -> E + AR(m) + # Rxn23: E + AR(r) -> E + AR(4p) + + # 4) Set values in params class + params.D[0] = De + params.D[1] = Di + params.D[2] = Dm + params.D[3] = Dr + params.D[4] = D4p + + params.mu[0] = mue + params.mu[1] = mui + params.mu[2] = mum + params.mu[3] = mur + params.mu[4] = mu4p + + params.A[:] = Ck[:] + params.B[:] = B[:] + params.C[:] = A[:] + + # Account for the 2/3 term to convert from electron temperature to electron energy + for i in range(len(params.A)): + params.A[i] *= (2/3)**(params.B[i]) + + params.dH[:] = dH[:] + params.dEps[:] = dEps[:] + params.qStar = qStar + params.alpha = alpha + params.ks = ks + params.gam = gam + params.kappaB = kappaB + params.nAronp0 = nAr / np0 + params.p0 = p0 + params.Tg0 = Tg0 + params.EC = 2.0 * me / mAr \ + * np.sqrt(16.0 * (me + mAr) * e0 * c**2 + / (3.0 * np.pi * me * mAr)) * se * nAr * tau + # params.EC = 2.0 * me / mAr * 3.8e9 * tau + + params.verticalShift = verticalShift / V0 + + # Parameters needed to compute the current with dimensions + params.V0Ltau = V0 / (L * tau) + params.V0L = V0 / L + params.LLV0tau = (L*L) / (V0*tau) + params.tauL = L / tau + params.np0 = np0 # "nominal" electron density [1/m^3] + params.qe = qe # unit charge [C] + params.eps0 = eps0 # unit charge [C] + params.eArea = electrodeArea # electrode area [m^2] + + + reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + f"{params.A[3]} * energy**{params.B[3]} * np.exp(-{params.C[3]} / energy)", + f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)", + f"{params.A[7]} * energy**{params.B[7]} * np.exp(-{params.C[7]} / energy)", + f"{params.A[8]} * energy**{params.B[8]} * np.exp(-{params.C[8]} / energy)", + f"{params.A[9]} * energy**{params.B[9]} * np.exp(-{params.C[9]} / energy)", + f"{params.A[10]} * energy**{params.B[10]} * np.exp(-{params.C[10]} / energy)", + f"{params.A[11]} * energy**{params.B[11]} * np.exp(-{params.C[11]} / energy)", + f"{params.A[12]} * energy**{params.B[12]} * np.exp(-{params.C[12]} / energy)", + f"{params.A[13]} * energy**{params.B[13]} * np.exp(-{params.C[13]} / energy)", + f"{params.A[14]} * energy**{params.B[14]} * np.exp(-{params.C[14]} / energy)", + f"{params.A[15]} * energy**{params.B[15]} * np.exp(-{params.C[15]} / energy)", + f"{params.A[16]} * energy**{params.B[16]} * np.exp(-{params.C[16]} / energy)", + f"{params.A[17]} * energy**{params.B[17]} * np.exp(-{params.C[17]} / energy)", + f"{params.A[18]} * energy**{params.B[18]} * np.exp(-{params.C[18]} / energy)", + f"{params.A[19]} * energy**{params.B[19]} * np.exp(-{params.C[19]} / energy)", + f"{params.A[20]} * energy**{params.B[20]} * np.exp(-{params.C[20]} / energy)", + f"{params.A[21]} * energy**{params.B[21]} * np.exp(-{params.C[21]} / energy)", + f"{params.A[22]} * energy**{params.B[22]} * np.exp(-{params.C[22]} / energy)"] + + reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + f"{params.A[3]} * (energy**({params.B[3]}-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]} + {params.C[3]}/energy)", + f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)", + f"{params.A[7]} * (energy**({params.B[7]}-1)) * np.exp(-{params.C[7]}/energy) * ({params.B[7]} + {params.C[7]}/energy)", + f"{params.A[8]} * (energy**({params.B[8]}-1)) * np.exp(-{params.C[8]}/energy) * ({params.B[8]} + {params.C[8]}/energy)", + f"{params.A[9]} * (energy**({params.B[9]}-1)) * np.exp(-{params.C[9]}/energy) * ({params.B[9]} + {params.C[9]}/energy)", + f"{params.A[10]} * (energy**({params.B[10]}-1)) * np.exp(-{params.C[10]}/energy) * ({params.B[10]} + {params.C[10]}/energy)", + f"{params.A[11]} * (energy**({params.B[11]}-1)) * np.exp(-{params.C[11]}/energy) * ({params.B[11]} + {params.C[11]}/energy)", + f"{params.A[12]} * (energy**({params.B[12]}-1)) * np.exp(-{params.C[12]}/energy) * ({params.B[12]} + {params.C[12]}/energy)", + f"{params.A[13]} * (energy**({params.B[13]}-1)) * np.exp(-{params.C[13]}/energy) * ({params.B[13]} + {params.C[13]}/energy)", + f"{params.A[14]} * (energy**({params.B[14]}-1)) * np.exp(-{params.C[14]}/energy) * ({params.B[14]} + {params.C[14]}/energy)", + f"{params.A[15]} * (energy**({params.B[15]}-1)) * np.exp(-{params.C[15]}/energy) * ({params.B[15]} + {params.C[15]}/energy)", + f"{params.A[16]} * (energy**({params.B[16]}-1)) * np.exp(-{params.C[16]}/energy) * ({params.B[16]} + {params.C[16]}/energy)", + f"{params.A[17]} * (energy**({params.B[17]}-1)) * np.exp(-{params.C[17]}/energy) * ({params.B[17]} + {params.C[17]}/energy)", + f"{params.A[18]} * (energy**({params.B[18]}-1)) * np.exp(-{params.C[18]}/energy) * ({params.B[18]} + {params.C[18]}/energy)", + f"{params.A[19]} * (energy**({params.B[19]}-1)) * np.exp(-{params.C[19]}/energy) * ({params.B[19]} + {params.C[19]}/energy)", + f"{params.A[20]} * (energy**({params.B[20]}-1)) * np.exp(-{params.C[20]}/energy) * ({params.B[20]} + {params.C[20]}/energy)", + f"{params.A[21]} * (energy**({params.B[21]}-1)) * np.exp(-{params.C[21]}/energy) * ({params.B[21]} + {params.C[21]}/energy)", + f"{params.A[22]} * (energy**({params.B[22]}-1)) * np.exp(-{params.C[22]}/energy) * ({params.B[22]} + {params.C[22]}/energy)"] + + + reactionExpressionTypelist = np.array([False,False,False,False,False,False,False,False,False, + True,True,True,True,False,True,False,False,False,False,False,False,False,False]) + + reactionsList = [] + LOGFilename = 'interpolationSample%s.log'%str(iSample) + f = open(LOGFilename, 'w') + # 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}/BOLSIGChemistry_6SpeciesRates/reaction300K_{1:d}.txt".format(root_dir, i), 'r') + temp_file = open("{0:s}/BOLSIGChemistry_6SpeciesRates/reaction300K_Te.txt".format(root_dir), 'r') + + #rateCoeff = np.fromfile(rate_file) + rateCoeff = np.genfromtxt(rate_file) + rate_file.close() + rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + + #Te = np.fromfile(temp_file) + Te = np.genfromtxt(temp_file) + temp_file.close() + Te = np.reshape(Te,[Nsample, N300]).T[:,iSample] + + # 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 + for k in range(len(Te)): + if Te[k] < 4.5 and indices[k] == False: + lastFalse = k + 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 == 9: + rateCoeffLog += - np.log(1.0/tau) + np.log(nAr) + elif i == 10: + rateCoeffLog += - np.log(1.0/tau) + np.log(nAr) + elif i == 11: + rateCoeffLog += - np.log(1.0/tau) + np.log(nAr) + elif i == 12: + 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) + + logging.basicConfig(filename=LOGFilename) + logging.warning('Interpolation info (Reaction %s):', i + 1) + logging.warning('First non-zero entry in the rate coefficient: %s', I[0][0]) + logging.warning('Monotonicity of the rate coefficient start from entry: %s', lastFalse) + logging.warning('Position of possible duplicates in mean energy array: %s', + TeDuplicateindsForLog[0]) + + else: + rxn = eval("lambda energy :" + reactionExpressionslist[i]) + rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) + + reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]], + rxnBolsig = reactionExpressionTypelist[i], + kf = rxn, kf_T = rxn_T) + reactionsList.append(reaction) + + params.reactionsList = reactionsList + + # 5) Dump to screen + params.print() diff --git a/psaapProperties_6Species_Nominal.py b/psaapProperties_6Species_Nominal.py new file mode 100755 index 000000000..9c6cc8338 --- /dev/null +++ b/psaapProperties_6Species_Nominal.py @@ -0,0 +1,461 @@ +import numpy as np +from scipy.interpolate import CubicSpline +import csv +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors +import h5py as h5 + +import logging + +class Reaction(object): + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) + +def setPsaapProperties_6Species_Nominal(gam, inputV0, inputVDC, params, Nr, iSample): + """Sets non-dimensional properties corresponding to Liu 2014 paper. + + Inputs: + gam : Secondary electron emission coefficient + params : chebSolver.modelParams class + + Outputs: None + params data is overwritten using values from Liu 2014. + """ + ################################################################### + # User specified parameters (you may change these if you wish to + # run a different scenario from Liu 2014) + ################################################################### + + # densities + nAr = 3.22e22 # background number density of Ar [1/m^3] (corresponds to p=100 mTorr) + np0 = 8e16 # "nominal" electron density [1/m^3] + + # masses + # me = 9.10938356e-31 # mass of an electron [kg] + # me = 5.489e-4 # mass of an electron [u] + me = 0.511e6 # mass of an electron [eV/c2] + # mAr = 39.948 # mass of an argon atom [u] + # mAr = 39.948 * 1.66054e-27 # mass of an argon atom [kg] + mAr = 37.2158e9 # mass of an argon atom [eV/c2] + # u = 931.4941e6 # eV/c2 + c = 299792458 # speed of light [m/s] + se = 40 # momentum cross section [A^2] + + # nominal electron energy + e0 = 1.0 # [eV] + + # pressure + p = 133.3224*1.5 # [J/m^3] *1.5 to convert it to energy (1 Torr) + + # gas energy at the wall + Tg0 = 0.038778 # 3/2*300K*kB ~ (p0 - nT[:,0])/ntot + + # characteristics of driving voltage + V0 = inputV0 # amplitude of driving voltage [V] + verticalShift = inputVDC # DC voltage (vertical shift in driving voltage) + tau = (1./13.56e6) # period of driving voltage [s] + L = 2.00*0.005 # half-gap-width [m] (gap width is 2 cm) + electrodeArea = np.pi*0.05**2 # electrode area [m^2] (electrode diameter = 0.1 m) + + # transport parameters + nmue = 9.66e21 # argon number density times electron mobility [1/(V*cm*s)] + #nmui = 4.65e19 # argon number density times ion mobility [1/(V*cm*s)] + #nmum = 1 / (np.sqrt(16.0 * (mAr + mAr) * 300 * 8.62e-5 * c**2 + # / (3.0 * np.pi * mAr * mAr)) * se * mAr * 1.6e-19 / c**2) + nmum = 0.0 + nmur = 0.0 + nmu4p = 0.0 + nmui = 8.0e19 + #nmum = 9.35e19 + nDe = 3.86e22 # argon number density times electron diffusivity [1/(cm*s)] + nDi = 2.07e18 # argon number density times ion diffusivity [1/(cm*s)] + nDm = 2.42e18 # argon number density times AR(m) diffusivity [1/(cm*s)] + nDr = 2.42e18 + nD4p = 2.42e18 + #nDm = 3.914e20 + + # reaction parameters (NB: k_i = Ck*Ee^B*exp(-A/Ee)) + # Ee = 3/2*Te (Te in eV) + # -> k_i = [Ck*(2/3)^B] * Ee^B * exp[-(3/2)*A/Ee] + # nominal + Ck = np.array([2.0e-7,2.1e-9,5.0e-10,6.4e-10,2.1e-15,1.0e5,3.2e7,3.0e7,3.0e7,0.0,0.0,0.0,0.0,4.3e-10,0.0,3.7e-8,8.9e-7,1.8e-7,3.0e-7,3.0e-7,4.3e-10,9.1e-7,8.9e-7]) # pre-exponential factors [cm^3/s] + B = np.array([0,0,0,0,0,0,0,0,0,0,0,0,0,0.74,0,0,0.51,0.61,0.51,0.51,0.74,0,0.51]) + A = np.array([0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1.59,2.61,0,0,0,0,1.59]) # activation temperature [eV] + dH = np.array([0.0,-7.412,-10.054,-7.336,0.0,0.0,0.0,0.0,0.0,11.548,11.624,12.907,15.76,-11.548,4.212,0.076,1.359,2.853,-1.283,-1.359,-11.624,-0.076,0.983]) # energy lost per electron due to ionization rxn [eV] + dEps = np.array([0.0,15.76,11.548,11.624,12.907,0.0]) # E, AR+, AR(m), AR(r), AR(4p), AR + + # BC parameters + # ks = 1.19e7 # electron recombination rate [cm/s] + ks = 1.366109824889323e7 # electron recombination rate [cm/s/eV] + + ################################################################### + # Constants of nature (probably shouldn't change unless you have + # root privileges on universe) + ################################################################### + qe = 1.6e-19 # unit charge [C] + eps0 = 8.86e-12 # permittivity of free space [F/m] + kB = 1.38e-23 # Boltzmann constant [J/K] + # kB = 8.62e-5 # Boltzmann constant [eV/K] + + + ################################################################### + # Calculate non-dimensional parameters + ################################################################### + + # 1) Convert input units to base SI (except eV) + nDe *= 100. # 1/(m*s) + nDi *= 100. # 1/(m*s) + nDm *= 100. + nDr *= 100. + nD4p *= 100. + nmue *= 100. # 1/(V*m*s) + nmui *= 100. # 1/(V*m*s) + Ck[0:5] *= 1e-6 # m^3/s + Ck[5:9] *= 1 # 1/s + Ck[9:22] *= 1e-6 # m^3/s + ks *= 0.01 # m/s + se *= 1.0e-20 # m^2 + + # 2) Compute "raw" transport parameters + De = nDe/nAr + Di = nDi/nAr + Dm = nDm/nAr + Dr = nDr/nAr + D4p = nD4p/nAr + + mue = nmue/nAr + mui = nmui/nAr + mum = nmum/nAr + mur = nmur/nAr + mu4p = nmu4p/nAr + + # 3) Compute non-dimensional properties required by solver + De = De*tau/(L*L) + Di = Di*tau/(L*L) + Dm = Dm*tau/(L*L) + Dr = Dr*tau/(L*L) + D4p = D4p*tau/(L*L) + + mue = mue*V0*tau/(L*L) + mui = mui*V0*tau/(L*L) + mum = mum*V0*tau/(L*L) + mur = mur*V0*tau/(L*L) + mu4p = mu4p*V0*tau/(L*L) + + Ck[0:4] *= tau*np0 + Ck[4] *= tau*nAr + Ck[5:9] *= tau + Ck[9:13] *= tau*nAr + Ck[13:] *= tau*np0 + A = A*1.5/e0 # 1.5 to convert from temperature to energy + dH = dH/e0 + qStar = V0/e0 # qe*V0/e0, since e0 in eV, need qe*V0 in eV, which is just V0 in V + alpha = qe*np0*L*L/(V0*eps0) + ks = ks*tau/L + p0 = p/qe/np0 + kappaB = 4.878171165833662*1.6129 # non-dimensional thermal conductivity of background specie + # (2/3)*tau/L**2*Kb/np0/kB, + # where Kb is the thermal conductivity of background specie + + params.beta = np.array([[0,1,1,1,0,0,0,0,0,1,1,1,2,1,2,1,1,2,1,1,1,1,1], # E + [0,1,1,1,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,0,0], # AR+ + [0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,1,0,1,0], # AR(m) + [0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,1,0,0,1,0,0,0,0], # AR(r) + [0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1], # AR(4p) + [2,1,1,1,2,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0]], dtype=np.int64) # AR + + params.alfa = np.array([[0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1], # E + [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], # AR+ + [2,1,0,2,1,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0], # AR(m) + [0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1], # AR(r) + [0,0,2,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,1,1,0,0,0], # AR(4p) + [0,0,0,0,1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0]], dtype=np.int64) # AR + # Rxn1: 2AR(m) -> 2AR + # Rxn2: AR(m) + AR(r) -> E + AR+ + AR + # Rxn3: 2AR(4p) -> E + AR+ + AR + # Rxn4: 2AR(m) -> E + AR+ + AR + # Rxn5: AR(m) + AR -> 2AR + # Rxn6: AR(r) -> AR + # Rxn7: AR(4p) -> AR + # Rxn8: AR(4p) -> AR(m) + # Rxn9: AR(4p) -> AR(r) + # Rxn10: E + AR -> E + AR(m) + # Rxn11: E + AR -> E + AR(r) + # Rxn12: E + AR -> E + AR(4p) + # Rxn13: E + AR -> 2E + AR+ + # Rxn14: E + AR(m) -> E + AR + # Rxn15: E + AR(m) -> 2E + AR+ + # Rxn16: E + AR(m) -> E + AR(r) + # Rxn17: E + AR(m) -> E + AR(4p) + # Rxn18: E + AR(4p) -> 2E + AR+ + # Rxn19: E + AR(4p) -> E + AR(r) + # Rxn20: E + AR(4p) -> E + AR(m) + # Rxn21: E + AR(r) -> E + AR + # Rxn22: E + AR(r) -> E + AR(m) + # Rxn23: E + AR(r) -> E + AR(4p) + + rxnNameDict = { 0: "2Ar(m) => 2Ar", + 1: "Ar(m) + Ar(r) => E + Ar+ + Ar", + 2: "2Ar(4p) => E + Ar+ + Ar", + 3: "2Ar(m) => E + Ar+ + Ar", + 4: "Ar(m) + Ar => 2Ar", + 5: "Ar(r) => Ar", + 6: "Ar(4p) => Ar", + 7: "Ar(4p) => Ar(m)", + 8: "Ar(4p) => Ar(r)", + 9: "lumped.metastable", + 10: "lumped.resonance", + 11: "lumped.2p", + 12: "ionization", + 13: "E + Ar(m) => E + Ar", + 14: "step_ionization", + 15: "E + Ar(m) => E + Ar(r)", + 16: "E + Ar(m) => E + Ar(4p)", + 17: "E + Ar(4p) => 2E + Ar+", + 18: "E + Ar(4p) => E + Ar(r)", + 19: "E + Ar(4p) => E + Ar(m)", + 20: "E + Ar(r) => E + Ar", + 21: "E + Ar(r) => E + Ar(m)", + 22: "E + Ar(r) => E + Ar(4p)"} + + # 4) Set values in params class + params.D[0] = De + params.D[1] = Di + params.D[2] = Dm + params.D[3] = Dr + params.D[4] = D4p + + params.mu[0] = mue + params.mu[1] = mui + params.mu[2] = mum + params.mu[3] = mur + params.mu[4] = mu4p + + params.A[:] = Ck[:] + params.B[:] = B[:] + params.C[:] = A[:] + + # Account for the 2/3 term to convert from electron temperature to electron energy + for i in range(len(params.A)): + params.A[i] *= (2/3)**(params.B[i]) + + params.dH[:] = dH[:] + params.dEps[:] = dEps[:] + params.qStar = qStar + params.alpha = alpha + params.ks = ks + params.gam = gam + params.kappaB = kappaB + params.nAronp0 = nAr / np0 + params.p0 = p0 + params.Tg0 = Tg0 + params.EC = 2.0 * me / mAr \ + * np.sqrt(16.0 * (me + mAr) * e0 * c**2 + / (3.0 * np.pi * me * mAr)) * se * nAr * tau + # params.EC = 2.0 * me / mAr * 3.8e9 * tau + + params.verticalShift = verticalShift / V0 + + # Parameters needed to compute the current with dimensions + params.V0Ltau = V0 / (L * tau) + params.V0L = V0 / L + params.LLV0tau = (L*L) / (V0*tau) + params.tauL = L / tau + params.np0 = np0 # "nominal" electron density [1/m^3] + params.qe = qe # unit charge [C] + params.eps0 = eps0 # unit charge [C] + params.eArea = electrodeArea # electrode area [m^2] + + + reactionExpressionslist = [f"{params.A[0]} * energy**{params.B[0]} * np.exp(-{params.C[0]} / energy)", + f"{params.A[1]} * energy**{params.B[1]} * np.exp(-{params.C[1]} / energy)", + f"{params.A[2]} * energy**{params.B[2]} * np.exp(-{params.C[2]} / energy)", + f"{params.A[3]} * energy**{params.B[3]} * np.exp(-{params.C[3]} / energy)", + f"{params.A[4]} * energy**{params.B[4]} * np.exp(-{params.C[4]} / energy)", + f"{params.A[5]} * energy**{params.B[5]} * np.exp(-{params.C[5]} / energy)", + f"{params.A[6]} * energy**{params.B[6]} * np.exp(-{params.C[6]} / energy)", + f"{params.A[7]} * energy**{params.B[7]} * np.exp(-{params.C[7]} / energy)", + f"{params.A[8]} * energy**{params.B[8]} * np.exp(-{params.C[8]} / energy)", + f"{params.A[9]} * energy**{params.B[9]} * np.exp(-{params.C[9]} / energy)", + f"{params.A[10]} * energy**{params.B[10]} * np.exp(-{params.C[10]} / energy)", + f"{params.A[11]} * energy**{params.B[11]} * np.exp(-{params.C[11]} / energy)", + f"{params.A[12]} * energy**{params.B[12]} * np.exp(-{params.C[12]} / energy)", + f"{params.A[13]} * energy**{params.B[13]} * np.exp(-{params.C[13]} / energy)", + f"{params.A[14]} * energy**{params.B[14]} * np.exp(-{params.C[14]} / energy)", + f"{params.A[15]} * energy**{params.B[15]} * np.exp(-{params.C[15]} / energy)", + f"{params.A[16]} * energy**{params.B[16]} * np.exp(-{params.C[16]} / energy)", + f"{params.A[17]} * energy**{params.B[17]} * np.exp(-{params.C[17]} / energy)", + f"{params.A[18]} * energy**{params.B[18]} * np.exp(-{params.C[18]} / energy)", + f"{params.A[19]} * energy**{params.B[19]} * np.exp(-{params.C[19]} / energy)", + f"{params.A[20]} * energy**{params.B[20]} * np.exp(-{params.C[20]} / energy)", + f"{params.A[21]} * energy**{params.B[21]} * np.exp(-{params.C[21]} / energy)", + f"{params.A[22]} * energy**{params.B[22]} * np.exp(-{params.C[22]} / energy)"] + + reactionTExpressionslist = [f"{params.A[0]} * (energy**({params.B[0]}-1)) * np.exp(-{params.C[0]}/energy) * ({params.B[0]} + {params.C[0]}/energy)", + f"{params.A[1]} * (energy**({params.B[1]}-1)) * np.exp(-{params.C[1]}/energy) * ({params.B[1]} + {params.C[1]}/energy)", + f"{params.A[2]} * (energy**({params.B[2]}-1)) * np.exp(-{params.C[2]}/energy) * ({params.B[2]} + {params.C[2]}/energy)", + f"{params.A[3]} * (energy**({params.B[3]}-1)) * np.exp(-{params.C[3]}/energy) * ({params.B[3]} + {params.C[3]}/energy)", + f"{params.A[4]} * (energy**({params.B[4]}-1)) * np.exp(-{params.C[4]}/energy) * ({params.B[4]} + {params.C[4]}/energy)", + f"{params.A[5]} * (energy**({params.B[5]}-1)) * np.exp(-{params.C[5]}/energy) * ({params.B[5]} + {params.C[5]}/energy)", + f"{params.A[6]} * (energy**({params.B[6]}-1)) * np.exp(-{params.C[6]}/energy) * ({params.B[6]} + {params.C[6]}/energy)", + f"{params.A[7]} * (energy**({params.B[7]}-1)) * np.exp(-{params.C[7]}/energy) * ({params.B[7]} + {params.C[7]}/energy)", + f"{params.A[8]} * (energy**({params.B[8]}-1)) * np.exp(-{params.C[8]}/energy) * ({params.B[8]} + {params.C[8]}/energy)", + f"{params.A[9]} * (energy**({params.B[9]}-1)) * np.exp(-{params.C[9]}/energy) * ({params.B[9]} + {params.C[9]}/energy)", + f"{params.A[10]} * (energy**({params.B[10]}-1)) * np.exp(-{params.C[10]}/energy) * ({params.B[10]} + {params.C[10]}/energy)", + f"{params.A[11]} * (energy**({params.B[11]}-1)) * np.exp(-{params.C[11]}/energy) * ({params.B[11]} + {params.C[11]}/energy)", + f"{params.A[12]} * (energy**({params.B[12]}-1)) * np.exp(-{params.C[12]}/energy) * ({params.B[12]} + {params.C[12]}/energy)", + f"{params.A[13]} * (energy**({params.B[13]}-1)) * np.exp(-{params.C[13]}/energy) * ({params.B[13]} + {params.C[13]}/energy)", + f"{params.A[14]} * (energy**({params.B[14]}-1)) * np.exp(-{params.C[14]}/energy) * ({params.B[14]} + {params.C[14]}/energy)", + f"{params.A[15]} * (energy**({params.B[15]}-1)) * np.exp(-{params.C[15]}/energy) * ({params.B[15]} + {params.C[15]}/energy)", + f"{params.A[16]} * (energy**({params.B[16]}-1)) * np.exp(-{params.C[16]}/energy) * ({params.B[16]} + {params.C[16]}/energy)", + f"{params.A[17]} * (energy**({params.B[17]}-1)) * np.exp(-{params.C[17]}/energy) * ({params.B[17]} + {params.C[17]}/energy)", + f"{params.A[18]} * (energy**({params.B[18]}-1)) * np.exp(-{params.C[18]}/energy) * ({params.B[18]} + {params.C[18]}/energy)", + f"{params.A[19]} * (energy**({params.B[19]}-1)) * np.exp(-{params.C[19]}/energy) * ({params.B[19]} + {params.C[19]}/energy)", + f"{params.A[20]} * (energy**({params.B[20]}-1)) * np.exp(-{params.C[20]}/energy) * ({params.B[20]} + {params.C[20]}/energy)", + f"{params.A[21]} * (energy**({params.B[21]}-1)) * np.exp(-{params.C[21]}/energy) * ({params.B[21]} + {params.C[21]}/energy)", + f"{params.A[22]} * (energy**({params.B[22]}-1)) * np.exp(-{params.C[22]}/energy) * ({params.B[22]} + {params.C[22]}/energy)"] + + + reactionExpressionTypelist = np.array([False,False,False,False,False,False,False,False,False, + True,True,True,True,False,True,True,True,False,True,True,False,True,True]) + + reactionsList = [] + LOGFilename = 'interpolationSample%s.log'%str(iSample) + f = open(LOGFilename, 'w') + # 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}/BOLSIGChemistry_6SpeciesRates/reaction300K_{1:d}.txt".format(root_dir, i), 'r') + #temp_file = open("{0:s}/BOLSIGChemistry_6SpeciesRates/reaction300K_Te.txt".format(root_dir), 'r') + + #rateCoeff = np.fromfile(rate_file) + #rateCoeff = np.genfromtxt(rate_file) + #rate_file.close() + #rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + + #Te = np.fromfile(temp_file) + #Te = np.genfromtxt(temp_file) + #temp_file.close() + #Te = np.reshape(Te,[Nsample, N300]).T[:,iSample] + + if (i < 15): + f = h5.File("../BOLSIGChemistry_6SpeciesRates/{0:s}.h5".format(rxnNameDict[i]), 'r') + dataset = f["table"] + else: + f = h5.File("../BOLSIGChemistry_6SpeciesRates/StepwiseExcitations.nominal.h5", 'r') + dataset = f[rxnNameDict[i]] + + Te = dataset[:,0] + Te /= 11604 + rateCoeff = dataset[:,1] + rateCoeff /= 6.022e23 + + # 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 + for k in range(len(Te)): + if Te[k] < 4.5 and indices[k] == False: + lastFalse = k + 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 < 13): + 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) + + logging.basicConfig(filename=LOGFilename) + logging.warning('Interpolation info (Reaction %s):', i + 1) + logging.warning('First non-zero entry in the rate coefficient: %s', I[0][0]) + logging.warning('Monotonicity of the rate coefficient start from entry: %s', lastFalse) + logging.warning('Position of possible duplicates in mean energy array: %s', + TeDuplicateindsForLog[0]) + + else: + rxn = eval("lambda energy :" + reactionExpressionslist[i]) + rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) + + reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]], + rxnBolsig = reactionExpressionTypelist[i], + kf = rxn, kf_T = rxn_T) + reactionsList.append(reaction) + + params.reactionsList = reactionsList + + # 5) Dump to screen + params.print() diff --git a/psaapProperties_6Species_Sampling.py b/psaapProperties_6Species_Sampling.py new file mode 100755 index 000000000..9abed73e7 --- /dev/null +++ b/psaapProperties_6Species_Sampling.py @@ -0,0 +1,411 @@ +import numpy as np +from scipy.interpolate import CubicSpline +import csv +import matplotlib.pyplot as plt +import matplotlib.colors as mcolors +import h5py as h5 +import logging + +class Reaction(object): + def __init__(self, *initial_data, **kwargs): + for dictionary in initial_data: + for key in dictionary: + setattr(self, key, dictionary[key]) + for key in kwargs: + setattr(self, key, kwargs[key]) + +def setPsaapProperties_6Species_Sampling(gam, inputV0, inputVDC, params, Nr, iSample): + """Sets non-dimensional properties corresponding to Liu 2014 paper. + + Inputs: + gam : Secondary electron emission coefficient + params : chebSolver.modelParams class + + Outputs: None + params data is overwritten using values from Liu 2014. + """ + ################################################################### + # User specified parameters (you may change these if you wish to + # run a different scenario from Liu 2014) + ################################################################### + + # densities + nAr = 3.22e22 # background number density of Ar [1/m^3] (corresponds to p=100 mTorr) + np0 = 8e16 # "nominal" electron density [1/m^3] + + # masses + # me = 9.10938356e-31 # mass of an electron [kg] + # me = 5.489e-4 # mass of an electron [u] + me = 0.511e6 # mass of an electron [eV/c2] + # mAr = 39.948 # mass of an argon atom [u] + # mAr = 39.948 * 1.66054e-27 # mass of an argon atom [kg] + mAr = 37.2158e9 # mass of an electron [eV/c2] + # u = 931.4941e6 # eV/c2 + c = 299792458 # speed of light [m/s] + se = 40 # momentum cross section [A^2] + + # nominal electron energy + e0 = 1.0 # [eV] + + # pressure + p = 133.3224*1.5 # [J/m^3] *1.5 to convert it to energy (1 Torr) + + # gas energy at the wall + Tg0 = 0.038778 # 3/2*300K*kB ~ (p0 - nT[:,0])/ntot + + # characteristics of driving voltage + V0 = inputV0 # amplitude of driving voltage [V] + verticalShift = inputVDC # DC voltage (vertical shift in driving voltage) + tau = (1./13.56e6) # period of driving voltage [s] + L = 2.00*0.005 # half-gap-width [m] (gap width is 2 cm) + electrodeArea = np.pi*0.05**2 # electrode area [m^2] (electrode diameter = 0.1 m) + + # transport parameters + nmue = 9.66e21 # argon number density times electron mobility [1/(V*cm*s)] + #nmui = 4.65e19 # argon number density times ion mobility [1/(V*cm*s)] + #nmum = 1 / (np.sqrt(16.0 * (mAr + mAr) * 300 * 8.62e-5 * c**2 + # / (3.0 * np.pi * mAr * mAr)) * se * mAr * 1.6e-19 / c**2) + nmum = 0.0 + nmui = 8.0e19 + #nmum = 9.35e19 + nDe = 3.86e22 # argon number density times electron diffusivity [1/(cm*s)] + nDi = 2.07e18 # argon number density times ion diffusivity [1/(cm*s)] + nDm = 2.42e18 # argon number density times metastable diffusivity [1/(cm*s)] + #nDm = 3.914e20 + + # reaction parameters (NB: k_i = Ck*Ee^B*exp(-A/Ee)) + # Ee = 3/2*Te (Te in eV) + # -> k_i = [Ck*(2/3)^B] * Ee^B * exp[-(3/2)*A/Ee] + # nominal + Ck = np.array([2.0e-7,2.1e-9,5.0e-10,6.4e-10,2.1e-15,1.0e-5,3.2e7,3.0e7,3.0e7,0.0,0.0,0.0,0.0,4.3e-10,0.0,3.7e-8,8.9e-7,1.8e-7,3.0e-7,3.0e-7,4.3e-10,9.1e-7,8.9e-7]) # pre-exponential factors [cm^3/s] + B = np.array([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.74,0.0,0.0,0.51,0.61,0.51,0.51,0.74,0.0,0.51]) + A = np.array([0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.59,2.61,0.0,0.0,0.0,0.0,1.59]) # activation temperature [eV] + dH = np.array([0.0,-7.412,-10.054,-7.336,0.0,0.0,0.0,0.0,0.0,11.548,11.624,12.907,15.76,-11.548,4.212,0.076,1.359,2.853,-1.283,-1.359,-11.624,-0.076,0.983]) # energy lost per electron due to ionization rxn [eV] + dEps = np.array([0.0,15.76,11.548,11.624,12.907,0.0]) + + # BC parameters + # ks = 1.19e7 # electron recombination rate [cm/s] + ks = 1.366109824889323e7 # electron recombination rate [cm/s/eV] + + ################################################################### + # Constants of nature (probably shouldn't change unless you have + # root privileges on universe) + ################################################################### + qe = 1.6e-19 # unit charge [C] + eps0 = 8.86e-12 # permittivity of free space [F/m] + kB = 1.38e-23 # Boltzmann constant [J/K] + # kB = 8.62e−5 # Boltzmann constant [eV/K] + + + ################################################################### + # Calculate non-dimensional parameters + ################################################################### + + # 1) Convert input units to base SI (except eV) + nDe *= 100. # 1/(m*s) + nDi *= 100. # 1/(m*s) + nDm *= 100. + nmue *= 100. # 1/(V*m*s) + nmui *= 100. # 1/(V*m*s) + Ck[0:5] *= 1e-6 # m^3/s + Ck[5:9] *= 1 # 1/s + Ck[9:22] *= 1e-6 # m^6/s + ks *= 0.01 # m/s + se *= 1.0e-20 # m^2 + + # 2) Compute "raw" transport parameters + De = nDe/nAr + Di = nDi/nAr + Dm = nDm/nAr + + mue = nmue/nAr + mui = nmui/nAr + mum = nmum/nAr + + # 3) Compute non-dimensional properties required by solver + De = De*tau/(L*L) + Di = Di*tau/(L*L) + Dm = Dm*tau/(L*L) + + mue = mue*V0*tau/(L*L) + mui = mui*V0*tau/(L*L) + mum = mum*V0*tau/(L*L) + + Ck[0:4] *= tau*np0 + Ck[4] *= tau*nAr + Ck[5:9] *= tau + Ck[9:13] *= tau*nAr + Ck[13:] *= tau*np0 + A = A*1.5/e0 # 1.5 to convert from temperature to energy + dH = dH/e0 + qStar = V0/e0 # qe*V0/e0, since e0 in eV, need qe*V0 in eV, which is just V0 in V + alpha = qe*np0*L*L/(V0*eps0) + ks = ks*tau/L + p0 = p/qe/np0 + kappaB = 4.878171165833662*1.6129 # non-dimensional thermal conductivity of background specie + # (2/3)*tau/L**2*Kb/np0/kB, + # where Kb is the thermal conductivity of background specie + + params.beta = np.array([[0,1,1,1,0,0,0,0,0,1,1,1,2,1,2,1,1,2,1,1,1,1,1], # E + [0,1,1,1,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,0,0,0], # AR+ + [0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,1,0,1,0], # AR(m) + [0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,1,0,0,1,0,0,0,0], # AR(r) + [0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,1], # AR(4p) + [2,1,1,1,2,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0]], dtype=np.int64) # AR + + params.alfa = np.array([[0,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1], # E + [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0], # AR+ + [2,1,0,2,1,0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0,0,0], # AR(m) + [0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1], # AR(r) + [0,0,2,0,0,0,1,1,1,0,0,0,0,0,0,0,0,1,1,1,0,0,0], # AR(4p) + [0,0,0,0,1,0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0,0,0]], dtype=np.int64) # AR + # Rxn1: 2AR(m) -> 2AR + # Rxn2: AR(m) + AR(r) -> E + AR+ + AR + # Rxn3: 2AR(4p) -> E + AR+ + AR + # Rxn4: 2AR(m) -> E + AR+ + AR + # Rxn5: AR(m) + AR -> 2AR + # Rxn6: AR(r) -> AR + # Rxn7: AR(4p) -> AR + # Rxn8: AR(4p) -> AR(m) + # Rxn9: AR(4p) -> AR(r) + # Rxn10: E + AR -> E + AR(m) + # Rxn11: E + AR -> E + AR(r) + # Rxn12: E + AR -> E + AR(4p) + # Rxn13: E + AR -> 2E + AR+ + # Rxn14: E + AR(m) -> E + AR + # Rxn15: E + AR(m) -> 2E + AR+ + # Rxn16: E + AR(m) -> E + AR(r) + # Rxn17: E + AR(m) -> E + AR(4p) + # Rxn18: E + AR(4p) -> 2E + AR+ + # Rxn19: E + AR(4p) -> E + AR(r) + # Rxn20: E + AR(4p) -> E + AR(m) + # Rxn21: E + AR(r) -> E + AR + # Rxn22: E + AR(r) -> E + AR(m) + # Rxn23: E + AR(r) -> E + AR(4p) + + rxnNameDict = {0: "2AR(m) => 2AR", + 1: "AR(m) + AR(r) => E + AR+ + AR", + 2: "2AR(4p) => E + AR+ + AR", + 3: "2AR(m) => E + AR+ + AR", + 4: "AR(m) + AR => 2AR", + 5: "AR(r) => AR", + 6: "AR(4p) => AR", + 7: "AR(4p) => AR(m)", + 8: "Ar(4p) => AR(r)", + 9: "1s-metastable", + 10: "1s-resonance", + 11: "2p-lumped", + 12: "Ionization", + 13: "E + AR(m) => E + AR", + 14: "StepIonization", + 15: "E + AR(m) => E + AR(r)", + 16: "E + AR(m) => E + AR(4p)", + 17: "E + AR(4p) => 2E + AR+", + 18: "E + AR(4p) => E + AR(r)", + 19: "E + AR(4p) => E + AR(m)", + 20: "E + AR(r) => E + AR", + 21: "E + AR(r) => E + AR(m)", + 22: "E + AR(r) => E + AR(4p)"} + + # 4) Set values in params class + params.D[0] = De + params.D[1] = Di + params.D[2] = Dm + + params.mu[0] = mue + params.mu[1] = mui + params.mu[2] = mum + + params.A[:] = Ck[:] + params.B[:] = B[:] + params.C[:] = A[:] + + # Account for the 2/3 term to convert from electron temperature to electron energy + for i in range(len(params.A)): + params.A[i] *= (2/3)**(params.B[i]) + + params.dH[:] = dH[:] + params.dEps[:] = dEps[:] + params.qStar = qStar + params.alpha = alpha + params.ks = ks + params.gam = gam + params.kappaB = kappaB + params.nAronp0 = nAr / np0 + params.p0 = p0 + params.Tg0 = Tg0 + params.EC = 2.0 * me / mAr \ + * np.sqrt(16.0 * (me + mAr) * e0 * c**2 + / (3.0 * np.pi * me * mAr)) * se * nAr * tau + # params.EC = 2.0 * me / mAr * 3.8e9 * tau + + params.verticalShift = verticalShift / V0 + + # Parameters needed to compute the current with dimensions + params.V0Ltau = V0 / (L * tau) + params.V0L = V0 / L + params.LLV0tau = (L*L) / (V0*tau) + params.tauL = L / tau + params.np0 = np0 # "nominal" electron density [1/m^3] + params.qe = qe # unit charge [C] + params.eps0 = eps0 # unit charge [C] + params.eArea = electrodeArea # electrode area [m^2] + + + reactionExpressionTypelist = np.array([False,False,False,False,False,False,False,False,False, + True,True,True,True,False,True,False,False,False,False,False,False,False,False]) + + reactionsList = [] + LOGFilename = 'interpolationSample%s.log'%str(iSample) + f = open(LOGFilename, 'w') + + for i in range(Nr): + Nsample = 20 + N300 = 200 + sample_root_dir = "../6species_SampleFiles" + if reactionExpressionTypelist[i]: + #Nsample = 1 + #N300 = 200 + fileString = sample_root_dir + "/" + rxnNameDict[i] + fileName = "%s.%08d.h5" % (fileString, iSample) + f = h5.File(fileName, 'r') + dataset = f["table"] + + #rateCoeff = np.fromfile(rate_file) + rateCoeff = dataset[:,1] + rateCoeff /= 6.022e23 + #rateCoeff = np.reshape(rateCoeff,[Nsample, N300]).T[:,iSample] + #rateCoeff = np.reshape(rateCoeff,[1, N300]).T[:,1] + + #Te = np.fromfile(temp_file) + Te = dataset[:,0] + Te /= 11604. + #Te = np.reshape(Te,[Nsample, N300]).T[:,iSample] + #Te = np.reshape(Te,[1, N300]).T[:,1] + + # 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 + for k in range(len(Te)): + if (Te[k] < 4.5 and indices[k] == False): + lasFalse = k + 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 == 9: + rateCoeffLog += - np.log(1.0/tau) + np.log(nAr) + elif i == 10: + rateCoeffLog += - np.log(1.0/tau) + np.log(nAr) + elif i == 11: + rateCoeffLog += - np.log(1.0/tau) + np.log(nAr) + elif i == 12: + 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) + + logging.basicConfig(filename=LOGFilename) + logging.warning('Interpolation info (Reaction %s):', i) + logging.warning('First non-zero entry in the rate coefficient: %s', I[0][0]) + logging.warning('Monotonicity of the rate coefficient start from entry: %s', lastFalse) + logging.warning('Position of possible duplicates in mean energy array: %s', + TeDuplicateindsForLog[0]) + + else: + fileName = "%s/Arrhenius.%08d.h5" % (sample_root_dir, iSample) + f = h5.File(fileName, 'r') + arrh_Coeffs = f[rxnNameDict[i]][...] + A, B, C = arrh_Coeffs + + # Non-dimensionalize the coefficients + if (i < 4): + A /= 6.022e23 + A *= tau*np0 + elif (i == 4): + A /= 6.022e23 + A *= tau*nAr + elif (i > 4 and i < 9): + A *= tau + elif (i > 8 and i < 13): + A /= 6.022e23 + A *= tau*nAr + else: + A /= 6.022e23 + A *= tau*np0 + + A *= 11604**B + C = C*1.5/e0 + + rxn = eval("lambda energy :" + f"{A} * energy**{B} * np.exp(-{C} / energy)") + rxn_T = eval("lambda energy :" + f"{A} * energy**({B}-1) * np.exp(-{C} / energy) * ({B} + {C} / energy)") + + #rxn = eval("lambda energy :" + reactionExpressionslist[i]) + #rxn_T = eval("lambda energy :" + reactionTExpressionslist[i]) + + reaction = Reaction(rxnAlfa = params.alfa[:,[i]], rxnBeta = params.beta[:,[i]], + rxnBolsig = reactionExpressionTypelist[i], + kf = rxn, kf_T = rxn_T) + reactionsList.append(reaction) + + params.reactionsList = reactionsList + + # 5) Dump to screen + params.print() diff --git a/reference_solns/Nominal_Np250_restart.npy b/reference_solns/Nominal_Np250_restart.npy new file mode 100644 index 000000000..de25b95b8 --- /dev/null +++ b/reference_solns/Nominal_Np250_restart.npy @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:80e433281f66962468805f299f8a05078790da53944450dbbdd8bb815196f21b +size 14128 diff --git a/reference_solns/Nominal_Np250_solution.npy b/reference_solns/Nominal_Np250_solution.npy new file mode 100644 index 000000000..59400c490 --- /dev/null +++ b/reference_solns/Nominal_Np250_solution.npy @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:460eb9c41df691b0ada34efb9146c74403422497f106b85fe008d5c008fbecbe +size 14128 diff --git a/runAllTests.sh b/runAllTests.sh new file mode 100755 index 000000000..33b15fbd2 --- /dev/null +++ b/runAllTests.sh @@ -0,0 +1,7 @@ +#!/bin/bash + +./testRunBase.sh || echo "testRunBase.sh failed" +./testRunCN.sh +./testInterpTrans.sh +./testRunBackground.sh +./test6Species.sh diff --git a/sampleRun.sh b/sampleRun.sh index 053ff5dc1..af1297760 100755 --- a/sampleRun.sh +++ b/sampleRun.sh @@ -6,8 +6,9 @@ error_exit() exit 1 } -EXE="python3 ./chebSolver.py" -NEWTEXE="python3 ./timePeriodicSolver.py" +module load python/3.8.2 +EXE="python3 ../chebSolver.py" +NEWTEXE="python3 ../timePeriodicSolver.py" # psaap 2 species case #Np=250 @@ -26,11 +27,11 @@ NEWTEXE="python3 ./timePeriodicSolver.py" #saveFile="newton_3spec_Np${Np}_fullsoln.npy" # Liu 4 species case -Np=250 +Np=150 Nt=25600 Nt1=128 dt=0.0078125 -scenario=2 +scenario=6 baseFile="restart_4spec_Np${Np}_" newtFile="newton_4spec_Np${Np}.npy" saveFile="newton_4spec_Np${Np}_fullsoln.npy" diff --git a/sampleRunCN.sh b/sampleRunCN.sh index 34df8591e..fb6a9f314 100755 --- a/sampleRunCN.sh +++ b/sampleRunCN.sh @@ -6,6 +6,7 @@ error_exit() exit 1 } +module load python/3.8.2 EXE="python3 ./chebSolver.py" NEWTEXE="python3 ./timePeriodicSolver.py" @@ -26,11 +27,11 @@ NEWTEXE="python3 ./timePeriodicSolver.py" #saveFile="newton_3spec_CN_Np${Np}_fullsoln.npy" # Liu 4 species case -Np=300 +Np=150 Nt=25600 Nt1=128 dt=0.0078125 -scenario=2 +scenario=6 baseFile="restart_4spec_CN_Np${Np}_" newtFile="newton_4spec_CN_Np${Np}.npy" saveFile="newton_4spec_CN_Np${Np}_fullsoln.npy" diff --git a/test6Species.sh b/test6Species.sh new file mode 100755 index 000000000..2a4e0ca72 --- /dev/null +++ b/test6Species.sh @@ -0,0 +1,42 @@ +#!/bin/bash + +error_exit() +{ + echo "$1" 1>&2 + exit 1 +} + +tmp_dir=$(mktemp -d -t tmp-test-XXXXXXXXXX --tmpdir=.) +echo "Running 6 species test in $tmp_dir" +cd $tmp_dir +EXE="python3 ../chebSolver.py" +NEWTEXE="python3 ../timePeriodicSolver.py" + +# 6 species + 23 rxn - Nominal Rates case +Np=250 +Nt=25600 +Nt1=128 +dt=0.0078125 +scenario=12 +baseFile="restart_6spec_CN_Np${Np}_" +newtFile="newton_6spec_CN_Np${Np}.npy" +saveFile="newton_6spec_CN_Np${Np}_fullsoln.npy" + +baseCmd="$EXE --Np $Np --Nt $Nt --dt $dt --scenario $scenario --EinsteinForm --elasticCollisionActivation --backgroundSpecieActivation" +newtCmd="$NEWTEXE --Np $Np --Nt $Nt1 --Nn 20 --scenario $scenario --tscheme CN --alpha0 0.25 --increaseFac 2.0 --EinsteinForm --elasticCollisionActivation --backgroundSpecieActivation" +saveCmd="$EXE --Np $Np --Nt $Nt1 --dt $dt --scenario $scenario --tscheme CN --EinsteinForm --elasticCollisionActivation --backgroundSpecieActivation" + +diffCmd="python3 ../diffSolns.py --reference ../reference_solns/Nominal_Np250_solution.npy --Np $Np --Nv 7" + + +screenOut="run6species.out" +rm -f $screenOut + +$newtCmd --V0 100 --VDC 0.0 --gam 0.01 --rtol 1e-8 --restart "../reference_solns/Nominal_Np250_restart.npy" \ + --outfile $newtFile >> $screenOut || error_exit "Shooting failed" + +$diffCmd --solution $newtFile >> $screenOut || error_exit "Solution differs from reference" + +# if test passed, delete tmp dir +cd .. +rm -rf $tmp_dir diff --git a/timePeriodicSolver.py b/timePeriodicSolver.py index 050861119..cb239002d 100644 --- a/timePeriodicSolver.py +++ b/timePeriodicSolver.py @@ -7,13 +7,13 @@ class timePeriodicSolver: def __init__(self, Ns, NT, Np, elasticCollisionActivationFactor, backgroundSpecieActivationFactor, EinsteinForm, gam, V0, VDC, restart=None, scenario=0, scheme='BE', - alpha0 = 1.0, increaseFac = 1.0): + alpha0 = 1.0, increaseFac = 1.0, iSample = 0): self.tds = cs.timeDomainCollocationSolver(Ns, NT, Np, elasticCollisionActivationFactor, backgroundSpecieActivationFactor, EinsteinForm, gam, V0, - VDC, scenario, scheme) + VDC, scenario, scheme, iSample) self.res = np.zeros((self.tds.Ndof,1)) self.jac = np.zeros((self.tds.Ndof,self.tds.Ndof)) @@ -141,7 +141,8 @@ def solveNewtonStep(self, Uic, Nt): type=float, help='Newton step under-relaxation factor') parser.add_argument('--increaseFac', metavar='increaseFac', default=1.0, type=float, help='Increase alpha by this factor each step') - + parser.add_argument('--iSample', metavar='iSample', default=0, + type=int, help='Sample index, if BOLSIG chemistry is used.') args = parser.parse_args() # Dump inputs to the screen for posterity @@ -181,11 +182,29 @@ def solveNewtonStep(self, Uic, Nt): elif(args.scenario==4): print("# Running scenario = 4 (4 species, 8 rxn, Liu 2017)") Ns = 4 + elif(args.scenario==5): + print("# Running scenario = 5 (4 species, 7 rxn, Bolsing and Lay, Moss et al, 2003)") + Ns = 4 + elif(args.scenario==6): + print("# Running scenario = 6 (4 species, 9 rxn, Juan's Mechanism)") + Ns = 4 + elif(args.scenario==7): + print("# Running scenario = 7 (4 species, 9 rxn, Nominal Reaction Rates)") + Ns = 4 + elif(args.scenario==8): + Ns = 4 + elif(args.scenario==9): + print("# Running scenario = 9 (6 species, 23 rxn)") + Ns = 6 + elif(args.scenario==10): + Ns = 6 + elif(args.scenario==12): + Ns = 6 elif(args.scenario==21): print("# Running scenario = 21 (4 species, 8 rxn, Liu 2017, interpolated transport)") Ns = 4 else: - print("ERROR: Scenario = {0:d} not recognized. Use --scenario i with i=0, 1, 2, or 3. Exiting.".format(args.scenario)) + print("ERROR: Scenario = {0:d} not recognized. Exiting.".format(args.scenario)) exit(-1) print("#") @@ -220,7 +239,8 @@ def solveNewtonStep(self, Uic, Nt): args.gam, args.V0, args.VDC, restart=args.restart, scenario=args.scenario, scheme=args.tscheme, - alpha0 = args.alpha0, increaseFac = args.increaseFac) + alpha0 = args.alpha0, increaseFac = args.increaseFac, + iSample = args.iSample) # Get the IC, for use in computing the residual below Uic = np.copy(tps.tds.U1)