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 p_opt = 2;
0012 %
0013 permission = test_permissions_yp;
0014 %
0015 if ~permission 
0016     disp('Please move the script to a local folder where you have write permission before to run it')
0017     return;
0018 end
0019 %
0020 % ***********************This part must be specified by the user, run make files if necessary) *****************************
0021 %
0022 id_simul = 'LH_karney_transp_bounce_nomag';%Simulation ID
0023 path_simul = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0024 %
0025 psin_S = [];%Normalized poloidal flux grid where calculations are performed (0 < psin_S < 1) (If one value: local calculation only, not used if empty)
0026 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)
0027 %
0028 id_path = '';%For all paths used by DKE solver
0029 path_path = '.';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0030 %
0031 id_equil = 'TScirc';%For plasma equilibrium
0032 path_equil = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0033 %
0034 id_dkeparam = 'BOUNCE3D_MUMPS';%For DKE code parameters
0035 path_dkeparam = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0036 %
0037 id_display = 'NO_DISPLAY';%For output code display
0038 path_display = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0039 %
0040 id_ohm = '';%For Ohmic electric contribution
0041 path_ohm = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0042 %
0043 ids_wave = {'LH_karney_offaxis_2'};%For RF waves contribution (put all the type of waves needed)
0044 paths_wave = {''};%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0045 %
0046 id_transpfaste = 'nomagneticturbulence';%For fast electron radial transport
0047 path_transpfaste = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0048 %
0049 id_ripple = '';%For fast electron magnetic ripple losses
0050 path_ripple = '';%if nothing is specified, the working directory is first used and then MatLab is looking in all the path
0051 %
0052 %************************************************************************************************************************************
0053 %************************************************************************************************************************************
0054 %************************************************************************************************************************************
0055 %
0056 [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);
0057 %
0058 %************************************************************************************************************************************
0059 %
0060 if exist('dmumpsmex');dkeparam.invproc = -2;end
0061 %
0062 dkeparam.psin_S = psin_S;
0063 dkeparam.rho_S = rho_S;
0064 %
0065 % Dpsipsi scan (0 -> 100 m2/s
0066 %
0067 transpfaste.Dr0 = 0;
0068 [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,[],[]);
0069 disp('D=0 m2/s done');
0070 transpfaste.Dr0 = 0.1;
0071 [Znorm01,Zcurr01,ZP0_01,dke_out01] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0072 disp('D=0.1 m2/s done');
0073 transpfaste.Dr0 = 0.2;
0074 %[Znorm02,Zcurr02,ZP0_02,dke_out02] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0075 %disp('D=0.2 m2/s done');
0076 %transpfaste.Dr0 = 0.5;
0077 %[Znorm05,Zcurr05,ZP0_05,dke_out05] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0078 %disp('D=0.5 m2/s done');
0079 transpfaste.Dr0 = 1;
0080 [Znorm1,Zcurr1,ZP0_1,dke_out1] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0081 disp('D=1 m2/s done');
0082 %transpfaste.Dr0 = 2;
0083 %[Znorm2,Zcurr2,ZP0_2,dke_out2] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0084 %disp('D=2 m2/s done');
0085 %transpfaste.Dr0 = 5;
0086 %[Znorm5,Zcurr5,ZP0_5,dke_out5] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0087 %disp('D=5 m2/s done');
0088 transpfaste.Dr0 = 10;
0089 [Znorm10,Zcurr10,ZP0_10,dke_out10] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0090 disp('D=10 m2/s done');
0091 transpfaste.Dr0= 20;
0092 [Znorm20,Zcurr20,ZP0_20,dke_out20] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0093 disp('D=20 m2/s done');
0094 %[Znorm50,Zcurr50,ZP0_50,dke_out50] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0095 %disp('D=50 m2/s done');
0096 %[Znorm100,Zcurr100,ZP0_100,dke_out100] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0097 %disp('D=100 m2/s done');%
0098 %
0099 % Nprocessor scan for Dpsipsi = 10 m2/s (1-> 10)
0100 %
0101 transpfaste.Dr0 = 10;
0102 %
0103 dkeparam.invproc = -1;
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 %
0111 %TODO: dkeparam.MUMPSparam.nproc not in operation
0112 %
0113 dkeparam.MUMPSparam.nproc = 1;%Number of processors
0114 [Znorm10_1,Zcurr10_1,ZP0_10_1,dke_out10_1] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0115 disp('Nproc = 1 done');
0116 dkeparam.MUMPSparam.nproc = 2;%Number of processors
0117 [Znorm10_2,Zcurr10_2,ZP0_10_2,dke_out10_2] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0118 disp('Nprocs = 2 done');
0119 dkeparam.MUMPSparam.nproc = 3;%Number of processors
0120 [Znorm10_3,Zcurr10_3,ZP0_10_3,dke_out10_3] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0121 disp('Nprocs = 3 done');
0122 dkeparam.MUMPSparam.nproc = 5;%Number of processors
0123 [Znorm10_5,Zcurr10_5,ZP0_10_5,dke_out10_5] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0124 disp('Nproc = 5 done');
0125 dkeparam.MUMPSparam.nproc = 6;%Number of processors
0126 [Znorm10_6,Zcurr10_6,ZP0_10_6,dke_out10_6] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0127 disp('Nproc = 6 done');
0128 dkeparam.MUMPSparam.nproc = 7;%Number of processors
0129 [Znorm10_7,Zcurr10_7,ZP0_10_7,dke_out10_7] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0130 disp('Nproc = 7 done');
0131 dkeparam.MUMPSparam.nproc = 8;%Number of processors
0132 [Znorm10_8,Zcurr10_8,ZP0_10_8,dke_out10_8] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0133 disp('Nproc = 8 done');
0134 dkeparam.MUMPSparam.nproc = 9;%Number of processors
0135 [Znorm10_9,Zcurr10_9,ZP0_10_9,dke_out10_9] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0136 disp('Nproc = 9 done');
0137 dkeparam.MUMPSparam.nproc = 10;%Number of processors
0138 [Znorm10_10,Zcurr10_10,ZP0_10_10,dke_out10_10] = main_dke_yp(id_simul,dkepath,equil,dkeparam,dkedisplay,ohm,waves,transpfaste,ripple,[],[]);
0139 disp('Nproc = 10 done');  
0140 %
0141 Dpsipsi=[0,0.1,1,10];
0142 tCPU = [dke_out0.time_out,dke_out01.time_out,dke_out1.time_out,dke_out10.time_out];
0143 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];
0144 %
0145 vparmin = 1/(waves{1}.rays{1}.sNpar(1)+imag(waves{1}.rays{1}.sdNpar(1))/2)/mksa.betath_ref;
0146 vparmax = 1/(waves{1}.rays{1}.sNpar(1)-imag(waves{1}.rays{1}.sdNpar(1))/2)/mksa.betath_ref;
0147 %
0148 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,'','','-','none','r',2,20,gca,1,0.6,0.6);
0149 axis([0.1,10,0,300])
0150 %
0151 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,'','','-','none','r',2,20,gca,1,0.6,0.6);
0152 axis([1,10,0,400])
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(p_opt,'tCPU_Dpsipsi','./figures',1)
0167 print_jd(p_opt,'tCPUDpsipsi10_Nproc','./figures',2)
0168 print_jd(p_opt,'Jprof_Dpsipsi','./figures',3)
0169 %
0170 %************************************************************************************************************************************
0171 eval(['save ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);
0172 info_dke_yp(2,['Data saved in ',path_simul,'DKE_RESULTS_',id_equil,'_',id_simul,'.mat']);