rundke

Script for running the DKE solver (can be modified by the user for specific simulations)
by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and Joan Decker MIT-RLE (jodecker@mit.edu)

0001 %Script for running the DKE solver (can be modified by the user for specific simulations)
0002 %by Y.Peysson CEA-DRFC <yves.peysson@cea.fr> and Joan Decker MIT-RLE (jodecker@mit.edu)
0003 %
0004 clear all
0005 clear mex
0006 clear functions
0007 close all
0008 warning off
0009 global nfig
0010 %
0011 permission = test_permissions_yp;
0012 %
0013 if ~permission 
0014     disp('Please move the script to a local folder where you have write permission before to run it')
0015     return;
0016 end
0017 %
0018 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0019 %
0020 id_simul = 'LH_karney_transp_bounce';%Simulation ID
0021 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0022 %
0023 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0024 rho_S = [0.05:0.02:0.2,0.22:0.03:0.6,0.62:0.05:0.95];%Normalized radial flux grid where calculations are performed (0 < rho_S < 1) (If one value: local calculation only, not used if empty)
0025 %
0026 id_path = '';%For all paths used by DKE solver
0027 path_path = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0028 %
0029 id_equil = 'TScirc';%For plasma equilibrium
0030 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0031 %
0032 id_dkeparam = 'BOUNCE3D_MUMPS';%For DKE code parameters
0033 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0034 %
0035 id_display = 'NO_DISPLAY';%For output code display
0036 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0037 %
0038 id_ohm = '';%For Ohmic electric contribution
0039 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0040 %
0041 ids_wave = {'LH_karney_offaxis_2'};%For RF waves contribution (put all the type of waves needed)
0042 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0043 %
0044 id_transpfaste = 'magneticturbulence';%For fast electron radial transport
0045 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0046 %
0047 id_ripple = '';%For fast electron magnetic ripple losses
0048 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0049 %
0050 %************************************************************************************************************************************
0051 %************************************************************************************************************************************
0052 %************************************************************************************************************************************
0053 %
0054 [dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple] = load_structures_yp('dkepath',id_path,path_path,'equil',id_equil,path_equil,'dkeparam',id_dkeparam,path_dkeparam,'dkedisplay',id_display,path_display,'ohm',id_ohm,path_ohm,'waves',ids_wave,paths_wave,'transpfaste',id_transpfaste,path_transpfaste,'ripple',id_ripple,path_ripple);
0055 %
0056 %************************************************************************************************************************************
0057 %
0058 if exist('dmumpsmex');dkeparam.invproc = -2;end
0059 %
0060 dkeparam.psin_S = psin_S;
0061 dkeparam.rho_S = rho_S;
0062 %
0063 % Dpsipsi scan (0 -> 20 m2/s)
0064 %
0065 transpfaste.Dr0 = 0;
0066 [Znorm0,Zcurr0,ZP0_0,dke_out0,radialDKE,equilDKE,momentumDKE,gridDKE,Zmomcoef,Zbouncecoef,Zmripple,mksa,XXsinksource] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0067 disp('D=0 m2/s done');
0068 transpfaste.Dr0 = 0.1;
0069 [Znorm01,Zcurr01,ZP0_01,dke_out01] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0070 disp('D=0.1 m2/s done');
0071 transpfaste.Dr0 = 0.2;
0072 %[Znorm02,Zcurr02,ZP0_02,dke_out02] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0073 %disp('D=0.2 m2/s done');
0074 %transpfaste.Dr0 = 0.5;
0075 %[Znorm05,Zcurr05,ZP0_05,dke_out05] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0076 %disp('D=0.5 m2/s done');
0077 transpfaste.Dr0 = 1;
0078 [Znorm1,Zcurr1,ZP0_1,dke_out1] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0079 disp('D=1 m2/s done');
0080 %transpfaste.Dr0 = 2;
0081 %[Znorm2,Zcurr2,ZP0_2,dke_out2] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0082 %disp('D=2 m2/s done');
0083 %transpfaste.Dr0 = 5;
0084 %[Znorm5,Zcurr5,ZP0_5,dke_out5] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0085 %disp('D=5 m2/s done');
0086 transpfaste.Dr0 = 10;
0087 [Znorm10,Zcurr10,ZP0_10,dke_out10] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0088 disp('D=10 m2/s done');
0089 transpfaste.Dr0= 20;
0090 [Znorm20,Zcurr20,ZP0_20,dke_out20] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0091 disp('D=20 m2/s done');
0092 %[Znorm50,Zcurr50,ZP0_50,dke_out50] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0093 %disp('D=50 m2/s done');
0094 %[Znorm100,Zcurr100,ZP0_100,dke_out100] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0095 %disp('D=100 m2/s done');%
0096 %
0097 % Nprocessor scan for Dpsipsi = 10 m2/s (1-> 10)
0098 %
0099 transpfaste.Dr0 = 10;
0100 %
0101 %Compare processing performances
0102 %
0103 dkeparam.invproc = -1;%MUMPS sequential or parallel.
0104 %
0105 Znorm10_4 = Znorm10;%Calculations already done
0106 Zcurr10_4 = Zcurr10;
0107 ZP0_10_4 = ZP0_10;
0108 dke_out10_4 = dke_out10;
0109 %
0110 dkeparam.MUMPSparam.nproc = 1;%Number of processors
0111 [Znorm10_1,Zcurr10_1,ZP0_10_1,dke_out10_1] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0112 disp('Nproc = 1 done');
0113 dkeparam.MUMPSparam.nproc = 2;%Number of processors
0114 [Znorm10_2,Zcurr10_2,ZP0_10_2,dke_out10_2] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0115 disp('Nprocs = 2 done');
0116 dkeparam.MUMPSparam.nproc = 3;%Number of processors
0117 [Znorm10_3,Zcurr10_3,ZP0_10_3,dke_out10_3] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0118 disp('Nprocs = 3 done');
0119 dkeparam.MUMPSparam.nproc = 5;%Number of processors
0120 [Znorm10_5,Zcurr10_5,ZP0_10_5,dke_out10_5] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0121 disp('Nproc = 5 done');
0122 dkeparam.MUMPSparam.nproc = 6;%Number of processors
0123 [Znorm10_6,Zcurr10_6,ZP0_10_6,dke_out10_6] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0124 disp('Nproc = 6 done');
0125 dkeparam.MUMPSparam.nproc = 7;%Number of processors
0126 [Znorm10_7,Zcurr10_7,ZP0_10_7,dke_out10_7] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0127 disp('Nproc = 7 done');
0128 dkeparam.MUMPSparam.nproc = 8;%Number of processors
0129 [Znorm10_8,Zcurr10_8,ZP0_10_8,dke_out10_8] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0130 disp('Nproc = 8 done');
0131 dkeparam.MUMPSparam.nproc = 9;%Number of processors
0132 [Znorm10_9,Zcurr10_9,ZP0_10_9,dke_out10_9] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0133 disp('Nproc = 9 done');
0134 dkeparam.MUMPSparam.nproc = 10;%Number of processors
0135 [Znorm10_10,Zcurr10_10,ZP0_10_10,dke_out10_10] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0136 disp('Nproc = 10 done'); 
0137 %
0138 Dpsipsi=[0,0.1,1,10];
0139 tCPU = [dke_out0.time_out,dke_out01.time_out,dke_out1.time_out,dke_out10.time_out];
0140 tCPU1 = [dke_out10_1.time_out,dke_out10_2.time_out,dke_out10_3.time_out,dke_out10_4.time_out,dke_out10_5.time_out,dke_out10_6.time_out,dke_out10_7.time_out,dke_out10_8.time_out,dke_out10_9.time_out,dke_out10_10.time_out];
0141 %
0142 vparmin = 1/(waves{1}.rays{1}.sNpar(1)+imag(waves{1}.rays{1}.sdNpar(1))/2)/mksa.betath_ref;
0143 vparmax = 1/(waves{1}.rays{1}.sNpar(1)-imag(waves{1}.rays{1}.sdNpar(1))/2)/mksa.betath_ref;
0144 %
0145 figure(1),clf,graph1D_jd(Dpsipsi,tCPU,1,0,'D_{\psi\psi} (m^{2}/s)','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,[0.1,10],[0,max(tCPU)*1.2],'-','none','r',2,20,gca,1,0.6,0.6);
0146 print_jd(2,'tCPU_Dpsipsi','./figures',1)
0147 %
0148 figure(2),clf,graph1D_jd([1:10],tCPU1,0,0,'Number of processors','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2),', D_{\psi\psi}  = 10 (m^{2}/s)'],NaN,[0,length(tCPU1)+1],[0,max(tCPU1)*1.2],'-','none','r',2,20,gca,1,0.6,0.6);
0149 graph1D_jd([1:10],dke_out10.time_out*ones(1,10),0,0,'Number of processors','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2),', D_{\psi\psi}  = 10 (m^{2}/s)'],NaN,[0,length(tCPU1)+1],[0,max(tCPU1)*1.2],'--','none','b',1,20,gca,1,0.6,0.6);
0150 leg = legend('MUMPS seq.','MUMPS mex')
0151 %
0152 print_jd(2,'tCPUDpsipsi10_Nproc','./figures',2)
0153 %
0154 figure(3),clf,graph1D_jd(equilDKE.xrho,Zcurr0.x_tot_fsav,0,0,'\rho','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,'','','-','none','r',2,20,gca,1,0.6,0.6);hold on
0155 graph1D_jd(equilDKE.xrho,Zcurr01.x_tot_fsav,0,0,'\rho','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,'','','-','none','b',2,20,gca,1,0.6,0.6);hold on
0156 %graph1D_jd(equilDKE.xrho,Zcurr02.x_tot_fsav,0,0,'\rho','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,'','','-','none','g',1,20,gca,1,0.6,0.6);hold on
0157 %graph1D_jd(equilDKE.xrho,Zcurr05.x_tot_fsav,0,0,'\rho','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,'','','-','none','c',1,20,gca,1,0.6,0.6);hold on
0158 graph1D_jd(equilDKE.xrho,Zcurr1.x_tot_fsav,0,0,'\rho','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,'','','-','none','g',2,20,gca,1,0.6,0.6);hold on
0159 %graph1D_jd(equilDKE.xrho,Zcurr2.x_tot_fsav,0,0,'\rho','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,'','','-','none','k',1,20,gca,1,0.6,0.6);hold on
0160 %graph1D_jd(equilDKE.xrho,Zcurr5.x_tot_fsav,0,0,'\rho','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,'','','--','none','r',1,20,gca,1,0.6,0.6);hold on
0161 graph1D_jd(equilDKE.xrho,Zcurr10.x_tot_fsav,0,0,'\rho','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,'','','-','none','m',2,20,gca,1,0.6,0.6);hold on
0162 graph1D_jd(equilDKE.xrho,Zcurr20.x_tot_fsav,0,0,'\rho','CPU time (s)',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,'','','-','none','k',2,20,gca,1,0.6,0.6);hold on
0163 graph1D_jd(equilDKE.xrho,besselj(0,2.4*equilDKE.xrho)*Zcurr20.x_tot_fsav(1),0,0,'\rho','current density',['v_{||min} = ',num2str(vparmin),', v_{||max} = ',num2str(vparmax),', D_{|| ||LH} = ',num2str(2)],NaN,'','','none','+','r',2,20,gca,1,0.6,0.6);
0164 legend('D_{\psi\psi} = 0 (m^{2}/s)','D_{\psi\psi} = 0.1 (m^{2}/s)','D_{\psi\psi} = 1 (m^{2}/s)','D_{\psi\psi} = 10 (m^{2}/s)','D_{\psi\psi} = 20 (m^{2}/s)','J_{0}(2.4*\rho)')
0165 %
0166 print_jd(2,'Jprof_Dpsipsi','./figures',3)
0167 %
0168 %************************************************************************************************************************************
0169 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0170 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);