diff --git a/MFC_load_range.asv b/MFC_load_range.asv new file mode 100644 index 0000000..27e07a7 --- /dev/null +++ b/MFC_load_range.asv @@ -0,0 +1,52 @@ +clear all +close all + +%% +mechanism={'MFC'}; +% fuel_sim={'base','v1','v2','v3' }; +fuel_sim={'modify'}; +fuel_name={'n_dodecane'}; +date = {'03_06_2017_'}; +equi =1; +% pressure=[20 40 ]; +pressure = [20 40]; +for k=1:length(pressure) + pressure_text{k}=['P',num2str(pressure(k)),'atm']; +end + +classnumb=[15 26 27]; +numbOfClass = length(classnumb); +class_numb_text = {}; +for k=1:numbOfClass + class_numb_text=[class_numb_text ['class',num2str(classnumb(k))]]; +end +num_cases=[26 11]; +num_cases_text = {}; +for k=1:num_cases(1) + num_cases_text=[num_cases_text ['numcases',num2str(k)]]; +end +%% + +for i=1:length(mechanism) + for j=1:length(fuel_sim) + for k = 1: length(classnumb) +% phi=1 simulation results + for m=1:length(pressure) + for q=1:num_cases(m) + location=['C:\Users\unghee\Dropbox\post_process','\',mechanism{i},'\',fuel_name{1},'_',num2str(pressure(m))... + ,'atm','_','phi',num2str(equi),'_',date{1},'\',fuel_sim{j},'\',class_numb_text{k},'\',num2str(q)]; + sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{m}).(num_cases_text{q})=read_ignition_delay(location,7); + end + end + + end + end +end + + +directory=['C:\Users\unghee\Dropbox\post_process','\',mechanism{1},'\',fuel_name{1},'_',num2str(pressure(1))... + ,'atm','_','phi',num2str(equi),'_',date{1}]; +cd(directory); +cd(fuel_sim{1}); +filename = ['simulation_result','_',fuel_sim{1},'_',num2str(pressure(1)),'atm','_','phi','_',num2str(equi),'.mat']; +save(filename,'sim') diff --git a/MFC_load_range.m b/MFC_load_range.m new file mode 100644 index 0000000..9a95a25 --- /dev/null +++ b/MFC_load_range.m @@ -0,0 +1,52 @@ +clear all +close all + +%% +mechanism={'MFC'}; +% fuel_sim={'base','v1','v2','v3' }; +fuel_sim={'modify'}; +fuel_name={'n_dodecane'}; +date = {'03_08_2017'}; +equi =1; +% pressure=[20 40 ]; +pressure = [20]; +for k=1:length(pressure) + pressure_text{k}=['P',num2str(pressure(k)),'atm']; +end + +classnumb=[15 22 26 27 28]; +numbOfClass = length(classnumb); +class_numb_text = {}; +for k=1:numbOfClass + class_numb_text=[class_numb_text ['class',num2str(classnumb(k))]]; +end +num_cases=[26 11]; +num_cases_text = {}; +for k=1:num_cases(1) + num_cases_text=[num_cases_text ['numcases',num2str(k)]]; +end +%% + +for i=1:length(mechanism) + for j=1:length(fuel_sim) + for k = 1: length(classnumb) +% phi=1 simulation results + for m=1:length(pressure) + for q=1:num_cases(m) + location=['C:\Users\unghee\Dropbox\post_process','\',mechanism{i},'\',fuel_name{1},'_',num2str(pressure(m))... + ,'atm','_','phi',num2str(equi),'_',date{1},'\',fuel_sim{j},'\',class_numb_text{k},'\',num2str(q)]; + sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{m}).(class_numb_text{k}).(num_cases_text{q})=read_ignition_delay(location,7); + end + end + + end + end +end + + +directory=['C:\Users\unghee\Dropbox\post_process','\',mechanism{1},'\',fuel_name{1},'_',num2str(pressure(1))... + ,'atm','_','phi',num2str(equi),'_',date{1}]; +cd(directory); +cd(fuel_sim{1}); +filename = ['simulation_result','_',fuel_sim{1},'_',num2str(pressure(1)),'atm','_','phi','_',num2str(equi),'.mat']; +save(filename,'sim') diff --git a/README.md b/README.md new file mode 100644 index 0000000..0528377 --- /dev/null +++ b/README.md @@ -0,0 +1,9 @@ +# Chemical Mechanisms Optimizer +5/08/2017 Jordan Lee (unghee@umich.edu)1 +#About Mechanisms Optimizer +The optimizer optimizes reaction rate parameters(pre-exponential and activation energyterms) in the selected reactions in the chemicalmechanism in Chemkin format. The optimizerruns in MATLAB framework. Chemkin homogeneous reactor simulationsare usedto calculate 0-D ignition delay time. + +# 2 Files description in the optimizer +The homogeneous reactor code in Chemkin is used. For a detailed description on this code, please refer to CHEMKIN ignition delay calculations manual. The flow chart of the entire process is as below.In whole, the optimizer is composed of sensitivity analysis sectionand + + diff --git a/Ra_Reitz_load_data.m b/Ra_Reitz_load_data.m new file mode 100644 index 0000000..3d0c7c8 --- /dev/null +++ b/Ra_Reitz_load_data.m @@ -0,0 +1,39 @@ +clear all +close all + +%% +mechanism={'MFC'}; + +fuel_sim={'afteroptimize_3'}; +fuel_name={'JetA'}; +date = {'05_03_2017'}; +equi =1; + +pressure = [40]; +for k=1:length(pressure) + pressure_text{k}=['P',num2str(pressure(k)),'atm']; +end +num_cases=25; + + +for i=1:length(mechanism) + for j=1:length(fuel_sim) + +% phi=1 simulation results + for k=1:length(pressure) + location=['C:\Users\unghee\Dropbox\post_process','\',mechanism{i},'\',fuel_name{1},'_',num2str(pressure(k)),'atm','_','phi',num2str(equi),'_',date{1},'\',fuel_sim{j}]; + + sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k})=read_ignition_delay(location,num_cases); + + end + + + end +end + +directory=['C:\Users\unghee\Dropbox\post_process','\',mechanism{1},'\',fuel_name{1},'_',num2str(pressure(1)),'atm','_','phi',num2str(equi),'_',date{1}]; + +cd(directory); +cd(fuel_sim{1}); +filename = ['simulation_result','_',fuel_sim{1},'_',num2str(pressure(1)),'atm','_','phi','_',num2str(equi),'.mat']; +save(filename,'sim') diff --git a/Ra_Reitz_plot_2.asv b/Ra_Reitz_plot_2.asv new file mode 100644 index 0000000..331e208 --- /dev/null +++ b/Ra_Reitz_plot_2.asv @@ -0,0 +1,187 @@ +clear all +% close all +%% Ra&Reitz + +%%%%%% setting %%%%%%%%%%%%%%%%%%%%%%%%%%% +mechanism={'MFC' }; + +pressure=[40]; +phi=1; +iteration_numb=3; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +for k=1:length(pressure) + pressure_text{k}=['P',num2str(pressure(k)),'atm']; + + simulation_result_beforeoptimize_text{k}=['simulation_result_beforeoptimize_',num2str(pressure(k)),'atm','_',... + 'phi','_',num2str(phi),'.mat']; + for i= 1 : length(iteration_numb) + simulation_result_afteroptimize_text{k}=['simulation_result_afteroptimize_',num2str(iteration_numb(i)),'_',... + num2str(pressure(k)),'atm','_','phi','_',num2str(phi),'.mat']; + end + fuel_sim_text = ['afteroptimize','_',num2str(iteration_numb(i))]; +end + +num_cases=25; + + +pure_component_ID; +real_fuel_ID; + +exp_markers={'ro' 'g*' 'b^'}; +sim_line={'k-' 'k-s','k-d','k-^','k-x','k-o','k-p'}; +marker_size=8; +line_width=2; + +i=1; %mechanism + + +num_cases=[26 11]; +num_cases_text = {}; +for k=1:num_cases(1) + num_cases_text=[num_cases_text ['numcases',num2str(k)]]; +end +% classnumb=[22 24 26 28]; +% numbOfClass = length(classnumb); +% class_numb_text = {}; +% for k=1:numbOfClass +% class_numb_text=[class_numb_text ['class',num2str(classnumb(k))]]; +% end +for k = 1: length(pressure) + +h3=figure('position',[20 50 580 480]); +set(gca,'Fontsize',13) +clear legend_text +% % legend_text{1}='Vasu et al.40atm, Shock Tube'; +% legend_text{1}='Shen et al.40atm, Shock Tube'; +% % semilogy(Vasu_40atm(:,2),Vasu_40atm(:,5)/1000,exp_markers{2},'markersize',marker_size) +% semilogy(Shen_dode_40atm(:,2),Shen_dode_40atm(:,5),exp_markers{2},'markersize',marker_size) + +% legend_text{1}=['ERC',pressure_text{k},'Dummy']; +% legend_text{2}=['ERC',pressure_text{k},'target(original)']; +% legend_text{3}=['ERC',pressure_text{k},'after optimize']; + +legend_text{1}=['MFC ',pressure_text{k},', ','beforeoptimize']; + if strcmp(pressure_text{k},'P20atm') +% legend_text{2}=['Vasu, ',pressure_text{k},', ','target']; + legend_text{3}=['exp, ',pressure_text{k},', ','target']; + + else +% legend_text{2}=['Shen, ',pressure_text{k},', ','target']; + legend_text{3}=['exp, ',pressure_text{k},', ','target']; + end + +legend_text{2}=['MFC ',pressure_text{k},', after optimize']; + + + load(simulation_result_beforeoptimize_text{1}) + fuel_sim={'beforeoptimize'}; + + for j=1:length(fuel_sim) + semilogy(sim.(mechanism{i}).(fuel_sim{1}).(pressure_text{k}).table.data(:,6),... + sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).table.data(:,10),... + 'k--','markersize',marker_size) + end + + + hold on + + + + + load(simulation_result_afteroptimize_text{1}) + fuel_sim={fuel_sim_text}; + + + + semilogy(sim.(mechanism{1}).(fuel_sim{1}).(pressure_text{k}).table.data(:,6),... + sim.(mechanism{1}).(fuel_sim{1}).(pressure_text{k}).table.data(:,10),... + 'b','markersize',marker_size) + + +ylim([0 10^5]); + + + +% hold on + + if strcmp(pressure_text{k},'P20atm') +% semilogy(Vasu_dode_20atm(:,2),Vasu_dode_20atm(:,5),'r*','markersize',marker_size) % pure dodecane + semilogy(Vasu_20atm(:,2),Vasu_20atm(:,5),'r^','markersize',marker_size,'markerface','r') % Jet A + hold on; + semilogy(Wang_20atm(:,2),Wang_20atm(:,5),'r^','markersize',marker_size,'markerface','r') % Jet A + hold on; + semilogy(Dooley_20atm(:,2),Dooley_20atm(:,5),'r^','markersize',marker_size,'markerface','r') % Jet A + else +% semilogy(Shen_dode_40atm(:,2),Shen_dode_40atm(:,5),'r*','markersize',marker_size)% pure dodecane + semilogy(Wang_40atm(:,2),Wang_40atm(:,5),'rd','markersize',marker_size,'markerface','r') % Jet A + end +%% upper bound & lowerbond +% The below shows the range of the variation + +% hold on +% +% +% load('simulation_result_modify_20atm_phi_1.mat') +% fuel_sim={'modify'}; +% +% +% % for j=1:length(fuel_sim) +% % for q= 1: num_cases(k) +% % semilogy(sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(num_cases_text{q}).table.data(:,6),... +% % sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(num_cases_text{q}).table.data(:,10),... +% % 'bo','markersize',marker_size); +% % end +% % end +% % hold on; +% +% x.(class_numb_text{1})=[]; y.(class_numb_text{1})=[]; +% for j=1:length(fuel_sim) +% for q= 1: num_cases(k) +% for m = 1: length(classnumb) +% x.(class_numb_text{m})(q)=sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(class_numb_text{m}).(num_cases_text{q}).table.data(1,6); +% y.(class_numb_text{m})(q)=sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(class_numb_text{m}).(num_cases_text{q}).table.data(4,10); +% end +% +% end +% end +% for m = 1: length(classnumb) +% % semilogy(x.(class_numb_text{m}),y.(class_numb_text{m}),'k','markersize',marker_size,'LineWidth',1.5); +% hold on; +% end +% x.(class_numb_text{1})=[]; y.(class_numb_text{1})=[]; +% for j=1:length(fuel_sim) +% for q= 1: num_cases(k) +% for m = 1: length(classnumb) +% x.(class_numb_text{m})(q)=sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(class_numb_text{m}).(num_cases_text{q}).table.data(1,6); +% y.(class_numb_text{m})(q)=sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(class_numb_text{m}).(num_cases_text{q}).table.data(7,10); +% end +% +% end +% end +% for m = 1: length(classnumb) +% % semilogy(x.(class_numb_text{m}),y.(class_numb_text{m}),'k','markersize',marker_size,'LineWidth',1.5); +% hold on; +% end +%% + legend(legend_text,'location','SouthEast','interpreter','none') + legend('boxoff') + xlabel('1000/T (1/K)') + ylabel('Ignition Delay Time (micros)') + % axis([0.7 1.6 20 40000]) + % clear legend_text + annotation(h3,'textbox',[0.413 0.38 0.279 0.05],... + 'String',{'JetA',... + '/Air', '\phi=1'},... + 'FontSize',13,... + 'FontName','Arial',... + 'FitBoxToText','off',... + 'LineStyle','none'); + annotation(h3,'textbox',[0.25 0.8 0.279 0.05],... + 'String',{pressure_text{k}},... + 'FontSize',13,... + 'FontName','Arial',... + 'FitBoxToText','off',... + 'LineStyle','none'); + + +end diff --git a/Ra_Reitz_plot_2.m b/Ra_Reitz_plot_2.m new file mode 100644 index 0000000..275370d --- /dev/null +++ b/Ra_Reitz_plot_2.m @@ -0,0 +1,173 @@ +clear all +% close all +% This function is used to plot the optimized ignition delay time results +% and compare with the target(experimental) data. + +%%%%%% setting %%%%%%%%%%%%%%%%%%%%%%%%%%% +mechanism={'MFC' }; +pressure=[40]; +phi=1; +iteration_numb=3; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +for k=1:length(pressure) + pressure_text{k}=['P',num2str(pressure(k)),'atm']; + + simulation_result_beforeoptimize_text{k}=['simulation_result_beforeoptimize_',num2str(pressure(k)),'atm','_',... + 'phi','_',num2str(phi),'.mat']; + for i= 1 : length(iteration_numb) + simulation_result_afteroptimize_text{k}=['simulation_result_afteroptimize_',num2str(iteration_numb(i)),'_',... + num2str(pressure(k)),'atm','_','phi','_',num2str(phi),'.mat']; + end + fuel_sim_text = ['afteroptimize','_',num2str(iteration_numb(i))]; +end + +num_cases=25; % 25 temperature points + +pure_component_ID; +real_fuel_ID; + +exp_markers={'ro' 'g*' 'b^'}; +sim_line={'k-' 'k-s','k-d','k-^','k-x','k-o','k-p'}; +marker_size=8; +line_width=2; + +i=1; %mechanism + +num_cases=[26 11]; +num_cases_text = {}; +for k=1:num_cases(1) + num_cases_text=[num_cases_text ['numcases',num2str(k)]]; +end + +for k = 1: length(pressure) + +h3=figure('position',[20 50 580 480]); +set(gca,'Fontsize',13) +clear legend_text + +%% defines the legend +legend_text{1}=['MFC ',pressure_text{k},', ','beforeoptimize']; + if strcmp(pressure_text{k},'P20atm') +% legend_text{2}=['Vasu, ',pressure_text{k},', ','target']; + legend_text{3}=['exp, ',pressure_text{k},', ','target']; + + else +% legend_text{2}=['Shen, ',pressure_text{k},', ','target']; + legend_text{3}=['exp, ',pressure_text{k},', ','target']; + end + +legend_text{2}=['MFC ',pressure_text{k},', after optimize']; +%% loading the ignition delay times + +% loads the base mechanism's ignition delay time result + load(simulation_result_beforeoptimize_text{1}) + fuel_sim={'beforeoptimize'}; + + for j=1:length(fuel_sim) + semilogy(sim.(mechanism{i}).(fuel_sim{1}).(pressure_text{k}).table.data(:,6),... + sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).table.data(:,10),... + 'k--','markersize',marker_size) + end + + + hold on + + +% loads the optimized mechanism's ignition delay time result + load(simulation_result_afteroptimize_text{1}) + fuel_sim={fuel_sim_text}; + + + + semilogy(sim.(mechanism{1}).(fuel_sim{1}).(pressure_text{k}).table.data(:,6),... + sim.(mechanism{1}).(fuel_sim{1}).(pressure_text{k}).table.data(:,10),... + 'b','markersize',marker_size) + +% fixes the axes +ylim([0 10^5]); + + +% loads the target(experimental data) ignition delay time result + if strcmp(pressure_text{k},'P20atm') +% semilogy(Vasu_dode_20atm(:,2),Vasu_dode_20atm(:,5),'r*','markersize',marker_size) % pure dodecane + semilogy(Vasu_20atm(:,2),Vasu_20atm(:,5),'r^','markersize',marker_size,'markerface','r') % Jet A + hold on; + semilogy(Wang_20atm(:,2),Wang_20atm(:,5),'r^','markersize',marker_size,'markerface','r') % Jet A + hold on; + semilogy(Dooley_20atm(:,2),Dooley_20atm(:,5),'r^','markersize',marker_size,'markerface','r') % Jet A + else +% semilogy(Shen_dode_40atm(:,2),Shen_dode_40atm(:,5),'r*','markersize',marker_size)% pure dodecane + semilogy(Wang_40atm(:,2),Wang_40atm(:,5),'rd','markersize',marker_size,'markerface','r') % Jet A + end +%% upper bound & lowerbond +% The below shows the range of the variation. This can be use to check the +% range of ignition time when the rate parameters are varied(7 variations). + + +% hold on +% +% +% load('simulation_result_modify_20atm_phi_1.mat') +% fuel_sim={'modify'}; +% +% +% % for j=1:length(fuel_sim) +% % for q= 1: num_cases(k) +% % semilogy(sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(num_cases_text{q}).table.data(:,6),... +% % sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(num_cases_text{q}).table.data(:,10),... +% % 'bo','markersize',marker_size); +% % end +% % end +% % hold on; +% +% x.(class_numb_text{1})=[]; y.(class_numb_text{1})=[]; +% for j=1:length(fuel_sim) +% for q= 1: num_cases(k) +% for m = 1: length(classnumb) +% x.(class_numb_text{m})(q)=sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(class_numb_text{m}).(num_cases_text{q}).table.data(1,6); +% y.(class_numb_text{m})(q)=sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(class_numb_text{m}).(num_cases_text{q}).table.data(4,10); +% end +% +% end +% end +% for m = 1: length(classnumb) +% % semilogy(x.(class_numb_text{m}),y.(class_numb_text{m}),'k','markersize',marker_size,'LineWidth',1.5); +% hold on; +% end +% x.(class_numb_text{1})=[]; y.(class_numb_text{1})=[]; +% for j=1:length(fuel_sim) +% for q= 1: num_cases(k) +% for m = 1: length(classnumb) +% x.(class_numb_text{m})(q)=sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(class_numb_text{m}).(num_cases_text{q}).table.data(1,6); +% y.(class_numb_text{m})(q)=sim.(mechanism{i}).(fuel_sim{j}).(pressure_text{k}).(class_numb_text{m}).(num_cases_text{q}).table.data(7,10); +% end +% +% end +% end +% for m = 1: length(classnumb) +% % semilogy(x.(class_numb_text{m}),y.(class_numb_text{m}),'k','markersize',marker_size,'LineWidth',1.5); +% hold on; +% end +%% setting the position of the figure + legend(legend_text,'location','SouthEast','interpreter','none') + legend('boxoff') + xlabel('1000/T (1/K)') + ylabel('Ignition Delay Time (micros)') + % axis([0.7 1.6 20 40000]) + % clear legend_text + annotation(h3,'textbox',[0.413 0.38 0.279 0.05],... + 'String',{'JetA',... + '/Air', '\phi=1'},... + 'FontSize',13,... + 'FontName','Arial',... + 'FitBoxToText','off',... + 'LineStyle','none'); + annotation(h3,'textbox',[0.25 0.8 0.279 0.05],... + 'String',{pressure_text{k}},... + 'FontSize',13,... + 'FontName','Arial',... + 'FitBoxToText','off',... + 'LineStyle','none'); + + +end diff --git a/bash_backup/afteroptimize.pbs b/bash_backup/afteroptimize.pbs new file mode 100644 index 0000000..11e74e5 --- /dev/null +++ b/bash_backup/afteroptimize.pbs @@ -0,0 +1,72 @@ +#PBS -A jmartz_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=00:20:00,pmem=2500mb,procs=1 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 1-25 +#PBS -joe + +cd $PBS_O_WORKDIR + + +##setting +pressure="40" +mechanism="MFC"; +fuel_name="n_dodecane"; +#class_name="class11_class15_class22_class24_class26_class27_class28"; +iteration_numb="3"; +equi="1"; +matlabdate="04_12_2017" +## + +## count number of file +#cd /scratch/engin_flux/unghee/chemkin/Ra_Reitz/${fuel_name}_${pressure}atm_phi${equi}_$matlabdate +#numfiles=(*) +#numfiles=${#numfiles[@]} +#echo "number of temperature points" +#echo $numfiles +## + +## access and make file +#cd /scratch/engin_flux/unghee/chemkin/Ra_Reitz +#now=$(date +"%m_%d_%Y") +#mkdir ${fuel_name}_${pressure} +#cd ${fuel_name}_${pressure} + + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + + +mkdir afteroptimize_${iteration_numb} +cd afteroptimize_${iteration_numb} +## + + + +##calculate chemkin + + echo "calcultating temp ${PBS_ARRAYID}" + mkdir ck_${PBS_ARRAYID} + cd ck_${PBS_ARRAYID} + + + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/result/${PBS_ARRAYID}/${mechanism}_${fuel_name}.inp ./ +# cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/output/${mechanism}_${fuel_name}_${class_name}_afteroptimize.inp ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/output/${mechanism}_afteroptimize_${iteration_numb}_iteration.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_therm.dat ./ + cp /scratch/engin_flux/unghee/chemkin/chemkindata.dtd ./ + + +# $CHEMKIN_BIN/chem -i ${mechanism}_${fuel_name}_afteroptimize_${iteration_numb}_iteration.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/chem -i ${mechanism}_afteroptimize_${iteration_numb}_iteration.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/CKReactorGenericClosed -i ${mechanism}_${fuel_name}.inp -o ${mechanism}_reg.out + $CHEMKIN_BIN/GetSolution + $CHEMKIN_BIN/CKSolnTranspose + cd .. + + + + diff --git a/bash_backup/beforeoptimize.pbs b/bash_backup/beforeoptimize.pbs new file mode 100644 index 0000000..2ba7440 --- /dev/null +++ b/bash_backup/beforeoptimize.pbs @@ -0,0 +1,57 @@ +#PBS -A engin_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=00:10:00,pmem=2500mb,procs=1 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 1-25 +#PBS -joe + +cd $PBS_O_WORKDIR + + +##setting +pressure="40" +mechanism="MFC"; +fuel_name="n_dodecane"; +equi="1"; +matlabdate="03_23_2017_1_iteration" +## + + +## access and make file + + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + + +mkdir beforeoptimize +cd beforeoptimize +## + + + +##calculate chemkin + + echo " calcultating temp ${PBS_ARRAYID} " + mkdir ck_${PBS_ARRAYID} + cd ck_${PBS_ARRAYID} + + + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/result/${PBS_ARRAYID}/${mechanism}_${fuel_name}.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_base.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_therm.dat ./ + cp /scratch/engin_flux/unghee/chemkin/chemkindata.dtd ./ + + + $CHEMKIN_BIN/chem -i ${mechanism}_base.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/CKReactorGenericClosed -i ${mechanism}_${fuel_name}.inp -o ${mechanism}_reg.out + $CHEMKIN_BIN/GetSolution + $CHEMKIN_BIN/CKSolnTranspose + cd .. + + + + diff --git a/bash_backup/beforeoptimize_dummy.pbs b/bash_backup/beforeoptimize_dummy.pbs new file mode 100644 index 0000000..28851ca --- /dev/null +++ b/bash_backup/beforeoptimize_dummy.pbs @@ -0,0 +1,58 @@ +#PBS -A jmartz_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=00:10:00,pmem=2500mb,procs=1 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 1-25 +#PBS -joe + +cd $PBS_O_WORKDIR + + +##setting +pressure="20" +mechanism="Ra_Reitz"; +fuel_name="n_heptane"; +class_name="class2_class4_class6_moderate"; +equi="1"; +matlabdate="01_30_2017" +## + + +## access and make file + + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + + +mkdir beforeoptimize_dummy +cd beforeoptimize_dummy +## + + + +##calculate chemkin + + echo " calcultating temp ${PBS_ARRAYID} " + mkdir ck_${PBS_ARRAYID} + cd ck_${PBS_ARRAYID} + + + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/result/${PBS_ARRAYID}/${mechanism}_${fuel_name}.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_${fuel_name}_${class_name}.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_therm.dat ./ + cp /scratch/engin_flux/unghee/chemkin/chemkindata.dtd ./ + + + $CHEMKIN_BIN/chem -i ${mechanism}_${fuel_name}_${class_name}.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/CKReactorGenericClosed -i ${mechanism}_${fuel_name}.inp -o ${mechanism}_reg.out + $CHEMKIN_BIN/GetSolution + $CHEMKIN_BIN/CKSolnTranspose + cd .. + + + + diff --git a/bash_backup/case_chemkinrun.pbs b/bash_backup/case_chemkinrun.pbs new file mode 100644 index 0000000..c2c17df --- /dev/null +++ b/bash_backup/case_chemkinrun.pbs @@ -0,0 +1,34 @@ +#PBS -A jmartz_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=6:20:00,pmem=2500mb,procs=1 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 1-7 +#PBS -joe + +#cd $PBS_O_WORKDIR +#echo $PBS_O_WORKDIR + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} +cd modify_sensitivity_${iteration_numb}_iteration +cd class${classnumb} +cd $casenumber + + + echo " calcultating v ${PBS_ARRAYID} " + mkdir ck_${PBS_ARRAYID} + cd ck_${PBS_ARRAYID} + + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/modify_sensitivity_${iteration_numb}_iteration/${mechanism}_${fuel_name}_class${classnumb}_v_${PBS_ARRAYID}.inp ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/result/$casenumber/${mechanism}_${fuel_name}.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_therm.dat ./ + cp /scratch/engin_flux/unghee/chemkin/chemkindata.dtd ./ + + $CHEMKIN_BIN/chem -i ${mechanism}_${fuel_name}_class${classnumb}_v_${PBS_ARRAYID}.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/CKReactorGenericClosed -i ${mechanism}_${fuel_name}.inp -o ${mechanism}_reg.out + $CHEMKIN_BIN/GetSolution + $CHEMKIN_BIN/CKSolnTranspose + cd .. \ No newline at end of file diff --git a/bash_backup/copy.sh b/bash_backup/copy.sh new file mode 100644 index 0000000..f6626eb --- /dev/null +++ b/bash_backup/copy.sh @@ -0,0 +1,34 @@ +#!/bin/bash +#PBS -A engin_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=00:01:00,pmem=2500mb,procs=1 +#PBS -l qos=flux +#PBS -q flux +#PBS -joe + +#cd $PBS_O_WORKDIR +#cd Ra_Reitz +# afterwards do it as a loop i=1:10 v{i}.. +#cd afteroptimize_test_40atm_phi1 +dir=$(pwd) +cd $dir +echo $dir + +for((i=1;i<=25;i++)) + do + mkdir id_$i + cd id_$i + + cp ../ck_$i/CKSoln_solution_point_value_vs_solution_number.csv ./ + cp ../ck_$i/*.inp ./ + cp ../ck_$i/*.out ./ + cp ../ck_$i/XMLdata.zip ./ + cp ../ck_$i/chem.asc ./ + cp ../ck_$i/chemkindata.dtd ./ + cd .. + + rm -rf ck_$i +done diff --git a/bash_backup/modify.pbs b/bash_backup/modify.pbs new file mode 100644 index 0000000..19cbab6 --- /dev/null +++ b/bash_backup/modify.pbs @@ -0,0 +1,86 @@ +#PBS -A engin_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=6:20:00,pmem=2500mb,procs=3 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 11,15,21,22,23,24,26,27,28 +#PBS -joe + +# 11,15,21,22,23,24,26,27,28 + +cd $PBS_O_WORKDIR +echo $PBS_O_WORKDIR +#DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" +#echo $DIR + +##setting +# check matlab setting as well!! +fuel_name="n_dodecane" +pressure="40" +equi="1" +matlabdate="03_23_2017_1_iteration" +#numbOfClass=2 +#classnumb="class6" +numbOfTemp=11; +mechanism="MFC"; +file_name="${mechanism}_base.inp"; +## + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + +mkdir modify +cd modify + +echo "running modifcation.m" +#cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_base.inp ./ +cp /scratch/engin_flux/unghee/chemkin/mechanisms/${file_name} ./ +#cp /scratch/engin_flux/unghee/chemkin/modification.m ./ +cp /scratch/engin_flux/unghee/chemkin/modification_func.m ./ + +echo $file_name + +#matlab -nodisplay -r "modification_func('${file_name}','${mechanism}','${fuel_name}')" +mkdir class${PBS_ARRAYID} +cd class${PBS_ARRAYID} + + + + +for j in $(seq 1 ${numbOfTemp}) +#for j in $(seq 1 1) +do + mkdir $j + + cd $j + + + + + + for i in {1..7} + do + echo " calcultating v $i " + mkdir ck_$i + cd ck_$i + +# cp ../$i/* ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/modify/${mechanism}_${fuel_name}_class${PBS_ARRAYID}_v_$i.inp ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/input/$j/${mechanism}_$fuel_name.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_therm.dat ./ + cp /scratch/engin_flux/unghee/chemkin/chemkindata.dtd ./ + + $CHEMKIN_BIN/chem -i ${mechanism}_${fuel_name}_class${PBS_ARRAYID}_v_$i.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/CKReactorGenericClosed -i ${mechanism}_${fuel_name}.inp -o ${mechanism}_reg.out + $CHEMKIN_BIN/GetSolution + $CHEMKIN_BIN/CKSolnTranspose + cd .. + + done + +cd .. + +done + diff --git a/bash_backup/modify_loop.pbs b/bash_backup/modify_loop.pbs new file mode 100644 index 0000000..c813958 --- /dev/null +++ b/bash_backup/modify_loop.pbs @@ -0,0 +1,84 @@ +#PBS -A jmartz_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=6:20:00,pmem=2500mb,procs=3 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 28 +#PBS -joe + + +cd $PBS_O_WORKDIR +echo $PBS_O_WORKDIR +#DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" +#echo $DIR + +##setting +# check matlab setting as well!! +fuel_name="n_dodecane" +pressure="20" +equi="1" +matlabdate="03_16_2017_2_iteration" +#numbOfClass=2 +#classnumb="class6" +class_name="class11_class15_class22_class24_class26_class27_class28"; +numbOfTemp=26; +mechanism="MFC"; +## + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + +mkdir modify +cd modify + + + + +mkdir class${PBS_ARRAYID} +cd class${PBS_ARRAYID} + + +echo "running modifcation.m" +#cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_base.inp ./ +cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/output/${mechanism}_${fuel_name}_${class_name}_afteroptimize.inp ./ +cp /scratch/engin_flux/unghee/chemkin/modification.m ./ +matlab -nodisplay -r modification + + +#for j in $(seq 1 ${numbOfTemp}) +for j in $(seq 26 26) +do + mkdir $j + + cd $j + + + + + + for i in {1..7} + do + echo " calcultating v $i " + mkdir ck_$i + cd ck_$i + +# cp ../$i/* ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/modify/class${PBS_ARRAYID}/${mechanism}_${fuel_name}_class${PBS_ARRAYID}_v_$i.inp ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/input/$j/${mechanism}_$fuel_name.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_therm.dat ./ + cp /scratch/engin_flux/unghee/chemkin/chemkindata.dtd ./ + + $CHEMKIN_BIN/chem -i ${mechanism}_${fuel_name}_class${PBS_ARRAYID}_v_$i.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/CKReactorGenericClosed -i ${mechanism}_${fuel_name}.inp -o ${mechanism}_reg.out + $CHEMKIN_BIN/GetSolution + $CHEMKIN_BIN/CKSolnTranspose + cd .. + + done + +cd .. + +done + diff --git a/bash_backup/modify_loop_2.pbs b/bash_backup/modify_loop_2.pbs new file mode 100644 index 0000000..c369a17 --- /dev/null +++ b/bash_backup/modify_loop_2.pbs @@ -0,0 +1,93 @@ +#PBS -A jmartz_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=6:20:00,pmem=2500mb,procs=4 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 1-7 +#PBS -joe + + +cd $PBS_O_WORKDIR +echo $PBS_O_WORKDIR +#DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" +#echo $DIR + +##setting +# check matlab setting as well!! +fuel_name="n_dodecane" +pressure="20" +equi="1" +matlabdate="03_16_2017_3_iteration" +class_name="class26_class27"; +mechanism="MFC"; +## + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + +mkdir modify +cd modify + +echo "running modifcation.m" +cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_base.inp ./ +cp /scratch/engin_flux/unghee/chemkin/modification.m ./ +cp /scratch/engin_flux/unghee/chemkin/modification.m ./ +matlab -nodisplay -r modification + + + +for i in 11 +do + +mkdir class$i +cd class$i + + +#echo "running modifcation.m" +#cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_base.inp ./ +#cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/output/${mechanism}_${fuel_name}_${class_name}_afteroptimize.inp ./ +#cp /scratch/engin_flux/unghee/chemkin/modification.m ./ +#matlab -nodisplay -r modification + + + for j in $(seq 1 ${numbOfTemp}) + do +#for j in $(seq 1 1) +#for j in 11 15 21 22 23 24 26 27 28 + + mkdir $j + + cd $j + + + + + +# for ${PBS_ARRAYID} in {1..7} +# do + echo " calcultating v ${PBS_ARRAYID} " + mkdir ck_${PBS_ARRAYID} + cd ck_${PBS_ARRAYID} + +# cp ../$i/* ./ + #cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/modify/class${PBS_ARRAYID}/${mechanism}_${fuel_name}_class${PBS_ARRAYID}_v_$i.inp ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/modify/${mechanism}_${fuel_name}_class$i_v_${PBS_ARRAYID}.inp ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/input/$j/${mechanism}_$fuel_name.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_therm.dat ./ + cp /scratch/engin_flux/unghee/chemkin/chemkindata.dtd ./ + + #$CHEMKIN_BIN/chem -i ${mechanism}_${fuel_name}_class${PBS_ARRAYID}_v_$i.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/chem -i ${mechanism}_${fuel_name}_class$i_v_${PBS_ARRAYID}.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/CKReactorGenericClosed -i ${mechanism}_${fuel_name}.inp -o ${mechanism}_reg.out + $CHEMKIN_BIN/GetSolution + $CHEMKIN_BIN/CKSolnTranspose + cd .. + +# done + done +cd .. + +done + diff --git a/bash_backup/modify_loop_matlab.pbs b/bash_backup/modify_loop_matlab.pbs new file mode 100644 index 0000000..19b39de --- /dev/null +++ b/bash_backup/modify_loop_matlab.pbs @@ -0,0 +1,51 @@ +#PBS -A jmartz_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=6:20:00,pmem=2500mb,procs=3 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 11,15,21,22,23,24,26,27,28 +#PBS -joe + + +cd $PBS_O_WORKDIR +echo $PBS_O_WORKDIR +#DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" +#echo $DIR + +##setting +# check matlab setting as well!! +fuel_name="n_dodecane" +pressure="20" +equi="1" +matlabdate="03_23_2017_1_iteration" +#numbOfClass=2 +#classnumb="class6" +#class_name="class26_class27"; +numbOfTemp=26; +mechanism="MFC"; +## + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + +mkdir modify +cd modify + +echo "running modifcation.m" +cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_base.inp ./ +#cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/output/${mechanism}_${fuel_name}_${class_name}_afteroptimize.inp ./ +cp /scratch/engin_flux/unghee/chemkin/modification.m ./ +matlab -nodisplay -r modification + + + + +mkdir class${PBS_ARRAYID} +cd class${PBS_ARRAYID} + + + + + diff --git a/bash_backup/modify_sensitivity.pbs b/bash_backup/modify_sensitivity.pbs new file mode 100644 index 0000000..5facef0 --- /dev/null +++ b/bash_backup/modify_sensitivity.pbs @@ -0,0 +1,80 @@ +#PBS -A jmartz_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=20:20:00,pmem=2500mb,procs=3 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 11,15,21,22,23,24,26,27,28 +#PBS -joe + + +cd $PBS_O_WORKDIR +echo $PBS_O_WORKDIR +#DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" +#echo $DIR + +##setting +# check matlab setting as well!! +fuel_name="n_dodecane" +pressure="20" +equi="1" +matlabdate="03_09_2017_3_iteration" +numbOfTemp=25; +mechanism="MFC"; +## + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + +mkdir modify_sensitivity +cd modify_sensitivity + + + + +mkdir class${PBS_ARRAYID} +cd class${PBS_ARRAYID} + + +echo "running modifcation.m" +cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_base.inp ./ +cp /scratch/engin_flux/unghee/chemkin/modification.m ./ +matlab -nodisplay -r modification + + +for j in $(seq 1 ${numbOfTemp}) +#for j in $(seq 10 10) +do + mkdir $j + + cd $j + + + + + + for i in {1..7} + do + echo " calcultating v $i " + mkdir ck_$i + cd ck_$i + +# cp ../$i/* ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/modify_sensitivity/class${PBS_ARRAYID}/${mechanism}_${fuel_name}_class${PBS_ARRAYID}_v_$i.inp ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/result/$j/${mechanism}_$fuel_name.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_therm.dat ./ + cp /scratch/engin_flux/unghee/chemkin/chemkindata.dtd ./ + + $CHEMKIN_BIN/chem -i ${mechanism}_${fuel_name}_class${PBS_ARRAYID}_v_$i.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/CKReactorGenericClosed -i ${mechanism}_${fuel_name}.inp -o ${mechanism}_reg.out + $CHEMKIN_BIN/GetSolution + $CHEMKIN_BIN/CKSolnTranspose + cd .. + + done + +cd .. + +done + diff --git a/bash_backup/modify_sensitivity_2.pbs b/bash_backup/modify_sensitivity_2.pbs new file mode 100644 index 0000000..4f660d1 --- /dev/null +++ b/bash_backup/modify_sensitivity_2.pbs @@ -0,0 +1,77 @@ +#PBS -A jmartz_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=6:20:00,pmem=2500mb,procs=1 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 1-7 +#PBS -joe + + +cd $PBS_O_WORKDIR +echo $PBS_O_WORKDIR +#DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" +#echo $DIR + +##setting +# check matlab setting as well!! +fuel_name="n_dodecane" +pressure="20" +equi="1" +matlabdate="03_23_2017_1_iteration" +numbOfTemp=25; +mechanism="MFC"; +## + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + +mkdir modify_sensitivity_2 +cd modify_sensitivity_2 + +echo "running modifcation.m" +cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_base.inp ./ +cp /scratch/engin_flux/unghee/chemkin/modification.m ./ +matlab -nodisplay -r modification + + + +for i in 23 24 +do + +mkdir class$i +cd class$i + + + + + for j in $(seq 1 ${numbOfTemp}) + do + + mkdir $j + + cd $j + + echo " calcultating v ${PBS_ARRAYID} " + mkdir ck_${PBS_ARRAYID} + cd ck_${PBS_ARRAYID} + + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/modify_sensitivity_2/${mechanism}_${fuel_name}_class${i}_v_${PBS_ARRAYID}.inp ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/result/$j/${mechanism}_${fuel_name}.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_therm.dat ./ + cp /scratch/engin_flux/unghee/chemkin/chemkindata.dtd ./ + + $CHEMKIN_BIN/chem -i ${mechanism}_${fuel_name}_class${i}_v_${PBS_ARRAYID}.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/CKReactorGenericClosed -i ${mechanism}_${fuel_name}.inp -o ${mechanism}_reg.out + $CHEMKIN_BIN/GetSolution + $CHEMKIN_BIN/CKSolnTranspose + cd .. + + cd .. + done + +cd .. + +done + diff --git a/bash_backup/modify_sensitivity_case.pbs b/bash_backup/modify_sensitivity_case.pbs new file mode 100644 index 0000000..7c1e5af --- /dev/null +++ b/bash_backup/modify_sensitivity_case.pbs @@ -0,0 +1,88 @@ +#PBS -A jmartz_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=8:20:00,pmem=2500mb,procs=1 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 1-25 +#PBS -joe + + +cd $PBS_O_WORKDIR +echo $PBS_O_WORKDIR + +##setting######################### + +fuel_name="n_dodecane" +pressure="20" +equi="1" +matlabdate="04_12_2017" +mechanism="MFC"; +iteration_numb="3"; +#class_name="class11_class15_class22_class24_class26_class27_class28"; +#file_name="${mechanism}_${fuel_name}_${class_name}_afteroptimize.inp"; +#file_name="${mechanism}_base.inp"; +file_name="${mechanism}_afteroptimize_${iteration_numb}_iteration.inp"; +#################################### + +##count numb of iteration############## + +#cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} +#cd output/ +#numfiles=(*) +#numfiles=${#numfiles[@]} +#echo "number of iteration" +#numfiles=$(expr $numfiles) +#echo $numfiles + +##################################### + + + + + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + +mkdir modify_sensitivity_${iteration_numb}_iteration +cd modify_sensitivity_${iteration_numb}_iteration + +#echo "running modifcation.m" +#cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_base.inp ./ +cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/output/${file_name} ./ + + +cp /scratch/engin_flux/unghee/chemkin/modification_func.m ./ +echo $file_name +matlab -nodisplay -r "modification_func('${file_name}','${mechanism}','${fuel_name}')" + + +for i in 11 15 21 22 23 24 26 27 28 +#for i in 28 +do + +mkdir class$i +cd class$i + + + mkdir ${PBS_ARRAYID} + + cd ${PBS_ARRAYID} + cp /scratch/engin_flux/unghee/chemkin/case_chemkinrun.pbs ./ + qsub -v casenumber=${PBS_ARRAYID},fuel_name=${fuel_name},pressure=${pressure},equi=${equi},matlabdate=${matlabdate},mechanism=${mechanism},classnumb=${i},iteration_numb=${iteration_numb} case_chemkinrun.pbs + cd .. + +cd .. + +done + + + +########### + +#matlab -nodisplay -r "modification_afteroptimize_github_func('${file_name}','${mechanism}','${fuel_name}','${numfiles}')" + + + + diff --git a/bash_backup/modify_sensitivity_case_past.pbs b/bash_backup/modify_sensitivity_case_past.pbs new file mode 100644 index 0000000..581d80c --- /dev/null +++ b/bash_backup/modify_sensitivity_case_past.pbs @@ -0,0 +1,74 @@ +#PBS -A engin_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=6:20:00,pmem=2500mb,procs=1 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 1-25 +#PBS -joe + + +cd $PBS_O_WORKDIR +echo $PBS_O_WORKDIR +#DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" +#echo $DIR + +##setting +# check matlab setting as well!! +fuel_name="n_dodecane" +pressure="20" +equi="1" +matlabdate="03_23_2017_2_iteration" +#numbOfTemp=25; +mechanism="MFC"; +class_name="class11_class15_class22_class24_class26_class27_class28"; +file_name="${mechanism}_${fuel_name}_${class_name}_afteroptimize.inp"; +#file_name="${mechanism}_base.inp"; +## + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + +mkdir modify_sensitivity_4 +cd modify_sensitivity_4 + +#echo "running modifcation.m" +#cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_base.inp ./ +cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/output/${file_name} ./ +#cp /scratch/engin_flux/unghee/chemkin/modification.m ./ +#matlab -nodisplay -r modification +cp /scratch/engin_flux/unghee/chemkin/modification_func.m ./ +echo $file_name +#matlab -nodisplay -r "modification_func('${file_name}','${mechanism}','${fuel_name}')" + + +for i in 11 15 21 22 23 24 26 27 28 +#for i in 11 +do + +mkdir class$i +cd class$i + + + + +# for j in $(seq 1 ${numbOfTemp}) +# for j in $(seq 1 1) +# do + +# mkdir $j + mkdir ${PBS_ARRAYID} +# cd $j + cd ${PBS_ARRAYID} + cp /scratch/engin_flux/unghee/chemkin/case_chemkinrun.pbs ./ +# qsub -v casenumber=$j fuel_name="${fuel_name}" pressure=$pressure equi=$equi matlabdate=$matlabdate mechanism=$mechanism case_chemkinrun.pbs +# qsub -v casenumber=$j,fuel_name=\"${fuel_name}\" case_chemkinrun.pbs + qsub -v casenumber=${PBS_ARRAYID},fuel_name=${fuel_name},pressure=${pressure},equi=${equi},matlabdate=${matlabdate},mechanism=${mechanism},classnumb=${i} case_chemkinrun.pbs + cd .. +# done + +cd .. + +done + diff --git a/bash_backup/modify_sensitivity_loop.pbs b/bash_backup/modify_sensitivity_loop.pbs new file mode 100644 index 0000000..0ed797f --- /dev/null +++ b/bash_backup/modify_sensitivity_loop.pbs @@ -0,0 +1,82 @@ +#PBS -A jmartz_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=20:20:00,pmem=2500mb,procs=3 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 11,15,21,22,23,24,26,27,28 +#PBS -joe + + +cd $PBS_O_WORKDIR +echo $PBS_O_WORKDIR +#DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" +#echo $DIR + +##setting +# check matlab setting as well!! +fuel_name="n_dodecane" +pressure="20" +equi="1" +matlabdate="03_16_2017_5_iteration" +numbOfTemp=25; +mechanism="MFC"; +class_name="class22_class24_class26_class27" +## + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + +mkdir modify_sensitivity +cd modify_sensitivity + + + + +mkdir class${PBS_ARRAYID} +cd class${PBS_ARRAYID} + + +echo "running modifcation.m" +cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/output/${mechanism}_${fuel_name}_${class_name}_afteroptimize.inp ./ +cp /scratch/engin_flux/unghee/chemkin/modification.m ./ +matlab -nodisplay -r modification + + +for j in $(seq 1 ${numbOfTemp}) +#for j in $(seq 21 21) +do + mkdir $j + + cd $j + + + + + + for i in {1..7} +# for i in 2 + do + echo " calcultating v $i " + mkdir ck_$i + cd ck_$i + +# cp ../$i/* ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/modify_sensitivity/class${PBS_ARRAYID}/${mechanism}_${fuel_name}_class${PBS_ARRAYID}_v_$i.inp ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/result/$j/${mechanism}_$fuel_name.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_therm.dat ./ + cp /scratch/engin_flux/unghee/chemkin/chemkindata.dtd ./ + + $CHEMKIN_BIN/chem -i ${mechanism}_${fuel_name}_class${PBS_ARRAYID}_v_$i.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/CKReactorGenericClosed -i ${mechanism}_${fuel_name}.inp -o ${mechanism}_reg.out + $CHEMKIN_BIN/GetSolution + $CHEMKIN_BIN/CKSolnTranspose + cd .. + + done + +cd .. + +done + diff --git a/bash_backup/modify_sensitivity_loop_2.pbs b/bash_backup/modify_sensitivity_loop_2.pbs new file mode 100644 index 0000000..9a72009 --- /dev/null +++ b/bash_backup/modify_sensitivity_loop_2.pbs @@ -0,0 +1,79 @@ +#PBS -A jmartz_flux +#PBS -M unghee@umich.edu +#PBS -m abe +#PBS -V +#PBS -N Ra_Reitz +#PBS -l walltime=6:20:00,pmem=2500mb,procs=1 +#PBS -l qos=flux +#PBS -q flux +#PBS -t 1-7 +#PBS -joe + + +cd $PBS_O_WORKDIR +echo $PBS_O_WORKDIR +#DIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" && pwd )" +#echo $DIR + +##setting +# check matlab setting as well!! +fuel_name="n_dodecane" +pressure="20" +equi="1" +matlabdate="03_16_2017_5_iteration" +class_name="class22_class24_class26_class27"; +numbOfTemp=25; +mechanism="MFC"; +## + +cd /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate} + +mkdir modify_sensitivity_2 +cd modify_sensitivity_2 + +echo "running modifcation.m" +cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_base.inp ./ +cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/output/${mechanism}_${fuel_name}_${class_name}_afteroptimize.inp ./ +cp /scratch/engin_flux/unghee/chemkin/modification.m ./ +matlab -nodisplay -r modification + + + +for i in 23 24 +do + +mkdir class$i +cd class$i + + + + + for j in $(seq 1 ${numbOfTemp}) + do + + mkdir $j + + cd $j + + echo " calcultating v ${PBS_ARRAYID} " + mkdir ck_${PBS_ARRAYID} + cd ck_${PBS_ARRAYID} + + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/modify_sensitivity_2/${mechanism}_${fuel_name}_class${i}_v_${PBS_ARRAYID}.inp ./ + cp /scratch/engin_flux/unghee/chemkin/${mechanism}/${fuel_name}_${pressure}atm_phi${equi}_${matlabdate}/result/$j/${mechanism}_${fuel_name}.inp ./ + cp /scratch/engin_flux/unghee/chemkin/mechanisms/${mechanism}_therm.dat ./ + cp /scratch/engin_flux/unghee/chemkin/chemkindata.dtd ./ + + $CHEMKIN_BIN/chem -i ${mechanism}_${fuel_name}_class${i}_v_${PBS_ARRAYID}.inp -o ${mechanism}_chem.out -d ${mechanism}_therm.dat + $CHEMKIN_BIN/CKReactorGenericClosed -i ${mechanism}_${fuel_name}.inp -o ${mechanism}_reg.out + $CHEMKIN_BIN/GetSolution + $CHEMKIN_BIN/CKSolnTranspose + cd .. + + cd .. + done + +cd .. + +done + diff --git a/chemkin_input_writer_Ra_Reitz.m b/chemkin_input_writer_Ra_Reitz.m new file mode 100644 index 0000000..440c119 --- /dev/null +++ b/chemkin_input_writer_Ra_Reitz.m @@ -0,0 +1,135 @@ +%% +clear all +close all + +addpath(pwd) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% change + +mechanism='MFC'; %mechanism + +fuel_name='JetA'; %surrogate fuel name + + +pressure_list=[40]; +for k=1:length(pressure_list) + pressure_text{k}=[num2str(pressure_list(k)),'atm']; + +end + +temp=700:25:1300; %K +% temp=975; +equi=1; %equivalence ratio +date = '05_03_2017'; + + + +% change +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +file_name=[mechanism,'_',fuel_name,'.inp']; + + +for k = 1:length(pressure_list) + + pressure = pressure_list(k); +directory=['C:\Users\unghee\Dropbox\post_process','\',mechanism,'\',fuel_name,'_',pressure_text{k},'_','phi',num2str(equi),'_',date]; + +mkdir(directory) +cd(directory) + + + +mkdir('result') +cd('result') + +%% problem type definition + +problem_type='CONV'; % problem type, CONV: constant volume, solve energy equation +energy_eq='ENRG'; % ENRG, solve gas energy equation +sol='TRAN'; % Transient solver + +%% physical property +% pressure=20; %atm +% temp=700:25:1300; %K +% +% equi=1.5; %equivalence ratio + +ifac=0.1; %ignition noise filtering factor +% klim='oh'; %species maximum fraction +if strcmp(mechanism,'MFC') + klim='OH'; %species maximum fraction +elseif strcmp(mechanism,'Ra_Reitz') + klim='oh'; + +end +dtig=400; %K, ignition delay definition - the time that the ambient I increases dtig +qlos=0; %heat loss + +%% species property +if strcmp(mechanism,'MFC') +comb_product={'CO2' 'H2O' 'N2'}; +elseif strcmp(mechanism,'Ra_Reitz') +comb_product={'co2' 'h2o' 'n2'}; +end + +load_surrogate_composition +% +% %UM1 +% fuel={'nc12h26' 'hmn' 'mch' 'c6h5ch3'}; +% fuel_composition=[0.3844 0.1484 0.2336 0.2336]; %mole fraction + +% %UM2 +% fuel={'nc12h26' 'hmn' 'decalin' 'c6h5ch3'}; +% fuel_composition=[0.2897 0.1424 0.3188 0.2491]; %mole fraction + + +% fuel={'nc12h26'}; +% fuel_composition=[1]; %mole fraction + +if strcmp(fuel_name,'n_heptane') + fuel={'nc7h16'}; + fuel_composition=[1]; %mole fraction + +elseif strcmp(fuel_name,'n_dodecane') + fuel={'NC12H26'}; + fuel_composition=[1]; %mole fraction + +elseif strcmp(fuel_name,'JetA') +fuel={'NC12H26' 'HMN' 'DECALIN' 'C6H5CH3'}; +fuel_composition=[0.4706 0.1669 0.2419 0.1206]; %mole fraction +end + + +if strcmp(mechanism,'MFC') +oxidizer={'N2' 'O2'}; +elseif strcmp(mechanism,'Ra_Reitz') +oxidizer={'n2' 'o2'}; +end +oxidizer_composition=[0.79 0.21]; %mole fraction + +%% solver control +time=1; %sec + +%% + + +% mkdir(directory); +% cd(directory); +for i=1:length(temp) + dir=num2str(i); + mkdir(dir); + cd(dir); + temperature=temp(i); + write_input_file; + cd .. +end +cd .. +cd .. + + + +end + diff --git a/find_rate_weighting.m b/find_rate_weighting.m new file mode 100644 index 0000000..b1f6bf1 --- /dev/null +++ b/find_rate_weighting.m @@ -0,0 +1,48 @@ +function differTotal = find_rate_weighting(X,coefs,Temp,numbOfClass,classnumb_text,Target_data,numbOftarget1,numbOftarget2) +%% This function automatically writes the quadratic function for defining objective function i.e. obj = (target-f(A,E)) + + +differTotal= 0; +timeTotal =zeros(length(Temp),1) ; +time = zeros(length(Temp),numbOfClass); +time_weight_differ = zeros(length(Temp),numbOfClass); +for j= 1: length(Temp) + + if j<=numbOftarget1 % normalizing weighting + W = 1/numbOftarget1; + else + W = 1/numbOftarget2; + end + + + for i = 1: numbOfClass + + A=X(2*i-1); + E=X(2*i); +% W=X(3*i); % setting weight as variable +% W=1; % equal weighting + +% automatically writing the second order polynomial func. + time(j,i)=coefs.(classnumb_text{i})(j,1)*log(A)+coefs.(classnumb_text{i})(j,2)*log(Temp(j))... + +coefs.(classnumb_text{i})(j,3)*E+coefs.(classnumb_text{i})(j,4)*log(A)*E... + +coefs.(classnumb_text{i})(j,5)*(log(A))^2+coefs.(classnumb_text{i})(j,6)*(E)^2; + + time_weight_differ(j,i)= W*abs(log10(Target_data(j))-log10(exp(time(j,i)))); + + timeTotal(j)=time_weight_differ(j,i)+timeTotal(j); % sum up class + + + end + + differTotal = timeTotal(j) + differTotal; % sum up Temp + + +end + + + + +end + + + diff --git a/flux_version/find_rate_weighting.m b/flux_version/find_rate_weighting.m new file mode 100644 index 0000000..172d434 --- /dev/null +++ b/flux_version/find_rate_weighting.m @@ -0,0 +1,49 @@ +function differTotal = find_rate_weighting(X,coefs,Temp,numbOfClass,classnumb_text,Target_data,numbOftarget1,numbOftarget2) +differTotal= 0; +timeTotal =zeros(length(Temp),1) ; +time = zeros(length(Temp),numbOfClass); +time_weight_differ = zeros(length(Temp),numbOfClass); +for j= 1: length(Temp) + + if j<=numbOftarget1 + W = 1/numbOftarget1; + else + W = 1/numbOftarget2; + end + + + for i = 1: numbOfClass + + A=X(2*i-1); + E=X(2*i); +% W=X(3*i); +% W=1; +% time(j,i)=coefs.(classnumb_text{i})(j,1)*log(X(2*i-1))+coefs.(classnumb_text{i})(j,2)*log(Temp(j))... +% +coefs.(classnumb_text{i})(j,3)*X(2*i)+coefs.(classnumb_text{i})(j,4)*log(X(2*i-1))*X(2*i)... +% +coefs.(classnumb_text{i})(j,5)*(log(X(2*i-1)))^2+coefs.(classnumb_text{i})(j,6)*(X(2*i))^2; + + time(j,i)=coefs.(classnumb_text{i})(j,1)*log(A)+coefs.(classnumb_text{i})(j,2)*log(Temp(j))... + +coefs.(classnumb_text{i})(j,3)*E+coefs.(classnumb_text{i})(j,4)*log(A)*E... + +coefs.(classnumb_text{i})(j,5)*(log(A))^2+coefs.(classnumb_text{i})(j,6)*(E)^2; + + time_weight_differ(j,i)= W*abs(log10(Target_data(j))-log10(exp(time(j,i)))); + + timeTotal(j)=time_weight_differ(j,i)+timeTotal(j); % sum up class +% timeTotal(j)= time(j,i) +timeTotal(j); +% X(3*i) + + end +% differ_tempj =W*abs(log10(Target_data(j))-log10(exp(timeTotal(j)))); + + differTotal = timeTotal(j) + differTotal; % sum up Temp + +% differTotal = differ_tempj + differTotal; % sum up Temp +end + + + + +end + + + diff --git a/flux_version/modification_afteroptimize_github.m b/flux_version/modification_afteroptimize_github.m new file mode 100644 index 0000000..cca1453 --- /dev/null +++ b/flux_version/modification_afteroptimize_github.m @@ -0,0 +1,221 @@ + +clear all; +%% Import data from text file. +% Script for importing data from the following text file: +% +% /home/jordan/Documents/Research/matlabfolder/mech.dat +% +% To extend the code to different selected data or a different text file, +% generate a function instead of a script. + +% Auto-generated by MATLAB on 2016/08/24 17:10:49 + +%% Initialize variables. +% filename = '/scratch/engin_flux/unghee/chemkin/mechanisms/mech_ERC-MultiChem+Bio_Brakora2012.inp'; +currentFolder = pwd; +cd(pwd); +% cd mechanisms/; +% filename = 'mech_ERC-MultiChem+Bio_Brakora2012_base.inp'; +filename = 'MFC_base.inp'; +delimiter = ' '; + +%% Read columns of data as strings: +% For more information, see the TEXTSCAN documentation. +% formatSpec = '%s%s%s%s%s%s%s%s%s%[^\n\r]'; +formatSpec = '%s%s%s%s%s%s%s%s%s%s%s%[^\n\r]'; +%% Open the text file. +fileID = fopen(filename,'r'); + +%% Read columns of data according to format string. +% This call is based on the structure of the file used to generate this +% code. If an error occurs for a different file, try regenerating the code +% from the Import Tool. +dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'ReturnOnError', false); + +%% Close the text file. +fclose(fileID); + +%% Convert the contents of columns containing numeric strings to numbers. +% Replace non-numeric strings with NaN. +raw = repmat({''},length(dataArray{1}),length(dataArray)-1); +for col=1:length(dataArray)-1 + raw(1:length(dataArray{col}),col) = dataArray{col}; +end +numericData = NaN(size(dataArray{1},1),size(dataArray,2)); + + +%% Split data into numeric and cell columns. +rawNumericColumns = {}; +rawCellColumns = raw(:, [1,2,3,4,5,6,7,8,9]); + + + %% Replace non-numeric cells with NaN + +[A,B]= find(ismember(rawCellColumns,'ELEMENTS')); + +locationClass = [A B]; + +[rowCellNumber,] = size(rawCellColumns); + +rawCellColumns2=rawCellColumns; + +ModStart = 0; + +mechanism='MFC'; +fuel_name='n_dodecane'; +date = '03_23_2017_1_iteration'; +iteration_numb = 1; +class_numb=[11 15 22 24 26 27 28]; +% class_numb=classnumb; +numbOfClass = length(class_numb); +class_numb_text = {}; +for k=1:numbOfClass +% classnumb_text{classnumb(k)}=['class',num2str(classnumb(k))]; + class_numb_text=[class_numb_text ['class',num2str(class_numb(k))]]; +end +pressure=[20 40]; +equi=1; +for k=1:length(pressure) + pressure_text{k}=[num2str(pressure(k)),'atm']; + +end +load final_result.mat + + +for i=1 : rowCellNumber + [ ~, columnCellNumber] = size(rawCellColumns(i,:)); +% for j= 1 : cellfun('length',rawCellColumns{i,:}) + for j= 1 : 1: columnCellNumber + + if strmatch('reactions',rawCellColumns{i,j})==1; + ModStart = 1; + end + + + for k = 1: numbOfClass +% if class_numb(k) == 1 & (strfind(rawCellColumns{i,j},'nc7h16+h=') == 1)... % n heptane erc +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1)... +% |(class_numb(k) == 2 & (strfind(rawCellColumns{i,j},'nc7h16+oh=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 3 & (strfind(rawCellColumns{i,j},'nc7h16+ho2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 4 & (strfind(rawCellColumns{i,j},'nc7h16+o2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 5 & (strfind(rawCellColumns{i,j},'c7h15-2+o2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% | class_numb(k) == 6 & (strfind(rawCellColumns{i,j},'c7h15o2+o2=') == 1)... % n heptane erc +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1)... +% |(class_numb(k) == 7 & (strfind(rawCellColumns{i,j},'c7ket12=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 8 & (strfind(rawCellColumns{i,j},'c5h11co=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 9 & (strfind(rawCellColumns{i,j},'c7h15-2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... + if class_numb(k) == 27 & (strfind(rawCellColumns{i,j},'C12OOH') ~= 0)... % n dodecane ske_361 + & (isempty(strfind(rawCellColumns{i,j},'O2=C12KET')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 21 & (strfind(rawCellColumns{i,j},'NC12H26+H=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+H=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 22 & (strfind(rawCellColumns{i,j},'NC12H26+OH=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+OH=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 23 & (strfind(rawCellColumns{i,j},'NC12H26+O=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+O=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 24 & (strfind(rawCellColumns{i,j},'NC12H26+HO2=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+HO2=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 11 & (strfind(rawCellColumns{i,j},'C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'+O2=C12H25O2')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 28 & (strfind(rawCellColumns{i,j},'C12KET') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'=OH+')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 26 & (strfind(rawCellColumns{i,j},'C12OOH') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'O2=C12OOH')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 15 & (strfind(rawCellColumns{i,j},'C12H25O2') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'=C12OOH')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i,j},'C12OOH6-3')) == 1)... + & (isempty(strfind(rawCellColumns{i,j},'C12OOH6-9')) == 1)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + + rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb(k)),'***','after Optimize']; + + rawCellColumns2{i,2} = final_result.(class_numb_text{k})(2,1); % Optimized A + rawCellColumns2{i,4} = final_result.(class_numb_text{k})(2,2); % Optimized E + + + end + end + end +end + +%% +% datename = '01_19_2017' +% save_filename=[mechanism,'_','output','_',datename]; + + +save_filename='output'; + +for k = 1 :length(pressure) +directory=['C:\Users\unghee\Dropbox\post_process','\',mechanism,'\',fuel_name,'_',pressure_text{k},'_','phi',num2str(equi),'_',date]; +% location=[currentloc,'\',mechanism,'\',directory]; +cd(directory) + +mkdir(save_filename); +cd(save_filename); + +class_numb_text = []; +for k=1:numbOfClass +% classnumb_text{classnumb(k)}=['class',num2str(classnumb(k))]; + class_numb_text=[class_numb_text '_' ['class',num2str(class_numb(k))]]; +end +% file_name=[mechanism,'_',fuel_name,class_numb_text,'_','afteroptimize','.inp']; +file_name=[mechanism,'_',num2str(iteration_numb),'_'class_numb_text,'_','afteroptimize','.inp']; +fileID = fopen(file_name,'w'); +output = cell(size(rawCellColumns2,1),size(rawCellColumns2,2)); + + + + for i = 1:size(rawCellColumns2,1) + for j = 1:size(rawCellColumns2,2) + if numel(rawCellColumns2{i,j}) == 0 + output{i,j} = ''; + % Check whether the content of cell i,j is + % numeric and convert numbers to strings. + elseif isnumeric(rawCellColumns2{i,j}) || islogical(rawCellColumns2{i,j}) + output{i,j} = num2str(rawCellColumns2{i,j}(1,1)); + + % If the cell already contains a string, nothing has to be done. + elseif ischar(rawCellColumns2{i,j}) + output{i,j} = rawCellColumns2{i,j}; + end; + + % Cell i,j is written to the output file. A delimiter is appended for + % all but the last element of each row. At the end of a row, a newline + % is written to the output file. + if j < size(rawCellColumns2,2) + fprintf(fileID,['%s',delimiter],output{i,j}); + else + fprintf(fileID,'%s\r\n',output{i,j}); + end + end; + end; + + +fclose(fileID); + +end +% exit +%type mech_ERC-MultiChem+Bio_Brakora2012_v1_test.inp + +%% Clear temporary variables +% clearvars filename delimiter formatSpec fileID dataArray ans raw col numericData rawCellColumns rawNumericColumns R; \ No newline at end of file diff --git a/flux_version/modification_afteroptimize_github_func.m b/flux_version/modification_afteroptimize_github_func.m new file mode 100644 index 0000000..c3f0329 --- /dev/null +++ b/flux_version/modification_afteroptimize_github_func.m @@ -0,0 +1,223 @@ +function modification_afteroptimize_github_func(filename,mechanism,fuel_name,iteration_numb) +%% Import data from text file. +% Script for importing data from the following text file: +% +% /home/jordan/Documents/Research/matlabfolder/mech.dat +% +% To extend the code to different selected data or a different text file, +% generate a function instead of a script. + +% Auto-generated by MATLAB on 2016/08/24 17:10:49 + +%% Initialize variables. +% filename = '/scratch/engin_flux/unghee/chemkin/mechanisms/mech_ERC-MultiChem+Bio_Brakora2012.inp'; +currentFolder = pwd; +cd(pwd); +% cd mechanisms/; +% filename = 'mech_ERC-MultiChem+Bio_Brakora2012_base.inp'; +% filename = 'MFC_base.inp'; +delimiter = ' '; + +%% Read columns of data as strings: +% For more information, see the TEXTSCAN documentation. +% formatSpec = '%s%s%s%s%s%s%s%s%s%[^\n\r]'; +formatSpec = '%s%s%s%s%s%s%s%s%s%s%s%[^\n\r]'; +%% Open the text file. +fileID = fopen(filename,'r'); + +%% Read columns of data according to format string. +% This call is based on the structure of the file used to generate this +% code. If an error occurs for a different file, try regenerating the code +% from the Import Tool. +dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'ReturnOnError', false); + +%% Close the text file. +fclose(fileID); + +%% Convert the contents of columns containing numeric strings to numbers. +% Replace non-numeric strings with NaN. +raw = repmat({''},length(dataArray{1}),length(dataArray)-1); +for col=1:length(dataArray)-1 + raw(1:length(dataArray{col}),col) = dataArray{col}; +end +numericData = NaN(size(dataArray{1},1),size(dataArray,2)); + + +%% Split data into numeric and cell columns. +rawNumericColumns = {}; +rawCellColumns = raw(:, [1,2,3,4,5,6,7,8,9]); + + + %% Replace non-numeric cells with NaN + +[A,B]= find(ismember(rawCellColumns,'ELEMENTS')); + +locationClass = [A B]; + +[rowCellNumber,] = size(rawCellColumns); + +rawCellColumns2=rawCellColumns; + +ModStart = 0; + +% mechanism='MFC'; +% fuel_name='n_dodecane'; +date = '03_23_2017_1_iteration'; +% iteration_numb = 1; +class_numb=[11 15 22 24 26 27 28]; +% class_numb=classnumb; +numbOfClass = length(class_numb); +class_numb_text = {}; +for k=1:numbOfClass +% classnumb_text{classnumb(k)}=['class',num2str(classnumb(k))]; + class_numb_text=[class_numb_text ['class',num2str(class_numb(k))]]; +end +pressure=[20 40]; +equi=1; +for k=1:length(pressure) + pressure_text{k}=[num2str(pressure(k)),'atm']; + +end +load final_result.mat + + +for i=1 : rowCellNumber + [ ~, columnCellNumber] = size(rawCellColumns(i,:)); +% for j= 1 : cellfun('length',rawCellColumns{i,:}) + for j= 1 : 1: columnCellNumber + + if strmatch('reactions',rawCellColumns{i,j})==1; + ModStart = 1; + end + + + for k = 1: numbOfClass +% if class_numb(k) == 1 & (strfind(rawCellColumns{i,j},'nc7h16+h=') == 1)... % n heptane erc +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1)... +% |(class_numb(k) == 2 & (strfind(rawCellColumns{i,j},'nc7h16+oh=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 3 & (strfind(rawCellColumns{i,j},'nc7h16+ho2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 4 & (strfind(rawCellColumns{i,j},'nc7h16+o2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 5 & (strfind(rawCellColumns{i,j},'c7h15-2+o2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% | class_numb(k) == 6 & (strfind(rawCellColumns{i,j},'c7h15o2+o2=') == 1)... % n heptane erc +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1)... +% |(class_numb(k) == 7 & (strfind(rawCellColumns{i,j},'c7ket12=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 8 & (strfind(rawCellColumns{i,j},'c5h11co=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 9 & (strfind(rawCellColumns{i,j},'c7h15-2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... + if class_numb(k) == 27 & (strfind(rawCellColumns{i,j},'C12OOH') ~= 0)... % n dodecane ske_361 + & (isempty(strfind(rawCellColumns{i,j},'O2=C12KET')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 21 & (strfind(rawCellColumns{i,j},'NC12H26+H=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+H=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 22 & (strfind(rawCellColumns{i,j},'NC12H26+OH=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+OH=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 23 & (strfind(rawCellColumns{i,j},'NC12H26+O=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+O=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 24 & (strfind(rawCellColumns{i,j},'NC12H26+HO2=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+HO2=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 11 & (strfind(rawCellColumns{i,j},'C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'+O2=C12H25O2')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 28 & (strfind(rawCellColumns{i,j},'C12KET') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'=OH+')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 26 & (strfind(rawCellColumns{i,j},'C12OOH') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'O2=C12OOH')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 15 & (strfind(rawCellColumns{i,j},'C12H25O2') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'=C12OOH')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i,j},'C12OOH6-3')) == 1)... + & (isempty(strfind(rawCellColumns{i,j},'C12OOH6-9')) == 1)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + + rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb(k)),'***','after Optimize']; + + rawCellColumns2{i,2} = final_result.(class_numb_text{k})(2,1); % Optimized A + rawCellColumns2{i,4} = final_result.(class_numb_text{k})(2,2); % Optimized E + + + end + end + end +end + +%% +% datename = '01_19_2017' +% save_filename=[mechanism,'_','output','_',datename]; + + +save_filename='output'; + +for k = 1 :length(pressure) +directory=['C:\Users\unghee\Dropbox\post_process','\',mechanism,'\',fuel_name,'_',pressure_text{k},'_','phi',num2str(equi),'_',date]; +% location=[currentloc,'\',mechanism,'\',directory]; +cd(directory) + +mkdir(save_filename); +cd(save_filename); + +class_numb_text = []; +for k=1:numbOfClass +% classnumb_text{classnumb(k)}=['class',num2str(classnumb(k))]; + class_numb_text=[class_numb_text '_' ['class',num2str(class_numb(k))]]; +end +% file_name=[mechanism,'_',fuel_name,class_numb_text,'_','afteroptimize','.inp']; +file_name=[mechanism,'_',num2str(iteration_numb),'_'class_numb_text,'_','afteroptimize','.inp']; +fileID = fopen(file_name,'w'); +output = cell(size(rawCellColumns2,1),size(rawCellColumns2,2)); + + + + for i = 1:size(rawCellColumns2,1) + for j = 1:size(rawCellColumns2,2) + if numel(rawCellColumns2{i,j}) == 0 + output{i,j} = ''; + % Check whether the content of cell i,j is + % numeric and convert numbers to strings. + elseif isnumeric(rawCellColumns2{i,j}) || islogical(rawCellColumns2{i,j}) + output{i,j} = num2str(rawCellColumns2{i,j}(1,1)); + + % If the cell already contains a string, nothing has to be done. + elseif ischar(rawCellColumns2{i,j}) + output{i,j} = rawCellColumns2{i,j}; + end; + + % Cell i,j is written to the output file. A delimiter is appended for + % all but the last element of each row. At the end of a row, a newline + % is written to the output file. + if j < size(rawCellColumns2,2) + fprintf(fileID,['%s',delimiter],output{i,j}); + else + fprintf(fileID,'%s\r\n',output{i,j}); + end + end; + end; + + +fclose(fileID); + +end + + +end +% exit +%type mech_ERC-MultiChem+Bio_Brakora2012_v1_test.inp + +%% Clear temporary variables +% clearvars filename delimiter formatSpec fileID dataArray ans raw col numericData rawCellColumns rawNumericColumns R; \ No newline at end of file diff --git a/flux_version/optimization_7_repro_multiple_2.m b/flux_version/optimization_7_repro_multiple_2.m new file mode 100644 index 0000000..6d2beb0 --- /dev/null +++ b/flux_version/optimization_7_repro_multiple_2.m @@ -0,0 +1,227 @@ +% clear all +% close all + +% classnumb=[15 22 26 27 28]; +% classnumb=[11 15 22 24 26 27 28]; +numbOfClass = length(classnumb); +class_numb_text = {}; +for k=1:numbOfClass + class_numb_text=[class_numb_text ['class',num2str(classnumb(k))]]; +end +fuel_sim={'modify'}; + +%% read modification ignition delay time +% mechanism={'MFC'}; +% mechanism={'MFC'}; +% date = {'03_23_2017_1_iteration'}; +% date = {'03_06_2017_increased'} +fuel_name = {'n_dodecane'}; +% fuel_name = {'n_heptane'}; +equi=1; + +currentloc = '/scratch/engin_flux/unghee/chemkin'; + + +% pressure=[20]; +for k=1:length(pressure) + pressure_text{k}=[num2str(pressure(k)),'atm']; + +end +m = 1; +directory=[fuel_name{1},'_',pressure_text{m},'_','phi',num2str(equi),'_',date{1}]; +% location_rateParam=[currentloc,'/',mechanism{1},'/',directory,'/',fuel_sim{1},'/',class_numb_text{1}]; +location_rateParam=[currentloc,'/',mechanism{1},'/',directory,'/',fuel_sim{1},'/']; +cd(location_rateParam) +load rateParam.mat; +clear m; +%% exp data + +% addpath('C:\Users\unghee\Dropbox\post_process'); +real_fuel_ID; +pure_component_ID; +Target_fuel1 = Vasu_dode_20atm; +Target_data1=Target_fuel1(:,5); +Temp1 = Target_fuel1(:,2); +numbOftarget1 =length(Target_data1); + + +Target_fuel2 = Shen_dode_40atm; +Target_data2=Target_fuel2(:,5); +Temp2 = Target_fuel2(:,2); +numbOftarget2 =length(Target_data2); + +Temp = [Temp1; Temp2]; +Target_data = [Target_data1; Target_data2]; + +% numbOftarget2=0; +% Temp = Temp1; +% Target_data = Target_data1; + +cd ../../.. + + + A = rateParam.(class_numb_text{1})(:,1); + E = rateParam.(class_numb_text{1})(:,2); +%% read modification ignition delay time + +currentloc = '/scratch/engin_flux/unghee/chemkin'; +num_cases_modification= size(A,1); + +m=1; % pressure 20atm +directory=[fuel_name{1},'_',pressure_text{m},'_','phi',num2str(equi),'_',date{1}]; + +for j = 1: numbOftarget1 + for k = 1: numbOfClass + + location_modification=[currentloc,'/',mechanism{1},'/',directory,'/',fuel_sim{1},'/',class_numb_text{k},'/',num2str(j)]; + time_struct_modification=read_ignition_delay(location_modification,num_cases_modification); + time_modification.(class_numb_text{k})(:,j)=time_struct_modification.table.data(:,10); + temp_modification.(class_numb_text{k})(:,j)=time_struct_modification.table.data(:,6); + end + +end +for k = 1: numbOfClass + + time_modification.(class_numb_text{k})=time_modification.(class_numb_text{k})(:,1: numbOftarget1); + temp_modification.(class_numb_text{k})=temp_modification.(class_numb_text{k})(:,1: numbOftarget1); +end + +m =2 ; % pressure 40atm +directory=[fuel_name{1},'_',pressure_text{m},'_','phi',num2str(equi),'_',date{1}]; + +for j = 1: numbOftarget2 + for k = 1: numbOfClass + + location_modification=[currentloc,'/',mechanism{1},'/',directory,'/',fuel_sim{1},'/',class_numb_text{k},'/',num2str(j)]; + time_struct_modification=read_ignition_delay(location_modification,num_cases_modification); + time_modification.(class_numb_text{k})(:,numbOftarget1+j)=time_struct_modification.table.data(:,10); + temp_modification.(class_numb_text{k})(:,numbOftarget1+j)=time_struct_modification.table.data(:,6); + end + +end + + + + + +%% coefficient +for j = 1 : size(Temp,1) + + Temp_current = Temp(j); + + + for k = 1: numbOfClass + time_current = time_modification.(class_numb_text{k})(:,j); + A = rateParam.(class_numb_text{k})(:,1); + E = rateParam.(class_numb_text{k})(:,2); + M = [log(A) log(Temp_current)*ones(7,1) E log(A).*(E) log(A).*log(A) E.^2]; + d = log(time_current); + coefs_inv = lsqlin(M,d); + coefs_element = coefs_inv'; + + + + coefs.(class_numb_text{k})(j,:)=coefs_element; + + prediction.(class_numb_text{k}){:,j}=M*coefs_element'-d; +% predictionreg.(class_numb_text{k}){:,j}=M*coefs_element'; +% plotregression(d,M*coefs_element') + end + + +end +% plot regression +% +% for k = 1: numbOfClass +% plotregression(log(time_modification.(class_numb_text{k})(:,1:10)),predictionreg.(class_numb_text{k})(1:10)) +% end + + + + +%% OPTIMIZER + +% Target Temp value + +numberOftempPoints = size(Target_data,1); +%% objective function + +ObjectiveFunction = @(X) find_rate_weighting(X,coefs,Temp,numbOfClass,class_numb_text,Target_data,numbOftarget1,numbOftarget2); +LB =[]; +UB =[]; +% class2 +if strcmp(mechanism{1},'Ra_Reitz') && ismember(2,classnumb) +LB =[LB rateParam.('class2')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class2')(1,1)*10 rateParam.('class2')(1,2)+2000 ]; + +end +% class4 +if strcmp(mechanism{1},'Ra_Reitz') && ismember(4,classnumb) +LB =[LB rateParam.('class4')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class4')(1,1)*10 rateParam.('class4')(1,2)+2000 ]; +end + +% class6 +if strcmp(mechanism{1},'Ra_Reitz') && ismember(6,classnumb) +LB =[LB rateParam.('class6')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class6')(1,1)*10 rateParam.('class6')(1,2)+2000 ]; +end + +% ndodecane +if strcmp(mechanism{1},'MFC') && ismember(11,classnumb) +LB =[LB rateParam.('class11')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class11')(1,1)*10 rateParam.('class11')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(15,classnumb) +LB =[LB rateParam.('class15')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class15')(1,1)*10 rateParam.('class15')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(21,classnumb) +LB =[LB rateParam.('class21')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class21')(1,1)*10 rateParam.('class21')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(22,classnumb) +LB =[LB rateParam.('class22')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class22')(1,1)*10 rateParam.('class22')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(23,classnumb) +LB =[LB rateParam.('class23')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class23')(1,1)*10 rateParam.('class23')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(24,classnumb) +LB =[LB rateParam.('class24')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class24')(1,1)*10 rateParam.('class24')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(26,classnumb) +LB =[LB rateParam.('class26')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class26')(1,1)*10 rateParam.('class26')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(27,classnumb) +LB =[LB rateParam.('class27')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class27')(1,1)*10 rateParam.('class27')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(28,classnumb) +LB =[LB rateParam.('class28')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class28')(1,1)*10 rateParam.('class28')(1,2)+2000 ]; +end + +nvars=2*numbOfClass; +options=gaoptimset('PopulationSize',500); +[result_ga,Fval,exitFlag,Output] = ga(ObjectiveFunction,nvars,[],[],[],[],LB,UB,[],options); +[result_fmin,Fval,exitFlag,Output] = fmincon(ObjectiveFunction,result_ga,[],[],[],[],LB,UB); + +X = result_ga; + +error=ObjectiveFunction(X) + +for i = 1 : numbOfClass +final_result.(class_numb_text{i})= [result_ga(1,2*i-1:2*i); result_fmin(1,2*i-1:2*i)]; +end + +% location_save=[currentloc,'\',mechanism{1},'\',directory]; +location_save='/scratch/engin_flux/unghee/chemkin/MFC'; +cd(location_save) +% cd ../ + +save('final_result.mat','final_result') + diff --git a/flux_version/optimization_driver_2.m b/flux_version/optimization_driver_2.m new file mode 100644 index 0000000..864d8df --- /dev/null +++ b/flux_version/optimization_driver_2.m @@ -0,0 +1,55 @@ +clear; +mechanism={'MFC'}; +fuel_name={'n_dodecane'}; +date = {'03_23_2017_1_iteration'}; +fuel_sim={'modify_sensitivity_3'}; +equi = 1; +classnumb=[11 15 21 22 23 24 26 27 28]; +pressure=[20 40]; +%% Target setting +addpath('/scratch/engin_flux/unghee/chemkin/MFC'); +real_fuel_ID; +pure_component_ID; +Target_fuel1 = Vasu_dode_20atm; + +Target_data1=Target_fuel1(:,5); +numbOftarget1=length(Target_data1); +Temp1=Target_fuel1(:,2); + +Target_fuel2 = Shen_dode_40atm; +Target_data2=Target_fuel2(:,5); +Temp2 = Target_fuel2(:,2); +numbOftarget2=length(Target_data2); + +Temp = [Temp1; Temp2]; +Target_data = [Target_data1; Target_data2]; + +%% sensitivity +% addpath('C:\Users\unghee\Dropbox\JPoptimization_github'); +sensitivity_analysis_3; +% sensitivity_analysis_3_modify; +% clearvars -except sensitivity classnumb classnumb_text pressure_text pressure mechanism fuel_name date fuel_sim equi + +m =1; +disp('looking at sensitivity at one pressure condition') +class_to_optimize.(pressure_text{m}) = []; +for k=1:length(classnumb) + + if (sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig_avg > 1)... + && (sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr_avg >= 0.38) + class_to_optimize.(pressure_text{m}) = [class_to_optimize.(pressure_text{m}) classnumb(k)]; + end +end + +classnumb = class_to_optimize.(pressure_text{m}) + + + +%% Optimization +% classnumb= [26 27] +clearvars -except classnumb Temp Temp1 Temp2 numberOftarget Target_data Target_data1 numbOftarget1 sensitivity pressure fuel_name equi fuel_name mechanism fuel_sim date +optimization_7_repro_multiple_2 + +%% call modification func. if the tolerance is larger than.. and loop back to optimization driver + +% if error < k , run optimization_driver_2 diff --git a/flux_version/pure_component_ID.m b/flux_version/pure_component_ID.m new file mode 100644 index 0000000..8728257 --- /dev/null +++ b/flux_version/pure_component_ID.m @@ -0,0 +1,290 @@ +Pfahl_dec_20atm=[ + 1297.87234 0.770491803 12 90.83238519 54.49943112 +1234.817814 0.809836066 12 189.938161 113.9628966 +1204.235463 0.830402385 12 251.9098497 151.1459098 +1095.689092 0.91266766 12 870.0784243 522.0470546 +1066.433566 0.937704918 12 791.2576649 474.7545989 +1013.289037 0.986885246 12 1680.364423 1008.218654 +936.1049107 1.068256334 12 1931.573337 1158.944002 +932.2033898 1.072727273 12 2816.325604 1689.795362 +915.4160982 1.092399404 12 3735.30403 2241.182418 +857.1793562 1.16661699 12 2765.596281 1659.357768 +836.6583541 1.195230999 12 2989.461373 1793.676824 +798.4293194 1.252459016 12 2471.89845 1483.13907 +764.0628558 1.308792846 12 2076.368423 1245.821054 +731.0960994 1.36780924 12 2464.440492 1478.664295 +697.3602162 1.433979136 12 3587.468327 2152.480996 +]; + +Shen_dec_20atm=[ + 1011 0.989119683 12 1135 681 +1041 0.960614793 11 881 484.55 +1085 0.921658986 11.7 615 359.775 +1088 0.919117647 11.4 513 292.41 +1145 0.873362445 10.4 359 186.68 +1175 0.85106383 10.9 311 169.495 +1263 0.791765637 10.5 108 56.7 +1276 0.78369906 9.9 114 56.43 +1378 0.725689405 10.6 40 21.2 +836 1.196172249 37.2 757 1408.02 +900 1.111111111 43.1 692 1491.26 +928 1.077586207 43.1 625 1346.875 +929 1.076426265 42.5 551 1170.875 +938 1.066098081 42.4 565 1197.8 +955 1.047120419 42.5 597 1268.625 +970 1.030927835 45.6 397 905.16 +976 1.024590164 43.3 471 1019.715 +1003 0.997008973 41.9 364 762.58 +1032 0.968992248 39.5 351 693.225 +1079 0.926784059 39.8 186 370.14 +1133 0.882612533 38.9 135 262.575 +1144 0.874125874 33.9 128 216.96 +1158 0.863557858 37.7 113 213.005 +1271 0.786782061 35.8 37 66.23 +]; + + +Fieweger_ioct_20atm=[ %% + 1208.053691 0.827777778 33.6 67.64376051 100.8592862 +1113.402062 0.898148148 33.6 293.7033517 437.9222885 +990.8256881 1.009259259 33.6 1148.263098 1712.101685 +980.9264305 1.019444444 33.6 1159.263587 1728.503812 +961.7097061 1.039814815 33.6 1251.147431 1865.505937 +926.2435678 1.07962963 33.6 2195.911873 3274.183788 +917.5870858 1.089814815 33.6 2632.024061 3924.442787 +800 1.25 33.6 10388.74645 15489.99558 +781.4761216 1.27962963 33.6 8107.793734 12089.01283 +729.7297297 1.37037037 33.6 9004.328811 13425.77896 +699.0291262 1.430555556 33.6 10792.60528 16092.16367 +1186.813187 0.842592593 16.8 191.2376943 167.2124265 +1055.718475 0.947222222 16.8 1105.295141 966.4364717 +980.0362976 1.02037037 16.8 2392.6655 2092.073979 +951.5418502 1.050925926 16.8 2813.684268 2460.199991 +]; + +Davidson_ioct_20atm=[ + 984 1.016260163 18.1 1511 1396.422513 +995 1.005025126 16.3 1535 1305.938997 +1043 0.958772771 17.1 927 819.0925291 +1077 0.928505107 18.4 604 565.4957199 +1109 0.901713255 15.9 516 430.4668708 +1159 0.86281277 14.9 214 169.5966341 +]; + +Shen_ioct_20atm=[ +1021 0.979431929 19.8 1059 1051.04753 +1032 0.968992248 25.5 850 1019.885845 +1045 0.956937799 26 645 785.2667863 +1058 0.945179584 24.5 704 819.7373542 +1072 0.932835821 22.8 509 561.5603303 +1078 0.927643785 23.9 439 501.7535446 +1084 0.922509225 24.2 494 569.9227281 +1085 0.921658986 24.1 429 493.3982188 +1088 0.919117647 25.9 449 545.0654663 +1118 0.894454383 23.5 295 332.9280721 +1132 0.883392226 23.2 220 245.9043281 +1142 0.875656743 19 220 211.6973321 +1145 0.873362445 21.3 207 217.0114013 +1205 0.829875519 20.7 104 106.7182271 +1223 0.817661488 20.4 91 92.36161562 +]; + +Oehlschlaeger_hmn_20atm=[ + 993 1.007049345 14 984 739.7305083 +972 1.028806584 31.2 813 1160.353253 +1039 0.962463908 31.8 299 433.3001983 +1074 0.931098696 27 193 245.3716164 +1083 0.923361034 25.4 183 221.5613826 +1179 0.848176421 28.5 63 83.63584343 +1195 0.836820084 26.5 62 77.65410026 +]; + +% Oehlschlaeger_hmn_20atm=[ +% 993 14 984 1.007049345 688.8 +% 1046 11.7 662 0.956022945 387.27 +% 1098 9.8 472 0.910746812 231.28 +% 1161 12.5 173 0.861326443 108.125 +% 1210 10.9 119 0.826446281 64.855 +% 1304 10.4 59 0.766871166 30.68 +% 1309 9.9 49 0.76394194 24.255 +% 953 46.9 1014 1.049317943 2377.83 +% 968 37.1 963 1.033057851 1786.365 +% 972 31.2 813 1.028806584 1268.28 +% 986 43.8 600 1.014198783 1314 +% 1039 31.8 299 0.962463908 475.41 +% 1074 27 193 0.931098696 260.55 +% 1083 25.4 183 0.923361034 232.41 +% 1179 28.5 63 0.848176421 89.775 +% 1195 26.5 62 0.836820084 82.15 +% ]; + +Oehlschlaeger_hmn_40atm=[ + 953 1.049317943 46.9 1014 1151.670483 +968 1.033057851 37.1 963 906.728824 +972 1.028806584 31.2 813 666.4479367 +986 1.014198783 43.8 600 645.1824505 +1039 0.962463908 31.8 299 248.8656125 +1074 0.931098696 27 193 140.928986 +1083 0.923361034 25.4 183 127.2535978 +1179 0.848176421 28.5 63 48.03617788 +1195 0.836820084 26.5 62 44.60056861 +]; + +Shen_dode_20atm=[ +924 1.082251082 17.8 1299 1156.11 +999 1.001001001 12.1 822 497.31 +1014 0.986193294 13.9 713 495.535 +1014 0.986193294 13.8 740 510.6 +1030 0.970873786 12.3 656 403.44 + +1103 0.906618314 12.8 343 219.52 + +1135 0.881057269 13.9 257 178.615 + +1177 0.849617672 14.9 162 120.69 +1178 0.848896435 12.4 194 120.28 +1210 0.826446281 13.4 123 82.41]; +Shen_dode_40atm=[ +877 1.140250855 45 801 901.125 +913 1.095290252 41.8 757 791.065 +945 1.058201058 43.2 737 795.96 +978 1.022494888 43.8 506 554.07 +1020 0.980392157 45.2 342 386.46 +1043 0.958772771 41 270 276.75 +1080 0.925925926 40.6 206 209.09 +1097 0.911577028 41.3 157 162.1025 +1102 0.907441016 44.2 155 171.275 +1114 0.897666068 38.6 129 124.485 +1122 0.891265597 42.6 120 127.8 +]; +Vasu_dode_20atm=[ + 727 1.375515818 27 809 1092.15 +773 1.293661061 28.7 556 797.86 +818 1.222493888 26.9 881 1184.945 +822 1.216545012 23.3 805 937.825 +855 1.169590643 25.9 875 1133.125 +869 1.150747986 22.5 1040 1170 +907 1.102535832 23.4 1081 1264.77 +942 1.061571125 21 1116 1171.8 +953 1.049317943 20 1141 1141 +957 1.044932079 19.8 1064 1053.36 +976 1.024590164 20.5 800 820 +978 1.022494888 22.9 912 1044.24 +987 1.013171226 22 699 768.9 +991 1.009081736 21.8 645 703.05 +992 1.008064516 20.4 739 753.78 +1008 0.992063492 22.4 570 638.4 +1015 0.985221675 22.3 508 566.42 +1036 0.965250965 22.1 432 477.36 +1087 0.919963201 19.4 403 390.91 +1102 0.907441016 20.8 268 278.72 +1107 0.903342367 23.3 298 347.17 +1109 0.901713255 23 278 319.7 +1118 0.894454383 33.7 183 308.355 +1123 0.89047195 24 262 314.4 +1125 0.888888889 18.7 188 175.78 +1135 0.881057269 20.6 178 183.34 +]; +Oehlschlaeger_decalin_20atm=[ +1037 0.964320154 15.2 867 699.9285351 +1049 0.953288847 13.5 871 641.0251926 +1057 0.946073794 12.5 825 571.7950114 +1063 0.940733772 13.9 627 472.0800943 +1074 0.931098696 12.2 580 394.4439188 +1105 0.904977376 14.1 371 282.4629047 +1134 0.881834215 10.5 315 190.5613454 +1164 0.859106529 11.4 200 129.0067348]; +Zhu_decalin_20atm = [ + 1202 0.831946755 16.6 109 90.47 +1193 0.838222967 18.8 106 99.64 +1145 0.873362445 20.5 157 160.925 +1100 0.909090909 18.9 284 268.38 +1090 0.917431193 19.9 261 259.695 +1068 0.936329588 19.2 381 365.76 +1059 0.944287063 20.5 396 405.9 +1051 0.951474786 19.3 440 424.6 +1014 0.986193294 26.8 553 741.02 +1012 0.988142292 23.3 788 918.02 +990 1.01010101 22.1 917 1013.285 +951 1.051524711 24.1 1271 1531.555 +914 1.094091904 19.5 2644 2577.9 +876 1.141552511 20.2 3291 3323.91 +830 1.204819277 24 3049 3658.8 +802 1.246882793 19.3 5341 5154.065 +769 1.300390117 15.5 8746 6778.15 +]; +Zhu_decalin_50atm = [ + 1043 0.958772771 46 264 242.88 +1011 0.989119683 50.4 413 416.304 +1006 0.994035785 46.4 443 411.104 +962 1.03950104 47.9 891 853.578 +961 1.040582726 50.1 804 805.608 +941 1.062699256 50.6 1024 1036.288 +912 1.096491228 48.9 903 883.134 +883 1.132502831 56.6 638 722.216 +852 1.17370892 37.8 1984 1499.904 +831 1.203369434 48.3 3411 3295.026 +]; + +Shen_tol_20atm=[ + 1129 0.885739593 13.4 1380 891.8680812 +1164 0.859106529 13.5 934 608.5384005 +1197 0.835421888 12.4 798 473.9252337 +1235 0.809716599 11.3 483 259.2268239 +1239 0.807102502 14.5 471 331.7335027 +1290 0.775193798 11.2 281 149.358956 +1371 0.729394602 13.2 137 87.10106333 +]; +Davidson_tol_20atm=[ + 971 1.029866117 17.2 1314 1142.033714 +995 1.005025126 16.5 1100 919.8030754 +1053 0.949667616 17.1 1288 1113.382395 +1067 0.937207123 16.1 1646 1345.302656 +1076 0.92936803 16.5 1174 981.6807368 +1097 0.911577028 16.6 1223 1028.416611 +1138 0.878734622 22.2 716 788.9752734 +1173 0.852514919 14.9 767 583.3117135 +]; + +Shen_tol_50atm = [ + 1021 0.979431929 53.8 1255 1350.38 +1041 0.960614793 47.6 1203 1145.256 +1082 0.924214418 46.9 757 710.066 +1085 0.921658986 51 546 556.92 +1116 0.896057348 49.5 492 487.08 +1117 0.895255148 46.2 459 424.116 +1135 0.881057269 50.9 346 352.228 +1218 0.821018062 51.4 117 120.276 +1232 0.811688312 57.5 97 111.55 +1256 0.796178344 61.5 61 75.03 +1305 0.766283525 60.7 44 53.416 +]; + +Davidson_tol_50atm = [ + 977 1.023541453 54.3 1173 1273.878 +990 1.01010101 50.5 952 961.52 +1005 0.995024876 50.4 870 876.96 +1068 0.936329588 48 605 580.8 +1111 0.900090009 46.7 397 370.798 +1166 0.857632933 44.6 276 246.192 +1179 0.848176421 41.5 273 226.59 +1183 0.845308538 46.2 191 176.484 +1242 0.805152979 42.4 121 102.608 +]; + +%n heptane +Shen_hep_40atm=[ +788 1.269035533 295 +797 1.254705144 379 +799 1.251564456 319 +834 1.199040767 284 +862 1.160092807 295 +899 1.112347052 381 +967 1.034126163 491 +991 1.009081736 447 +996 1.004016064 336 +1003 0.997008973 344 +1022 0.978473581 282 +1029 0.971817298 297 +1050 0.952380952 259]; + diff --git a/flux_version/read_ignition_delay.m b/flux_version/read_ignition_delay.m new file mode 100644 index 0000000..13b2cbe --- /dev/null +++ b/flux_version/read_ignition_delay.m @@ -0,0 +1,138 @@ +%updated for DTIG ouput - 3/28/2014 +%For MFL_2013, ignition delay prediction based on delta T of 400K is added +%to data output file. To accomodate this, code is updated. + + +function A = read_ignition_delay(location,num_cases) + + + if ~isempty(findstr(location,'MFL_2013')) || ~isempty(findstr(location,'MFL_2014')) || ~isempty(findstr(location,'MFC_2012'))||... + ~isempty(findstr(location,'Ra_Reitz')) || ~isempty(findstr(location,'MFC')) + % if the mechanism was wither MFL_2013, MFL_2014, or MFC_2012 + + for i=1:num_cases +% file_path=[location,'/id_',num2str(i),'/CKSoln_solution_point_value_vs_solution_number.csv']; + % location of the file +% cd([location,'/id_',num2str(i)]) + +% [data{i}.num data{i}.txt data{i}.raw]=csvread('CKSoln_solution_point_value_vs_solution_number.csv',1,0); + + + cd([location,'/id_',num2str(i)]) + + num=csvread('CKSoln_solution_point_value_vs_solution_number.csv',1,0); +% file_path=['/scratch/engin_flux/unghee/chemkin/MFC/n_dodecane_20atm_phi1_03_16_2017_5_iteration/modify_sensitivity/class11/1/id_1/','CKSoln_solution_point_value_vs_solution_number.csv']; + file_path=[location,'/id_',num2str(i),'/','CKSoln_solution_point_value_vs_solution_number.csv']; + formatSpec = '%s%s%s%s%s%s%s%s%s%s%s%[^\n\r]'; + delimiter = ','; + fileID=fopen(file_path,'r'); + dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'ReturnOnError', false); + for j = 1 : 12 + text{1,j} = dataArray{1,j}(1,1); + + rawdat{1,j} = dataArray{1,j}(2,1); + end + raw = [text{1,:}; rawdat{1,:}]; + + data{i}.num = num; + data{i}.txt = text; + data{i}.raw = raw; + % read the file + + fclose(fileID); + + + + + table.data(i,1:6)=data{i}.num(1:6); + % keep the first 6 columns + + [m n]=size(data{i}.num); + + if n<7 % when the data file is only 6 columns wide, igntion was not achieved (or + % both DT and OH peak ignition criteria were not met). + display(['Warning: no ID data for ',location,', case#',num2str(i)]); + table.data(i,7:8)=0; + % then send warning and ignition delay time is set to zero. + + else if n==7 && strcmp(data{i}.txt{7},' Ignition_time_1_by_ignition_T_(sec)') + % when the data file is 7 columns wide and the header + % of 7th column is as above, ignition delay time is + % from DTIG criteria only, and OH peak criteria is not + % met. + display(['Warning: Only DTIG ID data for ',location,', case#',num2str(i)]); + table.data(i,7)=data{i}.num(7); %delta T + table.data(i,8)=data{i}.num(7); %delta T + % Then, send warning and ignition delay time for OH + % peak is same as DTIG igintion delay time. + else + + % For the rest of the case, 7th column is DTIG ignition + % delay time, 8th ~ end columns are OH peak ignition + % delay time. OH peak ignition delay time can be + % multiple, since it looks for local minimum. + % Typically, the last one is responsible for the real + % ignition delay time, while + % + table.data(i,7)=data{i}.num(7); %delta T + table.data(i,8)=max(data{i}.num(8:end)); %oh peak + end + end + end + + + table.data=sortrows(table.data,6); %sort based on 1000/T + + temp=table.data(:,7:8); % to rearrange the table + + table.data(:,7)=temp(:,2); % oh peak to 7th and 8th column + table.data(:,8)=table.data(:,7)*1000000; + table.data(:,9)=temp(:,1); % delta T to 9th and 10th column + table.data(:,10)=table.data(:,9)*1000000; + + for j=1:6 + table.header{j}=data{1}.txt{j}; + end + table.header{7}='ignition_time_oh(s)'; + table.header{8}='ignition_time_oh(micros)'; + table.header{9}='ignition_time_dt(s)'; + table.header{10}='ignition_time_dt(micros)'; + + + + + + else % for rest of the mechanisms (previous) + +% for i=1:num_cases +% file_path=[location,'\id_',num2str(i),'\CKSoln_solution_point_value_vs_solution_number.csv']; +% [data{i}.num data{i}.txt data{i}.raw]=xlsread(file_path); +% table.data(i,1:6)=data{i}.num(1:6); +% [m n]=size(data{i}.num); +% if n>6 +% table.data(i,7)=max(data{i}.num(7:end)); +% else +% table.data(i,7)=0; +% +% display(['Warning: no ID data for ',location,', case#',num2str(i)]); +% +% end +% end +% table.data=sortrows(table.data,6); +% table.data(:,8)=table.data(:,7)*1000000; +% for j=1:6 +% table.header{j}=data{1}.txt{j}; +% end +% table.header{7}='ignition_time_(s)'; +% table.header{8}='ignition_time_(micros)'; + + end + + A.raw=data; + A.table=table; + +end + + + + diff --git a/flux_version/real_fuel_ID.m b/flux_version/real_fuel_ID.m new file mode 100644 index 0000000..52f426e --- /dev/null +++ b/flux_version/real_fuel_ID.m @@ -0,0 +1,226 @@ +% experimental data +% Jet-A 4658 +% Vasu -> Vasu et al. +% Wang -> Wang and Oehlschlaeger +% Dooley -> Dooley et al. + +% T 1000/T P ID_raw ID_scaled +% scaling factor of -1 was used + +Vasu_20atm=[ + 715 1.399 21.7 3286 3565.31 +753 1.328 21.8 2584 2816.56 +763 1.311 21.5 2423 2604.725 +774 1.292 19.4 2566 2489.02 +785 1.274 21.1 2484 2620.62 +786 1.272 19.7 2742 2700.87 +808 1.238 20.8 2628 2733.12 +874 1.144 23.3 3109 3621.985 +880 1.136 25.1 2355 2955.525 +952 1.05 24 1358 1629.6 +958 1.044 22.5 1287 1447.875 +963 1.038 21.7 1412 1532.02 +987 1.013 21.3 984 1047.96 +1035 0.966 21.9 539 590.205 +1042 0.96 22 486 534.6 +1130 0.885 24.8 138 171.12 +]; + +Vasu_40atm=[ +934 1.071 31.4 1230 965.55 +974 1.027 32.5 777 631.3125 +997 1.003 32.1 641 514.4025 +961 1.041 49.3 528 650.76 +1009 0.991 50.9 286 363.935 +1090 0.917 50.5 115 145.1875 +1124 0.89 47.8 80 95.6 +]; + +Wang_20atm=[ + 674 1.483679525 20.8 8974 9332.96 +680 1.470588235 21.8 6703 7306.27 +681 1.468428781 21.3 5821 6199.365 +700 1.428571429 23.5 3887 4567.225 +703 1.422475107 23.7 4216 4995.96 +704 1.420454545 23.4 3204 3748.68 +717 1.394700139 24.3 2737 3325.455 +723 1.383125864 20 2287 2287 +737 1.356852103 23.7 1871 2217.135 +750 1.333333333 20.9 1842 1924.89 +753 1.328021248 20.8 1894 1969.76 +755 1.324503311 21.2 1788 1895.28 +757 1.321003963 20.5 1757 1800.925 +797 1.254705144 23 1618 1860.7 +797 1.254705144 24.8 1618 2006.32 +802 1.246882793 21.5 1616 1737.2 +803 1.245330012 23.2 1625 1885 +815 1.226993865 22.1 1742 1924.91 +820 1.219512195 19.1 2102 2007.41 +826 1.210653753 21.2 1896 2009.76 +831 1.203369434 21.7 1789 1941.065 +836 1.196172249 20.5 2038 2088.95 +845 1.183431953 20.9 2138 2234.21 +849 1.177856302 21.2 2059 2182.54 +849 1.177856302 22.1 1897 2096.185 +851 1.175088132 21.8 1908 2079.72 +871 1.148105626 22.8 1843 2101.02 +882 1.133786848 21.6 2032 2194.56 +904 1.10619469 21.5 2031 2183.325 +931 1.074113856 20.9 1749 1827.705 +940 1.063829787 22.7 1675 1901.125 +955 1.047120419 21.1 1354 1428.47 +973 1.027749229 21.7 1236 1341.06 +984 1.016260163 19.9 1140 1134.3 +991 1.009081736 18.9 1050 992.25 +1025 0.975609756 19.2 816 783.36 +1043 0.958772771 16.7 706 589.51 +1046 0.956022945 21.6 563 608.04 +1050 0.952380952 21 458 480.9 +1072 0.932835821 19.3 464 447.76 +1083 0.923361034 19.1 375 358.125 +1092 0.915750916 18.4 318 292.56 +1100 0.909090909 16.9 355 299.975 +1106 0.904159132 18.2 261 237.51 +1109 0.901713255 18.1 249 225.345 +1117 0.895255148 17.7 225 199.125 +1126 0.888099467 19.2 185 177.6 +1138 0.878734622 18.7 181 169.235 +1143 0.874890639 18.1 176 159.28 +1145 0.873362445 17.2 184 158.24 +1177 0.849617672 18.1 104 94.12 +1187 0.842459983 16.3 116 94.54 +1197 0.835421888 16.5 95 78.375 +1222 0.818330606 16.8 77 64.68 +]; + +Wang_40atm=[ + 741 1.349527665 36.6 1885 1724.775 +763 1.31061599 40.3 1707 1719.8025 +812 1.231527094 38.5 1644 1582.35 +832 1.201923077 36 1661 1494.9 +848 1.179245283 40.1 1727 1731.3175 +852 1.17370892 36.9 1423 1312.7175 +854 1.170960187 36.9 1423 1312.7175 +855 1.169590643 35.7 1824 1627.92 +888 1.126126126 39.3 1550 1522.875 +896 1.116071429 38.1 1640 1562.1 +904 1.10619469 34.5 1689 1456.7625 +916 1.091703057 39.9 1483 1479.2925 +945 1.058201058 36.5 1192 1087.7 +964 1.037344398 41.1 819 841.5225 +978 1.022494888 40.2 821 825.105 +1007 0.993048659 38.8 542 525.74 +1008 0.992063492 40.3 689 694.1675 +1021 0.979431929 38.8 462 448.14 +1030 0.970873786 38.6 456 440.04 +1034 0.967117988 40.9 495 506.1375 +1054 0.948766603 38.5 298 286.825 +1103 0.906618314 40.2 169 169.845 +1132 0.883392226 39.6 107 105.93 +1161 0.861326443 40.3 92 92.69 +1162 0.860585198 40.9 79 80.7775 +1163 0.859845228 38.7 96 92.88 +1224 0.816993464 40.9 40 40.9 +1229 0.81366965 38.8 38 36.86 +]; + +Dooley_20atm=[ + 653.3 1.530690341 22 24100 26510 +663.4 1.507386192 22 18700 20570 +671.5 1.489203276 22 12980 14278 +680 1.470588235 22 8380 9218 +690.8 1.447596989 22 6600 7260 +699.4 1.429796969 22 4600 5060 +709.6 1.409244645 22 3480 3828 +]; + +Wang_20atm_rich=[ + 662 1.510574018 19.8 8311 8227.89 +670 1.492537313 20 6770 6770 +672 1.488095238 20.1 5285 5311.425 +684 1.461988304 21.1 3656 3857.08 +701 1.426533524 18.8 2811 2642.34 +706 1.416430595 23 2528 2907.2 +710 1.408450704 23.5 2331 2738.925 +717 1.394700139 19.7 1720 1694.2 +724 1.38121547 21.1 1973 2081.515 +743 1.34589502 18.7 1709 1597.915 +759 1.317523057 19.6 1832 1795.36 +774 1.291989664 20.6 1481 1525.43 +810 1.234567901 21 1678 1761.9 +813 1.2300123 18.4 1651 1518.92 +832 1.201923077 19.6 2052 2010.96 +837 1.19474313 19.1 1761 1681.755 +854 1.170960187 19.2 2064 1981.44 +880 1.136363636 22.2 1655 1837.05 +894 1.118568233 23.4 1662 1944.54 +901 1.109877913 20.4 1931 1969.62 +910 1.098901099 20.5 1612 1652.3 +972 1.028806584 20.2 1024 1034.24 +989 1.011122346 23.2 816 946.56 +1004 0.996015936 18.9 756 714.42 +1014 0.986193294 21.2 628 665.68 +1014 0.986193294 21.7 522 566.37 +1020 0.980392157 20.9 524 547.58 +1021 0.979431929 20.4 541 551.82 +1049 0.953288847 22.4 371 415.52 +1075 0.930232558 21.9 279 305.505 +1088 0.919117647 22.3 291 324.465 +1101 0.908265213 21.9 243 266.085 +1121 0.89206066 23.4 140 163.8 +1142 0.875656743 20.5 130 133.25 +1132 0.883392226 18.9 155 146.475 +1157 0.864304235 20.3 119 120.785 +1160 0.862068966 19.9 114 113.43 +1185 0.843881857 21.1 83 87.565 +]; + +Wang_IPK_20atm=[ +681 1.468428781 21.5 6642 7140.15 +687 1.455604076 20.4 5995 6114.9 +734 1.36239782 20.4 3060 3121.2 +805 1.242236025 20.6 3332 3431.96 +869 1.150747986 20 2686 2686 +914 1.094091904 20.2 1910 1929.1 +977 1.023541453 20.2 1162 1173.62 +989 1.011122346 20 1038 1038 +1033 0.968054211 20.4 563 574.26 +1080 0.925925926 21.5 303 325.725 +1151 0.868809731 20.3 135 137.025 +1213 0.824402308 20.5 70 71.75 +1252 0.798722045 20.8 45 46.8 +1323 0.755857899 21.9 27 29.565]; + +Wang_S8_20atm=[ +673 1.485884101 20.2 2676 2702.76 +700 1.428571429 22.5 1621 1823.625 +717 1.394700139 20.6 1334 1374.02 +732 1.366120219 20.1 1245 1251.225 +759 1.317523057 19.6 1453 1423.94 +804 1.243781095 20.3 1457 1478.855 +832 1.201923077 20 1295 1295 +857 1.166861144 19.8 1674 1657.26 +870 1.149425287 18.7 1547 1446.445 +926 1.079913607 18.7 1362 1273.47 +942 1.061571125 20.9 1142 1193.39 +956 1.046025105 20.1 1037 1042.185 +1016 0.984251969 19.6 593 581.14 +1044 0.957854406 24 423 507.6 +1046 0.956022945 20.8 387 402.48 +1066 0.938086304 20 343 343 +1086 0.920810313 19 292 277.4 +1119 0.893655049 20.3 221 224.315 +1142 0.875656743 19.8 161 159.39 +1179 0.848176421 19.8 103 101.97 +1180 0.847457627 20 83 83 +1260 0.793650794 20.8 39 40.56]; + + + + + + + + + + diff --git a/flux_version/sensitivity_analysis_3.m b/flux_version/sensitivity_analysis_3.m new file mode 100644 index 0000000..12e6f29 --- /dev/null +++ b/flux_version/sensitivity_analysis_3.m @@ -0,0 +1,159 @@ + +%1. matlab peaks func -> temp ????? sensitivity data (Before,after 25temps) - 4?? + +% pressure=[20 40]; + +numbOfPressure=length(pressure) ; +for k=1:length(pressure) + pressure_text{k}=['P',num2str(pressure(k)),'atm']; +end +% classnumb=[1 2 3 4 5 6 7 8 9]; +numbOfClass=length(classnumb) ; +for k=1:numbOfClass + classnumb_text{k}=['class',num2str(classnumb(k))]; +end + +exp_markers={'ro' 'g*' 'g^'}; +sim_line={'k-' 'k-s','k-d','k-^','k-x','k-o','k-p'}; +marker_size=8; +line_width=2; + + +%%%%setting%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % i=1 ; +% m=1 ; % pressure +% % j=1 ; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% sensitivity analysis +% mechanism={'Ra_Reitz'}; +% fuel_sim={'modify_sensitivity'}; +% fuel_name={'n_heptane'} +% equi = 1; +currentloc = '/scratch/engin_flux/unghee/chemkin/MFC'; +num_cases_modification= 3; + + +% for m = 1 : numbOfPressure +for m = 1 : 1 %% sensitivity analysis for only one pressure + for k = 1 :numbOfClass + for j = 1: 25 + location_modification=['/scratch/engin_flux/unghee/chemkin','/',mechanism{1},'/',fuel_name{1},'_',num2str(pressure(m)),... + 'atm','_','phi',num2str(equi),'_',date{1},'/',fuel_sim{1},'/',classnumb_text{k},'/',num2str(j)]; + time_struct_modification=read_ignition_delay(location_modification,num_cases_modification); + sensitivity.(pressure_text{m}).(classnumb_text{k}).data(:,j)=time_struct_modification.table.data(:,10); + + end + + + temp=700:25:1300; + h =25; + tbase = sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,:); % when base + tk1 = sensitivity.(pressure_text{m}).(classnumb_text{k}).data(2,:); % when k = 2 + tk2 = sensitivity.(pressure_text{m}).(classnumb_text{k}).data(3,:); % when k = 0.5 + d_log_tk1=diff(log10(tk1)); + d_log_tk2=diff(log10(tk2)); + dh = log(h); + k1 = 2; + k2 = 0.5; + + [peaks,locs] = findpeaks(tbase/1000); +% yvalue=tbase/1000; + locsValue = temp(locs); + [peaks2,locs2] = findpeaks(-tbase/1000); +% yvalue=tbase/1000; + locsValue2 = temp(locs2); + + p_tempPoints_ig =[1, locs2, locs, length(temp)]; + temp_p_tempPoints_ig=[temp(1),temp(locs2),temp(locs), temp(length(temp))]; + +%translational sensitivity +%change the normalize contant log10?e.x. log10^3 + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig = ones(1,length(p_tempPoints_ig)); + for i = 1 : length(p_tempPoints_ig) + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig(i)... + = 100*(log10(tk1(p_tempPoints_ig(i))) - log10(tk2(p_tempPoints_ig(i))))/(log10(tbase(i))*log10(k1/k2)); % translational shift + end + + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig_avg... + = mean(abs(sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig)); + + + p_tempPoints_gr =[1, locs2, locs, length(temp)]; % added end point + + temp_p_tempPoints_gr=[temp(locs2),temp(locs), temp(length(temp))]; + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr = ones(1,length(p_tempPoints_gr)-1); + %rotational sensitivity + for i = 1 : length(p_tempPoints_gr)-1 + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr(i)... + = 100*(d_log_tk1(p_tempPoints_gr(i))/dh - d_log_tk2(p_tempPoints_gr(i))/dh)/log10(k1/k2); + end + + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr_avg... + = mean(abs(sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr)); + end + +end + + +save('sensitivity.mat','sensitivity') + +%% plotting +% for m = 1 : numbOfPressure +% for k = 1 : numbOfClass +% h=figure('position',[20 50 1200 480]); +% subplot(1,2,1); +% set(gca,'Fontsize',13) +% semilogy(temp,sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,:)/1000,'k-','markersize',marker_size); +% hold on +% semilogy(temp,sensitivity.(pressure_text{m}).(classnumb_text{k}).data(2,:)/1000,'rs-','markersize',marker_size); +% hold on +% semilogy(temp,sensitivity.(pressure_text{m}).(classnumb_text{k}).data(3,:)/1000,'bx-','markersize',marker_size); +% +% +% semilogy([temp_p_tempPoints_ig(1),temp_p_tempPoints_ig(end)],... +% [sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,1)/1000,... +% sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,end)/1000]... +% ,'gp','markersize',12); +% hold on +% semilogy(locsValue,peaks,'gp','markersize',12); +% hold on +% semilogy(locsValue2,-peaks2,'gp','markersize',12); +% +% legend('baseline','k=2','k=1','extreme points') +% xlabel('1000/T (1/K)') +% ylabel('Ignition Delay Time (ms)') +% +% +% +% +% subplot(1,2,2); +% plot(temp_p_tempPoints_ig,sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig,'k^-','markersize',marker_size); +% ylim([-25 25]) +% hold on +% % figure +% plot(temp_p_tempPoints_gr,sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr,'ko-','markersize',marker_size); +% +% ylim([-30 25]) +% xlabel('1000/T (1/K)') +% ylabel('Ignition Delay Time (ms)') +% legend('ignition delay sensitivity','gradient sensitivity') +% +% annotation(h,'textbox',[0.213 0.38 0.279 0.05],... +% 'String',{fuel_name{1},... +% '/Air',classnumb_text{k}, '\phi=1',pressure_text{m}},... +% 'FontSize',13,... +% 'FontName','Arial',... +% 'FitBoxToText','off',... +% 'LineStyle','none'); +% +% % location_save=[currentloc,'\',mechanism{1},'\',directory]; +% % cd(location_save) +% +% mkdir('sensitivity'); +% cd('sensitivity'); +% saveas(h,classnumb_text{k},'fig') +% saveas(h,classnumb_text{k},'jpg') +% cd ../ +% end +% end \ No newline at end of file diff --git a/images/optimizer_flow.png b/images/optimizer_flow.png new file mode 100644 index 0000000..dd73630 Binary files /dev/null and b/images/optimizer_flow.png differ diff --git a/modification.m b/modification.m index f744b93..bb9394e 100644 --- a/modification.m +++ b/modification.m @@ -17,10 +17,11 @@ % cd mechanisms/; %% Settings %%%%%%%%%%%%%%%%%%%%%%%%%%% -fuel_name='n_heptane'; -mechanism='Ra_Reitz'; +fuel_name='n_dodecane'; +mechanism='MFC'; % class_numb =6; -numbOfClass = 9; +class_numb = [15]; +numbOfClass = length(class_numb); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %filename = [mechanism,'_',fuel_name,'_','class6.inp']; filename = [mechanism,'_','base.inp']; @@ -68,10 +69,10 @@ -for class_numb = 1 : numbOfClass +for k = 1 : numbOfClass rawCellColumns2=rawCellColumns; - class_numb_text{1} = ['class' num2str(class_numb)]; + class_numb_text{k} = ['class' num2str(class_numb(k))]; for variation_numb = 1:7; ModStart = 0; @@ -85,63 +86,95 @@ end %% classes modification - if class_numb == 1 & (strfind(rawCellColumns{i,j},'nc7h16+h=') == 1)... % n heptane erc - & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1)... - |(class_numb == 2 & (strfind(rawCellColumns{i,j},'nc7h16+oh=') == 1)... - & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... - |(class_numb == 3 & (strfind(rawCellColumns{i,j},'nc7h16+ho2=') == 1)... - & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... - |(class_numb == 4 & (strfind(rawCellColumns{i,j},'nc7h16+o2=') == 1)... - & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... - |(class_numb == 5 & (strfind(rawCellColumns{i,j},'c7h15-2+o2=') == 1)... - & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... - | class_numb == 6 & (strfind(rawCellColumns{i,j},'c7h15o2+o2=') == 1)... % n heptane erc - & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1)... - |(class_numb == 7 & (strfind(rawCellColumns{i,j},'c7ket12=') == 1)... - & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... - |(class_numb == 8 & (strfind(rawCellColumns{i,j},'c5h11co=') == 1)... - & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... - |(class_numb == 9 & (strfind(rawCellColumns{i,j},'c7h15-2=') == 1)... - & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% if class_numb == 1 & (strfind(rawCellColumns{i,j},'nc7h16+h=') == 1)... % n heptane erc +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1)... +% |(class_numb == 2 & (strfind(rawCellColumns{i,j},'nc7h16+oh=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb == 3 & (strfind(rawCellColumns{i,j},'nc7h16+ho2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb == 4 & (strfind(rawCellColumns{i,j},'nc7h16+o2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb == 5 & (strfind(rawCellColumns{i,j},'c7h15-2+o2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% | class_numb == 6 & (strfind(rawCellColumns{i,j},'c7h15o2+o2=') == 1)... % n heptane erc +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1)... +% |(class_numb == 7 & (strfind(rawCellColumns{i,j},'c7ket12=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb == 8 & (strfind(rawCellColumns{i,j},'c5h11co=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb == 9 & (strfind(rawCellColumns{i,j},'c7h15-2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... -% | class_numb == 27 & (strfind(rawCellColumns{i,j},'C12OOH') == 1)... % n dodecane ske_361 -% & (isempty(strfind(rawCellColumns{i,j},'O2=C12KET')) == 0)... -% & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1); + if class_numb(k) == 27 & (strfind(rawCellColumns{i,j},'C12OOH') ~= 0)... % n dodecane ske_361 + & (isempty(strfind(rawCellColumns{i,j},'O2=C12KET')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 21 & (strfind(rawCellColumns{i,j},'NC12H26+H=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+H=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 22 & (strfind(rawCellColumns{i,j},'NC12H26+OH=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+OH=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 23 & (strfind(rawCellColumns{i,j},'NC12H26+O=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+O=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 24 & (strfind(rawCellColumns{i,j},'NC12H26+HO2=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+HO2=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 11 & (strfind(rawCellColumns{i,j},'C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'+O2=C12H25O2')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 28 & (strfind(rawCellColumns{i,j},'C12KET') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'=OH+')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 26 & (strfind(rawCellColumns{i,j},'C12OOH') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'O2=C12OOH')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 15 & (strfind(rawCellColumns{i,j},'C12H25O2') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'=C12OOH')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i,j},'C12OOH6-3')) == 1)... + & (isempty(strfind(rawCellColumns{i,j},'C12OOH6-9')) == 1)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + %% if variation_numb == 1 - rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb),'***','v_',num2str(variation_numb)]; + rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb(k)),'***','v_',num2str(variation_numb)]; rawCellColumns2{i,2} = str2num(rawCellColumns{i,2}); % same lsqfit(variation_numb,1) = str2num(rawCellColumns{i,2}); lsqfit(variation_numb,2) = str2num(rawCellColumns{i,4}); % lsqfit(variation_numb+1,1) = str2num(rawCellColumns{i,2}); % lsqfit(variation_numb+1,2) = str2num(rawCellColumns{i,4}); elseif variation_numb == 2 - rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb),'***','v_',num2str(variation_numb)]; + rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb(k)),'***','v_',num2str(variation_numb)]; rawCellColumns2{i,2} = str2num(rawCellColumns{i,2})*2; % A*2 lsqfit(variation_numb,1) = str2num(rawCellColumns{i,2})*2; lsqfit(variation_numb,2) = str2num(rawCellColumns{i,4}); % lsqfit(variation_numb+1,1) = str2num(rawCellColumns{i,2})*2; % lsqfit(variation_numb+1,2) = str2num(rawCellColumns{i,4}); elseif variation_numb == 3 - rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb),'***','v_',num2str(variation_numb)]; + rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb(k)),'***','v_',num2str(variation_numb)]; rawCellColumns2{i,2} = str2num(rawCellColumns{i,2})*0.5; % A*0.5 lsqfit(variation_numb,1) = str2num(rawCellColumns{i,2})*0.5; lsqfit(variation_numb,2) = str2num(rawCellColumns{i,4}); elseif variation_numb == 4 - rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb),'***','v_',num2str(variation_numb)]; + rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb(k)),'***','v_',num2str(variation_numb)]; rawCellColumns2{i,2} = str2num(rawCellColumns{i,2})*7.5; % A>>1, A*7.5 rawCellColumns2{i,4} = str2num(rawCellColumns{i,4})+2000; % E+2000 lsqfit(variation_numb,1) = str2num(rawCellColumns{i,2})*7.5; lsqfit(variation_numb,2) = str2num(rawCellColumns{i,4})+2000; elseif variation_numb == 5 - rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb),'***','v_',num2str(variation_numb)]; + rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb(k)),'***','v_',num2str(variation_numb)]; rawCellColumns2{i,2} = str2num(rawCellColumns{i,2})*1.6;% A>1 rawCellColumns2{i,4} = str2num(rawCellColumns{i,4})+2000; % E+2000 lsqfit(variation_numb,1) = str2num(rawCellColumns{i,2})*1.6;% A>1 - lsqfit(variation_numb,2) = str2num(rawCellColumns{i,4})+2000; % E+2000 - + lsqfit(variation_numb,2) = str2num(rawCellColumns{i,4})+2000; % E+2000 elseif variation_numb == 6 - rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb),'***','v_',num2str(variation_numb)]; + rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb(k)),'***','v_',num2str(variation_numb)]; rawCellColumns2{i,2} = str2num(rawCellColumns{i,2})*0.6; % A<1 lsqfit(variation_numb,1) = str2num(rawCellColumns{i,2})*0.6; if str2num(rawCellColumns{i,4})-2000>0; @@ -152,7 +185,7 @@ lsqfit(variation_numb,2) = 0; end elseif variation_numb == 7 - rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb),'***','v_',num2str(variation_numb)]; + rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb(k)),'***','v_',num2str(variation_numb)]; rawCellColumns2{i,2} = str2num(rawCellColumns{i,2})*0.13; % A<<1 lsqfit(variation_numb,1) = str2num(rawCellColumns{i,2})*0.13; if str2num(rawCellColumns{i,4})-2000>0; @@ -168,7 +201,7 @@ %%saving files % save mech_ERC-MultiChem+Bio_Brakora2012_v1_test.mat rawCellColumns2 % file_name=[mechanism,'_',fuel_name,'_','class',class_numb,'_','v_',num2str(variation_numb),'.inp']; - file_name=[mechanism,'_',fuel_name,'_','class',num2str(class_numb),'_','v_',num2str(variation_numb),'.inp']; + file_name=[mechanism,'_',fuel_name,'_','class',num2str(class_numb(k)),'_','v_',num2str(variation_numb),'.inp']; % fileID = fopen('mech_ERC-MultiChem+Bio_Brakora2012_v1_test.inp','w'); fileID = fopen(file_name,'w'); output = cell(size(rawCellColumns2,1),size(rawCellColumns2,2)); @@ -219,7 +252,7 @@ end % clearvars filename delimiter formatSpec fileID dataArray ans raw col numericData rawCellColumns2 rawNumericColumns R; - rateParam.(class_numb_text{1})=lsqfit; + rateParam.(class_numb_text{k})=lsqfit; save('rateParam.mat','rateParam'); clear rawCellColumns2; end diff --git a/modification_afteroptimize_github.m b/modification_afteroptimize_github.m new file mode 100644 index 0000000..36557a1 --- /dev/null +++ b/modification_afteroptimize_github.m @@ -0,0 +1,222 @@ + +clear all; +%% Import data from text file. +% Script for importing data from the following text file: +% +% /home/jordan/Documents/Research/matlabfolder/mech.dat +% +% To extend the code to different selected data or a different text file, +% generate a function instead of a script. + +% Auto-generated by MATLAB on 2016/08/24 17:10:49 + +%% Initialize variables. +% filename = '/scratch/engin_flux/unghee/chemkin/mechanisms/mech_ERC-MultiChem+Bio_Brakora2012.inp'; +currentFolder = pwd; +cd(pwd); +% cd mechanisms/; +% filename = 'mech_ERC-MultiChem+Bio_Brakora2012_base.inp'; +filename = 'MFC_base.inp'; +delimiter = ' '; + +%% Read columns of data as strings: +% For more information, see the TEXTSCAN documentation. +% formatSpec = '%s%s%s%s%s%s%s%s%s%[^\n\r]'; +formatSpec = '%s%s%s%s%s%s%s%s%s%s%s%[^\n\r]'; +%% Open the text file. +fileID = fopen(filename,'r'); + +%% Read columns of data according to format string. +% This call is based on the structure of the file used to generate this +% code. If an error occurs for a different file, try regenerating the code +% from the Import Tool. +dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'ReturnOnError', false); + +%% Close the text file. +fclose(fileID); + +%% Convert the contents of columns containing numeric strings to numbers. +% Replace non-numeric strings with NaN. +raw = repmat({''},length(dataArray{1}),length(dataArray)-1); +for col=1:length(dataArray)-1 + raw(1:length(dataArray{col}),col) = dataArray{col}; +end +numericData = NaN(size(dataArray{1},1),size(dataArray,2)); + + +%% Split data into numeric and cell columns. +rawNumericColumns = {}; +rawCellColumns = raw(:, [1,2,3,4,5,6,7,8,9]); + + + %% Replace non-numeric cells with NaN + +[A,B]= find(ismember(rawCellColumns,'ELEMENTS')); + +locationClass = [A B]; + +[rowCellNumber,] = size(rawCellColumns); + +rawCellColumns2=rawCellColumns; + +ModStart = 0; +%%%%%%%%%%%%%%%%%%%%%%%% setting +mechanism='MFC'; +fuel_name='n_dodecane'; +date = '04_29_2017'; +iteration_numb = 1; +class_numb=[11 15 22 24 26 27 28]; % classes which have been optimized : this value should be changed manually by the user +pressure=[20]; +equi=1; +%%%%%%%%%%%%%%%%%%%%%%%% +numbOfClass = length(class_numb); +class_numb_text = {}; +for k=1:numbOfClass + class_numb_text=[class_numb_text ['class',num2str(class_numb(k))]]; +end + +for k=1:length(pressure) + pressure_text{k}=[num2str(pressure(k)),'atm']; + +end +load final_result.mat + + +for i=1 : rowCellNumber + [ ~, columnCellNumber] = size(rawCellColumns(i,:)); +% for j= 1 : cellfun('length',rawCellColumns{i,:}) + for j= 1 : 1: columnCellNumber + + if strmatch('reactions',rawCellColumns{i,j})==1; + ModStart = 1; + end + + + for k = 1: numbOfClass +%% User should add if condition in order to add more classes. +% if class_numb(k) == 1 & (strfind(rawCellColumns{i,j},'nc7h16+h=') == 1)... % n heptane erc +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1)... +% |(class_numb(k) == 2 & (strfind(rawCellColumns{i,j},'nc7h16+oh=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 3 & (strfind(rawCellColumns{i,j},'nc7h16+ho2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 4 & (strfind(rawCellColumns{i,j},'nc7h16+o2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 5 & (strfind(rawCellColumns{i,j},'c7h15-2+o2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% | class_numb(k) == 6 & (strfind(rawCellColumns{i,j},'c7h15o2+o2=') == 1)... % n heptane erc +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1)... +% |(class_numb(k) == 7 & (strfind(rawCellColumns{i,j},'c7ket12=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 8 & (strfind(rawCellColumns{i,j},'c5h11co=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +% |(class_numb(k) == 9 & (strfind(rawCellColumns{i,j},'c7h15-2=') == 1)... +% & isempty((strfind(rawCellColumns{i},'!')) == 0)& (ModStart == 1))... +%% + if class_numb(k) == 27 & (strfind(rawCellColumns{i,j},'C12OOH') ~= 0)... % n dodecane ske_361 + & (isempty(strfind(rawCellColumns{i,j},'O2=C12KET')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 21 & (strfind(rawCellColumns{i,j},'NC12H26+H=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+H=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 22 & (strfind(rawCellColumns{i,j},'NC12H26+OH=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+OH=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 23 & (strfind(rawCellColumns{i,j},'NC12H26+O=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+O=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 24 & (strfind(rawCellColumns{i,j},'NC12H26+HO2=C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'NC12H26+HO2=C12H25')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 11 & (strfind(rawCellColumns{i,j},'C12H25') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'+O2=C12H25O2')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 28 & (strfind(rawCellColumns{i,j},'C12KET') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'=OH+')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 26 & (strfind(rawCellColumns{i,j},'C12OOH') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'O2=C12OOH')) == 0)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + |class_numb(k) == 15 & (strfind(rawCellColumns{i,j},'C12H25O2') ~= 0)... + & (isempty(strfind(rawCellColumns{i,j},'=C12OOH')) == 0)... + & (isempty(strfind(rawCellColumns{i,j},'-1')) == 1)... %except -1 reaction + & (isempty(strfind(rawCellColumns{i,j},'C12OOH6-3')) == 1)... + & (isempty(strfind(rawCellColumns{i,j},'C12OOH6-9')) == 1)... + & (isempty(strfind(rawCellColumns{i},'!')) == 1)& (ModStart == 1)... + + rawCellColumns2{i,9} = ['!***ClASS',num2str(class_numb(k)),'***','after Optimize']; + + rawCellColumns2{i,2} = final_result.(class_numb_text{k})(2,1); % Optimized A + rawCellColumns2{i,4} = final_result.(class_numb_text{k})(2,2); % Optimized E + + + end + end + end +end + +%% + + + +save_filename='output'; + +for k = 1 :length(pressure) +directory=['C:\Users\unghee\Dropbox\post_process','\',mechanism,'\',fuel_name,'_',pressure_text{k},'_','phi',num2str(equi),'_',date]; + +cd(directory) + +mkdir(save_filename); +cd(save_filename); + +class_numb_text = []; +for k=1:numbOfClass + + class_numb_text=[class_numb_text '_' ['class',num2str(class_numb(k))]]; +end + +file_name=[mechanism,'_','afteroptimize','_',num2str(iteration_numb),'_','iteration','.inp']; +fileID = fopen(file_name,'w'); +output = cell(size(rawCellColumns2,1),size(rawCellColumns2,2)); + + + + for i = 1:size(rawCellColumns2,1) + for j = 1:size(rawCellColumns2,2) + if numel(rawCellColumns2{i,j}) == 0 + output{i,j} = ''; + % Check whether the content of cell i,j is + % numeric and convert numbers to strings. + elseif isnumeric(rawCellColumns2{i,j}) || islogical(rawCellColumns2{i,j}) + output{i,j} = num2str(rawCellColumns2{i,j}(1,1)); + + % If the cell already contains a string, nothing has to be done. + elseif ischar(rawCellColumns2{i,j}) + output{i,j} = rawCellColumns2{i,j}; + end; + + % Cell i,j is written to the output file. A delimiter is appended for + % all but the last element of each row. At the end of a row, a newline + % is written to the output file. + if j < size(rawCellColumns2,2) + fprintf(fileID,['%s',delimiter],output{i,j}); + else + fprintf(fileID,'%s\r\n',output{i,j}); + end + end; + end; + + +fclose(fileID); + +end +% exit %this should be uncommented when running in FLUX environment + + +%% Clear temporary variables +% clearvars filename delimiter formatSpec fileID dataArray ans raw col numericData rawCellColumns rawNumericColumns R; \ No newline at end of file diff --git a/optimization_7_repro_multiple_2.m b/optimization_7_repro_multiple_2.m index 241c849..f1a620a 100644 --- a/optimization_7_repro_multiple_2.m +++ b/optimization_7_repro_multiple_2.m @@ -1,206 +1,226 @@ -clear all -% close all +% clear all % comment out when you want to test optimizer alone +% close all % comment out when you want to test optimizer alone - -%test -classnumb=[4 6]; +% classnumb=[11 15 22 24 26 27 28]; % comment out when you want to test optimizer alone numbOfClass = length(classnumb); class_numb_text = {}; for k=1:numbOfClass -% classnumb_text{classnumb(k)}=['class',num2str(classnumb(k))]; class_numb_text=[class_numb_text ['class',num2str(classnumb(k))]]; end -fuel_sim={'modify_class4_class6_success'}; +fuel_sim={'modify'}; %% read modification ignition delay time -% mechanism={'MFC'}; -mechanism={'Ra_Reitz'}; -% currentloc = pwd; -date = {'01_30_2017'}; +% mechanism={'MFC'}; %comment out when you want to test optimizer alone +mechanism={'MFC'}; +%comment out when you want to test optimizer alone +% date = {'03_16_2017_1_iteration'}; % fuel_name = {'n_dodecane'}; -fuel_name = {'n_heptane'}; -equi=1; +% fuel_name = {'n_heptane'}; +% equi=1; currentloc = 'C:\Users\unghee\Dropbox\post_process'; - - - -pressure=[20 40]; +% pressure=[20]; for k=1:length(pressure) pressure_text{k}=[num2str(pressure(k)),'atm']; end -k = 2 -directory=[fuel_name{1},'_',pressure_text{k},'_','phi',num2str(equi),'_',date{1}]; -% location_rateParam = [currentloc,'\',mechanism{1},'\',fuel_name] -location_rateParam=[currentloc,'\',mechanism{1},'\',directory,'\',fuel_sim{1},'\',class_numb_text{1}]; +m = 1; +directory=[fuel_name{1},'_',pressure_text{m},'_','phi',num2str(equi),'_',date{1}]; +location_rateParam=[currentloc,'\',mechanism{1},'\',directory,'\',fuel_sim{1}]; cd(location_rateParam) load rateParam.mat; +clear m; +%% exp data : setting target for optimization. - -% load rateParam.mat; - +% addpath('C:\Users\unghee\Dropbox\post_process'); % real_fuel_ID; % pure_component_ID; -% Target_data = Shen_hep_40atm; -% % Target_data = Vasu_40atm; -% % Target_data = Wang_40atm; -% Temp = Target_data(:,2); -% Temp_un = Target_data(:,1); - -num_cases_target=25; -location_target=[currentloc,'\',mechanism{1},'\',directory,'\','beforeoptimize']; -cd(location_target) -addpath(currentloc) -time_struct_target=read_ignition_delay(location_target,num_cases_target); -Temp=time_struct_target.table.data(:,6)/1000; % do we have to include temp? - -Target_data=time_struct_target.table.data(:,10); -range = 13:20; -Temp = Temp(range); -Target_data = Target_data(range); - +% Target_fuel1 = Vasu_dode_20atm; % n_dodecane optimization at 20atm. +% Target_data1=Target_fuel1(:,5); +% Temp1 = Target_fuel1(:,2); +% numbOftarget1 =length(Target_data1); + + +% Target_fuel2 = Shen_dode_40atm; % n_dodecane optimization at 40atm. +% Target_data2=Target_fuel2(:,5); +% Temp2 = Target_fuel2(:,2); +% numbOftarget2 =length(Target_data2); + +% Temp = [Temp1; Temp2]; +% Target_data = [Target_data1; Target_data2]; +%% only optimizing one pressure condition : comment out when optimizing multiple pressure condition. +numbOftarget2=0; +Temp = Temp1; +Target_data = Target_data1; +%% cd ../../.. - -pressure=[20 40]; -for k=1:length(pressure) - pressure_text{k}=[num2str(pressure(k)),'atm']; - -end -equi=1; - -% class_numb_text{1} = ['class6']; - A = rateParam.(class_numb_text{1})(:,1); E = rateParam.(class_numb_text{1})(:,2); %% read modification ignition delay time -mechanism={'Ra_Reitz'}; +% loading ignition delay times with 7 variations and save into structure format +% this will be used in the regression model part + currentloc = 'C:\Users\unghee\Dropbox\post_process'; num_cases_modification= size(A,1); -date = {'01_30_2017'}; -fuel_name = {'n_heptane'}; -directory=[fuel_name{1},'_',pressure_text{k},'_','phi',num2str(equi),'_',date{1}]; -for j = 1: size(Temp,1) + +m=1; % pressure 20atm +directory=[fuel_name{1},'_',pressure_text{m},'_','phi',num2str(equi),'_',date{1}]; + +for j = 1: numbOftarget1 for k = 1: numbOfClass -% location_modification=[currentloc,'\',mechanism{1},'_','modify','_',date{1},'\',num2str(Temp_un(j))]; + location_modification=[currentloc,'\',mechanism{1},'\',directory,'\',fuel_sim{1},'\',class_numb_text{k},'\',num2str(j)]; time_struct_modification=read_ignition_delay(location_modification,num_cases_modification); - time_modification(:,j)=time_struct_modification.table.data(:,10); + time_modification.(class_numb_text{k})(:,j)=time_struct_modification.table.data(:,10); + temp_modification.(class_numb_text{k})(:,j)=time_struct_modification.table.data(:,6); end end +for k = 1: numbOfClass + time_modification.(class_numb_text{k})=time_modification.(class_numb_text{k})(:,1: numbOftarget1); + temp_modification.(class_numb_text{k})=temp_modification.(class_numb_text{k})(:,1: numbOftarget1); +end -%% coefficient -for k = 1 : size(Temp,1) -% for k = 1 : 1 - - Temp_current = Temp(k); - time_current = time_modification(:,k); -% type = fittype({'log(x)','log(Temp_current)','y','log(x)*y','(log(x))^2','y^2'}... % linear term3 -% ,'independent', {'x', 'y'},'dependent','z','problem','Temp_current'); -% options=fitoptions('Method','LinearLeastSquares') -% sf = fit([A,E],log(time_current),type,options, 'problem',Temp_current); +%% comment out for only one pressure condition optimization. +% m =2 ; % pressure 40atm +% directory=[fuel_name{1},'_',pressure_text{m},'_','phi',num2str(equi),'_',date{1}]; +% +% for j = 1: numbOftarget2 +% for k = 1: numbOfClass % -% coefs(k,:)=coeffvalues(sf); -% plot(sf,[A,E],log(time_current)); -% zlim([4 9]); - totalM = []; - for n = 1: numbOfClass - A = rateParam.(class_numb_text{n})(:,1); - E = rateParam.(class_numb_text{n})(:,2); - M = [log(A) log(Temp_current)*ones(7,1) E log(A).*(E) log(A).*log(A) E.^2]; - totalM = [totalM M]; +% location_modification=[currentloc,'\',mechanism{1},'\',directory,'\',fuel_sim{1},'\',class_numb_text{k},'\',num2str(j)]; +% time_struct_modification=read_ignition_delay(location_modification,num_cases_modification); +% time_modification.(class_numb_text{k})(:,numbOftarget1+j)=time_struct_modification.table.data(:,10); +% temp_modification.(class_numb_text{k})(:,numbOftarget1+j)=time_struct_modification.table.data(:,6); +% end +% +% end - end + +%% coefficient : defining regression model. +for j = 1 : size(Temp,1) + + Temp_current = Temp(j); % defines regression model for individual temp. + + + for k = 1: numbOfClass + time_current = time_modification.(class_numb_text{k})(:,j); + A = rateParam.(class_numb_text{k})(:,1); % preexponential + E = rateParam.(class_numb_text{k})(:,2); % activiation energy + M = [log(A) log(Temp_current)*ones(7,1) E log(A).*(E) log(A).*log(A) E.^2]; % form a matrix composed of 7 variation preexponential and activiation energy d = log(time_current); - coefs_inv = lsqlin(totalM,d); + coefs_inv = lsqlin(M,d); % by least square, define the coefficient of regression model coefs_element = coefs_inv'; - coefsTotal = []; - - for n = 1: length(class_numb_text) - coefs.(class_numb_text{n})(k,:)=coefs_element(1,6*n-5:6*n); - coefsTotal =[coefsTotal coefs.(class_numb_text{n})(k,:)]; - end + + coefs.(class_numb_text{k})(j,:)=coefs_element; % save into structure file + + prediction.(class_numb_text{k}){:,j}=M*coefs_element'-d; % check whether regression model matches with the target + +%% commented out for plotting regression. Does not affect the optimization calculation +% predictionreg.(class_numb_text{k}){:,j}=M*coefs_element'; +% plotregression(d,M*coefs_element') %plot the regression + end - prediction{:,k}=totalM*coefsTotal'-d; - -% Model_igtime=@(x,y) coefs(k,1)*log(x)+coefs(k,2)*log(Temp_current)+coefs(k,3)*y... -% + coefs(k,4)*log(x)*y + coefs(k,5)*(log(x))^2+coefs(k,6)*y^2; % linear term3 -% -% for i= 1: size(A,1) -% prediction{:,k}(i,1)=exp(Model_igtime(A(i),E(i))); -% end -% prediction{:,k}(:,2)=time_current; -% prediction{:,k}(:,3)= prediction{:,k}(:,1)-prediction{:,k}(:,2); -% prediction{:,k}(:,4)= mean(abs(prediction{:,k}(:,3))); end +%% plot regression +% for k = 1: numbOfClass +% plotregression(log(time_modification.(class_numb_text{k})(:,1:10)),predictionreg.(class_numb_text{k})(1:10)) +% end -% coefs=coefs.class4; %% OPTIMIZER -% Target Temp value - -% te(:,1) = Target_data(:,3); - -te(:,1) = Target_data; -numberOftempPoints = size(te(:,1),1); +numberOftempPoints = size(Target_data,1); %% objective function -% ObjectiveFunction=ObjectiveFunction_array5(numberOftempPoints,numbOfClass,te,Temp,coefs); -% ObjectiveFunction=ObjectiveFunction_array6(numberOftempPoints,numbOfClass,te,Temp,coefs,class_numb_text); -ObjectiveFunction = @(X) find_rate_3(X,coefs,Temp,numbOfClass,class_numb_text,Target_data); - +ObjectiveFunction = @(X) find_rate_weighting(X,coefs,Temp,numbOfClass,class_numb_text,Target_data,numbOftarget1,numbOftarget2); LB =[]; UB =[]; +%% setting boundry values +% If User want to add more classes, they have to include additional classes +% below. Currently the boundary values are set up for erc bio fuel +% mechansim and 361kp skeletal mechanism + +% erc bio fuel mechanism +% class2 +if strcmp(mechanism{1},'Ra_Reitz') && ismember(2,classnumb) +LB =[LB rateParam.('class2')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class2')(1,1)*10 rateParam.('class2')(1,2)+2000 ]; +end % class4 if strcmp(mechanism{1},'Ra_Reitz') && ismember(4,classnumb) -% LB =[0.706e+14*0.13 37904-2000 ]; -% LB =[0.706e+14*0.13 0 ]; -% UB =[0.706e+14*10 37904+2000]; -LB =[42360000000000*0.13 0 ]; -UB =[42360000000000*10 35904+2000]; +LB =[LB rateParam.('class4')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class4')(1,1)*10 rateParam.('class4')(1,2)+2000 ]; end % class6 if strcmp(mechanism{1},'Ra_Reitz') && ismember(6,classnumb) -LB =[LB 19955000000000*0.13 0]; -UB =[UB 19955000000000*10 16232.712+2000]; +LB =[LB rateParam.('class6')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class6')(1,1)*10 rateParam.('class6')(1,2)+2000 ]; +end + +% ndodecane 361kp skeletal mech +if strcmp(mechanism{1},'MFC') && ismember(11,classnumb) +LB =[LB rateParam.('class11')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class11')(1,1)*10 rateParam.('class11')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(15,classnumb) +LB =[LB rateParam.('class15')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class15')(1,1)*10 rateParam.('class15')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(21,classnumb) +LB =[LB rateParam.('class21')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class21')(1,1)*10 rateParam.('class21')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(22,classnumb) +LB =[LB rateParam.('class22')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class22')(1,1)*10 rateParam.('class22')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(23,classnumb) +LB =[LB rateParam.('class23')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class23')(1,1)*10 rateParam.('class23')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(24,classnumb) +LB =[LB rateParam.('class24')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class24')(1,1)*10 rateParam.('class24')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(26,classnumb) +LB =[LB rateParam.('class26')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class26')(1,1)*10 rateParam.('class26')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(27,classnumb) +LB =[LB rateParam.('class27')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class27')(1,1)*10 rateParam.('class27')(1,2)+2000 ]; +end +if strcmp(mechanism{1},'MFC') && ismember(28,classnumb) +LB =[LB rateParam.('class28')(1,1)*0.13 0 ]; +UB =[UB rateParam.('class28')(1,1)*10 rateParam.('class28')(1,2)+2000 ]; end nvars=2*numbOfClass; options=gaoptimset('PopulationSize',500); -[result_ga,Fval,exitFlag,Output] = ga(ObjectiveFunction,nvars,[],[],[],[],LB,UB,[],options); -[result_fmin,Fval,exitFlag,Output] = fmincon(ObjectiveFunction,result_ga,[],[],[],[],LB,UB); +[result_ga,Fval,exitFlag,Output] = ga(ObjectiveFunction,nvars,[],[],[],[],LB,UB,[],options); % save optimized rate params from ga +[result_fmin,Fval,exitFlag,Output] = fmincon(ObjectiveFunction,result_ga,[],[],[],[],LB,UB); % save optimized rate params from ga & fmincon X = result_ga; -error=ObjectiveFunction(X) +error=ObjectiveFunction(X) % outputs the error of the objective function. This shows how the predicted values deviate from the desired(target) value for i = 1 : numbOfClass -final_result.(class_numb_text{i})= [result_ga(1,2*i-1:2*i); result_fmin(1,2*i-1:2*i)]; +final_result.(class_numb_text{i})= [result_ga(1,2*i-1:2*i); result_fmin(1,2*i-1:2*i)]; % saves optimized value from ga,fmincon end location_save=[currentloc,'\',mechanism{1},'\',directory]; cd(location_save) +cd ../ save('final_result.mat','final_result') -%% read final ignition delay time after optimization plug into plot part. -% mechanism={'Ra_Reitz'}; -% fuel_sim={'afteroptimize_class6_temp_all_shen' }; -% pressure = 40; -% num_cases=25; -% currentloc = pwd; -% location=[currentloc,'\',mechanism{i},'\',fuel_sim{1},'_',num2str(pressure(1)),'atm_phi1']; -% -% % -% time_struct=read_ignition_delay(location,num_cases); -% time=time_struct.table.data(:,10)/1000; + diff --git a/optimization_driver_2.m b/optimization_driver_2.m new file mode 100644 index 0000000..b8fa069 --- /dev/null +++ b/optimization_driver_2.m @@ -0,0 +1,55 @@ +clear; + +%% initial setting +mechanism={'MFC'}; +fuel_name={'n_dodecane'}; +date = {'04_29_2017'}; +fuel_sim={'modify_sensitivity_0_iteration'}; % change the number of iteration +equi = 1; +classnumb=[11 15 21 22 23 24 26 27 28]; % this value can be changed based on what mechanism use. Currently, it is set up for 361kp skeletal mech. +pressure=[20]; % add more pressure for mutiple condition optimization +%% Target setting +addpath('C:\Users\unghee\Dropbox\post_process'); +real_fuel_ID; +pure_component_ID; +Target_fuel1 = Vasu_dode_20atm; % n dodecane target value at 20 atm + +Target_data1=Target_fuel1(:,5); +numbOftarget1=length(Target_data1); +Temp1=Target_fuel1(:,2); + +Target_fuel2 = Shen_dode_40atm; % n dodecane target value at 40 atm +Target_data2=Target_fuel2(:,5); +Temp2 = Target_fuel2(:,2); +numbOftarget2=length(Target_data2); + +Temp = [Temp1; Temp2]; +Target_data = [Target_data1; Target_data2]; + +%% sensitivity +addpath('C:\Users\unghee\Dropbox\JPoptimization_github'); +sensitivity_analysis_3; % executes sensitivity analysis + +m =1; +disp('looking at sensitivity at one pressure condition') +class_to_optimize.(pressure_text{m}) = []; % saves sensitivities for pressures. +for k=1:length(classnumb) + + if (sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig_avg > 1)... % thershold for translational sensitivity, user may change based on their needs + && (sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr_avg >= 0.38) % thershold for rotational sensitivity, user may change based on their needs + class_to_optimize.(pressure_text{m}) = [class_to_optimize.(pressure_text{m}) classnumb(k)]; + end +end + +classnumb = class_to_optimize.(pressure_text{m}) % classes to optimize after sensitivity analysis + + + + + +%% Optimization +% classnumb= [26 27] +clearvars -except classnumb Temp Temp1 Temp2 numberOftarget Target_data Target_data1 numbOftarget1 sensitivity pressure fuel_name equi fuel_name mechanism fuel_sim date +optimization_7_repro_multiple_2 % executes mechanism optimizer + + diff --git a/optimization_driver_2_lockdown.m b/optimization_driver_2_lockdown.m new file mode 100644 index 0000000..a87b209 --- /dev/null +++ b/optimization_driver_2_lockdown.m @@ -0,0 +1,111 @@ +clear; +close all; +mechanism={'Ra_Reitz'}; +fuel_name={'n_heptane'}; +date = {'01_30_2017'}; +fuel_sim={'modify_class2_class4_class6_moderate'}; +equi = 1; +pressure=[40]; +classnumb=[2 4 6]; + +sensitivity_analysis_3; + +clearvars -except sensitivity classnumb classnumb_text pressure_text pressure mechanism fuel_name date fuel_sim equi +% Determine classes with overall sensitivity +classnumb_optimize.entire.(pressure_text{1}) = []; +% classnumb_optimize.entire.(pressure_text{1}) = sturct([]); +pressure_optimize_text =[]; +pressure_optimize =[]; +for m = 1 : length(pressure) + for k=1:length(classnumb) + + if (sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig_avg > 1) && (sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr_avg > 0.5) + classnumb_optimize.entire.(pressure_text{m}) = [ classnumb_optimize.entire.(pressure_text{m}) classnumb(k)]; + + end + end + pressure_optimize_text = [pressure_optimize pressure_text{m}]; + pressure_optimize = [pressure_optimize pressure(m)]; +end + + for m = 1 : length(pressure_optimize) + classnumb=classnumb_optimize.entire.(pressure_text{m}); %have to change for multiple pressure in optimization.m + end +%% HTR optimize +% add sgr information for the future work +classnumb_optimize.HTR.(pressure_text{1}) = []; +pressure_optimize_text =[]; +pressure_optimize =[]; +for m = 1 : length(pressure) + for k=1:length(classnumb) + + if (abs(sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig(4))) > 1 + classnumb_optimize.HTR.(pressure_text{m}) = [classnumb_optimize.HTR.(pressure_text{m}) classnumb(k)]; + end + end + pressure_optimize_text = [pressure_optimize pressure_text{m}]; + pressure_optimize = [pressure_optimize pressure(m)]; +end +for m = 1 : length(pressure_optimize) +classnumb=classnumb_optimize.HTR.(pressure_text{m}); +end +pressure = pressure_optimize; + +clearvars -except classnumb pressure mechanism fuel_name date fuel_sim equi classnumb_optimize sensitivity +range = 1:10; +% optimization_10; +optimization_7_repro_multiple_2; +%% +% final result for HTR, and LockDouwn +clearvars -except final_result pressure mechanism fuel_name date fuel_sim equi classnumb_optimize sensitivity pressure_text + +numbOfPressure=length(pressure) ; +for k=1:length(pressure) + pressure_text{k}=['P',num2str(pressure(k)),'atm']; +end + +classnumb_optimize.NTC.(pressure_text{1})=[]; +pressure_optimize_text =[]; +pressure_optimize =[]; +Entire = classnumb_optimize.entire.(pressure_text{1}); +HTR = classnumb_optimize.HTR.(pressure_text{1}); + +classnumb_optimize.NTC.(pressure_text{1})... + = setxor(Entire,HTR); +numbOfClass_NTC = length(classnumb_optimize.NTC.(pressure_text{1})); + +class_numb_text_NTC = {}; +for k=1:numbOfClass_NTC +% classnumb_text{classnumb(k)}=['class',num2str(classnumb(k))]; + class_numb_text_NTC=[class_numb_text_NTC ['class',num2str(classnumb_optimize.NTC.(pressure_text{1})(k))]]; +end + +classnumbOfNTC = []; +for m = 1 : length(pressure) + for k=1:numbOfClass_NTC + + if (abs(sensitivity.(pressure_text{m}).(class_numb_text_NTC{k}).Sig(2)) > 1) + classnumbOfNTC = [classnumbOfNTC classnumb_optimize.NTC.(pressure_text{1})(k)]; + end + end + pressure_optimize_text = [pressure_optimize pressure_text{m}]; + pressure_optimize = [pressure_optimize pressure(m)]; +end +classnumb_optimize.NTC.(pressure_text{m}) = classnumbOfNTC; + +for m = 1 : length(pressure_optimize) +classnumb=classnumb_optimize.NTC.(pressure_text{m}); +end +pressure = pressure_optimize; + + +clearvars -except final_result classnumb pressure mechanism fuel_name date fuel_sim equi classnumb_optimize sensitivity +range = 12:20; + +optimization_7_repro_multiple_2; + +%modification after optimize + +%load + +% load -> output rightaway bash \ No newline at end of file diff --git a/sensitivity_analysis_3.m b/sensitivity_analysis_3.m new file mode 100644 index 0000000..32de71b --- /dev/null +++ b/sensitivity_analysis_3.m @@ -0,0 +1,164 @@ + + +% pressure=[20 40]; % comment out for not using driver file/ +numbOfPressure=length(pressure) ; +for k=1:length(pressure) + pressure_text{k}=['P',num2str(pressure(k)),'atm']; +end +% classnumb=[1 2 3 4 5 6 7 8 9]; % comment out for not using driver file +numbOfClass=length(classnumb) ; +for k=1:numbOfClass + classnumb_text{k}=['class',num2str(classnumb(k))]; +end + + +%%%%setting%%%%%%%%%%%%%%%%%%%%%%%%%%%% +exp_markers={'ro' 'g*' 'g^'}; +sim_line={'k-' 'k-s','k-d','k-^','k-x','k-o','k-p'}; +marker_size=8; +line_width=2; +% developer note +% currently it is set up to check the sensitivities at one pressure +% condition. Even performing the multiple pressure sensitivity analysis, it +% assumes that at both pressure condtions the sensitivities are similar. +% However, user may want to perform sensitivity for each pressure +% condition. The code works for it, but for the optimizier this is not +% implemented(you cannot send both sensitivities info to optimizer). and +% this does not make since, because there is only one mechanism optimized even for +% multiple pressure optimization. +% m=2 ; % pressure setting pressure for 40 atm if pressure = [20 40] + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% sensitivity analysis +% mechanism={'Ra_Reitz'}; +% fuel_sim={'modify_sensitivity'}; +% fuel_name={'n_heptane'} +% equi = 1; +currentloc = 'C:\Users\unghee\Dropbox\post_process'; +num_cases_modification= 3; % unlike 7 variations, sensitivity analysis only need 3 variations, which is base, k=2, k =0.5 + + +for m = 1 : numbOfPressure + + for k = 1 :numbOfClass + for j = 1: 25 + location_modification=['C:\Users\unghee\Dropbox\post_process','\',mechanism{1},'\',fuel_name{1},'_',num2str(pressure(m)),... + 'atm','_','phi',num2str(equi),'_',date{1},'\',fuel_sim{1},'\',classnumb_text{k},'\',num2str(j)]; + time_struct_modification=read_ignition_delay(location_modification,num_cases_modification); % load ignition time variation + sensitivity.(pressure_text{m}).(classnumb_text{k}).data(:,j)=time_struct_modification.table.data(:,10); + + end + + + temp=700:25:1300; + h =25; + tbase = sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,:); % when base + tk1 = sensitivity.(pressure_text{m}).(classnumb_text{k}).data(2,:); % when k = 2 + tk2 = sensitivity.(pressure_text{m}).(classnumb_text{k}).data(3,:); % when k = 0.5 + d_log_tk1=diff(log10(tk1)); % difference of ignition time at proximate temps when k =2 + d_log_tk2=diff(log10(tk2)); % difference of ignition time at proximate temps when k =0.5 + dh = log(h); + k1 = 2; + k2 = 0.5; + + [peaks,locs] = findpeaks(tbase/1000); % find two peak temp +% yvalue=tbase/1000; % when it is milisecond + locsValue = temp(locs); + [peaks2,locs2] = findpeaks(-tbase/1000); % find two peak temp +% yvalue=tbase/1000; % when it is milisecond + locsValue2 = temp(locs2); + + p_tempPoints_ig =[1, locs2, locs, length(temp)]; % find 4 temp points(2peak 2end point) for sensitivity analysis + temp_p_tempPoints_ig=[temp(1),temp(locs2),temp(locs), temp(length(temp))]; + +%translational sensitivity +% for definition, refer to Ra&Reitz,2011 CNF paper + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig = ones(1,length(p_tempPoints_ig)); + for i = 1 : length(p_tempPoints_ig) + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig(i)... + = 100*(log10(tk1(p_tempPoints_ig(i))) - log10(tk2(p_tempPoints_ig(i))))/(log10(tbase(i))*log10(k1/k2)); % translational shift + end + + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig_avg... + = mean(abs(sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig)); + + + p_tempPoints_gr =[1, locs2, locs, length(temp)]; % added end point + + temp_p_tempPoints_gr=[temp(locs2),temp(locs), temp(length(temp))]; + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr = ones(1,length(p_tempPoints_gr)-1); +%rotational sensitivity + for i = 1 : length(p_tempPoints_gr)-1 + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr(i)... + = 100*(d_log_tk1(p_tempPoints_gr(i))/dh - d_log_tk2(p_tempPoints_gr(i))/dh)/log10(k1/k2); + end + + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr_avg... + = mean(abs(sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr)); + end + +end + + +save('sensitivity.mat','sensitivity') + +%% plotting the sensitivities. Currently it is turned off for calculation speed +% for m = 1 : numbOfPressure +% for k = 1 : numbOfClass +% h=figure('position',[20 50 1200 480]); +% subplot(1,2,1); +% set(gca,'Fontsize',13) +% semilogy(temp,sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,:)/1000,'k-','markersize',marker_size); +% hold on +% semilogy(temp,sensitivity.(pressure_text{m}).(classnumb_text{k}).data(2,:)/1000,'rs-','markersize',marker_size); +% hold on +% semilogy(temp,sensitivity.(pressure_text{m}).(classnumb_text{k}).data(3,:)/1000,'bx-','markersize',marker_size); +% +% +% semilogy([temp_p_tempPoints_ig(1),temp_p_tempPoints_ig(end)],... +% [sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,1)/1000,... +% sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,end)/1000]... +% ,'gp','markersize',12); +% hold on +% semilogy(locsValue,peaks,'gp','markersize',12); +% hold on +% semilogy(locsValue2,-peaks2,'gp','markersize',12); +% +% legend('baseline','k=2','k=1','extreme points') +% xlabel('1000/T (1/K)') +% ylabel('Ignition Delay Time (ms)') +% +% +% +% +% subplot(1,2,2); +% plot(temp_p_tempPoints_ig,sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig,'k^-','markersize',marker_size); +% ylim([-25 25]) +% hold on +% % figure +% plot(temp_p_tempPoints_gr,sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr,'ko-','markersize',marker_size); +% +% ylim([-30 25]) +% xlabel('1000/T (1/K)') +% ylabel('Ignition Delay Time (ms)') +% legend('ignition delay sensitivity','gradient sensitivity') +% +% annotation(h,'textbox',[0.213 0.38 0.279 0.05],... +% 'String',{fuel_name{1},... +% '/Air',classnumb_text{k}, '\phi=1',pressure_text{m}},... +% 'FontSize',13,... +% 'FontName','Arial',... +% 'FitBoxToText','off',... +% 'LineStyle','none'); +% +% % location_save=[currentloc,'\',mechanism{1},'\',directory]; +% % cd(location_save) +% +% mkdir('sensitivity'); +% cd('sensitivity'); +% saveas(h,classnumb_text{k},'fig') +% saveas(h,classnumb_text{k},'jpg') +% cd ../ +% end +% end \ No newline at end of file diff --git a/sensitivity_analysis_3_modify.m b/sensitivity_analysis_3_modify.m new file mode 100644 index 0000000..aa1b5f6 --- /dev/null +++ b/sensitivity_analysis_3_modify.m @@ -0,0 +1,154 @@ + +%1. matlab peaks func -> temp ????? sensitivity data (Before,after 25temps) - 4?? + +% pressure=[20 40]; +numbOfPressure=length(pressure) ; +for k=1:length(pressure) + pressure_text{k}=['P',num2str(pressure(k)),'atm']; +end +% classnumb=[1 2 3 4 5 6 7 8 9]; +numbOfClass=length(classnumb) ; +for k=1:numbOfClass + classnumb_text{k}=['class',num2str(classnumb(k))]; +end + +exp_markers={'ro' 'g*' 'g^'}; +sim_line={'k-' 'k-s','k-d','k-^','k-x','k-o','k-p'}; +marker_size=8; +line_width=2; + + +%%%%setting%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % i=1 ; +% m=2 ; % pressure +% % j=1 ; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% sensitivity analysis +% mechanism={'Ra_Reitz'}; +% fuel_sim={'modify_sensitivity'}; +% fuel_name={'n_heptane'} +% equi = 1; +currentloc = 'C:\Users\unghee\Dropbox\post_process'; +num_cases_modification= 3; + + +for m = 1 : numbOfPressure + + for k = 1 :numbOfClass + for j = 1: numbOftarget1 + location_modification=['C:\Users\unghee\Dropbox\post_process','\',mechanism{1},'\',fuel_name{1},'_',num2str(pressure(m)),... + 'atm','_','phi',num2str(equi),'_',date{1},'\',fuel_sim{1},'\',classnumb_text{k},'\',num2str(j)]; + time_struct_modification=read_ignition_delay(location_modification,num_cases_modification); + sensitivity.(pressure_text{m}).(classnumb_text{k}).data(:,j)=time_struct_modification.table.data(:,10); + + end + + +% temp=700:25:1300; +% h =25; + temp=Temp1; + h = numbOftarget1; + + tbase = sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,:); % when base + tk1 = sensitivity.(pressure_text{m}).(classnumb_text{k}).data(2,:); % when k = 2 + tk2 = sensitivity.(pressure_text{m}).(classnumb_text{k}).data(3,:); % when k = 0.5 + d_log_tk1=diff(log10(tk1)); + d_log_tk2=diff(log10(tk2)); + dh = log(h); + k1 = 2; + k2 = 0.5; + + [peaks,locs] = findpeaks(tbase/1000); +% yvalue=tbase/1000; + locsValue = temp(locs); + [peaks2,locs2] = findpeaks(-tbase/1000); +% yvalue=tbase/1000; + locsValue2 = temp(locs2); + + p_tempPoints_ig =[1, locs2, locs, length(temp)]; + temp_p_tempPoints_ig=[temp(1),temp(locs2),temp(locs), temp(length(temp))]; + +%translational sensitivity +%change the normalize contant log10?e.x. log10^3 + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig = ones(1,length(p_tempPoints_ig)); + for i = 1 : length(p_tempPoints_ig) + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig(i)... + = 100*(log10(tk1(p_tempPoints_ig(i))) - log10(tk2(p_tempPoints_ig(i))))/(log10(tbase(i))*log10(k1/k2)); % translational shift + end + + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig_avg... + = mean(abs(sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig)); + + + p_tempPoints_gr =[1, locs2, locs, length(temp)]; % added end point + + temp_p_tempPoints_gr=[temp(locs2),temp(locs), temp(length(temp))]; + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr = ones(1,length(p_tempPoints_gr)-1); + %rotational sensitivity + for i = 1 : length(p_tempPoints_gr)-1 + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr(i)... + = 100*(d_log_tk1(p_tempPoints_gr(i))/dh - d_log_tk2(p_tempPoints_gr(i))/dh)/log10(k1/k2); + end + + sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr_avg... + = mean(abs(sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr)); + end + +end + + +save('sensitivity.mat','sensitivity') + +%% plotting +for m = 1 : numbOfPressure +for k = 1 : numbOfClass +h=figure('position',[20 50 1200 480]); +subplot(1,2,1); +set(gca,'Fontsize',13) +semilogy(temp,sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,:)/1000,'k-','markersize',marker_size); +hold on +semilogy(temp,sensitivity.(pressure_text{m}).(classnumb_text{k}).data(2,:)/1000,'rs-','markersize',marker_size); +hold on +semilogy(temp,sensitivity.(pressure_text{m}).(classnumb_text{k}).data(3,:)/1000,'bx-','markersize',marker_size); + + +semilogy([temp_p_tempPoints_ig(1),temp_p_tempPoints_ig(end)],... + [sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,1)/1000,... + sensitivity.(pressure_text{m}).(classnumb_text{k}).data(1,end)/1000]... + ,'gp','markersize',12); +hold on +semilogy(locsValue,peaks,'gp','markersize',12); +hold on +semilogy(locsValue2,-peaks2,'gp','markersize',12); + +legend('baseline','k=2','k=1','extreme points') +xlabel('1000/T (1/K)') +ylabel('Ignition Delay Time (ms)') + + + + +subplot(1,2,2); +plot(temp_p_tempPoints_ig,sensitivity.(pressure_text{m}).(classnumb_text{k}).Sig,'k^-','markersize',marker_size); +ylim([-25 25]) +hold on +% figure +plot(temp_p_tempPoints_gr,sensitivity.(pressure_text{m}).(classnumb_text{k}).Sgr,'ko-','markersize',marker_size); + +ylim([-30 25]) +xlabel('1000/T (1/K)') +ylabel('Ignition Delay Time (ms)') +legend('ignition delay sensitivity','gradient sensitivity') + +annotation(h,'textbox',[0.213 0.38 0.279 0.05],... + 'String',{fuel_name{1},... + '/Air',classnumb_text{k}, '\phi=1',pressure_text{m}},... + 'FontSize',13,... + 'FontName','Arial',... + 'FitBoxToText','off',... + 'LineStyle','none'); + + +end +end \ No newline at end of file